MatNeohookean.cpp

Go to the documentation of this file.

/**
 * @file MatNeohookean.cpp
 * @author your name (you@domain.com)
 * @brief 
 * @version 0.1
 * @date 2026-03-26
 * 
 * @copyright Copyright (c) 2026
 * 
 */
 
namespace MatOps {
template <typename M> struct MatNeohookeanOptions {
  static constexpr auto optionsPrefix = "neo_hookean_";
  static constexpr auto defaultLogTag = "Default Neohookean parameters:";
  static constexpr auto blockLogTag = "MatBlock for Neohookean";
};
 
template <int DIM> struct MatNeohookeanOptions<TopoMatElasticImpl<DIM>> {
  static constexpr auto optionsPrefix = "topo_neo_hookean_";
  static constexpr auto defaultLogTag = "Default Topo Neohookean parameters:";
  static constexpr auto blockLogTag = "MatBlock for Topo Neohookean";
};
 
template <typename M> struct MatNeohookeanGeneric : public M {
  using A = M;
  using A::A;
 
  MoFEMErrorCode getOptions(MoFEM::Interface *m_field_ptr = nullptr) override {
    MoFEMFunctionBegin;
 
    MOFEM_LOG_CHANNEL("WORLD");
 
    using Options = MatNeohookeanOptions<A>;
 
    PetscOptionsBegin(PETSC_COMM_WORLD, Options::optionsPrefix, "", "none");
    CHKERR PetscOptionsScalar("-c10", "C10", "", C10, &C10, PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-K", "Bulk modulus K", "", K, &K, PETSC_NULLPTR);
    PetscOptionsEnd();
 
    MOFEM_TAG_AND_LOG("WORLD", Sev::inform, Options::defaultLogTag)
        << " C10 = " << C10 << " K = " << K;
    defaultMaterialParameters = {C10, K};
 
    std::string block_name = "MAT_NEOHOOKEAN";
 
    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() < 2) {
        CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                          "Expected that block has two attribute");
      }
      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], block_data[1]}});
 
      MOFEM_TAG_AND_LOG("WORLD", Sev::inform, Options::blockLogTag)
          << *m << " C10 = " << block_data[0] << " K = " << block_data[1];
    }
 
    MoFEMFunctionReturn(0);
  };
 
protected:
  double C10 = 1;
  double K = 1;
  std::vector<double> defaultMaterialParameters = {C10, K};
};
 
template <int DIM>
struct MatNeohookean : public MatNeohookeanGeneric<MatElasticImpl<DIM>> {
  using A = MatNeohookeanGeneric<MatElasticImpl<DIM>>;
  using A::A;
 
  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, this->defaultMaterialParameters.size(),
                  this->defaultMaterialParameters.data());
    return 0;
  }
 
  MoFEMErrorCode recordTape() override {
    MoFEMFunctionBegin;
 
    A::matOpsDataPtr->insertCommonData("grad", MatrixDouble());
    A::matOpsDataPtr->insertCommonData("P", MatrixDouble());
    A::matOpsDataPtr->insertCommonData("P_dF", MatrixDouble());
 
    A::matOpsDataPtr->insertActiveData("F", MatrixDouble());
    A::matOpsDataPtr->insertDependentData("P", MatrixDouble());
    A::matOpsDataPtr->insertDependentDerivativesData("P_dF", MatrixDouble());
 
    A::matOpsDataPtr->getActiveDataPtr("F")->resize(DIM, DIM, false);
    A::matOpsDataPtr->getDependentDataPtr("P")->resize(DIM, DIM, false);
 
    auto t_F = getFTensor2FromPtr<DIM, DIM>(
        A::matOpsDataPtr->getActiveDataPtr("F")->data().data());
    auto t_P = getFTensor2FromPtr<DIM, DIM>(
        A::matOpsDataPtr->getDependentDataPtr("P")->data().data());
 
    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) = 0;
 
    FTensor::Tensor2<adouble, DIM, DIM> ta_F;
    FTensor::Tensor2<adouble, DIM, DIM> ta_P;
 
    adouble det_aF;
    FTensor::Tensor2<adouble, DIM, DIM> ta_invF;
 
    trace_on(A::tAg);
    auto p_c10 = mkparam(this->C10);
    auto p_K = mkparam(this->K);
 
    ta_F(i, J) <<= t_F(i, J);
    // assume that gradient from approximated displacement, that why we add diagonal
    if (!MatElastic::useDeformationGradient) {
      ta_F(i, J) += t_kd(i, J);
    }
 
    det_aF = determinantTensor(ta_F);
    CHKERR invertTensor(ta_F, det_aF, ta_invF);
 
    ta_P(i, I) = 2. * p_c10 * (ta_F(i, I) - ta_invF(i, I)) +
                 p_K * log(det_aF) * ta_invF(i, I);
 
    ta_P(i, I) >>= t_P(i, I);
 
    trace_off();
 
    MoFEMFunctionReturn(0);
  }
 
};
 
