MatElastic.hpp

Go to the documentation of this file.

/**
 * @file MatElastic.hpp
 * @brief 
 * @version 0.1
 * @date 2026-04-02
 * 
 * @copyright Copyright (c) 2026
 * 
 */
 
#ifndef MAT_ELASTIC_HPP
#define MAT_ELASTIC_HPP
 
#include <MatUmat.hpp>
#include <MatAxisymmetric.hpp>
 
namespace MatOps {
 
/**
 * @brief Conservative elastic response with a stored-energy density.
 *
 * Every concrete derived model records the scalar dependent datum "energy"
 * consistently with P = dW/dF.
 */
struct MatElastic : public MatPiolaResponse {
  using MatPiolaResponse::MatPiolaResponse;
 
  static constexpr auto energyDensityDataName = "energy";
 
  MatOpsData::MatrixPtr getEnergyDensityPtr() const {
    return matOpsDataPtr->getDependentDataPtr(energyDensityDataName);
  }
};
 
struct PIOLA_RESPONSE {
  PIOLA_RESPONSE() = delete;
  struct META {};
  struct VOLUMELENGTHQUALITY {};
};
 
struct ELASTICITY {
  ELASTICITY() = delete;
  // Behaviors
  struct NEOHOOKEAN {};
  struct MOONEYRIVLINWRIGGERSEQ63 {};
  struct STVENANTKIRCHHOFF {};
};
 
/** \name First-Piola Response Operator Factory */
/**@{*/
template <int MODEL_TYPE> struct OpMaterialFactory<PIOLA_RESPONSE, MODEL_TYPE> {
  OpMaterialFactory() = delete;
 
  static void
  addKinematics(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
                const std::string &field_name,
                boost::shared_ptr<MatrixDouble> m_grad) {
    constexpr int DIM = (MODEL_TYPE == MODEL_3D) ? 3 : 2;
    if constexpr (MODEL_TYPE == MODEL_AXISYMMETRIC) {
      auto m_disp = boost::make_shared<MatrixDouble>();
      auto m_grad_2d = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateVectorFieldValues<2>(field_name, m_disp));
      pip.push_back(
          new OpCalculateVectorFieldGradient<2, 2>(field_name, m_grad_2d));
      pip.push_back(
          new OpCalculateAxisymmetricGradient(m_disp, m_grad_2d, m_grad));
    } else {
      pip.push_back(
          new OpCalculateVectorFieldGradient<DIM, DIM>(field_name, m_grad));
    }
  }
 
  template <AssemblyType A, IntegrationType I, typename DomainEleOp>
  static MoFEMErrorCode
  opRhsFactory(MoFEM::Interface &m_field,
               boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
               std::string field_name,
               boost::shared_ptr<PhysicalEquations> physical_equations_ptr) {
    MoFEMFunctionBegin;
    constexpr int DIM = (MODEL_TYPE == MODEL_3D) ? 3 : 2;
    using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
        A>::template LinearForm<I>;
    using OpInternalForce = typename B::template OpGradTimesTensor<1, DIM, DIM>;
    auto m_grad =
        physical_equations_ptr->matOpsDataPtr->getCommonDataPtr("grad");
    auto m_P = physical_equations_ptr->matOpsDataPtr->getCommonDataPtr("P");
    addKinematics(pip, field_name, m_grad);
    pip.push_back(physical_equations_ptr->createOp(physical_equations_ptr, true,
                                                   false, false));
    if constexpr (MODEL_TYPE == MODEL_AXISYMMETRIC) {
      auto data = boost::make_shared<AxisymmetricAssemblyData>();
      auto jacobian = [](double r, double, double) { return 2 * M_PI * r; };
      pip.push_back(new OpCalculateAxisymmetricStress(m_P, data));
      pip.push_back(
          new OpInternalForce(field_name, data->inPlaneStress, jacobian));
      pip.push_back(new OpAxisymmetricRhs<A, DomainEleOp>(field_name, m_P));
    } else {
      pip.push_back(new OpInternalForce(field_name, m_P));
    }
    MoFEMFunctionReturn(0);
  }
 
  template <AssemblyType A, IntegrationType I, typename DomainEleOp>
  static MoFEMErrorCode
  opLhsFactory(MoFEM::Interface &m_field,
               boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
               std::string field_name,
               boost::shared_ptr<PhysicalEquations> physical_equations_ptr) {
    MoFEMFunctionBegin;
    constexpr int DIM = (MODEL_TYPE == MODEL_3D) ? 3 : 2;
    using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
        A>::template BiLinearForm<I>;
    using OpKPiola = typename B::template OpGradTensorGrad<1, DIM, DIM, -1>;
    auto m_grad =
        physical_equations_ptr->matOpsDataPtr->getCommonDataPtr("grad");
    auto m_P_dF =
        physical_equations_ptr->matOpsDataPtr->getCommonDataPtr("P_dF");
    addKinematics(pip, field_name, m_grad);
    pip.push_back(physical_equations_ptr->createOp(physical_equations_ptr,
                                                   false, true, false));
    if constexpr (MODEL_TYPE == MODEL_AXISYMMETRIC) {
      auto data = boost::make_shared<AxisymmetricAssemblyData>();
      auto jacobian = [](double r, double, double) { return 2 * M_PI * r; };
      pip.push_back(new OpCalculateAxisymmetricTangent(m_P_dF, data));
      pip.push_back(new OpKPiola(field_name, field_name, data->inPlaneTangent,
                                jacobian));
      pip.push_back(
          new OpAxisymmetricLhs<A, DomainEleOp>(field_name, m_P_dF));
    } else {
      pip.push_back(new OpKPiola(field_name, field_name, m_P_dF));
    }
    MoFEMFunctionReturn(0);
  }
 
  static MoFEMErrorCode
  opUpdateFactory(MoFEM::Interface &m_field,
                  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
                  std::string field_name,
                  boost::shared_ptr<PhysicalEquations> physical_equations_ptr) {
    MoFEMFunctionBegin;
    constexpr int DIM = (MODEL_TYPE == MODEL_3D) ? 3 : 2;
    auto m_grad =
        physical_equations_ptr->matOpsDataPtr->getCommonDataPtr("grad");
    auto m_P = physical_equations_ptr->matOpsDataPtr->getCommonDataPtr("P");
    addKinematics(pip, field_name, m_grad);
    pip.push_back(physical_equations_ptr->createOp(physical_equations_ptr,
                                                   false, false, true));
    MoFEMFunctionReturn(0);
  }
 
  static MoFEMErrorCode opPostProcFactory(
      MoFEM::Interface &m_field,
      boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
      std::string field_name,
      boost::shared_ptr<PhysicalEquations> physical_equations_ptr) {
    MoFEMFunctionBegin;
    constexpr int DIM = (MODEL_TYPE == MODEL_3D) ? 3 : 2;
    auto m_grad =
        physical_equations_ptr->matOpsDataPtr->getCommonDataPtr("grad");
    auto m_P = physical_equations_ptr->matOpsDataPtr->getCommonDataPtr("P");
    addKinematics(pip, field_name, m_grad);
    pip.push_back(physical_equations_ptr->createOp(physical_equations_ptr, true,
                                                   false, false));
    MoFEMFunctionReturn(0);
  }
};
/**@}*/
 
/** \name First-Piola Response and Elastic Model Creators */
/**@{*/
 
std::map<int, boost::shared_ptr<MatPiolaResponse>> &
getMetaPiolaResponseMap(boost::shared_ptr<MatPiolaResponse> meta_ptr);
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<PIOLA_RESPONSE::META, MODEL_3D>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag);
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<PIOLA_RESPONSE::META, MODEL_2D_PLANE_STRAIN>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag);
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<PIOLA_RESPONSE::META, MODEL_AXISYMMETRIC>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag);
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<ELASTICITY::NEOHOOKEAN, MODEL_3D>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag);
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<ELASTICITY::NEOHOOKEAN, MODEL_2D_PLANE_STRAIN>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag);
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<ELASTICITY::NEOHOOKEAN, MODEL_AXISYMMETRIC>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag);
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<ELASTICITY::MOONEYRIVLINWRIGGERSEQ63,
                                 MODEL_3D>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag);
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<ELASTICITY::STVENANTKIRCHHOFF, MODEL_3D>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag);
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<ELASTICITY::STVENANTKIRCHHOFF,
                                 MODEL_2D_PLANE_STRAIN>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag);
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<ELASTICITY::STVENANTKIRCHHOFF,
                                 MODEL_AXISYMMETRIC>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag);
 
enum VolumeLengthQualityType {
  QUALITY,
  BARRIER_AND_QUALITY,
  BARRIER_AND_CHANGE_QUALITY,
  BARRIER_AND_CHANGE_QUALITY_SCALED_BY_VOLUME,
  LASTOP_VOLUMELENGTHQUALITYTYPE
};
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<PIOLA_RESPONSE::VOLUMELENGTHQUALITY, MODEL_3D>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag);
 
/**@}*/
 
} // namespace MatOps
 
#endif // MAT_ELASTIC_HPP