MatMetaElastic.cpp

Go to the documentation of this file.

/**
 * @file MatMetaElastic.cpp
 * @brief 
 * @version 0.1
 * @date 2026-03-28
 * 
 * @copyright Copyright (c) 2026
 * 
 */
 
namespace MatOps {
 
template <int DIM> struct MatMetaElastic : public MatElasticImpl<DIM> {
  using MatElasticImpl<DIM>::MatElasticImpl;
 
  using A = MatElasticImpl<DIM>;
 
  MoFEMErrorCode getOptions(MoFEM::Interface *m_field_ptr = nullptr) override {
    MoFEMFunctionBegin;
 
    char umat_library[PETSC_MAX_PATH_LEN] = "./umat.so";
    CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR,
                                 "-umat_library", umat_library,
                                 sizeof(umat_library), PETSC_NULLPTR);
 
    if (std::filesystem::exists(std::filesystem::path(umat_library))) {
      metaElasticAtPts.insert(
          {MatOps::MatOpsTagsRegistry::setTagName("Umat"),
           createMatOpsPhysicalEquationsPtr<UMAT, modelType>(
               this->matOpsDataPtr,
               MatOps::MatOpsTagsRegistry::setTagName("Umat"))});
    }
 
    for (auto &m : metaElasticAtPts) {
      CHKERR m.second->getOptions(m_field_ptr);
    }
    A::tagVsRangePtr = boost::make_shared<std::map<int, Range>>();
    for (const auto &m : metaElasticAtPts) {
      A::tagVsRangePtr->insert(
          {m.first, Range()}); // Initialize with empty range
      for (auto &rp : m.second->paramVecByRange) {
        A::tagVsRangePtr->at(m.first).merge(
            rp.first); // Merge ranges for the same tag
      }
    }
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode setParams(FEMethod *fe_ptr, int gg) override {
    MoFEMFunctionBegin;
    const auto ent = fe_ptr->getFEEntityHandle();
 
    bool has_ranges = false;
    if (A::tagVsRangePtr) {
      for (const auto &tag_range : *A::tagVsRangePtr) {
        has_ranges = has_ranges || !tag_range.second.empty();
      }
    }
 
    if (!has_ranges) {
      for (auto &m : metaElasticAtPts) {
        CHKERR m.second->setParams(fe_ptr, gg);
      }
      MoFEMFunctionReturnHot(0);
    }
 
    for (auto &m : metaElasticAtPts) {
      if (A::tagVsRangePtr->at(m.first).find(ent) !=
          A::tagVsRangePtr->at(m.first).end()) {
        CHKERR m.second->setParams(fe_ptr, gg);
        MoFEMFunctionReturnHot(0);
      }
    }
 
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Meta elastic material parameters not found for entity %lu",
            static_cast<unsigned long>(ent));
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode recordTape() override {
    MoFEMFunctionBegin;
    for (auto &m : metaElasticAtPts) {
      CHKERR m.second->recordTape();
    }
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode evaluateVariable(int tag, EntityHandle entity,
                                  int gg) override {
    return metaElasticAtPts.at(tag)->evaluateVariable(tag, entity, gg);
  }
 
  MoFEMErrorCode evaluateDerivatives(int tag, EntityHandle entity,
                                     int gg) override {
    return metaElasticAtPts.at(tag)->evaluateDerivatives(tag, entity, gg);
  }
 
  MoFEMErrorCode updateState(int tag, EntityHandle entity, int gg) override {
    return metaElasticAtPts.at(tag)->updateState(tag, entity, gg);
  }
 
  std::map<int, boost::shared_ptr<PhysicalEquations>> &getMetaElasticMap() {
    return metaElasticAtPts;
  }
 
protected:
  static inline constexpr int modelType =
      (DIM == 3) ? MatOps::MODEL_3D : MatOps::MODEL_2D_PLANE_STRAIN;
 
  std::map<int, boost::shared_ptr<PhysicalEquations>> metaElasticAtPts = {
 
      {
 
          // First tag is always MetaElastic tag
          MatOps::MatOpsTagsRegistry::setTagName("MetaElastic"),
          createMatOpsPhysicalEquationsPtr<ELASTICITY::NEOHOOKEAN, modelType>(
              this->matOpsDataPtr,
              MatOps::MatOpsTagsRegistry::setTagName("MetaElastic"))
 
      },
 
      {
 
          MatOps::MatOpsTagsRegistry::setTagName("MooneyRivlinWriggersEq63"),
          createMatOpsPhysicalEquationsPtr<ELASTICITY::MOONEYRIVLINWRIGGERSEQ63,
                                          modelType>(
              this->matOpsDataPtr, MatOps::MatOpsTagsRegistry::setTagName(
                                      "MooneyRivlinWriggersEq63"))
 
      },
 
      {
 
          MatOps::MatOpsTagsRegistry::setTagName("StVenantKirchhoff"),
          createMatOpsPhysicalEquationsPtr<ELASTICITY::STVENANTKIRCHHOFF,
                                          modelType>(
              this->matOpsDataPtr,
              MatOps::MatOpsTagsRegistry::setTagName("StVenantKirchhoff"))
 
      },
 
  };
};
 
std::map<int, boost::shared_ptr<PhysicalEquations>> &
getMetaElasticMap(boost::shared_ptr<PhysicalEquations> meta_ptr) {
  auto meta_elastic_3d_ptr =
      boost::dynamic_pointer_cast<MatMetaElastic<3>>(meta_ptr);
  if (meta_elastic_3d_ptr) {
    return meta_elastic_3d_ptr->getMetaElasticMap();
  }
  auto meta_elastic_2d_ptr =
      boost::dynamic_pointer_cast<MatMetaElastic<2>>(meta_ptr);
  if (meta_elastic_2d_ptr) {
    return meta_elastic_2d_ptr->getMetaElasticMap();
  }
  CHK_THROW_MESSAGE(
      MOFEM_OPERATION_UNSUCCESSFUL,
      "Provided PhysicalEquations pointer is not of type MatMetaElastic");
  static std::map<int, boost::shared_ptr<PhysicalEquations>> empty_map;
  return empty_map;
}
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<ELASTICITY::META, MODEL_3D>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag) {
  return boost::make_shared<MatMetaElastic<3>>(mat_ops_data_ptr, tag);
}
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<ELASTICITY::META, MODEL_2D_PLANE_STRAIN>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag) {
  return boost::make_shared<MatMetaElastic<2>>(mat_ops_data_ptr, tag);
}
} // namespace MatOps