EshelbianMatCore.cpp

Go to the documentation of this file.

/**
 * @file EshelbianMatCore.cpp
 * @author your name (you@domain.com)
 * @brief
 * @version 0.1
 * @date 2026-04-14
 *
 * @copyright Copyright (c) 2026
 *
 */
 
namespace EshelbianPlasticity {
 
struct MatPhysicalEquations : public PhysicalEquations {
 
  MatPhysicalEquations(
      boost::shared_ptr<MatOps::MatPiolaResponse> mat_physical_equations_ptr,
      const MaterialModel material_model, const Features features)
      : PhysicalEquations(material_model, features),
        matPhysicalEquationsPtr(mat_physical_equations_ptr) {
    if (!matPhysicalEquationsPtr)
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Mat core physical equations are not allocated");
    CHK_THROW_MESSAGE(getOptions(), "get options failed");
  }
 
  MoFEMErrorCode getOptions() {
    MoFEMFunctionBegin;
    PetscOptionsBegin(PETSC_COMM_WORLD, "meta_", "", "none");
    alphaGradU = 0;
    CHKERR PetscOptionsScalar("-viscosity_alpha_grad_u",
                              "Logarithmic-stretch-gradient rate viscosity", "",
                              alphaGradU, &alphaGradU, PETSC_NULLPTR);
    PetscOptionsEnd();
    MOFEM_LOG("EP", Sev::inform)
        << "alphaGradU (-meta_viscosity_alpha_grad_u), "
           "logarithmic-stretch-gradient rate viscosity: "
        << alphaGradU;
    MoFEMFunctionReturn(0);
  }
 
  UserDataOperator *
  returnOpJacobian(const bool eval_rhs, const bool eval_lhs,
                   boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                   boost::shared_ptr<PhysicalEquations> physics_ptr) override {
    return matPhysicalEquationsPtr->createOp(matPhysicalEquationsPtr, eval_rhs,
                                             eval_lhs, false);
  }
 
  struct OpCalculateHelmholtzFreeEnergy : public VolUserDataOperator {
 
    OpCalculateHelmholtzFreeEnergy(
        boost::shared_ptr<DataAtIntegrationPts> data_ptr,
        boost::shared_ptr<MatOps::MatElastic> physics_ptr,
        boost::shared_ptr<double> total_helmholtz_free_energy_ptr)
        : VolUserDataOperator(NOSPACE, OPSPACE), dataAtPts(data_ptr),
          physicsPtr(physics_ptr),
          totalHelmholtzFreeEnergyPtr(total_helmholtz_free_energy_ptr) {
      if (!dataAtPts)
        CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                          "Mat core energy data is not allocated");
      if (!physicsPtr)
        CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                          "Mat core physical equations are not allocated");
    }
 
    MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
  private:
    boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
    boost::shared_ptr<MatOps::MatElastic> physicsPtr;
    boost::shared_ptr<double> totalHelmholtzFreeEnergyPtr;
  };
 
  VolUserDataOperator *
  returnOpSpatialPhysical(const std::string &field_name,
                          boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                          const double alpha) override {
    return new OpSpatialPhysical(field_name, matPhysicalEquationsPtr, alpha,
                                 alphaGradU);
  }
 
  VolUserDataOperator *returnOpSpatialPhysical_du_du(
      std::string row_field, std::string col_field,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      const double alpha) override {
    return new OpSpatialPhysical_du_du(
        row_field, col_field, matPhysicalEquationsPtr, alpha, alphaGradU);
  }
 
  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);
  }
 