template <int DIM>
struct TopoDerivativeMatElasticImpl
    : public MatNeohookeanGeneric<TopoMatElasticImpl<DIM>> {
  using A = MatNeohookeanGeneric<TopoMatElasticImpl<DIM>>;
  using A::A;
 
  MoFEMErrorCode setParams(FEMethod *fe_ptr, int gg) override {
    (void)gg;
    std::vector<double> adolc_c_param_vec;
    auto grad = A::matOpsDataPtr->getCommonDataPtr("grad");
    auto var_grad = A::matOpsDataPtr->getCommonDataPtr("var_grad");
    adolc_c_param_vec.reserve(grad->size2() + var_grad->size2() +
                              this->defaultMaterialParameters.size());
 
    for (auto cc = 0; cc < grad->size2(); ++cc) {
      adolc_c_param_vec.push_back((*grad)(gg, cc));
    }
    for (auto cc = 0; cc < var_grad->size2(); ++cc) {
      adolc_c_param_vec.push_back((*var_grad)(gg, cc));
    }
    const auto ent = fe_ptr->getFEEntityHandle();
    for (auto &[range, param_vec] : A::paramVecByRange) {
      if (std::find(range.begin(), range.end(), ent) != range.end()) {
        adolc_c_param_vec.insert(adolc_c_param_vec.end(), param_vec.begin(),
                                 param_vec.end());
        set_param_vec(A::tAg, adolc_c_param_vec.size(),
                      adolc_c_param_vec.data());
        return 0;
      }
    }
    adolc_c_param_vec.insert(adolc_c_param_vec.end(),
                             this->defaultMaterialParameters.begin(),
                             this->defaultMaterialParameters.end());
    set_param_vec(A::tAg, adolc_c_param_vec.size(), adolc_c_param_vec.data());
    return 0;
  }
 
  MoFEMErrorCode recordTape() override {
    MoFEMFunctionBegin;
 
    A::matOpsDataPtr->insertCommonData("Jac", MatrixDouble());
    A::matOpsDataPtr->insertCommonData("grad", MatrixDouble());
    A::matOpsDataPtr->insertCommonData("var_grad", MatrixDouble());
    A::matOpsDataPtr->insertActiveData("Jac", MatrixDouble(DIM, DIM));
    A::matOpsDataPtr->insertDependentData("Obj", MatrixDouble(1, 1));
    A::matOpsDataPtr->insertDependentDerivativesData("dObj",
                                                     MatrixDouble(DIM, DIM));
 
    auto t_jac = getFTensor2FromPtr<DIM, DIM>(
        A::matOpsDataPtr->getActiveDataPtr("Jac")->data().data());
    auto &obj = (*A::matOpsDataPtr->getDependentDataPtr("Obj"))(0, 0);
 
    auto t_kd = FTensor::Kronecker_Delta<int>();
 
    FTENSOR_INDEXES(DIM, i, J);
    t_jac(i, J) = t_kd(i, J);
 
    FTensor::Tensor2<adouble, DIM, DIM> tp_F;
    FTensor::Tensor2<adouble, DIM, DIM> tp_var_F;
    FTensor::Tensor2<adouble, DIM, DIM> ta_jac;
    FTensor::Tensor2<adouble, DIM, DIM> ta_F;
    FTensor::Tensor2<adouble, DIM, DIM> ta_P;
    adouble det_aF;
    FTensor::Tensor2<adouble, DIM, DIM> ta_invF;
 
    adouble ta_obj;
    adouble det;
 
    trace_on(A::tAg);
 
    for (int ii = 0; ii < DIM; ++ii) {
      for (int JJ = 0; JJ < DIM; ++JJ) {
        tp_F(ii, JJ) = mkparam(0.);
      }
    }
    for (int ii = 0; ii < DIM; ++ii) {
      for (int JJ = 0; JJ < DIM; ++JJ) {
        tp_var_F(ii, JJ) = mkparam(0.);
      }
    }
    auto p_c10 = mkparam(this->C10);
    auto p_K = mkparam(this->K);
 
    FTENSOR_INDEX(DIM, I);
    FTENSOR_INDEX(DIM, K);
 
    ta_jac(i, J) <<= t_jac(i, J);
    det = determinantTensor(ta_jac);
 
    if (MatElastic::useDeformationGradient) {
      ta_F(i, J) = ta_jac(J, K) * tp_F(i, K);
    } else {
      ta_F(i, J) = ta_jac(J, K) * (tp_F(i, K) + t_kd(i, K));
    }
 
    det_aF = determinantTensor(ta_F);
    CHKERR invertTensor(ta_F, det_aF, ta_invF);
    ta_P(i, I) = 2. * p_c10 * (ta_F(i, I) - ta_invF(i, I)) +
                 p_K * log(det_aF) * ta_invF(i, I);
    ta_obj = det * ta_P(i, I) * ta_jac(I, K) * tp_var_F(i, K);
    ta_obj >>= obj;
 
    trace_off();
 
    MoFEMFunctionReturn(0);
  }
 
};
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<ELASTICITY::NEOHOOKEAN, MODEL_3D>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag) {
  return boost::make_shared<MatNeohookean<3>>(mat_ops_data_ptr, tag);
}
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<ELASTICITY::NEOHOOKEAN, MODEL_2D_PLANE_STRAIN>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag) {
  return boost::make_shared<MatNeohookean<2>>(mat_ops_data_ptr, tag);
}
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<ELASTICITY::TOPO_NEOHOOKEAN, MODEL_3D>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag) {
  return boost::make_shared<TopoDerivativeMatElasticImpl<3>>(mat_ops_data_ptr,
                                                             tag);
}
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<ELASTICITY::TOPO_NEOHOOKEAN,
                                 MODEL_2D_PLANE_STRAIN>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag) {
  return boost::make_shared<TopoDerivativeMatElasticImpl<2>>(mat_ops_data_ptr,
                                                             tag);
}
} // namespace MatOps