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) {
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]);
969 &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3]);
979template <
int Tensor_Dim>
984 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
994 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
996 const int zero_side = 0)
1011 mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts,
false);
1020 const double *array;
1023 for (
int i = 0;
i != local_indices.size(); ++
i)
1024 if (local_indices[
i] != -1)
1032 const int nb_base_functions = data.
getN().size2();
1034 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1037 const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
1038 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1041 for (; bb != size; ++bb) {
1042 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
1046 for (; bb < nb_base_functions; ++bb)
1048 ++values_at_gauss_pts;
1073template <
int Tensor_Dim>
1078 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1092 constexpr auto symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
1094 mat.resize(symm_size, nb_gauss_pts,
false);
1098 const int nb_dofs = local_indices.size();
1106 <<
"In this case filed degrees of freedom are read from vector. "
1107 "That usually happens when time solver is used, and acces to "
1108 "first rates is needed. You probably not set "
1109 "ts_u_t and associated data structure data_ctx to CTX_SET_X_T "
1116 const double *array;
1118 for (
int i = 0;
i != local_indices.size(); ++
i)
1119 if (local_indices[
i] != -1)
1125 const int nb_base_functions = data.
getN().size2();
1128 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1131 const int size = nb_dofs / symm_size;
1132 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1133 auto field_data = getFTensorDotData<Tensor_Dim>();
1135 for (; bb != size; ++bb) {
1136 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
1140 for (; bb < nb_base_functions; ++bb)
1142 ++values_at_gauss_pts;
1155 static_assert(Dim || !Dim,
"not implemented");
1164 &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3], &dotVector[4],
1173 &dotVector[0], &dotVector[1], &dotVector[2]);
1187template <
int Tensor_Dim,
class T,
class L,
class A>
1192 Tensor_Dim,
T,
L,
A>::OpCalculateVectorFieldValues_General;
1199template <
int Tensor_Dim>
1203 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1206 Tensor_Dim,
double, ublas::row_major,
1225template <
int Tensor_Dim>
1227 Tensor_Dim,
double, ublas::row_major,
1232 const size_t nb_gauss_pts = this->getGaussPts().size2();
1233 auto &mat = *this->dataPtr;
1234 if (type == this->zeroType) {
1235 mat.resize(Tensor_Dim, nb_gauss_pts,
false);
1243 const int nb_base_functions = data.
getN().size2();
1245 auto gradients_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
1248 if (nb_dofs > nb_base_functions)
1250 "Number of base functions inconsistent with number of DOFs "
1252 nb_dofs, nb_base_functions);
1254 if (data.
getDiffN().size2() != nb_base_functions * Tensor_Dim)
1257 "Number of base functions inconsistent with number of derivatives "
1259 data.
getDiffN().size2(), nb_base_functions);
1261 if (data.
getDiffN().size1() != nb_gauss_pts)
1264 "Number of base functions inconsistent with number of integration "
1266 data.
getDiffN().size2(), nb_gauss_pts);
1271 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1274 for (; bb != nb_dofs; ++bb) {
1275 gradients_at_gauss_pts(
I) += field_data * diff_base_function(
I);
1277 ++diff_base_function;
1280 for (; bb < nb_base_functions; ++bb)
1281 ++diff_base_function;
1282 ++gradients_at_gauss_pts;
1295template <
int Tensor_Dim>
1298 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1300 Tensor_Dim,
double, ublas::row_major,
1308template <
int Tensor_Dim>
1311 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1314 Tensor_Dim,
double, ublas::row_major,
1328template <
int Tensor_Dim>
1333 const size_t nb_gauss_pts = this->getGaussPts().size2();
1334 auto &mat = *this->dataPtr;
1335 if (type == this->zeroType) {
1336 mat.resize(Tensor_Dim * Tensor_Dim, nb_gauss_pts,
false);
1344 const int nb_base_functions = data.
getN().size2();
1346 auto &hessian_base = data.
getN(BaseDerivatives::SecondDerivative);
1348 if (hessian_base.size1() != nb_gauss_pts) {
1350 "Wrong number of integration pts (%d != %d)",
1351 hessian_base.size1(), nb_gauss_pts);
1353 if (hessian_base.size2() != nb_base_functions * Tensor_Dim * Tensor_Dim) {
1355 "Wrong number of base functions (%d != %d)",
1356 hessian_base.size2() / (Tensor_Dim * Tensor_Dim),
1359 if (hessian_base.size2() < nb_dofs * Tensor_Dim * Tensor_Dim) {
1361 "Wrong number of base functions (%d < %d)",
1362 hessian_base.size2(), nb_dofs * Tensor_Dim * Tensor_Dim);
1366 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim, Tensor_Dim>(
1367 &*hessian_base.data().begin());
1369 auto hessian_at_gauss_pts =
1370 getFTensor2FromMat<Tensor_Dim, Tensor_Dim>(mat);
1374 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1377 for (; bb != nb_dofs; ++bb) {
1378 hessian_at_gauss_pts(
I,
J) +=
1379 field_data * t_diff2_base_function(
I,
J);
1381 ++t_diff2_base_function;
1384 for (; bb < nb_base_functions; ++bb) {
1385 ++t_diff2_base_function;
1388 ++hessian_at_gauss_pts;
1407template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
1413 Tensor_Dim0, Tensor_Dim1,
T,
L,
A>::OpCalculateTensor2FieldValues_General;
1416template <
int Tensor_Dim0,
int Tensor_Dim1>
1420 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1423 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1442template <
int Tensor_Dim0,
int Tensor_Dim1>
1444 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1450 "Data pointer not allocated");
1452 const size_t nb_gauss_pts = this->getGaussPts().size2();
1453 auto &mat = *this->dataPtr;
1454 if (type == this->zeroType) {
1455 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
1464 if (this->dataVec.use_count()) {
1465 this->dotVector.resize(nb_dofs,
false);
1466 const double *array;
1467 CHKERR VecGetArrayRead(this->dataVec, &array);
1469 for (
int i = 0;
i != local_indices.size(); ++
i)
1470 if (local_indices[
i] != -1)
1471 this->dotVector[
i] = array[local_indices[
i]];
1473 this->dotVector[
i] = 0;
1474 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1478 const int nb_base_functions = data.
getN().size2();
1480 auto gradients_at_gauss_pts =
1481 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1484 int size = nb_dofs / Tensor_Dim0;
1485 if (nb_dofs % Tensor_Dim0) {
1487 "Data inconsistency");
1489 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1492 for (; bb < size; ++bb) {
1493 gradients_at_gauss_pts(
I,
J) += field_data(
I) * diff_base_function(
J);
1495 ++diff_base_function;
1499 for (; bb != nb_base_functions; ++bb)
1500 ++diff_base_function;
1501 ++gradients_at_gauss_pts;
1504 if (this->dataVec.use_count()) {
1517template <
int Tensor_Dim0,
int Tensor_Dim1>
1520 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1523 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1532template <
int Tensor_Dim0,
int Tensor_Dim1>
1537 boost::shared_ptr<MatrixDouble> data_ptr,
1551 const int nb_dofs = local_indices.size();
1552 const int nb_gauss_pts = this->
getGaussPts().size2();
1556 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
1563 const double *array;
1565 for (
int i = 0;
i != local_indices.size(); ++
i)
1566 if (local_indices[
i] != -1)
1572 const int nb_base_functions = data.
getN().size2();
1574 auto gradients_at_gauss_pts =
1575 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1578 int size = nb_dofs / Tensor_Dim0;
1579 if (nb_dofs % Tensor_Dim0) {
1583 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1584 auto field_data = getFTensorDotData<Tensor_Dim0>();
1586 for (; bb < size; ++bb) {
1587 gradients_at_gauss_pts(
I,
J) += field_data(
I) * diff_base_function(
J);
1589 ++diff_base_function;
1593 for (; bb != nb_base_functions; ++bb)
1594 ++diff_base_function;
1595 ++gradients_at_gauss_pts;
1606 static_assert(Dim || !Dim,
"not implemented");
1614 &dotVector[0], &dotVector[1], &dotVector[2]);
1629template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
1635 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1642template <
int Tensor_Dim0,
int Tensor_Dim1>
1646 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1649 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1672template <
int Tensor_Dim0,
int Tensor_Dim1>
1674 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1680 "Data pointer not allocated");
1682 const size_t nb_gauss_pts = this->getGaussPts().size2();
1683 constexpr size_t msize = (Tensor_Dim0 * (Tensor_Dim0 + 1)) / 2;
1684 auto &mat = *this->dataPtr;
1685 if (type == this->zeroType) {
1686 mat.resize(msize * Tensor_Dim1, nb_gauss_pts,
false);
1695 const int nb_base_functions = data.
getN().size2();
1697 auto gradients_at_gauss_pts =
1698 getFTensor3DgFromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1702 int size = nb_dofs / msize;
1703 if (nb_dofs % msize) {
1705 "Data inconsistency");
1707 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1710 for (; bb < size; ++bb) {
1711 gradients_at_gauss_pts(
I,
J,
K) +=
1712 field_data(
I,
J) * diff_base_function(
K);
1714 ++diff_base_function;
1718 for (; bb != nb_base_functions; ++bb)
1719 ++diff_base_function;
1720 ++gradients_at_gauss_pts;
1732template <
int Tensor_Dim0,
int Tensor_Dim1>
1735 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1738 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1741 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1745template <
int Tensor_Dim0,
int Tensor_Dim1>
1748 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1751 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1770template <
int Tensor_Dim0,
int Tensor_Dim1>
1776 "Data pointer not allocated");
1778 const size_t nb_gauss_pts = this->getGaussPts().size2();
1779 auto &mat = *this->dataPtr;
1780 if (type == this->zeroType) {
1781 mat.resize(Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim1, nb_gauss_pts,
false);
1790 if (this->dataVec.use_count()) {
1791 this->dotVector.resize(nb_dofs,
false);
1792 const double *array;
1793 CHKERR VecGetArrayRead(this->dataVec, &array);
1795 for (
int i = 0;
i != local_indices.size(); ++
i)
1796 if (local_indices[
i] != -1)
1797 this->dotVector[
i] = array[local_indices[
i]];
1799 this->dotVector[
i] = 0;
1800 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1804 const int nb_base_functions = data.
getN().size2();
1806 auto &hessian_base = data.
getN(BaseDerivatives::SecondDerivative);
1808 if (hessian_base.size1() != nb_gauss_pts) {
1810 "Wrong number of integration pts (%d != %d)",
1811 hessian_base.size1(), nb_gauss_pts);
1813 if (hessian_base.size2() !=
1814 nb_base_functions * Tensor_Dim1 * Tensor_Dim1) {
1816 "Wrong number of base functions (%d != %d)",
1817 hessian_base.size2() / (Tensor_Dim1 * Tensor_Dim1),
1820 if (hessian_base.size2() < nb_dofs * Tensor_Dim1 * Tensor_Dim1) {
1822 "Wrong number of base functions (%d < %d)",
1823 hessian_base.size2(), nb_dofs * Tensor_Dim1 * Tensor_Dim1);
1827 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim1, Tensor_Dim1>(
1828 &*hessian_base.data().begin());
1830 auto t_hessian_at_gauss_pts =
1831 getFTensor3FromMat<Tensor_Dim0, Tensor_Dim1, Tensor_Dim1>(mat);
1837 int size = nb_dofs / Tensor_Dim0;
1839 if (nb_dofs % Tensor_Dim0) {
1841 "Data inconsistency");
1845 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1848 for (; bb < size; ++bb) {
1849 t_hessian_at_gauss_pts(
I,
J,
K) +=
1850 field_data(
I) * t_diff2_base_function(
J,
K);
1852 ++t_diff2_base_function;
1856 for (; bb != nb_base_functions; ++bb)
1857 ++t_diff2_base_function;
1858 ++t_hessian_at_gauss_pts;
1861 if (this->dataVec.use_count()) {
1883template <
int DIM_01,
int DIM_23,
int S = 0>
1891 boost::shared_ptr<MatrixDouble> in_mat,
1892 boost::shared_ptr<MatrixDouble> out_mat,
1893 boost::shared_ptr<MatrixDouble> d_mat)
1907 const size_t nb_gauss_pts =
getGaussPts().size2();
1908 auto t_D = getFTensor4DdgFromMat<DIM_01, DIM_23, S>(*(
dMat));
1909 auto t_in = getFTensor2SymmetricFromMat<DIM_01>(*(
inMat));
1910 outMat->resize((DIM_23 * (DIM_23 + 1)) / 2, nb_gauss_pts,
false);
1911 auto t_out = getFTensor2SymmetricFromMat<DIM_23>(*(
outMat));
1912 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1913 t_out(
i,
j) = t_D(
i,
j,
k,
l) * t_in(
k,
l);
1928 boost::shared_ptr<MatrixDouble>
dMat;
1938 boost::shared_ptr<MatrixDouble> in_mat,
1939 boost::shared_ptr<MatrixDouble> out_mat)
1951 const size_t nb_gauss_pts =
getGaussPts().size2();
1952 auto t_in = getFTensor2FromMat<DIM, DIM>(*(
inMat));
1953 outMat->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts,
false);
1954 auto t_out = getFTensor2SymmetricFromMat<DIM>(*(
outMat));
1955 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1956 t_out(
i,
j) = (t_in(
i,
j) || t_in(
j,
i)) / 2;
1976 boost::shared_ptr<MatrixDouble> in_mat,
1977 boost::shared_ptr<MatrixDouble> out_mat)
2010template <
int Base_Dim,
int Field_Dim,
class T,
class L,
class A>
2016template <
int Field_Dim>
2022 boost::shared_ptr<MatrixDouble> data_ptr,
2024 const int zero_side = 0)
2027 dataPtr(data_ptr), zeroType(
zero_type), zeroSide(zero_side) {
2048template <
int Field_Dim>
2050 3, Field_Dim,
double, ublas::row_major,
2054 const size_t nb_integration_points = this->getGaussPts().size2();
2055 if (type == zeroType && side == zeroSide) {
2056 dataPtr->resize(Field_Dim, nb_integration_points,
false);
2062 const size_t nb_base_functions = data.
getN().size2() / 3;
2065 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2066 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2069 for (; bb != nb_dofs; ++bb) {
2070 t_data(
i) += t_n_hdiv(
i) * t_dof;
2074 for (; bb != nb_base_functions; ++bb)
2085template <
int Base_Dim,
int Field_Dim = Base_Dim>
2088 Base_Dim, Field_Dim, double, ublas::row_major, DoubleAllocator> {
2090 Base_Dim, Field_Dim,
double, ublas::row_major,
2099template <
int Base_Dim,
int Field_Dim = Base_Dim>
2102template <
int Field_Dim>
2107 boost::shared_ptr<MatrixDouble> data_ptr,
2109 const int zero_side = 0)
2112 dataPtr(data_ptr), zeroType(
zero_type), zeroSide(zero_side) {
2133template <
int Field_Dim>
2138 const size_t nb_integration_points = this->getGaussPts().size2();
2139 if (type == zeroType && side == zeroSide) {
2140 dataPtr->resize(Field_Dim, nb_integration_points,
false);
2145 const size_t nb_dofs = local_indices.size();
2148 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2149 const double *array;
2150 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2151 for (
size_t i = 0;
i != nb_dofs; ++
i)
2152 if (local_indices[
i] != -1)
2153 dot_dofs_vector[
i] = array[local_indices[
i]];
2155 dot_dofs_vector[
i] = 0;
2156 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2158 const size_t nb_base_functions = data.
getN().size2() / 3;
2161 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2162 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2164 for (; bb != nb_dofs; ++bb) {
2165 t_data(
i) += t_n_hdiv(
i) * dot_dofs_vector[bb];
2168 for (; bb != nb_base_functions; ++bb)
2184template <
int BASE_DIM,
int SPACE_DIM>
2189 boost::shared_ptr<VectorDouble> data_ptr,
2191 const int zero_side = 0)
2202 const size_t nb_integration_points =
getGaussPts().size2();
2204 dataPtr->resize(nb_integration_points,
false);
2210 const size_t nb_base_functions = data.
getN().size2() /
BASE_DIM;
2215 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2218 for (; bb != nb_dofs; ++bb) {
2219 t_data += t_dof * t_n_diff_hdiv(
j,
j);
2223 for (; bb != nb_base_functions; ++bb)
2243template <
int BASE_DIM,
int SPACE_DIM>
2248 boost::shared_ptr<MatrixDouble> data_ptr,
2250 const int zero_side = 0)
2261 const size_t nb_integration_points =
getGaussPts().size2();
2269 const size_t nb_base_functions = data.
getN().size2() /
BASE_DIM;
2273 auto t_data = getFTensor2FromMat<BASE_DIM, SPACE_DIM>(*
dataPtr);
2274 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2277 for (; bb != nb_dofs; ++bb) {
2278 t_data(
i,
j) += t_dof * t_base_diff(
i,
j);
2282 for (; bb != nb_base_functions; ++bb)
2302template <
int BASE_DIM,
int SPACE_DIM>
2307 boost::shared_ptr<MatrixDouble> data_ptr,
2309 const int zero_side = 0)
2320 const size_t nb_integration_points =
getGaussPts().size2();
2330 const int nb_base_functions = data.
getN().size2() /
BASE_DIM;
2333 auto &hessian_base = data.
getN(BaseDerivatives::SecondDerivative);
2334 if (hessian_base.size1() != nb_integration_points) {
2336 "Wrong number of integration pts (%d != %d)",
2337 hessian_base.size1(), nb_integration_points);
2339 if (hessian_base.size2() !=
2342 "Wrong number of base functions (%d != %d)",
2348 "Wrong number of base functions (%d < %d)", hessian_base.size2(),
2359 auto t_data = getFTensor3FromMat<BASE_DIM, SPACE_DIM, SPACE_DIM>(*
dataPtr);
2360 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2363 for (; bb != nb_dofs; ++bb) {
2364 t_data(
i,
j,
k) += t_dof * t_base_diff2(
i,
j,
k);
2369 for (; bb != nb_base_functions; ++bb)
2388template <
int Tensor_Dim1,
int Tensor_Dim2>
2393 boost::shared_ptr<VectorDouble> data_ptr,
2395 const int zero_side = 0)
2406 const size_t nb_integration_points =
getGaussPts().size2();
2408 dataPtr->resize(nb_integration_points,
false);
2413 const int nb_dofs = local_indices.size();
2416 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2417 const double *array;
2419 for (
size_t i = 0;
i != local_indices.size(); ++
i)
2420 if (local_indices[
i] != -1)
2421 dot_dofs_vector[
i] = array[local_indices[
i]];
2423 dot_dofs_vector[
i] = 0;
2426 const size_t nb_base_functions = data.
getN().size2() / Tensor_Dim1;
2430 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2432 for (; bb != nb_dofs; ++bb) {
2434 for (
auto ii = 0; ii != Tensor_Dim2; ++ii)
2435 div += t_n_diff_hdiv(ii, ii);
2436 t_data += dot_dofs_vector[bb] * div;
2439 for (; bb != nb_base_functions; ++bb)
2461template <
int Base_Dim,
int Space_Dim>
2476 boost::shared_ptr<MatrixDouble> data_ptr,
2478 const int zero_side = 0);
2500 boost::shared_ptr<MatrixDouble> data_ptr,
2502 const int zero_side = 0);
2520template <
int Tensor_Dim0,
int Tensor_Dim1>
2525 boost::shared_ptr<MatrixDouble> data_ptr,
2527 const int zero_side = 0)
2538 const size_t nb_integration_points =
getGaussPts().size2();
2540 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points,
false);
2545 const size_t nb_base_functions = data.
getN().size2() / 3;
2549 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*
dataPtr);
2550 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2553 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2554 t_data(
i,
j) += t_dof(
i) * t_n_hvec(
j);
2558 for (; bb < nb_base_functions; ++bb)
2579template <
int Tensor_Dim0,
int Tensor_Dim1>
2584 boost::shared_ptr<MatrixDouble> data_ptr,
2586 const int zero_side = 0)
2597 const size_t nb_integration_points =
getGaussPts().size2();
2599 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points,
false);
2605 const size_t nb_base_functions =
2606 data.
getN().size2() / (Tensor_Dim0 * Tensor_Dim1);
2609 auto t_n_hten = data.
getFTensor2N<Tensor_Dim0, Tensor_Dim1>();
2610 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*
dataPtr);
2611 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2614 for (; bb != nb_dofs; ++bb) {
2615 t_data(
i,
j) += t_dof * t_n_hten(
i,
j);
2619 for (; bb < nb_base_functions; ++bb)
2639template <
int Tensor_Dim0,
int Tensor_Dim1,
2645 boost::shared_ptr<MatrixDouble> data_ptr,
2647 const int zero_side = 0)
2658 const size_t nb_integration_points =
getGaussPts().size2();
2659 if (type ==
zeroType && side == 0) {
2660 dataPtr->resize(Tensor_Dim0, nb_integration_points,
false);
2665 const size_t nb_base_functions = data.
getN().size2() / 3;
2669 auto t_data = getFTensor1FromMat<Tensor_Dim0>(*
dataPtr);
2672 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2675 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2676 double div = t_n_diff_hvec(
j,
j);
2677 t_data(
i) += t_dof(
i) * div;
2679 t_data(
i) += t_dof(
i) * (t_base(0) / t_coords(0));
2685 for (; bb < nb_base_functions; ++bb) {
2708template <
int Tensor_Dim,
typename OpBase>
2712 boost::shared_ptr<MatrixDouble> data_ptr,
2714 const int zero_side = 0)
2724 const size_t nb_integration_points = OpBase::getGaussPts().size2();
2725 if (type ==
zeroType && side == 0) {
2726 dataPtr->resize(Tensor_Dim, nb_integration_points,
false);
2731 auto t_normal = OpBase::getFTensor1NormalsAtGaussPts();
2732 const size_t nb_base_functions = data.
getN().size2() / 3;
2734 auto t_data = getFTensor1FromMat<Tensor_Dim>(*
dataPtr);
2735 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2737 t_normalized_normal(
j) = t_normal(
j);
2741 for (; bb != nb_dofs / Tensor_Dim; ++bb) {
2742 t_data(
i) += t_dof(
i) * (t_base(
j) * t_normalized_normal(
j));
2746 for (; bb < nb_base_functions; ++bb) {
2788 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
2791 template <
int D1,
int D2,
int J1,
int J2>
2795 static_assert(D2 == J2,
"Dimension of jacobian and dimension of <out> "
2796 "directive does not match");
2798 size_t nb_functions = diff_n.size2() / D1;
2800 size_t nb_gauss_pts = diff_n.size1();
2805 "Wrong number of Gauss Pts");
2806 if (diff_n.size2() % D1)
2808 "Number of direvatives of base functions and D1 dimension does "
2812 diffNinvJac.resize(diff_n.size1(), D2 * nb_functions,
false);
2816 auto t_diff_n = getFTensor1FromPtr<D2>(&*
diffNinvJac.data().begin());
2817 auto t_diff_n_ref = getFTensor1FromPtr<D1>(&*diff_n.data().begin());
2818 auto t_inv_jac = getFTensor2FromMat<J1, J2>(*
invJacPtr);
2819 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
2820 for (
size_t dd = 0; dd != nb_functions; ++dd) {
2821 t_diff_n(
i) = t_inv_jac(
k,
i) * t_diff_n_ref(
k);
2866template <
int DERIVARIVE = 1>
2873template <
int DERIVARIVE = 1>
2880template <
int DERIVARIVE = 1>
2884 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2888template <
int DERIVARIVE = 1>
2892 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2910 invJacPtr(inv_jac_ptr) {}
2964 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
3001 boost::shared_ptr<MatrixDouble> jac_ptr)
3018 2>::OpSetContravariantPiolaTransformOnFace2DImpl;
3202 boost::shared_ptr<VectorDouble> det_ptr,
3203 boost::shared_ptr<MatrixDouble> out_ptr)
3234 "Pointer for inPtr matrix not allocated");
3237 "Pointer for detPtr matrix not allocated");
3239 const auto nb_rows = inPtr->size1();
3240 const auto nb_integration_pts = inPtr->size2();
3244 detPtr->resize(nb_integration_pts,
false);
3245 auto t_in = getFTensor2FromMat<3, 3>(*inPtr);
3247 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
3256 outPtr->resize(nb_rows, nb_integration_pts,
false);
3257 auto t_in = getFTensor2FromMat<3, 3>(*inPtr);
3258 auto t_out = getFTensor2FromMat<3, 3>(*outPtr);
3260 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
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<2, 2>(*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<2, 2>(*inPtr);
3303 auto t_out = getFTensor2FromMat<2, 2>(*outPtr);
3305 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
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
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
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