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 EffectiveElasticParams {
  double youngModulus;
  double poissonRatio;
  double bulkModulusK;
  double shearModulusG;
  double lambda;
};
 
static auto calc_effective_elastic_params(double E, double nu,
                                         double diagonal_strain) {
  const double c10 = E / (4.0 * (1.0 + nu));
  // The coefficient of 0.5 * log(J)^2 is the Lame parameter lambda.
  const double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
  const double shear_modulus_G =
      2.0 * c10 * std::exp(-2.0 * diagonal_strain / 3.0);
  const double bulk_modulus_K = lambda + 2.0 * shear_modulus_G / 3.0;
  const double young_modulus =
      shear_modulus_G * (3.0 * lambda + 2.0 * shear_modulus_G) /
      (lambda + shear_modulus_G);
  const double poisson_ratio = lambda / (2.0 * (lambda + shear_modulus_G));
 
  return EffectiveElasticParams{young_modulus, poisson_ratio, bulk_modulus_K,
                                shear_modulus_G, lambda};
}
 
struct HMHHencky : public PhysicalEquations {
 
  HMHHencky(MoFEM::Interface &m_field, const double E, const double nu,
            const Features features)
      : PhysicalEquations(Hencky, features), mField(m_field), E(E), nu(nu) {
    CHK_THROW_MESSAGE(getOptions(), "getOptions failed");
    CHK_THROW_MESSAGE(extractBlockData(Sev::verbose),
                      "Can not get data from block");
    for (const auto &block : blockData) {
      if (block.matType != HenckyMatType::HOMOGENEOUS) {
        materialFeatures.set(NON_HOMOGENEOUS_MATERIAL);
        MOFEM_LOG("WORLD", Sev::verbose)
            << "Found non-homogeneous material block: " << block.blockName;
        break;
      }
    }
  }
 
  static constexpr int StrideMatD =
      DataAtIntegrationPts::SizeSymm * DataAtIntegrationPts::SizeSymm;
 
  template <int STRIDEMATD = 0> struct OpHenckyJacobian : public OpJacobian {
    OpHenckyJacobian(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                     boost::shared_ptr<HMHHencky> hencky_ptr)
        : OpJacobian(H1, OPLAST), dataAtGaussPts(data_ptr),
          henckyPtr(hencky_ptr) {
      std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
      doEntities[MBVERTEX] = true;
    }
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBegin;
      CHKERR henckyPtr->computeMaterialParamsAtPts<STRIDEMATD>(this, data,
                                                               dataAtGaussPts);
      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;
  };
 
  UserDataOperator *
  returnOpJacobian(const bool eval_rhs, const bool eval_lhs,
                   boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                   boost::shared_ptr<PhysicalEquations> physics_ptr) override {
 
    auto hencky_ptr = boost::dynamic_pointer_cast<HMHHencky>(physics_ptr);
 
    if (getFeatures().test(NON_HOMOGENEOUS_MATERIAL))
      return new OpHenckyJacobian<StrideMatD>(data_ptr, hencky_ptr);
    return new OpHenckyJacobian<0>(data_ptr, hencky_ptr);
  }
 
