12#include <adolc/adolc.h>
13#include <MethodForForceScaling.hpp>
15#include <MethodForForceScaling.hpp>
19#error "MoFEM need to be compiled with ADOL-C"
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);
522 int row_side,
int col_side, EntityType row_type, EntityType col_type,
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;
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),
801 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
807 if (row_type != MBVERTEX)
828 for (
int dd = 0; dd < 3; dd++) {
834 int nb_gauss_pts = row_data.
getN().size1();
839 int nb_active_vars = 0;
840 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
846 for (
int nn1 = 0; nn1 < 3; nn1++) {
851 for (
int nn1 = 0; nn1 < 3; nn1++) {
859 for (
int nn1 = 0; nn1 < 3; nn1++) {
860 for (
int nn2 = 0; nn2 < 3; nn2++) {
872 for (
int nn1 = 0; nn1 < 3; nn1++) {
880 for (
int nn1 = 0; nn1 < 3; nn1++) {
881 for (
int nn2 = 0; nn2 < 3; nn2++) {
912 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
914 t_F(
i,
j) = t_h(
i,
j);
917 t_dot_u(
i) = t_dot_w(
i) + t_F(
i,
j) * t_dot_W(
j);
918 t_a_res(
i) = t_v(
i) - t_dot_u(
i);
924 for (
int rr = 0; rr < 3; rr++) {
925 a_res[rr] >>= res[rr];
930 active.resize(nb_active_vars);
932 for (
int nn1 = 0; nn1 < 3; nn1++) {
936 for (
int nn1 = 0; nn1 < 3; nn1++) {
943 for (
int nn1 = 0; nn1 < 3; nn1++) {
944 for (
int nn2 = 0; nn2 < 3; nn2++) {
945 if (fieldDisp && nn1 == nn2) {
957 for (
int nn1 = 0; nn1 < 3; nn1++) {
964 for (
int nn1 = 0; nn1 < 3; nn1++) {
965 for (
int nn2 = 0; nn2 < 3; nn2++) {
978 r = ::function(
tAg, 3, nb_active_vars, &active[0], &res[0]);
981 "ADOL-C function evaluation with error");
984 double val = getVolume() * getGaussPts()(3, gg);
989 for (
int nn1 = 0; nn1 < 3; nn1++) {
993 r = jacobian(
tAg, 3, nb_active_vars, &active[0],
997 "ADOL-C function evaluation with error");
999 double val = getVolume() * getGaussPts()(3, gg);
1018 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1027 int nb_base_functions = row_data.
getN().size2();
1035 for (
unsigned int gg = 0; gg < row_data.
getN().size1(); gg++) {
1041 for (; dd < nb_dofs / 3; dd++) {
1042 t_nf(
i) += base * res(
i);
1046 for (; dd != nb_base_functions; dd++) {
1051 if (row_data.
getIndices().size() > 3 * row_data.
getN().size2()) {
1052 SETERRQ(PETSC_COMM_SELF, 1,
"data inconsistency");
1055 &row_data.
getIndices()[0], &nf[0], ADD_VALUES);
1072 double *base_ptr =
const_cast<double *
>(&col_data.
getN(gg)[0]);
1075 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1076 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1085 for (
int dd = 0; dd < nb_col / 3; dd++) {
1086 t_jac(
i,
j) += t_mass1(
i,
j) * base;
1106 double *base_ptr =
const_cast<double *
>(&col_data.
getN(gg)[0]);
1109 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1110 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1124 for (
int dd = 0; dd < nb_col / 3; dd++) {
1125 t_jac(
i,
j) += t_mass1(
i,
j) * base * getFEMethod()->ts_a;
1131 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
1133 const int s = 3 + 3;
1158 for (
int dd = 0; dd < nb_col / 3; dd++) {
1159 t_jac(
i,
j) += t_mass1(
i,
j) * base * getFEMethod()->ts_a;
1160 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
1181 double *base_ptr =
const_cast<double *
>(&col_data.
getN(gg)[0]);
1184 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
1187 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1188 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1189 const int u = 3 + 3 + 9;
1196 const int s = 3 + 3 + 9 + 3;
1224 for (
int dd = 0; dd < nb_col / 3; dd++) {
1225 t_jac(
i,
j) += t_mass1(
i,
j) * base * getFEMethod()->ts_a;
1226 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
1241 dAta(data),
commonData(common_data),
tAg(tag), jAcobian(jacobian),
1249 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1255 if (row_type != MBVERTEX)
1274 for (
int dd = 0; dd < 3; dd++) {
1279 int nb_gauss_pts = row_data.
getN().size1();
1284 int nb_active_vars = 0;
1285 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1291 for (
int nn1 = 0; nn1 < 3; nn1++) {
1298 for (
int nn1 = 0; nn1 < 3; nn1++) {
1303 for (
int nn1 = 0; nn1 < 3; nn1++) {
1304 for (
int nn2 = 0; nn2 < 3; nn2++) {
1311 for (
int nn1 = 0; nn1 < 3; nn1++) {
1312 for (
int nn2 = 0; nn2 < 3; nn2++) {
1325 for (
int nn1 = 0; nn1 < 3; nn1++) {
1326 for (
int nn2 = 0; nn2 < 3; nn2++) {
1357 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
1358 t_G(
i,
j) = t_g(
i,
k) * t_invH(
k,
j);
1359 t_a_T(
i) = t_F(
k,
i) * t_a(
k) + t_G(
k,
i) * t_v(
k);
1360 const auto rho0 = dAta.rho0;
1361 t_a_T(
i) = -rho0 * detH;
1364 for (
int nn = 0; nn < 3; nn++) {
1370 active.resize(nb_active_vars);
1372 for (
int nn1 = 0; nn1 < 3; nn1++) {
1378 for (
int nn1 = 0; nn1 < 3; nn1++) {
1382 for (
int nn1 = 0; nn1 < 3; nn1++) {
1383 for (
int nn2 = 0; nn2 < 3; nn2++) {
1389 for (
int nn1 = 0; nn1 < 3; nn1++) {
1390 for (
int nn2 = 0; nn2 < 3; nn2++) {
1391 if (fieldDisp && nn1 == nn2) {
1404 for (
int nn1 = 0; nn1 < 3; nn1++) {
1405 for (
int nn2 = 0; nn2 < 3; nn2++) {
1418 r = ::function(
tAg, 3, nb_active_vars, &active[0], &res[0]);
1421 "ADOL-C function evaluation with error r = %d", r);
1424 double val = getVolume() * getGaussPts()(3, gg);
1429 for (
int nn1 = 0; nn1 < 3; nn1++) {
1433 r = jacobian(
tAg, 3, nb_active_vars, &active[0],
1437 "ADOL-C function evaluation with error");
1439 double val = getVolume() * getGaussPts()(3, gg);
1444 }
catch (
const std::exception &ex) {
1445 std::ostringstream ss;
1446 ss <<
"throw in method: " << ex.what() << std::endl;
1457 Range *forcesonlyonentities_ptr)
1461 if (forcesonlyonentities_ptr != NULL) {
1471 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1485 int nb_base_functions = row_data.
getN().size2();
1488 for (
unsigned int gg = 0; gg < row_data.
getN().size1(); gg++) {
1494 for (; dd < nb_dofs / 3; dd++) {
1495 t_nf(
i) += base * res(
i);
1499 for (; dd != nb_base_functions; dd++) {
1504 if (row_data.
getIndices().size() > 3 * row_data.
getN().size2()) {
1505 SETERRQ(PETSC_COMM_SELF, 1,
"data inconsistency");
1507 if (!forcesOnlyOnEntities.empty()) {
1510 VectorDofs::iterator dit = dofs.begin();
1511 for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
1512 if (forcesOnlyOnEntities.find((*dit)->getEnt()) ==
1513 forcesOnlyOnEntities.end()) {
1520 &nf[0], ADD_VALUES);
1523 &row_data.
getIndices()[0], &nf[0], ADD_VALUES);
1526 }
catch (
const std::exception &ex) {
1527 std::ostringstream ss;
1528 ss <<
"throw in method: " << ex.what() << std::endl;
1540 Range *forcesonlyonentities_ptr)
1542 vel_field,
field_name, data, common_data, forcesonlyonentities_ptr) {}
1550 double *base_ptr =
const_cast<double *
>(&col_data.
getN(gg)[0]);
1553 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
1556 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1557 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1565 const int s = 3 + 3;
1593 for (
int dd = 0; dd < nb_col / 3; dd++) {
1594 t_jac(
i,
j) += t_mass1(
i,
j) * base * getFEMethod()->ts_a;
1595 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
1608 Range *forcesonlyonentities_ptr)
1610 vel_field,
field_name, data, common_data, forcesonlyonentities_ptr) {}
1619 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
1622 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1623 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1624 const int s = 3 + 3 + 9;
1652 for (
int dd = 0; dd < nb_col / 3; dd++) {
1653 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
1665 Range *forcesonlyonentities_ptr)
1667 vel_field,
field_name, data, common_data, forcesonlyonentities_ptr) {}
1676 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
1679 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1680 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1681 const int s = 3 + 3 + 9 + 9;
1709 for (
int dd = 0; dd < nb_col / 3; dd++) {
1710 t_jac(
i,
j) += t_mass3(
i,
j,
k) * diff(
k);
1719 const std::string spatial_position_field)
1720 :
mField(m_field), tS(_ts), velocityField(velocity_field),
1721 spatialPositionField(spatial_position_field), jacobianLag(-1) {}
1727 case CTX_TSSETIFUNCTION: {
1733 problemPtr,
COL, ts_u, INSERT_VALUES, SCATTER_REVERSE);
1735 problemPtr, velocityField,
"DOT_" + velocityField,
COL, ts_u_t,
1736 INSERT_VALUES, SCATTER_REVERSE);
1738 problemPtr, spatialPositionField,
"DOT_" + spatialPositionField,
COL,
1739 ts_u_t, INSERT_VALUES, SCATTER_REVERSE);
1742 case CTX_TSSETIJACOBIAN: {
1768 int id = it->getMeshsetId();
1774 CHKERR it->getAttributeDataStructure(mydata);
1786 CHKERR it->getAttributeDataStructure(mydata);
1787 if (mydata.
data.User1 == 0)
1792 tets = subtract(tets, added_tets);
1795 int id = it->getMeshsetId();
1810 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr) {
1813 if (!block_sets_ptr)
1815 "Pointer to block of sets is null");
1820 CHKERR it->getAttributeDataStructure(mydata);
1821 int id = it->getMeshsetId();
1822 auto &block_data = (*block_sets_ptr)[id];
1825 block_data.tEts,
true);
1826 block_data.rho0 = mydata.
data.density;
1827 block_data.a0.resize(3);
1828 block_data.a0[0] = mydata.
data.acceleration_x;
1829 block_data.a0[1] = mydata.
data.acceleration_y;
1830 block_data.a0[2] = mydata.
data.acceleration_z;
1837 string element_name,
string velocity_field_name,
1838 string spatial_position_field_name,
string material_position_field_name,
1846 velocity_field_name);
1848 velocity_field_name);
1850 velocity_field_name);
1852 element_name, spatial_position_field_name);
1854 element_name, spatial_position_field_name);
1856 element_name, spatial_position_field_name);
1860 element_name, material_position_field_name);
1862 element_name, material_position_field_name);
1864 element_name,
"DOT_" + material_position_field_name);
1867 element_name, material_position_field_name);
1870 element_name,
"DOT_" + velocity_field_name);
1872 element_name,
"DOT_" + spatial_position_field_name);
1880 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
1882 Range add_tets = sit->second.tEts;
1883 if (!tets.empty()) {
1884 add_tets = intersect(add_tets, tets);
1894 string element_name,
string velocity_field_name,
1895 string spatial_position_field_name,
string material_position_field_name,
1903 velocity_field_name);
1905 velocity_field_name);
1907 velocity_field_name);
1909 element_name, spatial_position_field_name);
1911 element_name, spatial_position_field_name);
1915 element_name, material_position_field_name);
1917 element_name,
"DOT_" + material_position_field_name);
1920 element_name, material_position_field_name);
1923 element_name,
"DOT_" + velocity_field_name);
1925 element_name,
"DOT_" + spatial_position_field_name);
1933 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
1935 Range add_tets = sit->second.tEts;
1936 if (!tets.empty()) {
1937 add_tets = intersect(add_tets, tets);
1947 string element_name,
string velocity_field_name,
1948 string spatial_position_field_name,
string material_position_field_name,
1956 velocity_field_name);
1958 velocity_field_name);
1960 element_name, spatial_position_field_name);
1962 element_name, spatial_position_field_name);
1966 element_name, material_position_field_name);
1968 element_name, material_position_field_name);
1970 element_name,
"DOT_" + material_position_field_name);
1973 element_name, material_position_field_name);
1976 element_name,
"DOT_" + velocity_field_name);
1978 element_name,
"DOT_" + spatial_position_field_name);
1985 if (intersected != NULL) {
1987 tets = *intersected;
1989 tets = intersect(*intersected, tets);
1993 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
1995 Range add_tets = sit->second.tEts;
1996 if (!tets.empty()) {
1997 add_tets = intersect(add_tets, tets);
2007 string velocity_field_name,
string spatial_position_field_name,
2008 string material_position_field_name,
bool ale,
bool linear) {
2024 "DOT_" + spatial_position_field_name,
commonData));
2030 "DOT_" + material_position_field_name,
commonData));
2035 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
2052 "DOT_" + spatial_position_field_name,
commonData));
2058 "DOT_" + material_position_field_name,
commonData));
2072 spatial_position_field_name, spatial_position_field_name, sit->second,
2077 spatial_position_field_name, material_position_field_name,
2105 string velocity_field_name,
string spatial_position_field_name,
2106 string material_position_field_name,
bool ale) {
2122 "DOT_" + spatial_position_field_name,
commonData));
2127 "DOT_" + material_position_field_name,
commonData));
2132 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
2149 "DOT_" + spatial_position_field_name,
commonData));
2154 "DOT_" + material_position_field_name,
commonData));
2164 velocity_field_name, velocity_field_name, sit->second,
commonData));
2166 velocity_field_name, spatial_position_field_name, sit->second,
2171 velocity_field_name, material_position_field_name, sit->second,
2183 string velocity_field_name,
string spatial_position_field_name,
2184 string material_position_field_name,
Range *forces_on_entities_ptr) {
2201 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
2205 material_position_field_name, sit->second,
commonData,
tAg,
false));
2207 material_position_field_name, sit->second,
commonData,
2208 forces_on_entities_ptr));
2226 material_position_field_name, sit->second,
commonData,
tAg));
2229 material_position_field_name, velocity_field_name, sit->second,
2233 material_position_field_name, spatial_position_field_name,
2234 sit->second,
commonData, forces_on_entities_ptr));
2237 material_position_field_name, material_position_field_name,
2238 sit->second,
commonData, forces_on_entities_ptr));
2245 string velocity_field_name,
string spatial_position_field_name,
2246 string material_position_field_name,
bool linear) {
2266 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
2293 spatial_position_field_name, spatial_position_field_name, sit->second,
2318 spatial_position_field_name, spatial_position_field_name, sit->second,
2349 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
2357#if PETSC_VERSION_GE(3, 5, 3)
2358 CHKERR MatCreateVecs(
K, &u, &Ku);
2359 CHKERR MatCreateVecs(M, &
v, &Mv);
2361 CHKERR MatGetVecs(
K, &u, &Ku);
2362 CHKERR MatGetVecs(M, &
v, &Mv);
2364 CHKERR MatDuplicate(
K, MAT_SHARE_NONZERO_PATTERN, &barK);
2378 CHKERR MatDestroy(&barK);
2379 iNitialized =
false;
2389 CHKERR PetscObjectGetComm((PetscObject)shellMat, &comm);
2390 CHKERR PCCreate(comm, &pC);
2413 if (
ts_ctx != CTX_TSSETIFUNCTION) {
2415 "It is used to residual of velocities");
2417 if (!shellMatCtx->iNitialized) {
2418 CHKERR shellMatCtx->iNit();
2421 CHKERR VecScatterBegin(shellMatCtx->scatterU, ts_u_t, shellMatCtx->u,
2422 INSERT_VALUES, SCATTER_FORWARD);
2423 CHKERR VecScatterEnd(shellMatCtx->scatterU, ts_u_t, shellMatCtx->u,
2424 INSERT_VALUES, SCATTER_FORWARD);
2425 CHKERR VecScatterBegin(shellMatCtx->scatterV, ts_u, shellMatCtx->v,
2426 INSERT_VALUES, SCATTER_FORWARD);
2427 CHKERR VecScatterEnd(shellMatCtx->scatterV, ts_u, shellMatCtx->v,
2428 INSERT_VALUES, SCATTER_FORWARD);
2429 CHKERR VecAXPY(shellMatCtx->v, -1, shellMatCtx->u);
2430 CHKERR VecScatterBegin(shellMatCtx->scatterV, shellMatCtx->v, ts_F,
2431 ADD_VALUES, SCATTER_REVERSE);
2432 CHKERR VecScatterEnd(shellMatCtx->scatterV, shellMatCtx->v, ts_F, ADD_VALUES,
2439#ifdef __DIRICHLET_HPP__
2441ConvectiveMassElement::ShellMatrixElement::ShellMatrixElement(
2445MoFEMErrorCode ConvectiveMassElement::ShellMatrixElement::preProcess() {
2448 if (
ts_ctx != CTX_TSSETIJACOBIAN) {
2450 "It is used to calculate shell matrix only");
2453 shellMatCtx->ts_a = ts_a;
2454 DirichletBcPtr->copyTs(*((
TSMethod *)
this));
2456 DirichletBcPtr->dIag = 1;
2457 DirichletBcPtr->ts_B = shellMatCtx->K;
2458 CHKERR MatZeroEntries(shellMatCtx->K);
2460 LoopsToDoType::iterator itk = loopK.begin();
2461 for (; itk != loopK.end(); itk++) {
2462 itk->second->copyTs(*((
TSMethod *)
this));
2463 itk->second->ts_B = shellMatCtx->K;
2466 LoopsToDoType::iterator itam = loopAuxM.begin();
2467 for (; itam != loopAuxM.end(); itam++) {
2468 itam->second->copyTs(*((
TSMethod *)
this));
2469 itam->second->ts_B = shellMatCtx->K;
2473 CHKERR MatAssemblyBegin(shellMatCtx->K, MAT_FINAL_ASSEMBLY);
2474 CHKERR MatAssemblyEnd(shellMatCtx->K, MAT_FINAL_ASSEMBLY);
2476 DirichletBcPtr->dIag = 0;
2477 DirichletBcPtr->ts_B = shellMatCtx->M;
2478 CHKERR MatZeroEntries(shellMatCtx->M);
2480 LoopsToDoType::iterator itm = loopM.begin();
2481 for (; itm != loopM.end(); itm++) {
2482 itm->second->copyTs(*((
TSMethod *)
this));
2483 itm->second->ts_B = shellMatCtx->M;
2487 CHKERR MatAssemblyBegin(shellMatCtx->M, MAT_FINAL_ASSEMBLY);
2488 CHKERR MatAssemblyEnd(shellMatCtx->M, MAT_FINAL_ASSEMBLY);
2491 CHKERR MatZeroEntries(shellMatCtx->barK);
2492 CHKERR MatCopy(shellMatCtx->K, shellMatCtx->barK, SAME_NONZERO_PATTERN);
2493 CHKERR MatAXPY(shellMatCtx->barK, ts_a, shellMatCtx->M, SAME_NONZERO_PATTERN);
2494 CHKERR MatAssemblyBegin(shellMatCtx->barK, MAT_FINAL_ASSEMBLY);
2495 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_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
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.
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
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.
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()
function is run at the end of loop
int getRule(int order)
it is used to calculate nb. of Gauss integration points
MoFEMErrorCode preProcess()
function is run at the beginning of loop
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()
function is run at the end of loop
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 field data coefficients.
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
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Elastic material data structure.
intrusive_ptr for managing petsc objects
data structure for TS (time stepping) context
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Volume finite element base.
std::string meshPositionsFieldName