protected:
  boost::shared_ptr<MatOps::MatPiolaResponse> matPhysicalEquationsPtr;
  double alphaGradU = 0;
 
  struct OpSpatialPhysical : public OpAssembleVolume {
 
    OpSpatialPhysical(const std::string &field_name,
                      boost::weak_ptr<MatOps::PhysicalEquations> physics_ptr,
                      const double alpha, const double alpha_grad_u)
        : OpAssembleVolume(field_name, nullptr, OPROW), alphaU(alpha),
          alphaGradU(alpha_grad_u), physicsPtr(physics_ptr) {}
 
    MoFEMErrorCode integrate(EntData &data);
 
  private:
    const double alphaU;
    const double alphaGradU;
    boost::shared_ptr<MatOps::PhysicalEquations>
        physicsPtr; ///< material physical equations
  };
 
  struct OpSpatialPhysical_du_du : public OpAssembleVolume {
    OpSpatialPhysical_du_du(
        std::string row_field, std::string col_field,
        boost::shared_ptr<MatOps::PhysicalEquations> physics_ptr,
        const double alpha, const double alpha_grad_u)
        : OpAssembleVolume(row_field, col_field, nullptr, OPROWCOL, false),
          alphaU(alpha), alphaGradU(alpha_grad_u), physicsPtr(physics_ptr) {
      sYmm = false;
    }
 
    MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
 
  private:
    const double alphaU;
    const double alphaGradU;
    boost::shared_ptr<MatOps::PhysicalEquations>
        physicsPtr; ///< material physical equations
  };
};
 
struct MatElasticPhysicalEquations : public MatPhysicalEquations {
 
  MatElasticPhysicalEquations(
      boost::shared_ptr<MatOps::MatElastic> mat_elastic_ptr,
      const MaterialModel material_model, const Features features)
      : MatPhysicalEquations(mat_elastic_ptr, material_model, features),
        matElasticPtr(mat_elastic_ptr) {
    if (!matElasticPtr)
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Mat core elastic equations are not allocated");
  }
 
  bool providesHelmholtzFreeEnergy() const override { return true; }
 
  VolUserDataOperator *
  returnOpCalculateHelmholtzFreeEnergy(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<double> total_helmholtz_free_energy_ptr) override {
    return new OpCalculateHelmholtzFreeEnergy(
        data_ptr, matElasticPtr, total_helmholtz_free_energy_ptr);
  }
 
private:
  boost::shared_ptr<MatOps::MatElastic> matElasticPtr;
};
 
MoFEMErrorCode
MatPhysicalEquations::OpCalculateHelmholtzFreeEnergy::doWork(
    int side, EntityType type, EntData &data) {
  MoFEMFunctionBegin;
  (void)side;
  (void)type;
  (void)data;
 
  auto mat_ops_data = physicsPtr->matOpsDataPtr;
  if (!mat_ops_data)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Mat core material data is not allocated");
 
  const int nb_integration_pts = getGaussPts().size2();
  if (nb_integration_pts <= 0)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Mat core energy operator has no integration points");
 
  auto mat_grad_ptr = mat_ops_data->getCommonDataPtr("grad");
  auto mat_F_ptr = mat_ops_data->getActiveDataPtr("F");
  auto mat_energy_ptr = physicsPtr->getEnergyDensityPtr();
  if (!mat_grad_ptr || !mat_F_ptr || !mat_energy_ptr)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Mat core energy evaluation has incomplete material data");
  if (mat_grad_ptr->size1() != nb_integration_pts ||
      mat_grad_ptr->size2() != SPACE_DIM * SPACE_DIM)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Mat core energy gradient data has size %zu x %zu; expected %d "
            "x %d",
            mat_grad_ptr->size1(), mat_grad_ptr->size2(), nb_integration_pts,
            SPACE_DIM * SPACE_DIM);
 
  mat_F_ptr->resize(SPACE_DIM, SPACE_DIM, false);
  mat_energy_ptr->resize(1, 1, false);
 
  const auto get_material_tag = [&]() {
    if (physicsPtr->tagVsRangePtr) {
      for (const auto &tag_range_pair : *physicsPtr->tagVsRangePtr) {
        if (tag_range_pair.second.find(getFEEntityHandle()) !=
            tag_range_pair.second.end())
          return tag_range_pair.first;
      }
    }
#ifndef NDEBUG
    if (MatOps::MatOpsTagsRegistry::getTagName(physicsPtr->tAg).empty())
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                        "ADOL-C tag not found " +
                            std::to_string(physicsPtr->tAg));
#endif
    return physicsPtr->tAg;
  };
 
  const int material_tag = get_material_tag();
  auto *fe_ptr = const_cast<FEMethod *>(getFEMethod());
  if (!fe_ptr)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Mat core energy operator has no finite-element method");
  const EntityHandle ent = fe_ptr->getFEEntityHandle();
  if (!ent)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Mat core energy operator has no finite-element entity");
 
  using GaussLayout = DataLayoutTraits<DataLayout::GaussByCoeffs>;
  auto get_grad_at_pts = MatrixSizeHelper<
      GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, GaussLayout>,
      GaussLayout>::get(*mat_grad_ptr, nb_integration_pts);
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, J);
  auto t_grad_at_pts = get_grad_at_pts();
  dataAtPts->energyAtPts.resize(nb_integration_pts, false);
  auto t_energy_at_pts = getFTensor0FromVec(dataAtPts->energyAtPts);
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    auto t_F =
        getFTensor2FromPtr<SPACE_DIM, SPACE_DIM>(mat_F_ptr->data().data());
    t_F(i, J) = t_grad_at_pts(i, J);
    CHKERR physicsPtr->setParams(fe_ptr, gg);
    CHKERR physicsPtr->evaluateVariable(material_tag, ent, gg);
 
    auto t_energy = getFTensor0FromMat(mat_energy_ptr);
    const double energy_density = t_energy;
    if (!std::isfinite(energy_density))
      SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP,
              "Mat core returned non-finite energy density at point %d on "
              "entity %lu",
              gg, static_cast<unsigned long>(ent));
    t_energy_at_pts = energy_density;
 
    ++t_grad_at_pts;
    ++t_energy_at_pts;
  }
 
  if (totalHelmholtzFreeEnergyPtr) {
    if (!std::isfinite(*totalHelmholtzFreeEnergyPtr))
      SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP,
              "Mat core accumulated energy is non-finite before integration");
    const double element_measure = getMeasure();
    if (!std::isfinite(element_measure) || element_measure < 0.)
      SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP,
              "Mat core element measure is invalid: %.16g", element_measure);
 
    auto t_w = getFTensor0IntegrationWeight();
    auto t_energy_at_pts_integral = getFTensor0FromVec(dataAtPts->energyAtPts);
    auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
    double local_energy = 0;
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      const double weight = t_w;
      const double energy_density = t_energy_at_pts_integral;
      const double det_plasticF = determinantTensor3by3(t_plasticF);
      if (!std::isfinite(weight) || !std::isfinite(det_plasticF) ||
          det_plasticF <= 0.)
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP,
                "Invalid Mat core energy quadrature data at point %d on "
                "entity %lu",
                gg, static_cast<unsigned long>(ent));
      local_energy += weight * det_plasticF * energy_density;
 
      ++t_w;
      ++t_energy_at_pts_integral;
      ++t_plasticF;
    }
 
    const double element_energy = element_measure * local_energy;
    const double accumulated_energy =
        *totalHelmholtzFreeEnergyPtr + element_energy;
    if (!std::isfinite(element_energy) || !std::isfinite(accumulated_energy))
      SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP,
              "Mat core integrated energy is non-finite on entity %lu",
              static_cast<unsigned long>(ent));
    *totalHelmholtzFreeEnergyPtr = accumulated_energy;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
MatPhysicalEquations::OpSpatialPhysical::integrate(EntData &data) {
  MoFEMFunctionBegin;
 
  auto physics_ptr = physicsPtr;
  auto mat_ops_data = physics_ptr->matOpsDataPtr;
 
  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_P = getFTensor2FromMat<3, 3, -1, DL>(
      mat_ops_data->getCommonDataPtr("PAtPts"));
  auto t_approx_P_adjoint_log_du = getFTensor1FromMat<size_symm, -1, DL>(
      mat_ops_data->getCommonDataPtr("adjointPdUAtPts"));
  auto t_diff_u = getFTensor4FromMat<3, 3, 3, 3, -1, DL>(
      mat_ops_data->getCommonDataPtr("diffStretchH1AtPts"));
  auto t_dot_log_u = getFTensor2SymmetricFromMat<3, -1, DL>(
      mat_ops_data->getCommonDataPtr("logStretchDotTensorAtPts"));
  auto t_grad_dot_log_u = getFTensor2FromMat<size_symm, 3, -1, DL>(
      mat_ops_data->getCommonDataPtr("gradLogStretchDotTensorAtPts"));
  auto t_plasticF = getFTensor2SymmetricFromMat<3, -1, DL>(
      mat_ops_data->getCommonDataPtr("plasticF"));
  auto t_invPlasticF = getFTensor2FromMat<3, 3, -1, DL>(
      mat_ops_data->getCommonDataPtr("invPlasticF"));
 
  FTENSOR_INDEX(3, i);
  FTENSOR_INDEX(3, j);
  FTENSOR_INDEX(3, k);
  FTENSOR_INDEX(3, l);
  FTENSOR_INDEX(3, m);
  FTENSOR_INDEX(3, n);
 
  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();
  auto t_row_grad_fun = data.getFTensor1DiffN<3>();
  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::Tensor3<double, 3, 3, size_symm> t_Ldiff_u;
    t_Ldiff_u(i, j, L) = t_diff_u(i, j, k, l) * t_L(k, l, L);
 
    FTensor::Tensor1<double, size_symm> t_residual;
    t_residual(L) =
        a * (t_approx_P_adjoint_log_du(L) - t_Ldiff_u(i, j, L) * t_P(i, j));
    t_residual(L) += (a * alphaU) * (t_dot_log_u(i, j) * t_L(i, j, L));
 
    FTensor::Tensor2<double, size_symm, 3> t_grad_dot_log_u_intermediate;
    t_grad_dot_log_u_intermediate(L, j) =
        t_grad_dot_log_u(L, i) * t_invPlasticF(i, j);
    FTensor::Tensor2<double, size_symm, 3> t_grad_residual;
    t_grad_residual(L, j) =
        (a * alphaGradU) * t_grad_dot_log_u_intermediate(L, j);
 
    int bb = 0;
    for (; bb != nb_dofs / size_symm; ++bb) {
      FTensor::Tensor1<double, 3> t_row_grad_intermediate;
      t_row_grad_intermediate(j) = t_row_grad_fun(i) * t_invPlasticF(i, j);
      t_nf(L) -= t_row_base_fun * t_residual(L);
      t_nf(L) += t_row_grad_intermediate(j) * t_grad_residual(L, j);
      ++t_nf;
      ++t_row_base_fun;
      ++t_row_grad_fun;
    }
    for (; bb != nb_base_functions; ++bb) {
      ++t_row_base_fun;
      ++t_row_grad_fun;
    }
 
    ++t_w;
    ++t_approx_P_adjoint_log_du;
    ++t_P;
    ++t_diff_u;
    ++t_dot_log_u;
    ++t_grad_dot_log_u;
    ++t_plasticF;
    ++t_invPlasticF;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
MatPhysicalEquations::OpSpatialPhysical_du_du::integrate(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>();
  constexpr auto t_kd_sym = FTensor::Kronecker_Delta_symmetric<int>();
 
  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(3, i);
  FTENSOR_INDEX(3, j);
  FTENSOR_INDEX(3, k);
  FTENSOR_INDEX(3, l);
  FTENSOR_INDEX(3, m);
  FTENSOR_INDEX(3, n);
 
  MatrixDouble dP;
 
  auto get_dP = [&]() {
    auto get_tensor_for_dP =
        MatrixSizeHelper<GetFTensor2FromMatType<size_symm, size_symm, -1, DL>,
                         DL>::size(dP, nb_integration_pts);
 
    auto mat_ops_data = physicsPtr->matOpsDataPtr;
 
    auto t_P = getFTensor2FromMat<3, 3, -1, DL>(
        mat_ops_data->getCommonDataPtr("PAtPts"));
    auto t_P_du = getFTensor4FromMat<3, 3, 3, 3, -1, DL>(
        mat_ops_data->getCommonDataPtr("PAtPts_du"));
    auto t_approx_P_adjoint_dstretch = getFTensor2FromMat<3, 3, -1, DL>(
        mat_ops_data->getCommonDataPtr("adjointPdstretchAtPts"));
    auto t_u = getFTensor2FromMat<3, 3, -1, DL>(
        mat_ops_data->getCommonDataPtr("stretchH1AtPts"));
    auto t_diff_u = getFTensor4FromMat<3, 3, 3, 3, -1, DL>(
        mat_ops_data->getCommonDataPtr("diffStretchH1AtPts"));
    auto t_grad_h1 = getFTensor2FromMat<3, 3, -1, DL>(
        mat_ops_data->getCommonDataPtr("wGradH1AtPts"));
    auto t_eigen_vals = getFTensor1FromMat<3, -1, DL>(
        mat_ops_data->getCommonDataPtr("eigenVals"));
    auto t_eigen_vecs = getFTensor2FromMat<3, 3, -1, DL>(
        mat_ops_data->getCommonDataPtr("eigenVecs"));
    auto t_invPlasticF = getFTensor2FromMat<3, 3, -1, DL>(
        mat_ops_data->getCommonDataPtr("invPlasticF"));
 
    auto ts_a = getTSa();
    auto t_kd = FTensor::Kronecker_Delta<double>();
 
    auto t_dP = get_tensor_for_dP();
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
 
      FTensor::Tensor2<double, 3, 3> t_h1;
 
      switch (EshelbianCore::gradApproximator) {
      case LARGE_ROT:
      case MODERATE_ROT:
        t_h1(i, j) = t_grad_h1(i, k) * t_invPlasticF(k, j) + t_kd(i, j);
        break;
      default:
        t_h1(i, j) = t_kd(i, j);
      }
 
      FTensor::Tensor3<double, 3, 3, size_symm> t_Ldiff_u;
      t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
      t_dP(L, J) =
          -t_Ldiff_u(i, j, L) * (t_P_du(i, j, k, l) * t_Ldiff_u(k, l, J));
      t_dP(L, J) += (ts_a * alphaU) *
                    (t_L(i, j, L) * (t_diff(i, j, k, l) * t_L(k, l, J)));
 
      if (EshelbianCore::stretchSelector > LINEAR &&
          EshelbianCore::gradApproximator > MODERATE_ROT) {
        FTensor::Tensor2<double, 3, 3> t_deltaP;
        t_deltaP(i, j) =
            t_approx_P_adjoint_dstretch(i, j) - t_P(i, k) * t_h1(j, k);
        FTensor::Tensor2_symmetric<double, 3> t_deltaP_sym;
        t_deltaP_sym(i, j) = (t_deltaP(i, j) || t_deltaP(j, i));
        t_deltaP_sym(i, j) /= 2.0;
        auto nb_uniq = getUniqNb<3>(getVectorAdaptor(&t_eigen_vals(0), 3),
                                    FTensor::Number<3>{});
        auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
            t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
            EshelbianCore::dd_f, t_deltaP_sym, nb_uniq);
        t_dP(L, J) += t_L(i, j, L) * (t_diff2_uP2(i, j, k, l) * t_L(k, l, J));
      }
 
      ++t_P;
      ++t_dP;
      ++t_P_du;
      ++t_approx_P_adjoint_dstretch;
      ++t_u;
      ++t_diff_u;
      ++t_grad_h1;
      ++t_eigen_vals;
      ++t_eigen_vecs;
      ++t_invPlasticF;
    }
 
    return get_tensor_for_dP();
  };
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
  auto t_row_grad_fun = row_data.getFTensor1DiffN<3>();
 
  auto t_dP = get_dP();
  auto v = getVolume();
  auto ts_a = getTSa();
  auto t_w = getFTensor0IntegrationWeight();
  auto mat_ops_data = physicsPtr->matOpsDataPtr;
  auto t_plasticF = getFTensor2SymmetricFromMat<3, -1, DL>(
      mat_ops_data->getCommonDataPtr("plasticF"));
  auto t_invPlasticF = getFTensor2FromMat<3, 3, -1, DL>(
      mat_ops_data->getCommonDataPtr("invPlasticF"));
 
  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 / size_symm; ++rr) {
      FTensor::Tensor1<double, 3> t_row_grad_intermediate;
      t_row_grad_intermediate(j) = t_row_grad_fun(i) * t_invPlasticF(i, j);
      auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
      auto t_col_grad_fun = col_data.getFTensor1DiffN<3>(gg, 0);
      auto t_m = get_ftensor2(K, size_symm * rr, 0);
      for (int cc = 0; cc != col_nb_dofs / size_symm; ++cc) {
        FTensor::Tensor1<double, 3> t_col_grad_intermediate;
        t_col_grad_intermediate(j) = t_col_grad_fun(i) * t_invPlasticF(i, j);
        const double b = a * t_row_base_fun * t_col_base_fun;
        t_m(L, J) -= b * t_dP(L, J);
        const double c = (a * alphaGradU * ts_a) * (t_row_grad_intermediate(j) *
                                                    t_col_grad_intermediate(j));
        t_m(L, J) += c * t_kd_sym(L, J);
        ++t_m;
        ++t_col_base_fun;
        ++t_col_grad_fun;
      }
      ++t_row_base_fun;
      ++t_row_grad_fun;
    }
 
    for (; rr != row_nb_base_functions; ++rr) {
      ++t_row_base_fun;
      ++t_row_grad_fun;
    }
 
    ++t_w;
    ++t_dP;
    ++t_plasticF;
    ++t_invPlasticF;
  }
  MoFEMFunctionReturn(0);
}
 
MatPhysicalEquations::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 MatPhysicalEquations::OpSpatialPhysicalExternalStrain::integrate(
    EntData &data) {
  MoFEMFunctionBegin;
 
  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)) {
    auto block_name = "(.*)ANALYTICAL_EXTERNALSTRAIN(.*)";
    std::regex reg_name(block_name);
    if (std::regex_match(ext_strain_block.blockName, reg_name)) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
              "Analytical external strain not implemented for Neo-Hookean "
              "material.");
    }
    // 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;
      }
 
      // get ExternalStrain block data
      double external_strain_val = scale * ext_strain_block.val;
      double K = ext_strain_block.bulkModulusK;
 
      FTensor::Index<'L', size_symm> L;
      auto t_L = FTensor::SymmLTensor<double, 3>();
      constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
 
      int nb_dofs = data.getIndices().size();
      int nb_integration_pts = data.getN().size1();
      auto vol = getVolume();
      auto t_w = getFTensor0IntegrationWeight();
      auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
 
      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);
 
      int nb_base_functions = data.getN().size2();
      auto t_row_base_fun = data.getFTensor0N();
 
      const double bulk_modulus = K;
      const double diag_val =
          EshelbianCore::inv_f(std::exp(external_strain_val));
      auto fun_neohookean_bulk = [](double K, double tr) { return K * tr; };
      const double sigma_J = fun_neohookean_bulk(bulk_modulus, 3 * diag_val);
 
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        const double det_plasticF = determinantTensor3by3(t_plasticF);
        const double a = vol * t_w * det_plasticF;
        ++t_w;
        ++t_plasticF;
 
        FTensor::Tensor1<double, size_symm> t_residual;
        t_residual(L) = (t_L(i, j, L) * t_kd(i, j)) * sigma_J;
        t_residual(L) *= a;
 
        auto t_nf = getFTensor1FromPtr<size_symm>(&*nF.data().begin());
        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;
        }
      }
    }
  }
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity