21                     boost::shared_ptr<HMHHencky> hencky_ptr)
 
   26                        "Can not get data from block");
 
 
   29    MoFEMErrorCode 
doWork(
int side, EntityType type,
 
   30                          EntitiesFieldData::EntData &data) {
 
   34        if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
 
   40                                       b.bulkModulusK, b.shearModulusG);
 
 
 
   71                   boost::shared_ptr<DataAtIntegrationPts> data_ptr,
 
   72                   boost::shared_ptr<PhysicalEquations> physics_ptr) {
 
   74        data_ptr, boost::dynamic_pointer_cast<HMHHencky>(physics_ptr)));
 
 
   80                      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
 
   81                      const double alpha_u);
 
 
   96                          boost::shared_ptr<DataAtIntegrationPts> data_ptr,
 
   97                          const double alpha_u) {
 
 
  104        boost::shared_ptr<DataAtIntegrationPts> data_ptr,
 
  105        boost::shared_ptr<ExternalStrainVec> &external_strain_vec_ptr,
 
  106        std::map<std::string, boost::shared_ptr<ScalingMethod>> smv);
 
 
  117      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
 
  118      boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
 
  119      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv)
 
  122        field_name, data_ptr, external_strain_vec_ptr, smv);
 
 
  128                            boost::shared_ptr<DataAtIntegrationPts> data_ptr,
 
 
  142      std::string row_field, std::string col_field,
 
  143      boost::shared_ptr<DataAtIntegrationPts> data_ptr, 
const double alpha) {
 
 
  163                      boost::shared_ptr<double> total_energy_ptr);
 
  164    MoFEMErrorCode 
doWork(
int side, EntityType type, 
EntData &data);
 
 
  173                          boost::shared_ptr<double> total_energy_ptr) {
 
 
  179        boost::shared_ptr<DataAtIntegrationPts> data_ptr,
 
  180        boost::shared_ptr<MatrixDouble> strain_ptr,
 
  181        boost::shared_ptr<MatrixDouble> stress_ptr,
 
  182        boost::shared_ptr<HMHHencky> hencky_ptr);
 
  183    MoFEMErrorCode 
doWork(
int side, EntityType type, 
EntData &data);
 
  186    boost::shared_ptr<DataAtIntegrationPts>
 
 
  194      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
 
  195      boost::shared_ptr<PhysicalEquations> physics_ptr) {
 
  197        data_ptr, data_ptr->getLogStretchTensorAtPts(),
 
  198        data_ptr->getApproxPAtPts(),
 
  199        boost::dynamic_pointer_cast<HMHHencky>(physics_ptr));
 
 
  203      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
 
  204      boost::shared_ptr<PhysicalEquations> physics_ptr) {
 
  206        data_ptr, data_ptr->getVarLogStreachPts(), data_ptr->getVarPiolaPts(),
 
  207        boost::dynamic_pointer_cast<HMHHencky>(physics_ptr));
 
 
  210  MoFEMErrorCode 
getOptions(boost::shared_ptr<DataAtIntegrationPts> data_ptr) {
 
  212    PetscOptionsBegin(PETSC_COMM_WORLD, 
"hencky_", 
"", 
"none");
 
  214    CHKERR PetscOptionsScalar(
"-young_modulus", 
"Young modulus", 
"", 
E, &
E,
 
  216    CHKERR PetscOptionsScalar(
"-poisson_ratio", 
"poisson ratio", 
"", 
nu, &
nu,
 
  222        << 
"Hencky: E = " << 
E << 
" nu = " << 
nu;
 
 
  235            (boost::format(
"%s(.*)") % 
"MAT_ELASTIC").str()
 
 
  247    for (
auto m : meshset_vec_ptr) {
 
  249      std::vector<double> block_data;
 
  250      CHKERR m->getAttributes(block_data);
 
  251      if (block_data.size() < 2) {
 
  253                "Expected that block has atleast two attributes");
 
  255      auto get_block_ents = [&]() {
 
 
  275  MoFEMErrorCode 
getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
 
  276                            boost::shared_ptr<MatrixDouble> mat_axiator_D_ptr,
 
  277                            boost::shared_ptr<MatrixDouble> mat_deviator_D_ptr,
 
  281    auto set_material_stiffness = [&]() {
 
  287      auto t_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
 
  288          &*(mat_D_ptr->data().begin()));
 
  289      auto t_axiator_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
 
  290          &*mat_axiator_D_ptr->data().begin());
 
  291      auto t_deviator_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
 
  292          &*mat_deviator_D_ptr->data().begin());
 
  295      t_deviator_D(
i, 
j, 
k, 
l) =
 
  297      t_D(
i, 
j, 
k, 
l) = t_axiator_D(
i, 
j, 
k, 
l) + t_deviator_D(
i, 
j, 
k, 
l);
 
  304    set_material_stiffness();
 
 
  313    auto set_material_compilance = [&]() {
 
  322      auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
 
  323          &*(mat_inv_D_ptr->data().begin()));
 
  324      t_inv_D(
i, 
j, 
k, 
l) =
 
  330    set_material_compilance();
 
 
 
  354    boost::shared_ptr<DataAtIntegrationPts> data_ptr, 
const double alpha_u)
 
  357  CHK_MOAB_THROW(PetscOptionsGetBool(PETSC_NULLPTR, 
"", 
"-poly_convex",
 
  359                 "get polyconvex option failed");
 
 
  365    CHKERR integratePolyconvexHencky(data);
 
  367    CHKERR integrateHencky(data);
 
 
  379  int nb_integration_pts = data.
getN().size1();
 
  380  auto v = getVolume();
 
  381  auto t_w = getFTensor0IntegrationWeight();
 
  382  auto t_approx_P_adjoint_log_du =
 
  383      getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
 
  384  auto t_log_stretch_h1 =
 
  385      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
 
  387      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
 
  389  auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
 
  395  auto get_ftensor2 = [](
auto &
v) {
 
  397        &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
 
  400  int nb_base_functions = data.
getN().size2();
 
  402  for (
int gg = 0; gg != nb_integration_pts; ++gg) {
 
  404    auto t_nf = get_ftensor2(nF);
 
  408        t_D(
i, 
j, 
k, 
l) * (t_log_stretch_h1(
k, 
l) + alphaU * t_dot_log_u(
k, 
l));
 
  411        a * (t_approx_P_adjoint_log_du(L) - t_L(
i, 
j, L) * t_T(
i, 
j));
 
  414    for (; bb != nb_dofs / 6; ++bb) {
 
  415      t_nf(L) -= t_row_base_fun * t_residual(L);
 
  419    for (; bb != nb_base_functions; ++bb)
 
  423    ++t_approx_P_adjoint_log_du;
 
 
  438  int nb_integration_pts = data.
getN().size1();
 
  439  auto v = getVolume();
 
  440  auto t_w = getFTensor0IntegrationWeight();
 
  441  auto t_approx_P_adjoint_log_du =
 
  442      getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
 
  443  auto t_log_stretch_h1 =
 
  444      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
 
  446      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
 
  448  auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
 
  454  auto get_ftensor2 = [](
auto &
v) {
 
  456        &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
 
  459  constexpr double nohat_k = 1. / 4;
 
  460  constexpr double hat_k = 1. / 8;
 
  461  double mu = dataAtPts->mu;
 
  462  double lambda = dataAtPts->lambda;
 
  464  constexpr double third = boost::math::constants::third<double>();
 
  468  int nb_base_functions = data.
getN().size2();
 
  470  for (
int gg = 0; gg != nb_integration_pts; ++gg) {
 
  472    auto t_nf = get_ftensor2(nF);
 
  474    double log_det = t_log_stretch_h1(
i, 
i);
 
  475    double log_det2 = log_det * log_det;
 
  478    double dev_norm2 = t_dev(
i, 
j) * t_dev(
i, 
j);
 
  481    auto A = 2 * 
mu * std::exp(nohat_k * dev_norm2);
 
  482    auto B = 
lambda * std::exp(hat_k * log_det2) * log_det;
 
  485        A * (t_dev(
k, 
l) * t_diff_deviator(
k, 
l, 
i, 
j))
 
  493        alphaU * t_D(
i, 
j, 
k, 
l) * t_dot_log_u(
k, 
l);
 
  497        a * (t_approx_P_adjoint_log_du(L) - t_L(
i, 
j, L) * t_T(
i, 
j));
 
  500    for (; bb != nb_dofs / 
size_symm; ++bb) {
 
  501      t_nf(L) -= t_row_base_fun * t_residual(L);
 
  505    for (; bb != nb_base_functions; ++bb)
 
  509    ++t_approx_P_adjoint_log_du;
 
 
  517    std::string row_field, std::string col_field,
 
  518    boost::shared_ptr<DataAtIntegrationPts> data_ptr, 
const double alpha)
 
  523  CHK_MOAB_THROW(PetscOptionsGetBool(PETSC_NULLPTR, 
"", 
"-poly_convex",
 
  525                 "get polyconvex option failed");
 
 
  530    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
 
  531    boost::shared_ptr<ExternalStrainVec> &external_strain_vec_ptr,
 
  532    std::map<std::string, boost::shared_ptr<ScalingMethod>> smv)
 
  534      externalStrainVecPtr(external_strain_vec_ptr), scalingMethodsMap(smv) {}
 
 
  550  for (
auto &ext_strain_block : (*externalStrainVecPtr)) {
 
  552    if (ext_strain_block.ents.find(fe_ent) != ext_strain_block.ents.end()) {
 
  555      if (scalingMethodsMap.find(ext_strain_block.blockName) !=
 
  556          scalingMethodsMap.end()) {
 
  558            scalingMethodsMap.at(ext_strain_block.blockName)->getScale(time);
 
  561            << 
"No scaling method found for " << ext_strain_block.blockName;
 
  565      double external_strain_val = 
scale * ext_strain_block.val;
 
  570      int nb_integration_pts = data.
getN().size1();
 
  571      auto v = getVolume();
 
  572      auto t_w = getFTensor0IntegrationWeight();
 
  580      auto get_ftensor2 = [](
auto &
v) {
 
  582            &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
 
  585      int nb_base_functions = data.
getN().size2();
 
  587      for (
int gg = 0; gg != nb_integration_pts; ++gg) {
 
  588        auto tr = 3.0 * external_strain_val;
 
  590        auto t_nf = get_ftensor2(nF);
 
  597        t_residual(L) = 
a * (t_L(
i, 
j, L) * t_T(
i, 
j));
 
  600        for (; bb != nb_dofs / 6; ++bb) {
 
  601          t_nf(L) += t_row_base_fun * t_residual(L);
 
  605        for (; bb != nb_base_functions; ++bb)
 
 
  622    CHKERR integratePolyconvexHencky(row_data, col_data);
 
  624    CHKERR integrateHencky(row_data, col_data);
 
 
  639  int nb_integration_pts = row_data.
getN().size1();
 
  640  int row_nb_dofs = row_data.
getIndices().size();
 
  641  int col_nb_dofs = col_data.
getIndices().size();
 
  643  auto get_ftensor2 = [](MatrixDouble &
m, 
const int r, 
const int c) {
 
  647        &
m(r + 0, 
c + 0), &
m(r + 0, 
c + 1), &
m(r + 0, 
c + 2), &
m(r + 0, 
c + 3),
 
  648        &
m(r + 0, 
c + 4), &
m(r + 0, 
c + 5),
 
  650        &
m(r + 1, 
c + 0), &
m(r + 1, 
c + 1), &
m(r + 1, 
c + 2), &
m(r + 1, 
c + 3),
 
  651        &
m(r + 1, 
c + 4), &
m(r + 1, 
c + 5),
 
  653        &
m(r + 2, 
c + 0), &
m(r + 2, 
c + 1), &
m(r + 2, 
c + 2), &
m(r + 2, 
c + 3),
 
  654        &
m(r + 2, 
c + 4), &
m(r + 2, 
c + 5),
 
  656        &
m(r + 3, 
c + 0), &
m(r + 3, 
c + 1), &
m(r + 3, 
c + 2), &
m(r + 3, 
c + 3),
 
  657        &
m(r + 3, 
c + 4), &
m(r + 3, 
c + 5),
 
  659        &
m(r + 4, 
c + 0), &
m(r + 4, 
c + 1), &
m(r + 4, 
c + 2), &
m(r + 4, 
c + 3),
 
  660        &
m(r + 4, 
c + 4), &
m(r + 4, 
c + 5),
 
  662        &
m(r + 5, 
c + 0), &
m(r + 5, 
c + 1), &
m(r + 5, 
c + 2), &
m(r + 5, 
c + 3),
 
  663        &
m(r + 5, 
c + 4), &
m(r + 5, 
c + 5)
 
  674  auto v = getVolume();
 
  675  auto t_w = getFTensor0IntegrationWeight();
 
  677  auto t_approx_P_adjoint__dstretch =
 
  678      getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
 
  679  auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
 
  680  auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
 
  682  int row_nb_base_functions = row_data.
getN().size2();
 
  685  auto get_dP = [&]() {
 
  687    auto ts_a = getTSa();
 
  689    auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
 
  691    t_dP_tmp(L, 
J) = -(1 + alphaU * ts_a) *
 
  693                      ((t_D(
i, 
j, 
m, 
n) * t_diff(
m, 
n, 
k, 
l)) * t_L(
k, 
l, 
J)));
 
  697      auto t_approx_P_adjoint__dstretch =
 
  698          getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
 
  699      auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
 
  700      auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
 
  701      auto &nbUniq = dataAtPts->nbUniq;
 
  703      auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
 
  704      for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
 
  709        t_sym(
i, 
j) = (t_approx_P_adjoint__dstretch(
i, 
j) ||
 
  710                       t_approx_P_adjoint__dstretch(
j, 
i));
 
  715        t_dP(L, 
J) = t_L(
i, 
j, L) *
 
  716                         ((t_diff2_uP2(
i, 
j, 
k, 
l) + t_diff2_uP2(
k, 
l, 
i, 
j)) *
 
  722        ++t_approx_P_adjoint__dstretch;
 
  727      auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
 
  728      for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
 
  729        t_dP(L, 
J) = t_dP_tmp(L, 
J);
 
  734    return getFTensor2FromMat<size_symm, size_symm>(dP);
 
  737  auto t_dP = get_dP();
 
  738  for (
int gg = 0; gg != nb_integration_pts; ++gg) {
 
  742    for (; rr != row_nb_dofs / 6; ++rr) {
 
  744      auto t_m = get_ftensor2(K, 6 * rr, 0);
 
  745      for (
int cc = 0; cc != col_nb_dofs / 6; ++cc) {
 
  746        const double b = 
a * t_row_base_fun * t_col_base_fun;
 
  747        t_m(L, 
J) -= b * t_dP(L, 
J);
 
  754    for (; rr != row_nb_base_functions; ++rr) {
 
 
  773  int nb_integration_pts = row_data.
getN().size1();
 
  774  int row_nb_dofs = row_data.
getIndices().size();
 
  775  int col_nb_dofs = col_data.
getIndices().size();
 
  777  auto get_ftensor2 = [](MatrixDouble &
m, 
const int r, 
const int c) {
 
  781        &
m(r + 0, 
c + 0), &
m(r + 0, 
c + 1), &
m(r + 0, 
c + 2), &
m(r + 0, 
c + 3),
 
  782        &
m(r + 0, 
c + 4), &
m(r + 0, 
c + 5),
 
  784        &
m(r + 1, 
c + 0), &
m(r + 1, 
c + 1), &
m(r + 1, 
c + 2), &
m(r + 1, 
c + 3),
 
  785        &
m(r + 1, 
c + 4), &
m(r + 1, 
c + 5),
 
  787        &
m(r + 2, 
c + 0), &
m(r + 2, 
c + 1), &
m(r + 2, 
c + 2), &
m(r + 2, 
c + 3),
 
  788        &
m(r + 2, 
c + 4), &
m(r + 2, 
c + 5),
 
  790        &
m(r + 3, 
c + 0), &
m(r + 3, 
c + 1), &
m(r + 3, 
c + 2), &
m(r + 3, 
c + 3),
 
  791        &
m(r + 3, 
c + 4), &
m(r + 3, 
c + 5),
 
  793        &
m(r + 4, 
c + 0), &
m(r + 4, 
c + 1), &
m(r + 4, 
c + 2), &
m(r + 4, 
c + 3),
 
  794        &
m(r + 4, 
c + 4), &
m(r + 4, 
c + 5),
 
  796        &
m(r + 5, 
c + 0), &
m(r + 5, 
c + 1), &
m(r + 5, 
c + 2), &
m(r + 5, 
c + 3),
 
  797        &
m(r + 5, 
c + 4), &
m(r + 5, 
c + 5)
 
  808  auto v = getVolume();
 
  809  auto t_w = getFTensor0IntegrationWeight();
 
  811  int row_nb_base_functions = row_data.
getN().size2();
 
  814  auto get_dP = [&]() {
 
  816    auto ts_a = getTSa();
 
  818    auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
 
  820    constexpr double nohat_k = 1. / 4;
 
  821    constexpr double hat_k = 1. / 8;
 
  822    double mu = dataAtPts->mu;
 
  823    double lambda = dataAtPts->lambda;
 
  825    constexpr double third = boost::math::constants::third<double>();
 
  829    auto t_approx_P_adjoint__dstretch =
 
  830        getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
 
  831    auto t_log_stretch_h1 =
 
  832        getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
 
  833    auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
 
  834    auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
 
  835    auto &nbUniq = dataAtPts->nbUniq;
 
  837    auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
 
  838    for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
 
  840      double log_det = t_log_stretch_h1(
i, 
i);
 
  841      double log_det2 = log_det * log_det;
 
  844      double dev_norm2 = t_dev(
i, 
j) * t_dev(
i, 
j);
 
  846      auto A = 2 * 
mu * std::exp(nohat_k * dev_norm2);
 
  847      auto B = 
lambda * std::exp(hat_k * log_det2) * log_det;
 
  851          (
A * 2 * nohat_k) * (t_dev(
k, 
l) * t_diff_deviator(
k, 
l, 
i, 
j));
 
  852      t_B_diff(
i, 
j) = (
B * 2 * hat_k) * log_det * 
t_kd(
i, 
j) +
 
  856          t_A_diff(
i, 
j) * (t_dev(
m, 
n) * t_diff_deviator(
m, 
n, 
k, 
l))
 
  860          A * t_diff_deviator(
m, 
n, 
i, 
j) * t_diff_deviator(
m, 
n, 
k, 
l)
 
  866      t_dP(L, 
J) = -t_L(
i, 
j, L) *
 
  873                        (alphaU * ts_a) * (t_D(
i, 
j, 
m, 
n) * t_diff(
m, 
n, 
k, 
l)
 
  883        t_sym(
i, 
j) = (t_approx_P_adjoint__dstretch(
i, 
j) ||
 
  884                       t_approx_P_adjoint__dstretch(
j, 
i));
 
  889        t_dP(L, 
J) += t_L(
i, 
j, L) *
 
  890                      ((t_diff2_uP2(
i, 
j, 
k, 
l) + t_diff2_uP2(
k, 
l, 
i, 
j)) *
 
  896      ++t_approx_P_adjoint__dstretch;
 
  902    return getFTensor2FromMat<size_symm, size_symm>(dP);
 
  905  auto t_dP = get_dP();
 
  906  for (
int gg = 0; gg != nb_integration_pts; ++gg) {
 
  910    for (; rr != row_nb_dofs / 6; ++rr) {
 
  912      auto t_m = get_ftensor2(K, 6 * rr, 0);
 
  913      for (
int cc = 0; cc != col_nb_dofs / 6; ++cc) {
 
  914        const double b = 
a * t_row_base_fun * t_col_base_fun;
 
  915        t_m(L, 
J) -= b * t_dP(L, 
J);
 
  922    for (; rr != row_nb_base_functions; ++rr) {
 
 
  933    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
 
  934    boost::shared_ptr<double> total_energy_ptr)
 
  936      totalEnergyPtr(total_energy_ptr) {
 
  940                      "dataAtPts is not allocated. Please set it before " 
  941                      "using this operator.");
 
 
  955  int nb_integration_pts = getGaussPts().size2();
 
  957      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
 
  960  auto &mat_d = dataAtPts->matD;
 
  963             "wrong matD size, should be 6 by 6 but is %d by %d", 
size_symm,
 
  967  auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(mat_d.data().data());
 
  969  dataAtPts->energyAtPts.resize(nb_integration_pts, 
false);
 
  970  auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
 
  972  for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
 
  974    t_energy = 0.5 * (t_log_u(
i, 
j) * (t_D(
i, 
j, 
k, 
l) * t_log_u(
k, 
l)));
 
  980  if (totalEnergyPtr) {
 
  981    auto t_w = getFTensor0IntegrationWeight();
 
  982    auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
 
  983    double loc_energy = 0;
 
  984    for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
 
  985      loc_energy += t_energy * t_w;
 
  989    *totalEnergyPtr += getMeasure() * loc_energy;
 
 
  996    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
 
  997    boost::shared_ptr<MatrixDouble> strain_ptr,
 
  998    boost::shared_ptr<MatrixDouble> stress_ptr,
 
  999    boost::shared_ptr<HMHHencky> hencky_ptr)
 
 1001      strainPtr(strain_ptr), stressPtr(stress_ptr), henckyPtr(hencky_ptr) {}
 
 
 1015  auto nb_integration_pts = stressPtr->size2();
 
 1017  if (nb_integration_pts != getGaussPts().size2()) {
 
 1019            "inconsistent number of integration points");
 
 1023  auto get_D = [&]() {
 
 1025    for (
auto &b : henckyPtr->blockData) {
 
 1027      if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
 
 1028        dataAtPts->mu = b.shearModulusG;
 
 1029        dataAtPts->lambda = b.bulkModulusK - 2 * b.shearModulusG / 3;
 
 1030        CHKERR henckyPtr->getMatDPtr(
 
 1031            dataAtPts->getMatDPtr(), dataAtPts->getMatAxiatorDPtr(),
 
 1032            dataAtPts->getMatDeviatorDPtr(), b.bulkModulusK, b.shearModulusG);
 
 1033        CHKERR henckyPtr->getInvMatDPtr(dataAtPts->getMatInvDPtr(),
 
 1034                                        b.bulkModulusK, b.shearModulusG);
 
 1039    const auto E = henckyPtr->E;
 
 1040    const auto nu = henckyPtr->nu;
 
 1048    CHKERR henckyPtr->getMatDPtr(
 
 1049        dataAtPts->getMatDPtr(), dataAtPts->getMatAxiatorDPtr(),
 
 1058  strainPtr->resize(
size_symm, nb_integration_pts, 
false);
 
 1059  auto t_strain = getFTensor2SymmetricFromMat<3>(*strainPtr);
 
 1060  auto t_stress = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*stressPtr);
 
 1061  auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
 
 1062      &*dataAtPts->matInvD.data().begin());
 
 1064  auto t_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
 
 1065      &*dataAtPts->matD.data().begin());
 
 1068  const double lambda = dataAtPts->lambda;
 
 1069  const double mu = dataAtPts->mu;
 
 1074  for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
 
 1075    t_strain(
i, 
j) = t_inv_D(
i, 
j, 
k, 
l) * t_stress(
k, 
l);
 
 1079    t_stress_symm_debug(
i, 
j) = (t_stress(
i, 
j) || t_stress(
j, 
i)) / 2;
 
 1081    t_stress_symm_debug_diff(
i, 
j) =
 
 1082        t_D(
i, 
j, 
k, 
l) * t_strain(
k, 
l) - t_stress_symm_debug(
i, 
j);
 
 1084        t_stress_symm_debug_diff(
i, 
j) * t_stress_symm_debug_diff(
i, 
j);
 
 1085    double nrm0 = t_stress_symm_debug(
i, 
j) * t_stress_symm_debug(
i, 
j) +
 
 1086                  std::numeric_limits<double>::epsilon();
 
 1087    constexpr double eps = 1e-10;
 
 1088    if (std::fabs(std::sqrt(nrm / nrm0)) > 
eps) {
 
 1090          << 
"Stress symmetry check failed: " << std::endl
 
 1091          << t_stress_symm_debug_diff << std::endl
 
 1094                        "Norm is too big: " + std::to_string(nrm / nrm0));
 
 
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define FTENSOR_INDEX(DIM, I)
static PetscErrorCode ierr
Kronecker Delta class symmetric.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#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 ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
auto get_D
Create deviatoric stress tensor.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'J', DIM1 > J
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto getDiffDiffMat(A &&t_val, B &&t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, C &&t_S, const int nb)
Get the Diff Diff Mat object.
auto diff_deviator(FTensor::Ddg< double, 3, 3 > &&t_diff_stress)
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static constexpr auto size_symm
double poisson_ratio
Poisson ratio.
double young_modulus
Young modulus.
constexpr auto field_name
FTensor::Index< 'm', 3 > m
static enum StretchSelector stretchSelector
static double dynamicTime
static enum RotSelector gradApproximator
static PetscBool dynamicRelaxation
static boost::function< double(const double)> f
static boost::function< double(const double)> dd_f
static boost::function< double(const double)> d_f
Calculate energy density for Hencky material model.
boost::shared_ptr< double > totalEnergyPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpCalculateEnergy(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< double > total_energy_ptr)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< MatrixDouble > strainPtr
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
boost::shared_ptr< HMHHencky > henckyPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< MatrixDouble > stressPtr
OpCalculateStretchFromStress(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< HMHHencky > hencky_ptr)
MoFEMErrorCode evaluateRhs(EntData &data)
boost::shared_ptr< HMHHencky > henckyPtr
MoFEMErrorCode evaluateLhs(EntData &data)
OpHenckyJacobian(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< HMHHencky > hencky_ptr)
boost::shared_ptr< DataAtIntegrationPts > dataAtGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
std::map< std::string, boost::shared_ptr< ScalingMethod > > scalingMethodsMap
boost::shared_ptr< ExternalStrainVec > externalStrainVecPtr
OpSpatialPhysicalExternalStrain(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< ExternalStrainVec > &external_strain_vec_ptr, std::map< std::string, boost::shared_ptr< ScalingMethod > > smv)
MoFEMErrorCode integrate(EntData &data)
OpSpatialPhysical_du_du(std::string row_field, std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha)
MoFEMErrorCode integrateHencky(EntData &row_data, EntData &col_data)
MoFEMErrorCode integratePolyconvexHencky(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integratePolyconvexHencky(EntData &data)
OpSpatialPhysical(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha_u)
MoFEMErrorCode integrateHencky(EntData &data)
MoFEMErrorCode integrate(EntData &data)
virtual VolUserDataOperator * returnOpSpatialPhysical(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha_u)
VolUserDataOperator * returnOpCalculateStretchFromStress(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
MoFEMErrorCode extractBlockData(Sev sev)
Sev getOptionsSeverityLevels
std::vector< BlockData > blockData
MoFEMErrorCode extractBlockData(std::vector< const CubitMeshSets * > meshset_vec_ptr, Sev sev)
HMHHencky(MoFEM::Interface &m_field, const double E, const double nu)
VolUserDataOperator * returnOpCalculateVarStretchFromStress(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
VolUserDataOperator * returnOpCalculateEnergy(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< double > total_energy_ptr)
MoFEMErrorCode getMatDPtr(boost::shared_ptr< MatrixDouble > mat_D_ptr, boost::shared_ptr< MatrixDouble > mat_axiator_D_ptr, boost::shared_ptr< MatrixDouble > mat_deviator_D_ptr, double bulk_modulus_K, double shear_modulus_G)
MoFEMErrorCode getOptions(boost::shared_ptr< DataAtIntegrationPts > data_ptr)
virtual UserDataOperator * returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
VolUserDataOperator * returnOpSpatialPhysical_du_du(std::string row_field, std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha)
MoFEM::Interface & mField
virtual VolUserDataOperator * returnOpSpatialPhysicalExternalStrain(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< ExternalStrainVec > external_strain_vec_ptr, std::map< std::string, boost::shared_ptr< ScalingMethod > > smv)
MoFEMErrorCode recordTape(const int tag, DTensor2Ptr *t_h)
MoFEMErrorCode getInvMatDPtr(boost::shared_ptr< MatrixDouble > mat_inv_D_ptr, double bulk_modulus_K, double shear_modulus_G)
virtual moab::Interface & get_moab()=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::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
PetscReal ts_t
Current time value.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.