MatElastic.cpp

Go to the documentation of this file.

/**
 * @file MatElastic.cpp
 * @author your name (you@domain.com)
 * @brief 
 * @version 0.1
 * @date 2026-03-26
 * 
 * @copyright Copyright (c) 2026
 * 
 */
 
#include <MoFEM.hpp>
 
using namespace MoFEM;
 
#include "MatOps.hpp"
#include "MatElastic.hpp"
 
namespace MatOps {
 
template <int DIM> struct MatElasticImpl : public MatElastic {
  using MatElastic::MatElastic;
  ForcesAndSourcesCore::UserDataOperator *
  createOp(boost::shared_ptr<PhysicalEquations> physical_ptr, bool eval_stress,
           bool eval_tangent, bool update) override;
 
protected:
};
 
template <int DIM> struct TopoMatElasticImpl : public MatElastic {
  using MatElastic::MatElastic;
  ForcesAndSourcesCore::UserDataOperator *
  createOp(boost::shared_ptr<PhysicalEquations> physical_ptr, bool eval_stress,
           bool eval_tangent, bool update) override;
 
protected:
};
 
template <int DIM, typename DomainEleOp>
struct OpEvalMatSimpleMaterialImpl : public DomainEleOp {
 
  using OP = DomainEleOp;
 
  OpEvalMatSimpleMaterialImpl(
      boost::shared_ptr<PhysicalEquations> physical_ptr, bool eval_stress,
      bool eval_tangent, bool update)
      : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), physicalPtr(physical_ptr),
        evalStress(eval_stress), evalTangent(eval_tangent),
        updateState(update) {}
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
protected:
  boost::shared_ptr<PhysicalEquations> physicalPtr;
  const bool evalStress;
  const bool evalTangent;
  const bool updateState;
};
 
template <int DIM, typename DomainEleOp>
MoFEMErrorCode OpEvalMatSimpleMaterialImpl<DIM, DomainEleOp>::doWork(
    int side, EntityType type, EntData &data) {
  MoFEMFunctionBegin;
 
  int nb_integration_pts = OP::getGaussPts().size2();
 
  auto get_tag = [&]() {
    if (physicalPtr->tagVsRangePtr) {
      for (const auto &tag_range_pair : *(physicalPtr->tagVsRangePtr)) {
        if (tag_range_pair.second.find(DomainEleOp::getFEEntityHandle()) !=
            tag_range_pair.second.end()) {
          return tag_range_pair.first;
        }
      }
    }
#ifndef NDEBUG
    if (MatOpsTagsRegistry::getTagName(physicalPtr->tAg).empty()) {
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                        "ADOL-C tag not found " +
                            std::to_string(physicalPtr->tAg));
    }
#endif
    return physicalPtr->tAg; // Default tag if no range matches
  };
 
  const int current_tag = get_tag();
  auto *fe_ptr = const_cast<FEMethod *>(this->getFEMethod());
  const auto ent = fe_ptr->getFEEntityHandle();
 
  physicalPtr->matOpsDataPtr->getActiveDataPtr("F")->resize(DIM, DIM, false);
  physicalPtr->matOpsDataPtr->getDependentDataPtr("P")->resize(DIM, DIM, false);
  physicalPtr->matOpsDataPtr->getDependentDerivativesDataPtr("P_dF")->resize(
      DIM * DIM, DIM * DIM, false);
  auto mat_grad_ptr = physicalPtr->matOpsDataPtr->getCommonDataPtr("grad");
  auto mat_P_ptr = physicalPtr->matOpsDataPtr->getCommonDataPtr("P");
  auto mat_P_dF_ptr = physicalPtr->matOpsDataPtr->getCommonDataPtr("P_dF");
 
#ifndef NDEBUG
  if (!mat_grad_ptr || !mat_P_ptr || !mat_P_dF_ptr) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "Missing common data for ADOL-C evaluation");
  }
  if (mat_grad_ptr->size2() != DIM * DIM) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "Inconsistent size of gradient matrix for ADOL-C evaluation");
  }
  if (mat_grad_ptr->size1() != nb_integration_pts) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "Inconsistent size of gradient matrix data for ADOL-C evaluation "
            "%zu != %d",
            mat_grad_ptr->size1(), nb_integration_pts);
  }
#endif
 
  using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
  auto get_grad_at_pts =
      MatrixSizeHelper<GetFTensor2FromMatType<DIM, DIM, -1, DL>, DL>::get(
          *mat_grad_ptr, nb_integration_pts);
  auto get_P_at_pts =
      MatrixSizeHelper<GetFTensor2FromMatType<DIM, DIM, -1, DL>, DL>::size(
          *mat_P_ptr, nb_integration_pts);
  auto get_P_dF_at_pts =
      MatrixSizeHelper<GetFTensor4FromMatType<DIM, DIM, DIM, DIM, -1, DL>,
                       DL>::size(*mat_P_dF_ptr, nb_integration_pts);
 
  if (evalStress) {
    FTENSOR_INDEX(DIM, i);
    FTENSOR_INDEX(DIM, J);
 
    auto t_grad_at_pts = get_grad_at_pts();
    auto t_P_at_pts = get_P_at_pts();
 
    auto next = [&]() {
      ++t_grad_at_pts;
      ++t_P_at_pts;
    };
 
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
      auto t_F = getFTensor2FromPtr<DIM, DIM>(
          physicalPtr->matOpsDataPtr->getActiveDataPtr("F")->data().data());
      t_F(i, J) = t_grad_at_pts(i, J);
      CHKERR physicalPtr->setParams(fe_ptr, gg);
      CHKERR physicalPtr->evaluateVariable(current_tag, ent, gg);
      auto t_P = getFTensor2FromPtr<DIM, DIM>(
          physicalPtr->matOpsDataPtr->getDependentDataPtr("P")->data().data());
      t_P_at_pts(i, J) = t_P(i, J);
      next();
    }
  }
 
  if (evalTangent) {
    FTENSOR_INDEX(DIM, i);
    FTENSOR_INDEX(DIM, J);
    FTENSOR_INDEX(DIM, k);
    FTENSOR_INDEX(DIM, L);
 
    auto t_grad_at_pts = get_grad_at_pts();
    auto t_P_dF_at_pts = get_P_dF_at_pts();
    auto next = [&]() {
      ++t_grad_at_pts;
      ++t_P_dF_at_pts;
    };
 
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
      auto t_F = getFTensor2FromPtr<DIM, DIM>(
          physicalPtr->matOpsDataPtr->getActiveDataPtr("F")->data().data());
      t_F(i, J) = t_grad_at_pts(i, J);
      CHKERR physicalPtr->setParams(fe_ptr, gg);
      CHKERR physicalPtr->evaluateDerivatives(current_tag, ent, gg);
      auto t_P_dF = getFTensor4FromPtr<DIM, DIM, DIM, DIM>(
          physicalPtr->matOpsDataPtr->getDependentDerivativesDataPtr("P_dF")
              ->data()
              .data());
      t_P_dF_at_pts(i, J, k, L) = t_P_dF(i, J, k, L);
      next();
    }
  }
 
  if (updateState) {
    FTENSOR_INDEX(DIM, i);
    FTENSOR_INDEX(DIM, J);
 
    auto t_grad_at_pts = get_grad_at_pts();
    auto next = [&]() { ++t_grad_at_pts; };
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      auto t_F = getFTensor2FromPtr<DIM, DIM>(
          physicalPtr->matOpsDataPtr->getActiveDataPtr("F")->data().data());
      t_F(i, J) = t_grad_at_pts(i, J);
      CHKERR physicalPtr->setParams(fe_ptr, gg);
      CHKERR physicalPtr->updateState(current_tag, ent, gg);
      next();
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM, typename DomainEleOp>
struct OpEvalTopoMatSimpleMaterialImpl : public DomainEleOp {
 
  using OP = DomainEleOp;
 
  OpEvalTopoMatSimpleMaterialImpl(
      boost::shared_ptr<PhysicalEquations> physical_ptr, bool eval_stress,
      bool eval_tangent, bool update)
      : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), physicalPtr(physical_ptr),
        evalStress(eval_stress), evalTangent(eval_tangent),
        updateState(update) {}
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
protected:
  boost::shared_ptr<PhysicalEquations> physicalPtr;
  const bool evalStress;
  const bool evalTangent;
  const bool updateState;
};
 
template <int DIM, typename DomainEleOp>
MoFEMErrorCode OpEvalTopoMatSimpleMaterialImpl<DIM, DomainEleOp>::doWork(
    int side, EntityType type, EntData &data) {
  MoFEMFunctionBegin;
 
  int nb_integration_pts = OP::getGaussPts().size2();
 
  auto get_tag = [&]() {
    if (physicalPtr->tagVsRangePtr) {
      for (const auto &tag_range_pair : *(physicalPtr->tagVsRangePtr)) {
        if (tag_range_pair.second.find(DomainEleOp::getFEEntityHandle()) !=
            tag_range_pair.second.end()) {
          return tag_range_pair.first;
        }
      }
    }
#ifndef NDEBUG
    if (MatOpsTagsRegistry::getTagName(physicalPtr->tAg).empty()) {
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                        "ADOL-C tag not found " +
                            std::to_string(physicalPtr->tAg));
    }
#endif
    return physicalPtr->tAg; // Default tag if no range matches
  };
 
  const int current_tag = get_tag();
  auto *fe_ptr = const_cast<FEMethod *>(this->getFEMethod());
  const auto ent = fe_ptr->getFEEntityHandle();
 
  physicalPtr->matOpsDataPtr->getActiveDataPtr("Jac")->resize(DIM, DIM, false);
  physicalPtr->matOpsDataPtr->getDependentDataPtr("Obj")->resize(1, false);
  physicalPtr->matOpsDataPtr->getDependentDerivativesDataPtr("dObj")->resize(
      DIM, DIM, false);
  auto mat_jac_ptr = physicalPtr->matOpsDataPtr->getCommonDataPtr("Jac");
  auto mat_grad_ptr = physicalPtr->matOpsDataPtr->getCommonDataPtr("grad");
  auto mat_var_grad_ptr = physicalPtr->matOpsDataPtr->getCommonDataPtr("var_grad");
  auto mat_obj_ptr = physicalPtr->matOpsDataPtr->getCommonDataPtr("Obj");
  auto mat_dObj_ptr = physicalPtr->matOpsDataPtr->getCommonDataPtr("dObj");
 
#ifndef NDEBUG
  if (!mat_jac_ptr || !mat_grad_ptr || !mat_var_grad_ptr || !mat_obj_ptr ||
      !mat_dObj_ptr) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "Missing common data for ADOL-C evaluation");
  }
