MatMetaPiolaResponse.cpp

Go to the documentation of this file.

/**
 * @file MatMetaPiolaResponse.cpp
 * @brief Dispatch first-Piola responses without implying stored energy
 * @version 0.1
 * @date 2026-03-28
 *
 * @copyright Copyright (c) 2026
 *
 */
 
namespace MatOps {
 
template <int DIM,
         int MODEL_TYPE =
             (DIM == 3) ? MatOps::MODEL_3D : MatOps::MODEL_2D_PLANE_STRAIN>
struct MatMetaPiolaResponse : public MatPiolaResponseImpl<DIM> {
  using A = MatPiolaResponseImpl<DIM>;
  using ResponsePtr = boost::shared_ptr<MatPiolaResponse>;
 
  MatMetaPiolaResponse(boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag)
      : A(mat_ops_data_ptr, tag) {
    const int default_tag = MatOpsTagsRegistry::setTagName("MetaPiolaResponse");
    metaPiolaResponses.emplace(
        default_tag,
        boost::static_pointer_cast<MatPiolaResponse>(
            createMatOpsPhysicalEquationsPtr<ELASTICITY::NEOHOOKEAN, modelType>(
                this->matOpsDataPtr, default_tag)));
 
    if constexpr (DIM == 3) {
      const int mooney_tag =
          MatOpsTagsRegistry::setTagName("MooneyRivlinWriggersEq63");
      metaPiolaResponses.emplace(
          mooney_tag, boost::static_pointer_cast<MatPiolaResponse>(
                          createMatOpsPhysicalEquationsPtr<
                              ELASTICITY::MOONEYRIVLINWRIGGERSEQ63, MODEL_3D>(
                              this->matOpsDataPtr, mooney_tag)));
    }
 
    const int stvenant_tag =
        MatOpsTagsRegistry::setTagName("StVenantKirchhoff");
    metaPiolaResponses.emplace(
        stvenant_tag,
        boost::static_pointer_cast<MatPiolaResponse>(
            createMatOpsPhysicalEquationsPtr<ELASTICITY::STVENANTKIRCHHOFF,
                                             modelType>(this->matOpsDataPtr,
                                                        stvenant_tag)));
  }
 
  MoFEMErrorCode getOptions(MoFEM::Interface *m_field_ptr = nullptr) override {
    MoFEMFunctionBegin;
 
    // Activate UMAT only for an explicit library or material block.
    char umat_library[PETSC_MAX_PATH_LEN] = "";
    CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-umat_library",
                                 umat_library, sizeof(umat_library),
                                 PETSC_NULLPTR);
 
    if (umat_library[0] ||
        (m_field_ptr &&
         !m_field_ptr->getInterface<MeshsetsManager>()
              ->getCubitMeshsetPtr(std::regex("MAT_UMAT(.*)"))
              .empty())) {
      const int umat_tag = MatOpsTagsRegistry::setTagName("Umat");
      metaPiolaResponses.insert(
          {umat_tag, boost::static_pointer_cast<MatPiolaResponse>(
                         createMatOpsPhysicalEquationsPtr<UMAT, modelType>(
                             this->matOpsDataPtr, umat_tag))});
    }
 
    for (auto &m : metaPiolaResponses) {
      CHKERR m.second->getOptions(m_field_ptr);
    }
    A::tagVsRangePtr = boost::make_shared<std::map<int, Range>>();
    for (const auto &m : metaPiolaResponses) {
      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 : metaPiolaResponses) {
        CHKERR m.second->setParams(fe_ptr, gg);
      }
      MoFEMFunctionReturnHot(0);
    }
 
    for (auto &m : metaPiolaResponses) {
      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 Piola-response parameters not found for entity %lu",
            static_cast<unsigned long>(ent));
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode recordTape() override {
    MoFEMFunctionBegin;
    for (auto &m : metaPiolaResponses) {
      CHKERR m.second->recordTape();
    }
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode evaluateVariable(int tag, EntityHandle entity,
                                  int gg) override {
    return metaPiolaResponses.at(tag)->evaluateVariable(tag, entity, gg);
  }
 
  MoFEMErrorCode evaluateDerivatives(int tag, EntityHandle entity,
                                     int gg) override {
    return metaPiolaResponses.at(tag)->evaluateDerivatives(tag, entity, gg);
  }
 
  MoFEMErrorCode updateState(int tag, EntityHandle entity, int gg) override {
    return metaPiolaResponses.at(tag)->updateState(tag, entity, gg);
  }
 
  std::map<int, ResponsePtr> &getMetaPiolaResponseMap() {
    return metaPiolaResponses;
  }
 
protected:
  static inline constexpr int modelType = MODEL_TYPE;
 
  std::map<int, ResponsePtr> metaPiolaResponses;
};
 
std::map<int, boost::shared_ptr<MatPiolaResponse>> &
getMetaPiolaResponseMap(boost::shared_ptr<MatPiolaResponse> meta_ptr) {
  auto meta_response_3d_ptr =
      boost::dynamic_pointer_cast<MatMetaPiolaResponse<3>>(meta_ptr);
  if (meta_response_3d_ptr) {
    return meta_response_3d_ptr->getMetaPiolaResponseMap();
  }
  auto meta_response_2d_ptr =
      boost::dynamic_pointer_cast<MatMetaPiolaResponse<2>>(meta_ptr);
  if (meta_response_2d_ptr) {
    return meta_response_2d_ptr->getMetaPiolaResponseMap();
  }
  auto meta_response_axisymmetric_ptr = boost::dynamic_pointer_cast<
      MatMetaPiolaResponse<3, MODEL_AXISYMMETRIC>>(meta_ptr);
  if (meta_response_axisymmetric_ptr) {
    return meta_response_axisymmetric_ptr->getMetaPiolaResponseMap();
  }
  CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                    "Provided response pointer is not a MatMetaPiolaResponse");
  static std::map<int, boost::shared_ptr<MatPiolaResponse>> empty_map;
  return empty_map;
}
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<PIOLA_RESPONSE::META, MODEL_3D>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag) {
  return boost::make_shared<MatMetaPiolaResponse<3>>(mat_ops_data_ptr, tag);
}
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<PIOLA_RESPONSE::META, MODEL_2D_PLANE_STRAIN>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag) {
  return boost::make_shared<MatMetaPiolaResponse<2>>(mat_ops_data_ptr, tag);
}
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<PIOLA_RESPONSE::META, MODEL_AXISYMMETRIC>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag) {
  return boost::make_shared<MatMetaPiolaResponse<3, MODEL_AXISYMMETRIC>>(
      mat_ops_data_ptr, tag);
}
} // namespace MatOps