96#ifndef __USER_DATA_OPERATORS_HPP__
97 #define __USER_DATA_OPERATORS_HPP__
123template <
class T,
class A>
136 boost::shared_ptr<ublas::vector<T, A>> data_ptr,
137 const EntityType zero_type = MBVERTEX)
154 boost::shared_ptr<ublas::vector<T, A>> data_ptr,
173 boost::shared_ptr<ublas::vector<T, A>>
dataPtr;
183template <
class T,
class A>
221 if (type ==
zeroType || vec.size() != nb_gauss_pts) {
222 vec.resize(nb_gauss_pts,
false);
235 for (
int i = 0;
i != local_indices.size(); ++
i)
236 if (local_indices[
i] != -1)
244 const size_t nb_base_functions = data.
getN().size2();
247 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
250 for (; bb != nb_dofs; ++bb) {
251 values_at_gauss_pts += field_data * base_function;
256 for (; bb < nb_base_functions; ++bb)
258 ++values_at_gauss_pts;
281template <PetscData::DataContext CTX>
293 const std::string
field_name, boost::shared_ptr<VectorDouble> data_ptr,
294 const EntityType zero_at_type = MBVERTEX)
309 if (type ==
zeroAtType || vec.size() != nb_gauss_pts) {
310 vec.resize(nb_gauss_pts,
false);
315 const size_t nb_dofs = local_indices.size();
320 auto get_array = [&](
const auto ctx,
auto vec) {
326 <<
"In this case field degrees of freedom are read from vector. "
327 "That usually happens when time solver is used, and acces to "
328 "first or second rates is needed. You probably not set ts_u, "
329 "ts_u_t, or ts_u_tt and associated data structure, i.e. "
330 "data_ctx to CTX_SET_X, CTX_SET_X_T, or CTX_SET_X_TT "
335 CHKERR VecGetArrayRead(vec, &array);
339 auto restore_array = [&](
auto vec) {
340 return VecRestoreArrayRead(vec, &array);
355 "That case is not implemented");
358 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
359 for (
int i = 0;
i != local_indices.size(); ++
i)
360 if (local_indices[
i] != -1)
361 dot_dofs_vector[
i] = array[local_indices[
i]];
363 dot_dofs_vector[
i] = 0;
377 "That case is not implemented");
380 const size_t nb_base_functions = data.
getN().size2();
384 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
386 for (; bb != nb_dofs; ++bb) {
387 values_at_gauss_pts += dot_dofs_vector[bb] * base_function;
392 for (; bb < nb_base_functions; ++bb)
394 ++values_at_gauss_pts;
434template <
int Tensor_Dim,
class T,
class L,
class A>
447 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
448 const EntityType zero_type = MBVERTEX)
467 boost::shared_ptr<ublas::matrix<T, L, A>>
dataPtr;
471template <
int Tensor_Dim,
class T,
class L,
class A>
477 "Not implemented for T = %s and dim = %d",
487template <
int Tensor_Dim>
493 boost::shared_ptr<MatrixDouble> data_ptr,
494 const EntityType zero_type = MBVERTEX)
503 boost::shared_ptr<MatrixDouble> data_ptr,
505 const EntityType zero_type = MBVERTEX)
527template <
int Tensor_Dim>
529 Tensor_Dim,
double, ublas::row_major,
534 const size_t nb_gauss_pts = getGaussPts().size2();
535 auto &mat = *dataPtr;
536 if (type == zeroType || mat.size2() != nb_gauss_pts) {
537 mat.resize(Tensor_Dim, nb_gauss_pts,
false);
544 if (dataVec.use_count()) {
545 dotVector.resize(nb_dofs,
false);
547 CHKERR VecGetArrayRead(dataVec, &array);
549 for (
int i = 0;
i != local_indices.size(); ++
i)
550 if (local_indices[
i] != -1)
551 dotVector[
i] = array[local_indices[
i]];
554 CHKERR VecRestoreArrayRead(dataVec, &array);
559 const size_t nb_base_functions = data.
getN().size2();
561 auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
563 const size_t size = nb_dofs / Tensor_Dim;
564 if (nb_dofs % Tensor_Dim) {
566 "Nb. of DOFs is inconsistent with Tensor_Dim");
568 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
572 if (field_data.l2() != field_data.l2()) {
573 MOFEM_LOG(
"SELF", Sev::error) <<
"field data: " << field_data;
575 "Wrong number in coefficients");
580 for (; bb != size; ++bb) {
584 if (base_function != base_function) {
585 MOFEM_LOG(
"SELF", Sev::error) <<
"base function: " << base_function;
587 "Wrong number number in base functions");
592 values_at_gauss_pts(
I) += field_data(
I) * base_function;
598 for (; bb < nb_base_functions; ++bb)
600 ++values_at_gauss_pts;
604 if (dataVec.use_count()) {
623template <
int Tensor_Dim>
626 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
629 Tensor_Dim,
double, ublas::row_major,
649template <
int Tensor_Dim, CoordinateTypes COORDINATE_SYSTEM = CARTESIAN>
661 const std::string
field_name, boost::shared_ptr<VectorDouble> data_ptr,
662 const EntityType zero_type = MBVERTEX)
675 if constexpr (COORDINATE_SYSTEM ==
POLAR || COORDINATE_SYSTEM ==
SPHERICAL)
677 "%s coordiante not implemented",
683 vec.resize(nb_gauss_pts,
false);
691 const size_t nb_base_functions = data.
getN().size2();
694 const size_t size = nb_dofs / Tensor_Dim;
696 if (nb_dofs % Tensor_Dim) {
698 "Number of dofs should multiple of dimensions");
703 if constexpr (COORDINATE_SYSTEM ==
CARTESIAN) {
705 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
708 for (; bb != size; ++bb) {
709 values_at_gauss_pts += field_data(
I) * diff_base_function(
I);
711 ++diff_base_function;
715 for (; bb < nb_base_functions; ++bb)
716 ++diff_base_function;
717 ++values_at_gauss_pts;
727 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
730 for (; bb != size; ++bb) {
731 values_at_gauss_pts += field_data(
I) * diff_base_function(
I);
732 values_at_gauss_pts +=
733 base_function * (field_data(0) / t_coords(0));
736 ++diff_base_function;
740 for (; bb < nb_base_functions; ++bb) {
742 ++diff_base_function;
744 ++values_at_gauss_pts;
765template <
int Tensor_Dim, PetscData::DataContext CTX>
770 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
771 const EntityType zero_at_type = MBVERTEX,
bool throw_error =
true)
784 const size_t nb_dofs = local_indices.size();
789 mat.resize(Tensor_Dim, nb_gauss_pts,
false);
803 auto get_array = [&](
const auto ctx,
auto vec) {
809 <<
"In this case field degrees of freedom are read from vector. "
810 "That usually happens when time solver is used, and access to "
811 "first or second rates is needed. You probably not set ts_u, "
812 "ts_u_t, or ts_u_tt and associated data structure, i.e. "
813 "data_ctx to CTX_SET_X, CTX_SET_DX, CTX_SET_X_T, or "
814 "CTX_SET_X_TT respectively";
818 CHKERR VecGetArrayRead(vec, &array);
822 auto restore_array = [&](
auto vec) {
823 return VecRestoreArrayRead(vec, &array);
841 "That case is not implemented");
845 for (
int i = 0;
i != local_indices.size(); ++
i)
846 if (local_indices[
i] != -1)
866 "That case is not implemented");
869 const size_t nb_base_functions = data.
getN().size2();
871 auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
874 const size_t size = nb_dofs / Tensor_Dim;
875 if (nb_dofs % Tensor_Dim) {
878 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
879 auto field_data = getFTensor1FromArray<Tensor_Dim, Tensor_Dim>(
dotVector);
881 for (; bb != size; ++bb) {
882 values_at_gauss_pts(
I) += field_data(
I) * base_function;
888 for (; bb < nb_base_functions; ++bb)
890 ++values_at_gauss_pts;
907template <
int Tensor_Dim>
917template <
int Tensor_Dim>
927template <
int Tensor_Dim>
941template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
947 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
948 const EntityType zero_type = MBVERTEX)
960 boost::shared_ptr<ublas::matrix<T, L, A>>
dataPtr;
964template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
970 "Not implemented for T = %s, dim0 = %d and dim1 = %d",
972 Tensor_Dim0, Tensor_Dim1);
976template <
int Tensor_Dim0,
int Tensor_Dim1>
984 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
986 const EntityType zero_type = MBVERTEX)
997 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
1011 boost::shared_ptr<ublas::matrix<double, ublas::row_major, DoubleAllocator>>
1018template <
int Tensor_Dim0,
int Tensor_Dim1>
1020 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1026 const size_t nb_gauss_pts = data.
getN().size1();
1027 if (type == zeroType) {
1028 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
1034 if (dataVec.use_count()) {
1035 dotVector.resize(nb_dofs,
false);
1036 const double *array;
1037 CHKERR VecGetArrayRead(dataVec, &array);
1039 for (
int i = 0;
i != local_indices.size(); ++
i)
1040 if (local_indices[
i] != -1)
1041 dotVector[
i] = array[local_indices[
i]];
1044 CHKERR VecRestoreArrayRead(dataVec, &array);
1049 const size_t nb_base_functions = data.
getN().size2();
1051 auto values_at_gauss_pts =
1052 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1055 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
1056 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1059 for (; bb != size; ++bb) {
1060 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
1064 for (; bb < nb_base_functions; ++bb)
1066 ++values_at_gauss_pts;
1069 if (dataVec.use_count()) {
1081template <
int Tensor_Dim0,
int Tensor_Dim1>
1084 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1087 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1096template <
int Tensor_Dim0,
int Tensor_Dim1>
1101 boost::shared_ptr<MatrixDouble> data_ptr,
1102 const EntityType zero_at_type = MBVERTEX)
1114 const size_t nb_gauss_pts =
getGaussPts().size2();
1117 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
1121 const size_t nb_dofs = local_indices.size();
1124 const double *array;
1126 for (
size_t i = 0;
i != local_indices.size(); ++
i)
1127 if (local_indices[
i] != -1)
1133 const size_t nb_base_functions = data.
getN().size2();
1136 auto values_at_gauss_pts =
1137 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1140 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
1141 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1143 getFTensor2FromPtr<Tensor_Dim0, Tensor_Dim1>(&*
dotVector.begin());
1145 for (; bb != size; ++bb) {
1146 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
1150 for (; bb < nb_base_functions; ++bb)
1152 ++values_at_gauss_pts;
1171template <
int Tensor_Dim>
1176 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1177 const EntityType zero_type = MBEDGE,
const int zero_side = 0)
1186 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1188 const int zero_side = 0)
1203 mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts,
false);
1212 const double *array;
1215 for (
int i = 0;
i != local_indices.size(); ++
i)
1216 if (local_indices[
i] != -1)
1224 const int nb_base_functions = data.
getN().size2();
1226 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1229 const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
1230 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1233 for (; bb != size; ++bb) {
1234 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
1238 for (; bb < nb_base_functions; ++bb)
1240 ++values_at_gauss_pts;
1265template <
int Tensor_Dim>
1270 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1271 const EntityType zero_type = MBEDGE,
const int zero_side = 0)
1284 constexpr auto symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
1286 mat.resize(symm_size, nb_gauss_pts,
false);
1290 const int nb_dofs = local_indices.size();
1298 <<
"In this case field degrees of freedom are read from vector. "
1299 "That usually happens when time solver is used, and acces to "
1300 "first rates is needed. You probably not set "
1301 "ts_u_t and associated data structure data_ctx to CTX_SET_X_T "
1308 const double *array;
1310 for (
int i = 0;
i != local_indices.size(); ++
i)
1311 if (local_indices[
i] != -1)
1317 const int nb_base_functions = data.
getN().size2();
1320 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1323 const int size = nb_dofs / symm_size;
1324 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1325 auto field_data = getFTensorDotData<Tensor_Dim>();
1327 for (; bb != size; ++bb) {
1328 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
1332 for (; bb < nb_base_functions; ++bb)
1334 ++values_at_gauss_pts;
1347 static_assert(Dim || !Dim,
"not implemented");
1356 &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3], &dotVector[4],
1365 &dotVector[0], &dotVector[1], &dotVector[2]);
1379template <
int Tensor_Dim,
class T,
class L,
class A>
1384 Tensor_Dim, T,
L,
A>::OpCalculateVectorFieldValues_General;
1391template <
int Tensor_Dim>
1395 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1398 Tensor_Dim,
double, ublas::row_major,
1417template <
int Tensor_Dim>
1419 Tensor_Dim,
double, ublas::row_major,
1424 const size_t nb_gauss_pts = this->getGaussPts().size2();
1425 auto &mat = *this->dataPtr;
1426 if (type == this->zeroType) {
1427 mat.resize(Tensor_Dim, nb_gauss_pts,
false);
1435 const int nb_base_functions = data.
getN().size2();
1437 auto gradients_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
1440 if (nb_dofs > nb_base_functions)
1442 "Number of base functions inconsistent with number of DOFs "
1444 nb_dofs, nb_base_functions);
1446 if (data.
getDiffN().size2() != nb_base_functions * Tensor_Dim)
1449 "Number of base functions inconsistent with number of derivatives "
1451 data.
getDiffN().size2(), nb_base_functions);
1453 if (data.
getDiffN().size1() != nb_gauss_pts)
1456 "Number of base functions inconsistent with number of integration "
1458 data.
getDiffN().size2(), nb_gauss_pts);
1463 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1466 for (; bb != nb_dofs; ++bb) {
1467 gradients_at_gauss_pts(
I) += field_data * diff_base_function(
I);
1469 ++diff_base_function;
1472 for (; bb < nb_base_functions; ++bb)
1473 ++diff_base_function;
1474 ++gradients_at_gauss_pts;
1487template <
int Tensor_Dim>
1490 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1492 Tensor_Dim,
double, ublas::row_major,
1500template <
int Tensor_Dim>
1503 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1506 Tensor_Dim,
double, ublas::row_major,
1520template <
int Tensor_Dim>
1525 const size_t nb_gauss_pts = this->getGaussPts().size2();
1526 auto &mat = *this->dataPtr;
1527 if (type == this->zeroType) {
1528 mat.resize(Tensor_Dim * Tensor_Dim, nb_gauss_pts,
false);
1536 const int nb_base_functions = data.
getN().size2();
1538 auto &hessian_base = data.
getN(BaseDerivatives::SecondDerivative);
1540 if (hessian_base.size1() != nb_gauss_pts) {
1542 "Wrong number of integration pts (%ld != %d)",
1543 hessian_base.size1(), nb_gauss_pts);
1545 if (hessian_base.size2() != nb_base_functions * Tensor_Dim * Tensor_Dim) {
1547 "Wrong number of base functions (%ld != %d)",
1548 hessian_base.size2() / (Tensor_Dim * Tensor_Dim),
1551 if (hessian_base.size2() < nb_dofs * Tensor_Dim * Tensor_Dim) {
1553 "Wrong number of base functions (%ld < %d)",
1554 hessian_base.size2(), nb_dofs * Tensor_Dim * Tensor_Dim);
1558 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim, Tensor_Dim>(
1559 &*hessian_base.data().begin());
1561 auto hessian_at_gauss_pts =
1562 getFTensor2FromMat<Tensor_Dim, Tensor_Dim>(mat);
1566 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1569 for (; bb != nb_dofs; ++bb) {
1570 hessian_at_gauss_pts(
I,
J) +=
1571 field_data * t_diff2_base_function(
I,
J);
1573 ++t_diff2_base_function;
1576 for (; bb < nb_base_functions; ++bb) {
1577 ++t_diff2_base_function;
1580 ++hessian_at_gauss_pts;
1606template <
int Tensor_Dim0,
int Tensor_Dim1,
int S,
class T,
class L,
class A>
1612 Tensor_Dim0, Tensor_Dim1, T,
L,
A>::OpCalculateTensor2FieldValues_General;
1615template <
int Tensor_Dim0,
int Tensor_Dim1,
int S>
1619 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1622 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1641template <
int Tensor_Dim0,
int Tensor_Dim1,
int S>
1643 Tensor_Dim0, Tensor_Dim1, S,
double, ublas::row_major,
1649 "Data pointer not allocated");
1651 const size_t nb_gauss_pts = this->getGaussPts().size2();
1652 auto &mat = *this->dataPtr;
1653 if (type == this->zeroType) {
1654 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
1663 if (this->dataVec.use_count()) {
1664 this->dotVector.resize(nb_dofs,
false);
1665 const double *array;
1666 CHKERR VecGetArrayRead(this->dataVec, &array);
1668 for (
int i = 0;
i != local_indices.size(); ++
i)
1669 if (local_indices[
i] != -1)
1670 this->dotVector[
i] = array[local_indices[
i]];
1672 this->dotVector[
i] = 0;
1673 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1677 const int nb_base_functions = data.
getN().size2();
1679 auto gradients_at_gauss_pts =
1680 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1683 int size = nb_dofs / Tensor_Dim0;
1684 if (nb_dofs % Tensor_Dim0) {
1686 "Data inconsistency");
1688 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1689 auto field_data = getFTensor1FromPtr<Tensor_Dim0, S>(
1693 if (field_data.l2() != field_data.l2()) {
1694 MOFEM_LOG(
"SELF", Sev::error) <<
"field data " << field_data;
1696 "Wrong number in coefficients");
1701 for (; bb < size; ++bb) {
1704 if (diff_base_function.l2() != diff_base_function.l2()) {
1706 <<
"diff_base_function: " << diff_base_function;
1708 "Wrong number number in base functions");
1713 gradients_at_gauss_pts(
I,
J) += field_data(
I) * diff_base_function(
J);
1715 ++diff_base_function;
1719 for (; bb != nb_base_functions; ++bb)
1720 ++diff_base_function;
1721 ++gradients_at_gauss_pts;
1724 if (this->dataVec.use_count()) {
1741template <
int Tensor_Dim0,
int Tensor_Dim1,
int S = Tensor_Dim0>
1744 double, ublas::row_major,
1748 Tensor_Dim0, Tensor_Dim1, S,
double, ublas::row_major,
1757template <
int Tensor_Dim0,
int Tensor_Dim1>
1762 boost::shared_ptr<MatrixDouble> data_ptr,
1763 const EntityType zero_at_type = MBVERTEX)
1776 const int nb_dofs = local_indices.size();
1777 const int nb_gauss_pts = this->
getGaussPts().size2();
1781 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts,
false);
1788 const double *array;
1790 for (
int i = 0;
i != local_indices.size(); ++
i)
1791 if (local_indices[
i] != -1)
1797 const int nb_base_functions = data.
getN().size2();
1799 auto gradients_at_gauss_pts =
1800 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1803 int size = nb_dofs / Tensor_Dim0;
1804 if (nb_dofs % Tensor_Dim0) {
1808 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1809 auto field_data = getFTensor1FromPtr<Tensor_Dim0>(&*
dotVector.begin());
1811 for (; bb < size; ++bb) {
1812 gradients_at_gauss_pts(
I,
J) += field_data(
I) * diff_base_function(
J);
1814 ++diff_base_function;
1818 for (; bb != nb_base_functions; ++bb)
1819 ++diff_base_function;
1820 ++gradients_at_gauss_pts;
1836template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
1842 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1843 const EntityType zero_type = MBVERTEX)
1849template <
int Tensor_Dim0,
int Tensor_Dim1>
1853 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1856 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1857 const EntityType zero_type = MBVERTEX)
1879template <
int Tensor_Dim0,
int Tensor_Dim1>
1881 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1887 "Data pointer not allocated");
1889 const size_t nb_gauss_pts = this->getGaussPts().size2();
1890 constexpr size_t msize = (Tensor_Dim0 * (Tensor_Dim0 + 1)) / 2;
1891 auto &mat = *this->dataPtr;
1892 if (type == this->zeroType) {
1893 mat.resize(msize * Tensor_Dim1, nb_gauss_pts,
false);
1902 const int nb_base_functions = data.
getN().size2();
1904 auto gradients_at_gauss_pts =
1905 getFTensor3DgFromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1909 int size = nb_dofs / msize;
1910 if (nb_dofs % msize) {
1912 "Data inconsistency");
1914 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
1917 for (; bb < size; ++bb) {
1918 gradients_at_gauss_pts(
I,
J,
K) +=
1919 field_data(
I,
J) * diff_base_function(
K);
1921 ++diff_base_function;
1925 for (; bb != nb_base_functions; ++bb)
1926 ++diff_base_function;
1927 ++gradients_at_gauss_pts;
1939template <
int Tensor_Dim0,
int Tensor_Dim1>
1942 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1945 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1946 const EntityType zero_type = MBVERTEX)
1948 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1952template <
int Tensor_Dim0,
int Tensor_Dim1>
1955 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1958 Tensor_Dim0, Tensor_Dim1,
double, ublas::row_major,
1977template <
int Tensor_Dim0,
int Tensor_Dim1>
1983 "Data pointer not allocated");
1985 const size_t nb_gauss_pts = this->getGaussPts().size2();
1986 auto &mat = *this->dataPtr;
1987 if (type == this->zeroType) {
1988 mat.resize(Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim1, nb_gauss_pts,
false);
1997 if (this->dataVec.use_count()) {
1998 this->dotVector.resize(nb_dofs,
false);
1999 const double *array;
2000 CHKERR VecGetArrayRead(this->dataVec, &array);
2002 for (
int i = 0;
i != local_indices.size(); ++
i)
2003 if (local_indices[
i] != -1)
2004 this->dotVector[
i] = array[local_indices[
i]];
2006 this->dotVector[
i] = 0;
2007 CHKERR VecRestoreArrayRead(this->dataVec, &array);
2011 const int nb_base_functions = data.
getN().size2();
2013 auto &hessian_base = data.
getN(BaseDerivatives::SecondDerivative);
2015 if (hessian_base.size1() != nb_gauss_pts) {
2017 "Wrong number of integration pts (%d != %d)",
2018 hessian_base.size1(), nb_gauss_pts);
2020 if (hessian_base.size2() !=
2021 nb_base_functions * Tensor_Dim1 * Tensor_Dim1) {
2023 "Wrong number of base functions (%d != %d)",
2024 hessian_base.size2() / (Tensor_Dim1 * Tensor_Dim1),
2027 if (hessian_base.size2() < nb_dofs * Tensor_Dim1 * Tensor_Dim1) {
2029 "Wrong number of base functions (%d < %d)",
2030 hessian_base.size2(), nb_dofs * Tensor_Dim1 * Tensor_Dim1);
2034 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim1, Tensor_Dim1>(
2035 &*hessian_base.data().begin());
2037 auto t_hessian_at_gauss_pts =
2038 getFTensor3FromMat<Tensor_Dim0, Tensor_Dim1, Tensor_Dim1>(mat);
2044 int size = nb_dofs / Tensor_Dim0;
2046 if (nb_dofs % Tensor_Dim0) {
2048 "Data inconsistency");
2052 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
2055 for (; bb < size; ++bb) {
2056 t_hessian_at_gauss_pts(
I,
J,
K) +=
2057 field_data(
I) * t_diff2_base_function(
J,
K);
2059 ++t_diff2_base_function;
2063 for (; bb != nb_base_functions; ++bb)
2064 ++t_diff2_base_function;
2065 ++t_hessian_at_gauss_pts;
2068 if (this->dataVec.use_count()) {
2090template <
int DIM_01,
int DIM_23,
int S = 0>
2102 boost::shared_ptr<MatrixDouble> in_mat,
2103 boost::shared_ptr<MatrixDouble> out_mat,
2104 boost::shared_ptr<MatrixDouble> d_mat)
2117 boost::shared_ptr<MatrixDouble> out_mat,
2118 boost::shared_ptr<MatrixDouble> d_mat)
2131 const size_t nb_gauss_pts =
getGaussPts().size2();
2132 auto t_D = getFTensor4DdgFromMat<DIM_01, DIM_23, S>(*(
dMat));
2133 auto t_in = getFTensor2SymmetricFromMat<DIM_01>(*(
inMat));
2134 outMat->resize((DIM_23 * (DIM_23 + 1)) / 2, nb_gauss_pts,
false);
2135 auto t_out = getFTensor2SymmetricFromMat<DIM_23>(*(
outMat));
2136 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
2137 t_out(
i,
j) = t_D(
i,
j,
k,
l) * t_in(
k,
l);
2152 boost::shared_ptr<MatrixDouble>
dMat;
2175 boost::shared_ptr<MatrixDouble> in_mat,
2176 boost::shared_ptr<MatrixDouble> out_mat)
2193 boost::shared_ptr<MatrixDouble> out_mat)
2204 const size_t nb_gauss_pts =
getGaussPts().size2();
2205 auto t_in = getFTensor2FromMat<DIM, DIM>(*(
inMat));
2206 outMat->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts,
false);
2207 auto t_out = getFTensor2SymmetricFromMat<DIM>(*(
outMat));
2208 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
2209 t_out(
i,
j) = (t_in(
i,
j) || t_in(
j,
i)) / 2;
2240 boost::shared_ptr<MatrixDouble> in_mat,
2241 boost::shared_ptr<MatrixDouble> out_mat)
2260 boost::shared_ptr<MatrixDouble> in_mat,
2261 boost::shared_ptr<MatrixDouble> out_mat)
2275 noalias(*
outMat) = (*scalePtr) * (*inMat);
2294template <
int Base_Dim,
int Field_Dim,
class T,
class L,
class A>
2300template <
int Field_Dim>
2306 boost::shared_ptr<MatrixDouble> data_ptr,
2308 const EntityType zero_type = MBEDGE,
2309 const int zero_side = 0)
2312 dataPtr(data_ptr), dataVec(data_vec), zeroType(zero_type),
2313 zeroSide(zero_side) {
2319 boost::shared_ptr<MatrixDouble> data_ptr,
2320 const EntityType zero_type = MBEDGE,
2321 const int zero_side = 0)
2343template <
int Field_Dim>
2345 3, Field_Dim,
double, ublas::row_major,
2349 const size_t nb_integration_points = this->getGaussPts().size2();
2350 if (type == zeroType && side == zeroSide) {
2351 dataPtr->resize(Field_Dim, nb_integration_points,
false);
2358 if (dataVec.use_count()) {
2359 dotVector.resize(nb_dofs,
false);
2360 const double *array;
2361 CHKERR VecGetArrayRead(dataVec, &array);
2363 for (
int i = 0;
i != local_indices.size(); ++
i)
2364 if (local_indices[
i] != -1)
2365 dotVector[
i] = array[local_indices[
i]];
2368 CHKERR VecRestoreArrayRead(dataVec, &array);
2372 const size_t nb_base_functions = data.
getN().size2() / 3;
2375 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2376 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2379 for (; bb != nb_dofs; ++bb) {
2380 t_data(
i) += t_n_hdiv(
i) * t_dof;
2384 for (; bb != nb_base_functions; ++bb)
2389 if (dataVec.use_count()) {
2399template <
int Base_Dim,
int Field_Dim = Base_Dim>
2402 Base_Dim, Field_Dim, double, ublas::row_major, DoubleAllocator> {
2404 Base_Dim, Field_Dim,
double, ublas::row_major,
2411template <
int Base_Dim,
int Field_Dim = Base_Dim>
2414template <
int Field_Dim>
2419 boost::shared_ptr<MatrixDouble> data_ptr,
2420 const EntityType zero_type = MBEDGE,
2421 const int zero_side = 0)
2424 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2445template <
int Field_Dim>
2450 const size_t nb_integration_points = this->getGaussPts().size2();
2451 if (type == zeroType && side == zeroSide) {
2452 dataPtr->resize(Field_Dim, nb_integration_points,
false);
2457 const size_t nb_dofs = local_indices.size();
2460 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2461 const double *array;
2462 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2463 for (
size_t i = 0;
i != nb_dofs; ++
i)
2464 if (local_indices[
i] != -1)
2465 dot_dofs_vector[
i] = array[local_indices[
i]];
2467 dot_dofs_vector[
i] = 0;
2468 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2470 const size_t nb_base_functions = data.
getN().size2() / 3;
2473 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2474 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2476 for (; bb != nb_dofs; ++bb) {
2477 t_data(
i) += t_n_hdiv(
i) * dot_dofs_vector[bb];
2480 for (; bb != nb_base_functions; ++bb)
2496template <
int BASE_DIM,
int SPACE_DIM>
2501 boost::shared_ptr<VectorDouble> data_ptr,
2502 const EntityType zero_type = MBEDGE,
2503 const int zero_side = 0)
2514 const size_t nb_integration_points =
getGaussPts().size2();
2516 dataPtr->resize(nb_integration_points,
false);
2522 const size_t nb_base_functions = data.
getN().size2() /
BASE_DIM;
2527 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2530 for (; bb != nb_dofs; ++bb) {
2531 t_data += t_dof * t_n_diff_hdiv(
j,
j);
2535 for (; bb != nb_base_functions; ++bb)
2555template <
int BASE_DIM,
int SPACE_DIM>
2560 boost::shared_ptr<MatrixDouble> data_ptr,
2561 const EntityType zero_type = MBEDGE,
2562 const int zero_side = 0)
2573 const size_t nb_integration_points =
getGaussPts().size2();
2581 const size_t nb_base_functions = data.
getN().size2() /
BASE_DIM;
2585 auto t_data = getFTensor2FromMat<BASE_DIM, SPACE_DIM>(*
dataPtr);
2586 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2589 for (; bb != nb_dofs; ++bb) {
2590 t_data(
i,
j) += t_dof * t_base_diff(
i,
j);
2594 for (; bb != nb_base_functions; ++bb)
2615template <
int BASE_DIM,
int FIELD_DIM,
int SPACE_DIM>
2628 boost::shared_ptr<MatrixDouble> data_ptr,
2629 const EntityType zero_type = MBEDGE,
2630 const int zero_side = 0)
2633 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2642 const size_t nb_integration_points = getGaussPts().size2();
2643 if (type == zeroType && side == zeroSide) {
2644 dataPtr->resize(27, nb_integration_points,
2653 "Data inconsistency, nb_dofs %% COEFF_DIM != 0, that is %ld %% %d "
2663 const size_t nb_base_functions = data.
getN().size2() / 3;
2669 auto t_data = getFTensor3FromMat<3, 3, 3>(*dataPtr);
2670 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2673 for (; bb != nb_dofs; ++bb) {
2674 t_data(
k,
i,
j) += t_base_diff(
i,
j) * t_dof(
k);
2678 for (; bb != nb_base_functions; ++bb)
2696 boost::shared_ptr<MatrixDouble> data_ptr,
2697 const EntityType zero_type = MBEDGE,
2698 const int zero_side = 0)
2701 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2710 const size_t nb_integration_points = getGaussPts().size2();
2711 if (type == zeroType && side == zeroSide) {
2712 dataPtr->resize(27, nb_integration_points,
false);
2720 const size_t nb_base_functions = data.
getN().size2() / 9;
2723 if (data.
getDiffN().size1() != nb_integration_points) {
2725 "Wrong number of integration pts (%ld != %ld)",
2726 data.
getDiffN().size1(), nb_integration_points);
2728 if (data.
getDiffN().size2() != nb_base_functions * 27) {
2730 "Wrong number of base functions (%ld != %ld)",
2731 data.
getDiffN().size2() / 27, nb_base_functions);
2733 if (nb_base_functions < nb_dofs) {
2735 "Wrong number of base functions (%ld < %ld)", nb_base_functions,
2745 auto t_data = getFTensor3FromMat<3, 3, 3>(*dataPtr);
2746 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2749 for (; bb != nb_dofs; ++bb) {
2750 t_data(
k,
i,
j) += t_base_diff(
k,
i,
j) * t_dof;
2754 for (; bb != nb_base_functions; ++bb)
2776template <
int BASE_DIM,
int SPACE_DIM>
2781 boost::shared_ptr<MatrixDouble> data_ptr,
2782 const EntityType zero_type = MBEDGE,
2783 const int zero_side = 0)
2794 const size_t nb_integration_points =
getGaussPts().size2();
2804 const int nb_base_functions = data.
getN().size2() /
BASE_DIM;
2807 auto &hessian_base = data.
getN(BaseDerivatives::SecondDerivative);
2808 if (hessian_base.size1() != nb_integration_points) {
2810 "Wrong number of integration pts (%d != %d)",
2811 hessian_base.size1(), nb_integration_points);
2813 if (hessian_base.size2() !=
2816 "Wrong number of base functions (%d != %d)",
2822 "Wrong number of base functions (%d < %d)", hessian_base.size2(),
2833 auto t_data = getFTensor3FromMat<BASE_DIM, SPACE_DIM, SPACE_DIM>(*
dataPtr);
2834 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2837 for (; bb != nb_dofs; ++bb) {
2838 t_data(
i,
j,
k) += t_dof * t_base_diff2(
i,
j,
k);
2843 for (; bb != nb_base_functions; ++bb)
2862template <
int Tensor_Dim1,
int Tensor_Dim2>
2867 boost::shared_ptr<VectorDouble> data_ptr,
2868 const EntityType zero_type = MBEDGE,
2869 const int zero_side = 0)
2880 const size_t nb_integration_points =
getGaussPts().size2();
2882 dataPtr->resize(nb_integration_points,
false);
2887 const int nb_dofs = local_indices.size();
2890 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2891 const double *array;
2893 for (
size_t i = 0;
i != local_indices.size(); ++
i)
2894 if (local_indices[
i] != -1)
2895 dot_dofs_vector[
i] = array[local_indices[
i]];
2897 dot_dofs_vector[
i] = 0;
2900 const size_t nb_base_functions = data.
getN().size2() / Tensor_Dim1;
2904 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
2906 for (; bb != nb_dofs; ++bb) {
2908 for (
auto ii = 0; ii != Tensor_Dim2; ++ii)
2909 div += t_n_diff_hdiv(ii, ii);
2910 t_data += dot_dofs_vector[bb] * div;
2913 for (; bb != nb_base_functions; ++bb)
2949 boost::shared_ptr<MatrixDouble> data_ptr,
2950 const EntityType zero_type = MBEDGE,
2951 const int zero_side = 0);
2973 boost::shared_ptr<MatrixDouble> data_ptr,
2974 const EntityType zero_type = MBVERTEX,
2975 const int zero_side = 0);
2998 boost::shared_ptr<MatrixDouble> data_ptr,
2999 const EntityType zero_type = MBVERTEX,
3000 const int zero_side = 0);
3018template <
int Tensor_Dim0,
int Tensor_Dim1>
3023 boost::shared_ptr<MatrixDouble> data_ptr,
3024 boost::shared_ptr<double> scale_ptr,
3026 const EntityType zero_type = MBEDGE,
3027 const int zero_side = 0)
3037 boost::shared_ptr<MatrixDouble> data_ptr,
3038 const EntityType zero_type = MBEDGE,
3039 const int zero_side = 0)
3047 const size_t nb_integration_points =
getGaussPts().size2();
3049 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points,
false);
3057 const double *array;
3060 for (
int i = 0;
i != local_indices.size(); ++
i)
3061 if (local_indices[
i] != -1)
3070 const size_t nb_base_functions = data.
getN().size2() / 3;
3074 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*
dataPtr);
3075 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
3078 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
3079 t_data(
i,
j) += (
scale * t_dof(
i)) * t_n_hvec(
j);
3083 for (; bb < nb_base_functions; ++bb)
3109template <
int Tensor_Dim0,
int Tensor_Dim1>
3116 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
3118 boost::shared_ptr<Range> broken_range_ptr =
nullptr,
3119 boost::shared_ptr<double> scale_ptr =
nullptr,
3120 const EntityType zero_type = MBEDGE,
const int zero_side = 0)
3148template <
int Tensor_Dim0,
int Tensor_Dim1>
3153 const size_t nb_integration_points = OP::getGaussPts().size2();
3154 if (type == zeroType) {
3155 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points,
false);
3162 if (dataVec.use_count()) {
3163 dotVector.resize(nb_dofs,
false);
3164 const double *array;
3165 CHKERR VecGetArrayRead(dataVec, &array);
3167 for (
int i = 0;
i != local_indices.size(); ++
i)
3168 if (local_indices[
i] != -1)
3169 dotVector[
i] = array[local_indices[
i]];
3172 CHKERR VecRestoreArrayRead(dataVec, &array);
3183 auto get_get_side_face_dofs = [&]() {
3184 auto fe_type = OP::getFEType();
3188 std::vector<int> side_face_dofs;
3189 side_face_dofs.reserve(data.
getIndices().size() / Tensor_Dim0);
3193 auto it = side_dof_map.get<1>().begin();
3194 it != side_dof_map.get<1>().end(); ++it
3197 if ((Tensor_Dim0 * it->dof) >= data.
getIndices().size()) {
3200 if (it->type == brokenType) {
3201 if (brokenRangePtr) {
3202 auto ent = OP::getSideEntity(it->side, brokenType);
3203 if (brokenRangePtr->find(ent) != brokenRangePtr->end()) {
3204 side_face_dofs.push_back(it->dof);
3207 side_face_dofs.push_back(it->dof);
3212 return side_face_dofs;
3215 auto side_face_dofs = get_get_side_face_dofs();
3219 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
3220 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
3221 for (
auto b : side_face_dofs) {
3223 auto t_dof = getFTensor1FromPtr<Tensor_Dim0>(data.
getFieldData().data() +
3225 t_data(
i,
j) += t_dof(
i) * t_row_base(
j);
3229 *dataPtr *= (scalePtr) ? *scalePtr : 1.0;
3231 if (dataVec.use_count()) {
3245template <
int Tensor_Dim0,
int Tensor_Dim1>
3250 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
3251 boost::shared_ptr<double> scale_ptr,
3253 const EntityType zero_type = MBEDGE,
const int zero_side = 0)
3263 boost::shared_ptr<MatrixDouble> data_ptr,
3264 const EntityType zero_type = MBEDGE,
3265 const int zero_side = 0)
3273 const size_t nb_integration_points =
getGaussPts().size2();
3275 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points,
false);
3284 const double *array;
3287 for (
int i = 0;
i != local_indices.size(); ++
i)
3288 if (local_indices[
i] != -1)
3297 const size_t nb_base_functions =
3298 data.
getN().size2() / (Tensor_Dim0 * Tensor_Dim1);
3301 auto t_n_hten = data.
getFTensor2N<Tensor_Dim0, Tensor_Dim1>();
3302 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*
dataPtr);
3303 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
3306 for (; bb != nb_dofs; ++bb) {
3307 t_data(
i,
j) += (
scale * t_dof) * t_n_hten(
i,
j);
3311 for (; bb < nb_base_functions; ++bb)
3339template <
int Tensor_Dim0,
int Tensor_Dim1,
3345 boost::shared_ptr<MatrixDouble> data_ptr,
3347 const EntityType zero_type = MBEDGE,
3348 const int zero_side = 0)
3358 boost::shared_ptr<MatrixDouble> data_ptr,
3359 const EntityType zero_type = MBEDGE,
3360 const int zero_side = 0)
3367 const size_t nb_integration_points =
getGaussPts().size2();
3369 dataPtr->resize(Tensor_Dim0, nb_integration_points,
false);
3377 const double *array;
3380 for (
int i = 0;
i != local_indices.size(); ++
i)
3381 if (local_indices[
i] != -1)
3389 const size_t nb_base_functions = data.
getN().size2() / 3;
3393 auto t_data = getFTensor1FromMat<Tensor_Dim0>(*
dataPtr);
3396 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
3399 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
3400 double div = t_n_diff_hvec(
j,
j);
3401 t_data(
i) += t_dof(
i) * div;
3403 t_data(
i) += t_base(0) * (t_dof(
i) / t_coords(0));
3409 for (; bb < nb_base_functions; ++bb) {
3442template <
int Tensor_Dim0,
int Tensor_Dim1,
3450 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
3452 boost::shared_ptr<Range> broken_range_ptr =
nullptr,
3453 boost::shared_ptr<double> scale_ptr =
nullptr,
3454 const EntityType zero_type = MBEDGE)
3466 const size_t nb_integration_points =
getGaussPts().size2();
3467 if (type ==
zeroType && side == 0) {
3468 dataPtr->resize(Tensor_Dim0, nb_integration_points,
false);
3477 const double *array;
3480 for (
int i = 0;
i != local_indices.size(); ++
i)
3481 if (local_indices[
i] != -1)
3496 auto get_get_side_face_dofs = [&]() {
3501 std::vector<int> side_face_dofs;
3502 side_face_dofs.reserve(data.
getIndices().size() / Tensor_Dim0);
3506 auto it = side_dof_map.get<1>().begin();
3507 it != side_dof_map.get<1>().end(); ++it
3510 if ((Tensor_Dim0 * it->dof) >= data.
getIndices().size()) {
3517 side_face_dofs.push_back(it->dof);
3520 side_face_dofs.push_back(it->dof);
3525 return side_face_dofs;
3528 auto side_face_dofs = get_get_side_face_dofs();
3532 auto t_data = getFTensor1FromMat<Tensor_Dim0>(*
dataPtr);
3534 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
3535 for (
auto b : side_face_dofs) {
3536 auto t_dof = getFTensor1FromPtr<Tensor_Dim0>(
3540 double div = t_diff_base(
j,
j);
3541 t_data(
i) += t_dof(
i) * div;
3543 t_data(
i) += t_base(0) * (t_dof(
i) / t_coords(0));
3574template <
int Tensor_Dim,
typename OpBase>
3578 boost::shared_ptr<MatrixDouble> data_ptr,
3579 boost::shared_ptr<double> scale_ptr,
3580 const EntityType zero_type = MBEDGE,
3581 const int zero_side = 0)
3589 boost::shared_ptr<MatrixDouble> data_ptr,
3590 const EntityType zero_type = MBEDGE,
3591 const int zero_side = 0)
3598 const size_t nb_integration_points = OpBase::getGaussPts().size2();
3599 if (type ==
zeroType && side == 0) {
3600 dataPtr->resize(Tensor_Dim, nb_integration_points,
false);
3606 auto t_normal = OpBase::getFTensor1NormalsAtGaussPts();
3607 const size_t nb_base_functions = data.
getN().size2() / 3;
3609 auto t_data = getFTensor1FromMat<Tensor_Dim>(*
dataPtr);
3610 for (
size_t gg = 0; gg != nb_integration_points; ++gg) {
3612 t_normalized_normal(
j) = t_normal(
j);
3616 for (; bb != nb_dofs / Tensor_Dim; ++bb) {
3618 (scale_val * t_dof(
i)) * (t_base(
j) * t_normalized_normal(
j));
3622 for (; bb < nb_base_functions; ++bb) {
3665 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
3679 template <
int D1,
int D2,
int J1,
int J2>
3683 static_assert(D2 == J2,
"Dimension of jacobian and dimension of <out> "
3684 "directive does not match");
3686 size_t nb_functions = diff_n.size2() / D1;
3688 size_t nb_gauss_pts = diff_n.size1();
3693 "Wrong number of Gauss Pts");
3694 if (diff_n.size2() % D1)
3696 "Number of directives of base functions and D1 dimension does "
3700 diffNinvJac.resize(diff_n.size1(), D2 * nb_functions,
false);
3704 auto t_diff_n = getFTensor1FromPtr<D2>(&*
diffNinvJac.data().begin());
3705 auto t_diff_n_ref = getFTensor1FromPtr<D1>(&*diff_n.data().begin());
3706 auto t_inv_jac = getFTensor2FromMat<J1, J2>(*
invJacPtr);
3707 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
3708 for (
size_t dd = 0; dd != nb_functions; ++dd) {
3709 t_diff_n(
i) = t_inv_jac(
k,
i) * t_diff_n_ref(
k);
3754template <
int DERIVARIVE = 1>
3761template <
int DERIVARIVE = 1>
3768template <
int DERIVARIVE = 1>
3772 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3776template <
int DERIVARIVE = 1>
3780 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3798 invJacPtr(inv_jac_ptr) {}
3852 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
3889 boost::shared_ptr<MatrixDouble> jac_ptr)
3906 2>::OpSetContravariantPiolaTransformOnFace2DImpl;
4106 boost::shared_ptr<VectorDouble> det_ptr,
4107 boost::shared_ptr<MatrixDouble> out_ptr)
4127 "Pointer for inPtr matrix not allocated");
4130 "Pointer for detPtr matrix not allocated");
4132 const auto nb_rows = inPtr->size1();
4133 const auto nb_integration_pts = inPtr->size2();
4137 detPtr->resize(nb_integration_pts,
false);
4138 auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
4140 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
4149 outPtr->resize(nb_rows, nb_integration_pts,
false);
4150 auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
4151 auto t_out = getFTensor2FromMat<DIM, DIM>(*outPtr);
4153 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
4188 boost::shared_ptr<VectorDouble> out_ptr)
4209 "Pointer for inPtr matrix not allocated");
4211 const auto nb_integration_pts = inPtr->size2();
4214 outPtr->resize(nb_integration_pts,
false);
4215 auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
4218 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
4241 boost::shared_ptr<VectorDouble> out_ptr)
4262 "Pointer for inPtr matrix not allocated");
4264 const auto nb_integration_pts = inPtr->size2();
4267 outPtr->resize(nb_integration_pts,
false);
4268 auto t_in = getFTensor2SymmetricFromMat<DIM>(*inPtr);
4271 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 CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#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', DIM1 > J
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VecAllocator< double > DoubleAllocator
implementation of Data Operators for Forces and Sources
static MoFEMErrorCode invertTensor(FTensor::Tensor2< T1, DIM, DIM > &t, T2 &det, FTensor::Tensor2< T3, DIM, DIM > &inv_t)
static auto determinantTensor(FTensor::Tensor2< T, DIM, DIM > &t)
Calculate the determinant of a tensor of rank DIM.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
constexpr IntegrationType I
constexpr auto field_name
multi_index_container< DofsSideMapData, indexed_by< ordered_non_unique< tag< TypeSide_mi_tag >, composite_key< DofsSideMapData, member< DofsSideMapData, EntityType, &DofsSideMapData::type >, member< DofsSideMapData, int, &DofsSideMapData::side > > >, ordered_unique< tag< EntDofIdx_mi_tag >, member< DofsSideMapData, int, &DofsSideMapData::dof > > > > DofsSideMap
Map entity stype and side to element/entity dof index.
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
const VectorFieldEntities & getFieldEntities() const
Get field entities (const version)
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 DOF values on entity.
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 degrees of freedom on entity.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from field data coefficients.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Return 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 field data coefficients.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2FieldData()
Return FTensor rank 2, i.e. matrix from field data coefficients.
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
FTensor::Tensor3< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 *Tensor_Dim2 >, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > getFTensor3DiffN(FieldApproximationBase base)
Get derivatives of base functions for tonsorial Hdiv space.
default operator for TRI element
default operator for Flat Prism element
default operator for Flat Prism element
FlatPrism finite element.
EntityType getFEType() const
Get dimension of finite element.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
EntityHandle getSideEntity(const int side_number, const EntityType type)
Get the side entity.
@ OPROW
operator doWork function is executed on FE rows
@ OPSPACE
operator do Work is execute on space data
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
structure to get information from mofem into EntitiesFieldData
Calculate divergence of tonsorial field using vectorial base.
boost::shared_ptr< double > scalePtr
boost::shared_ptr< Range > brokenRangePtr
SmartPetscObj< Vec > dataVec
OpCalculateBrokenHVecTensorDivergence(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, SmartPetscObj< Vec > data_vec, EntityType broken_type, boost::shared_ptr< Range > broken_range_ptr=nullptr, boost::shared_ptr< double > scale_ptr=nullptr, const EntityType zero_type=MBEDGE)
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
Get tensor field for H-div approximation.
boost::shared_ptr< MatrixDouble > dataPtr
boost::shared_ptr< Range > brokenRangePtr
boost::shared_ptr< double > scalePtr
OpCalculateBrokenHVecTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, SmartPetscObj< Vec > data_vec, EntityType broken_type, boost::shared_ptr< Range > broken_range_ptr=nullptr, boost::shared_ptr< double > scale_ptr=nullptr, 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
SmartPetscObj< Vec > dataVec
Calculate divergence of vector field at integration points.
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)
Constructor for vector field divergence calculation operator.
const EntityHandle zeroType
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
VectorDouble dotVector
Keeps temporary values of time derivatives.
OpCalculateHTensorTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, boost::shared_ptr< double > scale_ptr, SmartPetscObj< Vec > data_vec=SmartPetscObj< Vec >(), const EntityType zero_type=MBEDGE, const int zero_side=0)
SmartPetscObj< Vec > dataVec
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)
boost::shared_ptr< double > scalePtr
Calculate divergence of tonsorial field using vectorial base.
boost::shared_ptr< MatrixDouble > dataPtr
SmartPetscObj< Vec > dataVec
VectorDouble dotVector
Keeps temporary values of time derivatives.
const EntityHandle zeroType
OpCalculateHVecTensorDivergence(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)
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.
OpCalculateHVecTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, boost::shared_ptr< double > scale_ptr, SmartPetscObj< Vec > data_vec=SmartPetscObj< Vec >(), const EntityType zero_type=MBEDGE, const int zero_side=0)
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)
SmartPetscObj< Vec > dataVec
VectorDouble dotVector
Keeps temporary values of time derivatives.
boost::shared_ptr< double > scalePtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
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
OpCalculateHVecTensorGradient(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
const EntityHandle zeroType
OpCalculateHVecTensorGradient(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
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate gradient of tensor field.
Calculate trace of vector (Hdiv/Hcurl) space.
const EntityHandle zeroType
boost::shared_ptr< MatrixDouble > dataPtr
boost::shared_ptr< double > scalePtr
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
OpCalculateHVecTensorTrace(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, boost::shared_ptr< double > scale_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
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.
OpCalculateHVecVectorField_General(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)
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)
SmartPetscObj< Vec > dataVec
const EntityHandle zeroType
Get vector field for H-div approximation.
Get vector field for H-div approximation.
Calculate gradient of vector field.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
const EntityHandle zeroType
OpCalculateHVecVectorGradient(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate gradient of vector field.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
const EntityHandle zeroType
OpCalculateHVecVectorHessian(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
boost::shared_ptr< MatrixDouble > dataPtr
const EntityHandle zeroType
const EntityHandle zeroType
boost::shared_ptr< MatrixDouble > dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
const EntityHandle zeroType
Calculate curl of vector field.
Calculate divergence of vector field dot.
boost::shared_ptr< VectorDouble > dataPtr
OpCalculateHdivVectorDivergenceDot(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
const EntityHandle zeroType
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate divergence of vector field.
boost::shared_ptr< VectorDouble > dataPtr
OpCalculateHdivVectorDivergence(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
const EntityHandle zeroType
Calculate inverse of jacobian for face element.
OpCalculateInvJacForFatPrism(MatrixDouble &inv_jac)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateInvJacForFatPrism(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
const boost::shared_ptr< MatrixDouble > invJacPtr
Calculate inverse of jacobian for face element.
OpCalculateInvJacForFlatPrism(MatrixDouble &inv_jac_f3)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Evaluate field gradient values for scalar field, i.e. gradient is tensor rank 1 (vector)
Get field gradients at integration pts for scalar field 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
Calculate scalar field values from PETSc vector 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)
Constructor for PETSc vector-based scalar field calculation.
const EntityHandle zeroAtType
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)
Constructor with PETSc vector support.
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)
Constructor for scalar field values calculation operator.
Specialization for double precision scalar field values calculation.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of scalar field at integration points
Get time direvarive values at integration pts for tensor field 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 temporary values of time derivatives.
boost::shared_ptr< MatrixDouble > dataPtr
Data computed into this matrix.
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 values at integration pts for tensor field 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
Operator for calculating the trace of matrices at integration points.
boost::shared_ptr< MatrixDouble > inPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Index< 'i', DIM > i
boost::shared_ptr< VectorDouble > outPtr
OpCalculateTraceFromMat(boost::shared_ptr< MatrixDouble > in_ptr, boost::shared_ptr< VectorDouble > out_ptr)
Constructor for matrix trace calculation operator.
Calculates the trace of an input matrix.
FTensor::Index< 'i', DIM > i
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 > outPtr
boost::shared_ptr< MatrixDouble > inPtr
OpCalculateTraceFromSymmMat(boost::shared_ptr< MatrixDouble > in_ptr, boost::shared_ptr< VectorDouble > out_ptr)
Get field gradients time derivative at integration pts for scalar field 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 temporary values of time derivatives.
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 vector field, i.e. gradient is tensor rank 2.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of vector field at integration points
Approximate field values for given petsc vector.
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateVectorFieldValuesFromPetscVecImpl(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_type=MBVERTEX, bool throw_error=true)
const EntityHandle zeroAtType
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
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 tensor 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)
Constructor for vector field values calculation operator.
const EntityHandle zeroType
boost::shared_ptr< ublas::matrix< T, L, A > > dataPtr
Specialization for double precision vector field values calculation.
Operator for inverting matrices at integration points.
boost::shared_ptr< VectorDouble > detPtr
OpInvertMatrix(boost::shared_ptr< MatrixDouble > in_ptr, boost::shared_ptr< VectorDouble > det_ptr, boost::shared_ptr< MatrixDouble > out_ptr)
Constructor for matrix inversion operator.
boost::shared_ptr< MatrixDouble > outPtr
boost::shared_ptr< MatrixDouble > inPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Make Hdiv space from Hcurl space in 2d.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Operator for fat prism element updating integration weights in the volume.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpMultiplyDeterminantOfJacobianAndWeightsForFatPrisms()
Operator for scaling matrix values by a scalar factor.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > outMat
DEPRECATED OpScaleMatrix(const std::string field_name, const double scale, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
boost::shared_ptr< double > scalePtr
OpScaleMatrix(boost::shared_ptr< double > scale_ptr, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
Constructor for matrix scaling operator.
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)
Apply transformation to the input matrix.
OpSetInvJacToScalarBasesBasic(FieldSpace space, boost::shared_ptr< MatrixDouble > inv_jac_ptr)
boost::shared_ptr< MatrixDouble > invJacPtr
Operator for symmetrizing tensor fields.
FTensor::Index< 'i', DIM > i
boost::shared_ptr< MatrixDouble > outMat
OpSymmetrizeTensor(boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
Constructor for tensor symmetrization operator.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Index< 'j', DIM > j
DEPRECATED OpSymmetrizeTensor(const std::string field_name, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
boost::shared_ptr< MatrixDouble > inMat
FTensor::Index< 'i', DIM_01 > i
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Index< 'k', DIM_23 > k
boost::shared_ptr< MatrixDouble > dMat
OpTensorTimesSymmetricTensor(boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat, boost::shared_ptr< MatrixDouble > d_mat)
FTensor::Index< 'j', DIM_01 > j
boost::shared_ptr< MatrixDouble > inMat
FTensor::Index< 'l', DIM_23 > l
DEPRECATED OpTensorTimesSymmetricTensor(const std::string field_name, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat, boost::shared_ptr< MatrixDouble > d_mat)
boost::shared_ptr< MatrixDouble > outMat
static constexpr Switches CtxSetX
Solution vector switch.
static constexpr Switches CtxSetX_TT
Second time derivative switch.
std::bitset< 8 > Switches
Bitset type for context switches.
static constexpr Switches CtxSetX_T
First time derivative switch.
@ CTX_SET_X_T
Time derivative X_t is set.
@ CTX_SET_DX
Solution increment DX is set.
@ CTX_SET_X
Solution vector X is set.
@ CTX_SET_X_TT
Second time derivative X_tt is set.
intrusive_ptr for managing petsc objects