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)
constexpr IntegrationType G
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()
Post-processing function executed at loop completion.
int getRule(int order)
it is used to calculate nb. of Gauss integration points
MoFEMErrorCode preProcess()
Pre-processing function executed at loop initialization.
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()
Post-processing function executed at loop completion.
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 DOF values on entity.
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 DOF data structures (const version)
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
structure to get information from 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