10 #ifndef __ENTITIES_FIELD_DATA_HPP__
11 #define __ENTITIES_FIELD_DATA_HPP__
13 using namespace boost::numeric;
19 FEDofEntity *, std::allocator<FEDofEntity *>
23 using VectorDofs = ublas::vector<FEDofEntity *, DofsAllocator>;
44 std::bitset<LASTBASE>
bAse;
48 std::array<std::bitset<LASTSPACE>, MBMAXTYPE>
50 std::array<std::bitset<LASTBASE>, MBMAXTYPE>
52 std::array<std::bitset<LASTBASE>,
LASTSPACE>
54 std::array<std::bitset<LASTBASE>,
LASTSPACE>
56 std::array<boost::ptr_vector<EntData>, MBMAXTYPE>
115 const boost::shared_ptr<EntitiesFieldData> &data_ptr);
119 const boost::shared_ptr<EntitiesFieldData>
dataPtr;
143 EntData(
const bool allocate_base_matrices =
true);
154 virtual int getSense()
const;
160 inline const VectorInt &getIndices()
const;
166 inline const VectorInt &getLocalIndices()
const;
171 inline int &getSense();
192 inline const VectorDofs &getFieldDofs()
const;
204 virtual std::vector<BitRefLevel> &getEntDataBitRefLevel();
217 template <
int Tensor_Dim>
219 getFTensor1FieldData();
233 template <
int Tensor_Dim0,
int Tensor_Dim1>
235 Tensor_Dim0, Tensor_Dim1>
236 getFTensor2FieldData();
250 template <
int Tensor_Dim>
254 getFTensor2SymmetricFieldData();
286 virtual boost::shared_ptr<MatrixDouble> &
293 virtual boost::shared_ptr<MatrixDouble> &
299 virtual boost::shared_ptr<MatrixDouble> &
426 const int gg,
const int nb_base_functions);
439 inline const VectorAdaptor getN(
const int gg,
const int nb_base_functions);
455 const int nb_base_functions);
469 const int nb_base_functions);
494 template <
int DIM>
inline const MatrixAdaptor getVectorN(
const int gg);
503 template <
int DIM0,
int DIM1>
513 template <
int DIM0,
int DIM1>
523 template <
int DIM0,
int DIM1>
525 const int dof,
const int gg);
533 template <
int DIM0,
int DIM1>
534 inline const MatrixAdaptor getVectorDiffN(
const int dof,
const int gg);
593 getFTensor0N(
const int gg,
const int bb);
612 template <
int Tensor_Dim>
632 template <
int Tensor_Dim>
654 template <
int Tensor_Dim>
677 template <
int Tensor_Dim>
679 getFTensor1DiffN(
const int gg,
const int bb);
706 template <
int Tensor_Dim>
730 template <
int Tensor_Dim>
auto getFTensor1N();
734 template <
int Tensor_Dim0,
int Tensor_Dim1>
736 Tensor_Dim0, Tensor_Dim1>
742 template <
int Tensor_Dim0,
int Tensor_Dim1>
744 Tensor_Dim0, Tensor_Dim1>
749 template <
int Tensor_Dim0,
int Tensor_Dim1>
751 Tensor_Dim0, Tensor_Dim1>
753 return getFTensor2DiffN<Tensor_Dim0, Tensor_Dim1>(bAse);
759 template <
int Tensor_Dim0,
int Tensor_Dim1>
761 Tensor_Dim0, Tensor_Dim1>
763 return getFTensor2DiffN<Tensor_Dim0, Tensor_Dim1>(bAse, gg, bb);
768 template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2>
771 Tensor_Dim0, Tensor_Dim1, Tensor_Dim2>
776 template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2>
779 Tensor_Dim0, Tensor_Dim1, Tensor_Dim2>
781 return getFTensor3Diff2N<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2>(bAse);
798 template <
int Tensor_Dim>
816 template <
int Tensor_Dim>
817 inline auto getFTensor1N(
const int gg,
const int bb);
845 template <
int Tensor_Dim0,
int Tensor_Dim1>
847 Tensor_Dim0, Tensor_Dim1>
872 template <
int Tensor_Dim0,
int Tensor_Dim1>
auto getFTensor2N();
900 template <
int Tensor_Dim0,
int Tensor_Dim1>
902 Tensor_Dim0, Tensor_Dim1>
926 template <
int Tensor_Dim0,
int Tensor_Dim1>
927 auto getFTensor2N(
const int gg,
const int bb);
935 friend std::ostream &
operator<<(std::ostream &os,
964 virtual boost::shared_ptr<MatrixInt> &
965 getBBAlphaIndicesSharedPtr(
const std::string &
field_name);
970 virtual boost::shared_ptr<MatrixDouble> &
971 getBBNSharedPtr(
const std::string &
field_name);
976 virtual const boost::shared_ptr<MatrixDouble> &
977 getBBNSharedPtr(
const std::string &
field_name)
const;
982 virtual boost::shared_ptr<MatrixDouble> &
983 getBBDiffNSharedPtr(
const std::string &
field_name);
988 virtual const boost::shared_ptr<MatrixDouble> &
989 getBBDiffNSharedPtr(
const std::string &
field_name)
const;
991 virtual std::map<std::string, boost::shared_ptr<MatrixInt>> &
992 getBBAlphaIndicesMap();
999 virtual std::map<std::string, boost::shared_ptr<MatrixDouble>> &getBBNMap();
1007 virtual std::map<std::string, boost::shared_ptr<MatrixDouble>> &
1016 virtual boost::shared_ptr<MatrixInt> &
1017 getBBAlphaIndicesByOrderSharedPtr(
const size_t o);
1025 virtual boost::shared_ptr<MatrixDouble> &
1026 getBBNByOrderSharedPtr(
const size_t o);
1034 virtual boost::shared_ptr<MatrixDouble> &
1035 getBBDiffNByOrderSharedPtr(
const size_t o);
1039 virtual std::array<boost::shared_ptr<MatrixInt>, MaxBernsteinBezierOrder> &
1040 getBBAlphaIndicesByOrderArray();
1042 virtual std::array<boost::shared_ptr<MatrixDouble>, MaxBernsteinBezierOrder> &
1043 getBBNByOrderArray();
1045 virtual std::array<boost::shared_ptr<MatrixDouble>, MaxBernsteinBezierOrder> &
1046 getBBDiffNByOrderArray();
1078 std::array<std::array<boost::shared_ptr<MatrixDouble>,
LASTBASE>,
1083 std::array<boost::shared_ptr<MatrixDouble>,
LASTBASE>
1088 std::map<std::string, boost::shared_ptr<MatrixDouble>>
bbN;
1089 std::map<std::string, boost::shared_ptr<MatrixDouble>>
bbDiffN;
1090 std::map<std::string, boost::shared_ptr<MatrixInt>>
1093 std::array<boost::shared_ptr<MatrixDouble>, MaxBernsteinBezierOrder>
1095 std::array<boost::shared_ptr<MatrixDouble>, MaxBernsteinBezierOrder>
1097 std::array<boost::shared_ptr<MatrixInt>, MaxBernsteinBezierOrder>
1134 const boost::shared_ptr<EntitiesFieldData::EntData> &ent_data_ptr);
1136 int getSense()
const;
1139 std::vector<BitRefLevel> &getEntDataBitRefLevel();
1141 boost::shared_ptr<MatrixDouble> &
1145 boost::shared_ptr<MatrixDouble> &
1148 boost::shared_ptr<MatrixDouble> &
1151 const boost::shared_ptr<MatrixDouble> &
1154 const boost::shared_ptr<MatrixDouble> &
1157 inline boost::shared_ptr<MatrixDouble> &
1160 inline boost::shared_ptr<MatrixDouble> &
1163 boost::shared_ptr<MatrixInt> &
1164 getBBAlphaIndicesSharedPtr(
const std::string &
field_name);
1169 boost::shared_ptr<MatrixDouble> &
1170 getBBNSharedPtr(
const std::string &
field_name);
1175 const boost::shared_ptr<MatrixDouble> &
1176 getBBNSharedPtr(
const std::string &
field_name)
const;
1181 boost::shared_ptr<MatrixDouble> &
1182 getBBDiffNSharedPtr(
const std::string &
field_name);
1187 const boost::shared_ptr<MatrixDouble> &
1188 getBBDiffNSharedPtr(
const std::string &
field_name)
const;
1204 const VectorInt &EntitiesFieldData::EntData::getIndices()
const {
1209 EntitiesFieldData::EntData::getIndicesUpToOrder(
int order) {
1210 unsigned int size = 0;
1211 if (
auto dof = dOfs[0]) {
1212 size = dof->getOrderNbDofs(
order) * dof->getNbOfCoeffs();
1213 size = size < iNdices.size() ? size : iNdices.size();
1215 int *data = &*iNdices.data().begin();
1216 return VectorIntAdaptor(size, ublas::shallow_array_adaptor<int>(size, data));
1219 const VectorInt &EntitiesFieldData::EntData::getLocalIndices()
const {
1220 return localIndices;
1224 EntitiesFieldData::EntData::getLocalIndicesUpToOrder(
int order) {
1225 unsigned int size = 0;
1226 if (
auto dof = dOfs[0]) {
1227 size = dof->getOrderNbDofs(
order) * dof->getNbOfCoeffs();
1228 size = size < localIndices.size() ? size : localIndices.size();
1230 int *data = &*localIndices.data().begin();
1231 return VectorIntAdaptor(size, ublas::shallow_array_adaptor<int>(size, data));
1234 int &EntitiesFieldData::EntData::getSense() {
return sEnse; }
1238 VectorInt &EntitiesFieldData::EntData::getIndices() {
return iNdices; }
1240 VectorInt &EntitiesFieldData::EntData::getLocalIndices() {
1241 return localIndices;
1249 EntitiesFieldData::EntData::getFieldDataUpToOrder(
int order) {
1250 unsigned int size = 0;
1251 if (
auto dof = dOfs[0]) {
1252 size = dof->getOrderNbDofs(
order) * dof->getNbOfCoeffs();
1253 size = size < fieldData.size() ? size : fieldData.size();
1255 double *data = &*fieldData.data().begin();
1259 const VectorDofs &EntitiesFieldData::EntData::getFieldDofs()
const {
1263 VectorDofs &EntitiesFieldData::EntData::getFieldDofs() {
return dOfs; }
1265 VectorDouble &EntitiesFieldData::EntData::getFieldData() {
return fieldData; }
1268 return fieldEntities;
1272 EntitiesFieldData::EntData::getFieldEntities()
const {
1273 return fieldEntities;
1276 template <
int Tensor_Dim>
1278 EntitiesFieldData::EntData::getFTensor1FieldData() {
1279 std::stringstream s;
1280 s <<
"Not implemented for this dimension dim = " << Tensor_Dim;
1284 template <
int Tensor_Dim0,
int Tensor_Dim1>
1286 Tensor_Dim0, Tensor_Dim1>
1287 EntitiesFieldData::EntData::getFTensor2FieldData() {
1288 std::stringstream s;
1289 s <<
"Not implemented for this dimension dim0 = " << Tensor_Dim0;
1290 s <<
" and dim1 " << Tensor_Dim1;
1294 template <
int Tensor_Dim>
1296 FTensor::PackPtr<
double *, (Tensor_Dim * (Tensor_Dim + 1)) / 2>, Tensor_Dim>
1297 EntitiesFieldData::EntData::getFTensor2SymmetricFieldData() {
1298 std::stringstream s;
1299 s <<
"Not implemented for this dimension dim = " << Tensor_Dim;
1305 FieldSpace &EntitiesFieldData::EntData::getSpace() {
return sPace; }
1309 return *(getNSharedPtr(base));
1320 return *(getDiffNSharedPtr(base));
1327 if (!getNSharedPtr(base, derivative)) {
1329 "Ptr to base %s functions derivative %d is null",
1334 return *(getNSharedPtr(base, derivative));
1338 EntitiesFieldData::EntData::getDiffN(
const std::string &
field_name) {
1342 MatrixDouble &EntitiesFieldData::EntData::getDiffN() {
return getDiffN(bAse); }
1346 return getN(bAse, derivative);
1352 int size = getN(base).size2();
1353 double *data = &getN(base)(gg, 0);
1354 return VectorAdaptor(size, ublas::shallow_array_adaptor<double>(size, data));
1358 return getN(bAse, gg);
1367 if (getN(base).size1() == getDiffN(base).size1()) {
1368 int size = getN(base).size2();
1369 int dim = getDiffN(base).size2() / size;
1370 double *data = &getDiffN(base)(gg, 0);
1372 getN(base).size2(), dim,
1373 ublas::shallow_array_adaptor<double>(getDiffN(base).size2(), data));
1378 getN(base).size1(), getN(base).size2(),
1379 ublas::shallow_array_adaptor<double>(getDiffN(base).data().size(),
1380 &getDiffN(base).data()[0]));
1385 return getDiffN(bAse, gg);
1390 const int gg,
const int nb_base_functions) {
1391 (void)getN()(gg, nb_base_functions -
1393 double *data = &getN(base)(gg, 0);
1394 return VectorAdaptor(nb_base_functions, ublas::shallow_array_adaptor<double>(
1395 nb_base_functions, data));
1399 EntitiesFieldData::EntData::getN(
const int gg,
const int nb_base_functions) {
1400 return getN(bAse, gg, nb_base_functions);
1406 const int nb_base_functions) {
1410 if (getN(base).size1() == getDiffN(base).size1()) {
1411 (void)getN(base)(gg,
1414 int dim = getDiffN(base).size2() / getN(base).size2();
1415 double *data = &getDiffN(base)(gg, 0);
1417 nb_base_functions, dim,
1418 ublas::shallow_array_adaptor<double>(dim * nb_base_functions, data));
1423 getN(base).size1(), getN(base).size2(),
1424 ublas::shallow_array_adaptor<double>(getDiffN(base).data().size(),
1425 &getDiffN(base).data()[0]));
1430 EntitiesFieldData::EntData::getDiffN(
const int gg,
1431 const int nb_base_functions) {
1432 return getDiffN(bAse, gg, nb_base_functions);
1439 if (PetscUnlikely(getN(base).size2() % DIM)) {
1443 const int nb_base_functions = getN(base).size2() / DIM;
1444 double *data = &getN(base)(gg, 0);
1446 nb_base_functions, DIM,
1447 ublas::shallow_array_adaptor<double>(DIM * nb_base_functions, data));
1452 return getVectorN<DIM>(bAse, gg);
1455 template <
int DIM0,
int DIM1>
1459 if (PetscUnlikely(getDiffN(base).size2() % (DIM0 *
DIM1))) {
1463 const int nb_base_functions = getN(base).size2() / (DIM0 *
DIM1);
1464 double *data = &getN(base)(gg, 0);
1466 ublas::shallow_array_adaptor<double>(
1467 DIM0 *
DIM1 * nb_base_functions, data));
1470 template <
int DIM0,
int DIM1>
1472 return getVectorDiffN<DIM0, DIM1>(bAse, gg);
1475 template <
int DIM0,
int DIM1>
1478 const int dof,
const int gg) {
1480 &EntitiesFieldData::EntData::getDiffN(base)(gg, DIM0 *
DIM1 * dof);
1482 ublas::shallow_array_adaptor<double>(DIM0 *
DIM1, data));
1485 template <
int DIM0,
int DIM1>
1488 return getVectorDiffN<DIM0, DIM1>(bAse, dof, gg);
1493 double *ptr = &*getN(base).data().begin();
1498 EntitiesFieldData::EntData::getFTensor0N() {
1499 return getFTensor0N(bAse);
1504 const int gg,
const int bb) {
1505 double *ptr = &getN(base)(gg, bb);
1510 EntitiesFieldData::EntData::getFTensor0N(
const int gg,
const int bb) {
1511 return getFTensor0N(bAse, gg, bb);
1514 template <
int Tensor_Dim>
auto EntitiesFieldData::EntData::getFTensor1N() {
1515 return getFTensor1N<Tensor_Dim>(bAse);
1518 template <
int Tensor_Dim>
1519 auto EntitiesFieldData::EntData::getFTensor1N(
const int gg,
const int bb) {
1520 return getFTensor1N<Tensor_Dim>(bAse, gg, bb);
1523 template <
int Tensor_Dim0,
int Tensor_Dim1>
1524 auto EntitiesFieldData::EntData::getFTensor2N() {
1525 return getFTensor2N<Tensor_Dim0, Tensor_Dim1>(bAse);
1528 template <
int Tensor_Dim0,
int Tensor_Dim1>
1529 auto EntitiesFieldData::EntData::getFTensor2N(
const int gg,
const int bb) {
1530 return getFTensor2N<Tensor_Dim0, Tensor_Dim1>(bAse, gg, bb);
1537 VectorInt &EntitiesFieldData::EntData::getBBNodeOrder() {
return bbNodeOrder; }
1539 MatrixInt &EntitiesFieldData::EntData::getBBAlphaIndices() {
1540 return *getBBAlphaIndicesSharedPtr(bbFieldName);
1549 boost::shared_ptr<MatrixDouble> &
1550 DerivedEntitiesFieldData::DerivedEntData::getDerivedNSharedPtr(
1555 boost::shared_ptr<MatrixDouble> &
1556 DerivedEntitiesFieldData::DerivedEntData::getDerivedDiffNSharedPtr(
1578 template <
typename T = EntityStorage>
1581 const double *ptr, InsertMode iora) {
1582 static_assert(!std::is_same<T, T>::value,
1583 "VecSetValues value for this data storage is not implemented");
1590 const double *ptr, InsertMode iora) {
1610 template <
typename T = EntityStorage>
1614 return VecSetValues<T>(V, data, &*vec.data().begin(), iora);
1633 template <
typename T = EntityStorage>
1637 const double *ptr, InsertMode iora) {
1638 static_assert(!std::is_same<T, T>::value,
1639 "MatSetValues value for this data storage is not implemented");
1659 template <
typename T = EntityStorage>
1664 return MatSetValues<T>(
M, row_data, col_data, &*mat.data().begin(), iora);
1671 const double *ptr, InsertMode iora) {
1688 const int gg,
const int bb);
1702 EntitiesFieldData::EntData::getFTensor1DiffN<3>(
1706 EntitiesFieldData::EntData::getFTensor1DiffN<3>();
1710 EntitiesFieldData::EntData::getFTensor1DiffN<2>(
1714 EntitiesFieldData::EntData::getFTensor1DiffN<2>();
1722 const int gg,
const int bb);
1736 EntitiesFieldData::EntData::getFTensor1FieldData<3>();
1740 EntitiesFieldData::EntData::getFTensor1FieldData<2>();
1744 EntitiesFieldData::EntData::getFTensor1FieldData<1>();
1748 EntitiesFieldData::EntData::getFTensor2FieldData<1, 1>();
1752 EntitiesFieldData::EntData::getFTensor2FieldData<1, 2>();
1756 EntitiesFieldData::EntData::getFTensor2FieldData<1, 3>();
1760 EntitiesFieldData::EntData::getFTensor2FieldData<2, 2>();
1764 EntitiesFieldData::EntData::getFTensor2FieldData<3, 3>();
1768 EntitiesFieldData::EntData::getFTensor2SymmetricFieldData<3>();
1772 EntitiesFieldData::EntData::getFTensor2SymmetricFieldData<2>();
1788 #endif //__ENTITIES_FIELD_DATA_HPP__