MatHuHu.cpp

Go to the documentation of this file.

/**
 * @file MatHuHu.cpp
 * @author your name (you@domain.com)
 * @version 0.1
 * @date 2026-03-29
 * 
 * @copyright Copyright (c) 2026
 * 
 */
 
#include <MoFEM.hpp>
 
using namespace MoFEM;
 
#include "MatOps.hpp"
#include "MatHuHu.hpp"
 
namespace MatOps {
 
template <int DIM>
struct MatHuHu : public PhysicalEquations {
  using PhysicalEquations::PhysicalEquations;
  using A = PhysicalEquations;
  ForcesAndSourcesCore::UserDataOperator *
  createOp(boost::shared_ptr<PhysicalEquations> physical_ptr, bool eval_stress,
           bool eval_tangent, bool update) override;
 
  MoFEMErrorCode getOptions(MoFEM::Interface *m_field_ptr = nullptr) override {
    MoFEMFunctionBegin;
    MOFEM_LOG_CHANNEL("WORLD");
 
    PetscOptionsBegin(PETSC_COMM_WORLD, "huhu_", "", "none");
    CHKERR PetscOptionsScalar("-K", "Bulk modulus K", "", K, &K, PETSC_NULLPTR);
    PetscOptionsEnd();
 
    MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "Default HuHu parameters:")
        << " K = " << K;
    defaultMaterialParameters = {K};
 
    std::string block_name = "MAT_HUHU";
 
    for (auto &m :
 
         m_field_ptr->getInterface<MeshsetsManager>()->
 
         getCubitMeshsetPtr(
             std::regex((boost::format("%s(.*)") % block_name).str()))
 
    ) {
 
      std::vector<double> block_data;
      CHKERR m->getAttributes(block_data);
      if (block_data.size() < 1) {
        block_data.push_back(
            K); // Use default K if not provided in block attributes
      }
      auto get_block_ents = [&]() {
        Range ents;
        CHK_MOAB_THROW(m_field_ptr->get_moab().get_entities_by_handle(
                           m->meshset, ents, true),
                       "can not get block entities");
        return ents;
      };
 
      A::paramVecByRange.push_back({get_block_ents(), {block_data[0]}});
 
      MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "MatBlock for HuHu")
          << *m << " K = " << block_data[0];
    }
 
    MoFEMFunctionReturn(0);
  };
 
  MoFEMErrorCode setParams(FEMethod *fe_ptr, int gg) override {
    (void)gg;
    const auto ent = fe_ptr->getFEEntityHandle();
    for (auto &[range, param_vec] : A::paramVecByRange) {
      if (std::find(range.begin(), range.end(), ent) != range.end()) {
        set_param_vec(A::tAg, param_vec.size(), param_vec.data());
        return 0;
      }
    }
    set_param_vec(A::tAg, defaultMaterialParameters.size(),
                  defaultMaterialParameters.data());
    return 0;
  }
 
  MoFEMErrorCode recordTape() override {
    MoFEMFunctionBegin;
 
    A::matOpsDataPtr->insertCommonData("grad", MatrixDouble(DIM * DIM, 1));
    A::matOpsDataPtr->insertCommonData("k", MatrixDouble(1, 1));
    A::matOpsDataPtr->insertCommonData("k_dF", MatrixDouble(DIM * DIM, 1));
 
    A::matOpsDataPtr->insertActiveData("F", MatrixDouble(DIM, DIM));
    A::matOpsDataPtr->insertDependentData("k", MatrixDouble(1, 1));
    A::matOpsDataPtr->insertDependentDerivativesData("k_dF",
                                                    MatrixDouble(DIM, DIM));
 
    auto t_F = getFTensor2FromPtr<DIM, DIM>(
        A::matOpsDataPtr->getActiveDataPtr("F")->data().data());
    auto k = (*A::matOpsDataPtr->getDependentDataPtr("k"))(0, 0);
 
    auto t_kd = FTensor::Kronecker_Delta<int>();
 
    FTENSOR_INDEX(DIM, i);
    FTENSOR_INDEX(DIM, j);
    FTENSOR_INDEX(DIM, I);
    FTENSOR_INDEX(DIM, J);
 
    t_F(i, J) = t_kd(i, J);
 
    FTensor::Tensor2<adouble, DIM, DIM> ta_F;
    adouble ta_k;
 
    adouble det_aF;
    FTensor::Tensor2<adouble, DIM, DIM> ta_invF;
 
    trace_on(A::tAg);
    auto p_K = mkparam(K);
    ta_F(i, J) <<= t_F(i, J);
    ta_k = p_K;
    ta_k >>= k;
 
    trace_off();
 
    MoFEMFunctionReturn(0);
  }
 
