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