  template <int STRIDEMATD = 0>
  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);
 
  private:
    const double alphaU;
  };
 
  VolUserDataOperator *
  returnOpSpatialPhysical(const std::string &field_name,
                          boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                          const double alpha_u) override {
    if (getFeatures().test(NON_HOMOGENEOUS_MATERIAL)) {
      return new OpSpatialPhysical<StrideMatD>(field_name, data_ptr, alpha_u);
    } else {
      return new OpSpatialPhysical<0>(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;
  };
 
  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) override {
    return new OpSpatialPhysicalExternalStrain(field_name, data_ptr,
                                               external_strain_vec_ptr, smv);
  }
 
  template <int STRIDEMATD = 0>
  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);
 
  private:
    MatrixDouble dP;
  };
 
  VolUserDataOperator *returnOpSpatialPhysical_du_du(
      std::string row_field, std::string col_field,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      const double alpha) override {
 
    if (getFeatures().test(NON_HOMOGENEOUS_MATERIAL)) {
      return new OpSpatialPhysical_du_du<StrideMatD>(row_field, col_field,
                                                     data_ptr, alpha);
    } else {
      return new OpSpatialPhysical_du_du<0>(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.
   *
   */
  template <int STRIDEMATD>
  struct OpCalculateHelmholtzFreeEnergy : public VolUserDataOperator {
 
    OpCalculateHelmholtzFreeEnergy(
        boost::shared_ptr<DataAtIntegrationPts> data_ptr,
        boost::shared_ptr<double> total_helmholtz_free_energy_ptr);
    MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
  private:
    boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
    boost::shared_ptr<double> totalHelmholtzFreeEnergyPtr;
  };
 
  VolUserDataOperator *
  returnOpCalculateHelmholtzFreeEnergy(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<double> total_helmholtz_free_energy_ptr) override {
 
    if (getFeatures().test(NON_HOMOGENEOUS_MATERIAL)) {
      return new OpCalculateHelmholtzFreeEnergy<StrideMatD>(
          data_ptr, total_helmholtz_free_energy_ptr);
    } else {
      return new OpCalculateHelmholtzFreeEnergy<0>(
          data_ptr, total_helmholtz_free_energy_ptr);
    }
  }
 
  bool providesHelmholtzFreeEnergy() const override { return true; }
 
  template <int STRIDEMATD = 0>
  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,
      boost::shared_ptr<MatrixDouble> strain_ptr) override {
    return returnOpCalculateStretchFromStress(
        std::move(data_ptr), std::move(physics_ptr), std::move(strain_ptr),
        nullptr, nullptr);
  }
 
  VolUserDataOperator *returnOpCalculateStretchFromStress(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<PhysicalEquations> physics_ptr,
      boost::shared_ptr<MatrixDouble> strain_ptr,
      boost::shared_ptr<MatrixDouble> stress_ptr,
      VectorPtr external_pressure_ptr) override {
    auto hencky_ptr = boost::dynamic_pointer_cast<HMHHencky>(physics_ptr);
 
    if (getFeatures().test(NON_HOMOGENEOUS_MATERIAL)) {
      return new OpCalculateStretchFromStress<StrideMatD>(
          data_ptr,
          strain_ptr ? strain_ptr : data_ptr->getLogStretchTensorAtPts(),
          stress_ptr ? stress_ptr : data_ptr->getApproxPAtPts(), hencky_ptr);
    } else {
      return new OpCalculateStretchFromStress<0>(
          data_ptr,
          strain_ptr ? strain_ptr : data_ptr->getLogStretchTensorAtPts(),
          stress_ptr ? stress_ptr : data_ptr->getApproxPAtPts(), hencky_ptr);
    }
  }
 
  VolUserDataOperator *returnOpCalculateVarStretchFromStress(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<PhysicalEquations> physics_ptr) override {
    auto hencky_ptr = boost::dynamic_pointer_cast<HMHHencky>(physics_ptr);
 
    if (getFeatures().test(NON_HOMOGENEOUS_MATERIAL)) {
      return new OpCalculateStretchFromStress<StrideMatD>(
          data_ptr, data_ptr->getVarLogStreachPts(), data_ptr->getVarPiolaPts(),
          hencky_ptr);
    } else {
      return new OpCalculateStretchFromStress<0>(
          data_ptr, data_ptr->getVarLogStreachPts(), data_ptr->getVarPiolaPts(),
          hencky_ptr);
    }
  }
 
  MoFEMErrorCode getOptions() {
    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);
    CHKERR PetscOptionsBool("-effective_neohookean_stiffness",
                            "Use effective Neo-Hookean stiffness", "",
                            effectiveNehookeanStiffness,
                            &effectiveNehookeanStiffness, PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-effective_diagonal_strain",
                              "Diagonal logarithmic strain for effective "
                              "Neo-Hookean stiffness",
                              "", effectiveDiagonalStrain,
                              &effectiveDiagonalStrain, PETSC_NULLPTR);
 
    PetscOptionsEnd();
 
    MOFEM_LOG("EP", Sev::inform)
        << "Hencky: E = " << E << " nu = " << nu
        << " effective_neohookean_stiffness = "
        << (effectiveNehookeanStiffness ? "true" : "false")
        << " effective_diagonal_strain = " << effectiveDiagonalStrain;
 
    CHKERRG(ierr);
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode extractBlockData(Sev sev) {
    return extractBlockData(
 
        mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
            (boost::format("(.*)%s(.*)") % "_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::string block_name = m->getName();
 
      auto block_name_heterogeneous = "(.*)HETEROGENEOUS_ELASTIC(.*)";
      auto block_name_analytical = "(.*)ANALYTICAL_ELASTIC(.*)";
      std::regex reg_name_heterogeneous(block_name_heterogeneous);
      std::regex reg_name_analytical(block_name_analytical);
      const bool is_heterogeneous =
          std::regex_match(block_name, reg_name_heterogeneous);
      const bool is_analytical =
          std::regex_match(block_name, reg_name_analytical);
 
      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];
 
      HenckyMatType mat_type = HenckyMatType::HOMOGENEOUS;
      if (is_heterogeneous) {
        mat_type = HenckyMatType::HETEROGENEOUS;
      } else if (is_analytical) {
        mat_type = HenckyMatType::ANALYTICAL;
      }
 
      if (effectiveNehookeanStiffness &&
          mat_type == HenckyMatType::HOMOGENEOUS) {
        const auto effective_params = calc_effective_elastic_params(
            young_modulus, poisson_ratio, effectiveDiagonalStrain);
        young_modulus = effective_params.youngModulus;
        poisson_ratio = effective_params.poissonRatio;
      }
 
      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({block_name, young_modulus, poisson_ratio,
                           bulk_modulus_K, shear_modulus_G, get_block_ents(),
                           mat_type});
    }
    MOFEM_LOG_CHANNEL("WORLD");
    MoFEMFunctionReturn(0);
  }
 
  template <int STRIDEMATD, typename OP_PTR>
  MoFEMErrorCode computeMaterialParamsAtPts(
      OP_PTR op_ptr, EntitiesFieldData::EntData &data,
      boost::shared_ptr<DataAtIntegrationPts> dataAtGaussPts) {
    MoFEMFunctionBegin;
 
    auto getMaterialParams = [&](double E, double nu) {
      if (effectiveNehookeanStiffness) {
        return calc_effective_elastic_params(E, nu, effectiveDiagonalStrain);
      }
 
      const double bulk_modulus_K = E / (3 * (1 - 2 * nu));
      const double shear_modulus_G = E / (2 * (1 + nu));
      const double lambda = bulk_modulus_K - 2 * shear_modulus_G / 3;
      return EffectiveElasticParams{E, nu, bulk_modulus_K, shear_modulus_G,
                                    lambda};
    };
 
    auto fe_ent = op_ptr->getNumeredEntFiniteElementPtr()->getEnt();
    int nb_integration_pts = op_ptr->getGaussPts().size2();
 
    dataAtGaussPts->muAtPts.resize(nb_integration_pts, false);
    dataAtGaussPts->lambdaAtPts.resize(nb_integration_pts, false);
    dataAtGaussPts->muAtPts.clear();
    dataAtGaussPts->lambdaAtPts.clear();
 
    dataAtGaussPts->youngModulusAtPts.resize(nb_integration_pts, false);
    dataAtGaussPts->youngModulusAtPts.clear();
 
    auto t_young_modulus =
        getFTensor0FromVec(dataAtGaussPts->youngModulusAtPts);
    auto t_mu = getFTensor0FromVec(dataAtGaussPts->muAtPts);
    auto t_lambda = getFTensor0FromVec(dataAtGaussPts->lambdaAtPts);
 
    MatrixSizeHelper<
        GetFTensor4DdgFromMatType<SPACE_DIM, SPACE_DIM, STRIDEMATD, DL>,
        DL>::size(dataAtGaussPts->matD, nb_integration_pts);
    MatrixSizeHelper<
        GetFTensor4DdgFromMatType<SPACE_DIM, SPACE_DIM, STRIDEMATD, DL>,
        DL>::size(dataAtGaussPts->matAxiatorD, nb_integration_pts);
    MatrixSizeHelper<
        GetFTensor4DdgFromMatType<SPACE_DIM, SPACE_DIM, STRIDEMATD, DL>,
        DL>::size(dataAtGaussPts->matDeviatorD, nb_integration_pts);
    MatrixSizeHelper<
        GetFTensor4DdgFromMatType<SPACE_DIM, SPACE_DIM, STRIDEMATD, DL>,
        DL>::size(dataAtGaussPts->matInvD, nb_integration_pts);
 
    dataAtGaussPts->matD.clear();
    dataAtGaussPts->matAxiatorD.clear();
    dataAtGaussPts->matDeviatorD.clear();
    dataAtGaussPts->matInvD.clear();
 
    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 = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, STRIDEMATD>(
        dataAtGaussPts->matD);
    auto t_axiator_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, STRIDEMATD>(
        dataAtGaussPts->matAxiatorD);
    auto t_deviator_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, STRIDEMATD>(
        dataAtGaussPts->matDeviatorD);
    auto t_inv_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, STRIDEMATD>(
        dataAtGaussPts->matInvD);
 
    auto next = [&]() {
      ++t_young_modulus;
      ++t_mu;
      ++t_lambda;
      ++t_D;
      ++t_axiator_D;
      ++t_deviator_D;
      ++t_inv_D;
    };
 
    auto evalMatD = [&](double bulk_modulus_K, double shear_modulus_G) {
      MoFEMFunctionBegin;
      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);
      MoFEMFunctionReturn(0);
    };
 
    auto evalInvMatDPtr = [&](double bulk_modulus_K, double shear_modulus_G) {
      MoFEMFunctionBegin;
      const double A = 1. / (2. * shear_modulus_G);
      const double B =
          (1. / (9. * bulk_modulus_K)) - (1. / (6. * shear_modulus_G));
      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);
      MoFEMFunctionReturn(0);
    };
 
    // from block data (MAT_ELASTIC) or (ANALYTICAL_ELASTIC) if provided,
    // otherwise from command line options
    for (auto &b : this->blockData) {
      if (b.blockEnts.find(op_ptr->getFEEntityHandle()) != b.blockEnts.end()) {
 
        if (b.matType == HMHHencky::HenckyMatType::ANALYTICAL) {
          VectorDouble analytical_elastic;
          analytical_elastic = getAnalyticalElastic(op_ptr, b.blockName);
 
          auto t_analytical_elastic = getFTensor0FromVec(analytical_elastic);
 
          for (int gg = 0; gg != nb_integration_pts; ++gg) {
            const auto material_params =
                getMaterialParams(t_analytical_elastic, b.poissonRatio);
            t_young_modulus = material_params.youngModulus;
            t_mu = material_params.shearModulusG;
            t_lambda = material_params.lambda;
 
            CHKERR evalMatD(material_params.bulkModulusK,
                            material_params.shearModulusG);
            CHKERR evalInvMatDPtr(material_params.bulkModulusK,
                                  material_params.shearModulusG);
            ++t_analytical_elastic;
            next();
          }
 
        } else if (b.matType == HMHHencky::HenckyMatType::HETEROGENEOUS) {
          Tag tag_heterogenous_mat;
          CHKERR op_ptr->getPtrFE()->mField.get_moab().tag_get_handle(
              EshelbianCore::heterogeneousYoungModTagName.c_str(),
              tag_heterogenous_mat);
          int tag_length;
          CHKERR op_ptr->getPtrFE()->mField.get_moab().tag_get_length(
              tag_heterogenous_mat, tag_length);
          if (tag_length != 1) {
            SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                    "heterogeneous Young's modulus tag should be 1 but is %d",
                    tag_length);
          }
          if (EshelbianCore::meshTransferInterpOrder == 0) {
            // Constant interpolation (element-wise)
            double elem_young_mod = 0.0;
            CHKERR op_ptr->getPtrFE()->mField.get_moab().tag_get_data(
                tag_heterogenous_mat, &fe_ent, 1, &elem_young_mod);
 
            for (int gg = 0; gg != nb_integration_pts; ++gg) {
              const auto material_params =
                  getMaterialParams(elem_young_mod, b.poissonRatio);
              t_young_modulus = material_params.youngModulus;
              t_mu = material_params.shearModulusG;
              t_lambda = material_params.lambda;
 
              CHKERR evalMatD(material_params.bulkModulusK,
                              material_params.shearModulusG);
              CHKERR evalInvMatDPtr(material_params.bulkModulusK,
                                    material_params.shearModulusG);
              next();
            }
          } else if (EshelbianCore::meshTransferInterpOrder == 1) {
            // Linear interpolation (vertex-based)
            const EntityHandle *vert_conn;
            int vert_num;
            CHKERR op_ptr->getPtrFE()->mField.get_moab().get_connectivity(
                fe_ent, vert_conn, vert_num, true);
 
            VectorDouble vert_young_mod(vert_num);
            CHKERR op_ptr->getPtrFE()->mField.get_moab().tag_get_data(
                tag_heterogenous_mat, vert_conn, vert_num, &vert_young_mod[0]);
 
            auto t_shape_n = data.getFTensor0N();
            int nb_shape_fn = data.getN(NOBASE).size2();
 
            for (int gg = 0; gg != nb_integration_pts; ++gg) {
              t_young_modulus = 0; // Initialize to zero before accumulation
              auto t_vert_young_mod = getFTensor0FromVec(vert_young_mod);
              for (int bb = 0; bb != nb_shape_fn; ++bb) {
                t_young_modulus += t_vert_young_mod * t_shape_n;
                ++t_vert_young_mod;
                ++t_shape_n;
              }
              const auto material_params =
                  getMaterialParams(t_young_modulus, b.poissonRatio);
              t_young_modulus = material_params.youngModulus;
              t_mu = material_params.shearModulusG;
              t_lambda = material_params.lambda;
 
              CHKERR evalMatD(material_params.bulkModulusK,
                              material_params.shearModulusG);
              CHKERR evalInvMatDPtr(material_params.bulkModulusK,
                                    material_params.shearModulusG);
              next();
            }
          } else {
            SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
                    "Unsupported heterogeneous Young's modulus interpolation "
                    "order %d",
                    EshelbianCore::meshTransferInterpOrder);
          }
        } else {
          // MAT_ELASTIC block with homogeneous material properties
          for (int gg = 0; gg != nb_integration_pts; ++gg) {
            t_young_modulus = b.youngModulus;
            t_mu = b.shearModulusG;
            t_lambda = b.bulkModulusK - 2 * b.shearModulusG / 3;
 
            CHKERR evalMatD(b.bulkModulusK, b.shearModulusG);
            CHKERR evalInvMatDPtr(b.bulkModulusK, b.shearModulusG);
            next();
          }
        }
        MoFEMFunctionReturnHot(0);
      }
    }
 
    // From command line options if no block data is provided
    const auto material_params = getMaterialParams(this->E, this->nu);
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      t_young_modulus = material_params.youngModulus;
      t_mu = material_params.shearModulusG;
      t_lambda = material_params.lambda;
      CHKERR evalMatD(material_params.bulkModulusK,
                      material_params.shearModulusG);
      CHKERR evalInvMatDPtr(material_params.bulkModulusK,
                            material_params.shearModulusG);
      next();
    }
 
    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,
                          boost::shared_ptr<double> J_ptr);
 
    MoFEMErrorCode integrate(EntData &data) override;
 
    MoFEMErrorCode integrateHencky(EntData &data);
 
    MoFEMErrorCode assemble(int row_side, EntityType row_type,
                            EntData &data) override;
 
  private:
    const double alphaU;
    boost::shared_ptr<TopologicalData> topoDataPtr;
    SmartPetscObj<Vec> assembleVec;
    boost::shared_ptr<double> JPtr;
    double locJ;
  };
 
  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,
                              boost::shared_ptr<double> J_ptr) override {
    return new OpTopoSpatialPhysical(field_name, data_ptr, assemble_vec,
                                     topo_ptr, alpha_u, J_ptr);
  }
 
