HMHHencky.cpp

Go to the documentation of this file.

/**
 * @file Hencky.cpp
 * @brief Implementation of Hencky material
 * @date 2024-08-31
 *
 * @copyright Copyright (c) 2024
 *
 */
 
 
namespace EshelbianPlasticity {
  
struct HMHHencky : public PhysicalEquations {
 
  HMHHencky(MoFEM::Interface &m_field, const double E, const double nu)
      : PhysicalEquations(0, 0), mField(m_field), E(E), nu(nu) {}
 
  MoFEMErrorCode recordTape(const int tag, DTensor2Ptr *t_h) { return 0; }
 
  struct OpHenckyJacobian : public OpJacobian {
    OpHenckyJacobian(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                     boost::shared_ptr<HMHHencky> hencky_ptr)
        : OpJacobian(), dataAtGaussPts(data_ptr), henckyPtr(hencky_ptr) {
      CHK_THROW_MESSAGE(henckyPtr->getOptions(dataAtGaussPts),
                        "getOptions failed");
      CHK_THROW_MESSAGE(henckyPtr->extractBlockData(Sev::verbose),
                        "Can not get data from block");
    }
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBegin;
 
      for (auto &b : henckyPtr->blockData) {
        if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
          dataAtGaussPts->mu = b.shearModulusG;
          dataAtGaussPts->lambda = b.bulkModulusK - 2 * b.shearModulusG / 3;
          CHKERR henckyPtr->getMatDPtr(dataAtGaussPts->getMatDPtr(),
                                       dataAtGaussPts->getMatAxiatorDPtr(),
                                       dataAtGaussPts->getMatDeviatorDPtr(),
                                       b.bulkModulusK, b.shearModulusG);
          MoFEMFunctionReturnHot(0);
        }
      }
      const auto E = henckyPtr->E;
      const auto nu = henckyPtr->nu;
 
      double bulk_modulus_K = E / (3 * (1 - 2 * nu));
      double shear_modulus_G = E / (2 * (1 + nu));
 
      dataAtGaussPts->mu = shear_modulus_G;
      dataAtGaussPts->lambda = bulk_modulus_K - 2 * shear_modulus_G / 3;
 
      CHKERR henckyPtr->getMatDPtr(dataAtGaussPts->getMatDPtr(),
                                   dataAtGaussPts->getMatAxiatorDPtr(),
                                   dataAtGaussPts->getMatDeviatorDPtr(),
                                   bulk_modulus_K, shear_modulus_G);
 
      MoFEMFunctionReturn(0);
    }
 
    MoFEMErrorCode evaluateRhs(EntData &data) { return 0; }
    MoFEMErrorCode evaluateLhs(EntData &data) { return 0; }
 
