6 #ifndef __USER_DATA_OPERATORS_HPP__
7 #define __USER_DATA_OPERATORS_HPP__
22 template <
class T,
class A>
28 boost::shared_ptr<ublas::vector<T, A>> data_ptr,
29 const EntityType zero_type = MBVERTEX)
38 boost::shared_ptr<ublas::vector<T, A>> data_ptr,
57 boost::shared_ptr<ublas::vector<T, A>>
dataPtr;
67 template <
class T,
class A>
101 vec.resize(nb_gauss_pts,
false);
114 for (
int i = 0;
i != local_indices.size(); ++
i)
115 if (local_indices[
i] != -1)
123 const size_t nb_base_functions = data.
getN().size2();
126 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
129 for (; bb != nb_dofs; ++bb) {
130 values_at_gauss_pts += field_data * base_function;
135 for (; bb < nb_base_functions; ++bb)
137 ++values_at_gauss_pts;
155 template <PetscData::DataContext CTX>
160 const std::string
field_name, boost::shared_ptr<VectorDouble> data_ptr,
161 const EntityType zero_at_type = MBVERTEX)
177 vec.resize(nb_gauss_pts,
false);
182 const size_t nb_dofs = local_indices.size();
187 auto get_array = [&](
const auto ctx,
auto vec) {
193 <<
"In this case filed degrees of freedom are read from vector. "
194 "That usually happens when time solver is used, and acces to "
195 "first or second rates is needed. You probably not set ts_u, "
196 "ts_u_t, or ts_u_tt and associated data structure, i.e. "
197 "data_ctx to CTX_SET_X, CTX_SET_X_T, or CTX_SET_X_TT "
202 CHKERR VecGetArrayRead(vec, &array);
206 auto restore_array = [&](
auto vec) {
207 return VecRestoreArrayRead(vec, &array);
222 "That case is not implemented");
225 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
226 for (
int i = 0;
i != local_indices.size(); ++
i)
227 if (local_indices[
i] != -1)
228 dot_dofs_vector[
i] = array[local_indices[
i]];
230 dot_dofs_vector[
i] = 0;
244 "That case is not implemented");
247 const size_t nb_base_functions = data.
getN().size2();
251 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
253 for (; bb != nb_dofs; ++bb) {
254 values_at_gauss_pts += dot_dofs_vector[bb] * base_function;
259 for (; bb < nb_base_functions; ++bb)
261 ++values_at_gauss_pts;
292 template <
int Tensor_Dim,
class T,
class L,
class A>
298 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
299 const EntityType zero_type = MBVERTEX)
318 boost::shared_ptr<ublas::matrix<T, L, A>>
dataPtr;
322 template <
int Tensor_Dim,
class T,
class L,
class A>
328 "Not implemented for T = %s and dim = %d",
338 template <
int Tensor_Dim>
344 boost::shared_ptr<MatrixDouble> data_ptr,
345 const EntityType zero_type = MBVERTEX)
354 boost::shared_ptr<MatrixDouble> data_ptr,
356 const EntityType zero_type = MBVERTEX)
378 template <
int Tensor_Dim>
385 const size_t nb_gauss_pts = getGaussPts().size2();
386 auto &mat = *dataPtr;
387 if (
type == zeroType || mat.size2() != nb_gauss_pts) {
388 mat.resize(Tensor_Dim, nb_gauss_pts,
false);
395 if (dataVec.use_count()) {
396 dotVector.resize(nb_dofs,
false);
398 CHKERR VecGetArrayRead(dataVec, &array);
400 for (
int i = 0;
i != local_indices.size(); ++
i)
401 if (local_indices[
i] != -1)
402 dotVector[
i] = array[local_indices[
i]];
405 CHKERR VecRestoreArrayRead(dataVec, &array);
410 const size_t nb_base_functions = data.
getN().size2();
412 auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
414 const size_t size = nb_dofs / Tensor_Dim;
415 if (nb_dofs % Tensor_Dim) {
417 "Data inconsistency");
419 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
423 if (field_data.l2() != field_data.l2()) {
424 MOFEM_LOG(
"SELF", Sev::error) <<
"field data: " << field_data;
426 "Wrong number in coefficients");
431 for (; bb != size; ++bb) {
434 if (base_function != base_function) {
435 MOFEM_LOG(
"SELF", Sev::error) <<
"base finction: " << base_function;
437 "Wrong number number in base functions");
441 values_at_gauss_pts(
I) += field_data(
I) * base_function;
447 for (; bb < nb_base_functions; ++bb)
449 ++values_at_gauss_pts;
453 if (dataVec.use_count()) {
465 template <
int Tensor_Dim>
468 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
485 template <
int Tensor_Dim, CoordinateTypes COORDINATE_SYSTEM = CARTESIAN>
490 const std::string
field_name, boost::shared_ptr<VectorDouble> data_ptr,
491 const EntityType zero_type = MBVERTEX)
504 if constexpr (COORDINATE_SYSTEM ==
POLAR || COORDINATE_SYSTEM ==
SPHERICAL)
506 "%s coordiante not implemented",
512 vec.resize(nb_gauss_pts,
false);
520 const size_t nb_base_functions = data.
getN().size2();
523 const size_t size = nb_dofs / Tensor_Dim;
525 if (nb_dofs % Tensor_Dim) {
527 "Number of dofs should multiple of dimensions");
532 if constexpr (COORDINATE_SYSTEM ==
CARTESIAN) {
534 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
537 for (; bb != size; ++bb) {
538 values_at_gauss_pts += field_data(
I) * diff_base_function(
I);
540 ++diff_base_function;
544 for (; bb < nb_base_functions; ++bb)
545 ++diff_base_function;
546 ++values_at_gauss_pts;
556 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
559 for (; bb != size; ++bb) {
560 values_at_gauss_pts += field_data(
I) * diff_base_function(
I);
561 values_at_gauss_pts +=
562 base_function * (field_data(0) / t_coords(0));
565 ++diff_base_function;
569 for (; bb < nb_base_functions; ++bb) {
571 ++diff_base_function;
573 ++values_at_gauss_pts;
594 template <
int Tensor_Dim, PetscData::DataContext CTX>
599 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
600 const EntityType zero_at_type = MBVERTEX)
612 const size_t nb_dofs = local_indices.size();
617 mat.resize(Tensor_Dim, nb_gauss_pts,
false);
625 auto get_array = [&](
const auto ctx,
auto vec) {
631 <<
"In this case filed degrees of freedom are read from vector. "
632 "That usually happens when time solver is used, and acces to "
633 "first or second rates is needed. You probably not set ts_u, "
634 "ts_u_t, or ts_u_tt and associated data structure, i.e. "
635 "data_ctx to CTX_SET_X, CTX_SET_X_T, or CTX_SET_X_TT "
640 CHKERR VecGetArrayRead(vec, &array);
644 auto restore_array = [&](
auto vec) {
645 return VecRestoreArrayRead(vec, &array);
660 "That case is not implemented");
664 for (
int i = 0;
i != local_indices.size(); ++
i)
665 if (local_indices[
i] != -1)
682 "That case is not implemented");
685 const size_t nb_base_functions = data.
getN().size2();
687 auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
690 const size_t size = nb_dofs / Tensor_Dim;
691 if (nb_dofs % Tensor_Dim) {
694 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
695 auto field_data = getFTensor1FromArray<Tensor_Dim, Tensor_Dim>(
dotVector);
697 for (; bb != size; ++bb) {
698 values_at_gauss_pts(
I) += field_data(
I) * base_function;
704 for (; bb < nb_base_functions; ++bb)
706 ++values_at_gauss_pts;
722 template <
int Tensor_Dim>
732 template <
int Tensor_Dim>
746 template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
752 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
753 const EntityType zero_type = MBVERTEX)
765 boost::shared_ptr<ublas::matrix<T, L, A>>
dataPtr;
769 template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
775 "Not implemented for T = %s, dim0 = %d and dim1 = %d",
777 Tensor_Dim0, Tensor_Dim1);
781 template <
int Tensor_Dim0,
int Tensor_Dim1>
789 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
791 const EntityType zero_type = MBVERTEX)
802 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
816 boost::shared_ptr<ublas::matrix<double, ublas::row_major, DoubleAllocator>>
823 template <
int Tensor_Dim0,
int Tensor_Dim1>
831 const size_t nb_gauss_pts = data.
getN().size1();
832 if (
type == zeroType) {
833 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
839 if (dataVec.use_count()) {
840 dotVector.resize(nb_dofs,
false);
842 CHKERR VecGetArrayRead(dataVec, &array);
844 for (
int i = 0;
i != local_indices.size(); ++
i)
845 if (local_indices[
i] != -1)
846 dotVector[
i] = array[local_indices[
i]];
849 CHKERR VecRestoreArrayRead(dataVec, &array);
854 const size_t nb_base_functions = data.
getN().size2();
856 auto values_at_gauss_pts =
857 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
860 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
861 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
864 for (; bb != size; ++bb) {
865 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
869 for (; bb < nb_base_functions; ++bb)
871 ++values_at_gauss_pts;
874 if (dataVec.use_count()) {
886 template <
int Tensor_Dim0,
int Tensor_Dim1>
889 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
901 template <
int Tensor_Dim0,
int Tensor_Dim1>
906 boost::shared_ptr<MatrixDouble> data_ptr,
907 const EntityType zero_at_type = MBVERTEX)
922 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
926 const size_t nb_dofs = local_indices.size();
931 for (
size_t i = 0;
i != local_indices.size(); ++
i)
932 if (local_indices[
i] != -1)
938 const size_t nb_base_functions = data.
getN().size2();
941 auto values_at_gauss_pts =
942 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
945 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
946 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
948 getFTensor2FromPtr<Tensor_Dim0, Tensor_Dim1>(&*
dotVector.begin());
950 for (; bb != size; ++bb) {
951 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
955 for (; bb < nb_base_functions; ++bb)
957 ++values_at_gauss_pts;
977 template <
int Tensor_Dim>
982 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
983 const EntityType zero_type = MBEDGE,
const int zero_side = 0)
992 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
994 const int zero_side = 0)
1009 mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts,
false);
1018 const double *array;
1021 for (
int i = 0;
i != local_indices.size(); ++
i)
1022 if (local_indices[
i] != -1)
1030 const int nb_base_functions = data.
getN().size2();
1032 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1035 const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
1036 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1039 for (; bb != size; ++bb) {
1040 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
1044 for (; bb < nb_base_functions; ++bb)
1046 ++values_at_gauss_pts;
1071 template <
int Tensor_Dim>
1076 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1077 const EntityType zero_type = MBEDGE,
const int zero_side = 0)
1090 constexpr
auto symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
1092 mat.resize(symm_size, nb_gauss_pts,
false);
1096 const int nb_dofs = local_indices.size();
1104 <<
"In this case filed degrees of freedom are read from vector. "
1105 "That usually happens when time solver is used, and acces to "
1106 "first rates is needed. You probably not set "
1107 "ts_u_t and associated data structure data_ctx to CTX_SET_X_T "
1114 const double *array;
1116 for (
int i = 0;
i != local_indices.size(); ++
i)
1117 if (local_indices[
i] != -1)
1123 const int nb_base_functions = data.
getN().size2();
1126 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1129 const int size = nb_dofs / symm_size;
1130 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1131 auto field_data = getFTensorDotData<Tensor_Dim>();
1133 for (; bb != size; ++bb) {
1134 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
1138 for (; bb < nb_base_functions; ++bb)
1140 ++values_at_gauss_pts;
1153 static_assert(Dim || !Dim,
"not implemented");
1162 &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3], &dotVector[4],
1171 &dotVector[0], &dotVector[1], &dotVector[2]);
1185 template <
int Tensor_Dim,
class T,
class L,
class A>
1197 template <
int Tensor_Dim>
1201 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1223 template <
int Tensor_Dim>
1230 const size_t nb_gauss_pts = this->getGaussPts().size2();
1231 auto &mat = *this->dataPtr;
1232 if (
type == this->zeroType) {
1233 mat.resize(Tensor_Dim, nb_gauss_pts,
false);
1241 const int nb_base_functions = data.
getN().size2();
1243 auto gradients_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
1246 if (nb_dofs > nb_base_functions)
1248 "Number of base functions inconsistent with number of DOFs "
1250 nb_dofs, nb_base_functions);
1252 if (data.
getDiffN().size2() != nb_base_functions * Tensor_Dim)
1255 "Number of base functions inconsistent with number of derivatives "
1257 data.
getDiffN().size2(), nb_base_functions);
1259 if (data.
getDiffN().size1() != nb_gauss_pts)
1262 "Number of base functions inconsistent with number of integration "
1264 data.
getDiffN().size2(), nb_gauss_pts);
1269 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1272 for (; bb != nb_dofs; ++bb) {
1273 gradients_at_gauss_pts(
I) += field_data * diff_base_function(
I);
1275 ++diff_base_function;
1278 for (; bb < nb_base_functions; ++bb)
1279 ++diff_base_function;
1280 ++gradients_at_gauss_pts;
1293 template <
int Tensor_Dim>
1296 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1306 template <
int Tensor_Dim>
1309 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1326 template <
int Tensor_Dim>
1331 const size_t nb_gauss_pts = this->getGaussPts().size2();
1332 auto &mat = *this->dataPtr;
1333 if (
type == this->zeroType) {
1334 mat.resize(Tensor_Dim * Tensor_Dim, nb_gauss_pts,
false);
1342 const int nb_base_functions = data.
getN().size2();
1344 auto &hessian_base = data.
getN(BaseDerivatives::SecondDerivative);
1346 if (hessian_base.size1() != nb_gauss_pts) {
1348 "Wrong number of integration pts (%d != %d)",
1349 hessian_base.size1(), nb_gauss_pts);
1351 if (hessian_base.size2() != nb_base_functions * Tensor_Dim * Tensor_Dim) {
1353 "Wrong number of base functions (%d != %d)",
1354 hessian_base.size2() / (Tensor_Dim * Tensor_Dim),
1357 if (hessian_base.size2() < nb_dofs * Tensor_Dim * Tensor_Dim) {
1359 "Wrong number of base functions (%d < %d)",
1360 hessian_base.size2(), nb_dofs * Tensor_Dim * Tensor_Dim);
1364 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim, Tensor_Dim>(
1365 &*hessian_base.data().begin());
1367 auto hessian_at_gauss_pts =
1368 getFTensor2FromMat<Tensor_Dim, Tensor_Dim>(mat);
1372 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1375 for (; bb != nb_dofs; ++bb) {
1376 hessian_at_gauss_pts(
I,
J) +=
1377 field_data * t_diff2_base_function(
I,
J);
1379 ++t_diff2_base_function;
1382 for (; bb < nb_base_functions; ++bb) {
1383 ++t_diff2_base_function;
1386 ++hessian_at_gauss_pts;
1405 template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
1414 template <
int Tensor_Dim0,
int Tensor_Dim1>
1418 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1440 template <
int Tensor_Dim0,
int Tensor_Dim1>
1448 "Data pointer not allocated");
1450 const size_t nb_gauss_pts = this->getGaussPts().size2();
1451 auto &mat = *this->dataPtr;
1452 if (
type == this->zeroType) {
1453 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
1462 if (this->dataVec.use_count()) {
1463 this->dotVector.resize(nb_dofs,
false);
1464 const double *array;
1465 CHKERR VecGetArrayRead(this->dataVec, &array);
1467 for (
int i = 0;
i != local_indices.size(); ++
i)
1468 if (local_indices[
i] != -1)
1469 this->dotVector[
i] = array[local_indices[
i]];
1471 this->dotVector[
i] = 0;
1472 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1476 const int nb_base_functions = data.
getN().size2();
1478 auto gradients_at_gauss_pts =
1479 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1482 int size = nb_dofs / Tensor_Dim0;
1483 if (nb_dofs % Tensor_Dim0) {
1485 "Data inconsistency");
1487 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1491 if (field_data.l2() != field_data.l2()) {
1492 MOFEM_LOG(
"SELF", Sev::error) <<
"field data " << field_data;
1494 "Wrong number in coefficients");
1499 for (; bb < size; ++bb) {
1502 if (diff_base_function.l2() != diff_base_function.l2()) {
1504 <<
"diff_base_function: " << diff_base_function;
1506 "Wrong number number in base functions");
1511 gradients_at_gauss_pts(
I,
J) += field_data(
I) * diff_base_function(
J);
1513 ++diff_base_function;
1517 for (; bb != nb_base_functions; ++bb)
1518 ++diff_base_function;
1519 ++gradients_at_gauss_pts;
1522 if (this->dataVec.use_count()) {
1535 template <
int Tensor_Dim0,
int Tensor_Dim1>
1538 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1550 template <
int Tensor_Dim0,
int Tensor_Dim1>
1555 boost::shared_ptr<MatrixDouble> data_ptr,
1556 const EntityType zero_at_type = MBVERTEX)
1569 const int nb_dofs = local_indices.size();
1570 const int nb_gauss_pts = this->
getGaussPts().size2();
1574 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
1581 const double *array;
1583 for (
int i = 0;
i != local_indices.size(); ++
i)
1584 if (local_indices[
i] != -1)
1590 const int nb_base_functions = data.
getN().size2();
1592 auto gradients_at_gauss_pts =
1593 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1596 int size = nb_dofs / Tensor_Dim0;
1597 if (nb_dofs % Tensor_Dim0) {
1601 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1602 auto field_data = getFTensor1FromPtr<Tensor_Dim0>(&*
dotVector.begin());
1604 for (; bb < size; ++bb) {
1605 gradients_at_gauss_pts(
I,
J) += field_data(
I) * diff_base_function(
J);
1607 ++diff_base_function;
1611 for (; bb != nb_base_functions; ++bb)
1612 ++diff_base_function;
1613 ++gradients_at_gauss_pts;
1630 template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
1636 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1637 const EntityType zero_type = MBVERTEX)
1643 template <
int Tensor_Dim0,
int Tensor_Dim1>
1647 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1650 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1651 const EntityType zero_type = MBVERTEX)
1673 template <
int Tensor_Dim0,
int Tensor_Dim1>
1681 "Data pointer not allocated");
1683 const size_t nb_gauss_pts = this->getGaussPts().size2();
1684 constexpr
size_t msize = (Tensor_Dim0 * (Tensor_Dim0 + 1)) / 2;
1685 auto &mat = *this->dataPtr;
1686 if (
type == this->zeroType) {
1687 mat.resize(msize * Tensor_Dim1, nb_gauss_pts,
false);
1696 const int nb_base_functions = data.
getN().size2();
1698 auto gradients_at_gauss_pts =
1699 getFTensor3DgFromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1703 int size = nb_dofs / msize;
1704 if (nb_dofs % msize) {
1706 "Data inconsistency");
1708 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1711 for (; bb < size; ++bb) {
1712 gradients_at_gauss_pts(
I,
J,
K) +=
1713 field_data(
I,
J) * diff_base_function(
K);
1715 ++diff_base_function;
1719 for (; bb != nb_base_functions; ++bb)
1720 ++diff_base_function;
1721 ++gradients_at_gauss_pts;
1733 template <
int Tensor_Dim0,
int Tensor_Dim1>
1736 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1739 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1740 const EntityType zero_type = MBVERTEX)
1746 template <
int Tensor_Dim0,
int Tensor_Dim1>
1749 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1771 template <
int Tensor_Dim0,
int Tensor_Dim1>
1777 "Data pointer not allocated");
1779 const size_t nb_gauss_pts = this->getGaussPts().size2();
1780 auto &mat = *this->dataPtr;
1781 if (
type == this->zeroType) {
1782 mat.resize(Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim1, nb_gauss_pts,
false);
1791 if (this->dataVec.use_count()) {
1792 this->dotVector.resize(nb_dofs,
false);
1793 const double *array;
1794 CHKERR VecGetArrayRead(this->dataVec, &array);
1796 for (
int i = 0;
i != local_indices.size(); ++
i)
1797 if (local_indices[
i] != -1)
1798 this->dotVector[
i] = array[local_indices[
i]];
1800 this->dotVector[
i] = 0;
1801 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1805 const int nb_base_functions = data.
getN().size2();
1807 auto &hessian_base = data.
getN(BaseDerivatives::SecondDerivative);
1809 if (hessian_base.size1() != nb_gauss_pts) {
1811 "Wrong number of integration pts (%d != %d)",
1812 hessian_base.size1(), nb_gauss_pts);
1814 if (hessian_base.size2() !=
1815 nb_base_functions * Tensor_Dim1 * Tensor_Dim1) {
1817 "Wrong number of base functions (%d != %d)",
1818 hessian_base.size2() / (Tensor_Dim1 * Tensor_Dim1),
1821 if (hessian_base.size2() < nb_dofs * Tensor_Dim1 * Tensor_Dim1) {
1823 "Wrong number of base functions (%d < %d)",
1824 hessian_base.size2(), nb_dofs * Tensor_Dim1 * Tensor_Dim1);
1828 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim1, Tensor_Dim1>(
1829 &*hessian_base.data().begin());
1831 auto t_hessian_at_gauss_pts =
1832 getFTensor3FromMat<Tensor_Dim0, Tensor_Dim1, Tensor_Dim1>(mat);
1838 int size = nb_dofs / Tensor_Dim0;
1840 if (nb_dofs % Tensor_Dim0) {
1842 "Data inconsistency");
1846 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1849 for (; bb < size; ++bb) {
1850 t_hessian_at_gauss_pts(
I,
J,
K) +=
1851 field_data(
I) * t_diff2_base_function(
J,
K);
1853 ++t_diff2_base_function;
1857 for (; bb != nb_base_functions; ++bb)
1858 ++t_diff2_base_function;
1859 ++t_hessian_at_gauss_pts;
1862 if (this->dataVec.use_count()) {
1884 template <
int DIM_01,
int DIM_23,
int S = 0>
1895 boost::shared_ptr<MatrixDouble> in_mat,
1896 boost::shared_ptr<MatrixDouble> out_mat,
1897 boost::shared_ptr<MatrixDouble> d_mat)
1910 boost::shared_ptr<MatrixDouble> out_mat,
1911 boost::shared_ptr<MatrixDouble> d_mat)
1924 const size_t nb_gauss_pts =
getGaussPts().size2();
1925 auto t_D = getFTensor4DdgFromMat<DIM_01, DIM_23, S>(*(
dMat));
1926 auto t_in = getFTensor2SymmetricFromMat<DIM_01>(*(
inMat));
1927 outMat->resize((DIM_23 * (DIM_23 + 1)) / 2, nb_gauss_pts,
false);
1928 auto t_out = getFTensor2SymmetricFromMat<DIM_23>(*(
outMat));
1929 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1930 t_out(
i,
j) = t_D(
i,
j,
k,
l) * t_in(
k,
l);
1945 boost::shared_ptr<MatrixDouble>
dMat;
1958 boost::shared_ptr<MatrixDouble> in_mat,
1959 boost::shared_ptr<MatrixDouble> out_mat)
1970 boost::shared_ptr<MatrixDouble> out_mat)
1981 const size_t nb_gauss_pts =
getGaussPts().size2();
1982 auto t_in = getFTensor2FromMat<DIM, DIM>(*(
inMat));
1983 outMat->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts,
false);
1984 auto t_out = getFTensor2SymmetricFromMat<DIM>(*(
outMat));
1985 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1986 t_out(
i,
j) = (t_in(
i,
j) || t_in(
j,
i)) / 2;
2006 boost::shared_ptr<MatrixDouble> in_mat,
2007 boost::shared_ptr<MatrixDouble> out_mat)
2040 template <
int Base_Dim,
int Field_Dim,
class T,
class L,
class A>
2046 template <
int Field_Dim>
2052 boost::shared_ptr<MatrixDouble> data_ptr,
2053 const EntityType zero_type = MBEDGE,
2054 const int zero_side = 0)
2057 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2078 template <
int Field_Dim>
2084 const size_t nb_integration_points = this->getGaussPts().size2();
2085 if (
type == zeroType && side == zeroSide) {
2086 dataPtr->resize(Field_Dim, nb_integration_points,
false);
2092 const size_t nb_base_functions = data.
getN().size2() / 3;
2095 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2096 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2099 for (; bb != nb_dofs; ++bb) {
2100 t_data(
i) += t_n_hdiv(
i) * t_dof;
2104 for (; bb != nb_base_functions; ++bb)
2115 template <
int Base_Dim,
int Field_Dim = Base_Dim>
2118 Base_Dim, Field_Dim, double, ublas::row_major, DoubleAllocator> {
2129 template <
int Base_Dim,
int Field_Dim = Base_Dim>
2132 template <
int Field_Dim>
2137 boost::shared_ptr<MatrixDouble> data_ptr,
2138 const EntityType zero_type = MBEDGE,
2139 const int zero_side = 0)
2142 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2163 template <
int Field_Dim>
2168 const size_t nb_integration_points = this->getGaussPts().size2();
2169 if (
type == zeroType && side == zeroSide) {
2170 dataPtr->resize(Field_Dim, nb_integration_points,
false);
2175 const size_t nb_dofs = local_indices.size();
2178 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2179 const double *array;
2180 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2181 for (
size_t i = 0;
i != nb_dofs; ++
i)
2182 if (local_indices[
i] != -1)
2183 dot_dofs_vector[
i] = array[local_indices[
i]];
2185 dot_dofs_vector[
i] = 0;
2186 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2188 const size_t nb_base_functions = data.
getN().size2() / 3;
2191 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2192 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2194 for (; bb != nb_dofs; ++bb) {
2195 t_data(
i) += t_n_hdiv(
i) * dot_dofs_vector[bb];
2198 for (; bb != nb_base_functions; ++bb)
2214 template <
int BASE_DIM,
int SPACE_DIM>
2219 boost::shared_ptr<VectorDouble> data_ptr,
2220 const EntityType zero_type = MBEDGE,
2221 const int zero_side = 0)
2232 const size_t nb_integration_points =
getGaussPts().size2();
2234 dataPtr->resize(nb_integration_points,
false);
2240 const size_t nb_base_functions = data.
getN().size2() /
BASE_DIM;
2245 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2248 for (; bb != nb_dofs; ++bb) {
2249 t_data += t_dof * t_n_diff_hdiv(
j,
j);
2253 for (; bb != nb_base_functions; ++bb)
2273 template <
int BASE_DIM,
int SPACE_DIM>
2278 boost::shared_ptr<MatrixDouble> data_ptr,
2279 const EntityType zero_type = MBEDGE,
2280 const int zero_side = 0)
2291 const size_t nb_integration_points =
getGaussPts().size2();
2299 const size_t nb_base_functions = data.
getN().size2() /
BASE_DIM;
2303 auto t_data = getFTensor2FromMat<BASE_DIM, SPACE_DIM>(*
dataPtr);
2304 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2307 for (; bb != nb_dofs; ++bb) {
2308 t_data(
i,
j) += t_dof * t_base_diff(
i,
j);
2312 for (; bb != nb_base_functions; ++bb)
2332 template <
int BASE_DIM,
int SPACE_DIM>
2337 boost::shared_ptr<MatrixDouble> data_ptr,
2338 const EntityType zero_type = MBEDGE,
2339 const int zero_side = 0)
2350 const size_t nb_integration_points =
getGaussPts().size2();
2360 const int nb_base_functions = data.
getN().size2() /
BASE_DIM;
2363 auto &hessian_base = data.
getN(BaseDerivatives::SecondDerivative);
2364 if (hessian_base.size1() != nb_integration_points) {
2366 "Wrong number of integration pts (%d != %d)",
2367 hessian_base.size1(), nb_integration_points);
2369 if (hessian_base.size2() !=
2372 "Wrong number of base functions (%d != %d)",
2378 "Wrong number of base functions (%d < %d)", hessian_base.size2(),
2389 auto t_data = getFTensor3FromMat<BASE_DIM, SPACE_DIM, SPACE_DIM>(*
dataPtr);
2390 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2393 for (; bb != nb_dofs; ++bb) {
2394 t_data(
i,
j,
k) += t_dof * t_base_diff2(
i,
j,
k);
2399 for (; bb != nb_base_functions; ++bb)
2418 template <
int Tensor_Dim1,
int Tensor_Dim2>
2423 boost::shared_ptr<VectorDouble> data_ptr,
2424 const EntityType zero_type = MBEDGE,
2425 const int zero_side = 0)
2436 const size_t nb_integration_points =
getGaussPts().size2();
2438 dataPtr->resize(nb_integration_points,
false);
2443 const int nb_dofs = local_indices.size();
2446 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2447 const double *array;
2449 for (
size_t i = 0;
i != local_indices.size(); ++
i)
2450 if (local_indices[
i] != -1)
2451 dot_dofs_vector[
i] = array[local_indices[
i]];
2453 dot_dofs_vector[
i] = 0;
2456 const size_t nb_base_functions = data.
getN().size2() / Tensor_Dim1;
2460 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2462 for (; bb != nb_dofs; ++bb) {
2464 for (
auto ii = 0; ii != Tensor_Dim2; ++ii)
2465 div += t_n_diff_hdiv(ii, ii);
2466 t_data += dot_dofs_vector[bb] * div;
2469 for (; bb != nb_base_functions; ++bb)
2491 template <
int Base_Dim,
int Space_Dim>
2506 boost::shared_ptr<MatrixDouble> data_ptr,
2507 const EntityType zero_type = MBEDGE,
2508 const int zero_side = 0);
2530 boost::shared_ptr<MatrixDouble> data_ptr,
2531 const EntityType zero_type = MBVERTEX,
2532 const int zero_side = 0);
2555 boost::shared_ptr<MatrixDouble> data_ptr,
2556 const EntityType zero_type = MBVERTEX,
2557 const int zero_side = 0);
2575 template <
int Tensor_Dim0,
int Tensor_Dim1>
2580 boost::shared_ptr<MatrixDouble> data_ptr,
2581 const EntityType zero_type = MBEDGE,
2582 const int zero_side = 0)
2593 const size_t nb_integration_points =
getGaussPts().size2();
2595 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points,
false);
2600 const size_t nb_base_functions = data.
getN().size2() / 3;
2604 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*
dataPtr);
2605 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2608 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2609 t_data(
i,
j) += t_dof(
i) * t_n_hvec(
j);
2613 for (; bb < nb_base_functions; ++bb)
2634 template <
int Tensor_Dim0,
int Tensor_Dim1>
2639 boost::shared_ptr<MatrixDouble> data_ptr,
2640 const EntityType zero_type = MBEDGE,
2641 const int zero_side = 0)
2652 const size_t nb_integration_points =
getGaussPts().size2();
2654 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points,
false);
2660 const size_t nb_base_functions =
2661 data.
getN().size2() / (Tensor_Dim0 * Tensor_Dim1);
2664 auto t_n_hten = data.
getFTensor2N<Tensor_Dim0, Tensor_Dim1>();
2665 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*
dataPtr);
2666 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2669 for (; bb != nb_dofs; ++bb) {
2670 t_data(
i,
j) += t_dof * t_n_hten(
i,
j);
2674 for (; bb < nb_base_functions; ++bb)
2694 template <
int Tensor_Dim0,
int Tensor_Dim1,
2700 boost::shared_ptr<MatrixDouble> data_ptr,
2701 const EntityType zero_type = MBEDGE,
2702 const int zero_side = 0)
2713 const size_t nb_integration_points =
getGaussPts().size2();
2715 dataPtr->resize(Tensor_Dim0, nb_integration_points,
false);
2720 const size_t nb_base_functions = data.
getN().size2() / 3;
2724 auto t_data = getFTensor1FromMat<Tensor_Dim0>(*
dataPtr);
2727 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2730 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2731 double div = t_n_diff_hvec(
j,
j);
2732 t_data(
i) += t_dof(
i) * div;
2734 t_data(
i) += t_base(0) * (t_dof(
i) / t_coords(0));
2740 for (; bb < nb_base_functions; ++bb) {
2763 template <
int Tensor_Dim,
typename OpBase>
2767 boost::shared_ptr<MatrixDouble> data_ptr,
2768 const EntityType zero_type = MBEDGE,
2769 const int zero_side = 0)
2779 const size_t nb_integration_points = OpBase::getGaussPts().size2();
2781 dataPtr->resize(Tensor_Dim, nb_integration_points,
false);
2786 auto t_normal = OpBase::getFTensor1NormalsAtGaussPts();
2787 const size_t nb_base_functions = data.
getN().size2() / 3;
2789 auto t_data = getFTensor1FromMat<Tensor_Dim>(*
dataPtr);
2790 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2792 t_normalized_normal(
j) = t_normal(
j);
2796 for (; bb != nb_dofs / Tensor_Dim; ++bb) {
2797 t_data(
i) += t_dof(
i) * (t_base(
j) * t_normalized_normal(
j));
2801 for (; bb < nb_base_functions; ++bb) {
2843 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
2846 template <
int D1,
int D2,
int J1,
int J2>
2850 static_assert(D2 == J2,
"Dimension of jacobian and dimension of <out> "
2851 "directive does not match");
2853 size_t nb_functions = diff_n.size2() / D1;
2855 size_t nb_gauss_pts = diff_n.size1();
2860 "Wrong number of Gauss Pts");
2861 if (diff_n.size2() % D1)
2863 "Number of direvatives of base functions and D1 dimension does "
2867 diffNinvJac.resize(diff_n.size1(), D2 * nb_functions,
false);
2871 auto t_diff_n = getFTensor1FromPtr<D2>(&*
diffNinvJac.data().begin());
2872 auto t_diff_n_ref = getFTensor1FromPtr<D1>(&*diff_n.data().begin());
2873 auto t_inv_jac = getFTensor2FromMat<J1, J2>(*
invJacPtr);
2874 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
2875 for (
size_t dd = 0;
dd != nb_functions; ++
dd) {
2876 t_diff_n(
i) = t_inv_jac(
k,
i) * t_diff_n_ref(
k);
2921 template <
int DERIVARIVE = 1>
2928 template <
int DERIVARIVE = 1>
2935 template <
int DERIVARIVE = 1>
2939 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2943 template <
int DERIVARIVE = 1>
2947 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2965 invJacPtr(inv_jac_ptr) {}
3019 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
3056 boost::shared_ptr<MatrixDouble> jac_ptr)
3257 boost::shared_ptr<VectorDouble> det_ptr,
3258 boost::shared_ptr<MatrixDouble> out_ptr)
3279 "Pointer for inPtr matrix not allocated");
3282 "Pointer for detPtr matrix not allocated");
3284 const auto nb_rows = inPtr->size1();
3285 const auto nb_integration_pts = inPtr->size2();
3289 detPtr->resize(nb_integration_pts,
false);
3290 auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
3292 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
3301 outPtr->resize(nb_rows, nb_integration_pts,
false);
3302 auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
3303 auto t_out = getFTensor2FromMat<DIM, DIM>(*outPtr);
3305 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
3328 boost::shared_ptr<VectorDouble> out_ptr)
3350 "Pointer for inPtr matrix not allocated");
3352 const auto nb_integration_pts = inPtr->size2();
3355 outPtr->resize(nb_integration_pts,
false);
3356 auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
3359 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
3382 boost::shared_ptr<VectorDouble> out_ptr)
3404 "Pointer for inPtr matrix not allocated");
3406 const auto nb_integration_pts = inPtr->size2();
3409 outPtr->resize(nb_integration_pts,
false);
3410 auto t_in = getFTensor2SymmetricFromMat<DIM>(*inPtr);
3413 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
3426 #endif // __USER_DATA_OPERATORS_HPP__