MatMooneyRivlinWriggersEq63.cpp

Go to the documentation of this file.

/**
 * @file MatMooneyRivlinWriggersEq63.cpp
 * @author your name (you@domain.com)
 * @brief 
 * @version 0.1
 * @date 2026-03-26
 * 
 * @copyright Copyright (c) 2026
 * 
 */
 
namespace MatOps {
 
template <int DIM>
struct MatMooneyRivlinWriggersEq63 : public MatElasticImpl<DIM> {
  using MatElasticImpl<DIM>::MatElasticImpl;
 
  using A = MatElasticImpl<DIM>;
 
  MoFEMErrorCode getOptions(MoFEM::Interface *m_field_ptr = nullptr) override {
    MoFEMFunctionBegin;
 
    MOFEM_LOG_CHANNEL("WORLD");
 
    PetscOptionsBegin(PETSC_COMM_WORLD, "mooneyrivlin_", "", "none");
    CHKERR PetscOptionsScalar("-alpha", "Alpha", "", alpha, &alpha,
                              PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-beta", "Beta", "", beta, &beta, PETSC_NULLPTR);
 
    CHKERR PetscOptionsScalar("-lambda", "Lambda", "", lambda, &lambda,
                              PETSC_NULLPTR);
 
    CHKERR PetscOptionsScalar("-epsilon", "Epsilon", "", epsilon, &epsilon,
                              PETSC_NULLPTR);
    PetscOptionsEnd();
 
    if (!std::isfinite(epsilon) ||
        std::abs(epsilon) <= std::sqrt(std::numeric_limits<double>::epsilon()))
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
              "Mooney-Rivlin epsilon must be finite and nonzero");
 
    MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "Default Mooney-Rivlin parameters:")
        << " alpha = " << alpha << " beta = " << beta << " lambda = " << lambda
        << " epsilon = " << epsilon;
    defaultMaterialParameters = {alpha, beta, lambda, epsilon};
 
    std::string block_name = "MAT_MOONEY_RIVLIN";
 
    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);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  MoFEMErrorCode
  addBlockParameters(MoFEM::Interface &, const std::string &block_name, int,
                     const Range &block_entities,
                     const std::vector<double> &block_data) override {
    MoFEMFunctionBegin;
 
    if (block_data.size() < 4) {
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Expected that block has four attributes (alpha, "
                        "beta, lambda, epsilon), but given " +
                            std::to_string(block_data.size()));
    }
    if (!std::isfinite(block_data[3]) ||
        std::abs(block_data[3]) <=
            std::sqrt(std::numeric_limits<double>::epsilon()))
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
              "MAT_MOONEY_RIVLIN epsilon must be finite and nonzero");
    std::vector<double> params = {block_data[0], block_data[1], block_data[2],
                                  block_data[3]};
 
    A::paramVecByRange.push_back({block_entities, params});
 
    MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "MatBlock for Mooney-Rivlin")
        << block_name << " alpha = " << params[0] << " beta = " << params[1]
        << " lambda = " << params[2] << " epsilon = " << params[3];
 
    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());
    A::matOpsDataPtr->insertCommonData("P", MatrixDouble());
    A::matOpsDataPtr->insertCommonData("P_dF", MatrixDouble());
 
    A::matOpsDataPtr->insertActiveData("F", MatrixDouble());
    A::matOpsDataPtr->insertDependentData("P", MatrixDouble());
    A::matOpsDataPtr->insertDependentData("energy", MatrixDouble());
    A::matOpsDataPtr->insertDependentDerivativesData("P_dF", MatrixDouble());
    A::matOpsDataPtr->insertDependentDerivativesData("energy_dF",
                                                     MatrixDouble());
 
    A::matOpsDataPtr->getActiveDataPtr("F")->resize(DIM, DIM, false);
    A::matOpsDataPtr->getDependentDataPtr("P")->resize(DIM, DIM, false);
    A::matOpsDataPtr->getDependentDataPtr("energy")->resize(1, 1, false);
    A::matOpsDataPtr->getDependentDerivativesDataPtr("energy_dF")
        ->resize(1, 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 energy_ptr = A::matOpsDataPtr->getDependentDataPtr("energy");
 
    auto t_kd = FTensor::Kronecker_Delta<int>();
 
    FTENSOR_INDEX(DIM, i);
    FTENSOR_INDEX(DIM, j);
    FTENSOR_INDEX(DIM, I);
    FTENSOR_INDEX(DIM, J);
    FTENSOR_INDEX(DIM, k);
    FTENSOR_INDEX(DIM, K);
 
    t_F(i, J) = 0;
 
    FTensor::Tensor2<adouble, DIM, DIM> ta_F;
    FTensor::Tensor2<adouble, DIM, DIM> ta_P;
    FTensor::Tensor2<adouble, DIM, DIM> ta_Cof;
    FTensor::Tensor2<adouble, DIM, DIM> ta_BF;
    FTensor::Tensor2<adouble, DIM, DIM> ta_BCof;
    adouble ta_Bj, A, B, ta_energy;
 
    adouble det_aF;
    FTensor::Tensor2<adouble, DIM, DIM> ta_invF;
 
    trace_on(A::tAg);
    auto p_alpha = mkparam(alpha);
    auto p_beta = mkparam(beta);
    auto p_lambda = mkparam(lambda);
    auto p_epsilon = mkparam(epsilon);
    const auto ta_alpha = 1 * p_alpha;
    const auto ta_beta = 1 * p_beta;
    const auto ta_lambda = 1 * p_lambda;
    const auto ta_epsilon = 1 * p_epsilon;
 
    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_Cof(i, I) = det_aF * ta_invF(I, i);
 
    A = ta_F(k, K) * ta_F(k, K);
    B = ta_Cof(k, K) * ta_Cof(k, K);
 
    ta_BF(i, I) = 4 * ta_alpha * (A * ta_F(i, I));
    ta_BCof(i, I) = 4 * ta_beta * (B * ta_Cof(i, I));
    ta_Bj = (-12 * ta_alpha - 24 * ta_beta) / det_aF +
            0.5 * (ta_lambda / ta_epsilon) *
                (pow(det_aF, ta_epsilon - 1) - pow(det_aF, -ta_epsilon - 1));
 
    ta_P(i, I) = ta_BF(i, I);
    ta_P(i, I) += (levi_civita(i, j, k) * ta_BCof(j, J)) *
                  (levi_civita(I, J, K) * ta_F(k, K));
    ta_P(i, I) += ta_Cof(i, I) * ta_Bj;
 
    ta_energy = ta_alpha * A * A + ta_beta * B * B +
                (-12 * ta_alpha - 24 * ta_beta) * log(det_aF) +
                0.5 * (ta_lambda / (ta_epsilon * ta_epsilon)) *
                    (pow(det_aF, ta_epsilon) + pow(det_aF, -ta_epsilon) - 2.) -
                (DIM * DIM) * (ta_alpha + ta_beta);
 
    // Set dependent variables to ADOL-C
    ta_P(i, I) >>= t_P(i, I);
    ta_energy >>= (*energy_ptr)(0, 0);
 
    trace_off();
 
    MoFEMFunctionReturn(0);
  }
 
protected:
  double alpha = 1;
  double beta = 1;
  double lambda = 1;
  double epsilon = 2;
  std::vector<double> defaultMaterialParameters = {alpha, beta, lambda,
                                                   epsilon};
};
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<ELASTICITY::MOONEYRIVLINWRIGGERSEQ63,
                                 MODEL_3D>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag) {
  return boost::make_shared<MatMooneyRivlinWriggersEq63<3>>(mat_ops_data_ptr,
                                                            tag);
}
 
} // namespace MatOps