14#include <adolc/adolc.h> 
   21  auto create_vec = [&]() {
 
   22    constexpr int ghosts[] = {0};
 
 
   34  return 2 * (
order - 1) + addToRule;
 
 
   40  CHKERR VolumeElementForcesAndSourcesCore::preProcess();
 
   42  if (
A != PETSC_NULLPTR) {
 
   46  if (
F != PETSC_NULLPTR) {
 
 
   66    CHKERR VecAssemblyBegin(V);
 
   74  CHKERR VolumeElementForcesAndSourcesCore::postProcess();
 
 
   86    std::vector<VectorDouble> &values_at_gauss_pts,
 
   87    std::vector<MatrixDouble> &gardient_at_gauss_pts)
 
   90      valuesAtGaussPts(values_at_gauss_pts),
 
   91      gradientAtGaussPts(gardient_at_gauss_pts), zeroAtType(MBVERTEX) {}
 
 
   98  const int nb_base_functions = data.
getN().size2();
 
  102  const int nb_gauss_pts = data.
getN().size1();
 
  103  const int rank = data.
getFieldDofs()[0]->getNbOfCoeffs();
 
  106  if (type == zeroAtType) {
 
  107    valuesAtGaussPts.resize(nb_gauss_pts);
 
  108    gradientAtGaussPts.resize(nb_gauss_pts);
 
  109    for (
int gg = 0; gg != nb_gauss_pts; gg++) {
 
  110      valuesAtGaussPts[gg].resize(rank, 
false);
 
  111      gradientAtGaussPts[gg].resize(rank, 3, 
false);
 
  113    for (
int gg = 0; gg != nb_gauss_pts; gg++) {
 
  114      valuesAtGaussPts[gg].clear();
 
  115      gradientAtGaussPts[gg].clear();
 
  126    for (
int gg = 0; gg != nb_gauss_pts; gg++) {
 
  128      double &val = valuesAtGaussPts[gg][0];
 
  130                                         &gradientAtGaussPts[gg](0, 1),
 
  131                                         &gradientAtGaussPts[gg](0, 2));
 
  133      for (; bb != nb_dofs; bb++) {
 
  134        val += base_function * field_data;
 
  135        grad(
i) += diff_base_functions(
i) * field_data;
 
  136        ++diff_base_functions;
 
  140      for (; bb != nb_base_functions; bb++) {
 
  141        ++diff_base_functions;
 
  146  } 
else if (rank == 3) {
 
  148    for (
int gg = 0; gg != nb_gauss_pts; gg++) {
 
  151                                           &valuesAtGaussPts[gg][1],
 
  152                                           &valuesAtGaussPts[gg][2]);
 
  154          &gradientAtGaussPts[gg](0, 0), &gradientAtGaussPts[gg](0, 1),
 
  155          &gradientAtGaussPts[gg](0, 2), &gradientAtGaussPts[gg](1, 0),
 
  156          &gradientAtGaussPts[gg](1, 1), &gradientAtGaussPts[gg](1, 2),
 
  157          &gradientAtGaussPts[gg](2, 0), &gradientAtGaussPts[gg](2, 1),
 
  158          &gradientAtGaussPts[gg](2, 2));
 
  160      for (; bb != nb_dofs / 3; bb++) {
 
  161        values(
i) += base_function * field_data(
i);
 
  162        gradient(
i, 
j) += field_data(
i) * diff_base_functions(
j);
 
  163        ++diff_base_functions;
 
  167      for (; bb != nb_base_functions; bb++) {
 
  168        ++diff_base_functions;
 
  176    for (
int gg = 0; gg < nb_gauss_pts; gg++) {
 
  179      for (
int dd = 0; dd < nb_dofs / rank; dd++) {
 
  180        for (
int rr1 = 0; rr1 < rank; rr1++) {
 
  181          valuesAtGaussPts[gg][rr1] += 
N[dd] * values[rank * dd + rr1];
 
  182          for (
int rr2 = 0; rr2 < 3; rr2++) {
 
  183            gradientAtGaussPts[gg](rr1, rr2) +=
 
  184                diffN(dd, rr2) * values[rank * dd + rr1];
 
 
  202                                   int tag, 
bool jacobian, 
bool ale,
 
  206      dAta(data), 
commonData(common_data), 
tAg(tag), adlocReturnValue(0),
 
  207      jAcobian(jacobian), fUnction(!jacobian), aLe(ale), fieldDisp(field_disp) {
 
 
  216  CHKERR dAta.materialAdoublePtr->calculateP_PiolaKirchhoffI(
 
  217      dAta, getNumeredEntFiniteElementPtr());
 
  220    auto &t_P = dAta.materialAdoublePtr->t_P;
 
  221    auto &t_invH = dAta.materialAdoublePtr->t_invH;
 
  222    t_P(
i, 
j) = t_P(
i, 
k) * t_invH(
j, 
k);
 
  223    t_P(
i, 
j) *= dAta.materialAdoublePtr->detH;
 
  227  for (
int dd1 = 0; dd1 < 3; dd1++) {
 
  228    for (
int dd2 = 0; dd2 < 3; dd2++) {
 
  229      dAta.materialAdoublePtr->P(dd1, dd2) >>=
 
 
  244  dAta.materialAdoublePtr->F.resize(3, 3, 
false);
 
  248    nbActiveVariables = 0;
 
  249    for (
int dd1 = 0; dd1 < 3; dd1++) {
 
  250      for (
int dd2 = 0; dd2 < 3; dd2++) {
 
  251        dAta.materialAdoublePtr->F(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
 
  254            dAta.materialAdoublePtr->F(dd1, dd2) += 1;
 
  263    nbActiveVariables = 0;
 
  265    dAta.materialAdoublePtr->h.resize(3, 3, 
false);
 
  266    for (
int dd1 = 0; dd1 < 3; dd1++) {
 
  267      for (
int dd2 = 0; dd2 < 3; dd2++) {
 
  268        dAta.materialAdoublePtr->h(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
 
  273    dAta.materialAdoublePtr->H.resize(3, 3, 
false);
 
  274    for (
int dd1 = 0; dd1 < 3; dd1++) {
 
  275      for (
int dd2 = 0; dd2 < 3; dd2++) {
 
  276        dAta.materialAdoublePtr->H(dd1, dd2) <<= (*ptrH)[gg](dd1, dd2);
 
  282    dAta.materialAdoublePtr->invH.resize(3, 3, 
false);
 
  284                            dAta.materialAdoublePtr->detH,
 
  285                            dAta.materialAdoublePtr->invH);
 
  287    auto &t_F = dAta.materialAdoublePtr->t_F;
 
  288    auto &t_h = dAta.materialAdoublePtr->t_h;
 
  289    auto &t_invH = dAta.materialAdoublePtr->t_invH;
 
  291    t_F(
i, 
j) = t_h(
i, 
k) * t_invH(
k, 
j);
 
  295  CHKERR dAta.materialAdoublePtr->setUserActiveVariables(nbActiveVariables);
 
  296  CHKERR calculateStress(gg);
 
 
  312    r = ::function(
tAg, 9, nbActiveVariables, &activeVariables[0],
 
  314    if (r < adlocReturnValue) { 
 
  316              "ADOL-C function evaluation with error r = %d", r);
 
  322    double *jac_ptr[] = {
 
  329    r = jacobian(
tAg, 9, nbActiveVariables, &activeVariables[0], jac_ptr);
 
  330    if (r < adlocReturnValue) {
 
  332              "ADOL-C function evaluation with error");
 
 
  340    int row_side, EntityType row_type,
 
  345  if (row_type != MBVERTEX)
 
  348  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
 
  356  dAta.materialAdoublePtr->commonDataPtr = &
commonData;
 
  357  dAta.materialAdoublePtr->opPtr = 
this;
 
  359  int nb_gauss_pts = row_data.
getN().size1();
 
  368  for (
int gg = 0; gg != nb_gauss_pts; gg++) {
 
  370    dAta.materialAdoublePtr->gG = gg;
 
  373    if (recordTagForIntegrationPoint(gg)) {
 
  378    if (jAcobian || (!recordTagForIntegrationPoint(gg))) {
 
  379      activeVariables.resize(nbActiveVariables, 
false);
 
  381        for (
int dd1 = 0; dd1 < 3; dd1++) {
 
  382          for (
int dd2 = 0; dd2 < 3; dd2++) {
 
  383            activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
 
  387        for (
int dd1 = 0; dd1 < 3; dd1++) {
 
  388          for (
int dd2 = 0; dd2 < 3; dd2++) {
 
  389            activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
 
  392        for (
int dd1 = 0; dd1 < 3; dd1++) {
 
  393          for (
int dd2 = 0; dd2 < 3; dd2++) {
 
  394            activeVariables(9 + dd1 * 3 + dd2) = (*ptrH)[gg](dd1, dd2);
 
  398      CHKERR dAta.materialAdoublePtr->setUserActiveVariables(activeVariables);
 
  401      if (jAcobian || (!recordTagForIntegrationPoint(gg))) {
 
 
  414    bool hessian, 
bool ale, 
bool field_disp)
 
  417      dAta(data), 
commonData(common_data), 
tAg(tag), gRadient(gradient),
 
  418      hEssian(hessian), aLe(ale), fieldDisp(field_disp) {}
 
 
  423  CHKERR dAta.materialAdoublePtr->calculateElasticEnergy(
 
  424      dAta, getNumeredEntFiniteElementPtr());
 
 
  437    nbActiveVariables = 0;
 
  438    for (
int dd1 = 0; dd1 < 3; dd1++) {
 
  439      for (
int dd2 = 0; dd2 < 3; dd2++) {
 
  440        dAta.materialAdoublePtr->F(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
 
  443            dAta.materialAdoublePtr->F(dd1, dd2) += 1;
 
  452    nbActiveVariables = 0;
 
  454    dAta.materialAdoublePtr->h.resize(3, 3, 
false);
 
  455    for (
int dd1 = 0; dd1 < 3; dd1++) {
 
  456      for (
int dd2 = 0; dd2 < 3; dd2++) {
 
  457        dAta.materialAdoublePtr->h(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
 
  462    dAta.materialAdoublePtr->H.resize(3, 3, 
false);
 
  463    for (
int dd1 = 0; dd1 < 3; dd1++) {
 
  464      for (
int dd2 = 0; dd2 < 3; dd2++) {
 
  465        dAta.materialAdoublePtr->H(dd1, dd2) <<= (*ptrH)[gg](dd1, dd2);
 
  471    dAta.materialAdoublePtr->invH.resize(3, 3, 
false);
 
  473                            dAta.materialAdoublePtr->detH,
 
  474                            dAta.materialAdoublePtr->invH);
 
  476    auto &t_F = dAta.materialAdoublePtr->t_F;
 
  477    auto &t_h = dAta.materialAdoublePtr->t_h;
 
  478    auto &t_invH = dAta.materialAdoublePtr->t_invH;
 
  480    t_F(
i, 
j) = t_h(
i, 
k) * t_invH(
k, 
j);
 
  484  CHKERR dAta.materialAdoublePtr->setUserActiveVariables(nbActiveVariables);
 
 
  498    int r = ::gradient(
tAg, nbActiveVariables, &activeVariables[0],
 
  504              "ADOL-C function evaluation with error");
 
  511    double *
H[nbActiveVariables];
 
  512    for (
int n = 0; 
n != nbActiveVariables; 
n++) {
 
  515    int r = ::hessian(
tAg, nbActiveVariables, &*activeVariables.begin(), 
H);
 
  520              "ADOL-C function evaluation with error");
 
 
  528    int row_side, EntityType row_type,
 
  533  if (row_type != MBVERTEX)
 
  536  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
 
  544  dAta.materialAdoublePtr->commonDataPtr = &
commonData;
 
  545  dAta.materialAdoublePtr->opPtr = 
this;
 
  547  int nb_gauss_pts = row_data.
getN().size1();
 
  556  for (
int gg = 0; gg != nb_gauss_pts; gg++) {
 
  558    dAta.materialAdoublePtr->gG = gg;
 
  561    if (recordTagForIntegrationPoint(gg)) {
 
  565    activeVariables.resize(nbActiveVariables, 
false);
 
  567      for (
int dd1 = 0; dd1 < 3; dd1++) {
 
  568        for (
int dd2 = 0; dd2 < 3; dd2++) {
 
  569          activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
 
  573      for (
int dd1 = 0; dd1 < 3; dd1++) {
 
  574        for (
int dd2 = 0; dd2 < 3; dd2++) {
 
  575          activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
 
  578      for (
int dd1 = 0; dd1 < 3; dd1++) {
 
  579        for (
int dd2 = 0; dd2 < 3; dd2++) {
 
  580          activeVariables(9 + dd1 * 3 + dd2) = (*ptrH)[gg](dd1, dd2);
 
  584    CHKERR dAta.materialAdoublePtr->setUserActiveVariables(activeVariables);
 
 
  597      dAta(data), 
commonData(common_data), aLe(false) {}
 
 
  600    int row_side, EntityType row_type,
 
  606  if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
 
  607    iNdices.resize(nb_dofs, 
false);
 
  609    indices_ptr = &iNdices[0];
 
  611    VectorDofs::iterator dit = dofs.begin();
 
  612    for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
 
  613      if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
 
  614          dAta.forcesOnlyOnEntitiesRow.end()) {
 
 
  625    int row_side, EntityType row_type,
 
  629  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
 
  634  const int nb_dofs = row_data.
getIndices().size();
 
  637  if ((
unsigned int)nb_dofs > 3 * row_data.
getN().size2()) {
 
  638    SETERRQ(PETSC_COMM_SELF, 1, 
"data inconsistency");
 
  640  const int nb_base_functions = row_data.
getN().size2();
 
  641  const int nb_gauss_pts = row_data.
getN().size1();
 
  643  nf.resize(nb_dofs, 
false);
 
  650  for (
int gg = 0; gg != nb_gauss_pts; gg++) {
 
  651    double val = getVolume() * getGaussPts()(3, gg);
 
  654        &stress(0, 0), &stress(0, 1), &stress(0, 2), &stress(1, 0),
 
  655        &stress(1, 1), &stress(1, 2), &stress(2, 0), &stress(2, 1),
 
  659    for (; bb != nb_dofs / 3; bb++) {
 
  660      rhs(
i) += val * t3(
i, 
j) * diff_base_functions(
j);
 
  662      ++diff_base_functions;
 
  664    for (; bb != nb_base_functions; bb++) {
 
  665      ++diff_base_functions;
 
  669  CHKERR aSemble(row_side, row_type, row_data);
 
 
  681      dAta(data), 
commonData(common_data), ghostVec(ghost_vec, true),
 
  682      fieldDisp(field_disp) {}
 
 
  685    int row_side, EntityType row_type,
 
  689  if (row_type != MBVERTEX)
 
  691  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
 
  696  std::vector<MatrixDouble> &
F =
 
  698  dAta.materialDoublePtr->F.resize(3, 3, 
false);
 
  702  for (
unsigned int gg = 0; gg != row_data.
getN().size1(); ++gg) {
 
  703    double val = getVolume() * getGaussPts()(3, gg);
 
  704    noalias(dAta.materialDoublePtr->F) = 
F[gg];
 
  706      for (
int dd = 0; dd < 3; dd++) {
 
  707        dAta.materialDoublePtr->F(dd, dd) += 1;
 
  710    int nb_active_variables = 0;
 
  711    CHKERR dAta.materialDoublePtr->setUserActiveVariables(nb_active_variables);
 
  712    CHKERR dAta.materialDoublePtr->calculateElasticEnergy(
 
  713        dAta, getNumeredEntFiniteElementPtr());
 
  714    energy += val * dAta.materialDoublePtr->eNergy;
 
  717  CHKERR VecSetValue(ghostVec, 0, energy, ADD_VALUES);
 
 
  726      dAta(data), 
commonData(common_data), aLe(false) {}
 
 
  739      const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
 
  743      &jac_stress(3 * 0 + 0, S + 0), &jac_stress(3 * 0 + 0, S + 1),
 
  744      &jac_stress(3 * 0 + 0, S + 2), &jac_stress(3 * 0 + 1, S + 0),
 
  745      &jac_stress(3 * 0 + 1, S + 1), &jac_stress(3 * 0 + 1, S + 2),
 
  746      &jac_stress(3 * 0 + 2, S + 0), &jac_stress(3 * 0 + 2, S + 1),
 
  747      &jac_stress(3 * 0 + 2, S + 2), &jac_stress(3 * 1 + 0, S + 0),
 
  748      &jac_stress(3 * 1 + 0, S + 1), &jac_stress(3 * 1 + 0, S + 2),
 
  749      &jac_stress(3 * 1 + 1, S + 0), &jac_stress(3 * 1 + 1, S + 1),
 
  750      &jac_stress(3 * 1 + 1, S + 2), &jac_stress(3 * 1 + 2, S + 0),
 
  751      &jac_stress(3 * 1 + 2, S + 1), &jac_stress(3 * 1 + 2, S + 2),
 
  752      &jac_stress(3 * 2 + 0, S + 0), &jac_stress(3 * 2 + 0, S + 1),
 
  753      &jac_stress(3 * 2 + 0, S + 2), &jac_stress(3 * 2 + 1, S + 0),
 
  754      &jac_stress(3 * 2 + 1, S + 1), &jac_stress(3 * 2 + 1, S + 2),
 
  755      &jac_stress(3 * 2 + 2, S + 0), &jac_stress(3 * 2 + 2, S + 1),
 
  756      &jac_stress(3 * 2 + 2, S + 2));
 
  758      &jac_stress(3 * 0 + 0, S + 3), &jac_stress(3 * 0 + 0, S + 4),
 
  759      &jac_stress(3 * 0 + 0, S + 5), &jac_stress(3 * 0 + 1, S + 3),
 
  760      &jac_stress(3 * 0 + 1, S + 4), &jac_stress(3 * 0 + 1, S + 5),
 
  761      &jac_stress(3 * 0 + 2, S + 3), &jac_stress(3 * 0 + 2, S + 4),
 
  762      &jac_stress(3 * 0 + 2, S + 5), &jac_stress(3 * 1 + 0, S + 3),
 
  763      &jac_stress(3 * 1 + 0, S + 4), &jac_stress(3 * 1 + 0, S + 5),
 
  764      &jac_stress(3 * 1 + 1, S + 3), &jac_stress(3 * 1 + 1, S + 4),
 
  765      &jac_stress(3 * 1 + 1, S + 5), &jac_stress(3 * 1 + 2, S + 3),
 
  766      &jac_stress(3 * 1 + 2, S + 4), &jac_stress(3 * 1 + 2, S + 5),
 
  767      &jac_stress(3 * 2 + 0, S + 3), &jac_stress(3 * 2 + 0, S + 4),
 
  768      &jac_stress(3 * 2 + 0, S + 5), &jac_stress(3 * 2 + 1, S + 3),
 
  769      &jac_stress(3 * 2 + 1, S + 4), &jac_stress(3 * 2 + 1, S + 5),
 
  770      &jac_stress(3 * 2 + 2, S + 3), &jac_stress(3 * 2 + 2, S + 4),
 
  771      &jac_stress(3 * 2 + 2, S + 5));
 
  773      &jac_stress(3 * 0 + 0, S + 6), &jac_stress(3 * 0 + 0, S + 7),
 
  774      &jac_stress(3 * 0 + 0, S + 8), &jac_stress(3 * 0 + 1, S + 6),
 
  775      &jac_stress(3 * 0 + 1, S + 7), &jac_stress(3 * 0 + 1, S + 8),
 
  776      &jac_stress(3 * 0 + 2, S + 6), &jac_stress(3 * 0 + 2, S + 7),
 
  777      &jac_stress(3 * 0 + 2, S + 8), &jac_stress(3 * 1 + 0, S + 6),
 
  778      &jac_stress(3 * 1 + 0, S + 7), &jac_stress(3 * 1 + 0, S + 8),
 
  779      &jac_stress(3 * 1 + 1, S + 6), &jac_stress(3 * 1 + 1, S + 7),
 
  780      &jac_stress(3 * 1 + 1, S + 8), &jac_stress(3 * 1 + 2, S + 6),
 
  781      &jac_stress(3 * 1 + 2, S + 7), &jac_stress(3 * 1 + 2, S + 8),
 
  782      &jac_stress(3 * 2 + 0, S + 6), &jac_stress(3 * 2 + 0, S + 7),
 
  783      &jac_stress(3 * 2 + 0, S + 8), &jac_stress(3 * 2 + 1, S + 6),
 
  784      &jac_stress(3 * 2 + 1, S + 7), &jac_stress(3 * 2 + 1, S + 8),
 
  785      &jac_stress(3 * 2 + 2, S + 6), &jac_stress(3 * 2 + 2, S + 7),
 
  786      &jac_stress(3 * 2 + 2, S + 8));
 
  790      &jac(0, 0), &jac(1, 0), &jac(2, 0), &jac(3, 0), &jac(4, 0), &jac(5, 0),
 
  791      &jac(6, 0), &jac(7, 0), &jac(8, 0));
 
  793      &jac(0, 1), &jac(1, 1), &jac(2, 1), &jac(3, 1), &jac(4, 1), &jac(5, 1),
 
  794      &jac(6, 1), &jac(7, 1), &jac(8, 1));
 
  796      &jac(0, 2), &jac(1, 2), &jac(2, 2), &jac(3, 2), &jac(4, 2), &jac(5, 2),
 
  797      &jac(6, 2), &jac(7, 2), &jac(8, 2));
 
  799      diff_ptr, &diff_ptr[1], &diff_ptr[2]);
 
  800  for (
int dd = 0; dd != nb_col / 3; ++dd) {
 
  801    t2_1_0(
i, 
j) += t3_1_0(
i, 
j, 
k) * diff(
k);
 
  802    t2_1_1(
i, 
j) += t3_1_1(
i, 
j, 
k) * diff(
k);
 
  803    t2_1_2(
i, 
j) += t3_1_2(
i, 
j, 
k) * diff(
k);
 
 
  818    int row_side, 
int col_side, EntityType row_type, EntityType col_type,
 
  826  int *row_indices_ptr = &row_data.
getIndices()[0];
 
  827  int *col_indices_ptr = &col_data.
getIndices()[0];
 
  829  if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
 
  830    rowIndices.resize(nb_row, 
false);
 
  832    row_indices_ptr = &rowIndices[0];
 
  834    VectorDofs::iterator dit = dofs.begin();
 
  835    for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
 
  836      if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
 
  837          dAta.forcesOnlyOnEntitiesRow.end()) {
 
  843  if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
 
  844    colIndices.resize(nb_col, 
false);
 
  846    col_indices_ptr = &colIndices[0];
 
  848    VectorDofs::iterator dit = dofs.begin();
 
  849    for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
 
  850      if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
 
  851          dAta.forcesOnlyOnEntitiesCol.end()) {
 
  858                      col_indices_ptr, &
k(0, 0), ADD_VALUES);
 
  861  if (row_side != col_side || row_type != col_type) {
 
  866    if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
 
  867      rowIndices.resize(nb_row, 
false);
 
  869      row_indices_ptr = &rowIndices[0];
 
  871      VectorDofs::iterator dit = dofs.begin();
 
  872      for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
 
  873        if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
 
  874            dAta.forcesOnlyOnEntitiesCol.end()) {
 
  880    if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
 
  881      colIndices.resize(nb_col, 
false);
 
  883      col_indices_ptr = &colIndices[0];
 
  885      VectorDofs::iterator dit = dofs.begin();
 
  886      for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
 
  887        if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
 
  888            dAta.forcesOnlyOnEntitiesRow.end()) {
 
  894    trans_k.resize(nb_col, nb_row, 
false);
 
  895    noalias(trans_k) = trans(
k);
 
  897                        row_indices_ptr, &trans_k(0, 0), ADD_VALUES);
 
 
  904    int row_side, 
int col_side, EntityType row_type, EntityType col_type,
 
  916  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
 
  922  const int nb_gauss_pts = row_data.
getN().size1();
 
  928  k.resize(nb_row, nb_col, 
false);
 
  930  jac.resize(9, nb_col, 
false);
 
  932  for (
int gg = 0; gg != nb_gauss_pts; gg++) {
 
  933    CHKERR getJac(col_data, gg);
 
  934    double val = getVolume() * getGaussPts()(3, gg);
 
  936        &jac(3 * 0 + 0, 0), &jac(3 * 0 + 0, 1), &jac(3 * 0 + 0, 2),
 
  937        &jac(3 * 0 + 1, 0), &jac(3 * 0 + 1, 1), &jac(3 * 0 + 1, 2),
 
  938        &jac(3 * 0 + 2, 0), &jac(3 * 0 + 2, 1), &jac(3 * 0 + 2, 2),
 
  939        &jac(3 * 1 + 0, 0), &jac(3 * 1 + 0, 1), &jac(3 * 1 + 0, 2),
 
  940        &jac(3 * 1 + 1, 0), &jac(3 * 1 + 1, 1), &jac(3 * 1 + 1, 2),
 
  941        &jac(3 * 1 + 2, 0), &jac(3 * 1 + 2, 1), &jac(3 * 1 + 2, 2),
 
  942        &jac(3 * 2 + 0, 0), &jac(3 * 2 + 0, 1), &jac(3 * 2 + 0, 2),
 
  943        &jac(3 * 2 + 1, 0), &jac(3 * 2 + 1, 1), &jac(3 * 2 + 1, 2),
 
  944        &jac(3 * 2 + 2, 0), &jac(3 * 2 + 2, 1), &jac(3 * 2 + 2, 2));
 
  945    for (
int cc = 0; cc != nb_col / 3; cc++) {
 
  948          &
k(0, 3 * cc + 0), &
k(0, 3 * cc + 1), &
k(0, 3 * cc + 2),
 
  949          &
k(1, 3 * cc + 0), &
k(1, 3 * cc + 1), &
k(1, 3 * cc + 2),
 
  950          &
k(2, 3 * cc + 0), &
k(2, 3 * cc + 1), &
k(2, 3 * cc + 2), 3 * nb_col);
 
  951      for (
int rr = 0; rr != nb_row / 3; rr++) {
 
  952        lhs(
i, 
j) += val * t3_1(
i, 
m, 
j) * diff_base_functions(
m);
 
  953        ++diff_base_functions;
 
  960  CHKERR aSemble(row_side, col_side, row_type, col_type, row_data, col_data);
 
 
  978    int row_side, 
int col_side, EntityType row_type, EntityType col_type,
 
  986  int *row_indices_ptr = &row_data.
getIndices()[0];
 
  987  if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
 
  988    rowIndices.resize(nb_row, 
false);
 
  990    row_indices_ptr = &rowIndices[0];
 
  992    VectorDofs::iterator dit = dofs.begin();
 
  993    for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
 
  994      if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
 
  995          dAta.forcesOnlyOnEntitiesRow.end()) {
 
 1001  int *col_indices_ptr = &col_data.
getIndices()[0];
 
 1002  if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
 
 1003    colIndices.resize(nb_col, 
false);
 
 1005    col_indices_ptr = &colIndices[0];
 
 1007    VectorDofs::iterator dit = dofs.begin();
 
 1008    for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
 
 1009      if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
 
 1010          dAta.forcesOnlyOnEntitiesCol.end()) {
 
 1011        colIndices[ii] = -1;
 
 1025                      col_indices_ptr, &
k(0, 0), ADD_VALUES);
 
 
 1032    int tag, 
bool jacobian, 
bool ale)
 
 1034                                     jacobian, ale, false) {}
 
 
 1041  CHKERR dAta.materialAdoublePtr->calculatesIGma_EshelbyStress(
 
 1042      dAta, getNumeredEntFiniteElementPtr());
 
 1044    auto &t_sIGma = dAta.materialAdoublePtr->t_sIGma;
 
 1045    auto &t_invH = dAta.materialAdoublePtr->t_invH;
 
 1046    t_sIGma(
i, 
j) = t_sIGma(
i, 
k) * t_invH(
j, 
k);
 
 1047    t_sIGma(
i, 
j) *= dAta.materialAdoublePtr->detH;
 
 1051  for (
int dd1 = 0; dd1 < 3; dd1++) {
 
 1052    for (
int dd2 = 0; dd2 < 3; dd2++) {
 
 1053      dAta.materialAdoublePtr->sIGma(dd1, dd2) >>=
 
 
 1089        materialAdoublePtr) {
 
 1092  if (!materialDoublePtr) {
 
 1094            "Pointer for materialDoublePtr not allocated");
 
 1096  if (!materialAdoublePtr) {
 
 1098            "Pointer for materialAdoublePtr not allocated");
 
 1104    CHKERR it->getAttributeDataStructure(mydata);
 
 1105    int id = it->getMeshsetId();
 
 1112    setOfBlocks[id].materialDoublePtr = materialDoublePtr;
 
 1113    setOfBlocks[id].materialAdoublePtr = materialAdoublePtr;
 
 
 1120    const std::string element_name,
 
 1121    const std::string spatial_position_field_name,
 
 1122    const std::string material_position_field_name, 
const bool ale) {
 
 1127      element_name, spatial_position_field_name);
 
 1129      element_name, spatial_position_field_name);
 
 1131      element_name, spatial_position_field_name);
 
 1135          element_name, material_position_field_name);
 
 1137          element_name, material_position_field_name);
 
 1140        element_name, material_position_field_name);
 
 1143  std::map<int, BlockData>::iterator sit = 
setOfBlocks.begin();
 
 
 1153    const std::string spatial_position_field_name,
 
 1154    const std::string material_position_field_name, 
const bool ale,
 
 1155    const bool field_disp) {
 
 1168  std::map<int, BlockData>::iterator sit = 
setOfBlocks.begin();
 
 1171        spatial_position_field_name, sit->second, 
commonData, 
tAg, 
false, ale,
 
 1174        spatial_position_field_name, sit->second, 
commonData));
 
 1201        spatial_position_field_name, sit->second, 
commonData, 
tAg, 
true, ale,
 
 1204        spatial_position_field_name, spatial_position_field_name, sit->second,
 
 
static MoFEMErrorCode calculateEnergy(DM dm, boost::shared_ptr< map< int, BlockData > > block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, SmartPetscObj< Vec > &v_energy_ptr)
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static MoFEMErrorCode get_jac(EntitiesFieldData::EntData &col_data, int gg, MatrixDouble &jac_stress, MatrixDouble &jac)
Operators and data structures for non-linear elastic analysis.
#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 ...
@ 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
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
const double n
refractive index of diffusive medium
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorShallowArrayAdaptor< double > VectorAdaptor
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
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.
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
FTensor::Index< 'm', 3 > m
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
bool sYmm
If true assume that matrix is symmetric structure.
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.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Return scalar files as a FTensor of rank 0.
const VectorDofs & getFieldDofs() const
Get DOF data structures (const version)
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Elastic material data structure.
intrusive_ptr for managing petsc objects
FTensor::Tensor2< double *, 3, 3 > & getJac()
get element Jacobian
Volume finite element base.
data for calculation heat conductivity and heat capacity elements
common data used by volume elements
std::vector< VectorDouble > jacEnergy
std::vector< MatrixDouble > jacStress
this is simply material tangent operator
std::vector< double > eNergy
std::vector< VectorDouble > hessianEnergy
std::vector< MatrixDouble3by3 > sTress
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
Implementation of elastic (non-linear) St. Kirchhoff equation.
MoFEMErrorCode preProcess()
Pre-processing function executed at loop initialization.
MoFEMErrorCode postProcess()
Post-processing function executed at loop completion.
MyVolumeFE(MoFEM::Interface &m_field)
int getRule(int order)
it is used to calculate nb. of Gauss integration points
OpEnergy(const std::string field_name, BlockData &data, CommonData &common_data, SmartPetscObj< Vec > ghost_vec, bool field_disp)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpGetCommonDataAtGaussPts(const std::string field_name, CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
operator calculating deformation gradient
OpGetDataAtGaussPts(const std::string field_name, std::vector< VectorDouble > &values_at_gauss_pts, std::vector< MatrixDouble > &gradient_at_gauss_pts)
virtual MoFEMErrorCode playTag(const int gg)
Play ADOL-C tape.
OpJacobianEnergy(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool gradient, bool hessian, bool ale, bool field_disp)
virtual MoFEMErrorCode recordTag(const int gg)
Record ADOL-C tape.
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
virtual MoFEMErrorCode calculateEnergy(const int gg)
Calculate Paola-Kirchhoff I stress.
MoFEMErrorCode calculateStress(const int gg)
Calculate Paola-Kirchhoff I stress.
OpJacobianEshelbyStress(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian, bool ale)
Operator performs automatic differentiation.
virtual MoFEMErrorCode recordTag(const int gg)
Record ADOL-C tape.
virtual MoFEMErrorCode playTag(const int gg)
Play ADOL-C tape.
OpJacobianPiolaKirchhoffStress(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian, bool ale, bool field_disp)
Construct operator to calculate Piola-Kirchhoff stress or its derivatives over gradient deformation.
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Calculate stress or jacobian at gauss points.
virtual MoFEMErrorCode calculateStress(const int gg)
Calculate Paola-Kirchhoff I stress.
OpLhsEshelby_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpLhsEshelby_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpLhsPiolaKirchhoff_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
virtual MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpLhsPiolaKirchhoff_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpRhsEshelbyStress(const std::string field_name, BlockData &data, CommonData &common_data)
virtual MoFEMErrorCode aSemble(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpRhsPiolaKirchhoff(const std::string field_name, BlockData &data, CommonData &common_data)
MyVolumeFE feRhs
calculate right hand side for tetrahedral elements
NonlinearElasticElement(MoFEM::Interface &m_field, short int tag)
MoFEMErrorCode addElement(const std::string element_name, const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false)
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
FTensor::Index< 'k', 3 > k
FTensor::Index< 'j', 3 > j
MoFEM::Interface & mField
MoFEMErrorCode setBlocks(boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double > > materialDoublePtr, boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble > > materialAdoublePtr)
MyVolumeFE feEnergy
calculate elastic energy
FTensor::Index< 'i', 3 > i
MoFEMErrorCode setOperators(const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false, const bool field_disp=false)
Set operators to calculate left hand tangent matrix and right hand residual.