private:
  MoFEM::Interface &mField;
 
  enum HenckyMatType { HOMOGENEOUS, HETEROGENEOUS, ANALYTICAL };
  struct BlockData {
    std::string blockName;
    double youngModulus;
    double poissonRatio;
    double bulkModulusK;
    double shearModulusG;
    Range blockEnts;
    HenckyMatType matType;
  };
  std::vector<BlockData> blockData;
 
  double E;
  double nu;
  PetscBool effectiveNehookeanStiffness = PETSC_FALSE;
  double effectiveDiagonalStrain = 0;
 
  // Set verbosity level, it verbile can changes that informatno is pronated
  // only once at particular level
};
 
template <int STRIDEMATD>
HMHHencky::OpSpatialPhysical<STRIDEMATD>::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) {}
 
template <int STRIDEMATD>
MoFEMErrorCode
HMHHencky::OpSpatialPhysical<STRIDEMATD>::integrate(EntData &data) {
  MoFEMFunctionBegin;
  CHKERR integrateHencky(data);
  MoFEMFunctionReturn(0);
}
 
template <int STRIDEMATD>
MoFEMErrorCode
HMHHencky::OpSpatialPhysical<STRIDEMATD>::integrateHencky(EntData &data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
  auto t_L = FTensor::SymmLTensor<double, 3>();
 
  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 =
      dataAtPts->getFTensorAdjointPdU(nb_integration_pts);
  auto t_log_stretch_h1 =
      dataAtPts->getFTensorLogStretchTotal(nb_integration_pts);
  auto t_dot_log_u = dataAtPts->getFTensorLogStretchDot(nb_integration_pts);
  auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
 
  auto t_D = getFTensor4DdgFromMat<3, 3, STRIDEMATD>(dataAtPts->matD);
 
  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) {
    const double det_plasticF = determinantTensor3by3(t_plasticF);
    const double a = v * t_w * det_plasticF;
    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_D;
    ++t_w;
    ++t_approx_P_adjoint_log_du;
    ++t_dot_log_u;
    ++t_log_stretch_h1;
    ++t_plasticF;
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int STRIDEMATD>
HMHHencky::OpSpatialPhysical_du_du<STRIDEMATD>::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;
}
 
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::physicalTimeFlg) {
    time = EshelbianCore::currentPhysicalTime;
  }
  // 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();
      auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
      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 = FTensor::SymmLTensor<double, 3>();
 
      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;
        const double det_plasticF = determinantTensor3by3(t_plasticF);
        const double a = v * t_w * det_plasticF;
        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;
        ++t_plasticF;
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int STRIDEMATD>
MoFEMErrorCode
HMHHencky::OpSpatialPhysical_du_du<STRIDEMATD>::integrate(EntData &row_data,
                                                          EntData &col_data) {
  MoFEMFunctionBegin;
  CHKERR integrateHencky(row_data, col_data);
  MoFEMFunctionReturn(0);
}
 
template <int STRIDEMATD>
MoFEMErrorCode HMHHencky::OpSpatialPhysical_du_du<STRIDEMATD>::integrateHencky(
    EntData &row_data, EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
  FTensor::Index<'J', size_symm> J;
  auto t_L = FTensor::SymmLTensor<double, 3>();
  auto t_diff = FTensor::DiffTensor<double>();
 
  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_INDEXES(SPACE_DIM, i, j, k, l, m, 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 = [&]() {
    auto get_stress =
        MatrixSizeHelper<GetFTensor2FromMatType<size_symm, size_symm, -1, DL>,
                         DL>::size(dP, nb_integration_pts);
    auto ts_a = getTSa();
 
    auto t_D = getFTensor4DdgFromMat<3, 3, STRIDEMATD>(dataAtPts->matD);
    FTensor::Tensor2<double, size_symm, size_symm> t_dP_tmp;
    if constexpr (!STRIDEMATD) {
      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)));
    }
    // allocate FTensors
    FTensor::Tensor3<double, 3, 3, size_symm> L_left;
    L_left(i, j, L) = t_L(i, j, L);
    FTensor::Tensor3<double, 3, 3, size_symm> L_right;
    L_right(k, l, J) = t_L(k, l, J);
    FTensor::Tensor4<double, 3, 3, 3, 3> temp;
 
    if (EshelbianCore::stretchSelector > LINEAR &&
        EshelbianCore::gradApproximator > MODERATE_ROT) {
      auto t_approx_P_adjoint__dstretch =
          dataAtPts->getFTensorAdjointPdstretch(nb_integration_pts);
      auto t_eigen_vals = dataAtPts->getFTensorEigenVals(nb_integration_pts);
      auto t_eigen_vecs = dataAtPts->getFTensorEigenVecs(nb_integration_pts);
      auto &nbUniq = dataAtPts->nbUniq;
 
      auto t_dP = get_stress();
      for (auto gg = 0; gg != nb_integration_pts; ++gg) {
        if constexpr (STRIDEMATD) {
          temp(i, j, k, l) = t_D(i, j, m, n) * t_diff(m, n, k, l);
 
          t_dP_tmp(L, J) =
              -(1 + alphaU * ts_a) *
              (L_left(i, j, L) * (temp(i, j, k, l) * L_right(k, l, J)));
 
          ++t_D;
        }
 
        // 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 = get_stress();
      for (auto gg = 0; gg != nb_integration_pts; ++gg) {
        if constexpr (STRIDEMATD) {
          temp(i, j, k, l) = t_D(i, j, m, n) * t_diff(m, n, k, l);
 
          t_dP_tmp(L, J) =
              -(1 + alphaU * ts_a) *
              (L_left(i, j, L) * (temp(i, j, k, l) * L_right(k, l, J)));
 
          ++t_D;
        }
        t_dP(L, J) = t_dP_tmp(L, J);
 
        ++t_dP;
      }
    }
 
    return get_stress();
  };
 
  auto t_dP = get_dP();
  auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    const double det_plasticF = determinantTensor3by3(t_plasticF);
    const double a = v * t_w * det_plasticF;
 
    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;
    ++t_plasticF;
  }
  MoFEMFunctionReturn(0);
}
 