#endif
 
  using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
  auto get_jac_at_pts =
      MatrixSizeHelper<GetFTensor2FromMatType<DIM, DIM, -1, DL>, DL>::get(
          *mat_jac_ptr, nb_integration_pts);
  [[maybe_unused]] auto get_grad_at_pts =
      MatrixSizeHelper<GetFTensor2FromMatType<DIM, DIM, -1, DL>, DL>::get(
          *mat_grad_ptr, nb_integration_pts);
  [[maybe_unused]] auto get_var_grad_at_pts =
      MatrixSizeHelper<GetFTensor2FromMatType<DIM, DIM, -1, DL>, DL>::get(
          *mat_var_grad_ptr, nb_integration_pts);
  auto get_obj_at_pts =
      MatrixSizeHelper<GetFTensor1FromMatType<1, -1, DL>, DL>::size(
          *mat_obj_ptr, nb_integration_pts);
  auto get_dObj_at_pts =
      MatrixSizeHelper<GetFTensor2FromMatType<DIM, DIM, -1, DL>, DL>::size(
          *mat_dObj_ptr, nb_integration_pts);
 
  if (evalStress) {
    FTENSOR_INDEX(DIM, i);
    FTENSOR_INDEX(DIM, J);
 
    auto t_jac_at_pts = get_jac_at_pts();
    auto t_obj_at_pts = get_obj_at_pts();
 
    auto next = [&]() {
      ++t_jac_at_pts;
      ++t_obj_at_pts;
    };
 
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
      auto t_jac = getFTensor2FromPtr<DIM, DIM>(
          physicalPtr->matOpsDataPtr->getActiveDataPtr("Jac")->data().data());
      t_jac(i, J) = t_jac_at_pts(i, J);
      CHKERR physicalPtr->setParams(fe_ptr, gg);
      CHKERR physicalPtr->evaluateVariable(current_tag, ent, gg);
      auto t_obj = getFTensor1FromPtr<1, 1>(
          physicalPtr->matOpsDataPtr->getDependentDataPtr("Obj")
              ->data()
              .data());
      t_obj_at_pts(0) = t_obj(0);
      next();
    }
  }
 
  if (evalTangent) {
    FTENSOR_INDEX(DIM, i);
    FTENSOR_INDEX(DIM, J);
    FTENSOR_INDEX(DIM, k);
    FTENSOR_INDEX(DIM, L);
 
    auto t_jac_at_pts = get_jac_at_pts();
    auto t_dObj_at_pts = get_dObj_at_pts();
    auto next = [&]() {
      ++t_jac_at_pts;
      ++t_dObj_at_pts;
    };
 
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
      auto t_jac = getFTensor2FromPtr<DIM, DIM>(
          physicalPtr->matOpsDataPtr->getActiveDataPtr("Jac")->data().data());
      t_jac(i, J) = t_jac_at_pts(i, J);
      CHKERR physicalPtr->setParams(fe_ptr, gg);
      CHKERR physicalPtr->evaluateDerivatives(current_tag, ent, gg);
      auto t_dObj = getFTensor2FromPtr<DIM, DIM>(
          physicalPtr->matOpsDataPtr->getDependentDerivativesDataPtr("dObj")
              ->data()
              .data());
      t_dObj_at_pts(i, J) = t_dObj(i, J);
      next();
    }
  }
 
  if (updateState) {
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM, typename EleOp>
static EleOp *createOpImpl(boost::shared_ptr<PhysicalEquations> physical_ptr,
                           bool eval_stress, bool eval_tangent, bool update) {
  return new OpEvalMatSimpleMaterialImpl<DIM, EleOp>(
      physical_ptr, eval_stress, eval_tangent, update);
}
 
template <>
ForcesAndSourcesCore::UserDataOperator *
MatElasticImpl<3>::createOp(boost::shared_ptr<PhysicalEquations> physical_ptr,
                              bool eval_stress, bool eval_tangent,
                              bool update) {
  using EleOp = ForcesAndSourcesCore::UserDataOperator;
  return createOpImpl<3, EleOp>(physical_ptr, eval_stress, eval_tangent,
                                update);
}
 
template <>
ForcesAndSourcesCore::UserDataOperator *
MatElasticImpl<2>::createOp(boost::shared_ptr<PhysicalEquations> physical_ptr,
                              bool eval_stress, bool eval_tangent,
                              bool update) {
  using EleOp = ForcesAndSourcesCore::UserDataOperator;
  return createOpImpl<2, EleOp>(physical_ptr, eval_stress, eval_tangent,
                                update);
}
 
template <int DIM, typename EleOp>
static EleOp *
createTopoOpImpl(boost::shared_ptr<PhysicalEquations> physical_ptr,
                 bool eval_stress, bool eval_tangent, bool update) {
  return new OpEvalTopoMatSimpleMaterialImpl<DIM, EleOp>(
      physical_ptr, eval_stress, eval_tangent, update);
}
 
template <>
ForcesAndSourcesCore::UserDataOperator *TopoMatElasticImpl<3>::createOp(
    boost::shared_ptr<PhysicalEquations> physical_ptr, bool eval_stress,
    bool eval_tangent, bool update) {
  using EleOp = ForcesAndSourcesCore::UserDataOperator;
  return createTopoOpImpl<3, EleOp>(physical_ptr, eval_stress, eval_tangent,
                                update);
}
 
template <>
ForcesAndSourcesCore::UserDataOperator *TopoMatElasticImpl<2>::createOp(
    boost::shared_ptr<PhysicalEquations> physical_ptr, bool eval_stress,
    bool eval_tangent, bool update) {
  using EleOp = ForcesAndSourcesCore::UserDataOperator;
  return createTopoOpImpl<2, EleOp>(physical_ptr, eval_stress, eval_tangent,
                                    update);
}
 
} // namespace MatOps
 
#include "MatNeohookean.cpp"
#include "MatMetaElastic.cpp"
#include "MatVolumeLengthQuality.cpp"
#include "MatMooneyRivlinWriggersEq63.cpp"
#include "MatStVenantKirchhoff.cpp"
#include "MatGenericElastic.cpp"