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>
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,
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) {
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) {
422 for (; bb != size; ++bb) {
423 values_at_gauss_pts(
I) += field_data(
I) * base_function;
429 for (; bb < nb_base_functions; ++bb)
431 ++values_at_gauss_pts;
435 if (dataVec.use_count()) {
447template <
int Tensor_Dim>
450 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
453 Tensor_Dim,
double, ublas::row_major,
467template <
int Tensor_Dim, CoordinateTypes COORDINATE_SYSTEM = CARTESIAN>
472 const std::string
field_name, boost::shared_ptr<VectorDouble> data_ptr,
486 if constexpr (COORDINATE_SYSTEM ==
POLAR || COORDINATE_SYSTEM ==
SPHERICAL)
488 "%s coordiante not implemented",
494 vec.resize(nb_gauss_pts,
false);
502 const size_t nb_base_functions = data.
getN().size2();
505 const size_t size = nb_dofs / Tensor_Dim;
507 if (nb_dofs % Tensor_Dim) {
509 "Number of dofs should multiple of dimensions");
514 if constexpr (COORDINATE_SYSTEM ==
CARTESIAN) {
516 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
519 for (; bb != size; ++bb) {
520 values_at_gauss_pts += field_data(
I) * diff_base_function(
I);
522 ++diff_base_function;
526 for (; bb < nb_base_functions; ++bb)
527 ++diff_base_function;
528 ++values_at_gauss_pts;
538 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
541 for (; bb != size; ++bb) {
542 values_at_gauss_pts += field_data(
I) * diff_base_function(
I);
543 values_at_gauss_pts +=
544 field_data(0) * base_function / t_coords(0);
547 ++diff_base_function;
551 for (; bb < nb_base_functions; ++bb) {
553 ++diff_base_function;
555 ++values_at_gauss_pts;
576template <
int Tensor_Dim, PetscData::DataContext CTX>
581 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
594 const size_t nb_dofs = local_indices.size();
599 mat.resize(Tensor_Dim, nb_gauss_pts,
false);
607 auto get_array = [&](
const auto ctx,
auto vec) {
613 <<
"In this case filed degrees of freedom are read from vector. "
614 "That usually happens when time solver is used, and acces to "
615 "first or second rates is needed. You probably not set ts_u, "
616 "ts_u_t, or ts_u_tt and associated data structure, i.e. "
617 "data_ctx to CTX_SET_X, CTX_SET_X_T, or CTX_SET_X_TT "
622 CHKERR VecGetArrayRead(vec, &array);
626 auto restore_array = [&](
auto vec) {
627 return VecRestoreArrayRead(vec, &array);
642 "That case is not implemented");
646 for (
int i = 0;
i != local_indices.size(); ++
i)
647 if (local_indices[
i] != -1)
664 "That case is not implemented");
667 const size_t nb_base_functions = data.
getN().size2();
669 auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
672 const size_t size = nb_dofs / Tensor_Dim;
673 if (nb_dofs % Tensor_Dim) {
676 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
677 auto field_data = getFTensor1FromArray<Tensor_Dim, Tensor_Dim>(
dotVector);
679 for (; bb != size; ++bb) {
680 values_at_gauss_pts(
I) += field_data(
I) * base_function;
686 for (; bb < nb_base_functions; ++bb)
688 ++values_at_gauss_pts;
704template <
int Tensor_Dim>
714template <
int Tensor_Dim>
728template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
734 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
747 boost::shared_ptr<ublas::matrix<T, L, A>>
dataPtr;
751template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
757 "Not implemented for T = %s, dim0 = %d and dim1 = %d",
759 Tensor_Dim0, Tensor_Dim1);
763template <
int Tensor_Dim0,
int Tensor_Dim1>
771 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
784 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
798 boost::shared_ptr<ublas::matrix<double, ublas::row_major, DoubleAllocator>>
805template <
int Tensor_Dim0,
int Tensor_Dim1>
807 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
813 const size_t nb_gauss_pts = data.
getN().size1();
814 if (type == zeroType) {
815 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
821 if (dataVec.use_count()) {
822 dotVector.resize(nb_dofs,
false);
824 CHKERR VecGetArrayRead(dataVec, &array);
826 for (
int i = 0;
i != local_indices.size(); ++
i)
827 if (local_indices[
i] != -1)
828 dotVector[
i] = array[local_indices[
i]];
831 CHKERR VecRestoreArrayRead(dataVec, &array);
836 const size_t nb_base_functions = data.
getN().size2();
838 auto values_at_gauss_pts =
839 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
842 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
843 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
846 for (; bb != size; ++bb) {
847 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
851 for (; bb < nb_base_functions; ++bb)
853 ++values_at_gauss_pts;
856 if (dataVec.use_count()) {
868template <
int Tensor_Dim0,
int Tensor_Dim1>
871 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
874 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
883template <
int Tensor_Dim0,
int Tensor_Dim1>
888 boost::shared_ptr<MatrixDouble> data_ptr,
904 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
908 const size_t nb_dofs = local_indices.size();
913 for (
size_t i = 0;
i != local_indices.size(); ++
i)
914 if (local_indices[
i] != -1)
920 const size_t nb_base_functions = data.
getN().size2();
923 auto values_at_gauss_pts =
924 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
927 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
928 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
929 auto field_data = getFTensorDotData<Tensor_Dim0, Tensor_Dim1>();
931 for (; bb != size; ++bb) {
932 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
936 for (; bb < nb_base_functions; ++bb)
938 ++values_at_gauss_pts;
950 static_assert(Dim0 || !Dim0 || Dim1 || !Dim1,
"not implemented");
958 &dotVector[0], &dotVector[1], &dotVector[2],
960 &dotVector[3], &dotVector[4], &dotVector[5],
962 &dotVector[6], &dotVector[7], &dotVector[8]);
972template <
int Tensor_Dim>
977 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
987 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
989 const int zero_side = 0)
1004 mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts,
false);
1013 const double *array;
1016 for (
int i = 0;
i != local_indices.size(); ++
i)
1017 if (local_indices[
i] != -1)
1025 const int nb_base_functions = data.
getN().size2();
1027 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1030 const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
1031 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1034 for (; bb != size; ++bb) {
1035 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
1039 for (; bb < nb_base_functions; ++bb)
1041 ++values_at_gauss_pts;
1066template <
int Tensor_Dim>
1071 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1085 constexpr auto symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
1087 mat.resize(symm_size, nb_gauss_pts,
false);
1091 const int nb_dofs = local_indices.size();
1099 <<
"In this case filed degrees of freedom are read from vector. "
1100 "That usually happens when time solver is used, and acces to "
1101 "first rates is needed. You probably not set "
1102 "ts_u_t and associated data structure data_ctx to CTX_SET_X_T "
1109 const double *array;
1111 for (
int i = 0;
i != local_indices.size(); ++
i)
1112 if (local_indices[
i] != -1)
1118 const int nb_base_functions = data.
getN().size2();
1121 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1124 const int size = nb_dofs / symm_size;
1125 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1126 auto field_data = getFTensorDotData<Tensor_Dim>();
1128 for (; bb != size; ++bb) {
1129 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
1133 for (; bb < nb_base_functions; ++bb)
1135 ++values_at_gauss_pts;
1148 static_assert(Dim || !Dim,
"not implemented");
1157 &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3], &dotVector[4],
1166 &dotVector[0], &dotVector[1], &dotVector[2]);
1180template <
int Tensor_Dim,
class T,
class L,
class A>
1185 Tensor_Dim,
T,
L,
A>::OpCalculateVectorFieldValues_General;
1192template <
int Tensor_Dim>
1196 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1199 Tensor_Dim,
double, ublas::row_major,
1218template <
int Tensor_Dim>
1220 Tensor_Dim,
double, ublas::row_major,
1225 const size_t nb_gauss_pts = this->getGaussPts().size2();
1226 auto &mat = *this->dataPtr;
1227 if (type == this->zeroType) {
1228 mat.resize(Tensor_Dim, nb_gauss_pts,
false);
1236 const int nb_base_functions = data.
getN().size2();
1238 auto gradients_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
1241 if (nb_dofs > nb_base_functions)
1243 "Number of base functions inconsistent with number of DOFs "
1245 nb_dofs, nb_base_functions);
1247 if (data.
getDiffN().size2() != nb_base_functions * Tensor_Dim)
1250 "Number of base functions inconsistent with number of derivatives "
1252 data.
getDiffN().size2(), nb_base_functions);
1254 if (data.
getDiffN().size1() != nb_gauss_pts)
1257 "Number of base functions inconsistent with number of integration "
1259 data.
getDiffN().size2(), nb_gauss_pts);
1264 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1267 for (; bb != nb_dofs; ++bb) {
1268 gradients_at_gauss_pts(
I) += field_data * diff_base_function(
I);
1270 ++diff_base_function;
1273 for (; bb < nb_base_functions; ++bb)
1274 ++diff_base_function;
1275 ++gradients_at_gauss_pts;
1288template <
int Tensor_Dim>
1291 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1293 Tensor_Dim,
double, ublas::row_major,
1301template <
int Tensor_Dim>
1304 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1307 Tensor_Dim,
double, ublas::row_major,
1321template <
int Tensor_Dim>
1326 const size_t nb_gauss_pts = this->getGaussPts().size2();
1327 auto &mat = *this->dataPtr;
1328 if (type == this->zeroType) {
1329 mat.resize(Tensor_Dim * Tensor_Dim, nb_gauss_pts,
false);
1337 const int nb_base_functions = data.
getN().size2();
1339 auto &hessian_base = data.
getN(BaseDerivatives::SecondDerivative);
1341 if (hessian_base.size1() != nb_gauss_pts) {
1343 "Wrong number of integration pts (%d != %d)",
1344 hessian_base.size1(), nb_gauss_pts);
1346 if (hessian_base.size2() != nb_base_functions * Tensor_Dim * Tensor_Dim) {
1348 "Wrong number of base functions (%d != %d)",
1349 hessian_base.size2() / (Tensor_Dim * Tensor_Dim),
1352 if (hessian_base.size2() < nb_dofs * Tensor_Dim * Tensor_Dim) {
1354 "Wrong number of base functions (%d < %d)",
1355 hessian_base.size2(), nb_dofs * Tensor_Dim * Tensor_Dim);
1359 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim, Tensor_Dim>(
1360 &*hessian_base.data().begin());
1362 auto hessian_at_gauss_pts =
1363 getFTensor2FromMat<Tensor_Dim, Tensor_Dim>(mat);
1367 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1370 for (; bb != nb_dofs; ++bb) {
1371 hessian_at_gauss_pts(
I,
J) +=
1372 field_data * t_diff2_base_function(
I,
J);
1374 ++t_diff2_base_function;
1377 for (; bb < nb_base_functions; ++bb) {
1378 ++t_diff2_base_function;
1381 ++hessian_at_gauss_pts;
1400template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
1406 Tensor_Dim0, Tensor_Dim1,
T,
L,
A>::OpCalculateTensor2FieldValues_General;
1409template <
int Tensor_Dim0,
int Tensor_Dim1>
1413 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1416 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1435template <
int Tensor_Dim0,
int Tensor_Dim1>
1437 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1443 "Data pointer not allocated");
1445 const size_t nb_gauss_pts = this->getGaussPts().size2();
1446 auto &mat = *this->dataPtr;
1447 if (type == this->zeroType) {
1448 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
1457 if (this->dataVec.use_count()) {
1458 this->dotVector.resize(nb_dofs,
false);
1459 const double *array;
1460 CHKERR VecGetArrayRead(this->dataVec, &array);
1462 for (
int i = 0;
i != local_indices.size(); ++
i)
1463 if (local_indices[
i] != -1)
1464 this->dotVector[
i] = array[local_indices[
i]];
1466 this->dotVector[
i] = 0;
1467 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1471 const int nb_base_functions = data.
getN().size2();
1473 auto gradients_at_gauss_pts =
1474 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1477 int size = nb_dofs / Tensor_Dim0;
1478 if (nb_dofs % Tensor_Dim0) {
1480 "Data inconsistency");
1482 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1485 for (; bb < size; ++bb) {
1486 gradients_at_gauss_pts(
I,
J) += field_data(
I) * diff_base_function(
J);
1488 ++diff_base_function;
1492 for (; bb != nb_base_functions; ++bb)
1493 ++diff_base_function;
1494 ++gradients_at_gauss_pts;
1497 if (this->dataVec.use_count()) {
1510template <
int Tensor_Dim0,
int Tensor_Dim1>
1513 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1516 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1525template <
int Tensor_Dim0,
int Tensor_Dim1>
1530 boost::shared_ptr<MatrixDouble> data_ptr,
1544 const int nb_dofs = local_indices.size();
1545 const int nb_gauss_pts = this->
getGaussPts().size2();
1549 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
1556 const double *array;
1558 for (
int i = 0;
i != local_indices.size(); ++
i)
1559 if (local_indices[
i] != -1)
1565 const int nb_base_functions = data.
getN().size2();
1567 auto gradients_at_gauss_pts =
1568 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1571 int size = nb_dofs / Tensor_Dim0;
1572 if (nb_dofs % Tensor_Dim0) {
1576 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1577 auto field_data = getFTensorDotData<Tensor_Dim0>();
1579 for (; bb < size; ++bb) {
1580 gradients_at_gauss_pts(
I,
J) += field_data(
I) * diff_base_function(
J);
1582 ++diff_base_function;
1586 for (; bb != nb_base_functions; ++bb)
1587 ++diff_base_function;
1588 ++gradients_at_gauss_pts;
1599 static_assert(Dim || !Dim,
"not implemented");
1607 &dotVector[0], &dotVector[1], &dotVector[2]);
1622template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
1628 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1635template <
int Tensor_Dim0,
int Tensor_Dim1>
1639 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1642 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1665template <
int Tensor_Dim0,
int Tensor_Dim1>
1667 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1673 "Data pointer not allocated");
1675 const size_t nb_gauss_pts = this->getGaussPts().size2();
1676 constexpr size_t msize = (Tensor_Dim0 * (Tensor_Dim0 + 1)) / 2;
1677 auto &mat = *this->dataPtr;
1678 if (type == this->zeroType) {
1679 mat.resize(msize * Tensor_Dim1, nb_gauss_pts,
false);
1688 const int nb_base_functions = data.
getN().size2();
1690 auto gradients_at_gauss_pts =
1691 getFTensor3DgFromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1695 int size = nb_dofs / msize;
1696 if (nb_dofs % msize) {
1698 "Data inconsistency");
1700 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1703 for (; bb < size; ++bb) {
1704 gradients_at_gauss_pts(
I,
J,
K) +=
1705 field_data(
I,
J) * diff_base_function(
K);
1707 ++diff_base_function;
1711 for (; bb != nb_base_functions; ++bb)
1712 ++diff_base_function;
1713 ++gradients_at_gauss_pts;
1725template <
int Tensor_Dim0,
int Tensor_Dim1>
1728 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1731 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1734 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1738template <
int Tensor_Dim0,
int Tensor_Dim1>
1741 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1744 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1763template <
int Tensor_Dim0,
int Tensor_Dim1>
1769 "Data pointer not allocated");
1771 const size_t nb_gauss_pts = this->getGaussPts().size2();
1772 auto &mat = *this->dataPtr;
1773 if (type == this->zeroType) {
1774 mat.resize(Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim1, nb_gauss_pts,
false);
1783 if (this->dataVec.use_count()) {
1784 this->dotVector.resize(nb_dofs,
false);
1785 const double *array;
1786 CHKERR VecGetArrayRead(this->dataVec, &array);
1788 for (
int i = 0;
i != local_indices.size(); ++
i)
1789 if (local_indices[
i] != -1)
1790 this->dotVector[
i] = array[local_indices[
i]];
1792 this->dotVector[
i] = 0;
1793 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1797 const int nb_base_functions = data.
getN().size2();
1799 auto &hessian_base = data.
getN(BaseDerivatives::SecondDerivative);
1801 if (hessian_base.size1() != nb_gauss_pts) {
1803 "Wrong number of integration pts (%d != %d)",
1804 hessian_base.size1(), nb_gauss_pts);
1806 if (hessian_base.size2() !=
1807 nb_base_functions * Tensor_Dim1 * Tensor_Dim1) {
1809 "Wrong number of base functions (%d != %d)",
1810 hessian_base.size2() / (Tensor_Dim1 * Tensor_Dim1),
1813 if (hessian_base.size2() < nb_dofs * Tensor_Dim1 * Tensor_Dim1) {
1815 "Wrong number of base functions (%d < %d)",
1816 hessian_base.size2(), nb_dofs * Tensor_Dim1 * Tensor_Dim1);
1820 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim1, Tensor_Dim1>(
1821 &*hessian_base.data().begin());
1823 auto t_hessian_at_gauss_pts =
1824 getFTensor3FromMat<Tensor_Dim0, Tensor_Dim1, Tensor_Dim1>(mat);
1830 int size = nb_dofs / Tensor_Dim0;
1832 if (nb_dofs % Tensor_Dim0) {
1834 "Data inconsistency");
1838 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1841 for (; bb < size; ++bb) {
1842 t_hessian_at_gauss_pts(
I,
J,
K) +=
1843 field_data(
I) * t_diff2_base_function(
J,
K);
1845 ++t_diff2_base_function;
1849 for (; bb != nb_base_functions; ++bb)
1850 ++t_diff2_base_function;
1851 ++t_hessian_at_gauss_pts;
1854 if (this->dataVec.use_count()) {
1876template <
int DIM_01,
int DIM_23,
int S = 0>
1884 boost::shared_ptr<MatrixDouble> in_mat,
1885 boost::shared_ptr<MatrixDouble> out_mat,
1886 boost::shared_ptr<MatrixDouble> d_mat)
1900 const size_t nb_gauss_pts =
getGaussPts().size2();
1901 auto t_D = getFTensor4DdgFromMat<DIM_01, DIM_23, S>(*(
dMat));
1902 auto t_in = getFTensor2SymmetricFromMat<DIM_01>(*(
inMat));
1903 outMat->resize((DIM_23 * (DIM_23 + 1)) / 2, nb_gauss_pts,
false);
1904 auto t_out = getFTensor2SymmetricFromMat<DIM_23>(*(
outMat));
1905 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1906 t_out(
i,
j) = t_D(
i,
j,
k,
l) * t_in(
k,
l);
1921 boost::shared_ptr<MatrixDouble>
dMat;
1931 boost::shared_ptr<MatrixDouble> in_mat,
1932 boost::shared_ptr<MatrixDouble> out_mat)
1944 const size_t nb_gauss_pts =
getGaussPts().size2();
1945 auto t_in = getFTensor2FromMat<DIM, DIM>(*(
inMat));
1946 outMat->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts,
false);
1947 auto t_out = getFTensor2SymmetricFromMat<DIM>(*(
outMat));
1948 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1949 t_out(
i,
j) = (t_in(
i,
j) || t_in(
j,
i)) / 2;
1969 boost::shared_ptr<MatrixDouble> in_mat,
1970 boost::shared_ptr<MatrixDouble> out_mat)
2003template <
int Base_Dim,
int Field_Dim,
class T,
class L,
class A>
2009template <
int Field_Dim>
2015 boost::shared_ptr<MatrixDouble> data_ptr,
2017 const int zero_side = 0)
2020 dataPtr(data_ptr), zeroType(
zero_type), zeroSide(zero_side) {
2041template <
int Field_Dim>
2043 3, Field_Dim,
double, ublas::row_major,
2047 const size_t nb_integration_points = this->getGaussPts().size2();
2048 if (type == zeroType && side == zeroSide) {
2049 dataPtr->resize(Field_Dim, nb_integration_points,
false);
2055 const size_t nb_base_functions = data.
getN().size2() / 3;
2058 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2059 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2062 for (; bb != nb_dofs; ++bb) {
2063 t_data(
i) += t_n_hdiv(
i) * t_dof;
2067 for (; bb != nb_base_functions; ++bb)
2078template <
int Base_Dim,
int Field_Dim = Base_Dim>
2081 Base_Dim, Field_Dim, double, ublas::row_major, DoubleAllocator> {
2083 Base_Dim, Field_Dim,
double, ublas::row_major,
2092template <
int Base_Dim,
int Field_Dim = Base_Dim>
2095template <
int Field_Dim>
2100 boost::shared_ptr<MatrixDouble> data_ptr,
2102 const int zero_side = 0)
2105 dataPtr(data_ptr), zeroType(
zero_type), zeroSide(zero_side) {
2126template <
int Field_Dim>
2131 const size_t nb_integration_points = this->getGaussPts().size2();
2132 if (type == zeroType && side == zeroSide) {
2133 dataPtr->resize(Field_Dim, nb_integration_points,
false);
2138 const size_t nb_dofs = local_indices.size();
2141 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2142 const double *array;
2143 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2144 for (
size_t i = 0;
i != nb_dofs; ++
i)
2145 if (local_indices[
i] != -1)
2146 dot_dofs_vector[
i] = array[local_indices[
i]];
2148 dot_dofs_vector[
i] = 0;
2149 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2151 const size_t nb_base_functions = data.
getN().size2() / 3;
2154 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2155 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2157 for (; bb != nb_dofs; ++bb) {
2158 t_data(
i) += t_n_hdiv(
i) * dot_dofs_vector[bb];
2161 for (; bb != nb_base_functions; ++bb)
2177template <
int BASE_DIM,
int SPACE_DIM>
2182 boost::shared_ptr<VectorDouble> data_ptr,
2184 const int zero_side = 0)
2195 const size_t nb_integration_points =
getGaussPts().size2();
2197 dataPtr->resize(nb_integration_points,
false);
2203 const size_t nb_base_functions = data.
getN().size2() /
BASE_DIM;
2208 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2211 for (; bb != nb_dofs; ++bb) {
2212 t_data += t_dof * t_n_diff_hdiv(
j,
j);
2216 for (; bb != nb_base_functions; ++bb)
2236template <
int BASE_DIM,
int SPACE_DIM>
2241 boost::shared_ptr<MatrixDouble> data_ptr,
2243 const int zero_side = 0)
2254 const size_t nb_integration_points =
getGaussPts().size2();
2262 const size_t nb_base_functions = data.
getN().size2() /
BASE_DIM;
2266 auto t_data = getFTensor2FromMat<BASE_DIM, SPACE_DIM>(*
dataPtr);
2267 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2270 for (; bb != nb_dofs; ++bb) {
2271 t_data(
i,
j) += t_dof * t_base_diff(
i,
j);
2275 for (; bb != nb_base_functions; ++bb)
2295template <
int BASE_DIM,
int SPACE_DIM>
2300 boost::shared_ptr<MatrixDouble> data_ptr,
2302 const int zero_side = 0)
2313 const size_t nb_integration_points =
getGaussPts().size2();
2323 const int nb_base_functions = data.
getN().size2() /
BASE_DIM;
2324 auto &hessian_base = data.
getN(BaseDerivatives::SecondDerivative);
2327 if (hessian_base.size1() != nb_integration_points) {
2329 "Wrong number of integration pts (%d != %d)",
2330 hessian_base.size1(), nb_integration_points);
2332 if (hessian_base.size2() !=
2335 "Wrong number of base functions (%d != %d)",
2341 "Wrong number of base functions (%d < %d)", hessian_base.size2(),
2352 auto t_data = getFTensor3FromMat<BASE_DIM, SPACE_DIM, SPACE_DIM>(*
dataPtr);
2353 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2356 for (; bb != nb_dofs; ++bb) {
2357 t_data(
i,
j,
k) += t_dof * t_base_diff2(
i,
j,
k);
2362 for (; bb != nb_base_functions; ++bb)
2381template <
int Tensor_Dim1,
int Tensor_Dim2>
2386 boost::shared_ptr<VectorDouble> data_ptr,
2388 const int zero_side = 0)
2399 const size_t nb_integration_points =
getGaussPts().size2();
2401 dataPtr->resize(nb_integration_points,
false);
2406 const int nb_dofs = local_indices.size();
2409 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2410 const double *array;
2412 for (
size_t i = 0;
i != local_indices.size(); ++
i)
2413 if (local_indices[
i] != -1)
2414 dot_dofs_vector[
i] = array[local_indices[
i]];
2416 dot_dofs_vector[
i] = 0;
2419 const size_t nb_base_functions = data.
getN().size2() / Tensor_Dim1;
2423 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2425 for (; bb != nb_dofs; ++bb) {
2427 for (
auto ii = 0; ii != Tensor_Dim2; ++ii)
2428 div += t_n_diff_hdiv(ii, ii);
2429 t_data += dot_dofs_vector[bb] * div;
2432 for (; bb != nb_base_functions; ++bb)
2454template <
int Base_Dim,
int Space_Dim>
2469 boost::shared_ptr<MatrixDouble> data_ptr,
2471 const int zero_side = 0);
2493 boost::shared_ptr<MatrixDouble> data_ptr,
2495 const int zero_side = 0);
2513template <
int Tensor_Dim0,
int Tensor_Dim1>
2518 boost::shared_ptr<MatrixDouble> data_ptr,
2520 const int zero_side = 0)
2531 const size_t nb_integration_points =
getGaussPts().size2();
2533 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points,
false);
2538 const size_t nb_base_functions = data.
getN().size2() / 3;
2542 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*
dataPtr);
2543 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2546 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2547 t_data(
i,
j) += t_dof(
i) * t_n_hvec(
j);
2551 for (; bb < nb_base_functions; ++bb)
2572template <
int Tensor_Dim0,
int Tensor_Dim1>
2577 boost::shared_ptr<MatrixDouble> data_ptr,
2579 const int zero_side = 0)
2590 const size_t nb_integration_points =
getGaussPts().size2();
2592 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points,
false);
2598 const size_t nb_base_functions =
2599 data.
getN().size2() / (Tensor_Dim0 * Tensor_Dim1);
2602 auto t_n_hten = data.
getFTensor2N<Tensor_Dim0, Tensor_Dim1>();
2603 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*
dataPtr);
2604 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2607 for (; bb != nb_dofs; ++bb) {
2608 t_data(
i,
j) += t_dof * t_n_hten(
i,
j);
2612 for (; bb < nb_base_functions; ++bb)
2632template <
int Tensor_Dim0,
int Tensor_Dim1>
2637 boost::shared_ptr<MatrixDouble> data_ptr,
2639 const int zero_side = 0)
2650 const size_t nb_integration_points =
getGaussPts().size2();
2651 if (type ==
zeroType && side == 0) {
2652 dataPtr->resize(Tensor_Dim0, nb_integration_points,
false);
2657 const size_t nb_base_functions = data.
getN().size2() / 3;
2661 auto t_data = getFTensor1FromMat<Tensor_Dim0>(*
dataPtr);
2662 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2665 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2666 double div = t_n_diff_hvec(
j,
j);
2667 t_data(
i) += t_dof(
i) * div;
2671 for (; bb < nb_base_functions; ++bb)
2691template <
int Tensor_Dim,
typename OpBase>
2695 boost::shared_ptr<MatrixDouble> data_ptr,
2697 const int zero_side = 0)
2707 const size_t nb_integration_points = OpBase::getGaussPts().size2();
2708 if (type ==
zeroType && side == 0) {
2709 dataPtr->resize(Tensor_Dim, nb_integration_points,
false);
2714 auto t_normal = OpBase::getFTensor1Normal();
2715 t_normal(
i) /= sqrt(t_normal(
j) * t_normal(
j));
2716 const size_t nb_base_functions = data.
getN().size2() / 3;
2718 auto t_data = getFTensor1FromMat<Tensor_Dim>(*
dataPtr);
2719 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2722 for (; bb != nb_dofs / Tensor_Dim; ++bb) {
2723 t_data(
i) += t_dof(
i) * (t_base(
j) * t_normal(
j));
2727 for (; bb < nb_base_functions; ++bb)
2749 boost::shared_ptr<MatrixDouble> data_ptr,
2751 const int zero_side = 0)
2761 const size_t nb_integration_points = getGaussPts().size2();
2762 if (type ==
zeroType && side == 0) {
2763 dataPtr->resize(3, nb_integration_points,
false);
2768 auto t_normal = getFTensor1Normal();
2769 t_normal(
i) /= sqrt(t_normal(
j) * t_normal(
j));
2770 const size_t nb_base_functions = data.
getN().size2() / 3;
2772 auto t_data = getFTensor1FromMat<3>(*
dataPtr);
2773 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2775 if (getNormalsAtGaussPts().size1() == nb_integration_points) {
2777 auto t_n = getFTensor1FromPtr<3>(&*
n.data().begin());
2778 t_normal(
i) = t_n(
i) / sqrt(t_n(
j) * t_n(
j));
2781 for (; bb != nb_dofs / 3; ++bb) {
2782 t_data(
i) += t_dof(
i) * (t_base(
j) * t_normal(
j));
2786 for (; bb < nb_base_functions; ++bb)
2826 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
2829 template <
int D1,
int D2,
int J1,
int J2>
2833 static_assert(D2 == J2,
"Dimension of jacobian and dimension of <out> "
2834 "directive does not match");
2836 size_t nb_functions = diff_n.size2() / D1;
2838 size_t nb_gauss_pts = diff_n.size1();
2843 "Wrong number of Gauss Pts");
2844 if (diff_n.size2() % D1)
2846 "Number of direvatives of base functions and D1 dimension does "
2850 diffNinvJac.resize(diff_n.size1(), D2 * nb_functions,
false);
2854 auto t_diff_n = getFTensor1FromPtr<D2>(&*
diffNinvJac.data().begin());
2855 auto t_diff_n_ref = getFTensor1FromPtr<D1>(&*diff_n.data().begin());
2856 auto t_inv_jac = getFTensor2FromMat<J1, J2>(*
invJacPtr);
2857 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
2858 for (
size_t dd = 0; dd != nb_functions; ++dd) {
2859 t_diff_n(
i) = t_inv_jac(
k,
i) * t_diff_n_ref(
k);
2904template <
int DERIVARIVE = 1>
2911template <
int DERIVARIVE = 1>
2918template <
int DERIVARIVE = 1>
2922 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2926template <
int DERIVARIVE = 1>
2930 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2948 invJacPtr(inv_jac_ptr) {}
3002 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
3039 boost::shared_ptr<MatrixDouble> jac_ptr)
3056 2>::OpSetContravariantPiolaTransformOnFace2DImpl;
3240 boost::shared_ptr<VectorDouble> det_ptr,
3241 boost::shared_ptr<MatrixDouble> out_ptr)
3272 "Pointer for inPtr matrix not allocated");
3275 "Pointer for detPtr matrix not allocated");
3277 const auto nb_rows = inPtr->size1();
3278 const auto nb_integration_pts = inPtr->size2();
3282 detPtr->resize(nb_integration_pts,
false);
3283 auto t_in = getFTensor2FromMat<3, 3>(*inPtr);
3285 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
3294 outPtr->resize(nb_rows, nb_integration_pts,
false);
3295 auto t_in = getFTensor2FromMat<3, 3>(*inPtr);
3296 auto t_out = getFTensor2FromMat<3, 3>(*outPtr);
3298 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
3317 "Pointer for inPtr matrix not allocated");
3320 "Pointer for detPtr matrix not allocated");
3322 const auto nb_rows = inPtr->size1();
3323 const auto nb_integration_pts = inPtr->size2();
3327 detPtr->resize(nb_integration_pts,
false);
3328 auto t_in = getFTensor2FromMat<2, 2>(*inPtr);
3330 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
3339 outPtr->resize(nb_rows, nb_integration_pts,
false);
3340 auto t_in = getFTensor2FromMat<2, 2>(*inPtr);
3341 auto t_out = getFTensor2FromMat<2, 2>(*outPtr);
3343 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
ForcesAndSourcesCore::UserDataOperator UserDataOperator
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
#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.
FTensor::Index< 'n', SPACE_DIM > n
#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
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
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
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_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.
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
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
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
OpCalculateHVecTensorDivergence(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.
const EntityHandle zeroType
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.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Index< 'j', 3 > j
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateHVecTensorTrace(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
const EntityHandle zeroType
FTensor::Index< 'i', 3 > i
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
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
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 doWorkImpl(int side, EntityType type, EntitiesFieldData::EntData &data, const FTensor::Number< 3 > &)
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
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSymmetrizeTensor(const std::string field_name, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
FTensor::Index< 'j', DIM > j
boost::shared_ptr< MatrixDouble > inMat
FTensor::Index< 'i', DIM_01 > i
OpTensorTimesSymmetricTensor(const std::string field_name, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat, boost::shared_ptr< MatrixDouble > d_mat)
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
FTensor::Index< 'j', DIM_01 > j
boost::shared_ptr< MatrixDouble > inMat
FTensor::Index< 'l', DIM_23 > l
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