template <int STRIDEMATD>
HMHHencky::OpCalculateHelmholtzFreeEnergy<
    STRIDEMATD>::OpCalculateHelmholtzFreeEnergy(
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    boost::shared_ptr<double> total_helmholtz_free_energy_ptr)
    : VolUserDataOperator(NOSPACE, OPSPACE), dataAtPts(data_ptr),
      totalHelmholtzFreeEnergyPtr(total_helmholtz_free_energy_ptr) {
 
  if (!dataAtPts) {
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "dataAtPts is not allocated. Please set it before "
                      "using this operator.");
  }
}
 
template <int STRIDEMATD>
MoFEMErrorCode
HMHHencky::OpCalculateHelmholtzFreeEnergy<STRIDEMATD>::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 = dataAtPts->getFTensorLogStretchTotal(nb_integration_pts);
 
#ifndef NDEBUG
  auto &mat_d = dataAtPts->matD;
  if (mat_d.size2() != size_symm * size_symm) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "wrong matD size, number of columns should be %d but is %zu",
            size_symm * size_symm, mat_d.size2());
  }
  if constexpr (STRIDEMATD != 0) {
    if (mat_d.size1() != nb_integration_pts) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "wrong matD size, number of rows should be %d but is %zu",
              nb_integration_pts, mat_d.size1());
    }
  }
