12#include <adolc/adolc.h>
19#error "MoFEM need to be compiled with ADOL-C"
24 meshPositionsFieldName =
"NoNE";
26 auto create_vec = [&]() {
27 constexpr int ghosts[] = {0};
43 CHKERR VolumeElementForcesAndSourcesCore::preProcess();
59 CHKERR VolumeElementForcesAndSourcesCore::postProcess();
64 CHKERR VecAssemblyBegin(V);
83 std::vector<VectorDouble> &values_at_gauss_pts,
84 std::vector<MatrixDouble> &gardient_at_gauss_pts)
87 valuesAtGaussPts(values_at_gauss_pts),
88 gradientAtGaussPts(gardient_at_gauss_pts), zeroAtType(MBVERTEX) {}
98 int nb_gauss_pts = data.
getN().size1();
99 int nb_base_functions = data.
getN().size2();
103 valuesAtGaussPts.resize(nb_gauss_pts);
104 gradientAtGaussPts.resize(nb_gauss_pts);
105 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
106 valuesAtGaussPts[gg].resize(3);
107 gradientAtGaussPts[gg].resize(3, 3);
110 if (type == zeroAtType) {
111 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
112 valuesAtGaussPts[gg].clear();
113 gradientAtGaussPts[gg].clear();
122 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
125 &valuesAtGaussPts[gg][1],
126 &valuesAtGaussPts[gg][2]);
128 &gradientAtGaussPts[gg](0, 0), &gradientAtGaussPts[gg](0, 1),
129 &gradientAtGaussPts[gg](0, 2), &gradientAtGaussPts[gg](1, 0),
130 &gradientAtGaussPts[gg](1, 1), &gradientAtGaussPts[gg](1, 2),
131 &gradientAtGaussPts[gg](2, 0), &gradientAtGaussPts[gg](2, 1),
132 &gradientAtGaussPts[gg](2, 2));
134 for (; bb != nb_dofs / 3; bb++) {
135 values(
i) += base_function * field_data(
i);
136 gradient(
i,
j) += field_data(
i) * diff_base_functions(
j);
137 ++diff_base_functions;
141 for (; bb != nb_base_functions; bb++) {
142 ++diff_base_functions;
156 boost::ptr_vector<MethodForForceScaling> &methods_op,
int tag,
160 dAta(data),
commonData(common_data),
tAg(tag), jAcobian(jacobian),
167 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
173 if (row_type != MBVERTEX)
184 dot_W.resize(3,
false);
185 a_res.resize(3,
false);
186 g.resize(3, 3,
false);
187 G.resize(3, 3,
false);
188 h.resize(3, 3,
false);
189 H.resize(3, 3,
false);
190 invH.resize(3, 3,
false);
191 F.resize(3, 3,
false);
194 std::fill(dot_W.begin(), dot_W.end(), 0);
195 std::fill(
H.data().begin(),
H.data().end(), 0);
196 std::fill(invH.data().begin(), invH.data().end(), 0);
197 for (
int ii = 0; ii != 3; ii++) {
202 int nb_gauss_pts = row_data.
getN().size1();
207 const std::vector<VectorDouble> &dot_spacial_vel =
210 const std::vector<MatrixDouble> &spatial_positions_grad =
213 const std::vector<MatrixDouble> &spatial_velocities_grad =
216 const std::vector<VectorDouble> &meshpos_vel =
219 const std::vector<MatrixDouble> &mesh_positions_gradient =
222 int nb_active_vars = 0;
223 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
229 for (
int nn1 = 0; nn1 < 3; nn1++) {
231 a[nn1] <<= dot_spacial_vel[gg][nn1];
234 for (
int nn1 = 0; nn1 < 3; nn1++) {
235 for (
int nn2 = 0; nn2 < 3; nn2++) {
237 h(nn1, nn2) <<= spatial_positions_grad[gg](nn1, nn2);
248 for (
int nn1 = 0; nn1 < 3; nn1++) {
249 for (
int nn2 = 0; nn2 < 3; nn2++) {
251 g(nn1, nn2) <<= spatial_velocities_grad[gg](nn1, nn2);
255 for (
int nn1 = 0; nn1 < 3; nn1++) {
257 dot_W(nn1) <<= meshpos_vel[gg][nn1];
260 for (
int nn1 = 0; nn1 < 3; nn1++) {
261 for (
int nn2 = 0; nn2 < 3; nn2++) {
263 H(nn1, nn2) <<= mesh_positions_gradient[gg](nn1, nn2);
284 const double rho0 = dAta.rho0;
289 t_G(
i,
j) = t_g(
i,
k) * t_invH(
k,
j);
290 t_a_res(
i) = t_a(
i) - t_a0(
i) + t_G(
i,
j) * t_dotW(
j);
296 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
297 t_a_res(
i) *= rho0 * detH;
302 t_a_res(
i) *= rho0 * detH;
308 for (
int rr = 0; rr < 3; rr++) {
309 a_res[rr] >>= res[rr];
315 active.resize(nb_active_vars);
317 for (
int nn1 = 0; nn1 < 3; nn1++) {
318 active[aa++] = dot_spacial_vel[gg][nn1];
320 for (
int nn1 = 0; nn1 < 3; nn1++) {
321 for (
int nn2 = 0; nn2 < 3; nn2++) {
322 if (fieldDisp && nn1 == nn2) {
323 active[aa++] = spatial_positions_grad[gg](nn1, nn2) + 1;
325 active[aa++] = spatial_positions_grad[gg](nn1, nn2);
331 for (
int nn1 = 0; nn1 < 3; nn1++) {
332 for (
int nn2 = 0; nn2 < 3; nn2++) {
333 active[aa++] = spatial_velocities_grad[gg](nn1, nn2);
336 for (
int nn1 = 0; nn1 < 3; nn1++) {
337 active[aa++] = meshpos_vel[gg][nn1];
339 for (
int nn1 = 0; nn1 < 3; nn1++) {
340 for (
int nn2 = 0; nn2 < 3; nn2++) {
341 active[aa++] = mesh_positions_gradient[gg](nn1, nn2);
351 r = ::function(
tAg, 3, nb_active_vars, &active[0], &res[0]);
354 "ADOL-C function evaluation with error r = %d", r);
357 double val = getVolume() * getGaussPts()(3, gg);
363 for (
int nn1 = 0; nn1 < 3; nn1++) {
368 r = jacobian(
tAg, 3, nb_active_vars, &active[0],
372 "ADOL-C function evaluation with error");
374 double val = getVolume() * getGaussPts()(3, gg);
395 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
404 int nb_base_functions = row_data.
getN().size2();
413 for (
unsigned int gg = 0; gg < row_data.
getN().size1(); gg++) {
419 for (; dd < nb_dofs / 3; dd++) {
420 t_nf(
i) += base * res(
i);
424 for (; dd != nb_base_functions; dd++) {
429 if ((
unsigned int)nb_dofs > 3 * row_data.
getN().size2()) {
430 SETERRQ(PETSC_COMM_SELF, 1,
"data inconsistency");
447 if (forcesonlyonentities_ptr != NULL) {
463 &jac(1, 0), &jac(1, 1), &jac(1, 2),
464 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
471 double *base_ptr =
const_cast<double *
>(&col_data.
getN(gg)[0]);
475 for (
int dd = 0; dd < nb_col / 3; dd++) {
476 t_jac(
i,
j) += t_mass1(
i,
j) * base * getFEMethod()->ts_a;
508 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
510 for (
int dd = 0; dd < nb_col / 3; dd++) {
511 t_jac(
i,
j) += t_mass1(
i,
j) * base * getFEMethod()->ts_a;
512 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
527 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
540 int nb_base_functions = row_data.
getN().size2();
544 k.resize(nb_row, nb_col);
546 jac.resize(3, nb_col);
548 for (
unsigned int gg = 0; gg < row_data.
getN().size1(); gg++) {
551 CHKERR getJac(col_data, gg);
552 }
catch (
const std::exception &ex) {
553 std::ostringstream ss;
554 ss <<
"throw in method: " << ex.what() << std::endl;
555 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
564 for (; dd1 < nb_row / 3; dd1++) {
566 &jac(0, 0), &jac(0, 1), &jac(0, 2), &jac(1, 0), &jac(1, 1),
567 &jac(1, 2), &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
568 for (
int dd2 = 0; dd2 < nb_col / 3; dd2++) {
570 &
k(3 * dd1 + 0, 3 * dd2 + 0), &
k(3 * dd1 + 0, 3 * dd2 + 1),
571 &
k(3 * dd1 + 0, 3 * dd2 + 2), &
k(3 * dd1 + 1, 3 * dd2 + 0),
572 &
k(3 * dd1 + 1, 3 * dd2 + 1), &
k(3 * dd1 + 1, 3 * dd2 + 2),
573 &
k(3 * dd1 + 2, 3 * dd2 + 0), &
k(3 * dd1 + 2, 3 * dd2 + 1),
574 &
k(3 * dd1 + 2, 3 * dd2 + 2));
575 t_k(
i,
j) += base * t_jac(
i,
j);
585 for (; dd1 != nb_base_functions; dd1++) {
591 if (!forcesOnlyOnEntities.empty()) {
594 VectorDofs::iterator dit = dofs.begin();
595 for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
596 if (forcesOnlyOnEntities.find((*dit)->getEnt()) ==
597 forcesOnlyOnEntities.end()) {
626 &jac(1, 0), &jac(1, 1), &jac(1, 2),
627 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
654 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
656 for (
int dd = 0; dd < nb_col / 3; dd++) {
657 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
674 double *base_ptr =
const_cast<double *
>(&col_data.
getN(gg)[0]);
677 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
680 &jac(1, 0), &jac(1, 1), &jac(1, 2),
681 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
682 const int u = 3 + 9 + 9;
689 const int s = 3 + 9 + 9 + 3;
717 for (
int dd = 0; dd < nb_col / 3; dd++) {
718 t_jac(
i,
j) += t_mass1(
i,
j) * base * getFEMethod()->ts_a;
719 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
741 if (row_type != MBVERTEX) {
744 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
751 for (
unsigned int gg = 0; gg < row_data.
getN().size1(); gg++) {
752 double val = getVolume() * getGaussPts()(3, gg);
753 double rho0 = dAta.rho0;
770 noalias(
F) = prod(
h, invH);
780 energy += 0.5 * (
rho * val) * inner_prod(
v,
v);
782 CHKERR VecSetValue(V, 0, energy, ADD_VALUES);
790 int tag,
bool jacobian)
793 dAta(data),
commonData(common_data),
tAg(tag), jAcobian(jacobian),
800 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
806 if (row_type != MBVERTEX)
827 for (
int dd = 0; dd < 3; dd++) {
833 int nb_gauss_pts = row_data.
getN().size1();
838 int nb_active_vars = 0;
839 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
845 for (
int nn1 = 0; nn1 < 3; nn1++) {
850 for (
int nn1 = 0; nn1 < 3; nn1++) {
858 for (
int nn1 = 0; nn1 < 3; nn1++) {
859 for (
int nn2 = 0; nn2 < 3; nn2++) {
871 for (
int nn1 = 0; nn1 < 3; nn1++) {
879 for (
int nn1 = 0; nn1 < 3; nn1++) {
880 for (
int nn2 = 0; nn2 < 3; nn2++) {
911 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
913 t_F(
i,
j) = t_h(
i,
j);
916 t_dot_u(
i) = t_dot_w(
i) + t_F(
i,
j) * t_dot_W(
j);
917 t_a_res(
i) = t_v(
i) - t_dot_u(
i);
923 for (
int rr = 0; rr < 3; rr++) {
924 a_res[rr] >>= res[rr];
929 active.resize(nb_active_vars);
931 for (
int nn1 = 0; nn1 < 3; nn1++) {
935 for (
int nn1 = 0; nn1 < 3; nn1++) {
942 for (
int nn1 = 0; nn1 < 3; nn1++) {
943 for (
int nn2 = 0; nn2 < 3; nn2++) {
944 if (fieldDisp && nn1 == nn2) {
956 for (
int nn1 = 0; nn1 < 3; nn1++) {
963 for (
int nn1 = 0; nn1 < 3; nn1++) {
964 for (
int nn2 = 0; nn2 < 3; nn2++) {
977 r = ::function(
tAg, 3, nb_active_vars, &active[0], &res[0]);
980 "ADOL-C function evaluation with error");
983 double val = getVolume() * getGaussPts()(3, gg);
988 for (
int nn1 = 0; nn1 < 3; nn1++) {
992 r = jacobian(
tAg, 3, nb_active_vars, &active[0],
996 "ADOL-C function evaluation with error");
998 double val = getVolume() * getGaussPts()(3, gg);
1017 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1026 int nb_base_functions = row_data.
getN().size2();
1034 for (
unsigned int gg = 0; gg < row_data.
getN().size1(); gg++) {
1040 for (; dd < nb_dofs / 3; dd++) {
1041 t_nf(
i) += base * res(
i);
1045 for (; dd != nb_base_functions; dd++) {
1050 if (row_data.
getIndices().size() > 3 * row_data.
getN().size2()) {
1051 SETERRQ(PETSC_COMM_SELF, 1,
"data inconsistency");
1054 &row_data.
getIndices()[0], &nf[0], ADD_VALUES);
1071 double *base_ptr =
const_cast<double *
>(&col_data.
getN(gg)[0]);
1074 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1075 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1084 for (
int dd = 0; dd < nb_col / 3; dd++) {
1085 t_jac(
i,
j) += t_mass1(
i,
j) * base;
1105 double *base_ptr =
const_cast<double *
>(&col_data.
getN(gg)[0]);
1108 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1109 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1123 for (
int dd = 0; dd < nb_col / 3; dd++) {
1124 t_jac(
i,
j) += t_mass1(
i,
j) * base * getFEMethod()->ts_a;
1130 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
1132 const int s = 3 + 3;
1157 for (
int dd = 0; dd < nb_col / 3; dd++) {
1158 t_jac(
i,
j) += t_mass1(
i,
j) * base * getFEMethod()->ts_a;
1159 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
1180 double *base_ptr =
const_cast<double *
>(&col_data.
getN(gg)[0]);
1183 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
1186 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1187 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1188 const int u = 3 + 3 + 9;
1195 const int s = 3 + 3 + 9 + 3;
1223 for (
int dd = 0; dd < nb_col / 3; dd++) {
1224 t_jac(
i,
j) += t_mass1(
i,
j) * base * getFEMethod()->ts_a;
1225 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
1240 dAta(data),
commonData(common_data),
tAg(tag), jAcobian(jacobian),
1248 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1254 if (row_type != MBVERTEX)
1273 for (
int dd = 0; dd < 3; dd++) {
1278 int nb_gauss_pts = row_data.
getN().size1();
1283 int nb_active_vars = 0;
1284 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1290 for (
int nn1 = 0; nn1 < 3; nn1++) {
1297 for (
int nn1 = 0; nn1 < 3; nn1++) {
1302 for (
int nn1 = 0; nn1 < 3; nn1++) {
1303 for (
int nn2 = 0; nn2 < 3; nn2++) {
1310 for (
int nn1 = 0; nn1 < 3; nn1++) {
1311 for (
int nn2 = 0; nn2 < 3; nn2++) {
1324 for (
int nn1 = 0; nn1 < 3; nn1++) {
1325 for (
int nn2 = 0; nn2 < 3; nn2++) {
1356 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
1357 t_G(
i,
j) = t_g(
i,
k) * t_invH(
k,
j);
1358 t_a_T(
i) = t_F(
k,
i) * t_a(
k) + t_G(
k,
i) * t_v(
k);
1359 const auto rho0 = dAta.rho0;
1360 t_a_T(
i) = -rho0 * detH;
1363 for (
int nn = 0; nn < 3; nn++) {
1369 active.resize(nb_active_vars);
1371 for (
int nn1 = 0; nn1 < 3; nn1++) {
1377 for (
int nn1 = 0; nn1 < 3; nn1++) {
1381 for (
int nn1 = 0; nn1 < 3; nn1++) {
1382 for (
int nn2 = 0; nn2 < 3; nn2++) {
1388 for (
int nn1 = 0; nn1 < 3; nn1++) {
1389 for (
int nn2 = 0; nn2 < 3; nn2++) {
1390 if (fieldDisp && nn1 == nn2) {
1403 for (
int nn1 = 0; nn1 < 3; nn1++) {
1404 for (
int nn2 = 0; nn2 < 3; nn2++) {
1417 r = ::function(
tAg, 3, nb_active_vars, &active[0], &res[0]);
1420 "ADOL-C function evaluation with error r = %d", r);
1423 double val = getVolume() * getGaussPts()(3, gg);
1428 for (
int nn1 = 0; nn1 < 3; nn1++) {
1432 r = jacobian(
tAg, 3, nb_active_vars, &active[0],
1436 "ADOL-C function evaluation with error");
1438 double val = getVolume() * getGaussPts()(3, gg);
1443 }
catch (
const std::exception &ex) {
1444 std::ostringstream ss;
1445 ss <<
"throw in method: " << ex.what() << std::endl;
1446 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
1456 Range *forcesonlyonentities_ptr)
1460 if (forcesonlyonentities_ptr != NULL) {
1470 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1484 int nb_base_functions = row_data.
getN().size2();
1487 for (
unsigned int gg = 0; gg < row_data.
getN().size1(); gg++) {
1493 for (; dd < nb_dofs / 3; dd++) {
1494 t_nf(
i) += base * res(
i);
1498 for (; dd != nb_base_functions; dd++) {
1503 if (row_data.
getIndices().size() > 3 * row_data.
getN().size2()) {
1504 SETERRQ(PETSC_COMM_SELF, 1,
"data inconsistency");
1506 if (!forcesOnlyOnEntities.empty()) {
1509 VectorDofs::iterator dit = dofs.begin();
1510 for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
1511 if (forcesOnlyOnEntities.find((*dit)->getEnt()) ==
1512 forcesOnlyOnEntities.end()) {
1519 &nf[0], ADD_VALUES);
1522 &row_data.
getIndices()[0], &nf[0], ADD_VALUES);
1525 }
catch (
const std::exception &ex) {
1526 std::ostringstream ss;
1527 ss <<
"throw in method: " << ex.what() << std::endl;
1528 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
1539 Range *forcesonlyonentities_ptr)
1541 vel_field,
field_name, data, common_data, forcesonlyonentities_ptr) {}
1549 double *base_ptr =
const_cast<double *
>(&col_data.
getN(gg)[0]);
1552 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
1555 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1556 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1564 const int s = 3 + 3;
1592 for (
int dd = 0; dd < nb_col / 3; dd++) {
1593 t_jac(
i,
j) += t_mass1(
i,
j) * base * getFEMethod()->ts_a;
1594 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
1607 Range *forcesonlyonentities_ptr)
1609 vel_field,
field_name, data, common_data, forcesonlyonentities_ptr) {}
1618 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
1621 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1622 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1623 const int s = 3 + 3 + 9;
1651 for (
int dd = 0; dd < nb_col / 3; dd++) {
1652 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
1664 Range *forcesonlyonentities_ptr)
1666 vel_field,
field_name, data, common_data, forcesonlyonentities_ptr) {}
1675 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
1678 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1679 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1680 const int s = 3 + 3 + 9 + 9;
1708 for (
int dd = 0; dd < nb_col / 3; dd++) {
1709 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
1718 const std::string spatial_position_field)
1719 :
mField(m_field), tS(_ts), velocityField(velocity_field),
1720 spatialPositionField(spatial_position_field), jacobianLag(-1) {}
1726 case CTX_TSSETIFUNCTION: {
1732 problemPtr,
COL, ts_u, INSERT_VALUES, SCATTER_REVERSE);
1734 problemPtr, velocityField,
"DOT_" + velocityField,
COL, ts_u_t,
1735 INSERT_VALUES, SCATTER_REVERSE);
1737 problemPtr, spatialPositionField,
"DOT_" + spatialPositionField,
COL,
1738 ts_u_t, INSERT_VALUES, SCATTER_REVERSE);
1741 case CTX_TSSETIJACOBIAN: {
1767 int id = it->getMeshsetId();
1773 CHKERR it->getAttributeDataStructure(mydata);
1785 CHKERR it->getAttributeDataStructure(mydata);
1786 if (mydata.
data.User1 == 0)
1791 tets = subtract(tets, added_tets);
1794 int id = it->getMeshsetId();
1809 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr) {
1812 if (!block_sets_ptr)
1814 "Pointer to block of sets is null");
1819 CHKERR it->getAttributeDataStructure(mydata);
1820 int id = it->getMeshsetId();
1821 auto &block_data = (*block_sets_ptr)[id];
1824 block_data.tEts,
true);
1825 block_data.rho0 = mydata.
data.density;
1826 block_data.a0.resize(3);
1827 block_data.a0[0] = mydata.
data.acceleration_x;
1828 block_data.a0[1] = mydata.
data.acceleration_y;
1829 block_data.a0[2] = mydata.
data.acceleration_z;
1836 string element_name,
string velocity_field_name,
1837 string spatial_position_field_name,
string material_position_field_name,
1845 velocity_field_name);
1847 velocity_field_name);
1849 velocity_field_name);
1851 element_name, spatial_position_field_name);
1853 element_name, spatial_position_field_name);
1855 element_name, spatial_position_field_name);
1859 element_name, material_position_field_name);
1861 element_name, material_position_field_name);
1863 element_name,
"DOT_" + material_position_field_name);
1866 element_name, material_position_field_name);
1869 element_name,
"DOT_" + velocity_field_name);
1871 element_name,
"DOT_" + spatial_position_field_name);
1879 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
1881 Range add_tets = sit->second.tEts;
1882 if (!tets.empty()) {
1883 add_tets = intersect(add_tets, tets);
1893 string element_name,
string velocity_field_name,
1894 string spatial_position_field_name,
string material_position_field_name,
1902 velocity_field_name);
1904 velocity_field_name);
1906 velocity_field_name);
1908 element_name, spatial_position_field_name);
1910 element_name, spatial_position_field_name);
1914 element_name, material_position_field_name);
1916 element_name,
"DOT_" + material_position_field_name);
1919 element_name, material_position_field_name);
1922 element_name,
"DOT_" + velocity_field_name);
1924 element_name,
"DOT_" + spatial_position_field_name);
1932 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
1934 Range add_tets = sit->second.tEts;
1935 if (!tets.empty()) {
1936 add_tets = intersect(add_tets, tets);
1946 string element_name,
string velocity_field_name,
1947 string spatial_position_field_name,
string material_position_field_name,
1955 velocity_field_name);
1957 velocity_field_name);
1959 element_name, spatial_position_field_name);
1961 element_name, spatial_position_field_name);
1965 element_name, material_position_field_name);
1967 element_name, material_position_field_name);
1969 element_name,
"DOT_" + material_position_field_name);
1972 element_name, material_position_field_name);
1975 element_name,
"DOT_" + velocity_field_name);
1977 element_name,
"DOT_" + spatial_position_field_name);
1984 if (intersected != NULL) {
1986 tets = *intersected;
1988 tets = intersect(*intersected, tets);
1992 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
1994 Range add_tets = sit->second.tEts;
1995 if (!tets.empty()) {
1996 add_tets = intersect(add_tets, tets);
2006 string velocity_field_name,
string spatial_position_field_name,
2007 string material_position_field_name,
bool ale,
bool linear) {
2023 "DOT_" + spatial_position_field_name,
commonData));
2029 "DOT_" + material_position_field_name,
commonData));
2031 feMassRhs.meshPositionsFieldName = material_position_field_name;
2034 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
2051 "DOT_" + spatial_position_field_name,
commonData));
2057 "DOT_" + material_position_field_name,
commonData));
2059 feMassLhs.meshPositionsFieldName = material_position_field_name;
2071 spatial_position_field_name, spatial_position_field_name, sit->second,
2076 spatial_position_field_name, material_position_field_name,
2079 feMassLhs.meshPositionsFieldName = material_position_field_name;
2085 feEnergy.getOpPtrVector().push_back(
2087 feEnergy.getOpPtrVector().push_back(
2092 feEnergy.meshPositionsFieldName = material_position_field_name;
2104 string velocity_field_name,
string spatial_position_field_name,
2105 string material_position_field_name,
bool ale) {
2113 feVelRhs.getOpPtrVector().push_back(
2115 feVelRhs.getOpPtrVector().push_back(
2117 feVelRhs.getOpPtrVector().push_back(
2121 "DOT_" + spatial_position_field_name,
commonData));
2126 "DOT_" + material_position_field_name,
commonData));
2128 feVelRhs.meshPositionsFieldName = material_position_field_name;
2131 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
2135 feVelRhs.getOpPtrVector().push_back(
2140 feVelLhs.getOpPtrVector().push_back(
2142 feVelLhs.getOpPtrVector().push_back(
2144 feVelLhs.getOpPtrVector().push_back(
2148 "DOT_" + spatial_position_field_name,
commonData));
2153 "DOT_" + material_position_field_name,
commonData));
2155 feVelLhs.meshPositionsFieldName = material_position_field_name;
2163 velocity_field_name, velocity_field_name, sit->second,
commonData));
2165 velocity_field_name, spatial_position_field_name, sit->second,
2170 velocity_field_name, material_position_field_name, sit->second,
2173 feVelLhs.meshPositionsFieldName = material_position_field_name;
2182 string velocity_field_name,
string spatial_position_field_name,
2183 string material_position_field_name,
Range *forces_on_entities_ptr) {
2191 feTRhs.getOpPtrVector().push_back(
2193 feTRhs.getOpPtrVector().push_back(
2195 feTRhs.getOpPtrVector().push_back(
2197 feTRhs.getOpPtrVector().push_back(
2200 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
2202 feTRhs.getOpPtrVector().push_back(
2204 material_position_field_name, sit->second,
commonData,
tAg,
false));
2206 material_position_field_name, sit->second,
commonData,
2207 forces_on_entities_ptr));
2211 feTLhs.getOpPtrVector().push_back(
2213 feTLhs.getOpPtrVector().push_back(
2215 feTLhs.getOpPtrVector().push_back(
2223 feTLhs.getOpPtrVector().push_back(
2225 material_position_field_name, sit->second,
commonData,
tAg));
2226 feTLhs.getOpPtrVector().push_back(
2228 material_position_field_name, velocity_field_name, sit->second,
2230 feTLhs.getOpPtrVector().push_back(
2232 material_position_field_name, spatial_position_field_name,
2233 sit->second,
commonData, forces_on_entities_ptr));
2234 feTLhs.getOpPtrVector().push_back(
2236 material_position_field_name, material_position_field_name,
2237 sit->second,
commonData, forces_on_entities_ptr));
2244 string velocity_field_name,
string spatial_position_field_name,
2245 string material_position_field_name,
bool linear) {
2263 feMassRhs.meshPositionsFieldName = material_position_field_name;
2265 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
2284 feMassLhs.meshPositionsFieldName = material_position_field_name;
2292 spatial_position_field_name, spatial_position_field_name, sit->second,
2295 feMassLhs.meshPositionsFieldName = material_position_field_name;
2309 feMassAuxLhs.meshPositionsFieldName = material_position_field_name;
2317 spatial_position_field_name, spatial_position_field_name, sit->second,
2320 feMassAuxLhs.meshPositionsFieldName = material_position_field_name;
2325 feEnergy.getOpPtrVector().push_back(
2327 feEnergy.getOpPtrVector().push_back(
2332 feEnergy.meshPositionsFieldName = material_position_field_name;
2348 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
2356#if PETSC_VERSION_GE(3, 5, 3)
2357 CHKERR MatCreateVecs(
K, &u, &Ku);
2360 CHKERR MatGetVecs(
K, &u, &Ku);
2363 CHKERR MatDuplicate(
K, MAT_SHARE_NONZERO_PATTERN, &barK);
2377 CHKERR MatDestroy(&barK);
2378 iNitialized =
false;
2388 CHKERR PetscObjectGetComm((PetscObject)shellMat, &comm);
2389 CHKERR PCCreate(comm, &pC);
2412 if (
ts_ctx != CTX_TSSETIFUNCTION) {
2414 "It is used to residual of velocities");
2416 if (!shellMatCtx->iNitialized) {
2417 CHKERR shellMatCtx->iNit();
2420 CHKERR VecScatterBegin(shellMatCtx->scatterU, ts_u_t, shellMatCtx->u,
2421 INSERT_VALUES, SCATTER_FORWARD);
2422 CHKERR VecScatterEnd(shellMatCtx->scatterU, ts_u_t, shellMatCtx->u,
2423 INSERT_VALUES, SCATTER_FORWARD);
2424 CHKERR VecScatterBegin(shellMatCtx->scatterV, ts_u, shellMatCtx->v,
2425 INSERT_VALUES, SCATTER_FORWARD);
2426 CHKERR VecScatterEnd(shellMatCtx->scatterV, ts_u, shellMatCtx->v,
2427 INSERT_VALUES, SCATTER_FORWARD);
2428 CHKERR VecAXPY(shellMatCtx->v, -1, shellMatCtx->u);
2429 CHKERR VecScatterBegin(shellMatCtx->scatterV, shellMatCtx->v, ts_F,
2430 ADD_VALUES, SCATTER_REVERSE);
2431 CHKERR VecScatterEnd(shellMatCtx->scatterV, shellMatCtx->v, ts_F, ADD_VALUES,
2438#ifdef __DIRICHLET_HPP__
2440ConvectiveMassElement::ShellMatrixElement::ShellMatrixElement(
2444MoFEMErrorCode ConvectiveMassElement::ShellMatrixElement::preProcess() {
2447 if (
ts_ctx != CTX_TSSETIJACOBIAN) {
2449 "It is used to calculate shell matrix only");
2452 shellMatCtx->ts_a = ts_a;
2453 DirichletBcPtr->copyTs(*((
TSMethod *)
this));
2455 DirichletBcPtr->dIag = 1;
2456 DirichletBcPtr->ts_B = shellMatCtx->K;
2457 CHKERR MatZeroEntries(shellMatCtx->K);
2459 LoopsToDoType::iterator itk = loopK.begin();
2460 for (; itk != loopK.end(); itk++) {
2461 itk->second->copyTs(*((
TSMethod *)
this));
2462 itk->second->ts_B = shellMatCtx->K;
2465 LoopsToDoType::iterator itam = loopAuxM.begin();
2466 for (; itam != loopAuxM.end(); itam++) {
2467 itam->second->copyTs(*((
TSMethod *)
this));
2468 itam->second->ts_B = shellMatCtx->K;
2472 CHKERR MatAssemblyBegin(shellMatCtx->K, MAT_FINAL_ASSEMBLY);
2473 CHKERR MatAssemblyEnd(shellMatCtx->K, MAT_FINAL_ASSEMBLY);
2475 DirichletBcPtr->dIag = 0;
2476 DirichletBcPtr->ts_B = shellMatCtx->M;
2477 CHKERR MatZeroEntries(shellMatCtx->M);
2479 LoopsToDoType::iterator itm = loopM.begin();
2480 for (; itm != loopM.end(); itm++) {
2481 itm->second->copyTs(*((
TSMethod *)
this));
2482 itm->second->ts_B = shellMatCtx->M;
2486 CHKERR MatAssemblyBegin(shellMatCtx->M, MAT_FINAL_ASSEMBLY);
2487 CHKERR MatAssemblyEnd(shellMatCtx->M, MAT_FINAL_ASSEMBLY);
2490 CHKERR MatZeroEntries(shellMatCtx->barK);
2491 CHKERR MatCopy(shellMatCtx->K, shellMatCtx->barK, SAME_NONZERO_PATTERN);
2492 CHKERR MatAXPY(shellMatCtx->barK, ts_a, shellMatCtx->M, SAME_NONZERO_PATTERN);
2493 CHKERR MatAssemblyBegin(shellMatCtx->barK, MAT_FINAL_ASSEMBLY);
2494 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 add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
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 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_row(const std::string &fe_name, const std::string name_row)=0
set field row 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
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
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.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr, const size_t cc=0)
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
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.
intrusive_ptr for managing petsc objects
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.