  private:
    boost::shared_ptr<DataAtIntegrationPts> dataAtGaussPts;
    boost::shared_ptr<HMHHencky> henckyPtr;
  };
 
  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) {
    return (new OpHenckyJacobian(
        data_ptr, boost::dynamic_pointer_cast<HMHHencky>(physics_ptr)));
  }
 
  struct OpSpatialPhysical : public OpAssembleVolume {
 
    OpSpatialPhysical(const std::string &field_name,
                      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                      const double alpha_u);
 
    MoFEMErrorCode integrate(EntData &data);
 
    MoFEMErrorCode integrateHencky(EntData &data);
 
    MoFEMErrorCode integratePolyconvexHencky(EntData &data);
 
  private:
    const double alphaU;
    PetscBool polyConvex = PETSC_FALSE;
  };
 
  virtual VolUserDataOperator *
  returnOpSpatialPhysical(const std::string &field_name,
                          boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                          const double alpha_u) {
    return new OpSpatialPhysical(field_name, data_ptr, alpha_u);
  }
 
  struct OpSpatialPhysicalExternalStrain : public OpAssembleVolume {
    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);
 
  private:
    boost::shared_ptr<ExternalStrainVec> externalStrainVecPtr;
    std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
  };
 
  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)
       {
    return new OpSpatialPhysicalExternalStrain(
        field_name, data_ptr, external_strain_vec_ptr, smv);
  }
 
  struct OpSpatialPhysical_du_du : public OpAssembleVolume {
    const double alphaU;
    OpSpatialPhysical_du_du(std::string row_field, std::string col_field,
                            boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                            const double alpha);
    MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
    MoFEMErrorCode integrateHencky(EntData &row_data, EntData &col_data);
    MoFEMErrorCode integratePolyconvexHencky(EntData &row_data,
                                             EntData &col_data);
 
  private:
    PetscBool polyConvex = PETSC_FALSE;
 
    MatrixDouble dP;
  };
 
  VolUserDataOperator *returnOpSpatialPhysical_du_du(
      std::string row_field, std::string col_field,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha) {
    return new OpSpatialPhysical_du_du(row_field, col_field, data_ptr, alpha);
  }
 
  /**
   * @brief Calculate energy density for Hencky material model
   *
   *
   * \f[
   *
   * \Psi(\log{\mathbf{U}}) = \frac{1}{2} U_{IJ} D_{IJKL} U_{KL} = \frac{1}{2}
   * U_{IJ} T_{IJ}
   *
   * \f]
   * where \f$T_{IJ} = D_{IJKL} U_{KL}\f$ is a a Hencky stress.
   *
   */
  struct OpCalculateEnergy : public VolUserDataOperator {
 
    OpCalculateEnergy(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                      boost::shared_ptr<double> total_energy_ptr);
    MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
  private:
    boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
    boost::shared_ptr<double> totalEnergyPtr;
  };
 
  VolUserDataOperator *
  returnOpCalculateEnergy(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                          boost::shared_ptr<double> total_energy_ptr) {
    return new OpCalculateEnergy(data_ptr, total_energy_ptr);
  }
 
  struct OpCalculateStretchFromStress : public VolUserDataOperator {
    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 doWork(int side, EntityType type, EntData &data);
 
  private:
    boost::shared_ptr<DataAtIntegrationPts>
        dataAtPts; ///< data at integration pts
    boost::shared_ptr<MatrixDouble> strainPtr;
    boost::shared_ptr<MatrixDouble> stressPtr;
    boost::shared_ptr<HMHHencky> henckyPtr;
  };
 
  VolUserDataOperator *returnOpCalculateStretchFromStress(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<PhysicalEquations> physics_ptr) {
    return new OpCalculateStretchFromStress(
        data_ptr, data_ptr->getLogStretchTensorAtPts(),
        data_ptr->getApproxPAtPts(),
        boost::dynamic_pointer_cast<HMHHencky>(physics_ptr));
  }
 
  VolUserDataOperator *returnOpCalculateVarStretchFromStress(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<PhysicalEquations> physics_ptr) {
    return new OpCalculateStretchFromStress(
        data_ptr, data_ptr->getVarLogStreachPts(), data_ptr->getVarPiolaPts(),
        boost::dynamic_pointer_cast<HMHHencky>(physics_ptr));
  }
 
  MoFEMErrorCode getOptions(boost::shared_ptr<DataAtIntegrationPts> data_ptr) {
    MoFEMFunctionBegin;
    PetscOptionsBegin(PETSC_COMM_WORLD, "hencky_", "", "none");
 
    CHKERR PetscOptionsScalar("-young_modulus", "Young modulus", "", E, &E,
                              PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-poisson_ratio", "poisson ratio", "", nu, &nu,
                              PETSC_NULLPTR);
 
    PetscOptionsEnd();
 
    MOFEM_LOG("EP", getOptionsSeverityLevels)
        << "Hencky: E = " << E << " nu = " << nu;
    getOptionsSeverityLevels = Sev::verbose;
 
    CHKERRG(ierr);
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode extractBlockData(Sev sev) {
    return extractBlockData(
 
        mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
            (boost::format("%s(.*)") % "MAT_ELASTIC").str()
 
                )),
 
        sev);
  }
 
  MoFEMErrorCode
  extractBlockData(std::vector<const CubitMeshSets *> meshset_vec_ptr,
                   Sev sev) {
    MoFEMFunctionBegin;
 
    for (auto m : meshset_vec_ptr) {
      MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
      std::vector<double> block_data;
      CHKERR m->getAttributes(block_data);
      if (block_data.size() < 2) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Expected that block has atleast two attributes");
      }
      auto get_block_ents = [&]() {
        Range ents;
        CHKERR mField.get_moab().get_entities_by_handle(m->meshset, ents, true);
        return ents;
      };
 
      double young_modulus = block_data[0];
      double poisson_ratio = block_data[1];
      double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
      double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
 
      MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
          << "E = " << young_modulus << " nu = " << poisson_ratio;
 
      blockData.push_back({bulk_modulus_K, shear_modulus_G, get_block_ents()});
    }
    MOFEM_LOG_CHANNEL("WORLD");
    MoFEMFunctionReturn(0);
  }
 
  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) {
    MoFEMFunctionBegin;
    //! [Calculate elasticity tensor]
    auto set_material_stiffness = [&]() {
      FTENSOR_INDEX(SPACE_DIM, i);
      FTENSOR_INDEX(SPACE_DIM, j);
      FTENSOR_INDEX(SPACE_DIM, k);
      FTENSOR_INDEX(SPACE_DIM, l);
      constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
      auto t_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
          &*(mat_D_ptr->data().begin()));
      auto t_axiator_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
          &*mat_axiator_D_ptr->data().begin());
      auto t_deviator_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
          &*mat_deviator_D_ptr->data().begin());
      t_axiator_D(i, j, k, l) = (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
                                t_kd(i, j) * t_kd(k, l);
      t_deviator_D(i, j, k, l) =
          2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
      t_D(i, j, k, l) = t_axiator_D(i, j, k, l) + t_deviator_D(i, j, k, l);
    };
    //! [Calculate elasticity tensor]
    constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
    mat_D_ptr->resize(size_symm, size_symm, false);
    mat_axiator_D_ptr->resize(size_symm, size_symm, false);
    mat_deviator_D_ptr->resize(size_symm, size_symm, false);
    set_material_stiffness();
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode
  getInvMatDPtr(boost::shared_ptr<MatrixDouble> mat_inv_D_ptr,
                double bulk_modulus_K, double shear_modulus_G) {
    MoFEMFunctionBegin;
    //! [Calculate elasticity tensor]
    auto set_material_compilance = [&]() {
      FTENSOR_INDEX(SPACE_DIM, i);
      FTENSOR_INDEX(SPACE_DIM, j);
      FTENSOR_INDEX(SPACE_DIM, k);
      FTENSOR_INDEX(SPACE_DIM, l);
      constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
      const double A = 1. / (2. * shear_modulus_G);
      const double B =
          (1. / (9. * bulk_modulus_K)) - (1. / (6. * shear_modulus_G));
      auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
          &*(mat_inv_D_ptr->data().begin()));
      t_inv_D(i, j, k, l) =
          A * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) + B * t_kd(i, j) * t_kd(k, l);
    };
    //! [Calculate elasticity tensor]
    constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
    mat_inv_D_ptr->resize(size_symm, size_symm, false);
    set_material_compilance();
    MoFEMFunctionReturn(0);
  }
 
  struct OpTopoSpatialPhysical : public OpAssembleVolume {
 
    OpTopoSpatialPhysical(const std::string &field_name,
                          boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                          SmartPetscObj<Vec> assemble_vec,
                          boost::shared_ptr<TopologicalData> topo_ptr,
                          const double alpha_u);
 
    MoFEMErrorCode integrate(EntData &data);
 
    MoFEMErrorCode integrateHencky(EntData &data);
 
    MoFEMErrorCode integratePolyconvexHencky(EntData &data);
 
    MoFEMErrorCode assemble(EntData &data) override;
 
  private:
    const double alphaU;
    PetscBool polyConvex = PETSC_FALSE;
    boost::shared_ptr<TopologicalData> topoDataPtr;
    SmartPetscObj<Vec> assembleVec;
  };
 
  virtual VolUserDataOperator *
  returnOpTopoSpatialPhysical(const std::string &field_name,
                              boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                              SmartPetscObj<Vec> assemble_vec,
                              boost::shared_ptr<TopologicalData> topo_ptr,
                              const double alpha_u) {
    return new OpTopoSpatialPhysical(field_name, data_ptr, assemble_vec,
                                     topo_ptr, alpha_u);
  }
 
