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 "Nb. of DOFs is inconsistent with Tensor_Dim");
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) {
435 if (base_function != base_function) {
436 MOFEM_LOG(
"SELF", Sev::error) <<
"base finction: " << base_function;
438 "Wrong number number in base functions");
443 values_at_gauss_pts(
I) += field_data(
I) * base_function;
449 for (; bb < nb_base_functions; ++bb)
451 ++values_at_gauss_pts;
455 if (dataVec.use_count()) {
467 template <
int Tensor_Dim>
470 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
487 template <
int Tensor_Dim, CoordinateTypes COORDINATE_SYSTEM = CARTESIAN>
492 const std::string
field_name, boost::shared_ptr<VectorDouble> data_ptr,
493 const EntityType zero_type = MBVERTEX)
506 if constexpr (COORDINATE_SYSTEM ==
POLAR || COORDINATE_SYSTEM ==
SPHERICAL)
508 "%s coordiante not implemented",
514 vec.resize(nb_gauss_pts,
false);
522 const size_t nb_base_functions = data.
getN().size2();
525 const size_t size = nb_dofs / Tensor_Dim;
527 if (nb_dofs % Tensor_Dim) {
529 "Number of dofs should multiple of dimensions");
534 if constexpr (COORDINATE_SYSTEM ==
CARTESIAN) {
536 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
539 for (; bb != size; ++bb) {
540 values_at_gauss_pts += field_data(
I) * diff_base_function(
I);
542 ++diff_base_function;
546 for (; bb < nb_base_functions; ++bb)
547 ++diff_base_function;
548 ++values_at_gauss_pts;
558 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
561 for (; bb != size; ++bb) {
562 values_at_gauss_pts += field_data(
I) * diff_base_function(
I);
563 values_at_gauss_pts +=
564 base_function * (field_data(0) / t_coords(0));
567 ++diff_base_function;
571 for (; bb < nb_base_functions; ++bb) {
573 ++diff_base_function;
575 ++values_at_gauss_pts;
596 template <
int Tensor_Dim, PetscData::DataContext CTX>
601 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
602 const EntityType zero_at_type = MBVERTEX,
bool throw_error =
true)
615 const size_t nb_dofs = local_indices.size();
620 mat.resize(Tensor_Dim, nb_gauss_pts,
false);
634 auto get_array = [&](
const auto ctx,
auto vec) {
640 <<
"In this case filed degrees of freedom are read from vector. "
641 "That usually happens when time solver is used, and access to "
642 "first or second rates is needed. You probably not set ts_u, "
643 "ts_u_t, or ts_u_tt and associated data structure, i.e. "
644 "data_ctx to CTX_SET_X, CTX_SET_DX, CTX_SET_X_T, or "
645 "CTX_SET_X_TT respectively";
649 CHKERR VecGetArrayRead(vec, &array);
653 auto restore_array = [&](
auto vec) {
654 return VecRestoreArrayRead(vec, &array);
672 "That case is not implemented");
676 for (
int i = 0;
i != local_indices.size(); ++
i)
677 if (local_indices[
i] != -1)
697 "That case is not implemented");
700 const size_t nb_base_functions = data.
getN().size2();
702 auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
705 const size_t size = nb_dofs / Tensor_Dim;
706 if (nb_dofs % Tensor_Dim) {
709 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
710 auto field_data = getFTensor1FromArray<Tensor_Dim, Tensor_Dim>(
dotVector);
712 for (; bb != size; ++bb) {
713 values_at_gauss_pts(
I) += field_data(
I) * base_function;
719 for (; bb < nb_base_functions; ++bb)
721 ++values_at_gauss_pts;
738 template <
int Tensor_Dim>
748 template <
int Tensor_Dim>
758 template <
int Tensor_Dim>
772 template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
778 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
779 const EntityType zero_type = MBVERTEX)
791 boost::shared_ptr<ublas::matrix<T, L, A>>
dataPtr;
795 template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
801 "Not implemented for T = %s, dim0 = %d and dim1 = %d",
803 Tensor_Dim0, Tensor_Dim1);
807 template <
int Tensor_Dim0,
int Tensor_Dim1>
815 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
817 const EntityType zero_type = MBVERTEX)
828 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
842 boost::shared_ptr<ublas::matrix<double, ublas::row_major, DoubleAllocator>>
849 template <
int Tensor_Dim0,
int Tensor_Dim1>
857 const size_t nb_gauss_pts = data.
getN().size1();
858 if (
type == zeroType) {
859 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
865 if (dataVec.use_count()) {
866 dotVector.resize(nb_dofs,
false);
868 CHKERR VecGetArrayRead(dataVec, &array);
870 for (
int i = 0;
i != local_indices.size(); ++
i)
871 if (local_indices[
i] != -1)
872 dotVector[
i] = array[local_indices[
i]];
875 CHKERR VecRestoreArrayRead(dataVec, &array);
880 const size_t nb_base_functions = data.
getN().size2();
882 auto values_at_gauss_pts =
883 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
886 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
887 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
890 for (; bb != size; ++bb) {
891 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
895 for (; bb < nb_base_functions; ++bb)
897 ++values_at_gauss_pts;
900 if (dataVec.use_count()) {
912 template <
int Tensor_Dim0,
int Tensor_Dim1>
915 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
927 template <
int Tensor_Dim0,
int Tensor_Dim1>
932 boost::shared_ptr<MatrixDouble> data_ptr,
933 const EntityType zero_at_type = MBVERTEX)
948 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
952 const size_t nb_dofs = local_indices.size();
957 for (
size_t i = 0;
i != local_indices.size(); ++
i)
958 if (local_indices[
i] != -1)
964 const size_t nb_base_functions = data.
getN().size2();
967 auto values_at_gauss_pts =
968 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
971 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
972 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
974 getFTensor2FromPtr<Tensor_Dim0, Tensor_Dim1>(&*
dotVector.begin());
976 for (; bb != size; ++bb) {
977 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
981 for (; bb < nb_base_functions; ++bb)
983 ++values_at_gauss_pts;
1002 template <
int Tensor_Dim>
1007 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1008 const EntityType zero_type = MBEDGE,
const int zero_side = 0)
1017 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1019 const int zero_side = 0)
1034 mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts,
false);
1043 const double *array;
1046 for (
int i = 0;
i != local_indices.size(); ++
i)
1047 if (local_indices[
i] != -1)
1055 const int nb_base_functions = data.
getN().size2();
1057 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1060 const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
1061 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1064 for (; bb != size; ++bb) {
1065 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
1069 for (; bb < nb_base_functions; ++bb)
1071 ++values_at_gauss_pts;
1096 template <
int Tensor_Dim>
1101 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1102 const EntityType zero_type = MBEDGE,
const int zero_side = 0)
1115 constexpr
auto symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
1117 mat.resize(symm_size, nb_gauss_pts,
false);
1121 const int nb_dofs = local_indices.size();
1129 <<
"In this case filed degrees of freedom are read from vector. "
1130 "That usually happens when time solver is used, and acces to "
1131 "first rates is needed. You probably not set "
1132 "ts_u_t and associated data structure data_ctx to CTX_SET_X_T "
1139 const double *array;
1141 for (
int i = 0;
i != local_indices.size(); ++
i)
1142 if (local_indices[
i] != -1)
1148 const int nb_base_functions = data.
getN().size2();
1151 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1154 const int size = nb_dofs / symm_size;
1155 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1156 auto field_data = getFTensorDotData<Tensor_Dim>();
1158 for (; bb != size; ++bb) {
1159 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
1163 for (; bb < nb_base_functions; ++bb)
1165 ++values_at_gauss_pts;
1178 static_assert(Dim || !Dim,
"not implemented");
1187 &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3], &dotVector[4],
1196 &dotVector[0], &dotVector[1], &dotVector[2]);
1210 template <
int Tensor_Dim,
class T,
class L,
class A>
1222 template <
int Tensor_Dim>
1226 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1248 template <
int Tensor_Dim>
1255 const size_t nb_gauss_pts = this->getGaussPts().size2();
1256 auto &mat = *this->dataPtr;
1257 if (
type == this->zeroType) {
1258 mat.resize(Tensor_Dim, nb_gauss_pts,
false);
1266 const int nb_base_functions = data.
getN().size2();
1268 auto gradients_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
1271 if (nb_dofs > nb_base_functions)
1273 "Number of base functions inconsistent with number of DOFs "
1275 nb_dofs, nb_base_functions);
1277 if (data.
getDiffN().size2() != nb_base_functions * Tensor_Dim)
1280 "Number of base functions inconsistent with number of derivatives "
1282 data.
getDiffN().size2(), nb_base_functions);
1284 if (data.
getDiffN().size1() != nb_gauss_pts)
1287 "Number of base functions inconsistent with number of integration "
1289 data.
getDiffN().size2(), nb_gauss_pts);
1294 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1297 for (; bb != nb_dofs; ++bb) {
1298 gradients_at_gauss_pts(
I) += field_data * diff_base_function(
I);
1300 ++diff_base_function;
1303 for (; bb < nb_base_functions; ++bb)
1304 ++diff_base_function;
1305 ++gradients_at_gauss_pts;
1318 template <
int Tensor_Dim>
1321 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1331 template <
int Tensor_Dim>
1334 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1351 template <
int Tensor_Dim>
1356 const size_t nb_gauss_pts = this->getGaussPts().size2();
1357 auto &mat = *this->dataPtr;
1358 if (
type == this->zeroType) {
1359 mat.resize(Tensor_Dim * Tensor_Dim, nb_gauss_pts,
false);
1367 const int nb_base_functions = data.
getN().size2();
1369 auto &hessian_base = data.
getN(BaseDerivatives::SecondDerivative);
1371 if (hessian_base.size1() != nb_gauss_pts) {
1373 "Wrong number of integration pts (%d != %d)",
1374 hessian_base.size1(), nb_gauss_pts);
1376 if (hessian_base.size2() != nb_base_functions * Tensor_Dim * Tensor_Dim) {
1378 "Wrong number of base functions (%d != %d)",
1379 hessian_base.size2() / (Tensor_Dim * Tensor_Dim),
1382 if (hessian_base.size2() < nb_dofs * Tensor_Dim * Tensor_Dim) {
1384 "Wrong number of base functions (%d < %d)",
1385 hessian_base.size2(), nb_dofs * Tensor_Dim * Tensor_Dim);
1389 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim, Tensor_Dim>(
1390 &*hessian_base.data().begin());
1392 auto hessian_at_gauss_pts =
1393 getFTensor2FromMat<Tensor_Dim, Tensor_Dim>(mat);
1397 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1400 for (; bb != nb_dofs; ++bb) {
1401 hessian_at_gauss_pts(
I,
J) +=
1402 field_data * t_diff2_base_function(
I,
J);
1404 ++t_diff2_base_function;
1407 for (; bb < nb_base_functions; ++bb) {
1408 ++t_diff2_base_function;
1411 ++hessian_at_gauss_pts;
1430 template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
1439 template <
int Tensor_Dim0,
int Tensor_Dim1>
1443 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1465 template <
int Tensor_Dim0,
int Tensor_Dim1>
1473 "Data pointer not allocated");
1475 const size_t nb_gauss_pts = this->getGaussPts().size2();
1476 auto &mat = *this->dataPtr;
1477 if (
type == this->zeroType) {
1478 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
1487 if (this->dataVec.use_count()) {
1488 this->dotVector.resize(nb_dofs,
false);
1489 const double *array;
1490 CHKERR VecGetArrayRead(this->dataVec, &array);
1492 for (
int i = 0;
i != local_indices.size(); ++
i)
1493 if (local_indices[
i] != -1)
1494 this->dotVector[
i] = array[local_indices[
i]];
1496 this->dotVector[
i] = 0;
1497 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1501 const int nb_base_functions = data.
getN().size2();
1503 auto gradients_at_gauss_pts =
1504 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1507 int size = nb_dofs / Tensor_Dim0;
1508 if (nb_dofs % Tensor_Dim0) {
1510 "Data inconsistency");
1512 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1516 if (field_data.l2() != field_data.l2()) {
1517 MOFEM_LOG(
"SELF", Sev::error) <<
"field data " << field_data;
1519 "Wrong number in coefficients");
1524 for (; bb < size; ++bb) {
1527 if (diff_base_function.l2() != diff_base_function.l2()) {
1529 <<
"diff_base_function: " << diff_base_function;
1531 "Wrong number number in base functions");
1536 gradients_at_gauss_pts(
I,
J) += field_data(
I) * diff_base_function(
J);
1538 ++diff_base_function;
1542 for (; bb != nb_base_functions; ++bb)
1543 ++diff_base_function;
1544 ++gradients_at_gauss_pts;
1547 if (this->dataVec.use_count()) {
1560 template <
int Tensor_Dim0,
int Tensor_Dim1>
1563 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1575 template <
int Tensor_Dim0,
int Tensor_Dim1>
1580 boost::shared_ptr<MatrixDouble> data_ptr,
1581 const EntityType zero_at_type = MBVERTEX)
1594 const int nb_dofs = local_indices.size();
1595 const int nb_gauss_pts = this->
getGaussPts().size2();
1599 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
1606 const double *array;
1608 for (
int i = 0;
i != local_indices.size(); ++
i)
1609 if (local_indices[
i] != -1)
1615 const int nb_base_functions = data.
getN().size2();
1617 auto gradients_at_gauss_pts =
1618 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1621 int size = nb_dofs / Tensor_Dim0;
1622 if (nb_dofs % Tensor_Dim0) {
1626 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1627 auto field_data = getFTensor1FromPtr<Tensor_Dim0>(&*
dotVector.begin());
1629 for (; bb < size; ++bb) {
1630 gradients_at_gauss_pts(
I,
J) += field_data(
I) * diff_base_function(
J);
1632 ++diff_base_function;
1636 for (; bb != nb_base_functions; ++bb)
1637 ++diff_base_function;
1638 ++gradients_at_gauss_pts;
1654 template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
1660 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1661 const EntityType zero_type = MBVERTEX)
1667 template <
int Tensor_Dim0,
int Tensor_Dim1>
1671 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1674 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1675 const EntityType zero_type = MBVERTEX)
1697 template <
int Tensor_Dim0,
int Tensor_Dim1>
1705 "Data pointer not allocated");
1707 const size_t nb_gauss_pts = this->getGaussPts().size2();
1708 constexpr
size_t msize = (Tensor_Dim0 * (Tensor_Dim0 + 1)) / 2;
1709 auto &mat = *this->dataPtr;
1710 if (
type == this->zeroType) {
1711 mat.resize(msize * Tensor_Dim1, nb_gauss_pts,
false);
1720 const int nb_base_functions = data.
getN().size2();
1722 auto gradients_at_gauss_pts =
1723 getFTensor3DgFromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1727 int size = nb_dofs / msize;
1728 if (nb_dofs % msize) {
1730 "Data inconsistency");
1732 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1735 for (; bb < size; ++bb) {
1736 gradients_at_gauss_pts(
I,
J,
K) +=
1737 field_data(
I,
J) * diff_base_function(
K);
1739 ++diff_base_function;
1743 for (; bb != nb_base_functions; ++bb)
1744 ++diff_base_function;
1745 ++gradients_at_gauss_pts;
1757 template <
int Tensor_Dim0,
int Tensor_Dim1>
1760 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1763 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1764 const EntityType zero_type = MBVERTEX)
1770 template <
int Tensor_Dim0,
int Tensor_Dim1>
1773 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1795 template <
int Tensor_Dim0,
int Tensor_Dim1>
1801 "Data pointer not allocated");
1803 const size_t nb_gauss_pts = this->getGaussPts().size2();
1804 auto &mat = *this->dataPtr;
1805 if (
type == this->zeroType) {
1806 mat.resize(Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim1, nb_gauss_pts,
false);
1815 if (this->dataVec.use_count()) {
1816 this->dotVector.resize(nb_dofs,
false);
1817 const double *array;
1818 CHKERR VecGetArrayRead(this->dataVec, &array);
1820 for (
int i = 0;
i != local_indices.size(); ++
i)
1821 if (local_indices[
i] != -1)
1822 this->dotVector[
i] = array[local_indices[
i]];
1824 this->dotVector[
i] = 0;
1825 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1829 const int nb_base_functions = data.
getN().size2();
1831 auto &hessian_base = data.
getN(BaseDerivatives::SecondDerivative);
1833 if (hessian_base.size1() != nb_gauss_pts) {
1835 "Wrong number of integration pts (%d != %d)",
1836 hessian_base.size1(), nb_gauss_pts);
1838 if (hessian_base.size2() !=
1839 nb_base_functions * Tensor_Dim1 * Tensor_Dim1) {
1841 "Wrong number of base functions (%d != %d)",
1842 hessian_base.size2() / (Tensor_Dim1 * Tensor_Dim1),
1845 if (hessian_base.size2() < nb_dofs * Tensor_Dim1 * Tensor_Dim1) {
1847 "Wrong number of base functions (%d < %d)",
1848 hessian_base.size2(), nb_dofs * Tensor_Dim1 * Tensor_Dim1);
1852 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim1, Tensor_Dim1>(
1853 &*hessian_base.data().begin());
1855 auto t_hessian_at_gauss_pts =
1856 getFTensor3FromMat<Tensor_Dim0, Tensor_Dim1, Tensor_Dim1>(mat);
1862 int size = nb_dofs / Tensor_Dim0;
1864 if (nb_dofs % Tensor_Dim0) {
1866 "Data inconsistency");
1870 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1873 for (; bb < size; ++bb) {
1874 t_hessian_at_gauss_pts(
I,
J,
K) +=
1875 field_data(
I) * t_diff2_base_function(
J,
K);
1877 ++t_diff2_base_function;
1881 for (; bb != nb_base_functions; ++bb)
1882 ++t_diff2_base_function;
1883 ++t_hessian_at_gauss_pts;
1886 if (this->dataVec.use_count()) {
1908 template <
int DIM_01,
int DIM_23,
int S = 0>
1920 boost::shared_ptr<MatrixDouble> in_mat,
1921 boost::shared_ptr<MatrixDouble> out_mat,
1922 boost::shared_ptr<MatrixDouble> d_mat)
1935 boost::shared_ptr<MatrixDouble> out_mat,
1936 boost::shared_ptr<MatrixDouble> d_mat)
1949 const size_t nb_gauss_pts =
getGaussPts().size2();
1950 auto t_D = getFTensor4DdgFromMat<DIM_01, DIM_23, S>(*(
dMat));
1951 auto t_in = getFTensor2SymmetricFromMat<DIM_01>(*(
inMat));
1952 outMat->resize((DIM_23 * (DIM_23 + 1)) / 2, nb_gauss_pts,
false);
1953 auto t_out = getFTensor2SymmetricFromMat<DIM_23>(*(
outMat));
1954 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1955 t_out(
i,
j) = t_D(
i,
j,
k,
l) * t_in(
k,
l);
1970 boost::shared_ptr<MatrixDouble>
dMat;
1983 boost::shared_ptr<MatrixDouble> in_mat,
1984 boost::shared_ptr<MatrixDouble> out_mat)
1995 boost::shared_ptr<MatrixDouble> out_mat)
2006 const size_t nb_gauss_pts =
getGaussPts().size2();
2007 auto t_in = getFTensor2FromMat<DIM, DIM>(*(
inMat));
2008 outMat->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts,
false);
2009 auto t_out = getFTensor2SymmetricFromMat<DIM>(*(
outMat));
2010 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
2011 t_out(
i,
j) = (t_in(
i,
j) || t_in(
j,
i)) / 2;
2031 boost::shared_ptr<MatrixDouble> in_mat,
2032 boost::shared_ptr<MatrixDouble> out_mat)
2065 template <
int Base_Dim,
int Field_Dim,
class T,
class L,
class A>
2071 template <
int Field_Dim>
2077 boost::shared_ptr<MatrixDouble> data_ptr,
2079 const EntityType zero_type = MBEDGE,
2080 const int zero_side = 0)
2083 dataPtr(data_ptr), dataVec(data_vec), zeroType(zero_type),
2084 zeroSide(zero_side) {
2090 boost::shared_ptr<MatrixDouble> data_ptr,
2091 const EntityType zero_type = MBEDGE,
2092 const int zero_side = 0)
2114 template <
int Field_Dim>
2120 const size_t nb_integration_points = this->getGaussPts().size2();
2121 if (
type == zeroType && side == zeroSide) {
2122 dataPtr->resize(Field_Dim, nb_integration_points,
false);
2129 if (dataVec.use_count()) {
2130 dotVector.resize(nb_dofs,
false);
2131 const double *array;
2132 CHKERR VecGetArrayRead(dataVec, &array);
2134 for (
int i = 0;
i != local_indices.size(); ++
i)
2135 if (local_indices[
i] != -1)
2136 dotVector[
i] = array[local_indices[
i]];
2139 CHKERR VecRestoreArrayRead(dataVec, &array);
2143 const size_t nb_base_functions = data.
getN().size2() / 3;
2146 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2147 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2150 for (; bb != nb_dofs; ++bb) {
2151 t_data(
i) += t_n_hdiv(
i) * t_dof;
2155 for (; bb != nb_base_functions; ++bb)
2160 if (dataVec.use_count()) {
2170 template <
int Base_Dim,
int Field_Dim = Base_Dim>
2173 Base_Dim, Field_Dim, double, ublas::row_major, DoubleAllocator> {
2182 template <
int Base_Dim,
int Field_Dim = Base_Dim>
2185 template <
int Field_Dim>
2190 boost::shared_ptr<MatrixDouble> data_ptr,
2191 const EntityType zero_type = MBEDGE,
2192 const int zero_side = 0)
2195 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2216 template <
int Field_Dim>
2221 const size_t nb_integration_points = this->getGaussPts().size2();
2222 if (
type == zeroType && side == zeroSide) {
2223 dataPtr->resize(Field_Dim, nb_integration_points,
false);
2228 const size_t nb_dofs = local_indices.size();
2231 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2232 const double *array;
2233 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2234 for (
size_t i = 0;
i != nb_dofs; ++
i)
2235 if (local_indices[
i] != -1)
2236 dot_dofs_vector[
i] = array[local_indices[
i]];
2238 dot_dofs_vector[
i] = 0;
2239 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2241 const size_t nb_base_functions = data.
getN().size2() / 3;
2244 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2245 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2247 for (; bb != nb_dofs; ++bb) {
2248 t_data(
i) += t_n_hdiv(
i) * dot_dofs_vector[bb];
2251 for (; bb != nb_base_functions; ++bb)
2267 template <
int BASE_DIM,
int SPACE_DIM>
2272 boost::shared_ptr<VectorDouble> data_ptr,
2273 const EntityType zero_type = MBEDGE,
2274 const int zero_side = 0)
2285 const size_t nb_integration_points =
getGaussPts().size2();
2287 dataPtr->resize(nb_integration_points,
false);
2293 const size_t nb_base_functions = data.
getN().size2() /
BASE_DIM;
2298 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2301 for (; bb != nb_dofs; ++bb) {
2302 t_data += t_dof * t_n_diff_hdiv(
j,
j);
2306 for (; bb != nb_base_functions; ++bb)
2326 template <
int BASE_DIM,
int SPACE_DIM>
2331 boost::shared_ptr<MatrixDouble> data_ptr,
2332 const EntityType zero_type = MBEDGE,
2333 const int zero_side = 0)
2344 const size_t nb_integration_points =
getGaussPts().size2();
2352 const size_t nb_base_functions = data.
getN().size2() /
BASE_DIM;
2356 auto t_data = getFTensor2FromMat<BASE_DIM, SPACE_DIM>(*
dataPtr);
2357 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2360 for (; bb != nb_dofs; ++bb) {
2361 t_data(
i,
j) += t_dof * t_base_diff(
i,
j);
2365 for (; bb != nb_base_functions; ++bb)
2385 template <
int BASE_DIM,
int SPACE_DIM>
2390 boost::shared_ptr<MatrixDouble> data_ptr,
2391 const EntityType zero_type = MBEDGE,
2392 const int zero_side = 0)
2403 const size_t nb_integration_points =
getGaussPts().size2();
2413 const int nb_base_functions = data.
getN().size2() /
BASE_DIM;
2416 auto &hessian_base = data.
getN(BaseDerivatives::SecondDerivative);
2417 if (hessian_base.size1() != nb_integration_points) {
2419 "Wrong number of integration pts (%d != %d)",
2420 hessian_base.size1(), nb_integration_points);
2422 if (hessian_base.size2() !=
2425 "Wrong number of base functions (%d != %d)",
2431 "Wrong number of base functions (%d < %d)", hessian_base.size2(),
2442 auto t_data = getFTensor3FromMat<BASE_DIM, SPACE_DIM, SPACE_DIM>(*
dataPtr);
2443 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2446 for (; bb != nb_dofs; ++bb) {
2447 t_data(
i,
j,
k) += t_dof * t_base_diff2(
i,
j,
k);
2452 for (; bb != nb_base_functions; ++bb)
2471 template <
int Tensor_Dim1,
int Tensor_Dim2>
2476 boost::shared_ptr<VectorDouble> data_ptr,
2477 const EntityType zero_type = MBEDGE,
2478 const int zero_side = 0)
2489 const size_t nb_integration_points =
getGaussPts().size2();
2491 dataPtr->resize(nb_integration_points,
false);
2496 const int nb_dofs = local_indices.size();
2499 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2500 const double *array;
2502 for (
size_t i = 0;
i != local_indices.size(); ++
i)
2503 if (local_indices[
i] != -1)
2504 dot_dofs_vector[
i] = array[local_indices[
i]];
2506 dot_dofs_vector[
i] = 0;
2509 const size_t nb_base_functions = data.
getN().size2() / Tensor_Dim1;
2513 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2515 for (; bb != nb_dofs; ++bb) {
2517 for (
auto ii = 0; ii != Tensor_Dim2; ++ii)
2518 div += t_n_diff_hdiv(ii, ii);
2519 t_data += dot_dofs_vector[bb] * div;
2522 for (; bb != nb_base_functions; ++bb)
2558 boost::shared_ptr<MatrixDouble> data_ptr,
2559 const EntityType zero_type = MBEDGE,
2560 const int zero_side = 0);
2582 boost::shared_ptr<MatrixDouble> data_ptr,
2583 const EntityType zero_type = MBVERTEX,
2584 const int zero_side = 0);
2607 boost::shared_ptr<MatrixDouble> data_ptr,
2608 const EntityType zero_type = MBVERTEX,
2609 const int zero_side = 0);
2627 template <
int Tensor_Dim0,
int Tensor_Dim1>
2632 boost::shared_ptr<MatrixDouble> data_ptr,
2633 boost::shared_ptr<double> scale_ptr,
2635 const EntityType zero_type = MBEDGE,
2636 const int zero_side = 0)
2646 boost::shared_ptr<MatrixDouble> data_ptr,
2647 const EntityType zero_type = MBEDGE,
2648 const int zero_side = 0)
2656 const size_t nb_integration_points =
getGaussPts().size2();
2658 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points,
false);
2666 const double *array;
2669 for (
int i = 0;
i != local_indices.size(); ++
i)
2670 if (local_indices[
i] != -1)
2679 const size_t nb_base_functions = data.
getN().size2() / 3;
2683 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*
dataPtr);
2684 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2687 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2688 t_data(
i,
j) += (
scale * t_dof(
i)) * t_n_hvec(
j);
2692 for (; bb < nb_base_functions; ++bb)
2718 template <
int Tensor_Dim0,
int Tensor_Dim1>
2725 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
2727 boost::shared_ptr<Range> broken_range_ptr =
nullptr,
2728 boost::shared_ptr<double> scale_ptr =
nullptr,
2729 const EntityType zero_type = MBEDGE,
const int zero_side = 0)
2757 template <
int Tensor_Dim0,
int Tensor_Dim1>
2762 const size_t nb_integration_points = OP::getGaussPts().size2();
2763 if (
type == zeroType) {
2764 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points,
false);
2771 if (dataVec.use_count()) {
2772 dotVector.resize(nb_dofs,
false);
2773 const double *array;
2774 CHKERR VecGetArrayRead(dataVec, &array);
2776 for (
int i = 0;
i != local_indices.size(); ++
i)
2777 if (local_indices[
i] != -1)
2778 dotVector[
i] = array[local_indices[
i]];
2781 CHKERR VecRestoreArrayRead(dataVec, &array);
2792 auto get_get_side_face_dofs = [&]() {
2793 auto fe_type = OP::getFEType();
2797 std::vector<int> side_face_dofs;
2798 side_face_dofs.reserve(data.
getIndices().size() / Tensor_Dim0);
2802 auto it = side_dof_map.get<1>().begin();
2803 it != side_dof_map.get<1>().end(); ++it
2806 if ((Tensor_Dim0 * it->dof) >= data.
getIndices().size()) {
2809 if (it->type == brokenType) {
2810 if (brokenRangePtr) {
2811 auto ent = OP::getSideEntity(it->side, brokenType);
2812 if (brokenRangePtr->find(ent) != brokenRangePtr->end()) {
2813 side_face_dofs.push_back(it->dof);
2816 side_face_dofs.push_back(it->dof);
2821 return side_face_dofs;
2824 auto side_face_dofs = get_get_side_face_dofs();
2828 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
2829 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2830 for (
auto b : side_face_dofs) {
2832 auto t_dof = getFTensor1FromPtr<Tensor_Dim0>(data.
getFieldData().data() +
2834 t_data(
i,
j) += t_dof(
i) * t_row_base(
j);
2838 *dataPtr *= (scalePtr) ? *scalePtr : 1.0;
2840 if (dataVec.use_count()) {
2854 template <
int Tensor_Dim0,
int Tensor_Dim1>
2859 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
2860 boost::shared_ptr<double> scale_ptr,
2862 const EntityType zero_type = MBEDGE,
const int zero_side = 0)
2872 boost::shared_ptr<MatrixDouble> data_ptr,
2873 const EntityType zero_type = MBEDGE,
2874 const int zero_side = 0)
2882 const size_t nb_integration_points =
getGaussPts().size2();
2884 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points,
false);
2893 const double *array;
2896 for (
int i = 0;
i != local_indices.size(); ++
i)
2897 if (local_indices[
i] != -1)
2906 const size_t nb_base_functions =
2907 data.
getN().size2() / (Tensor_Dim0 * Tensor_Dim1);
2910 auto t_n_hten = data.
getFTensor2N<Tensor_Dim0, Tensor_Dim1>();
2911 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*
dataPtr);
2912 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2915 for (; bb != nb_dofs; ++bb) {
2916 t_data(
i,
j) += (
scale * t_dof) * t_n_hten(
i,
j);
2920 for (; bb < nb_base_functions; ++bb)
2948 template <
int Tensor_Dim0,
int Tensor_Dim1,
2954 boost::shared_ptr<MatrixDouble> data_ptr,
2956 const EntityType zero_type = MBEDGE,
2957 const int zero_side = 0)
2967 boost::shared_ptr<MatrixDouble> data_ptr,
2968 const EntityType zero_type = MBEDGE,
2969 const int zero_side = 0)
2976 const size_t nb_integration_points =
getGaussPts().size2();
2978 dataPtr->resize(Tensor_Dim0, nb_integration_points,
false);
2986 const double *array;
2989 for (
int i = 0;
i != local_indices.size(); ++
i)
2990 if (local_indices[
i] != -1)
2998 const size_t nb_base_functions = data.
getN().size2() / 3;
3002 auto t_data = getFTensor1FromMat<Tensor_Dim0>(*
dataPtr);
3005 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
3008 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
3009 double div = t_n_diff_hvec(
j,
j);
3010 t_data(
i) += t_dof(
i) * div;
3012 t_data(
i) += t_base(0) * (t_dof(
i) / t_coords(0));
3018 for (; bb < nb_base_functions; ++bb) {
3051 template <
int Tensor_Dim0,
int Tensor_Dim1,
3059 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
3061 boost::shared_ptr<Range> broken_range_ptr =
nullptr,
3062 boost::shared_ptr<double> scale_ptr =
nullptr,
3063 const EntityType zero_type = MBEDGE)
3075 const size_t nb_integration_points =
getGaussPts().size2();
3077 dataPtr->resize(Tensor_Dim0, nb_integration_points,
false);
3086 const double *array;
3089 for (
int i = 0;
i != local_indices.size(); ++
i)
3090 if (local_indices[
i] != -1)
3105 auto get_get_side_face_dofs = [&]() {
3110 std::vector<int> side_face_dofs;
3111 side_face_dofs.reserve(data.
getIndices().size() / Tensor_Dim0);
3115 auto it = side_dof_map.get<1>().begin();
3116 it != side_dof_map.get<1>().end(); ++it
3119 if ((Tensor_Dim0 * it->dof) >= data.
getIndices().size()) {
3126 side_face_dofs.push_back(it->dof);
3129 side_face_dofs.push_back(it->dof);
3134 return side_face_dofs;
3137 auto side_face_dofs = get_get_side_face_dofs();
3141 auto t_data = getFTensor1FromMat<Tensor_Dim0>(*
dataPtr);
3143 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
3144 for (
auto b : side_face_dofs) {
3145 auto t_dof = getFTensor1FromPtr<Tensor_Dim0>(
3149 double div = t_diff_base(
j,
j);
3150 t_data(
i) += t_dof(
i) * div;
3152 t_data(
i) += t_base(0) * (t_dof(
i) / t_coords(0));
3183 template <
int Tensor_Dim,
typename OpBase>
3187 boost::shared_ptr<MatrixDouble> data_ptr,
3188 boost::shared_ptr<double> scale_ptr,
3189 const EntityType zero_type = MBEDGE,
3190 const int zero_side = 0)
3198 boost::shared_ptr<MatrixDouble> data_ptr,
3199 const EntityType zero_type = MBEDGE,
3200 const int zero_side = 0)
3207 const size_t nb_integration_points = OpBase::getGaussPts().size2();
3209 dataPtr->resize(Tensor_Dim, nb_integration_points,
false);
3215 auto t_normal = OpBase::getFTensor1NormalsAtGaussPts();
3216 const size_t nb_base_functions = data.
getN().size2() / 3;
3218 auto t_data = getFTensor1FromMat<Tensor_Dim>(*
dataPtr);
3219 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
3221 t_normalized_normal(
j) = t_normal(
j);
3225 for (; bb != nb_dofs / Tensor_Dim; ++bb) {
3227 (scale_val * t_dof(
i)) * (t_base(
j) * t_normalized_normal(
j));
3231 for (; bb < nb_base_functions; ++bb) {
3274 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
3277 template <
int D1,
int D2,
int J1,
int J2>
3281 static_assert(D2 == J2,
"Dimension of jacobian and dimension of <out> "
3282 "directive does not match");
3284 size_t nb_functions = diff_n.size2() / D1;
3286 size_t nb_gauss_pts = diff_n.size1();
3291 "Wrong number of Gauss Pts");
3292 if (diff_n.size2() % D1)
3294 "Number of directives of base functions and D1 dimension does "
3298 diffNinvJac.resize(diff_n.size1(), D2 * nb_functions,
false);
3302 auto t_diff_n = getFTensor1FromPtr<D2>(&*
diffNinvJac.data().begin());
3303 auto t_diff_n_ref = getFTensor1FromPtr<D1>(&*diff_n.data().begin());
3304 auto t_inv_jac = getFTensor2FromMat<J1, J2>(*
invJacPtr);
3305 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
3306 for (
size_t dd = 0;
dd != nb_functions; ++
dd) {
3307 t_diff_n(
i) = t_inv_jac(
k,
i) * t_diff_n_ref(
k);
3352 template <
int DERIVARIVE = 1>
3359 template <
int DERIVARIVE = 1>
3366 template <
int DERIVARIVE = 1>
3370 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3374 template <
int DERIVARIVE = 1>
3378 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3396 invJacPtr(inv_jac_ptr) {}
3450 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
3487 boost::shared_ptr<MatrixDouble> jac_ptr)
3687 boost::shared_ptr<VectorDouble> det_ptr,
3688 boost::shared_ptr<MatrixDouble> out_ptr)
3708 "Pointer for inPtr matrix not allocated");
3711 "Pointer for detPtr matrix not allocated");
3713 const auto nb_rows = inPtr->size1();
3714 const auto nb_integration_pts = inPtr->size2();
3718 detPtr->resize(nb_integration_pts,
false);
3719 auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
3721 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
3730 outPtr->resize(nb_rows, nb_integration_pts,
false);
3731 auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
3732 auto t_out = getFTensor2FromMat<DIM, DIM>(*outPtr);
3734 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
3757 boost::shared_ptr<VectorDouble> out_ptr)
3778 "Pointer for inPtr matrix not allocated");
3780 const auto nb_integration_pts = inPtr->size2();
3783 outPtr->resize(nb_integration_pts,
false);
3784 auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
3787 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
3810 boost::shared_ptr<VectorDouble> out_ptr)
3831 "Pointer for inPtr matrix not allocated");
3833 const auto nb_integration_pts = inPtr->size2();
3836 outPtr->resize(nb_integration_pts,
false);
3837 auto t_in = getFTensor2SymmetricFromMat<DIM>(*inPtr);
3840 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
3852 #endif // __USER_DATA_OPERATORS_HPP__