MatGenericElasticSimpleDamage.hpp

Go to the documentation of this file.

/**
 * @file MatGenericElasticSimpleDamage.hpp
 * @version 0.1
 * @date 2026-04-19
 * 
 * @copyright Copyright (c) 2026
 * 
 */
 
/*
  @note: This is only example. This model is not real damage model, and do not
  have any physical meaning. It is only for demonstration of how to implement 
  a simple damage model and update state variables. 
*/
 
#ifndef ADOLC_GENERIC_ELASTIC_SIMPLE_DAMAGE_HPP
#define ADOLC_GENERIC_ELASTIC_SIMPLE_DAMAGE_HPP
 
#include <algorithm>
#include <cmath>
 
namespace GenericElasticExamples {
 
inline auto makeSimpleDamageStiffness(double bulk_modulus_K,
                                      double shear_modulus_G) {
  constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
  auto mat_D_ptr = boost::make_shared<MatrixDouble>(1, size_symm * size_symm);
 
  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>();
 
  double A = 1.;
  if (SPACE_DIM == 2) {
    A = 2 * shear_modulus_G /
        (bulk_modulus_K + (4. / 3.) * shear_modulus_G);
  }
 
  auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
  t_D(i, j, k, l) = 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
                    A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
                        t_kd(i, j) * t_kd(k, l);
 
  return mat_D_ptr;
}
 
inline MoFEMErrorCode getSimpleDamageModuliFromBlock(
    MoFEM::Interface &m_field, double &bulk_modulus_K,
    double &shear_modulus_G) {
  MoFEMFunctionBegin;
 
  const auto blocks =
      m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
          std::regex("MAT_ELASTIC.*"));
 
  if (blocks.empty()) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "No MAT_ELASTIC block found for simple damage example");
  }
  // Note: we expect only one block with elastic properties, but we read the first
  std::vector<double> block_data;
  CHKERR blocks.front()->getAttributes(block_data);
  if (block_data.size() != 2) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Expected MAT_ELASTIC block with two attributes");
  }
 
  const double young_modulus = block_data[0];
  const double poisson_ratio = block_data[1];
  bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
  shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
 
  MoFEMFunctionReturn(0);
}
 
inline MoFEMErrorCode configureSimpleDamage(
    MoFEM::Interface &m_field,
    boost::shared_ptr<MatOps::GenericElastic> physical_ptr) {
  MoFEMFunctionBegin;
 
  double max_damage = 0.999;
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-generic_max_damage",
                               &max_damage, PETSC_NULLPTR);
  max_damage = std::min(0.999, std::max(0., max_damage));
 
  double bulk_modulus_K = 0;
  double shear_modulus_G = 0;
  CHKERR getSimpleDamageModuliFromBlock(m_field, bulk_modulus_K,
                                        shear_modulus_G);
  auto mat_D_ptr =
      makeSimpleDamageStiffness(bulk_modulus_K, shear_modulus_G);
 
  CHKERR physical_ptr->matOpsDataPtr->bindStateTag(m_field, "damage", 1);
  CHKERR physical_ptr->matOpsDataPtr->bindStateTag(m_field, "F",
                                                  SPACE_DIM * SPACE_DIM);
  CHKERR physical_ptr->matOpsDataPtr->setupStateData();
 
  physical_ptr->hookEvaluateVariable =
      [mat_D_ptr](boost::shared_ptr<MatOps::MatOpsData> mat_ops_data_ptr,
                  int, EntityHandle ent, int gg) {
        MoFEMFunctionBegin;
 
        auto t_grad = getFTensor2FromPtr<SPACE_DIM, SPACE_DIM>(
            mat_ops_data_ptr->getActiveDataPtr("F")->data().data());
        auto t_P = getFTensor2FromPtr<SPACE_DIM, SPACE_DIM>(
            mat_ops_data_ptr->getDependentDataPtr("P")->data().data());
        auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
        auto t_damage = getFTensor0FromPtr(
            &(*mat_ops_data_ptr->getStateDataPtr("damage", ent, gg))(0, 0));
 
        FTENSOR_INDEX(SPACE_DIM, i);
        FTENSOR_INDEX(SPACE_DIM, J);
        FTENSOR_INDEX(SPACE_DIM, k);
        FTENSOR_INDEX(SPACE_DIM, L);
 
        t_P(i, J) = (1. - t_damage) * t_D(i, J, k, L) * t_grad(k, L);
 
        MoFEMFunctionReturn(0);
      };
 
  physical_ptr->hookEvaluateDerivatives =
      [mat_D_ptr](boost::shared_ptr<MatOps::MatOpsData> mat_ops_data_ptr,
                  int, EntityHandle ent, int gg) {
        MoFEMFunctionBegin;
 
        auto t_dP = getFTensor4FromPtr<SPACE_DIM, SPACE_DIM, SPACE_DIM,
                                       SPACE_DIM>(mat_ops_data_ptr
                                                      ->getDependentDerivativesDataPtr(
                                                          "P_dF")
                                                      ->data()
                                                      .data());
        auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
        auto t_damage = getFTensor0FromPtr(
            &(*mat_ops_data_ptr->getStateDataPtr("damage", ent, gg))(0, 0));
 
        FTENSOR_INDEX(SPACE_DIM, i);
        FTENSOR_INDEX(SPACE_DIM, J);
        FTENSOR_INDEX(SPACE_DIM, k);
        FTENSOR_INDEX(SPACE_DIM, L);
 
        t_dP(i, J, k, L) = (1. - t_damage) * t_D(i, J, k, L);
 
        MoFEMFunctionReturn(0);
      };
 
  physical_ptr->hookUpdateState =
      [max_damage](boost::shared_ptr<MatOps::MatOpsData> mat_ops_data_ptr,
                   int, EntityHandle ent, int gg) {
        MoFEMFunctionBegin;
 
        auto t_F = getFTensor2FromPtr<SPACE_DIM, SPACE_DIM>(
            mat_ops_data_ptr->getActiveDataPtr("F")->data().data());
        // state variables damage and F (for testing only)
        auto damage_ptr = mat_ops_data_ptr->getStateDataPtr("damage", ent, gg);
        auto t_damage = getFTensor0FromPtr(&(*damage_ptr)(0, 0));
        auto F_state_ptr = mat_ops_data_ptr->getStateDataPtr("F", ent, gg);
        auto t_F_state = getFTensor2FromPtr<SPACE_DIM, SPACE_DIM>(
            &(*F_state_ptr)(0, 0));
 
        FTENSOR_INDEX(SPACE_DIM, i);
        FTENSOR_INDEX(SPACE_DIM, J);
 
        const double grad_norm = std::sqrt(t_F(i, J) * t_F(i, J));
        const double old_damage = t_damage;
        
        t_damage = std::min(max_damage, std::max(old_damage, grad_norm));
        t_F_state(i, J) = t_F(i, J);
 
        MoFEMFunctionReturn(0);
      };
 
  MoFEMFunctionReturn(0);
}
 
} // namespace GenericElasticExamples
 
#endif // ADOLC_GENERIC_ELASTIC_SIMPLE_DAMAGE_HPP