HookeOps.hpp

Go to the documentation of this file.

/**
 * @file HookeOps.hpp
 * @brief Implementation of Hookes operator Hookes for linear elastic problems
 * in MoFEM.
 * @version 0.14.0
 * @date 2025-04-10
 */
 
#ifndef __HOOKE_OPS_HPP__
#define __HOOKE_OPS_HPP__
 
namespace HookeOps {
// Shared data structure for sotting .
struct CommonData : public boost::enable_shared_from_this<CommonData> {
  boost::shared_ptr<MatrixDouble> matGradPtr;
  boost::shared_ptr<MatrixDouble> matDPtr;
  MatrixDouble matCauchyStressPtr;
  MatrixDouble matStrainPtr;
 
  /**
   * @brief Get shared Cauchy stress
   * @return Shared Cauchy stress matrix
   */
  inline auto getMatCauchyStress() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &matCauchyStressPtr);
  }
  /**
   * @brief Get shared strain
   * @return Shared strain matrix
   */
  inline auto getMatStrain() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matStrainPtr);
  }
};
 
/**
 * @brief // Add operator to compute material stiffness tensor D based on block
 * attributes.
 *
 * Extracts block-wise material properties and initializes elasticity tensor
 * for each finite element block.
 *
 * @tparam DIM Problem dimension
 * @param pip Operator pipeline container
 * @param m_field MoFEM interface reference
 * @param mat_D_Ptr Elasticity tensor pointer
 * @param sev Logging Verbosity level
 * @param scale Young's modulus scale
 */
template <int DIM>
MoFEMErrorCode
addMatBlockOps(MoFEM::Interface &m_field,
               boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
               std::string block_name,
               boost::shared_ptr<MatrixDouble> mat_D_Ptr, Sev sev,
               double scale = 1) {
  MoFEMFunctionBegin;
 
  struct OpMatBlocks : public DomainEleOp {
    OpMatBlocks(boost::shared_ptr<MatrixDouble> m, double bulk_modulus_K,
                double shear_modulus_G, bool is_plane_strain,
                MoFEM::Interface &m_field, Sev sev,
                std::vector<const CubitMeshSets *> meshset_vec_ptr,
                double scale)
        : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), matDPtr(m),
          bulkModulusKDefault(bulk_modulus_K),
          shearModulusGDefault(shear_modulus_G), isPlaneStrain(is_plane_strain),
          scaleYoungModulus(scale) {
      CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
                        "Can not get data from block");
    }
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBegin;
 
      for (auto &b : blockData) {
 
        if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
          CHKERR getMatDPtr(matDPtr, b.bulkModulusK * scaleYoungModulus,
                            b.shearModulusG * scaleYoungModulus, isPlaneStrain);
          MoFEMFunctionReturnHot(0);
        }
      }
 
      CHKERR getMatDPtr(matDPtr, bulkModulusKDefault * scaleYoungModulus,
                        shearModulusGDefault * scaleYoungModulus,
                        isPlaneStrain);
      MoFEMFunctionReturn(0);
    }
 
  private:
    boost::shared_ptr<MatrixDouble> matDPtr;
    const double scaleYoungModulus;
 
    struct BlockData {
      double bulkModulusK;  ///< Bulk modulus
      double shearModulusG; ///< Shear modulus
      Range blockEnts;      ///< Entities in the block
    };
 
    double bulkModulusKDefault;
    double shearModulusGDefault;
    std::vector<BlockData> blockData;
    bool isPlaneStrain; ///< Plane strain flag
 
    MoFEMErrorCode
    extractBlockData(MoFEM::Interface &m_field,
                     std::vector<const CubitMeshSets *> meshset_vec_ptr,
                     Sev sev) {
      MoFEMFunctionBegin;
 
      for (auto m : meshset_vec_ptr) {
        MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *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
          m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
          return ents;
        };
 
        double young_modulus = block_data[0];
        double poisson_ratio = block_data[1];
        double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
        double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
 
        MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
            << "E = " << young_modulus << " nu = " << poisson_ratio;
 
        blockData.push_back(
            {bulk_modulus_K, shear_modulus_G, get_block_ents()});
      }
      MOFEM_LOG_CHANNEL("WORLD");
      MoFEMFunctionReturn(0);
    }
    //! [Calculate elasticity tensor]
    // Compute elasticity tensor D
    MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
                              double bulk_modulus_K, double shear_modulus_G,
                              bool is_plane_strain) {
      MoFEMFunctionBegin;
 
      auto set_material_stiffness = [&]() {
        FTensor::Index<'i', DIM> i;
        FTensor::Index<'j', DIM> j;
        FTensor::Index<'k', DIM> k;
        FTensor::Index<'l', DIM> l;
        constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
        double A = 1.;
        /**
         * @brief Toggle for plane strain.
         * 
         * @note If `DIM == 2` and `is_plane_strain` is true, the code uses the
         *       plane strain elasticity tensor.
         * @note If `DIM == 2` and `is_plane_strain` is false, the code uses the
         *       plane stress elasticity tensor.
         */
        if (DIM == 2 && !is_plane_strain) {
          A = 2 * shear_modulus_G /
              (bulk_modulus_K + (4. / 3.) * shear_modulus_G);
        }
        auto t_D = getFTensor4DdgFromMat<DIM, 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);
      };
 
      constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
      mat_D_ptr->resize(size_symm * size_symm, 1);
      set_material_stiffness();
      MoFEMFunctionReturn(0);
    }
    //! [Calculate elasticity tensor]
  };
 
  double E = 1.0;
  double nu = 0.3;
 
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-young_modulus", &E,
                               PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-poisson_ratio", &nu,
                               PETSC_NULLPTR);
 
  PetscBool is_plane_strain = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-plane_strain",
                             &is_plane_strain, PETSC_NULLPTR);
 
  double bulk_modulus_K = E / (3 * (1 - 2 * nu));
  double shear_modulus_G = E / (2 * (1 + nu));
  pip.push_back(new OpMatBlocks(
      mat_D_Ptr, bulk_modulus_K, shear_modulus_G, is_plane_strain, m_field, sev,
 
      // Get blockset using regular expression
      m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
          (boost::format("%s(.*)") % block_name).str()
 
              )),
      scale
 
      ));
 
  MoFEMFunctionReturn(0);
}
// Factory to initialize common data to get strain, stress and D
template <int DIM, IntegrationType I, typename DomainEleOp>
auto commonDataFactory(
    MoFEM::Interface &m_field,
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
    std::string field_name, std::string block_name, Sev sev, double scale = 1) {
 
  auto common_ptr = boost::make_shared<HookeOps::CommonData>();
  common_ptr->matDPtr = boost::make_shared<MatrixDouble>();
  common_ptr->matGradPtr = boost::make_shared<MatrixDouble>();
 
  CHK_THROW_MESSAGE(addMatBlockOps<DIM>(m_field, pip, block_name,
                                        common_ptr->matDPtr, sev, scale),
                    "addMatBlockOps");
 
  pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(
      field_name, common_ptr->matGradPtr));
  pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(common_ptr->matGradPtr,
                                                  common_ptr->getMatStrain()));
  pip.push_back(new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
      common_ptr->getMatStrain(), common_ptr->getMatCauchyStress(),
      common_ptr->matDPtr));
 
  return common_ptr;
}
 
/**
@brief Factory function to create and push internal force operators for the
domain RHS. (RHS first overload)
 * Assemble domain right-hand-side vector- internal force vector.
 * @param m_field MoFEM interface reference
 * @param field_name Field name string
 * @param common_ptr Hookeian Material data pointer
 * @param is_non_linear option for solving linear nonlinear
 */
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
MoFEMErrorCode opFactoryDomainRhs(
    MoFEM::Interface &m_field,
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
    std::string field_name, boost::shared_ptr<HookeOps::CommonData> common_ptr,
    Sev sev, bool is_non_linear = false) {
  MoFEMFunctionBegin;
  //! [Push Internal forces]
  using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
      A>::template LinearForm<I>;
  using OpInternalForce =
      typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
 
  pip.push_back(
      new OpInternalForce(field_name, common_ptr->getMatCauchyStress(),
                          [is_non_linear](double, double, double) constexpr {
                            return (is_non_linear ? 1 : -1);
                          }));
  //! [Push Internal forces]
  MoFEMFunctionReturn(0);
}
 
/**
 * @brief Overloaded factory for domain RHS for material  scaling  (RHS Second overload)
 *
 * Initializes and pushes operators to assemble the domain RHS vector, 
 * representing internal forces, with optional scaling.
 * @tparam DIM Spatial dimension
 * @tparam A Assembly type
 * @tparam I Integration type
 
 * @param pip Operator container
 * @param field_name Field variable string
 * @param block_name Mesh block name
 * @param sev Logging severity level
 * @param is_non_linear Linear or Nonlinear problem flag
 * @param scale Yonuns modulas scaling factor
 
 */
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
MoFEMErrorCode opFactoryDomainRhs(
    MoFEM::Interface &m_field,
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
    std::string field_name, std::string block_name, Sev sev,
    bool is_non_linear = false, double scale = 1) {
  MoFEMFunctionBegin;
 
  auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
      m_field, pip, field_name, block_name, sev, scale);
  CHKERR opFactoryDomainRhs<DIM, A, I, DomainEleOp>(
      m_field, pip, field_name, common_ptr, sev, is_non_linear);
 
  MoFEMFunctionReturn(0);
}
/**
 * @brief  Assemble domain LHS K factory (LHS first overload)
 * Initializes and pushes operators to assemble the LHS 
 
 */
 
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
MoFEMErrorCode opFactoryDomainLhs(
    MoFEM::Interface &m_field,
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
    std::string field_name, boost::shared_ptr<HookeOps::CommonData> common_ptr,
    Sev sev) {
  MoFEMFunctionBegin;
  //! [OpK]
  using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
      A>::template BiLinearForm<I>;
 
  using OpKCauchy = typename B::template OpGradSymTensorGrad<1, DIM, DIM, 0>;
  //! [OpK]
  // Assumes constant D matrix per entity
  pip.push_back(new OpKCauchy(field_name, field_name, common_ptr->matDPtr));
 
  MoFEMFunctionReturn(0);
}
/**
 * @brief Overloaded LHS factory with optional scaling (LHS Second overload)
 * 
 * @param scale Young’s modulus scale
 */
 
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
MoFEMErrorCode opFactoryDomainLhs(
    MoFEM::Interface &m_field,
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
    std::string field_name, std::string block_name, Sev sev, double scale = 1) {
  MoFEMFunctionBegin;
 
  auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
      m_field, pip, field_name, block_name, sev, scale);
  CHKERR opFactoryDomainLhs<DIM, A, I, DomainEleOp>(m_field, pip, field_name,
                                                    common_ptr, sev);
 
  MoFEMFunctionReturn(0);
}
} // namespace HookeOps
 
#endif // __HOOKE_OPS_HPP__