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_OFF_INVERSE
// #define NEOHOOKEAN_SCALING
struct HMHNeohookean : public PhysicalEquations {
 
  static constexpr int numberOfActiveVariables = 0;
  static constexpr int numberOfDependentVariables = 0;
 
  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 (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);
  }
 
  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; }
  };
 
  virtual EshelbianPlasticity::OpJacobian *
  returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs,
                   boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                   boost::shared_ptr<PhysicalEquations> physics_ptr) {
    return (new OpJacobian(tag, eval_rhs, eval_lhs, data_ptr, physics_ptr));
  }
 
  MoFEMErrorCode getOptions() {
    MoFEMFunctionBegin;
    CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "neo_hookean_", "", "none");
 
    c10_default = 1;
    CHKERR PetscOptionsScalar("-c10", "C10", "", c10_default, &c10_default,
                              PETSC_NULL);
    K_default = 1;
    CHKERR PetscOptionsScalar("-K", "Bulk modulus K", "", K_default, &K_default,
                              PETSC_NULL);
 
    alphaGradU = 0;
    CHKERR PetscOptionsScalar("-viscosity_alpha_grad_u", "viscosity", "",
                              alphaGradU, &alphaGradU, PETSC_NULL);
 
    ierr = PetscOptionsEnd();
    CHKERRG(ierr);
 
    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;
    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 two attribute");
      }
      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;
    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) {
    return new OpSpatialPhysical(field_name, data_ptr, alpha_u);
  }
 
  struct OpSpatialPhysical_du_du : public OpAssembleVolume {
    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) {
    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.;
 
};
 
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) {}
 
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_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;
 
      FTensor::Dg<double, SPACE_DIM, size_symm> t_Ldiff_u;
      t_Ldiff_u(i, j, L) = t_diff_u(i, j, k, l) * t_L(k, l, L);
      ++t_diff_u;
 
      FTensor::Tensor2_symmetric<double, SPACE_DIM> t_Sigma_u;
      t_Sigma_u(i, j) = 2.0 * c10 * t_u(i, j);
      const double tr = t_log_u(i, i);
      const double J = EshelbianCore::f(tr);
      const double Sigma_J = -2.0 * c10 + K * (J * J - J);
 
      FTensor::Tensor2_symmetric<double, SPACE_DIM> t_viscous_P;
      t_viscous_P(i, j) =
          alpha_u * (t_dot_log_u(i, j) -
                     t_dot_log_u(k, l) * t_kd_sym(k, l) * t_kd_sym(i, j));
 
      FTensor::Tensor1<double, size_symm> t_residual;
      t_residual(L) = t_approx_P_adjoint_log_du(L);
      t_residual(L) -= t_Ldiff_u(i, j, L) * t_Sigma_u(i, j);
      t_residual(L) -= (t_L(i, j, L) * t_kd(i, j)) * Sigma_J;
      t_residual(L) -= t_L(i, j, L) * t_viscous_P(i, j);
      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_u;
      ++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;
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      double a = v * t_w;
      ++t_w;
  
      FTensor::Tensor3<double, SPACE_DIM, SPACE_DIM, size_symm> t_Ldiff_u;
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_h1;
      t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
      t_Ldiff_u(i, j, L) = (t_diff_u(i, m, k, l) * t_h1(m, j)) * t_L(k, l, L);
      ++t_diff_u;
      ++t_grad_h1;
  
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_Sigma_u;
      const double tr_h1 = t_log_u2_h1(i, j) * t_kd_sym(i, j) / 2;
      const double tr_u = t_log_u(i, j) * t_kd_sym(i, j);
      const double tr = tr_h1 + tr_u;
      const double J = EshelbianCore::f(tr);
 
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_total_u;
      t_total_u(i, j) = t_u(i, m) * t_h1(m, j);
#ifndef NEOHOOKEAN_OFF_INVERSE
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_inv_total_u;
      invertTensor3by3(t_total_u, J, t_inv_total_u);
#endif // NEOHOOKEAN_OFF_INVERSE
 
#ifndef NEOHOOKEAN_OFF_INVERSE
      if (tr > 0) {
        t_Sigma_u(i, j) =
            2.0 * c10 * t_total_u(i, m) *
            (t_kd_sym(m, j) - t_inv_total_u(m, k) * t_inv_total_u(k, j));
      } else {
        t_Sigma_u(i, j) = -2.0 * c10 * t_inv_total_u(i, m) *
                          (t_kd_sym(m, j) - t_total_u(m, k) * t_total_u(k, j));
      }
#else
      t_Sigma_u(i, j) = 2.0 * c10 * t_total_u(i, j);
#endif // NEOHOOKEAN_OFF_INVERSE
 
#ifndef NEOHOOKEAN_OFF_INVERSE
      const double Sigma_J = K * (J * J - J);
#else
      const double Sigma_J = -2.0 * c10 + K * (J * J - J);
#endif // NEOHOOKEAN_OFF_INVERSE
 
      FTensor::Tensor2_symmetric<double, SPACE_DIM> t_viscous_P;
      t_viscous_P(i, j) =
          alpha_u * (t_dot_log_u(i, j) -
                     t_dot_log_u(k, l) * t_kd_sym(k, l) * t_kd_sym(i, j));
  
      FTensor::Tensor1<double, size_symm> t_residual;
      t_residual(L) = t_approx_P_adjoint_log_du(L);
      t_residual(L) -= t_Ldiff_u(i, j, L) * t_Sigma_u(i, j);
      t_residual(L) -= (t_L(i, j, L) * t_kd(i, j)) * Sigma_J;
      t_residual(L) -= t_L(i, j, L) * t_viscous_P(i, j);
      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_u;
      ++t_log_u;
      ++t_log_u2_h1;
      ++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);
  };
 
  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::OpSpatialPhysical_du_du::OpSpatialPhysical_du_du(
    std::string row_field, std::string col_field,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha)
    : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false),
      alphaU(alpha) {
  sYmm = false;
}
 
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_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
  auto r_P_du = getFTensor4FromMat<3, 3, 3, 3>(dataAtPts->P_du);
 
  auto no_h1 = [&]() {
    MoFEMFunctionBegin;
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      double a = v * t_w;
      ++t_w;
 
      FTensor::Dg<double, SPACE_DIM, size_symm> t_Ldiff_u;
      t_Ldiff_u(i, j, L) = t_diff_u(i, j, k, l) * t_L(k, l, L);
      ++t_diff_u;
 
      const double tr = t_log_u(i, i);
      const double det_u = EshelbianCore::f(tr);
      ++t_log_u;
 
      FTensor::Tensor1<double, size_symm> Sigma_J_dtr;
      Sigma_J_dtr(L) =
          (lambda * (2 * det_u * det_u - det_u)) * (t_kd(i, j) * t_L(i, j, L));
 
 
      FTensor::Tensor2<double, size_symm, size_symm> t_dP;
      t_dP(L, J) = (-2.0 * c10) * (t_Ldiff_u(i, j, L) * t_Ldiff_u(i, j, J));
      t_dP(L, J) -= (t_L(i, j, L) * t_kd(i, j)) * Sigma_J_dtr(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) -
                                     (t_diff(m, n, k, l) * t_kd_sym(m, n)) *
                                         t_kd_sym(i, j) * t_L(k, l, J)));
 
      FTensor::Tensor2_symmetric<double, SPACE_DIM> t_Sigma_u;
      t_Sigma_u(i, j) = 2.0 * c10 * t_u(i, j);
      ++t_u;
 
      if constexpr (1) {
        FTensor::Tensor2<double, 3, 3> t_deltaP;
        t_deltaP(i, j) = t_approx_P_adjoint_dstretch(i, j) - t_Sigma_u(i, j);
        auto t_diff2_uP = EigenMatrix::getDiffDiffMat(
            t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
            EshelbianCore::dd_f, t_deltaP, nbUniq[gg]);
        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_eigen_vals;
      ++t_eigen_vecs;
 
      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;
      }
 
      ++t_P;
      ++r_P_du;
    }
    MoFEMFunctionReturn(0);
  };
 
  auto large = [&]() {
    MoFEMFunctionBegin;
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      double a = v * t_w;
      ++t_w;
 
      FTensor::Tensor3<double, SPACE_DIM, SPACE_DIM, size_symm> t_Ldiff_u;
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_h1;
 
      switch (EshelbianCore::gradApproximator) {
      case NO_H1_CONFIGURATION:
        t_h1(i, j) = t_kd(i, j);
        t_Ldiff_u(i, j, L) = t_diff_u(i, j, k, l) * t_L(k, l, L);
        break;
      case LARGE_ROT:
      case MODERATE_ROT:
        t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
        t_Ldiff_u(i, j, L) = (t_diff_u(i, m, k, l) * t_h1(m, j)) * t_L(k, l, L);
        break;
      default:
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "gradApproximator not handled");
      };
      ++t_diff_u;
      ++t_grad_h1;
 
      const double tr_h1 = t_log_u2_h1(i, j) * t_kd_sym(i, j) / 2;
      const double tr_u = t_log_u(i, j) * t_kd_sym(i, j);
      const double tr = tr_h1 + tr_u;
      const double det_u = EshelbianCore::f(tr);
      ++t_log_u;
      ++t_log_u2_h1;
 
      FTensor::Tensor1<double, size_symm> Sigma_J_dtr;
      Sigma_J_dtr(L) =
          (lambda * (2 * det_u * det_u - det_u)) * (t_kd(i, j) * t_L(i, j, L));
 
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_total_u;
      t_total_u(i, j) = t_u(i, m) * t_h1(m, j);
#ifndef NEOHOOKEAN_OFF_INVERSE
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_inv_total_u;
      invertTensor3by3(t_total_u, det_u, t_inv_total_u);
#endif // NEOHOOKEAN_OFF_INVERSE
      ++t_u;
 
      FTensor::Tensor2<double, size_symm, size_symm> t_dP;
#ifndef NEOHOOKEAN_OFF_INVERSE
      FTensor::Tensor3<double, SPACE_DIM, SPACE_DIM, size_symm> t_Sigma_u_du;
      t_Sigma_u_du(i, j, J) = t_Ldiff_u(i, j, J) + t_inv_total_u(i, k) *
                                                       t_Ldiff_u(k, l, J) *
                                                       t_inv_total_u(l, j);
      t_dP(L, J) = (-2.0 * c10) * (t_Ldiff_u(i, j, L) * t_Sigma_u_du(i, j, J));
#else
      t_dP(L, J) = (-2.0 * c10) * (t_Ldiff_u(i, j, L) * t_Ldiff_u(i, j, J));
#endif // NEOHOOKEAN_OFF_INVERSE
 
      t_dP(L, J) -= (t_L(i, j, L) * t_kd(i, j)) * Sigma_J_dtr(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) -
                                     (t_diff(m, n, k, l) * t_kd_sym(m, n)) *
                                         t_kd_sym(i, j) * t_L(k, l, J)));
 
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_Sigma_u;
#ifndef NEOHOOKEAN_OFF_INVERSE
      if (tr > 0) {
        t_Sigma_u(i, j) =
            2.0 * c10 * t_total_u(i, m) *
            (t_kd_sym(m, j) - t_inv_total_u(m, k) * t_inv_total_u(k, j));
      } else {
        t_Sigma_u(i, j) = -2.0 * c10 * t_inv_total_u(i, m) *
                          (t_kd_sym(m, j) - t_total_u(m, k) * t_total_u(k, j));
      }
#else
      t_Sigma_u(i, j) = 2.0 * c10 * t_total_u(i, j);
#endif // NEOHOOKEAN_OFF_INVERSE
 
      FTensor::Tensor2<double, 3, 3> t_deltaP;
      switch (EshelbianCore::gradApproximator) {
      case NO_H1_CONFIGURATION:
      case LARGE_ROT:
        t_deltaP(i, j) =
            t_approx_P_adjoint_dstretch(i, j) - t_h1(j, n) * t_Sigma_u(i, n);
        break;
      case MODERATE_ROT:
        t_deltaP(i, j) = -t_h1(j, n) * t_Sigma_u(i, n);
        break;
      default:
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "gradApproximator not handled");
      };
      ++t_approx_P_adjoint_dstretch;
 
      if constexpr (1) {
        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 t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
            t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
            EshelbianCore::dd_f, t_deltaP_sym, nbUniq[gg]);
        t_dP(L, J) += t_L(i, j, L) * (t_diff2_uP2(i, j, k, l) * t_L(k, l, J));
      }
      ++t_eigen_vals;
      ++t_eigen_vecs;
 
      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;
      }
 
      ++t_P;
      ++r_P_du;
    }
    MoFEMFunctionReturn(0);
  };
 
  switch (EshelbianCore::gradApproximator) {
    case NO_H1_CONFIGURATION:
      CHKERR large();
      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