#endif
 
  auto t_D = getFTensor4DdgFromMat<3, 3, STRIDEMATD>(dataAtPts->matD);
 
  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_D;
    ++t_log_u;
    ++t_energy;
  }
 
  if (totalHelmholtzFreeEnergyPtr) {
    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;
    }
    *totalHelmholtzFreeEnergyPtr += getMeasure() * loc_energy;
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int STRIDEMATD>
HMHHencky::OpCalculateStretchFromStress<STRIDEMATD>::
    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(H1, OPLAST), dataAtPts(data_ptr),
      strainPtr(strain_ptr), stressPtr(stress_ptr), henckyPtr(hencky_ptr) {
  std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
  doEntities[MBVERTEX] = true;
}
 
template <int STRIDEMATD>
MoFEMErrorCode HMHHencky::OpCalculateStretchFromStress<STRIDEMATD>::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->size1();
#ifndef NDEBUG
  if (nb_integration_pts != getGaussPts().size2()) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "inconsistent number of integration points");
  }
#endif // NDEBUG
 
  CHKERR henckyPtr->computeMaterialParamsAtPts<STRIDEMATD>(this, data,
                                                           dataAtPts);
 
  auto get_strain =
      MatrixSizeHelper<GetFTensor2SymmetricFromMatType<3, -1, DL>, DL>::size(
          *strainPtr, nb_integration_pts);
  auto t_strain = get_strain();
  auto t_stress = getFTensor2FromMat<SPACE_DIM, SPACE_DIM, -1, DL>(*stressPtr);
  auto t_inv_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, STRIDEMATD>(
      dataAtPts->matInvD);
#ifndef NDEBUG
  auto t_D =
      getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, STRIDEMATD>(dataAtPts->matD);
#endif
 
  // 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::physicalTimeFlg) {
    ts_time = EshelbianCore::currentPhysicalTime;
    ts_time_step = EshelbianCore::physicalDt;
  }
 
  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);
};
 
template <typename OP_PTR>
VectorDouble getAnalyticalElastic(OP_PTR op_ptr, 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::physicalTimeFlg) {
    ts_time = EshelbianCore::currentPhysicalTime;
    ts_time_step = EshelbianCore::physicalStepNumber;
  }
 
  MatrixDouble m_ref_coords = op_ptr->getCoordsAtGaussPts();
 
  auto v_analytical_expr = analytical_elastic_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 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,
    boost::shared_ptr<double> J_ptr)
    : OpAssembleVolume(field_name, data_ptr, OPROW), alphaU(alpha_u),
      JPtr(J_ptr), topoDataPtr(topo_ptr), assembleVec(assemble_vec) {}
 
MoFEMErrorCode HMHHencky::OpTopoSpatialPhysical::integrate(EntData &data) {
  MoFEMFunctionBegin;
  CHKERR integrateHencky(data);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
HMHHencky::OpTopoSpatialPhysical::integrateHencky(EntData &data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
  auto t_L = FTensor::SymmLTensor<double, 3>();
 
  int nb_dofs = data.getIndices().size();
  int nb_integration_pts = data.getN().size1();
  auto v = getVolume();
 
  FTENSOR_INDEXES(SPACE_DIM, i, j, k, l, I, J);
 
  auto get_ftensor1 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, SPACE_DIM>, SPACE_DIM>(
        &v[0], &v[1], &v[2]);
  };
 
  locJ = 0;
  int nb_base_functions = data.getN().size2();
 
  auto integrate = [&](auto t_D) {
    MoFEMFunctionBegin;
 
    auto t_w = getFTensor0IntegrationWeight();
    auto t_det = topoDataPtr->getFTensorDetJacobian(nb_integration_pts);
    auto t_inv_jac = topoDataPtr->getFTensorInvJacobian(nb_integration_pts);
 
    auto t_var_log_u = dataAtPts->getFTensorVarLogStreach(nb_integration_pts);
    auto t_approx_P = dataAtPts->getFTensorApproxP(nb_integration_pts);
    auto t_approx_P_adjoint_log_du =
        dataAtPts->getFTensorAdjointPdU(nb_integration_pts);
    auto t_h_dlog_u =
        dataAtPts->getFTensorSmallHdLogStretch(nb_integration_pts);
    auto t_log_stretch_h1 =
        dataAtPts->getFTensorLogStretchTotal(nb_integration_pts);
 
    // Swithc off rate effects for now, as we don't have the dot_log_u at the
    // pts; we would need to compute it from the solution, which is not
    // straightforward
    // auto t_dot_log_u =
    // dataAtPts->getFTensorLogStretchDot(nb_integration_pts);
 
    auto t_diff_base = data.getFTensor1DiffN<SPACE_DIM>();
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      const double a = v * t_w;
 
      FTensor::Tensor2_symmetric<double, SPACE_DIM> 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_stress_residual;
      t_stress_residual(L) = t_L(i, j, L) * t_T(i, j);
 
      FTensor::Tensor1<double, size_symm> t_residual;
      t_residual(L) = t_approx_P_adjoint_log_du(L) - t_stress_residual(L);
 
      locJ -= (a * t_det) * t_residual(L) * t_var_log_u(L);
 
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_cof;
      t_cof(I, J) = t_det * t_inv_jac(J, I);
 
      const double var_stress_residual = t_var_log_u(L) * t_stress_residual(L);
 
      // detJ cancels the 1 / detJ from the Piola transform; the remaining
      // derivative is from J(k, n) * P(i, n).
      FTensor::Tensor3<double, size_symm, SPACE_DIM, SPACE_DIM>
          t_approx_P_adjoint_log_du_dX;
      t_approx_P_adjoint_log_du_dX(L, I, J) =
          t_h_dlog_u(i, I, L) * t_approx_P(i, J);
 
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_residual_dX;
      t_residual_dX(I, J) =
          t_var_log_u(L) * t_approx_P_adjoint_log_du_dX(L, I, J) -
          var_stress_residual * t_cof(I, J);
 
      auto t_nf = get_ftensor1(nF);
      int bb = 0;
      for (; bb != nb_dofs / SPACE_DIM; ++bb) {
        t_nf(i) -= a * t_residual_dX(i, j) * t_diff_base(j);
        ++t_nf;
        ++t_diff_base;
      }
      for (; bb != nb_base_functions; ++bb)
        ++t_diff_base;
 
      ++t_w;
      ++t_det;
      ++t_inv_jac;
      ++t_var_log_u;
      ++t_approx_P;
      ++t_approx_P_adjoint_log_du;
      ++t_h_dlog_u;
      ++t_log_stretch_h1;
      // ++t_dot_log_u;
      ++t_D;
    }
    MoFEMFunctionReturn(0);
  };
 
  if (!dataAtPts->physicsPtr->getFeatures().test(
          PhysicalEquations::NON_HOMOGENEOUS_MATERIAL)) {
    CHKERR integrate(
        getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(dataAtPts->matD));
  } else {
    CHKERR integrate(
        getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->matD));
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode HMHHencky::OpTopoSpatialPhysical::assemble(int,
                                                          EntityType row_type,
                                                          EntData &data) {
  MoFEMFunctionBegin;
  if (assembleVec) {
    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);
  }
  if (row_type == MBVERTEX) {
    if (JPtr) {
      *JPtr += locJ;
    }
  }
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity