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<LASTSPACE>
sPace;
45 std::bitset<LASTBASE>
bAse;
49 std::array<std::bitset<LASTSPACE>, MBMAXTYPE>
51 std::array<std::bitset<LASTBASE>, MBMAXTYPE>
53 std::array<std::bitset<LASTBASE>,
LASTSPACE>
55 std::array<boost::ptr_vector<EntData>, MBMAXTYPE>
114 const boost::shared_ptr<EntitiesFieldData> &data_ptr);
118 const boost::shared_ptr<EntitiesFieldData>
dataPtr;
142 EntData(
const bool allocate_base_matrices =
true);
153 virtual int getSense()
const;
159 inline const VectorInt &getIndices()
const;
165 inline const VectorInt &getLocalIndices()
const;
170 inline int &getSense();
191 inline const VectorDofs &getFieldDofs()
const;
203 virtual std::vector<BitRefLevel> &getEntDataBitRefLevel();
216 template <
int Tensor_Dim>
218 getFTensor1FieldData();
232 template <
int Tensor_Dim0,
int Tensor_Dim1>
234 Tensor_Dim0, Tensor_Dim1>
235 getFTensor2FieldData();
249 template <
int Tensor_Dim>
253 getFTensor2SymmetricFieldData();
285 virtual boost::shared_ptr<MatrixDouble> &
292 virtual boost::shared_ptr<MatrixDouble> &
298 virtual boost::shared_ptr<MatrixDouble> &
425 const int gg,
const int nb_base_functions);
438 inline const VectorAdaptor getN(
const int gg,
const int nb_base_functions);
454 const int nb_base_functions);
468 const int nb_base_functions);
493 template <
int DIM>
inline const MatrixAdaptor getVectorN(
const int gg);
502 template <
int DIM0,
int DIM1>
512 template <
int DIM0,
int DIM1>
522 template <
int DIM0,
int DIM1>
524 const int dof,
const int gg);
532 template <
int DIM0,
int DIM1>
533 inline const MatrixAdaptor getVectorDiffN(
const int dof,
const int gg);
592 getFTensor0N(
const int gg,
const int bb);
611 template <
int Tensor_Dim>
631 template <
int Tensor_Dim>
653 template <
int Tensor_Dim>
676 template <
int Tensor_Dim>
678 getFTensor1DiffN(
const int gg,
const int bb);
705 template <
int Tensor_Dim>
729 template <
int Tensor_Dim>
auto getFTensor1N();
733 template <
int Tensor_Dim0,
int Tensor_Dim1>
735 Tensor_Dim0, Tensor_Dim1>
741 template <
int Tensor_Dim0,
int Tensor_Dim1>
743 Tensor_Dim0, Tensor_Dim1>
748 template <
int Tensor_Dim0,
int Tensor_Dim1>
750 Tensor_Dim0, Tensor_Dim1>
752 return getFTensor2DiffN<Tensor_Dim0, Tensor_Dim1>(bAse);
758 template <
int Tensor_Dim0,
int Tensor_Dim1>
760 Tensor_Dim0, Tensor_Dim1>
762 return getFTensor2DiffN<Tensor_Dim0, Tensor_Dim1>(bAse, gg, bb);
767 template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2>
770 Tensor_Dim0, Tensor_Dim1, Tensor_Dim2>
775 template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2>
778 Tensor_Dim0, Tensor_Dim1, Tensor_Dim2>
780 return getFTensor3Diff2N<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2>(bAse);
797 template <
int Tensor_Dim>
815 template <
int Tensor_Dim>
816 inline auto getFTensor1N(
const int gg,
const int bb);
844 template <
int Tensor_Dim0,
int Tensor_Dim1>
846 Tensor_Dim0, Tensor_Dim1>
871 template <
int Tensor_Dim0,
int Tensor_Dim1>
auto getFTensor2N();
899 template <
int Tensor_Dim0,
int Tensor_Dim1>
901 Tensor_Dim0, Tensor_Dim1>
925 template <
int Tensor_Dim0,
int Tensor_Dim1>
926 auto getFTensor2N(
const int gg,
const int bb);
934 friend std::ostream &
operator<<(std::ostream &os,
963 virtual boost::shared_ptr<MatrixInt> &
964 getBBAlphaIndicesSharedPtr(
const std::string &
field_name);
969 virtual boost::shared_ptr<MatrixDouble> &
970 getBBNSharedPtr(
const std::string &
field_name);
975 virtual const boost::shared_ptr<MatrixDouble> &
976 getBBNSharedPtr(
const std::string &
field_name)
const;
981 virtual boost::shared_ptr<MatrixDouble> &
982 getBBDiffNSharedPtr(
const std::string &
field_name);
987 virtual const boost::shared_ptr<MatrixDouble> &
988 getBBDiffNSharedPtr(
const std::string &
field_name)
const;
990 virtual std::map<std::string, boost::shared_ptr<MatrixInt>> &
991 getBBAlphaIndicesMap();
998 virtual std::map<std::string, boost::shared_ptr<MatrixDouble>> &getBBNMap();
1006 virtual std::map<std::string, boost::shared_ptr<MatrixDouble>> &
1015 virtual boost::shared_ptr<MatrixInt> &
1016 getBBAlphaIndicesByOrderSharedPtr(
const size_t o);
1024 virtual boost::shared_ptr<MatrixDouble> &
1025 getBBNByOrderSharedPtr(
const size_t o);
1033 virtual boost::shared_ptr<MatrixDouble> &
1034 getBBDiffNByOrderSharedPtr(
const size_t o);
1038 virtual std::array<boost::shared_ptr<MatrixInt>, MaxBernsteinBezierOrder> &
1039 getBBAlphaIndicesByOrderArray();
1041 virtual std::array<boost::shared_ptr<MatrixDouble>, MaxBernsteinBezierOrder> &
1042 getBBNByOrderArray();
1044 virtual std::array<boost::shared_ptr<MatrixDouble>, MaxBernsteinBezierOrder> &
1045 getBBDiffNByOrderArray();
1077 std::array<std::array<boost::shared_ptr<MatrixDouble>,
LASTBASE>,
1082 std::array<boost::shared_ptr<MatrixDouble>,
LASTBASE>
1087 std::map<std::string, boost::shared_ptr<MatrixDouble>>
bbN;
1088 std::map<std::string, boost::shared_ptr<MatrixDouble>>
bbDiffN;
1089 std::map<std::string, boost::shared_ptr<MatrixInt>>
1092 std::array<boost::shared_ptr<MatrixDouble>, MaxBernsteinBezierOrder>
1094 std::array<boost::shared_ptr<MatrixDouble>, MaxBernsteinBezierOrder>
1096 std::array<boost::shared_ptr<MatrixInt>, MaxBernsteinBezierOrder>
1131 const boost::shared_ptr<EntitiesFieldData::EntData> &ent_data_ptr);
1133 int getSense()
const;
1136 std::vector<BitRefLevel> &getEntDataBitRefLevel();
1138 boost::shared_ptr<MatrixDouble> &
1142 boost::shared_ptr<MatrixDouble> &
1145 boost::shared_ptr<MatrixDouble> &
1148 const boost::shared_ptr<MatrixDouble> &
1151 const boost::shared_ptr<MatrixDouble> &
1154 inline boost::shared_ptr<MatrixDouble> &
1157 inline boost::shared_ptr<MatrixDouble> &
1160 boost::shared_ptr<MatrixInt> &
1161 getBBAlphaIndicesSharedPtr(
const std::string &
field_name);
1166 boost::shared_ptr<MatrixDouble> &
1167 getBBNSharedPtr(
const std::string &
field_name);
1172 const boost::shared_ptr<MatrixDouble> &
1173 getBBNSharedPtr(
const std::string &
field_name)
const;
1178 boost::shared_ptr<MatrixDouble> &
1179 getBBDiffNSharedPtr(
const std::string &
field_name);
1184 const boost::shared_ptr<MatrixDouble> &
1185 getBBDiffNSharedPtr(
const std::string &
field_name)
const;
1201 const VectorInt &EntitiesFieldData::EntData::getIndices()
const {
1206 EntitiesFieldData::EntData::getIndicesUpToOrder(
int order) {
1207 unsigned int size = 0;
1208 if (
auto dof = dOfs[0]) {
1209 size = dof->getOrderNbDofs(
order) * dof->getNbOfCoeffs();
1210 size = size < iNdices.size() ? size : iNdices.size();
1212 int *data = &*iNdices.data().begin();
1213 return VectorIntAdaptor(size, ublas::shallow_array_adaptor<int>(size, data));
1216 const VectorInt &EntitiesFieldData::EntData::getLocalIndices()
const {
1217 return localIndices;
1221 EntitiesFieldData::EntData::getLocalIndicesUpToOrder(
int order) {
1222 unsigned int size = 0;
1223 if (
auto dof = dOfs[0]) {
1224 size = dof->getOrderNbDofs(
order) * dof->getNbOfCoeffs();
1225 size = size < localIndices.size() ? size : localIndices.size();
1227 int *data = &*localIndices.data().begin();
1228 return VectorIntAdaptor(size, ublas::shallow_array_adaptor<int>(size, data));
1231 int &EntitiesFieldData::EntData::getSense() {
return sEnse; }
1235 VectorInt &EntitiesFieldData::EntData::getIndices() {
return iNdices; }
1237 VectorInt &EntitiesFieldData::EntData::getLocalIndices() {
1238 return localIndices;
1246 EntitiesFieldData::EntData::getFieldDataUpToOrder(
int order) {
1247 unsigned int size = 0;
1248 if (
auto dof = dOfs[0]) {
1249 size = dof->getOrderNbDofs(
order) * dof->getNbOfCoeffs();
1250 size = size < fieldData.size() ? size : fieldData.size();
1252 double *data = &*fieldData.data().begin();
1256 const VectorDofs &EntitiesFieldData::EntData::getFieldDofs()
const {
1260 VectorDofs &EntitiesFieldData::EntData::getFieldDofs() {
return dOfs; }
1262 VectorDouble &EntitiesFieldData::EntData::getFieldData() {
return fieldData; }
1265 return fieldEntities;
1269 EntitiesFieldData::EntData::getFieldEntities()
const {
1270 return fieldEntities;
1273 template <
int Tensor_Dim>
1275 EntitiesFieldData::EntData::getFTensor1FieldData() {
1276 std::stringstream s;
1277 s <<
"Not implemented for this dimension dim = " << Tensor_Dim;
1281 template <
int Tensor_Dim0,
int Tensor_Dim1>
1283 Tensor_Dim0, Tensor_Dim1>
1284 EntitiesFieldData::EntData::getFTensor2FieldData() {
1285 std::stringstream s;
1286 s <<
"Not implemented for this dimension dim0 = " << Tensor_Dim0;
1287 s <<
" and dim1 " << Tensor_Dim1;
1291 template <
int Tensor_Dim>
1293 FTensor::PackPtr<
double *, (Tensor_Dim * (Tensor_Dim + 1)) / 2>, Tensor_Dim>
1294 EntitiesFieldData::EntData::getFTensor2SymmetricFieldData() {
1295 std::stringstream s;
1296 s <<
"Not implemented for this dimension dim = " << Tensor_Dim;
1302 FieldSpace &EntitiesFieldData::EntData::getSpace() {
return sPace; }
1306 return *(getNSharedPtr(base));
1317 return *(getDiffNSharedPtr(base));
1324 if (!getNSharedPtr(base, derivative)) {
1326 "Ptr to base %s functions derivative %d is null",
1331 return *(getNSharedPtr(base, derivative));
1335 EntitiesFieldData::EntData::getDiffN(
const std::string &
field_name) {
1339 MatrixDouble &EntitiesFieldData::EntData::getDiffN() {
return getDiffN(bAse); }
1343 return getN(bAse, derivative);
1349 int size = getN(base).size2();
1350 double *data = &getN(base)(gg, 0);
1351 return VectorAdaptor(size, ublas::shallow_array_adaptor<double>(size, data));
1355 return getN(bAse, gg);
1364 if (getN(base).size1() == getDiffN(base).size1()) {
1365 int size = getN(base).size2();
1366 int dim = getDiffN(base).size2() / size;
1367 double *data = &getDiffN(base)(gg, 0);
1369 getN(base).size2(), dim,
1370 ublas::shallow_array_adaptor<double>(getDiffN(base).size2(), data));
1375 getN(base).size1(), getN(base).size2(),
1376 ublas::shallow_array_adaptor<double>(getDiffN(base).data().size(),
1377 &getDiffN(base).data()[0]));
1382 return getDiffN(bAse, gg);
1387 const int gg,
const int nb_base_functions) {
1388 (void)getN()(gg, nb_base_functions -
1390 double *data = &getN(base)(gg, 0);
1391 return VectorAdaptor(nb_base_functions, ublas::shallow_array_adaptor<double>(
1392 nb_base_functions, data));
1396 EntitiesFieldData::EntData::getN(
const int gg,
const int nb_base_functions) {
1397 return getN(bAse, gg, nb_base_functions);
1403 const int nb_base_functions) {
1407 if (getN(base).size1() == getDiffN(base).size1()) {
1408 (void)getN(base)(gg,
1411 int dim = getDiffN(base).size2() / getN(base).size2();
1412 double *data = &getDiffN(base)(gg, 0);
1414 nb_base_functions, dim,
1415 ublas::shallow_array_adaptor<double>(dim * nb_base_functions, data));
1420 getN(base).size1(), getN(base).size2(),
1421 ublas::shallow_array_adaptor<double>(getDiffN(base).data().size(),
1422 &getDiffN(base).data()[0]));
1427 EntitiesFieldData::EntData::getDiffN(
const int gg,
1428 const int nb_base_functions) {
1429 return getDiffN(bAse, gg, nb_base_functions);
1436 if (PetscUnlikely(getN(base).size2() % DIM)) {
1440 const int nb_base_functions = getN(base).size2() / DIM;
1441 double *data = &getN(base)(gg, 0);
1443 nb_base_functions, DIM,
1444 ublas::shallow_array_adaptor<double>(DIM * nb_base_functions, data));
1449 return getVectorN<DIM>(bAse, gg);
1452 template <
int DIM0,
int DIM1>
1456 if (PetscUnlikely(getDiffN(base).size2() % (DIM0 *
DIM1))) {
1460 const int nb_base_functions = getN(base).size2() / (DIM0 *
DIM1);
1461 double *data = &getN(base)(gg, 0);
1463 ublas::shallow_array_adaptor<double>(
1464 DIM0 *
DIM1 * nb_base_functions, data));
1467 template <
int DIM0,
int DIM1>
1469 return getVectorDiffN<DIM0, DIM1>(bAse, gg);
1472 template <
int DIM0,
int DIM1>
1475 const int dof,
const int gg) {
1477 &EntitiesFieldData::EntData::getDiffN(base)(gg, DIM0 *
DIM1 * dof);
1479 ublas::shallow_array_adaptor<double>(DIM0 *
DIM1, data));
1482 template <
int DIM0,
int DIM1>
1485 return getVectorDiffN<DIM0, DIM1>(bAse, dof, gg);
1490 double *ptr = &*getN(base).data().begin();
1495 EntitiesFieldData::EntData::getFTensor0N() {
1496 return getFTensor0N(bAse);
1501 const int gg,
const int bb) {
1502 double *ptr = &getN(base)(gg, bb);
1507 EntitiesFieldData::EntData::getFTensor0N(
const int gg,
const int bb) {
1508 return getFTensor0N(bAse, gg, bb);
1511 template <
int Tensor_Dim>
auto EntitiesFieldData::EntData::getFTensor1N() {
1512 return getFTensor1N<Tensor_Dim>(bAse);
1515 template <
int Tensor_Dim>
1516 auto EntitiesFieldData::EntData::getFTensor1N(
const int gg,
const int bb) {
1517 return getFTensor1N<Tensor_Dim>(bAse, gg, bb);
1520 template <
int Tensor_Dim0,
int Tensor_Dim1>
1521 auto EntitiesFieldData::EntData::getFTensor2N() {
1522 return getFTensor2N<Tensor_Dim0, Tensor_Dim1>(bAse);
1525 template <
int Tensor_Dim0,
int Tensor_Dim1>
1526 auto EntitiesFieldData::EntData::getFTensor2N(
const int gg,
const int bb) {
1527 return getFTensor2N<Tensor_Dim0, Tensor_Dim1>(bAse, gg, bb);
1534 VectorInt &EntitiesFieldData::EntData::getBBNodeOrder() {
return bbNodeOrder; }
1536 MatrixInt &EntitiesFieldData::EntData::getBBAlphaIndices() {
1537 return *getBBAlphaIndicesSharedPtr(bbFieldName);
1546 boost::shared_ptr<MatrixDouble> &
1547 DerivedEntitiesFieldData::DerivedEntData::getDerivedNSharedPtr(
1552 boost::shared_ptr<MatrixDouble> &
1553 DerivedEntitiesFieldData::DerivedEntData::getDerivedDiffNSharedPtr(
1575 template <
typename T = EntityStorage>
1578 const double *ptr, InsertMode iora) {
1579 static_assert(!std::is_same<T, T>::value,
1580 "VecSetValues value for this data storage is not implemented");
1587 const double *ptr, InsertMode iora) {
1607 template <
typename T = EntityStorage>
1611 return VecSetValues<T>(V, data, &*vec.data().begin(), iora);
1630 template <
typename T = EntityStorage>
1634 const double *ptr, InsertMode iora) {
1635 static_assert(!std::is_same<T, T>::value,
1636 "MatSetValues value for this data storage is not implemented");
1656 template <
typename T = EntityStorage>
1661 return MatSetValues<T>(
M, row_data, col_data, &*mat.data().begin(), iora);
1668 const double *ptr, InsertMode iora) {
1685 const int gg,
const int bb);
1699 EntitiesFieldData::EntData::getFTensor1DiffN<3>(
1703 EntitiesFieldData::EntData::getFTensor1DiffN<3>();
1707 EntitiesFieldData::EntData::getFTensor1DiffN<2>(
1711 EntitiesFieldData::EntData::getFTensor1DiffN<2>();
1719 const int gg,
const int bb);
1733 EntitiesFieldData::EntData::getFTensor1FieldData<3>();
1737 EntitiesFieldData::EntData::getFTensor1FieldData<2>();
1741 EntitiesFieldData::EntData::getFTensor1FieldData<1>();
1745 EntitiesFieldData::EntData::getFTensor2FieldData<1, 1>();
1749 EntitiesFieldData::EntData::getFTensor2FieldData<1, 2>();
1753 EntitiesFieldData::EntData::getFTensor2FieldData<1, 3>();
1757 EntitiesFieldData::EntData::getFTensor2FieldData<2, 2>();
1761 EntitiesFieldData::EntData::getFTensor2FieldData<3, 3>();
1765 EntitiesFieldData::EntData::getFTensor2SymmetricFieldData<3>();
1769 EntitiesFieldData::EntData::getFTensor2SymmetricFieldData<2>();
1785 #endif //__ENTITIES_FIELD_DATA_HPP__