MatNeohookean.cpp

Go to the documentation of this file.

/**
 * @file MatNeohookean.cpp
 * @author your name (you@domain.com)
 * @brief Abaqus-compatible compressible Neo-Hookean material
 * @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;
 
  static inline double getShearModulus(const double c10) {
    return 2. * c10;
  }
 
  static inline double getAbaqusD1(const double K) {
    return 2. / K;
  }
 
  static MoFEMErrorCode validateMaterialParameters(const double c10,
                                                   const double K,
                                                   const char *source) {
    MoFEMFunctionBegin;
    if (!std::isfinite(c10) || c10 <= 0.) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
              "Neohookean C10 must be finite and positive in %s", source);
    }
    if (!std::isfinite(K) || K <= 0.) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
              "Neohookean bulk modulus K must be finite and positive in %s "
              "(received C10=%.16g, K=%.16g)",
              source, c10, K);
    }
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode getOptions(MoFEM::Interface *m_field_ptr = nullptr) override {
    MoFEMFunctionBegin;
 
    MOFEM_LOG_CHANNEL("WORLD");
 
    if (!A::paramVecByRange.empty()) {
      for (const auto &[range, parameters] : A::paramVecByRange) {
        (void)range;
        if (parameters.size() != 2) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "Neohookean material range needs exactly C10 and K");
        }
        CHKERR validateMaterialParameters(parameters[0], parameters[1],
                                          "preconfigured material range");
      }
      MoFEMFunctionReturnHot(0);
    }
 
    PetscOptionsBegin(PETSC_COMM_WORLD, "neo_hookean_", "", "none");
    CHKERR PetscOptionsScalar("-c10", "C10", "", C10, &C10, PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-K", "Bulk modulus K", "", K, &K,
                              PETSC_NULLPTR);
    PetscOptionsEnd();
    CHKERR validateMaterialParameters(C10, K, "command-line parameters");
    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);
      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;
      };
 
      CHKERR addBlockParameters(*m_field_ptr, m->getName(),
                                m->getMeshsetId(), get_block_ents(),
                                block_data);
    }
 
    if (A::paramVecByRange.empty()) {
      commandLineParameters = {C10, K};
      MOFEM_TAG_AND_LOG("WORLD", Sev::inform,
                        "Command line Neohookean parameters")
          << "C10 = " << C10 << " K = " << K
          << " mu = " << getShearModulus(C10)
          << " Abaqus D1 = " << getAbaqusD1(K);
    }
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode
  addBlockParameters(MoFEM::Interface &m_field, const std::string &block_name,
                     int block_id, const Range &block_entities,
                     const std::vector<double> &block_data) override {
    MoFEMFunctionBegin;
 
    const auto json_params =
        m_field.getInterface<JsonConfigManager>()->getParamsFromBlockset(
            "MAT_NEOHOOKEAN", block_id);
    std::vector<double> params;
    if (!json_params.empty()) {
      if (json_params.size() != 2 ||
          json_params.find("c10") == json_params.end() ||
          json_params.find("k") == json_params.end()) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Neohookean JSON block %s needs exactly C10 and K attributes",
                block_name.c_str());
      }
      params = {json_params.at("c10"), json_params.at("k")};
    } else {
      if (block_data.size() < 2) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Neohookean block %s has %d attributes; expected C10 and K",
                block_name.c_str(), static_cast<int>(block_data.size()));
      }
      params = {block_data[0], block_data[1]};
    }
 
    CHKERR validateMaterialParameters(params[0], params[1],
                                      block_name.c_str());
 
    A::paramVecByRange.push_back({block_entities, params});
 
    MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "MatBlock for Neohookean")
        << block_name << " C10 = " << params[0] << " K = " << params[1]
        << " mu = " << getShearModulus(params[0])
        << " Abaqus D1 = " << getAbaqusD1(params[1]);
 
    MoFEMFunctionReturn(0);
  }
 
protected:
  double C10 = 1;
  double K = 2;
  std::vector<double> defaultMaterialParameters = {C10, K};
  std::vector<double> commandLineParameters;
};
 
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;
      }
    }
    // command line overrides block parameters
    if (!this->commandLineParameters.empty()) {
      set_param_vec(A::tAg, this->commandLineParameters.size(),
                    this->commandLineParameters.data());
      return 0;
    }
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Neohookean material parameters not found for entity %lu",
            static_cast<unsigned long>(ent));
  }
 
  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, I1;
    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);
    I1 = ta_F(i, I) * ta_F(i, I);
    if constexpr (DIM == 2) {
      // Abaqus plane strain uses the three-dimensional potential with F33=1.
      I1 += 1.;
    }
    CHKERR invertTensor(ta_F, det_aF, ta_invF);
 
    // Abaqus compressible Neo-Hookean potential (reduced polynomial N=1):
    //
    //   W = C10 * (I1_bar - 3) + (1 / D1) * (J - 1)^2,
    //   I1_bar = J^(-2/3) * I1,  mu0 = 2*C10,  D1 = 2/K0.
    //
    // Its first Piola stress is
    //
    //   P = 2*C10*J^(-2/3)*(F - I1/3*F^(-T))
    //       + K0*J*(J - 1)*F^(-T).
    const auto jacobian_to_minus_two_thirds = pow(det_aF, -2. / 3.);
    ta_P(i, I) =
        2. * p_c10 * jacobian_to_minus_two_thirds *
            (ta_F(i, I) - (I1 / 3.) * ta_invF(I, i)) +
        p_K * det_aF * (det_aF - 1.) * ta_invF(I, i);
 
    ta_P(i, I) >>= t_P(i, I);
 
    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);
}
 
} // namespace MatOps