6#ifndef __USER_DATA_OPERATORS_HPP__
7#define __USER_DATA_OPERATORS_HPP__
22template <
class T,
class A>
28 boost::shared_ptr<ublas::vector<T, A>> data_ptr,
38 boost::shared_ptr<ublas::vector<T, A>> data_ptr,
57 boost::shared_ptr<ublas::vector<T, A>>
dataPtr;
67template <
class T,
class A>
100 if (type ==
zeroType || vec.size() != nb_gauss_pts) {
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;
155template <PetscData::DataContext CTX>
160 const std::string
field_name, boost::shared_ptr<VectorDouble> data_ptr,
176 if (type ==
zeroAtType || vec.size() != nb_gauss_pts) {
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;
292template <
int Tensor_Dim,
class T,
class L,
class A>
298 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
318 boost::shared_ptr<ublas::matrix<T, L, A>>
dataPtr;
322template <
int Tensor_Dim,
class T,
class L,
class A>
328 "Not implemented for T = %s and dim = %d",
338template <
int Tensor_Dim>
344 boost::shared_ptr<MatrixDouble> data_ptr,
354 boost::shared_ptr<MatrixDouble> data_ptr,
378template <
int Tensor_Dim>
380 Tensor_Dim,
double, ublas::row_major,
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()) {
465template <
int Tensor_Dim>
468 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
471 Tensor_Dim,
double, ublas::row_major,
485template <
int Tensor_Dim, CoordinateTypes COORDINATE_SYSTEM = CARTESIAN>
490 const std::string
field_name, boost::shared_ptr<VectorDouble> data_ptr,
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 field_data(0) * base_function / t_coords(0);
565 ++diff_base_function;
569 for (; bb < nb_base_functions; ++bb) {
571 ++diff_base_function;
573 ++values_at_gauss_pts;
594template <
int Tensor_Dim, PetscData::DataContext CTX>
599 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
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;
722template <
int Tensor_Dim>
732template <
int Tensor_Dim>
746template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
752 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
765 boost::shared_ptr<ublas::matrix<T, L, A>>
dataPtr;
769template <
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);
781template <
int Tensor_Dim0,
int Tensor_Dim1>
789 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
802 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
816 boost::shared_ptr<ublas::matrix<double, ublas::row_major, DoubleAllocator>>
823template <
int Tensor_Dim0,
int Tensor_Dim1>
825 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
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()) {
886template <
int Tensor_Dim0,
int Tensor_Dim1>
889 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
892 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
901template <
int Tensor_Dim0,
int Tensor_Dim1>
906 boost::shared_ptr<MatrixDouble> data_ptr,
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;
977template <
int Tensor_Dim>
982 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
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;
1071template <
int Tensor_Dim>
1076 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
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]);
1185template <
int Tensor_Dim,
class T,
class L,
class A>
1197template <
int Tensor_Dim>
1201 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1204 Tensor_Dim,
double, ublas::row_major,
1223template <
int Tensor_Dim>
1225 Tensor_Dim,
double, ublas::row_major,
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;
1293template <
int Tensor_Dim>
1296 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1298 Tensor_Dim,
double, ublas::row_major,
1306template <
int Tensor_Dim>
1309 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1312 Tensor_Dim,
double, ublas::row_major,
1326template <
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;
1405template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
1414template <
int Tensor_Dim0,
int Tensor_Dim1>
1418 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1421 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1440template <
int Tensor_Dim0,
int Tensor_Dim1>
1442 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
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()) {
1535template <
int Tensor_Dim0,
int Tensor_Dim1>
1538 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1541 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1550template <
int Tensor_Dim0,
int Tensor_Dim1>
1555 boost::shared_ptr<MatrixDouble> data_ptr,
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;
1630template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
1636 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1643template <
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,
1673template <
int Tensor_Dim0,
int Tensor_Dim1>
1675 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
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;
1733template <
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,
1742 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1746template <
int Tensor_Dim0,
int Tensor_Dim1>
1749 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1752 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1771template <
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()) {
1884template <
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)
2040template <
int Base_Dim,
int Field_Dim,
class T,
class L,
class A>
2046template <
int Field_Dim>
2052 boost::shared_ptr<MatrixDouble> data_ptr,
2054 const int zero_side = 0)
2057 dataPtr(data_ptr), zeroType(
zero_type), zeroSide(zero_side) {
2078template <
int Field_Dim>
2080 3, Field_Dim,
double, ublas::row_major,
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)
2115template <
int Base_Dim,
int Field_Dim = Base_Dim>
2118 Base_Dim, Field_Dim, double, ublas::row_major, DoubleAllocator> {
2120 Base_Dim, Field_Dim,
double, ublas::row_major,
2129template <
int Base_Dim,
int Field_Dim = Base_Dim>
2132template <
int Field_Dim>
2137 boost::shared_ptr<MatrixDouble> data_ptr,
2139 const int zero_side = 0)
2142 dataPtr(data_ptr), zeroType(
zero_type), zeroSide(zero_side) {
2163template <
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)
2214template <
int BASE_DIM,
int SPACE_DIM>
2219 boost::shared_ptr<VectorDouble> data_ptr,
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)
2273template <
int BASE_DIM,
int SPACE_DIM>
2278 boost::shared_ptr<MatrixDouble> data_ptr,
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)
2332template <
int BASE_DIM,
int SPACE_DIM>
2337 boost::shared_ptr<MatrixDouble> data_ptr,
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)
2418template <
int Tensor_Dim1,
int Tensor_Dim2>
2423 boost::shared_ptr<VectorDouble> data_ptr,
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)
2491template <
int Base_Dim,
int Space_Dim>
2506 boost::shared_ptr<MatrixDouble> data_ptr,
2508 const int zero_side = 0);
2530 boost::shared_ptr<MatrixDouble> data_ptr,
2532 const int zero_side = 0);
2555 boost::shared_ptr<MatrixDouble> data_ptr,
2557 const int zero_side = 0);
2575template <
int Tensor_Dim0,
int Tensor_Dim1>
2580 boost::shared_ptr<MatrixDouble> data_ptr,
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)
2634template <
int Tensor_Dim0,
int Tensor_Dim1>
2639 boost::shared_ptr<MatrixDouble> data_ptr,
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)
2694template <
int Tensor_Dim0,
int Tensor_Dim1,
2700 boost::shared_ptr<MatrixDouble> data_ptr,
2702 const int zero_side = 0)
2713 const size_t nb_integration_points =
getGaussPts().size2();
2714 if (type ==
zeroType && side == 0) {
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_dof(
i) * (t_base(0) / t_coords(0));
2740 for (; bb < nb_base_functions; ++bb) {
2763template <
int Tensor_Dim,
typename OpBase>
2767 boost::shared_ptr<MatrixDouble> data_ptr,
2769 const int zero_side = 0)
2779 const size_t nb_integration_points = OpBase::getGaussPts().size2();
2780 if (type ==
zeroType && side == 0) {
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);
2921template <
int DERIVARIVE = 1>
2928template <
int DERIVARIVE = 1>
2935template <
int DERIVARIVE = 1>
2939 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2943template <
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)
3073 2>::OpSetContravariantPiolaTransformOnFace2DImpl;
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)
3333 DataForcesAndSourcesCore::EntData &data);
3345 DataForcesAndSourcesCore::EntData &data) {
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)
3387 DataForcesAndSourcesCore::EntData &data);
3399 DataForcesAndSourcesCore::EntData &data) {
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) {
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Tensor1< T, Tensor_Dim > normalize()
static const char *const CoordinateTypesNames[]
Coordinate system names.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
CoordinateTypes
Coodinate system.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VecAllocator< double > DoubleAllocator
implementation of Data Operators for Forces and Sources
static MoFEMErrorCode invertTensor(FTensor::Tensor2< T1, DIM, DIM > &t, T2 &det, FTensor::Tensor2< T3, DIM, DIM > &inv_t)
static auto determinantTensor(FTensor::Tensor2< T, DIM, DIM > &t)
Calculate the determinant of a tensor of rank DIM.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
constexpr IntegrationType I
constexpr auto field_name
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
default operator for EDGE element
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *,(Tensor_Dim *(Tensor_Dim+1))/2 >, Tensor_Dim getFTensor2SymmetricFieldData)()
Return symmetric FTensor rank 2, i.e. matrix from filed data coeffinects.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor3< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 *Tensor_Dim2 >, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > getFTensor3Diff2N(FieldApproximationBase base)
Get second derivatives of base functions for Hvec space.
const VectorInt & getLocalIndices() const
get local indices of dofs on entity
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2FieldData()
Return FTensor rank 2, i.e. matrix from filed data coeffinects.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
default operator for TRI element
default operator for Flat Prism element
default operator for Flat Prism element
FlatPrism finite element.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
@ OPROW
operator doWork function is executed on FE rows
@ OPSPACE
operator do Work is execute on space data
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
structure to get information form mofem into EntitiesFieldData
friend class UserDataOperator
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
boost::shared_ptr< VectorDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateDivergenceVectorFieldValues(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_type=MBVERTEX)
const EntityHandle zeroType
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
boost::shared_ptr< MatrixDouble > dataPtr
const EntityHandle zeroType
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHTensorTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate divergence of tonsorial field using vectorial base.
boost::shared_ptr< MatrixDouble > dataPtr
const EntityHandle zeroType
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHVecTensorDivergence(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
const EntityHandle zeroType
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateHVecTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate trace of vector (Hdiv/Hcurl) space.
const EntityHandle zeroType
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpCalculateHVecTensorTrace(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
FTensor::Index< 'j', Tensor_Dim > j
FTensor::Index< 'i', Tensor_Dim > i
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateHVecVectorField_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
const EntityHandle zeroType
Get vector field for H-div approximation.
const EntityHandle zeroType
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateHVecVectorFieldDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Get vector field for H-div approximation.
Get vector field for H-div approximation.
Calculate gradient of vector field.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
const EntityHandle zeroType
OpCalculateHVecVectorGradient(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate gradient of vector field.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
const EntityHandle zeroType
OpCalculateHVecVectorHessian(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
boost::shared_ptr< MatrixDouble > dataPtr
const EntityHandle zeroType
const EntityHandle zeroType
boost::shared_ptr< MatrixDouble > dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
const EntityHandle zeroType
Calculate curl of vector field.
Calculate divergence of vector field dot.
boost::shared_ptr< VectorDouble > dataPtr
OpCalculateHdivVectorDivergenceDot(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
const EntityHandle zeroType
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate divergence of vector field.
boost::shared_ptr< VectorDouble > dataPtr
OpCalculateHdivVectorDivergence(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
const EntityHandle zeroType
Calculate inverse of jacobian for face element.
OpCalculateInvJacForFatPrism(MatrixDouble &inv_jac)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateInvJacForFatPrism(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
const boost::shared_ptr< MatrixDouble > invJacPtr
Calculate inverse of jacobian for face element.
OpCalculateInvJacForFlatPrism(MatrixDouble &inv_jac_f3)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Evaluate field gradient values for scalar field, i.e. gradient is tensor rank 1 (vector)
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Evaluate field gradient values for scalar field, i.e. gradient is tensor rank 1 (vector),...
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate gradient values of scalar field at integration points
Scalar field values at integration points.
OpCalculateScalarFieldValues_General(const std::string field_name, boost::shared_ptr< ublas::vector< T, A > > data_ptr, SmartPetscObj< Vec > data_vec, const EntityType zero_type=MBVERTEX)
boost::shared_ptr< ublas::vector< T, A > > dataPtr
SmartPetscObj< Vec > dataVec
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of scalar field at integration points
const EntityHandle zeroType
OpCalculateScalarFieldValues_General(const std::string field_name, boost::shared_ptr< ublas::vector< T, A > > data_ptr, const EntityType zero_type=MBVERTEX)
Get rate of scalar field at integration points.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< VectorDouble > dataPtr
OpCalculateScalarFieldValuesFromPetscVecImpl(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
const EntityHandle zeroAtType
Get value at integration points for scalar field.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of scalar field at integration points
SmartPetscObj< Vec > dataVec
OpCalculateTensor2FieldValues_General(const std::string field_name, boost::shared_ptr< ublas::matrix< double, ublas::row_major, DoubleAllocator > > data_ptr, SmartPetscObj< Vec > data_vec, const EntityType zero_type=MBVERTEX)
OpCalculateTensor2FieldValues_General(const std::string field_name, boost::shared_ptr< ublas::matrix< double, ublas::row_major, DoubleAllocator > > data_ptr, const EntityType zero_type=MBVERTEX)
boost::shared_ptr< ublas::matrix< double, ublas::row_major, DoubleAllocator > > dataPtr
const EntityHandle zeroType
Calculate field values for tenor field rank 2.
const EntityHandle zeroType
OpCalculateTensor2FieldValues_General(const std::string field_name, boost::shared_ptr< ublas::matrix< T, L, A > > data_ptr, const EntityType zero_type=MBVERTEX)
boost::shared_ptr< ublas::matrix< T, L, A > > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Get time direvarive values at integration pts for tensor filed rank 2, i.e. matrix field.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateTensor2FieldValuesDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
EntityType zeroAtType
Zero values at Gauss point at this type.
VectorDouble dotVector
Keeps temoorary values of time directives.
boost::shared_ptr< MatrixDouble > dataPtr
Data computed into this matrix.
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
OpCalculateTensor2SymmetricFieldGradient_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Evaluate field gradient values for symmetric 2nd order tensor field, i.e. gradient is tensor rank 3.
OpCalculateTensor2SymmetricFieldGradient_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Get field gradients at integration pts for symmetric tensorial field rank 2.
OpCalculateTensor2SymmetricFieldGradient(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Calculate symmetric tensor field rates ant integratio pts.
const EntityHandle zeroType
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateTensor2SymmetricFieldValuesDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate symmetric tensor field values at integration pts.
OpCalculateTensor2SymmetricFieldValues(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, SmartPetscObj< Vec > data_vec, const EntityType zero_type=MBEDGE, const int zero_side=0)
OpCalculateTensor2SymmetricFieldValues(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
SmartPetscObj< Vec > dataVec
boost::shared_ptr< MatrixDouble > dataPtr
const EntityHandle zeroType
Calculates the trace of an input matrix.
boost::shared_ptr< MatrixDouble > inPtr
FTensor::Index< 'i', DIM > i
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
boost::shared_ptr< VectorDouble > outPtr
OpCalculateTraceFromMat(boost::shared_ptr< MatrixDouble > in_ptr, boost::shared_ptr< VectorDouble > out_ptr)
Calculates the trace of an input matrix.
FTensor::Index< 'i', DIM > i
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
boost::shared_ptr< VectorDouble > outPtr
boost::shared_ptr< MatrixDouble > inPtr
OpCalculateTraceFromSymmMat(boost::shared_ptr< MatrixDouble > in_ptr, boost::shared_ptr< VectorDouble > out_ptr)
Evaluate field gradient values for vector field, i.e. gradient is tensor rank 2.
Get field gradients time derivative at integration pts for scalar filed rank 0, i....
boost::shared_ptr< MatrixDouble > dataPtr
Data computed into this matrix.
OpCalculateVectorFieldGradientDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
EntityType zeroAtType
Zero values at Gauss point at this type.
VectorDouble dotVector
Keeps temoorary values of time directives.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of vector field at integration points
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateVectorFieldValues_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
OpCalculateVectorFieldValues_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, SmartPetscObj< Vec > data_vec, const EntityType zero_type=MBVERTEX)
SmartPetscObj< Vec > dataVec
const EntityHandle zeroType
Calculate field values for tenor field rank 1, i.e. vector field.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of vector field at integration points
OpCalculateVectorFieldValues_General(const std::string field_name, boost::shared_ptr< ublas::matrix< T, L, A > > data_ptr, const EntityType zero_type=MBVERTEX)
const EntityHandle zeroType
boost::shared_ptr< ublas::matrix< T, L, A > > dataPtr
Approximate field values for given petsc vector.
boost::shared_ptr< MatrixDouble > dataPtr
const EntityHandle zeroAtType
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateVectorFieldValuesFromPetscVecImpl(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
Get values at integration pts for tensor filed rank 1, i.e. vector field.
boost::shared_ptr< VectorDouble > detPtr
OpInvertMatrix(boost::shared_ptr< MatrixDouble > in_ptr, boost::shared_ptr< VectorDouble > det_ptr, boost::shared_ptr< MatrixDouble > out_ptr)
boost::shared_ptr< MatrixDouble > outPtr
boost::shared_ptr< MatrixDouble > inPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Make Hdiv space from Hcurl space in 2d.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Operator for fat prism element updating integration weights in the volume.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpMultiplyDeterminantOfJacobianAndWeightsForFatPrisms()
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpScaleMatrix(const std::string field_name, const double scale, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
boost::shared_ptr< MatrixDouble > outMat
boost::shared_ptr< MatrixDouble > inMat
OpSetInvJacH1ForFaceEmbeddedIn3DSpace(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
OpSetInvJacH1ForFace(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Transform local reference derivatives of shape functions to global derivatives.
const boost::shared_ptr< MatrixDouble > invJacPtr
OpSetInvJacH1ForFatPrism(MatrixDouble &inv_jac)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSetInvJacH1ForFatPrism(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Transform local reference derivatives of shape functions to global derivatives.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSetInvJacH1ForFlatPrism(MatrixDouble &inv_jac_f3)
boost::shared_ptr< MatrixDouble > invJacPtr
MatrixDouble diffHcurlInvJac
OpSetInvJacHcurlFaceImpl(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Transform local reference derivatives of shape function to global derivatives for face.
OpSetInvJacL2ForFaceEmbeddedIn3DSpace(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
OpSetInvJacL2ForFace(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Transform local reference derivatives of shape functions to global derivatives.
MoFEMErrorCode applyTransform(MatrixDouble &diff_n)
OpSetInvJacToScalarBasesBasic(FieldSpace space, boost::shared_ptr< MatrixDouble > inv_jac_ptr)
boost::shared_ptr< MatrixDouble > invJacPtr
FTensor::Index< 'i', DIM > i
boost::shared_ptr< MatrixDouble > outMat
OpSymmetrizeTensor(boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Index< 'j', DIM > j
DEPRECATED OpSymmetrizeTensor(const std::string field_name, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
boost::shared_ptr< MatrixDouble > inMat
FTensor::Index< 'i', DIM_01 > i
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Index< 'k', DIM_23 > k
boost::shared_ptr< MatrixDouble > dMat
OpTensorTimesSymmetricTensor(boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat, boost::shared_ptr< MatrixDouble > d_mat)
FTensor::Index< 'j', DIM_01 > j
boost::shared_ptr< MatrixDouble > inMat
FTensor::Index< 'l', DIM_23 > l
DEPRECATED OpTensorTimesSymmetricTensor(const std::string field_name, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat, boost::shared_ptr< MatrixDouble > d_mat)
boost::shared_ptr< MatrixDouble > outMat
static constexpr Switches CtxSetX
static constexpr Switches CtxSetX_TT
static constexpr Switches CtxSetX_T
intrusive_ptr for managing petsc objects