HMHNeohookean.cpp

Go to the documentation of this file.

/**
 * @file HMHNeohookean.cpp
 * @brief Implementation of NeoHookean material
 * @date 2024-08-31
 *
 * @copyright Copyright (c) 2024
 *
 */
 
namespace EshelbianPlasticity {
 
// #define NEOHOOKEAN_SCALING
struct HMHNeohookean : public PhysicalEquations {
 
  /**
   * @brief Definition of Neo-hookean function
   * 
   * P(i, I) = 2. * c10_default * (tF(i, I) - tInvF(i, I)) +
   *               K_default * log(detF) * tInvF(i, I)
   * tCof(i, I) = detF * tInvF(I, i)
   * Psi = c10* (trace(tCof) + 2ln(det(F)) - 3) + (K/2.) * (log(detF))^2
   * 
   * @param c10
   * @param v
   * @return double
   */
  static inline double fun_neohookean(double c10, double v) {
    return (2.0 * c10) * (-1. + EshelbianCore::f(2 * v));
  }
 
  /**
   * @brief Definition of derivative of Neo-hookean function
   * 
   * @param c10
   * @param v
   * @return * double
   */
  static inline double fun_d_neohookean(double c10, double v) {
    return
 
        (4.0 * c10) * EshelbianCore::d_f(2 * v);
  }
 
  /**
   * @brief Definition of axiator of Neo-hookean function
   * 
   * Psi = (K/2.) * (log(J))^2 = (K/2.) * (log(exp tr H))^2 = (K/2.) * (tr H)^2
   *
   * @param c10
   * @param v
   * @return double
   */
  static inline double fun_neohookean_bulk(double K, double tr) {
    return K * tr; 
  }
 
  /**
   * @brief Definition of derivative of axiator of Neo-hookean function
   *
   * @param K
   * @param tr
   * @return double
   */
  static inline double fun_diff_neohookean_bulk(double K, double tr) {
    return K;//K * exp(tr) * (1. + tr); // K*2*ln(J)/J
  }
 
  static constexpr int numberOfActiveVariables = 9;
  static constexpr int numberOfDependentVariables = 9;
 
  HMHNeohookean(MoFEM::Interface &m_field, const double c10, const double K)
      : PhysicalEquations(numberOfActiveVariables, numberOfDependentVariables),
        mField(m_field), c10_default(c10), K_default(K) {
 
    CHK_THROW_MESSAGE(getOptions(), "get options failed");
    CHK_THROW_MESSAGE(extractBlockData(Sev::inform),
                      "extract block data failed");
 
#ifdef NEOHOOKEAN_SCALING
    if (blockData.size()) {
      double ave_K = 0;
      for (auto b : blockData) {
        ave_K += b.K;
      }
      eqScaling = ave_K / blockData.size();
    }
#endif
 
    MOFEM_LOG("EP", Sev::inform) << "Neo-Hookean scaling = " << eqScaling;
 
    if (adolCOn == PETSC_FALSE) {
      if (EshelbianCore::stretchSelector < StretchSelector::LOG) {
        CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                          "Stretch selector is not equal to LOG");
      } else {
        if (EshelbianCore::exponentBase != exp(1)) {
          CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                            "Exponent base is not equal to exp(1)");
        }
      }
    }
  }
 
  auto getMaterialParameters(EntityHandle ent) {
    for (auto &b : blockData) {
      if (b.blockEnts.find(ent) != b.blockEnts.end()) {
        return std::make_pair(b.c10, b.K);
      }
    }
    if (blockData.size() != 0)
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Block not found for entity handle. If you mat set "
                        "block, set it to all elements");
    return std::make_pair(c10_default, K_default);
  }
 
  MoFEMErrorCode setParams(const int tag, EntityHandle ent, int gg) {
    MoFEMFunctionBegin;
    auto [def_c10, def_K] = getMaterialParameters(ent);
    std::vector<double> params = {def_c10, def_K};
    set_param_vec(tag, params.size(), params.data());
    MoFEMFunctionReturn(0);
  };
 
  struct OpGetScale : public VolUserDataOperator {
    OpGetScale(boost::shared_ptr<double> scale_ptr,
               boost::shared_ptr<PhysicalEquations> physics_ptr)
        : VolUserDataOperator(NOSPACE, VolUserDataOperator::OPSPACE),
          scalePtr(scale_ptr), physicsPtr(physics_ptr) {}
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBegin;
      auto neohookean_ptr =
          boost::dynamic_pointer_cast<HMHNeohookean>(physicsPtr);
      *(scalePtr) = neohookean_ptr->eqScaling;
      MoFEMFunctionReturn(0);
    }
 
  private:
    boost::shared_ptr<double> scalePtr;
    boost::shared_ptr<PhysicalEquations> physicsPtr;
  };
 
  VolUserDataOperator *
  returnOpSetScale(boost::shared_ptr<double> scale_ptr,
                   boost::shared_ptr<PhysicalEquations> physics_ptr) {
 
    return new OpGetScale(scale_ptr, physics_ptr);
  };
 
  struct OpJacobian : public EshelbianPlasticity::OpJacobian {
    using EshelbianPlasticity::OpJacobian::OpJacobian;
    MoFEMErrorCode evaluateRhs(EntData &data) { return 0; }
    MoFEMErrorCode evaluateLhs(EntData &data) { return 0; }
  };
 
  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) {
    if (adolCOn) {
      return PhysicalEquations::returnOpJacobian(tag, eval_rhs, eval_lhs,
                                                 data_ptr, physics_ptr);
    } else {
      return (new OpJacobian(tag, eval_rhs, eval_lhs, data_ptr, physics_ptr));
    }
  }
 
  MoFEMErrorCode getOptions() {
    MoFEMFunctionBegin;
    PetscOptionsBegin(PETSC_COMM_WORLD, "neo_hookean_", "", "none");
 
    c10_default = 1;
    CHKERR PetscOptionsScalar("-c10", "C10", "", c10_default, &c10_default,
                              PETSC_NULLPTR);
    K_default = 1;
    CHKERR PetscOptionsScalar("-K", "Bulk modulus K", "", K_default, &K_default,
                              PETSC_NULLPTR);
 
    alphaGradU = 0;
    CHKERR PetscOptionsScalar("-viscosity_alpha_grad_u", "viscosity", "",
                              alphaGradU, &alphaGradU, PETSC_NULLPTR);
 
    CHKERR PetscOptionsBool("-adolc_on", "adolc ON", "", adolCOn, &adolCOn,
                            PETSC_NULLPTR);
 
    PetscOptionsEnd();
 
    MOFEM_LOG_CHANNEL("WORLD");
    MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "MatBlock Neo-Hookean (default)")
        << "c10 = " << c10_default << " K = " << K_default
        << " grad alpha u = " << alphaGradU;
    MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "MatBlock Neo-Hookean (default)")
            << "ADOL-C On " << (adolCOn)
        ? "Yes"
        : "No";
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode extractBlockData(Sev sev) {
    return extractBlockData(
 
        mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
            (boost::format("%s(.*)") % "MAT_NEOHOOKEAN").str()
 
                )),
 
        sev);
  }
 
  MoFEMErrorCode
  extractBlockData(std::vector<const CubitMeshSets *> meshset_vec_ptr,
                   Sev sev) {
    MoFEMFunctionBegin;
 
    for (auto m : meshset_vec_ptr) {
      MOFEM_LOG("EP", sev) << *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 c10 = block_data[0];
      double K = block_data[1];
 
      blockData.push_back({c10, K, get_block_ents()});
 
      MOFEM_LOG("EP", sev) << "MatBlock Neo-Hookean c10 = "
                           << blockData.back().c10
                           << " K = " << blockData.back().K << " nb ents. = "
                           << blockData.back().blockEnts.size();
    }
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr) {
    MoFEMFunctionBegin;
 
    FTENSOR_INDEX(3, i);
    FTENSOR_INDEX(3, j);
    FTENSOR_INDEX(3, I);
    FTENSOR_INDEX(3, J);
 
    auto large = [&]() {
      MoFEMFunctionBegin;
 
      constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
 
      auto ih = get_h();
      if (t_h_ptr)
        ih(i, j) = (*t_h_ptr)(i, j);
      else {
        ih(i, j) = t_kd(i, j);
      }
 
      auto r_P = get_P();
 
      enableMinMaxUsingAbs();
      trace_on(tape);
 
      // Set active variables to ADOL-C
      th(i, j) <<= get_h()(i, j);
 
      tF(i, I) = th(i, I);
      CHKERR determinantTensor3by3(tF, detF);
      CHKERR invertTensor3by3(tF, detF, tInvF);
 
      // tCof(i, I) = detF * tInvF(I, i);
 
      // Psi = c10* (trace(tCof) + 2ln(det(F)) - 3) + (K/2.) * (log(detF))^2;
 
      auto p_c10 = mkparam(c10_default);
      auto p_K = mkparam(K_default);
      tP(i, I) =
          2. * p_c10 * (tF(i, I) - tInvF(i, I)) + p_K * log(detF) * tInvF(i, I);
 
      // Set dependent variables to ADOL-C
      tP(i, j) >>= r_P(i, j);
 
      trace_off();
 
      MoFEMFunctionReturn(0);
    };
 
    switch (EshelbianCore::gradApproximator) {
    case NO_H1_CONFIGURATION:
    case LARGE_ROT:
    case MODERATE_ROT:
      CHKERR large();
      break;
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "gradApproximator not handled");
    };
 
    MoFEMFunctionReturn(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);
 
  private:
    const double alphaU;
  };
 
  virtual VolUserDataOperator *
  returnOpSpatialPhysical(const std::string &field_name,
                          boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                          const double alpha_u) {
    if (adolCOn) {
      return PhysicalEquations::returnOpSpatialPhysical(field_name, data_ptr,
                                                        alpha_u);
    } else {
      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 OpAssembleVolumePositiveDefine {
    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);
 
  private:
    const double alphaU;
  };
 
  VolUserDataOperator *returnOpSpatialPhysical_du_du(
      std::string row_field, std::string col_field,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha) {
    if (adolCOn) {
      return PhysicalEquations::returnOpSpatialPhysical_du_du(
          row_field, col_field, data_ptr, alpha);
    } else {
      return new OpSpatialPhysical_du_du(row_field, col_field, data_ptr, alpha);
    }
  }
 
private:
  MoFEM::Interface &mField;
 
  double c10_default;
  double K_default;
  double alphaGradU;
  struct BlockData {
    double c10;
    double K;
    Range blockEnts;
  };
  std::vector<BlockData> blockData;
 
  double eqScaling = 1.;
 
  // ADOL-C
  PetscBool adolCOn = PETSC_FALSE;
 
  ATensor2 th;
  ATensor2 tF;
  adouble detF;
  ATensor2 tInvF;
  ATensor2 tCof;
  ATensor2 tP;
};
 
HMHNeohookean::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) {}
 
extern "C" {
void tetcircumcenter_tp(double a[3], double b[3], double c[3], double d[3],
                        double circumcenter[3], double *xi, double *eta,
                        double *zeta);
}
 
MoFEMErrorCode HMHNeohookean::OpSpatialPhysical::integrate(EntData &data) {
  MoFEMFunctionBegin;
 
  auto neohookean_ptr =
      boost::dynamic_pointer_cast<HMHNeohookean>(dataAtPts->physicsPtr);
  if (!neohookean_ptr) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Pointer to HMHNeohookean is null");
  }
  auto [def_c10, def_K] =
      neohookean_ptr->getMaterialParameters(getFEEntityHandle());
 
  double c10 = def_c10 / neohookean_ptr->eqScaling;
  double alpha_u = alphaU / neohookean_ptr->eqScaling;
  double K = def_K / neohookean_ptr->eqScaling;
 
  double alpha_grad_u =
      neohookean_ptr->alphaGradU / neohookean_ptr->eqScaling;
 
  FTensor::Index<'L', size_symm> L;
  auto t_L = symm_L_tensor();
 
  constexpr auto t_kd_sym = FTensor::Kronecker_Delta_symmetric<int>();
  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
 
  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_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
  auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
  auto t_dot_log_u =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
  auto t_diff_u =
      getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
  auto t_log_u =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
  auto t_grad_log_u =
      getFTensor2FromMat<6, 3>(dataAtPts->gradLogStretchDotTensorAtPts);
  auto t_log_u2_h1 =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
 
  auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
  auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
  auto &nbUniq = dataAtPts->nbUniq;
  auto t_nb_uniq =
      FTensor::Tensor0<FTensor::PackPtr<int *, 1>>(nbUniq.data().data());
 
  auto t_diff = diff_tensor();
 
  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 get_ftensor2 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, size_symm>, size_symm>(
        &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_grad_base_fun = data.getFTensor1DiffN<3>();
 
  auto no_h1 = [&]() {
    MoFEMFunctionBegin;
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      double a = v * t_w;
      ++t_w;
 
      auto neohookean = [c10](auto v) { return fun_neohookean(c10, v); };
      auto t_neohookean_hencky =
          EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, neohookean);
      ++t_eigen_vals;
      ++t_eigen_vecs;
      ++t_nb_uniq;
 
      const double tr = t_log_u(i, i);
      FTensor::Tensor1<double, size_symm> t_P;
      t_P(L) = t_L(i, j, L) * (t_neohookean_hencky(i, j) +
                               t_kd(i, j) * fun_neohookean_bulk(K, tr));
      FTensor::Tensor1<double, size_symm> t_viscous_P;
      t_viscous_P(L) = alpha_u * (t_L(i, j, L) * t_dot_log_u(i, j));
 
      FTensor::Tensor1<double, size_symm> t_residual;
      t_residual(L) = t_approx_P_adjoint_log_du(L) - t_P(L) - t_viscous_P(L);
      t_residual(L) *= a;
 
      FTensor::Tensor2<double, size_symm, SPACE_DIM> t_grad_residual;
      t_grad_residual(L, i) = alpha_grad_u * t_grad_log_u(L, i);
      t_grad_residual(L, i) *= a;
 
      ++t_approx_P_adjoint_log_du;
      ++t_log_u;
      ++t_dot_log_u;
      ++t_grad_log_u;
 
      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(L) += t_grad_base_fun(i) * t_grad_residual(L, i);
        ++t_nf;
        ++t_row_base_fun;
        ++t_grad_base_fun;
      }
      for (; bb != nb_base_functions; ++bb) {
        ++t_row_base_fun;
        ++t_grad_base_fun;
      }
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto large = [&]() {
    MoFEMFunctionBegin;
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Not implemented for Neo-Hookean (used ADOL-C)");
    MoFEMFunctionReturn(0);
  };
 
  switch (EshelbianCore::gradApproximator) {
  case NO_H1_CONFIGURATION:
    CHKERR no_h1();
    break;
  case LARGE_ROT:
  case MODERATE_ROT:
    CHKERR large();
    break;
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "gradApproximator not handled");
  };
 
   MoFEMFunctionReturn(0);
}
 
HMHNeohookean::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
HMHNeohookean::OpSpatialPhysicalExternalStrain::integrate(EntData &data) {
  MoFEMFunctionBegin;
 
  auto neohookean_ptr =
      boost::dynamic_pointer_cast<HMHNeohookean>(dataAtPts->physicsPtr);
  if (!neohookean_ptr) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Pointer to HMHNeohookean is null");
  }
 
  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)) {
    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 = symm_L_tensor();
      constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
 
      int nb_dofs = data.getIndices().size();
      int nb_integration_pts = data.getN().size1();
      auto v = getVolume();
      auto t_w = getFTensor0IntegrationWeight();
 
      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 tr = 3.0 * external_strain_val;
      const double sigma_J = K * tr;
      
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        double a = v * t_w;
        ++t_w;
 
        FTensor::Tensor1<double, size_symm> t_residual;
        t_residual(L) = 0.0;
        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);
}
 
HMHNeohookean::OpSpatialPhysical_du_du::OpSpatialPhysical_du_du(
    std::string row_field, std::string col_field,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha)
    : OpAssembleVolumePositiveDefine(row_field, col_field, data_ptr, OPROWCOL,
                                     false),
      alphaU(alpha) {
  sYmm = false;
}
 
MoFEMErrorCode
HMHNeohookean::OpSpatialPhysical_du_du::integrate(EntData &row_data,
                                                  EntData &col_data) {
  MoFEMFunctionBegin;
 
  auto neohookean_ptr =
      boost::dynamic_pointer_cast<HMHNeohookean>(dataAtPts->physicsPtr);
  if (!neohookean_ptr) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Pointer to HMHNeohookean is null");
  }
  auto [def_c10, def_K] =
      neohookean_ptr->getMaterialParameters(getFEEntityHandle());
 
  double c10 = def_c10 / neohookean_ptr->eqScaling;
  double alpha_u = alphaU / neohookean_ptr->eqScaling;
  double lambda = def_K / neohookean_ptr->eqScaling;
 
  double alpha_grad_u =
      neohookean_ptr->alphaGradU / neohookean_ptr->eqScaling;
 
  FTENSOR_INDEX(size_symm, L);
  FTENSOR_INDEX(size_symm, J);
 
  constexpr auto t_kd_sym = FTensor::Kronecker_Delta_symmetric<int>();
  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
 
  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(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 v = getVolume();
  auto ts_a = getTSa();
  auto t_w = getFTensor0IntegrationWeight();
 
  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_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
  auto t_diff_u =
      getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
  auto t_log_u =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
  auto t_log_u2_h1 =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
  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_nb_uniq =
      FTensor::Tensor0<FTensor::PackPtr<int *, 1>>(nbUniq.data().data());
 
  auto no_h1 = [&]() {
    MoFEMFunctionBegin;
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      double a = v * t_w;
      ++t_w;
 
      auto neohookean = [c10](auto v) { return fun_neohookean(c10, v); };
      auto d_neohookean = [c10, lambda](auto v) {
        return fun_d_neohookean(c10, v);
      };
 
      auto t_diff_neohookean = EigenMatrix::getDiffMat(
          t_eigen_vals, t_eigen_vecs, neohookean, d_neohookean, t_nb_uniq);
 
      const auto tr = t_log_u(i, i);
      FTensor::Tensor2<double, size_symm, size_symm> t_dP;
      t_dP(L, J) = -t_L(i, j, L) * ((t_diff_neohookean(i, j, k, l) +
                                     fun_diff_neohookean_bulk(lambda, tr) *
                                         t_kd_sym(i, j) * t_kd_sym(k, l)) *
                                    t_L(k, l, J));
      t_dP(L, J) -= (alpha_u * ts_a) *
                    (t_L(i, j, L) * (t_diff(i, j, k, l) * t_L(k, l, J)));
 
      if constexpr (1) {
        FTensor::Tensor2_symmetric<double, 3> t_deltaP;
        t_deltaP(i, j) = (t_approx_P_adjoint__dstretch(i, j) ||
                          t_approx_P_adjoint__dstretch(j, i)) /
                         2.;
        auto t_diff2_uP = EigenMatrix::getDiffDiffMat(
            t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
            EshelbianCore::dd_f, t_deltaP, t_nb_uniq);
        t_dP(L, J) += t_L(i, j, L) * (t_diff2_uP(i, j, k, l) * t_L(k, l, J));
      }
      ++t_approx_P_adjoint__dstretch;
      ++t_log_u;
      ++t_eigen_vals;
      ++t_eigen_vecs;
      ++t_nb_uniq;
 
      int rr = 0;
      for (; rr != row_nb_dofs / size_symm; ++rr) {
        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, 6 * rr, 0);
        for (int cc = 0; cc != col_nb_dofs / size_symm; ++cc) {
          double b = a * t_row_base_fun * t_col_base_fun;
          double c = (a * alpha_grad_u * ts_a) *
                     (t_row_grad_fun(i) * t_col_grad_fun(i));
          t_m(L, J) -= b * t_dP(L, 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;
      }
 
    }
    MoFEMFunctionReturn(0);
  };
 
  auto large = [&]() {
    MoFEMFunctionBegin;
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Not implemented for Neo-Hookean (used ADOL-C)");
    MoFEMFunctionReturn(0);
  };
 
  switch (EshelbianCore::gradApproximator) {
    case NO_H1_CONFIGURATION:
      CHKERR no_h1();
      break;
    case LARGE_ROT:
    case MODERATE_ROT:
      CHKERR large();
      break;
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "gradApproximator not handled");
    };
 
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity