MatOpsExampleGeneric.hpp

Go to the documentation of this file.

/**
 * @file MatOpsExampleGeneric.hpp
 * @author your name (you@domain.com)
 * @brief 
 * @version 0.1
 * @date 2026-04-12
 * 
 * @copyright Copyright (c) 2026
 * 
 */
 
struct MatOpsExampleGeneric {
  MatOpsExampleGeneric() {
    constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
    matDPtr = boost::make_shared<MatrixDouble>(1, size_symm * size_symm);
    auto set_material_stiffness = [this]() {
      MoFEMFunctionBegin;
      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.;
      /**
         * @brief Toggle for plane strain.
         * 
         * @note If `SPACE_DIM == 2` and `is_plane_strain` is true, the code uses the
         *       plane strain elasticity tensor.
         * @note If `SPACE_DIM == 2` and `is_plane_strain` is false, the code uses the
         *       plane stress elasticity tensor.
         */
 
      constexpr double bulk_modulus_K = 10;
      constexpr double shear_modulus_G = 1;
      if (SPACE_DIM == 2 /*&& !is_plane_strain*/) {
        A = 2 * shear_modulus_G /
            (bulk_modulus_K + (4. / 3.) * shear_modulus_G);
      }
      auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*matDPtr);
      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);
      MoFEMFunctionReturn(0);
    };
    CHK_THROW_MESSAGE(set_material_stiffness(),
                      "Failed to set material stiffness");
  }
 
  auto getMaterialStiffness() { return matDPtr; }
 
  static MoFEMErrorCode
  hookEvaluateVariable(boost::shared_ptr<MatOps::MatOpsData> mat_ops_data_ptr,
                       int, EntityHandle, int) {
    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());
    FTENSOR_INDEX(SPACE_DIM, i);
    FTENSOR_INDEX(SPACE_DIM, J);
    FTENSOR_INDEX(SPACE_DIM, k);
    FTENSOR_INDEX(SPACE_DIM, L);
    auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*matDPtr);
    t_P(i, J) = t_D(i, J, k, L) * t_grad(k, L);
    return 0;
  }
 
  static MoFEMErrorCode
  hookEvaluateDerivatives(boost::shared_ptr<MatOps::MatOpsData> mat_ops_data_ptr,
                          int, EntityHandle, int) {
    auto t_dP = getFTensor4FromPtr<SPACE_DIM, SPACE_DIM, SPACE_DIM, SPACE_DIM>(
        mat_ops_data_ptr->getDependentDerivativesDataPtr("P_dF")->data().data());
    FTENSOR_INDEX(SPACE_DIM, i);
    FTENSOR_INDEX(SPACE_DIM, J);
    FTENSOR_INDEX(SPACE_DIM, k);
    FTENSOR_INDEX(SPACE_DIM, L);
    auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*matDPtr);
    t_dP(i, J, k, L) = t_D(i, J, k, L);
    return 0;
  };
 
  static MoFEMErrorCode hookUpdateState(boost::shared_ptr<MatOps::MatOpsData>,
                                        int, EntityHandle, int) {
    return 0;
  };
 
protected:
  inline static boost::shared_ptr<MatrixDouble> matDPtr;
};