private:
  MoFEM::Interface &mField;
 
  struct BlockData {
    double bulkModulusK;
    double shearModulusG;
    Range blockEnts;
  };
  std::vector<BlockData> blockData;
 
  double E;
  double nu;
 
  // Set verbosity level, it verbile can changes that informatno is pronated
  // only once at particular level
  Sev getOptionsSeverityLevels = Sev::inform;
};
 
HMHHencky::OpSpatialPhysical::OpSpatialPhysical(
    const std::string &field_name,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha_u)
    : OpAssembleVolume(field_name, data_ptr, OPROW), alphaU(alpha_u) {
 
  CHK_MOAB_THROW(PetscOptionsGetBool(PETSC_NULLPTR, "", "-poly_convex",
                                     &polyConvex, PETSC_NULLPTR),
                 "get polyconvex option failed");
}
 
MoFEMErrorCode HMHHencky::OpSpatialPhysical::integrate(EntData &data) {
  MoFEMFunctionBegin;
  if (polyConvex) {
    CHKERR integratePolyconvexHencky(data);
  } else {
    CHKERR integrateHencky(data);
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode HMHHencky::OpSpatialPhysical::integrateHencky(EntData &data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
  auto t_L = symm_L_tensor();
 
  int nb_dofs = data.getIndices().size();
  int nb_integration_pts = data.getN().size1();
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_approx_P_adjoint_log_du =
      getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
  auto t_log_stretch_h1 =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
  auto t_dot_log_u =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
 
  auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  auto get_ftensor2 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 6>, 6>(
        &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
  };
 
  int nb_base_functions = data.getN().size2();
  auto t_row_base_fun = data.getFTensor0N();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
    auto t_nf = get_ftensor2(nF);
 
    FTensor::Tensor2_symmetric<double, 3> t_T;
    t_T(i, j) =
        t_D(i, j, k, l) * (t_log_stretch_h1(k, l) + alphaU * t_dot_log_u(k, l));
    FTensor::Tensor1<double, size_symm> t_residual;
    t_residual(L) =
        a * (t_approx_P_adjoint_log_du(L) - t_L(i, j, L) * t_T(i, j));
 
    int bb = 0;
    for (; bb != nb_dofs / 6; ++bb) {
      t_nf(L) -= t_row_base_fun * t_residual(L);
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb)
      ++t_row_base_fun;
 
    ++t_w;
    ++t_approx_P_adjoint_log_du;
    ++t_dot_log_u;
    ++t_log_stretch_h1;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
HMHHencky::OpSpatialPhysical::integratePolyconvexHencky(EntData &data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
  auto t_L = symm_L_tensor();
 
  int nb_dofs = data.getIndices().size();
  int nb_integration_pts = data.getN().size1();
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_approx_P_adjoint_log_du =
      getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
  auto t_log_stretch_h1 =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
  auto t_dot_log_u =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
 
  auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  auto get_ftensor2 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 6>, 6>(
        &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
  };
 
  constexpr double nohat_k = 1. / 4;
  constexpr double hat_k = 1. / 8;
  double mu = dataAtPts->mu;
  double lambda = dataAtPts->lambda;
 
  constexpr double third = boost::math::constants::third<double>();
  auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
  auto t_diff_deviator = diff_deviator(diff_tensor());
 
  int nb_base_functions = data.getN().size2();
  auto t_row_base_fun = data.getFTensor0N();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
    auto t_nf = get_ftensor2(nF);
 
    double log_det = t_log_stretch_h1(i, i);
    double log_det2 = log_det * log_det;
    FTensor::Tensor2_symmetric<double, 3> t_dev;
    t_dev(i, j) = t_log_stretch_h1(i, j) - t_kd(i, j) * (third * log_det);
    double dev_norm2 = t_dev(i, j) * t_dev(i, j);
 
    FTensor::Tensor2_symmetric<double, 3> t_T;
    auto A = 2 * mu * std::exp(nohat_k * dev_norm2);
    auto B = lambda * std::exp(hat_k * log_det2) * log_det;
    t_T(i, j) =
 
        A * (t_dev(k, l) * t_diff_deviator(k, l, i, j))
 
        +
 
        B * t_kd(i, j)
 
        +
 
        alphaU * t_D(i, j, k, l) * t_dot_log_u(k, l);
 
    FTensor::Tensor1<double, size_symm> t_residual;
    t_residual(L) =
        a * (t_approx_P_adjoint_log_du(L) - t_L(i, j, L) * t_T(i, j));
 
    int bb = 0;
    for (; bb != nb_dofs / size_symm; ++bb) {
      t_nf(L) -= t_row_base_fun * t_residual(L);
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb)
      ++t_row_base_fun;
 
    ++t_w;
    ++t_approx_P_adjoint_log_du;
    ++t_dot_log_u;
    ++t_log_stretch_h1;
  }
  MoFEMFunctionReturn(0);
}
 
HMHHencky::OpSpatialPhysical_du_du::OpSpatialPhysical_du_du(
    std::string row_field, std::string col_field,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha)
    : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false),
      alphaU(alpha) {
  sYmm = false;
 
  CHK_MOAB_THROW(PetscOptionsGetBool(PETSC_NULLPTR, "", "-poly_convex",
                                     &polyConvex, PETSC_NULLPTR),
                 "get polyconvex option failed");
}
 
HMHHencky::OpSpatialPhysicalExternalStrain::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)
    : OpAssembleVolume(field_name, data_ptr, OPROW),
      externalStrainVecPtr(external_strain_vec_ptr), scalingMethodsMap(smv) {}
 
MoFEMErrorCode
HMHHencky::OpSpatialPhysicalExternalStrain::integrate(EntData &data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
 
  double time = OpAssembleVolume::getFEMethod()->ts_t;
  if (EshelbianCore::dynamicRelaxation) {
    time = EshelbianCore::dynamicTime;
  }
  // get entity of tet
  EntityHandle fe_ent = OpAssembleVolume::getFEEntityHandle();
  // iterate over all block data
 
  for (auto &ext_strain_block : (*externalStrainVecPtr)) {
    // check if finite element entity is part of the EXTERNALSTRAIN block
    if (ext_strain_block.ents.find(fe_ent) != ext_strain_block.ents.end()) {
 
      double scale = 1;
      if (scalingMethodsMap.find(ext_strain_block.blockName) !=
          scalingMethodsMap.end()) {
        scale *=
            scalingMethodsMap.at(ext_strain_block.blockName)->getScale(time);
      } else {
        MOFEM_LOG("SELF", Sev::warning)
            << "No scaling method found for " << ext_strain_block.blockName;
      }
 
      int nb_dofs = data.getIndices().size();
      int nb_integration_pts = data.getN().size1();
      auto v = getVolume();
      auto t_w = getFTensor0IntegrationWeight();
      constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
 
      double external_strain_val;
      VectorDouble v_external_strain;
      auto block_name = "(.*)ANALYTICAL_EXTERNALSTRAIN(.*)";
      std::regex reg_name(block_name);
      if (std::regex_match(ext_strain_block.blockName, reg_name)) {
        VectorDouble analytical_external_strain;
        std::string block_name_tmp;
        std::tie(block_name_tmp, v_external_strain) =
            getAnalyticalExternalStrain(this, analytical_external_strain,
                                        ext_strain_block.blockName);
      } else {
        // get ExternalStrain data from block
        external_strain_val = scale * ext_strain_block.val;
        // fill with same scalar value for all integration points
        v_external_strain.resize(nb_integration_pts);
        std::fill(v_external_strain.begin(), v_external_strain.end(),
                  external_strain_val);
      }
      auto t_external_strain = getFTensor0FromVec(v_external_strain);
      double bulk_modulus_K = ext_strain_block.bulkModulusK;
      auto t_L = symm_L_tensor();
 
      FTensor::Index<'i', 3> i;
      FTensor::Index<'j', 3> j;
      FTensor::Index<'k', 3> k;
      FTensor::Index<'l', 3> l;
      auto get_ftensor2 = [](auto &v) {
        return FTensor::Tensor1<FTensor::PackPtr<double *, 6>, 6>(
            &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
      };
 
      int nb_base_functions = data.getN().size2();
      auto t_row_base_fun = data.getFTensor0N();
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        auto tr = 3.0 * t_external_strain;
        double a = v * t_w;
        auto t_nf = get_ftensor2(nF);
 
        FTensor::Tensor2_symmetric<double, 3> t_T;
 
        t_T(i, j) = -bulk_modulus_K * tr * t_kd(i, j);
 
        FTensor::Tensor1<double, size_symm> t_residual;
        t_residual(L) = a * (t_L(i, j, L) * t_T(i, j));
 
        int bb = 0;
        for (; bb != nb_dofs / 6; ++bb) {
          t_nf(L) += t_row_base_fun * t_residual(L);
          ++t_nf;
          ++t_row_base_fun;
        }
        for (; bb != nb_base_functions; ++bb)
          ++t_row_base_fun;
        ++t_external_strain;
        ++t_w;
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
HMHHencky::OpSpatialPhysical_du_du::integrate(EntData &row_data,
                                              EntData &col_data) {
  MoFEMFunctionBegin;
 
  if (polyConvex) {
    CHKERR integratePolyconvexHencky(row_data, col_data);
  } else {
    CHKERR integrateHencky(row_data, col_data);
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
HMHHencky::OpSpatialPhysical_du_du::integrateHencky(EntData &row_data,
                                                    EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
  FTensor::Index<'J', size_symm> J;
  auto t_L = symm_L_tensor();
  auto t_diff = diff_tensor();
 
  int nb_integration_pts = row_data.getN().size1();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
 
  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor2<FTensor::PackPtr<double *, 6>, size_symm,
                            size_symm>(
 
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
        &m(r + 0, c + 4), &m(r + 0, c + 5),
 
        &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
        &m(r + 1, c + 4), &m(r + 1, c + 5),
 
        &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
        &m(r + 2, c + 4), &m(r + 2, c + 5),
 
        &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
        &m(r + 3, c + 4), &m(r + 3, c + 5),
 
        &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
        &m(r + 4, c + 4), &m(r + 4, c + 5),
 
        &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
        &m(r + 5, c + 4), &m(r + 5, c + 5)
 
    );
  };
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  FTensor::Index<'m', 3> m;
  FTensor::Index<'n', 3> n;
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
 
  auto t_approx_P_adjoint__dstretch =
      getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
  auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
  auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
 
  auto get_dP = [&]() {
    dP.resize(size_symm * size_symm, nb_integration_pts, false);
    auto ts_a = getTSa();
 
    auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
    FTensor::Tensor2<double, size_symm, size_symm> t_dP_tmp;
    t_dP_tmp(L, J) = -(1 + alphaU * ts_a) *
                     (t_L(i, j, L) *
                      ((t_D(i, j, m, n) * t_diff(m, n, k, l)) * t_L(k, l, J)));
 
    if (EshelbianCore::stretchSelector > LINEAR &&
        EshelbianCore::gradApproximator > MODERATE_ROT) {
      auto t_approx_P_adjoint__dstretch =
          getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
      auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
      auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
      auto &nbUniq = dataAtPts->nbUniq;
 
      auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
      for (auto gg = 0; gg != nb_integration_pts; ++gg) {
 
        // Work of symmetric tensor on undefined tensor is equal to the work
        // of the symmetric part of it
        FTensor::Tensor2_symmetric<double, 3> t_sym;
        t_sym(i, j) = (t_approx_P_adjoint__dstretch(i, j) ||
                       t_approx_P_adjoint__dstretch(j, i));
        t_sym(i, j) /= 2.0;
        auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
            t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
            EshelbianCore::dd_f, t_sym, nbUniq[gg]);
        t_dP(L, J) = t_L(i, j, L) *
                         ((t_diff2_uP2(i, j, k, l) + t_diff2_uP2(k, l, i, j)) *
                          t_L(k, l, J)) /
                         2. +
                     t_dP_tmp(L, J);
 
        ++t_dP;
        ++t_approx_P_adjoint__dstretch;
        ++t_eigen_vals;
        ++t_eigen_vecs;
      }
    } else {
      auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
      for (auto gg = 0; gg != nb_integration_pts; ++gg) {
        t_dP(L, J) = t_dP_tmp(L, J);
        ++t_dP;
      }
    }
 
    return getFTensor2FromMat<size_symm, size_symm>(dP);
  };
 
  auto t_dP = get_dP();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
 
    int rr = 0;
    for (; rr != row_nb_dofs / 6; ++rr) {
      auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
      auto t_m = get_ftensor2(K, 6 * rr, 0);
      for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
        const double b = a * t_row_base_fun * t_col_base_fun;
        t_m(L, J) -= b * t_dP(L, J);
        ++t_m;
        ++t_col_base_fun;
      }
      ++t_row_base_fun;
    }
 
    for (; rr != row_nb_base_functions; ++rr) {
      ++t_row_base_fun;
    }
 
    ++t_w;
    ++t_dP;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode HMHHencky::OpSpatialPhysical_du_du::integratePolyconvexHencky(
    EntData &row_data, EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
  FTensor::Index<'J', size_symm> J;
  auto t_L = symm_L_tensor();
  auto t_diff = diff_tensor();
 
  int nb_integration_pts = row_data.getN().size1();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
 
  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor2<FTensor::PackPtr<double *, 6>, size_symm,
                            size_symm>(
 
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
        &m(r + 0, c + 4), &m(r + 0, c + 5),
 
        &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
        &m(r + 1, c + 4), &m(r + 1, c + 5),
 
        &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
        &m(r + 2, c + 4), &m(r + 2, c + 5),
 
        &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
        &m(r + 3, c + 4), &m(r + 3, c + 5),
 
        &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
        &m(r + 4, c + 4), &m(r + 4, c + 5),
 
        &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
        &m(r + 5, c + 4), &m(r + 5, c + 5)
 
    );
  };
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  FTensor::Index<'m', 3> m;
  FTensor::Index<'n', 3> n;
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
 
  auto get_dP = [&]() {
    dP.resize(size_symm * size_symm, nb_integration_pts, false);
    auto ts_a = getTSa();
 
    auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
 
    constexpr double nohat_k = 1. / 4;
    constexpr double hat_k = 1. / 8;
    double mu = dataAtPts->mu;
    double lambda = dataAtPts->lambda;
 
    constexpr double third = boost::math::constants::third<double>();
    auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
    auto t_diff_deviator = diff_deviator(diff_tensor());
 
    auto t_approx_P_adjoint__dstretch =
        getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
    auto t_log_stretch_h1 =
        getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
    auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
    auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
    auto &nbUniq = dataAtPts->nbUniq;
 
    auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
 
      double log_det = t_log_stretch_h1(i, i);
      double log_det2 = log_det * log_det;
      FTensor::Tensor2_symmetric<double, 3> t_dev;
      t_dev(i, j) = t_log_stretch_h1(i, j) - t_kd(i, j) * (third * log_det);
      double dev_norm2 = t_dev(i, j) * t_dev(i, j);
 
      auto A = 2 * mu * std::exp(nohat_k * dev_norm2);
      auto B = lambda * std::exp(hat_k * log_det2) * log_det;
 
      FTensor::Tensor2_symmetric<double, 3> t_A_diff, t_B_diff;
      t_A_diff(i, j) =
          (A * 2 * nohat_k) * (t_dev(k, l) * t_diff_deviator(k, l, i, j));
      t_B_diff(i, j) = (B * 2 * hat_k) * log_det * t_kd(i, j) +
                       lambda * std::exp(hat_k * log_det2) * t_kd(i, j);
      FTensor::Ddg<double, 3, 3> t_dT;
      t_dT(i, j, k, l) =
          t_A_diff(i, j) * (t_dev(m, n) * t_diff_deviator(m, n, k, l))
 
          +
 
          A * t_diff_deviator(m, n, i, j) * t_diff_deviator(m, n, k, l)
 
          +
 
          t_B_diff(i, j) * t_kd(k, l);
 
      t_dP(L, J) = -t_L(i, j, L) *
                   ((
 
                        t_dT(i, j, k, l)
 
                        +
 
                        (alphaU * ts_a) * (t_D(i, j, m, n) * t_diff(m, n, k, l)
 
                                               )) *
                    t_L(k, l, J));
 
      // Work of symmetric tensor on undefined tensor is equal to the work
      // of the symmetric part of it
      if (EshelbianCore::stretchSelector > LINEAR &&
          EshelbianCore::gradApproximator > MODERATE_ROT) {
        FTensor::Tensor2_symmetric<double, 3> t_sym;
        t_sym(i, j) = (t_approx_P_adjoint__dstretch(i, j) ||
                       t_approx_P_adjoint__dstretch(j, i));
        t_sym(i, j) /= 2.0;
        auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
            t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
            EshelbianCore::dd_f, t_sym, nbUniq[gg]);
        t_dP(L, J) += t_L(i, j, L) *
                      ((t_diff2_uP2(i, j, k, l) + t_diff2_uP2(k, l, i, j)) *
                       t_L(k, l, J)) /
                      2.;
      }
 
      ++t_dP;
      ++t_approx_P_adjoint__dstretch;
      ++t_log_stretch_h1;
      ++t_eigen_vals;
      ++t_eigen_vecs;
    }
 
    return getFTensor2FromMat<size_symm, size_symm>(dP);
  };
 
  auto t_dP = get_dP();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
 
    int rr = 0;
    for (; rr != row_nb_dofs / 6; ++rr) {
      auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
      auto t_m = get_ftensor2(K, 6 * rr, 0);
      for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
        const double b = a * t_row_base_fun * t_col_base_fun;
        t_m(L, J) -= b * t_dP(L, J);
        ++t_m;
        ++t_col_base_fun;
      }
      ++t_row_base_fun;
    }
 
    for (; rr != row_nb_base_functions; ++rr) {
      ++t_row_base_fun;
    }
 
    ++t_w;
    ++t_dP;
  }
  MoFEMFunctionReturn(0);
}
 
HMHHencky::OpCalculateEnergy::OpCalculateEnergy(
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    boost::shared_ptr<double> total_energy_ptr)
    : VolUserDataOperator(NOSPACE, OPSPACE), dataAtPts(data_ptr),
      totalEnergyPtr(total_energy_ptr) {
 
  if (!dataAtPts) {
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "dataAtPts is not allocated. Please set it before "
                      "using this operator.");
  }
  
}
 
MoFEMErrorCode HMHHencky::OpCalculateEnergy::doWork(int side, EntityType type,
                                                    EntData &data) {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, l);
 
  int nb_integration_pts = getGaussPts().size2();
  auto t_log_u =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
 
 
  auto &mat_d = dataAtPts->matD;
  if (mat_d.size1() != size_symm || mat_d.size2() != size_symm) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
             "wrong matD size, should be 6 by 6 but is %d by %d", size_symm,
             size_symm);
  }
 
  auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(mat_d.data().data());
 
  dataAtPts->energyAtPts.resize(nb_integration_pts, false);
  auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
 
  for (auto gg = 0; gg != nb_integration_pts; ++gg) {
 
    t_energy = 0.5 * (t_log_u(i, j) * (t_D(i, j, k, l) * t_log_u(k, l)));
 
    ++t_log_u;
    ++t_energy;
  }
 
  if (totalEnergyPtr) {
    auto t_w = getFTensor0IntegrationWeight();
    auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
    double loc_energy = 0;
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
      loc_energy += t_energy * t_w;
      ++t_w;
      ++t_energy;
    }
    *totalEnergyPtr += getMeasure() * loc_energy;
  }
 
  MoFEMFunctionReturn(0);
}
 
HMHHencky::OpCalculateStretchFromStress::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)
    : VolUserDataOperator(NOSPACE, OPSPACE), dataAtPts(data_ptr),
      strainPtr(strain_ptr), stressPtr(stress_ptr), henckyPtr(hencky_ptr) {}
 
MoFEMErrorCode HMHHencky::OpCalculateStretchFromStress::doWork(int side,
                                                               EntityType type,
                                                               EntData &data) {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, l);
  FTENSOR_INDEX(SPACE_DIM, m);
  FTENSOR_INDEX(SPACE_DIM, n);
 
  auto nb_integration_pts = stressPtr->size2();
#ifndef NDEBUG
  if (nb_integration_pts != getGaussPts().size2()) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "inconsistent number of integration points");
  }
#endif // NDEBUG
 
  auto get_D = [&]() {
    MoFEMFunctionBegin;
    for (auto &b : henckyPtr->blockData) {
 
      if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
        dataAtPts->mu = b.shearModulusG;
        dataAtPts->lambda = b.bulkModulusK - 2 * b.shearModulusG / 3;
        CHKERR henckyPtr->getMatDPtr(
            dataAtPts->getMatDPtr(), dataAtPts->getMatAxiatorDPtr(),
            dataAtPts->getMatDeviatorDPtr(), b.bulkModulusK, b.shearModulusG);
        CHKERR henckyPtr->getInvMatDPtr(dataAtPts->getMatInvDPtr(),
                                        b.bulkModulusK, b.shearModulusG);
        MoFEMFunctionReturnHot(0);
      }
    }
 
    const auto E = henckyPtr->E;
    const auto nu = henckyPtr->nu;
 
    double bulk_modulus_K = E / (3 * (1 - 2 * nu));
    double shear_modulus_G = E / (2 * (1 + nu));
 
    dataAtPts->mu = shear_modulus_G;
    dataAtPts->lambda = bulk_modulus_K - 2 * shear_modulus_G / 3;
 
    CHKERR henckyPtr->getMatDPtr(
        dataAtPts->getMatDPtr(), dataAtPts->getMatAxiatorDPtr(),
        dataAtPts->getMatDeviatorDPtr(), bulk_modulus_K, shear_modulus_G);
    CHKERR henckyPtr->getInvMatDPtr(dataAtPts->getMatInvDPtr(), bulk_modulus_K,
                                    shear_modulus_G);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR get_D();
 
  strainPtr->resize(size_symm, nb_integration_pts, false);
  auto t_strain = getFTensor2SymmetricFromMat<3>(*strainPtr);
  auto t_stress = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*stressPtr);
  auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
      &*dataAtPts->matInvD.data().begin());
#ifndef NDEBUG
  auto t_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
      &*dataAtPts->matD.data().begin());
#endif
 
  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
 
  // note: add rotation, so we can extract rigid body motion, work then with
  // symmetric part.
  for (auto gg = 0; gg != nb_integration_pts; ++gg) {
    t_strain(i, j) = t_inv_D(i, j, k, l) * t_stress(k, l);
 
#ifndef NDEBUG
    FTensor::Tensor2_symmetric<double, 3> t_stress_symm_debug;
    t_stress_symm_debug(i, j) = (t_stress(i, j) || t_stress(j, i)) / 2;
    FTensor::Tensor2_symmetric<double, 3> t_stress_symm_debug_diff;
    t_stress_symm_debug_diff(i, j) =
        t_D(i, j, k, l) * t_strain(k, l) - t_stress_symm_debug(i, j);
    double nrm =
        t_stress_symm_debug_diff(i, j) * t_stress_symm_debug_diff(i, j);
    double nrm0 = t_stress_symm_debug(i, j) * t_stress_symm_debug(i, j) +
                  std::numeric_limits<double>::epsilon();
    constexpr double eps = 1e-10;
    if (std::fabs(std::sqrt(nrm / nrm0)) > eps) {
      MOFEM_LOG("SELF", Sev::error)
          << "Stress symmetry check failed: " << std::endl
          << t_stress_symm_debug_diff << std::endl
          << t_stress;
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Norm is too big: " + std::to_string(nrm / nrm0));
    }
    ++t_D;
#endif
 
    ++t_strain;
    ++t_stress;
    ++t_inv_D;
  }
 
  MoFEMFunctionReturn(0);
}
 
template <typename OP_PTR>
std::tuple<std::string, VectorDouble>
getAnalyticalExternalStrain(OP_PTR op_ptr, VectorDouble &analytical_expr,
                            const std::string block_name) {
 
  auto nb_gauss_pts = op_ptr->getGaussPts().size2();
 
  auto ts_time = op_ptr->getTStime();
  auto ts_time_step = op_ptr->getTStimeStep();
  if (EshelbianCore::dynamicRelaxation) {
    ts_time = EshelbianCore::dynamicTime;
    ts_time_step = EshelbianCore::dynamicStep;
  }
 
  MatrixDouble m_ref_coords = op_ptr->getCoordsAtGaussPts();
 
  auto v_analytical_expr = analytical_externalstrain_function(
      ts_time_step, ts_time, nb_gauss_pts, m_ref_coords, block_name);
 
#ifndef NDEBUG
  if (v_analytical_expr.size() != nb_gauss_pts)
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "Wrong number of integration pts");
#endif // NDEBUG
 
  return std::make_tuple(block_name, v_analytical_expr);
};
 
// Topo
 
HMHHencky::OpTopoSpatialPhysical::OpTopoSpatialPhysical(
    const std::string &field_name,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    SmartPetscObj<Vec> assemble_vec,
    boost::shared_ptr<TopologicalData> topo_ptr, const double alpha_u)
    : OpAssembleVolume(field_name, data_ptr, OPROW), alphaU(alpha_u),
      topoDataPtr(topo_ptr), assembleVec(assemble_vec) {
 
  CHK_MOAB_THROW(PetscOptionsGetBool(PETSC_NULLPTR, "", "-poly_convex",
                                     &polyConvex, PETSC_NULLPTR),
                 "get polyconvex option failed");
}
 
MoFEMErrorCode HMHHencky::OpTopoSpatialPhysical::integrate(EntData &data) {
  MoFEMFunctionBegin;
  if (polyConvex) {
    CHKERR integratePolyconvexHencky(data);
  } else {
    CHKERR integrateHencky(data);
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode HMHHencky::OpTopoSpatialPhysical::integrateHencky(EntData &data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
  auto t_L = symm_L_tensor();
 
  int nb_dofs = data.getIndices().size();
  int nb_integration_pts = data.getN().size1();
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_approx_P_adjoint_log_du =
      getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
  auto t_log_stretch_h1 =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
  auto t_dot_log_u =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
  auto t_var_log_u = getFTensor1FromMat<SPACE_DIM>(dataAtPts->varLogStreach);
 
  auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
 
  auto t_det = getFTensor0FromVec(topoDataPtr->detJacobianAtPts);
  auto t_inv_jac =
      getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(topoDataPtr->invJacobianAtPtr);
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  auto get_ftensor2 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 6>, 6>(
        &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
  };
 
  int nb_base_functions = data.getN().size2();
  auto t_diff_base = data.getFTensor1DiffN<3>();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
    // Calculate the variation of the gradient due to geometry change
    FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_cof;
    t_cof(i, j) = t_det * t_inv_jac(j, i);
 
 
    double a = v * t_w;
    auto t_nf = get_ftensor2(nF);
 
    FTensor::Tensor2_symmetric<double, 3> t_T;
    t_T(i, j) =
        t_D(i, j, k, l) * (t_log_stretch_h1(k, l) + alphaU * t_dot_log_u(k, l));
    FTensor::Tensor1<double, size_symm> t_residual;
    t_residual(L) =
        a * (t_approx_P_adjoint_log_du(L) - t_L(i, j, L) * t_T(i, j));
 
    int bb = 0;
    for (; bb != nb_dofs / 6; ++bb) {
      t_nf(i) -=
          t_var_log_u(L) * t_residual(L) * t_cof(i, j) * t_diff_base(j);
      ++t_nf;
      ++t_diff_base;
    }
    for (; bb != nb_base_functions; ++bb)
      ++t_diff_base;
 
    ++t_w;
    ++t_approx_P_adjoint_log_du;
    ++t_dot_log_u;
    ++t_log_stretch_h1;
    ++t_var_log_u;
 
    ++t_det;
    ++t_inv_jac;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
HMHHencky::OpTopoSpatialPhysical::integratePolyconvexHencky(EntData &data) {
  MoFEMFunctionBegin;
  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
          "Polyconvex Hencky is not implemented yet");
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode HMHHencky::OpTopoSpatialPhysical::assemble(EntData &data) {
  MoFEMFunctionBegin;
  double *vec_ptr = nF.data().data();
  const int nb_dofs = data.getIndices().size();
  int *ind_ptr = data.getIndices().data().data();
  CHKERR VecSetValues(assembleVec, nb_dofs, ind_ptr, vec_ptr, ADD_VALUES);
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity