20#ifndef __USER_DATA_OPERATORS_HPP__
21#define __USER_DATA_OPERATORS_HPP__
36template <
class T,
class A>
42 boost::shared_ptr<ublas::vector<T, A>> data_ptr,
52 boost::shared_ptr<ublas::vector<T, A>> data_ptr,
71 boost::shared_ptr<ublas::vector<T, A>>
dataPtr;
81template <
class T,
class A>
115 vec.resize(nb_gauss_pts,
false);
128 for (
int i = 0;
i != local_indices.size(); ++
i)
129 if (local_indices[
i] != -1)
137 const size_t nb_base_functions = data.
getN().size2();
140 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
143 for (; bb != nb_dofs; ++bb) {
144 values_at_gauss_pts += field_data * base_function;
149 for (; bb != nb_base_functions; ++bb)
151 ++values_at_gauss_pts;
169template <PetscData::DataContext CTX>
174 const std::string
field_name, boost::shared_ptr<VectorDouble> data_ptr,
191 vec.resize(nb_gauss_pts,
false);
196 const size_t nb_dofs = local_indices.size();
201 auto get_array = [&](
const auto ctx,
auto vec) {
205 CHKERR VecGetArrayRead(vec, &array);
209 auto restore_array = [&](
auto vec) {
210 return VecRestoreArrayRead(vec, &array);
225 "That case is not implemented");
228 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
229 for (
int i = 0;
i != local_indices.size(); ++
i)
230 if (local_indices[
i] != -1)
231 dot_dofs_vector[
i] = array[local_indices[
i]];
233 dot_dofs_vector[
i] = 0;
247 "That case is not implemented");
250 const size_t nb_base_functions = data.
getN().size2();
254 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
256 for (; bb != nb_dofs; ++bb) {
257 values_at_gauss_pts += dot_dofs_vector[bb] * base_function;
262 for (; bb != nb_base_functions; ++bb)
264 ++values_at_gauss_pts;
295template <
int Tensor_Dim,
class T,
class L,
class A>
301 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
321 boost::shared_ptr<ublas::matrix<T, L, A>>
dataPtr;
325template <
int Tensor_Dim,
class T,
class L,
class A>
331 "Not implemented for T = %s and dim = %d",
341template <
int Tensor_Dim>
347 boost::shared_ptr<MatrixDouble> data_ptr,
368template <
int Tensor_Dim>
375 const size_t nb_gauss_pts = getGaussPts().size2();
376 auto &mat = *dataPtr;
377 if (
type == zeroType) {
378 mat.resize(Tensor_Dim, nb_gauss_pts,
false);
386 const size_t nb_base_functions = data.
getN().size2();
388 auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
390 const size_t size = nb_dofs / Tensor_Dim;
391 if (nb_dofs % Tensor_Dim) {
393 "Data inconsistency");
395 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
398 for (; bb != size; ++bb) {
399 values_at_gauss_pts(
I) += field_data(
I) * base_function;
405 for (; bb != nb_base_functions; ++bb)
407 ++values_at_gauss_pts;
419template <
int Tensor_Dim>
422 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
439template <
int Tensor_Dim, CoordinateTypes COORDINATE_SYSTEM = CARTESIAN>
444 const std::string
field_name, boost::shared_ptr<VectorDouble> data_ptr,
458 if constexpr (COORDINATE_SYSTEM ==
POLAR || COORDINATE_SYSTEM ==
SPHERICAL)
460 "%s coordiante not implemented",
466 vec.resize(nb_gauss_pts,
false);
474 const size_t nb_base_functions = data.
getN().size2();
477 const size_t size = nb_dofs / Tensor_Dim;
479 if (nb_dofs % Tensor_Dim) {
481 "Number of dofs should multiple of dimensions");
486 if constexpr (COORDINATE_SYSTEM ==
CARTESIAN) {
488 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
491 for (; bb != size; ++bb) {
492 values_at_gauss_pts += field_data(
I) * diff_base_function(
I);
494 ++diff_base_function;
498 for (; bb != nb_base_functions; ++bb)
499 ++diff_base_function;
500 ++values_at_gauss_pts;
510 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
513 for (; bb != size; ++bb) {
514 values_at_gauss_pts += field_data(
I) * diff_base_function(
I);
515 values_at_gauss_pts +=
516 field_data(0) * base_function / t_coords(0);
519 ++diff_base_function;
523 for (; bb != nb_base_functions; ++bb) {
525 ++diff_base_function;
527 ++values_at_gauss_pts;
548template <
int Tensor_Dim, PetscData::DataContext CTX>
553 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
566 const size_t nb_dofs = local_indices.size();
571 mat.resize(Tensor_Dim, nb_gauss_pts,
false);
579 auto get_array = [&](
const auto ctx,
auto vec) {
583 CHKERR VecGetArrayRead(vec, &array);
587 auto restore_array = [&](
auto vec) {
588 return VecRestoreArrayRead(vec, &array);
603 "That case is not implemented");
607 for (
int i = 0;
i != local_indices.size(); ++
i)
608 if (local_indices[
i] != -1)
625 "That case is not implemented");
628 const size_t nb_base_functions = data.
getN().size2();
630 auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
633 const size_t size = nb_dofs / Tensor_Dim;
634 if (nb_dofs % Tensor_Dim) {
637 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
638 auto field_data = getFTensor1FromArray<Tensor_Dim, Tensor_Dim>(
dotVector);
640 for (; bb != size; ++bb) {
641 values_at_gauss_pts(
I) += field_data(
I) * base_function;
647 for (; bb != nb_base_functions; ++bb)
649 ++values_at_gauss_pts;
665template <
int Tensor_Dim>
675template <
int Tensor_Dim>
689template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
695 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
708 boost::shared_ptr<ublas::matrix<T, L, A>>
dataPtr;
712template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
718 "Not implemented for T = %s, dim0 = %d and dim1 = %d",
720 Tensor_Dim0, Tensor_Dim1);
724template <
int Tensor_Dim0,
int Tensor_Dim1>
732 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
745 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
759 boost::shared_ptr<ublas::matrix<double, ublas::row_major, DoubleAllocator>>
766template <
int Tensor_Dim0,
int Tensor_Dim1>
774 const size_t nb_gauss_pts = data.
getN().size1();
775 if (
type == zeroType) {
776 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
782 if (dataVec.use_count()) {
783 dotVector.resize(nb_dofs,
false);
785 CHKERR VecGetArrayRead(dataVec, &array);
787 for (
int i = 0;
i != local_indices.size(); ++
i)
788 if (local_indices[
i] != -1)
789 dotVector[
i] = array[local_indices[
i]];
792 CHKERR VecRestoreArrayRead(dataVec, &array);
797 const size_t nb_base_functions = data.
getN().size2();
799 auto values_at_gauss_pts =
800 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
803 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
804 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
807 for (; bb != size; ++bb) {
808 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
812 for (; bb != nb_base_functions; ++bb)
814 ++values_at_gauss_pts;
817 if (dataVec.use_count()) {
829template <
int Tensor_Dim0,
int Tensor_Dim1>
832 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
844template <
int Tensor_Dim0,
int Tensor_Dim1>
849 boost::shared_ptr<MatrixDouble> data_ptr,
865 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
869 const size_t nb_dofs = local_indices.size();
874 for (
size_t i = 0;
i != local_indices.size(); ++
i)
875 if (local_indices[
i] != -1)
881 const size_t nb_base_functions = data.
getN().size2();
884 auto values_at_gauss_pts =
885 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
888 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
889 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
890 auto field_data = getFTensorDotData<Tensor_Dim0, Tensor_Dim1>();
892 for (; bb != size; ++bb) {
893 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
897 for (; bb != nb_base_functions; ++bb)
899 ++values_at_gauss_pts;
911 static_assert(Dim0 || !Dim0 || Dim1 || !Dim1,
"not implemented");
919 &dotVector[0], &dotVector[1], &dotVector[2],
921 &dotVector[3], &dotVector[4], &dotVector[5],
923 &dotVector[6], &dotVector[7], &dotVector[8]);
933template <
int Tensor_Dim>
938 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
948 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
950 const int zero_side = 0)
965 mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts,
false);
977 for (
int i = 0;
i != local_indices.size(); ++
i)
978 if (local_indices[
i] != -1)
986 const int nb_base_functions = data.
getN().size2();
988 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
991 const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
992 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
995 for (; bb != size; ++bb) {
996 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
1000 for (; bb != nb_base_functions; ++bb)
1002 ++values_at_gauss_pts;
1027template <
int Tensor_Dim>
1032 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1047 mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts,
false);
1051 const int nb_dofs = local_indices.size();
1056 const double *array;
1058 for (
int i = 0;
i != local_indices.size(); ++
i)
1059 if (local_indices[
i] != -1)
1065 const int nb_base_functions = data.
getN().size2();
1068 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1071 const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
1072 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1073 auto field_data = getFTensorDotData<Tensor_Dim>();
1075 for (; bb != size; ++bb) {
1076 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
1080 for (; bb != nb_base_functions; ++bb)
1082 ++values_at_gauss_pts;
1095 static_assert(Dim || !Dim,
"not implemented");
1104 &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3], &dotVector[4],
1113 &dotVector[0], &dotVector[1], &dotVector[2]);
1127template <
int Tensor_Dim,
class T,
class L,
class A>
1139template <
int Tensor_Dim>
1143 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1165template <
int Tensor_Dim>
1172 const size_t nb_gauss_pts = this->getGaussPts().size2();
1173 auto &mat = *this->dataPtr;
1174 if (
type == this->zeroType) {
1175 mat.resize(Tensor_Dim, nb_gauss_pts,
false);
1183 const int nb_base_functions = data.
getN().size2();
1185 auto gradients_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
1188 if (nb_dofs > nb_base_functions)
1190 "Number of base functions inconsistent with number of DOFs "
1192 nb_dofs, nb_base_functions);
1194 if (data.
getDiffN().size2() != nb_base_functions * Tensor_Dim)
1197 "Number of base functions inconsistent with number of derivatives "
1199 data.
getDiffN().size2(), nb_base_functions);
1201 if (data.
getDiffN().size1() != nb_gauss_pts)
1204 "Number of base functions inconsistent with number of integration "
1206 data.
getDiffN().size2(), nb_gauss_pts);
1211 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1214 for (; bb != nb_dofs; ++bb) {
1215 gradients_at_gauss_pts(
I) += field_data * diff_base_function(
I);
1217 ++diff_base_function;
1220 for (; bb < nb_base_functions; ++bb)
1221 ++diff_base_function;
1222 ++gradients_at_gauss_pts;
1235template <
int Tensor_Dim>
1238 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1248template <
int Tensor_Dim>
1251 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1268template <
int Tensor_Dim>
1273 const size_t nb_gauss_pts = this->getGaussPts().size2();
1274 auto &mat = *this->dataPtr;
1275 if (
type == this->zeroType) {
1276 mat.resize(Tensor_Dim * Tensor_Dim, nb_gauss_pts,
false);
1284 const int nb_base_functions = data.
getN().size2();
1286 auto &hessian_base = data.
getN(BaseDerivatives::SecondDerivative);
1288 if (hessian_base.size1() != nb_gauss_pts) {
1290 "Wrong number of integration pts (%d != %d)",
1291 hessian_base.size1(), nb_gauss_pts);
1293 if (hessian_base.size2() != nb_base_functions * Tensor_Dim * Tensor_Dim) {
1295 "Wrong number of base functions (%d != %d)",
1296 hessian_base.size2() / (Tensor_Dim * Tensor_Dim),
1299 if (hessian_base.size2() < nb_dofs * Tensor_Dim * Tensor_Dim) {
1301 "Wrong number of base functions (%d < %d)",
1302 hessian_base.size2(), nb_dofs * Tensor_Dim * Tensor_Dim);
1306 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim, Tensor_Dim>(
1307 &*hessian_base.data().begin());
1309 auto hessian_at_gauss_pts =
1310 getFTensor2FromMat<Tensor_Dim, Tensor_Dim>(mat);
1314 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1317 for (; bb != nb_dofs; ++bb) {
1318 hessian_at_gauss_pts(
I,
J) +=
1319 field_data * t_diff2_base_function(
I,
J);
1321 ++t_diff2_base_function;
1324 for (; bb < nb_base_functions; ++bb) {
1325 ++t_diff2_base_function;
1328 ++hessian_at_gauss_pts;
1347template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
1356template <
int Tensor_Dim0,
int Tensor_Dim1>
1360 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1382template <
int Tensor_Dim0,
int Tensor_Dim1>
1390 "Data pointer not allocated");
1392 const size_t nb_gauss_pts = this->getGaussPts().size2();
1393 auto &mat = *this->dataPtr;
1394 if (
type == this->zeroType) {
1395 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
1404 if (this->dataVec.use_count()) {
1405 this->dotVector.resize(nb_dofs,
false);
1406 const double *array;
1407 CHKERR VecGetArrayRead(this->dataVec, &array);
1409 for (
int i = 0;
i != local_indices.size(); ++
i)
1410 if (local_indices[
i] != -1)
1411 this->dotVector[
i] = array[local_indices[
i]];
1413 this->dotVector[
i] = 0;
1414 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1418 const int nb_base_functions = data.
getN().size2();
1420 auto gradients_at_gauss_pts =
1421 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1424 int size = nb_dofs / Tensor_Dim0;
1425 if (nb_dofs % Tensor_Dim0) {
1427 "Data inconsistency");
1429 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1432 for (; bb < size; ++bb) {
1433 gradients_at_gauss_pts(
I,
J) += field_data(
I) * diff_base_function(
J);
1435 ++diff_base_function;
1439 for (; bb != nb_base_functions; ++bb)
1440 ++diff_base_function;
1441 ++gradients_at_gauss_pts;
1444 if (this->dataVec.use_count()) {
1457template <
int Tensor_Dim0,
int Tensor_Dim1>
1460 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1472template <
int Tensor_Dim0,
int Tensor_Dim1>
1477 boost::shared_ptr<MatrixDouble> data_ptr,
1491 const int nb_dofs = local_indices.size();
1492 const int nb_gauss_pts = this->
getGaussPts().size2();
1496 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
1503 const double *array;
1505 for (
int i = 0;
i != local_indices.size(); ++
i)
1506 if (local_indices[
i] != -1)
1512 const int nb_base_functions = data.
getN().size2();
1514 auto gradients_at_gauss_pts =
1515 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1518 int size = nb_dofs / Tensor_Dim0;
1519 if (nb_dofs % Tensor_Dim0) {
1523 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1524 auto field_data = getFTensorDotData<Tensor_Dim0>();
1526 for (; bb < size; ++bb) {
1527 gradients_at_gauss_pts(
I,
J) += field_data(
I) * diff_base_function(
J);
1529 ++diff_base_function;
1533 for (; bb != nb_base_functions; ++bb)
1534 ++diff_base_function;
1535 ++gradients_at_gauss_pts;
1546 static_assert(Dim || !Dim,
"not implemented");
1554 &dotVector[0], &dotVector[1], &dotVector[2]);
1569template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
1575 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1582template <
int Tensor_Dim0,
int Tensor_Dim1>
1586 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1589 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1612template <
int Tensor_Dim0,
int Tensor_Dim1>
1620 "Data pointer not allocated");
1622 const size_t nb_gauss_pts = this->getGaussPts().size2();
1623 constexpr size_t msize = (Tensor_Dim0 * (Tensor_Dim0 + 1)) / 2;
1624 auto &mat = *this->dataPtr;
1625 if (
type == this->zeroType) {
1626 mat.resize(msize * Tensor_Dim1, nb_gauss_pts,
false);
1635 const int nb_base_functions = data.
getN().size2();
1637 auto gradients_at_gauss_pts =
1638 getFTensor3DgFromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1642 int size = nb_dofs / msize;
1643 if (nb_dofs % msize) {
1645 "Data inconsistency");
1647 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1650 for (; bb < size; ++bb) {
1651 gradients_at_gauss_pts(
I,
J,
K) +=
1652 field_data(
I,
J) * diff_base_function(
K);
1654 ++diff_base_function;
1658 for (; bb != nb_base_functions; ++bb)
1659 ++diff_base_function;
1660 ++gradients_at_gauss_pts;
1672template <
int Tensor_Dim0,
int Tensor_Dim1>
1675 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1678 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1685template <
int Tensor_Dim0,
int Tensor_Dim1>
1688 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1710template <
int Tensor_Dim0,
int Tensor_Dim1>
1716 "Data pointer not allocated");
1718 const size_t nb_gauss_pts = this->getGaussPts().size2();
1719 auto &mat = *this->dataPtr;
1720 if (
type == this->zeroType) {
1721 mat.resize(Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim1, nb_gauss_pts,
false);
1730 if (this->dataVec.use_count()) {
1731 this->dotVector.resize(nb_dofs,
false);
1732 const double *array;
1733 CHKERR VecGetArrayRead(this->dataVec, &array);
1735 for (
int i = 0;
i != local_indices.size(); ++
i)
1736 if (local_indices[
i] != -1)
1737 this->dotVector[
i] = array[local_indices[
i]];
1739 this->dotVector[
i] = 0;
1740 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1744 const int nb_base_functions = data.
getN().size2();
1746 auto &hessian_base = data.
getN(BaseDerivatives::SecondDerivative);
1748 if (hessian_base.size1() != nb_gauss_pts) {
1750 "Wrong number of integration pts (%d != %d)",
1751 hessian_base.size1(), nb_gauss_pts);
1753 if (hessian_base.size2() !=
1754 nb_base_functions * Tensor_Dim1 * Tensor_Dim1) {
1756 "Wrong number of base functions (%d != %d)",
1757 hessian_base.size2() / (Tensor_Dim1 * Tensor_Dim1),
1760 if (hessian_base.size2() < nb_dofs * Tensor_Dim1 * Tensor_Dim1) {
1762 "Wrong number of base functions (%d < %d)",
1763 hessian_base.size2(), nb_dofs * Tensor_Dim1 * Tensor_Dim1);
1767 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim1, Tensor_Dim1>(
1768 &*hessian_base.data().begin());
1770 auto t_hessian_at_gauss_pts =
1771 getFTensor3FromMat<Tensor_Dim0, Tensor_Dim1, Tensor_Dim1>(mat);
1777 int size = nb_dofs / Tensor_Dim0;
1779 if (nb_dofs % Tensor_Dim0) {
1781 "Data inconsistency");
1785 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1788 for (; bb < size; ++bb) {
1789 t_hessian_at_gauss_pts(
I,
J,
K) +=
1790 field_data(
I) * t_diff2_base_function(
J,
K);
1792 ++t_diff2_base_function;
1796 for (; bb != nb_base_functions; ++bb)
1797 ++t_diff2_base_function;
1798 ++t_hessian_at_gauss_pts;
1801 if (this->dataVec.use_count()) {
1823template <
int DIM_01,
int DIM_23,
int S = 0>
1831 boost::shared_ptr<MatrixDouble> in_mat,
1832 boost::shared_ptr<MatrixDouble> out_mat,
1833 boost::shared_ptr<MatrixDouble> d_mat)
1847 const size_t nb_gauss_pts =
getGaussPts().size2();
1848 auto t_D = getFTensor4DdgFromMat<DIM_01, DIM_23, S>(*(
dMat));
1849 auto t_in = getFTensor2SymmetricFromMat<DIM_01>(*(
inMat));
1850 outMat->resize((DIM_23 * (DIM_23 + 1)) / 2, nb_gauss_pts,
false);
1851 auto t_out = getFTensor2SymmetricFromMat<DIM_23>(*(
outMat));
1852 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1853 t_out(
i,
j) = t_D(
i,
j,
k,
l) * t_in(
k,
l);
1868 boost::shared_ptr<MatrixDouble>
dMat;
1878 boost::shared_ptr<MatrixDouble> in_mat,
1879 boost::shared_ptr<MatrixDouble> out_mat)
1891 const size_t nb_gauss_pts =
getGaussPts().size2();
1892 auto t_in = getFTensor2FromMat<DIM, DIM>(*(
inMat));
1893 outMat->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts,
false);
1894 auto t_out = getFTensor2SymmetricFromMat<DIM>(*(
outMat));
1895 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1896 t_out(
i,
j) = (t_in(
i,
j) || t_in(
j,
i)) / 2;
1916 boost::shared_ptr<MatrixDouble> in_mat,
1917 boost::shared_ptr<MatrixDouble> out_mat)
1950template <
int Tensor_Dim0,
class T,
class L,
class A>
1956 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
1976 boost::shared_ptr<ublas::matrix<T, L, A>>
dataPtr;
1981template <
int Tensor_Dim,
class T,
class L,
class A>
1986 "Not implemented for T = %s and dim = %d",
1995template <
int Tensor_Dim>
2001 boost::shared_ptr<MatrixDouble> data_ptr,
2003 const int zero_side = 0)
2027template <
int Tensor_Dim>
2033 const size_t nb_integration_points = this->getGaussPts().size2();
2034 if (
type == zeroType && side == zeroSide) {
2035 dataPtr->resize(Tensor_Dim, nb_integration_points,
false);
2041 const size_t nb_base_functions = data.
getN().size2() / Tensor_Dim;
2044 auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
2045 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2048 for (; bb != nb_dofs; ++bb) {
2049 t_data(
i) += t_n_hdiv(
i) * t_dof;
2053 for (; bb != nb_base_functions; ++bb)
2064template <
int Tensor_Dim>
2067 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
2076template <
int Tensor_Dim>
2081 boost::shared_ptr<MatrixDouble> data_ptr,
2083 const int zero_side = 0)
2107template <
int Tensor_Dim>
2111 const size_t nb_integration_points = this->getGaussPts().size2();
2112 if (
type == zeroType && side == zeroSide) {
2113 dataPtr->resize(Tensor_Dim, nb_integration_points,
false);
2118 const size_t nb_dofs = local_indices.size();
2121 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2122 const double *array;
2123 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2124 for (
size_t i = 0;
i != nb_dofs; ++
i)
2125 if (local_indices[
i] != -1)
2126 dot_dofs_vector[
i] = array[local_indices[
i]];
2128 dot_dofs_vector[
i] = 0;
2129 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2131 const size_t nb_base_functions = data.
getN().size2() / 3;
2134 auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
2135 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2137 for (; bb != nb_dofs; ++bb) {
2138 t_data(
i) += t_n_hdiv(
i) * dot_dofs_vector[bb];
2141 for (; bb != nb_base_functions; ++bb)
2156template <
int BASE_DIM,
int SPACE_DIM>
2161 boost::shared_ptr<VectorDouble> data_ptr,
2163 const int zero_side = 0)
2174 const size_t nb_integration_points =
getGaussPts().size2();
2176 dataPtr->resize(nb_integration_points,
false);
2182 const size_t nb_base_functions = data.
getN().size2() /
BASE_DIM;
2187 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2190 for (; bb != nb_dofs; ++bb) {
2191 t_data += t_dof * t_n_diff_hdiv(
j,
j);
2195 for (; bb != nb_base_functions; ++bb)
2215template <
int BASE_DIM,
int SPACE_DIM>
2220 boost::shared_ptr<MatrixDouble> data_ptr,
2222 const int zero_side = 0)
2233 const size_t nb_integration_points =
getGaussPts().size2();
2241 const size_t nb_base_functions = data.
getN().size2() /
BASE_DIM;
2245 auto t_data = getFTensor2FromMat<BASE_DIM, SPACE_DIM>(*
dataPtr);
2246 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2249 for (; bb != nb_dofs; ++bb) {
2250 t_data(
i,
j) += t_dof * t_base_diff(
i,
j);
2254 for (; bb != nb_base_functions; ++bb)
2274template <
int BASE_DIM,
int SPACE_DIM>
2279 boost::shared_ptr<MatrixDouble> data_ptr,
2281 const int zero_side = 0)
2292 const size_t nb_integration_points =
getGaussPts().size2();
2302 const int nb_base_functions = data.
getN().size2() /
BASE_DIM;
2303 auto &hessian_base = data.
getN(BaseDerivatives::SecondDerivative);
2306 if (hessian_base.size1() != nb_integration_points) {
2308 "Wrong number of integration pts (%d != %d)",
2309 hessian_base.size1(), nb_integration_points);
2311 if (hessian_base.size2() !=
2314 "Wrong number of base functions (%d != %d)",
2320 "Wrong number of base functions (%d < %d)", hessian_base.size2(),
2330 auto t_data = getFTensor3FromMat<BASE_DIM, SPACE_DIM, SPACE_DIM>(*
dataPtr);
2331 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2334 for (; bb != nb_dofs; ++bb) {
2335 t_data(
i,
j,
k) += t_dof * t_base_diff2(
i,
j,
k);
2340 for (; bb != nb_base_functions; ++bb)
2359template <
int Tensor_Dim1,
int Tensor_Dim2>
2364 boost::shared_ptr<VectorDouble> data_ptr,
2366 const int zero_side = 0)
2377 const size_t nb_integration_points =
getGaussPts().size2();
2379 dataPtr->resize(nb_integration_points,
false);
2384 const int nb_dofs = local_indices.size();
2387 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2388 const double *array;
2390 for (
size_t i = 0;
i != local_indices.size(); ++
i)
2391 if (local_indices[
i] != -1)
2392 dot_dofs_vector[
i] = array[local_indices[
i]];
2394 dot_dofs_vector[
i] = 0;
2397 const size_t nb_base_functions = data.
getN().size2() / Tensor_Dim1;
2401 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2403 for (; bb != nb_dofs; ++bb) {
2405 for (
auto ii = 0; ii != Tensor_Dim2; ++ii)
2406 div += t_n_diff_hdiv(ii, ii);
2407 t_data += dot_dofs_vector[bb] * div;
2410 for (; bb != nb_base_functions; ++bb)
2430template <
int Tensor_Dim>
2435 boost::shared_ptr<MatrixDouble> data_ptr,
2437 const int zero_side = 0)
2448 const size_t nb_integration_points =
getGaussPts().size2();
2450 dataPtr->resize(Tensor_Dim, nb_integration_points,
false);
2459 const size_t nb_base_functions = data.
getN().size2() / Tensor_Dim;
2464 auto t_data = getFTensor1FromMat<Tensor_Dim>(*
dataPtr);
2465 for (
int gg = 0; gg != nb_integration_points; ++gg) {
2469 for (; bb != nb_dofs; ++bb) {
2476 for (; bb != nb_base_functions; ++bb)
2497template <
int Tensor_Dim0,
int Tensor_Dim1>
2502 boost::shared_ptr<MatrixDouble> data_ptr,
2504 const int zero_side = 0)
2515 const size_t nb_integration_points =
getGaussPts().size2();
2517 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points,
false);
2522 const size_t nb_base_functions = data.
getN().size2() / 3;
2526 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*
dataPtr);
2527 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2530 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2531 t_data(
i,
j) += t_dof(
i) * t_n_hvec(
j);
2535 for (; bb != nb_base_functions; ++bb)
2556template <
int Tensor_Dim0,
int Tensor_Dim1>
2561 boost::shared_ptr<MatrixDouble> data_ptr,
2563 const int zero_side = 0)
2574 const size_t nb_integration_points =
getGaussPts().size2();
2576 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points,
false);
2582 const size_t nb_base_functions =
2583 data.
getN().size2() / (Tensor_Dim0 * Tensor_Dim1);
2586 auto t_n_hten = data.
getFTensor2N<Tensor_Dim0, Tensor_Dim1>();
2587 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*
dataPtr);
2588 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2591 for (; bb != nb_dofs; ++bb) {
2592 t_data(
i,
j) += t_dof * t_n_hten(
i,
j);
2596 for (; bb != nb_base_functions; ++bb)
2616template <
int Tensor_Dim0,
int Tensor_Dim1>
2621 boost::shared_ptr<MatrixDouble> data_ptr,
2623 const int zero_side = 0)
2634 const size_t nb_integration_points =
getGaussPts().size2();
2636 dataPtr->resize(Tensor_Dim0, nb_integration_points,
false);
2641 const size_t nb_base_functions = data.
getN().size2() / 3;
2645 auto t_data = getFTensor1FromMat<Tensor_Dim0>(*
dataPtr);
2646 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2649 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2650 double div = t_n_diff_hvec(
j,
j);
2651 t_data(
i) += t_dof(
i) * div;
2655 for (; bb < nb_base_functions; ++bb)
2675template <
int Tensor_Dim,
typename OpBase>
2679 boost::shared_ptr<MatrixDouble> data_ptr,
2681 const int zero_side = 0)
2691 const size_t nb_integration_points = OpBase::getGaussPts().size2();
2693 dataPtr->resize(Tensor_Dim, nb_integration_points,
false);
2698 auto t_normal = OpBase::getFTensor1Normal();
2699 t_normal(
i) /= sqrt(t_normal(
j) * t_normal(
j));
2700 const size_t nb_base_functions = data.
getN().size2() / 3;
2702 auto t_data = getFTensor1FromMat<Tensor_Dim>(*
dataPtr);
2703 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2706 for (; bb != nb_dofs / Tensor_Dim; ++bb) {
2707 t_data(
i) += t_dof(
i) * (t_base(
j) * t_normal(
j));
2711 for (; bb < nb_base_functions; ++bb)
2733 boost::shared_ptr<MatrixDouble> data_ptr,
2735 const int zero_side = 0)
2745 const size_t nb_integration_points = getGaussPts().size2();
2747 dataPtr->resize(3, nb_integration_points,
false);
2752 auto t_normal = getFTensor1Normal();
2753 t_normal(
i) /= sqrt(t_normal(
j) * t_normal(
j));
2754 const size_t nb_base_functions = data.
getN().size2() / 3;
2756 auto t_data = getFTensor1FromMat<3>(*
dataPtr);
2757 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2759 if (getNormalsAtGaussPts().size1() == nb_integration_points) {
2762 t_normal(
i) = t_n(
i) / sqrt(t_n(
j) * t_n(
j));
2765 for (; bb != nb_dofs / 3; ++bb) {
2766 t_data(
i) += t_dof(
i) * (t_base(
j) * t_normal(
j));
2770 for (; bb < nb_base_functions; ++bb)
2810 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
2813 template <
int D1,
int D2,
int J1,
int J2>
2816 size_t nb_functions = diff_n.size2() / D1;
2818 size_t nb_gauss_pts = diff_n.size1();
2823 "Wrong number of Gauss Pts");
2824 if (diff_n.size2() % D1)
2826 "Number of direvatives of base functions and D1 dimension does "
2830 diffNinvJac.resize(diff_n.size1(), diff_n.size2(),
false);
2834 auto t_diff_n = getFTensor1FromPtr<D2>(&*
diffNinvJac.data().begin());
2835 auto t_diff_n_ref = getFTensor1FromPtr<D1>(&*diff_n.data().begin());
2836 auto t_inv_jac = getFTensor2FromMat<J1, J2>(*
invJacPtr);
2837 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
2838 for (
size_t dd = 0;
dd != nb_functions; ++
dd) {
2839 t_diff_n(
i) = t_inv_jac(
k,
i) * t_diff_n_ref(
k);
2884template <
int DERIVARIVE = 1>
2891template <
int DERIVARIVE = 1>
2898template <
int DERIVARIVE = 1>
2902 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2906template <
int DERIVARIVE = 1>
2910 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2928 invJacPtr(inv_jac_ptr) {}
2984 boost::shared_ptr<MatrixDouble> jac_ptr)
3001 2>::OpSetContravariantPiolaTransformOnFace2DImpl;
3187 boost::shared_ptr<VectorDouble> det_ptr,
3188 boost::shared_ptr<MatrixDouble> out_ptr)
3219 "Pointer for inPtr matrix not allocated");
3222 "Pointer for detPtr matrix not allocated");
3224 const auto nb_rows = inPtr->size1();
3225 const auto nb_integration_pts = inPtr->size2();
3229 detPtr->resize(nb_integration_pts,
false);
3230 auto t_in = getFTensor2FromMat<3, 3>(*inPtr);
3232 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
3241 outPtr->resize(nb_rows, nb_integration_pts,
false);
3242 auto t_in = getFTensor2FromMat<3, 3>(*inPtr);
3243 auto t_out = getFTensor2FromMat<3, 3>(*outPtr);
3245 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
3264 "Pointer for inPtr matrix not allocated");
3267 "Pointer for detPtr matrix not allocated");
3269 const auto nb_rows = inPtr->size1();
3270 const auto nb_integration_pts = inPtr->size2();
3274 detPtr->resize(nb_integration_pts,
false);
3275 auto t_in = getFTensor2FromMat<2, 2>(*inPtr);
3277 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
3286 outPtr->resize(nb_rows, nb_integration_pts,
false);
3287 auto t_in = getFTensor2FromMat<2, 2>(*inPtr);
3288 auto t_out = getFTensor2FromMat<2, 2>(*outPtr);
3290 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
ForcesAndSourcesCore::UserDataOperator UserDataOperator
EntitiesFieldData::EntData EntData
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< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
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
FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > getFTensor1FromPtr< 3 >(double *ptr)
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
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
friend class UserDataOperator
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
boost::shared_ptr< VectorDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateDivergenceVectorFieldValues(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_type=MBVERTEX)
const EntityHandle zeroType
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
boost::shared_ptr< MatrixDouble > dataPtr
const EntityHandle zeroType
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHTensorTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate divergence of tonsorial field using vectorial base.
boost::shared_ptr< MatrixDouble > dataPtr
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.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of vector field at integration points
OpCalculateHVecVectorField_General(const std::string field_name, boost::shared_ptr< ublas::matrix< T, L, A > > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
const EntityHandle zeroType
boost::shared_ptr< ublas::matrix< T, L, A > > dataPtr
Get vector field for H-div approximation.
OpCalculateHVecVectorFieldDot(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)
Calculate values of vector field at integration points.
const EntityHandle zeroType
boost::shared_ptr< MatrixDouble > dataPtr
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)
Calculate curl of vector field.
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateHcurlVectorCurl(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 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)
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 valuse 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