protected:
  double K = 1e-6;
  std::vector<double> defaultMaterialParameters = {K};
};
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<HUHU, MODEL_3D>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag) {
  return boost::make_shared<MatHuHu<3>>(mat_ops_data_ptr, tag);
}
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<HUHU, MODEL_2D_PLANE_STRAIN>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag) {
  return boost::make_shared<MatHuHu<2>>(mat_ops_data_ptr, tag);
}
 
template <int DIM, typename DomainEleOp>
struct OpEvalMatHuHuImpl : public DomainEleOp {
 
  OpEvalMatHuHuImpl(boost::shared_ptr<PhysicalEquations> physical_ptr,
                      bool eval_stress, bool eval_tangent)
      : DomainEleOp("U", DomainEleOp::OPROW), 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 OpEvalMatHuHuImpl<DIM, DomainEleOp>::doWork(int side,
                                                             EntityType type,
                                                             EntData &data) {
  MoFEMFunctionBegin;
  int nb_dofs = data.getIndices().size();
  if (!nb_dofs)
    MoFEMFunctionReturnHot(0);
  int nb_integration_pts = data.getN().size1();
  auto *fe_ptr = const_cast<FEMethod *>(this->getFEMethod());
  const auto ent = fe_ptr->getFEEntityHandle();
 
  auto get_tag = [&]() {
    if (physicalPtr->tagVsRangePtr) {
      for (const auto &tag_range_pair : *physicalPtr->tagVsRangePtr) {
        if (tag_range_pair.second.find(ent) != tag_range_pair.second.end()) {
          return tag_range_pair.first;
        }
      }
    }
    return physicalPtr->tAg; // Default tag if no range matches
  };
 
  const int current_tag = get_tag();
 
  auto mat_grad_ptr = physicalPtr->matOpsDataPtr->getCommonDataPtr("grad");
  auto mat_k_ptr = physicalPtr->matOpsDataPtr->getCommonDataPtr("k");
  auto mat_k_dF_ptr = physicalPtr->matOpsDataPtr->getCommonDataPtr("k_dF");
 
#ifndef NDEBUG
  if (!mat_grad_ptr || !mat_k_ptr || !mat_k_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");
  }
#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_k_at_pts =
      MatrixSizeHelper<GetFTensor1FromMatType<1, -1, DL>, DL>::size(
          *mat_k_ptr, nb_integration_pts);
  auto get_k_dF_at_pts =
      MatrixSizeHelper<GetFTensor2FromMatType<DIM, DIM, -1, DL>, DL>::size(
          *mat_k_dF_ptr, nb_integration_pts);
 
  if (evalStress) {
    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_k_at_pts = get_k_at_pts();
 
    auto next = [&]() {
      ++t_grad_at_pts;
      ++t_k_at_pts;
    };
 
    auto t_kd = FTensor::Kronecker_Delta<int>();
 
    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) + t_kd(i, J);
      CHKERR physicalPtr->setParams(fe_ptr, gg);
      CHKERR physicalPtr->evaluateVariable(current_tag, ent, gg);
      t_k_at_pts(0) =
          (*physicalPtr->matOpsDataPtr->getDependentDataPtr("k"))(0, 0);
      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_k_dF_at_pts = get_k_dF_at_pts();
    auto t_kd = FTensor::Kronecker_Delta<int>();
 
    auto next = [&]() {
      ++t_grad_at_pts;
      ++t_k_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) + t_kd(i, J);
      CHKERR physicalPtr->setParams(fe_ptr, gg);
      CHKERR physicalPtr->evaluateDerivatives(current_tag, ent, gg);
      auto t_k_dF = getFTensor2FromPtr<DIM, DIM>(
          physicalPtr->matOpsDataPtr->getDependentDerivativesDataPtr("k_dF")
              ->data()
              .data());
      t_k_dF_at_pts(i, J) = t_k_dF(i, J);
      next();
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM, typename EleOp>
static EleOp *createOpImpl(boost::shared_ptr<PhysicalEquations> physical_ptr,
                           bool eval_stress, bool eval_tangent, bool update) {
  (void)update;
  return new OpEvalMatHuHuImpl<DIM, EleOp>(physical_ptr, eval_stress,
                                             eval_tangent);
}
 
template <>
ForcesAndSourcesCore::UserDataOperator *
MatHuHu<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 *
MatHuHu<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);
}
 
} // namespace MatOps