AdolCElastic.cpp

Go to the documentation of this file.

/**
 * @file AdolCElastic.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 "AdolCOps.hpp"
#include "AdolCElastic.hpp"
 
namespace AdolCOps {
 
template <int DIM> struct AdolCElasticImpl : public PhysicalEquations {
  using PhysicalEquations::PhysicalEquations;
 
  ForcesAndSourcesCore::UserDataOperator *
  createOp(boost::shared_ptr<PhysicalEquations> physical_ptr, bool eval_stress,
           bool eval_tangent) override;
 
protected:
};
 
template <int DIM, typename DomainEleOp>
struct OpEvalAdolCSimpleMaterialImpl : public DomainEleOp {
 
  using OP = DomainEleOp;
 
  OpEvalAdolCSimpleMaterialImpl(
      boost::shared_ptr<PhysicalEquations> physical_ptr, bool eval_stress,
      bool eval_tangent, bool add_diagonal = false)
      : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), physicalPtr(physical_ptr),
        evalStress(eval_stress), evalTangent(eval_tangent) {}
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
protected:
  boost::shared_ptr<PhysicalEquations> physicalPtr;
  const bool evalStress;
  const bool evalTangent;
};
 
template <int DIM, typename DomainEleOp>
MoFEMErrorCode OpEvalAdolCSimpleMaterialImpl<DIM, DomainEleOp>::doWork(
    int side, EntityType type, EntData &data) {
  MoFEMFunctionBegin;
 
  int nb_integration_pts = OP::getGaussPts().size2();
 
  constexpr int adolc_return_value = 0;
 
  auto physics_ptr = physicalPtr;
 
  auto get_tag = [&]() {
    if (physics_ptr->tagVsRangePtr) {
      for (const auto &tag_range_pair : *physics_ptr->tagVsRangePtr) {
        if (tag_range_pair.second.find(DomainEleOp::getFEEntityHandle()) !=
            tag_range_pair.second.end()) {
          return tag_range_pair.first;
        }
      }
    }
    return physics_ptr->tAg; // Default tag if no range matches
  };
 
  auto mat_grad_ptr = physics_ptr->adolcDataPtr->getCommonDataPtr("grad");
  auto mat_P_ptr = physics_ptr->adolcDataPtr->getCommonDataPtr("P");
  auto mat_P_dF_ptr = physics_ptr->adolcDataPtr->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->size1() != DIM * DIM) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "Inconsistent size of gradient matrix for ADOL-C evaluation");
  }
  if (mat_grad_ptr->size2() != nb_integration_pts) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "Inconsistent size of gradient matrix data for ADOL-C evaluation "
            "%lu != %lu",
            mat_grad_ptr->size2(), nb_integration_pts);
  }
#endif
 
  if (evalStress) {
    FTENSOR_INDEX(DIM, i);
    FTENSOR_INDEX(DIM, J);
 
    auto t_grad_at_pts = getFTensor2FromMat<DIM, DIM>(*mat_grad_ptr);
    mat_P_ptr->resize(DIM * DIM, nb_integration_pts, false);
    auto t_P_at_pts = getFTensor2FromMat<DIM, DIM>(*mat_P_ptr);
 
    auto next = [&]() {
      ++t_grad_at_pts;
      ++t_P_at_pts;
    };
 
    auto t_kd = FTensor::Kronecker_Delta<int>();
 
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
      auto t_F = getFTensor2FromPtr<DIM, DIM>(
          physics_ptr->adolcDataPtr->getActiveDataPtr("F")->data().data());
      t_F(i, J) = t_grad_at_pts(i, J);
      CHKERR physics_ptr->setActiveContinuousVector();
      CHKERR physics_ptr->setDependentContinuousVector();
      CHKERR physics_ptr->setParams(this->getFEEntityHandle(), gg);
      int r = ::function(get_tag(), physics_ptr->dependentVariables.size(),
                         physics_ptr->activeVariables.size(),
                         physics_ptr->activeVariables.data(),
                         physics_ptr->dependentVariables.data());
      if (PetscUnlikely(r < adolc_return_value)) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                "ADOL-C function evaluation with error");
      }
      CHKERR physics_ptr->getDependentContinuousVector();
      auto t_P = getFTensor2FromPtr<DIM, DIM>(
          physics_ptr->adolcDataPtr->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);
 
    const int number_of_active_variables = physics_ptr->activeVariables.size();
    const int number_of_dependent_variables =
        physics_ptr->dependentVariables.size();
    CHKERR physics_ptr->setDependentDerivativesContinuousVector();
    std::vector<double *> jac_ptr(number_of_dependent_variables);
    for (unsigned int n = 0; n != number_of_dependent_variables; ++n) {
      jac_ptr[n] = &(
          physics_ptr
              ->dependentVariablesDerivatives[n * number_of_active_variables]);
    }
    auto t_grad_at_pts = getFTensor2FromMat<DIM, DIM>(*mat_grad_ptr);
    mat_P_dF_ptr->resize(DIM * DIM * DIM * DIM, nb_integration_pts, false);
    auto t_P_dF_at_pts =
        getFTensor4FromMat<DIM, DIM, DIM, DIM, 1>(*mat_P_dF_ptr);
    auto t_kd = FTensor::Kronecker_Delta<int>();
 
    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>(
          physics_ptr->adolcDataPtr->getActiveDataPtr("F")->data().data());
      t_F(i, J) = t_grad_at_pts(i, J);
      CHKERR physics_ptr->setActiveContinuousVector();
      CHKERR physics_ptr->setDependentContinuousVector();
      CHKERR physics_ptr->setParams(this->getFEEntityHandle(), gg);
      int r = ::jacobian(get_tag(), number_of_dependent_variables,
                         number_of_active_variables,
                         physics_ptr->activeVariables.data(), jac_ptr.data());
      if (PetscUnlikely(r < adolc_return_value)) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                "ADOL-C function evaluation with error");
      }
      CHKERR physics_ptr->getDependentDerivativesContinuousVector();
      auto t_P_dF = getFTensor4FromPtr<DIM, DIM, DIM, DIM>(
          physics_ptr->adolcDataPtr->getDependentDerivativesDataPtr("P_dF")
              ->data()
              .data());
      t_P_dF_at_pts(i, J, k, L) = t_P_dF(i, J, k, L);
      next();
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM, typename EleOp>
static EleOp *createOpImpl(boost::shared_ptr<PhysicalEquations> physical_ptr,
                           bool eval_stress, bool eval_tangent) {
  return new OpEvalAdolCSimpleMaterialImpl<DIM, EleOp>(
      physical_ptr, eval_stress, eval_tangent);
}
 
template <>
ForcesAndSourcesCore::UserDataOperator *
AdolCElasticImpl<3>::createOp(boost::shared_ptr<PhysicalEquations> physical_ptr,
                              bool eval_stress, bool eval_tangent) {
  using EleOp = ForcesAndSourcesCore::UserDataOperator;
  return createOpImpl<3, EleOp>(physical_ptr, eval_stress, eval_tangent);
}
 
template <>
ForcesAndSourcesCore::UserDataOperator *
AdolCElasticImpl<2>::createOp(boost::shared_ptr<PhysicalEquations> physical_ptr,
                              bool eval_stress, bool eval_tangent) {
  using EleOp = ForcesAndSourcesCore::UserDataOperator;
  return createOpImpl<2, EleOp>(physical_ptr, eval_stress, eval_tangent);
}
 
} // namespace AdolCOps
 
#include "AdolCNeohookean.cpp"
#include "AdolCMetaElastic.cpp"
#include "AdolCVolumeLengthQuality.cpp"