14#include <adolc/adolc.h>
21#error "MoFEM need to be compiled with ADOL-C"
26 meshPositionsFieldName =
"NoNE";
28 auto create_vec = [&]() {
29 constexpr int ghosts[] = {0};
45 CHKERR VolumeElementForcesAndSourcesCore::preProcess();
61 CHKERR VolumeElementForcesAndSourcesCore::postProcess();
66 CHKERR VecAssemblyBegin(V);
85 std::vector<VectorDouble> &values_at_gauss_pts,
86 std::vector<MatrixDouble> &gardient_at_gauss_pts)
89 valuesAtGaussPts(values_at_gauss_pts),
90 gradientAtGaussPts(gardient_at_gauss_pts), zeroAtType(MBVERTEX) {}
100 int nb_gauss_pts = data.
getN().size1();
101 int nb_base_functions = data.
getN().size2();
105 valuesAtGaussPts.resize(nb_gauss_pts);
106 gradientAtGaussPts.resize(nb_gauss_pts);
107 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
108 valuesAtGaussPts[gg].resize(3);
109 gradientAtGaussPts[gg].resize(3, 3);
112 if (
type == zeroAtType) {
113 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
114 valuesAtGaussPts[gg].clear();
115 gradientAtGaussPts[gg].clear();
124 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
127 &valuesAtGaussPts[gg][1],
128 &valuesAtGaussPts[gg][2]);
130 &gradientAtGaussPts[gg](0, 0), &gradientAtGaussPts[gg](0, 1),
131 &gradientAtGaussPts[gg](0, 2), &gradientAtGaussPts[gg](1, 0),
132 &gradientAtGaussPts[gg](1, 1), &gradientAtGaussPts[gg](1, 2),
133 &gradientAtGaussPts[gg](2, 0), &gradientAtGaussPts[gg](2, 1),
134 &gradientAtGaussPts[gg](2, 2));
136 for (; bb != nb_dofs / 3; bb++) {
137 values(
i) += base_function * field_data(
i);
138 gradient(
i,
j) += field_data(
i) * diff_base_functions(
j);
139 ++diff_base_functions;
143 for (; bb != nb_base_functions; bb++) {
144 ++diff_base_functions;
158 boost::ptr_vector<MethodForForceScaling> &methods_op,
int tag,
162 dAta(data),
commonData(common_data),
tAg(tag), jAcobian(jacobian),
170 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
176 if (row_type != MBVERTEX)
187 dot_W.resize(3,
false);
188 a_res.resize(3,
false);
189 g.resize(3, 3,
false);
190 G.resize(3, 3,
false);
191 h.resize(3, 3,
false);
192 H.resize(3, 3,
false);
193 invH.resize(3, 3,
false);
194 F.resize(3, 3,
false);
200 for (
int dd = 0;
dd < 3;
dd++) {
205 int nb_gauss_pts = row_data.
getN().size1();
210 const std::vector<VectorDouble> &dot_spacial_vel =
213 const std::vector<MatrixDouble> &spatial_positions_grad =
216 const std::vector<MatrixDouble> &spatial_velocities_grad =
219 const std::vector<VectorDouble> &meshpos_vel =
222 const std::vector<MatrixDouble> &mesh_positions_gradient =
225 int nb_active_vars = 0;
226 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
232 for (
int nn1 = 0; nn1 < 3; nn1++) {
234 a[nn1] <<= dot_spacial_vel[gg][nn1];
237 for (
int nn1 = 0; nn1 < 3; nn1++) {
238 for (
int nn2 = 0; nn2 < 3; nn2++) {
240 h(nn1, nn2) <<= spatial_positions_grad[gg](nn1, nn2);
251 for (
int nn1 = 0; nn1 < 3; nn1++) {
252 for (
int nn2 = 0; nn2 < 3; nn2++) {
254 g(nn1, nn2) <<= spatial_velocities_grad[gg](nn1, nn2);
258 for (
int nn1 = 0; nn1 < 3; nn1++) {
260 dot_W(nn1) <<= meshpos_vel[gg][nn1];
263 for (
int nn1 = 0; nn1 < 3; nn1++) {
264 for (
int nn2 = 0; nn2 < 3; nn2++) {
266 H(nn1, nn2) <<= mesh_positions_gradient[gg](nn1, nn2);
287 const double rho0 = dAta.rho0;
292 t_G(
i,
j) = t_g(
i,
k) * t_invH(
k,
j);
293 t_a_res(
i) = t_a(
i) - t_a0(
i) + t_G(
i,
j) * t_dotW(
j);
297 t_F(
i,
j) = t_h(
i,
k)*t_invH(
k,
j);
298 t_a_res(
i) *= rho0 * detH;
303 t_a_res(
i) *= rho0 * detH;
310 for (
int rr = 0; rr < 3; rr++) {
311 a_res[rr] >>= res[rr];
317 active.resize(nb_active_vars);
319 for (
int nn1 = 0; nn1 < 3; nn1++) {
320 active[aa++] = dot_spacial_vel[gg][nn1];
322 for (
int nn1 = 0; nn1 < 3; nn1++) {
323 for (
int nn2 = 0; nn2 < 3; nn2++) {
324 if (fieldDisp && nn1 == nn2) {
325 active[aa++] = spatial_positions_grad[gg](nn1, nn2) + 1;
327 active[aa++] = spatial_positions_grad[gg](nn1, nn2);
333 for (
int nn1 = 0; nn1 < 3; nn1++) {
334 for (
int nn2 = 0; nn2 < 3; nn2++) {
335 active[aa++] = spatial_velocities_grad[gg](nn1, nn2);
338 for (
int nn1 = 0; nn1 < 3; nn1++) {
339 active[aa++] = meshpos_vel[gg][nn1];
341 for (
int nn1 = 0; nn1 < 3; nn1++) {
342 for (
int nn2 = 0; nn2 < 3; nn2++) {
343 active[aa++] = mesh_positions_gradient[gg](nn1, nn2);
353 r = ::function(
tAg, 3, nb_active_vars, &active[0], &res[0]);
356 "ADOL-C function evaluation with error r = %d",
r);
359 double val = getVolume() * getGaussPts()(3, gg);
364 for (
int nn1 = 0; nn1 < 3; nn1++) {
369 r = jacobian(
tAg, 3, nb_active_vars, &active[0],
373 "ADOL-C function evaluation with error");
375 double val = getVolume() * getGaussPts()(3, gg);
396 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
405 int nb_base_functions = row_data.
getN().size2();
414 for (
unsigned int gg = 0; gg < row_data.
getN().size1(); gg++) {
420 for (;
dd < nb_dofs / 3;
dd++) {
421 t_nf(
i) += base * res(
i);
425 for (;
dd != nb_base_functions;
dd++) {
430 if ((
unsigned int)nb_dofs > 3 * row_data.
getN().size2()) {
431 SETERRQ(PETSC_COMM_SELF, 1,
"data inconsistency");
448 if (forcesonlyonentities_ptr != NULL) {
464 &jac(1, 0), &jac(1, 1), &jac(1, 2),
465 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
472 double *base_ptr =
const_cast<double *
>(&col_data.
getN(gg)[0]);
476 for (
int dd = 0;
dd < nb_col / 3;
dd++) {
477 t_jac(
i,
j) += t_mass1(
i,
j) * base * getFEMethod()->ts_a;
509 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
511 for (
int dd = 0;
dd < nb_col / 3;
dd++) {
512 t_jac(
i,
j) += t_mass1(
i,
j) * base * getFEMethod()->ts_a;
513 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
528 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
541 int nb_base_functions = row_data.
getN().size2();
545 k.resize(nb_row, nb_col);
547 jac.resize(3, nb_col);
549 for (
unsigned int gg = 0; gg < row_data.
getN().size1(); gg++) {
552 CHKERR getJac(col_data, gg);
553 }
catch (
const std::exception &ex) {
554 std::ostringstream ss;
555 ss <<
"throw in method: " << ex.what() << std::endl;
556 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
565 for (; dd1 < nb_row / 3; dd1++) {
567 &jac(0, 0), &jac(0, 1), &jac(0, 2), &jac(1, 0), &jac(1, 1),
568 &jac(1, 2), &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
569 for (
int dd2 = 0; dd2 < nb_col / 3; dd2++) {
571 &
k(3 * dd1 + 0, 3 * dd2 + 0), &
k(3 * dd1 + 0, 3 * dd2 + 1),
572 &
k(3 * dd1 + 0, 3 * dd2 + 2), &
k(3 * dd1 + 1, 3 * dd2 + 0),
573 &
k(3 * dd1 + 1, 3 * dd2 + 1), &
k(3 * dd1 + 1, 3 * dd2 + 2),
574 &
k(3 * dd1 + 2, 3 * dd2 + 0), &
k(3 * dd1 + 2, 3 * dd2 + 1),
575 &
k(3 * dd1 + 2, 3 * dd2 + 2));
576 t_k(
i,
j) += base * t_jac(
i,
j);
586 for (; dd1 != nb_base_functions; dd1++) {
592 if (!forcesOnlyOnEntities.empty()) {
595 VectorDofs::iterator dit = dofs.begin();
596 for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
597 if (forcesOnlyOnEntities.find((*dit)->getEnt()) ==
598 forcesOnlyOnEntities.end()) {
627 &jac(1, 0), &jac(1, 1), &jac(1, 2),
628 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
655 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
657 for (
int dd = 0;
dd < nb_col / 3;
dd++) {
658 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
675 double *base_ptr =
const_cast<double *
>(&col_data.
getN(gg)[0]);
678 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
681 &jac(1, 0), &jac(1, 1), &jac(1, 2),
682 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
683 const int u = 3 + 9 + 9;
690 const int s = 3 + 9 + 9 + 3;
718 for (
int dd = 0;
dd < nb_col / 3;
dd++) {
719 t_jac(
i,
j) += t_mass1(
i,
j) * base * getFEMethod()->ts_a;
720 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
742 if (row_type != MBVERTEX) {
745 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
752 for (
unsigned int gg = 0; gg < row_data.
getN().size1(); gg++) {
753 double val = getVolume() * getGaussPts()(3, gg);
754 double rho0 = dAta.rho0;
771 noalias(F) = prod(
h, invH);
781 energy += 0.5 * (
rho * val) * inner_prod(
v,
v);
783 CHKERR VecSetValue(V, 0, energy, ADD_VALUES);
791 int tag,
bool jacobian)
794 dAta(data),
commonData(common_data),
tAg(tag), jAcobian(jacobian),
802 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
808 if (row_type != MBVERTEX)
829 for (
int dd = 0;
dd < 3;
dd++) {
835 int nb_gauss_pts = row_data.
getN().size1();
840 int nb_active_vars = 0;
841 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
847 for (
int nn1 = 0; nn1 < 3; nn1++) {
852 for (
int nn1 = 0; nn1 < 3; nn1++) {
860 for (
int nn1 = 0; nn1 < 3; nn1++) {
861 for (
int nn2 = 0; nn2 < 3; nn2++) {
873 for (
int nn1 = 0; nn1 < 3; nn1++) {
881 for (
int nn1 = 0; nn1 < 3; nn1++) {
882 for (
int nn2 = 0; nn2 < 3; nn2++) {
913 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
915 t_F(
i,
j) = t_h(
i,
j);
918 t_dot_u(
i) = t_dot_w(
i) + t_F(
i,
j) * t_dot_W(
j);
919 t_a_res(
i) = t_v(
i) - t_dot_u(
i);
925 for (
int rr = 0; rr < 3; rr++) {
926 a_res[rr] >>= res[rr];
931 active.resize(nb_active_vars);
933 for (
int nn1 = 0; nn1 < 3; nn1++) {
937 for (
int nn1 = 0; nn1 < 3; nn1++) {
944 for (
int nn1 = 0; nn1 < 3; nn1++) {
945 for (
int nn2 = 0; nn2 < 3; nn2++) {
946 if (fieldDisp && nn1 == nn2) {
958 for (
int nn1 = 0; nn1 < 3; nn1++) {
965 for (
int nn1 = 0; nn1 < 3; nn1++) {
966 for (
int nn2 = 0; nn2 < 3; nn2++) {
979 r = ::function(
tAg, 3, nb_active_vars, &active[0], &res[0]);
982 "ADOL-C function evaluation with error");
985 double val = getVolume() * getGaussPts()(3, gg);
990 for (
int nn1 = 0; nn1 < 3; nn1++) {
994 r = jacobian(
tAg, 3, nb_active_vars, &active[0],
998 "ADOL-C function evaluation with error");
1000 double val = getVolume() * getGaussPts()(3, gg);
1020 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1029 int nb_base_functions = row_data.
getN().size2();
1037 for (
unsigned int gg = 0; gg < row_data.
getN().size1(); gg++) {
1043 for (;
dd < nb_dofs / 3;
dd++) {
1044 t_nf(
i) += base * res(
i);
1048 for (;
dd != nb_base_functions;
dd++) {
1053 if (row_data.
getIndices().size() > 3 * row_data.
getN().size2()) {
1054 SETERRQ(PETSC_COMM_SELF, 1,
"data inconsistency");
1057 &row_data.
getIndices()[0], &nf[0], ADD_VALUES);
1074 double *base_ptr =
const_cast<double *
>(&col_data.
getN(gg)[0]);
1077 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1078 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1087 for (
int dd = 0;
dd < nb_col / 3;
dd++) {
1088 t_jac(
i,
j) += t_mass1(
i,
j) * base;
1108 double *base_ptr =
const_cast<double *
>(&col_data.
getN(gg)[0]);
1111 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1112 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1126 for (
int dd = 0;
dd < nb_col / 3;
dd++) {
1127 t_jac(
i,
j) += t_mass1(
i,
j) * base * getFEMethod()->ts_a;
1133 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
1135 const int s = 3 + 3;
1160 for (
int dd = 0;
dd < nb_col / 3;
dd++) {
1161 t_jac(
i,
j) += t_mass1(
i,
j) * base * getFEMethod()->ts_a;
1162 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
1183 double *base_ptr =
const_cast<double *
>(&col_data.
getN(gg)[0]);
1186 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
1189 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1190 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1191 const int u = 3 + 3 + 9;
1198 const int s = 3 + 3 + 9 + 3;
1226 for (
int dd = 0;
dd < nb_col / 3;
dd++) {
1227 t_jac(
i,
j) += t_mass1(
i,
j) * base * getFEMethod()->ts_a;
1228 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
1243 dAta(data),
commonData(common_data),
tAg(tag), jAcobian(jacobian),
1252 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1258 if (row_type != MBVERTEX)
1277 for (
int dd = 0;
dd < 3;
dd++) {
1282 int nb_gauss_pts = row_data.
getN().size1();
1287 int nb_active_vars = 0;
1288 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1294 for (
int nn1 = 0; nn1 < 3; nn1++) {
1301 for (
int nn1 = 0; nn1 < 3; nn1++) {
1306 for (
int nn1 = 0; nn1 < 3; nn1++) {
1307 for (
int nn2 = 0; nn2 < 3; nn2++) {
1314 for (
int nn1 = 0; nn1 < 3; nn1++) {
1315 for (
int nn2 = 0; nn2 < 3; nn2++) {
1328 for (
int nn1 = 0; nn1 < 3; nn1++) {
1329 for (
int nn2 = 0; nn2 < 3; nn2++) {
1360 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
1361 t_G(
i,
j) = t_g(
i,
k) * t_invH(
k,
j);
1362 t_a_T(
i) = t_F(
k,
i) * t_a(
k) + t_G(
k,
i) * t_v(
k);
1363 const auto rho0 = dAta.rho0;
1364 t_a_T(
i) = -rho0 * detH;
1367 for (
int nn = 0; nn < 3; nn++) {
1373 active.resize(nb_active_vars);
1375 for (
int nn1 = 0; nn1 < 3; nn1++) {
1381 for (
int nn1 = 0; nn1 < 3; nn1++) {
1385 for (
int nn1 = 0; nn1 < 3; nn1++) {
1386 for (
int nn2 = 0; nn2 < 3; nn2++) {
1392 for (
int nn1 = 0; nn1 < 3; nn1++) {
1393 for (
int nn2 = 0; nn2 < 3; nn2++) {
1394 if (fieldDisp && nn1 == nn2) {
1407 for (
int nn1 = 0; nn1 < 3; nn1++) {
1408 for (
int nn2 = 0; nn2 < 3; nn2++) {
1421 r = ::function(
tAg, 3, nb_active_vars, &active[0], &res[0]);
1424 "ADOL-C function evaluation with error r = %d",
r);
1427 double val = getVolume() * getGaussPts()(3, gg);
1432 for (
int nn1 = 0; nn1 < 3; nn1++) {
1436 r = jacobian(
tAg, 3, nb_active_vars, &active[0],
1440 "ADOL-C function evaluation with error");
1442 double val = getVolume() * getGaussPts()(3, gg);
1447 }
catch (
const std::exception &ex) {
1448 std::ostringstream ss;
1449 ss <<
"throw in method: " << ex.what() << std::endl;
1450 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
1460 Range *forcesonlyonentities_ptr)
1464 if (forcesonlyonentities_ptr != NULL) {
1475 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1489 int nb_base_functions = row_data.
getN().size2();
1492 for (
unsigned int gg = 0; gg < row_data.
getN().size1(); gg++) {
1498 for (;
dd < nb_dofs / 3;
dd++) {
1499 t_nf(
i) += base * res(
i);
1503 for (;
dd != nb_base_functions;
dd++) {
1508 if (row_data.
getIndices().size() > 3 * row_data.
getN().size2()) {
1509 SETERRQ(PETSC_COMM_SELF, 1,
"data inconsistency");
1511 if (!forcesOnlyOnEntities.empty()) {
1514 VectorDofs::iterator dit = dofs.begin();
1515 for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
1516 if (forcesOnlyOnEntities.find((*dit)->getEnt()) ==
1517 forcesOnlyOnEntities.end()) {
1524 &nf[0], ADD_VALUES);
1527 &row_data.
getIndices()[0], &nf[0], ADD_VALUES);
1530 }
catch (
const std::exception &ex) {
1531 std::ostringstream ss;
1532 ss <<
"throw in method: " << ex.what() << std::endl;
1533 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
1544 Range *forcesonlyonentities_ptr)
1546 vel_field,
field_name, data, common_data, forcesonlyonentities_ptr) {}
1554 double *base_ptr =
const_cast<double *
>(&col_data.
getN(gg)[0]);
1557 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
1560 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1561 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1569 const int s = 3 + 3;
1597 for (
int dd = 0;
dd < nb_col / 3;
dd++) {
1598 t_jac(
i,
j) += t_mass1(
i,
j) * base * getFEMethod()->ts_a;
1599 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
1612 Range *forcesonlyonentities_ptr)
1614 vel_field,
field_name, data, common_data, forcesonlyonentities_ptr) {}
1623 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
1626 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1627 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1628 const int s = 3 + 3 + 9;
1656 for (
int dd = 0;
dd < nb_col / 3;
dd++) {
1657 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
1669 Range *forcesonlyonentities_ptr)
1671 vel_field,
field_name, data, common_data, forcesonlyonentities_ptr) {}
1680 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
1683 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1684 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1685 const int s = 3 + 3 + 9 + 9;
1713 for (
int dd = 0;
dd < nb_col / 3;
dd++) {
1714 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
1723 const std::string spatial_position_field)
1724 :
mField(m_field), tS(_ts), velocityField(velocity_field),
1725 spatialPositionField(spatial_position_field), jacobianLag(-1) {}
1731 case CTX_TSSETIFUNCTION: {
1737 problemPtr,
COL, ts_u, INSERT_VALUES, SCATTER_REVERSE);
1739 problemPtr, velocityField,
"DOT_" + velocityField,
COL, ts_u_t,
1740 INSERT_VALUES, SCATTER_REVERSE);
1742 problemPtr, spatialPositionField,
"DOT_" + spatialPositionField,
COL,
1743 ts_u_t, INSERT_VALUES, SCATTER_REVERSE);
1746 case CTX_TSSETIJACOBIAN: {
1772 int id = it->getMeshsetId();
1773 EntityHandle meshset = it->getMeshset();
1778 CHKERR it->getAttributeDataStructure(mydata);
1790 CHKERR it->getAttributeDataStructure(mydata);
1791 if (mydata.
data.User1 == 0)
1794 EntityHandle meshset = it->getMeshset();
1796 tets = subtract(tets, added_tets);
1799 int id = it->getMeshsetId();
1814 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr) {
1817 if (!block_sets_ptr)
1819 "Pointer to block of sets is null");
1824 CHKERR it->getAttributeDataStructure(mydata);
1825 int id = it->getMeshsetId();
1826 auto &block_data = (*block_sets_ptr)[id];
1827 EntityHandle meshset = it->getMeshset();
1829 block_data.tEts,
true);
1830 block_data.rho0 = mydata.
data.density;
1831 block_data.a0.resize(3);
1832 block_data.a0[0] = mydata.
data.acceleration_x;
1833 block_data.a0[1] = mydata.
data.acceleration_y;
1834 block_data.a0[2] = mydata.
data.acceleration_z;
1841 string element_name,
string velocity_field_name,
1842 string spatial_position_field_name,
string material_position_field_name,
1850 velocity_field_name);
1852 velocity_field_name);
1854 velocity_field_name);
1856 element_name, spatial_position_field_name);
1858 element_name, spatial_position_field_name);
1860 element_name, spatial_position_field_name);
1864 element_name, material_position_field_name);
1866 element_name, material_position_field_name);
1868 element_name,
"DOT_" + material_position_field_name);
1871 element_name, material_position_field_name);
1874 element_name,
"DOT_" + velocity_field_name);
1876 element_name,
"DOT_" + spatial_position_field_name);
1884 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
1886 Range add_tets = sit->second.tEts;
1887 if (!tets.empty()) {
1888 add_tets = intersect(add_tets, tets);
1898 string element_name,
string velocity_field_name,
1899 string spatial_position_field_name,
string material_position_field_name,
1907 velocity_field_name);
1909 velocity_field_name);
1911 velocity_field_name);
1913 element_name, spatial_position_field_name);
1915 element_name, spatial_position_field_name);
1919 element_name, material_position_field_name);
1921 element_name,
"DOT_" + material_position_field_name);
1924 element_name, material_position_field_name);
1927 element_name,
"DOT_" + velocity_field_name);
1929 element_name,
"DOT_" + spatial_position_field_name);
1937 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
1939 Range add_tets = sit->second.tEts;
1940 if (!tets.empty()) {
1941 add_tets = intersect(add_tets, tets);
1951 string element_name,
string velocity_field_name,
1952 string spatial_position_field_name,
string material_position_field_name,
1960 velocity_field_name);
1962 velocity_field_name);
1964 element_name, spatial_position_field_name);
1966 element_name, spatial_position_field_name);
1970 element_name, material_position_field_name);
1972 element_name, material_position_field_name);
1974 element_name,
"DOT_" + material_position_field_name);
1977 element_name, material_position_field_name);
1980 element_name,
"DOT_" + velocity_field_name);
1982 element_name,
"DOT_" + spatial_position_field_name);
1989 if (intersected != NULL) {
1991 tets = *intersected;
1993 tets = intersect(*intersected, tets);
1997 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
1999 Range add_tets = sit->second.tEts;
2000 if (!tets.empty()) {
2001 add_tets = intersect(add_tets, tets);
2011 string velocity_field_name,
string spatial_position_field_name,
2012 string material_position_field_name,
bool ale,
bool linear) {
2028 "DOT_" + spatial_position_field_name,
commonData));
2034 "DOT_" + material_position_field_name,
commonData));
2036 feMassRhs.meshPositionsFieldName = material_position_field_name;
2039 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
2056 "DOT_" + spatial_position_field_name,
commonData));
2062 "DOT_" + material_position_field_name,
commonData));
2064 feMassLhs.meshPositionsFieldName = material_position_field_name;
2076 spatial_position_field_name, spatial_position_field_name, sit->second,
2081 spatial_position_field_name, material_position_field_name,
2084 feMassLhs.meshPositionsFieldName = material_position_field_name;
2090 feEnergy.getOpPtrVector().push_back(
2092 feEnergy.getOpPtrVector().push_back(
2097 feEnergy.meshPositionsFieldName = material_position_field_name;
2109 string velocity_field_name,
string spatial_position_field_name,
2110 string material_position_field_name,
bool ale) {
2118 feVelRhs.getOpPtrVector().push_back(
2120 feVelRhs.getOpPtrVector().push_back(
2122 feVelRhs.getOpPtrVector().push_back(
2126 "DOT_" + spatial_position_field_name,
commonData));
2131 "DOT_" + material_position_field_name,
commonData));
2133 feVelRhs.meshPositionsFieldName = material_position_field_name;
2136 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
2140 feVelRhs.getOpPtrVector().push_back(
2145 feVelLhs.getOpPtrVector().push_back(
2147 feVelLhs.getOpPtrVector().push_back(
2149 feVelLhs.getOpPtrVector().push_back(
2153 "DOT_" + spatial_position_field_name,
commonData));
2158 "DOT_" + material_position_field_name,
commonData));
2160 feVelLhs.meshPositionsFieldName = material_position_field_name;
2168 velocity_field_name, velocity_field_name, sit->second,
commonData));
2170 velocity_field_name, spatial_position_field_name, sit->second,
2175 velocity_field_name, material_position_field_name, sit->second,
2178 feVelLhs.meshPositionsFieldName = material_position_field_name;
2187 string velocity_field_name,
string spatial_position_field_name,
2188 string material_position_field_name,
Range *forces_on_entities_ptr) {
2196 feTRhs.getOpPtrVector().push_back(
2198 feTRhs.getOpPtrVector().push_back(
2200 feTRhs.getOpPtrVector().push_back(
2202 feTRhs.getOpPtrVector().push_back(
2205 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
2207 feTRhs.getOpPtrVector().push_back(
2209 material_position_field_name, sit->second,
commonData,
tAg,
false));
2211 material_position_field_name, sit->second,
commonData,
2212 forces_on_entities_ptr));
2216 feTLhs.getOpPtrVector().push_back(
2218 feTLhs.getOpPtrVector().push_back(
2220 feTLhs.getOpPtrVector().push_back(
2228 feTLhs.getOpPtrVector().push_back(
2230 material_position_field_name, sit->second,
commonData,
tAg));
2231 feTLhs.getOpPtrVector().push_back(
2233 material_position_field_name, velocity_field_name, sit->second,
2235 feTLhs.getOpPtrVector().push_back(
2237 material_position_field_name, spatial_position_field_name,
2238 sit->second,
commonData, forces_on_entities_ptr));
2239 feTLhs.getOpPtrVector().push_back(
2241 material_position_field_name, material_position_field_name,
2242 sit->second,
commonData, forces_on_entities_ptr));
2249 string velocity_field_name,
string spatial_position_field_name,
2250 string material_position_field_name,
bool linear) {
2268 feMassRhs.meshPositionsFieldName = material_position_field_name;
2270 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
2289 feMassLhs.meshPositionsFieldName = material_position_field_name;
2297 spatial_position_field_name, spatial_position_field_name, sit->second,
2300 feMassLhs.meshPositionsFieldName = material_position_field_name;
2314 feMassAuxLhs.meshPositionsFieldName = material_position_field_name;
2322 spatial_position_field_name, spatial_position_field_name, sit->second,
2325 feMassAuxLhs.meshPositionsFieldName = material_position_field_name;
2330 feEnergy.getOpPtrVector().push_back(
2332 feEnergy.getOpPtrVector().push_back(
2337 feEnergy.meshPositionsFieldName = material_position_field_name;
2353 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
2361#if PETSC_VERSION_GE(3, 5, 3)
2362 CHKERR MatCreateVecs(
K, &u, &Ku);
2365 CHKERR MatGetVecs(
K, &u, &Ku);
2368 CHKERR MatDuplicate(
K, MAT_SHARE_NONZERO_PATTERN, &barK);
2382 CHKERR MatDestroy(&barK);
2383 iNitialized =
false;
2393 CHKERR PetscObjectGetComm((PetscObject)shellMat, &comm);
2394 CHKERR PCCreate(comm, &pC);
2417 if (ts_ctx != CTX_TSSETIFUNCTION) {
2419 "It is used to residual of velocities");
2421 if (!shellMatCtx->iNitialized) {
2422 CHKERR shellMatCtx->iNit();
2425 CHKERR VecScatterBegin(shellMatCtx->scatterU, ts_u_t, shellMatCtx->u,
2426 INSERT_VALUES, SCATTER_FORWARD);
2427 CHKERR VecScatterEnd(shellMatCtx->scatterU, ts_u_t, shellMatCtx->u,
2428 INSERT_VALUES, SCATTER_FORWARD);
2429 CHKERR VecScatterBegin(shellMatCtx->scatterV, ts_u, shellMatCtx->v,
2430 INSERT_VALUES, SCATTER_FORWARD);
2431 CHKERR VecScatterEnd(shellMatCtx->scatterV, ts_u, shellMatCtx->v,
2432 INSERT_VALUES, SCATTER_FORWARD);
2433 CHKERR VecAXPY(shellMatCtx->v, -1, shellMatCtx->u);
2434 CHKERR VecScatterBegin(shellMatCtx->scatterV, shellMatCtx->v, ts_F,
2435 ADD_VALUES, SCATTER_REVERSE);
2436 CHKERR VecScatterEnd(shellMatCtx->scatterV, shellMatCtx->v, ts_F, ADD_VALUES,
2444#ifdef __DIRICHLET_HPP__
2446ConvectiveMassElement::ShellMatrixElement::ShellMatrixElement(
2450MoFEMErrorCode ConvectiveMassElement::ShellMatrixElement::preProcess() {
2453 if (ts_ctx != CTX_TSSETIJACOBIAN) {
2455 "It is used to calculate shell matrix only");
2458 shellMatCtx->ts_a = ts_a;
2459 DirichletBcPtr->copyTs(*((
TSMethod *)
this));
2461 DirichletBcPtr->dIag = 1;
2462 DirichletBcPtr->ts_B = shellMatCtx->K;
2463 CHKERR MatZeroEntries(shellMatCtx->K);
2465 LoopsToDoType::iterator itk = loopK.begin();
2466 for (; itk != loopK.end(); itk++) {
2467 itk->second->copyTs(*((
TSMethod *)
this));
2468 itk->second->ts_B = shellMatCtx->K;
2471 LoopsToDoType::iterator itam = loopAuxM.begin();
2472 for (; itam != loopAuxM.end(); itam++) {
2473 itam->second->copyTs(*((
TSMethod *)
this));
2474 itam->second->ts_B = shellMatCtx->K;
2478 CHKERR MatAssemblyBegin(shellMatCtx->K, MAT_FINAL_ASSEMBLY);
2479 CHKERR MatAssemblyEnd(shellMatCtx->K, MAT_FINAL_ASSEMBLY);
2481 DirichletBcPtr->dIag = 0;
2482 DirichletBcPtr->ts_B = shellMatCtx->M;
2483 CHKERR MatZeroEntries(shellMatCtx->M);
2485 LoopsToDoType::iterator itm = loopM.begin();
2486 for (; itm != loopM.end(); itm++) {
2487 itm->second->copyTs(*((
TSMethod *)
this));
2488 itm->second->ts_B = shellMatCtx->M;
2492 CHKERR MatAssemblyBegin(shellMatCtx->M, MAT_FINAL_ASSEMBLY);
2493 CHKERR MatAssemblyEnd(shellMatCtx->M, MAT_FINAL_ASSEMBLY);
2496 CHKERR MatZeroEntries(shellMatCtx->barK);
2497 CHKERR MatCopy(shellMatCtx->K, shellMatCtx->barK, SAME_NONZERO_PATTERN);
2498 CHKERR MatAXPY(shellMatCtx->barK, ts_a, shellMatCtx->M, SAME_NONZERO_PATTERN);
2499 CHKERR MatAssemblyBegin(shellMatCtx->barK, MAT_FINAL_ASSEMBLY);
2500 CHKERR MatAssemblyEnd(shellMatCtx->barK, MAT_FINAL_ASSEMBLY);
Operators and data structures for mass and convective mass element.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ BODYFORCESSET
block name is "BODY_FORCES"
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
@ MOFEM_OPERATION_UNSUCCESSFUL
@ 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 ...
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual bool check_field(const std::string &name) const =0
check if field is in database
virtual MoFEMErrorCode problem_basic_method_postProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
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)
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
UBlasVector< int > VectorInt
implementation of Data Operators for Forces and Sources
auto createSmartVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr)
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 determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
constexpr auto VecSetValues
constexpr auto MatSetValues
const double r
rate factor
constexpr auto field_name
data for calculation inertia forces
common data used by volume elements
std::vector< std::vector< double * > > jacTRowPtr
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
std::vector< std::vector< double * > > jacVelRowPtr
std::vector< VectorDouble > valVel
std::vector< MatrixDouble > jacMass
std::vector< VectorDouble > valT
std::vector< MatrixDouble > jacVel
std::vector< VectorDouble > valMass
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
std::vector< MatrixDouble > jacT
std::vector< std::vector< double * > > jacMassRowPtr
MoFEMErrorCode postProcess()
int getRule(int order)
it is used to calculate nb. of Gauss integration points
MoFEMErrorCode preProcess()
MyVolumeFE(MoFEM::Interface &m_field)
OpEnergy(const std::string field_name, BlockData &data, CommonData &common_data, SmartPetscObj< Vec > v)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpEshelbyDynamicMaterialMomentumJacobian(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian=true)
OpEshelbyDynamicMaterialMomentumLhs_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpEshelbyDynamicMaterialMomentumLhs_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpEshelbyDynamicMaterialMomentumLhs_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
Range forcesOnlyOnEntities
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpEshelbyDynamicMaterialMomentumRhs(const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
OpGetCommonDataAtGaussPts(const std::string field_name, CommonData &common_data)
OpGetDataAtGaussPts(const std::string field_name, std::vector< VectorDouble > &values_at_gauss_pts, std::vector< MatrixDouble > &gardient_at_gauss_pts)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
operator calculating deformation gradient
OpMassJacobian(const std::string field_name, BlockData &data, CommonData &common_data, boost::ptr_vector< MethodForForceScaling > &methods_op, int tag, bool linear=false)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpMassLhs_dM_dX(const std::string field_name, const std::string col_field, BlockData &data, CommonData &common_data)
MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpMassLhs_dM_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr=NULL)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
Range forcesOnlyOnEntities
MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpMassLhs_dM_dx(const std::string field_name, const std::string col_field, BlockData &data, CommonData &common_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpMassRhs(const std::string field_name, BlockData &data, CommonData &common_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpVelocityJacobian(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian=true)
OpVelocityLhs_dV_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpVelocityLhs_dV_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpVelocityLhs_dV_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpVelocityRhs(const std::string field_name, BlockData &data, CommonData &common_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
MoFEMErrorCode preProcess()
Calculate inconsistency between approximation of velocities and velocities calculated from displaceme...
ShellResidualElement(MoFEM::Interface &m_field)
UpdateAndControl(MoFEM::Interface &m_field, TS _ts, const std::string velocity_field, const std::string spatial_position_field)
MoFEMErrorCode postProcess()
MoFEMErrorCode preProcess()
Scatter values from t_u_dt on the fields.
structure grouping operators and data used for calculation of mass (convective) element \ nonlinear_e...
ConvectiveMassElement(MoFEM::Interface &m_field, short int tag)
MoFEMErrorCode setVelocityOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false)
MyVolumeFE feEnergy
calculate kinetic energy
MoFEMErrorCode setShellMatrixMassOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool linear=false)
MyVolumeFE feVelRhs
calculate right hand side for tetrahedral elements
MoFEMErrorCode addVelocityElement(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel())
MoFEMErrorCode addConvectiveMassElement(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel())
MoFEMErrorCode setKinematicEshelbyOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", Range *forces_on_entities_ptr=NULL)
MoFEMErrorCode addEshelbyDynamicMaterialMomentum(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel(), Range *intersected=NULL)
MyVolumeFE feTRhs
calculate right hand side for tetrahedral elements
MyVolumeFE feMassRhs
calculate right hand side for tetrahedral elements
MyVolumeFE feTLhs
calculate left hand side for tetrahedral elements
MoFEMErrorCode setConvectiveMassOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, bool linear=false)
boost::ptr_vector< MethodForForceScaling > methodsOp
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
MoFEM::Interface & mField
MoFEMErrorCode setBlocks()
MyVolumeFE feVelLhs
calculate left hand side for tetrahedral elements
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
Body force data structure.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
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::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
structure to get information form mofem into EntitiesFieldData
Elastic material data structure.
data structure for TS (time stepping) context
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Volume finite element base.