AdolCEP.hpp

Go to the documentation of this file.

/**
 * @file AdolCStreach.hpp
 * @author your name (you@domain.com)
 * @brief Tested by stretch of Log of streach
 * @date 2026-04-04
 * 
 * @copyright Copyright (c) 2026
 * 
 */
 
#ifndef ADOLC_OPS_STRETCH_HPP
#define ADOLC_OPS_STRETCH_HPP
 
using namespace MoFEM;
using namespace AdolCOps;
 
namespace EshelbianPlasticity {
namespace AdolCOps {    
 
struct EP {
  EP() = delete;
  // Behaviors
  struct META {};
  struct NEOHOOKEAN {};
  struct KIRCHHOFF_SAINT_VENANT {};
  struct VOLUMELENGTHQUALITY {};
};
 
template <int FIELD_DIM, int SPACE_DIM, AssemblyType A, IntegrationType I,
          typename OpBase>
struct OpRhsU;
 
template <int FIELD_DIM, int SPACE_DIM, AssemblyType A, IntegrationType I,
          typename OpBase>
struct OpLhsUU;
 
template <int FIELD_DIM, int SPACE_DIM, AssemblyType A, IntegrationType I,
          typename OpBase>
struct OpLhsUR;
 
template <int FIELD_DIM, int SPACE_DIM, AssemblyType A, IntegrationType I,
          typename OpBase>
struct OpLhsUP;
 
template <int MODEL_TYPE> struct OpMaterialFactory<HUHU, MODEL_TYPE> {
  OpMaterialFactory() = delete;
 
  template <AssemblyType A, IntegrationType I, typename DomainEleOp>
  static MoFEMErrorCode
  opRhsFactory(MoFEM::Interface &m_field,
               boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
               std::string fe_name, std::string field_name,
               boost::shared_ptr<PhysicalEquations> physical_equations_ptr,
               Sev sev = Sev::noisy) {
    MoFEMFunctionBegin;
    constexpr int DIM = (MODEL_TYPE == MODEL_3D) ? 3 : 2;
    auto op_this = getPipThis<DomainEleOp>(m_field, pip, fe_name, field_name,
                                           physical_equations_ptr, sev);
    auto m_grad_grad = boost::make_shared<MatrixDouble>();
    op_this->getOpPtrVector().push_back(
        new OpCalculateVectorFieldHessian<DIM, DIM>(field_name, m_grad_grad));
    op_this->getOpPtrVector().push_back(
        physical_equations_ptr->createOp(physical_equations_ptr, true, false));
    auto m_k = physical_equations_ptr->adolcDataPtr->getCommonDataPtr("k");
    op_this->getOpPtrVector().push_back(
        new OpRhsHu<DIM, DIM, A, I, DomainEleOp>(field_name, m_grad_grad, m_k));
    MoFEMFunctionReturn(0);
  }
 
  template <AssemblyType A, IntegrationType I, typename DomainEleOp>
  static MoFEMErrorCode
  opLhsFactory(MoFEM::Interface &m_field,
               boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
               std::string fe_name, std::string field_name,
               boost::shared_ptr<PhysicalEquations> physical_equations_ptr,
               Sev sev = Sev::noisy) {
    MoFEMFunctionBegin;
    constexpr int DIM = (MODEL_TYPE == MODEL_3D) ? 3 : 2;
    auto op_this = getPipThis<DomainEleOp>(m_field, pip, fe_name, field_name,
                                           physical_equations_ptr, sev);
    auto m_grad_grad = boost::make_shared<MatrixDouble>();
    op_this->getOpPtrVector().push_back(
        new OpCalculateVectorFieldHessian<DIM, DIM>(field_name, m_grad_grad));
    op_this->getOpPtrVector().push_back(
        physical_equations_ptr->createOp(physical_equations_ptr, true, true));
    auto m_k = physical_equations_ptr->adolcDataPtr->getCommonDataPtr("k");
    auto m_diff_k =
        physical_equations_ptr->adolcDataPtr->getCommonDataPtr("k_dF");
    op_this->getOpPtrVector().push_back(
        new OpLhsHuHu<DIM, DIM, A, I, DomainEleOp>(field_name, m_k));
    op_this->getOpPtrVector().push_back(
        new OpLhsHuGrad<DIM, DIM, A, I, DomainEleOp>(field_name, m_grad_grad,
                                                     m_diff_k));
    MoFEMFunctionReturn(0);
  }
 
private:
  template <typename DomainEleOp>
  static auto
  getPipThis(MoFEM::Interface &m_field,
             boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
             std::string fe_name, std::string field_name,
             boost::shared_ptr<PhysicalEquations> physical_equations_ptr,
             Sev sev) {
 
    auto &param_vec_by_range = physical_equations_ptr->paramVecByRange;
    // range only if parameters are set via blockset
    auto r = boost::make_shared<Range>();
    for (auto &p : param_vec_by_range) {
      r->merge(p.first);
    }
    MOFEM_LOG("WORLD", Sev::inform) << "HuHu number of entities " << r->size();
 
    constexpr int DIM = (MODEL_TYPE == MODEL_3D) ? 3 : 2;
 
    auto field_structure = m_field.get_field_structure(field_name);
    auto base = field_structure->getApproxBase();
    auto space = field_structure->getSpace();
 
    using DomainEle = typename parent_of<DomainEleOp>::type;
    auto op_this = new OpLoopThis<DomainEle>(m_field, fe_name, r, sev);
    pip.push_back(op_this);
    auto this_fe_ptr = op_this->getThisFEPtr();
    auto &this_pip = op_this->getOpPtrVector();
    this_fe_ptr->getRuleHook = [](int order_row, int order_col,
                                  int order_data) {
      return 2 * (order_data - 1);
    };
 
    auto base_mass = boost::make_shared<MatrixDouble>();
    auto data_l2 = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
    auto jac_ptr = boost::make_shared<MatrixDouble>();
    auto det_ptr = boost::make_shared<VectorDouble>();
    auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
 
    // calculate jacobian at integration points
    this_pip.push_back(new OpCalculateHOJac<DIM>(jac_ptr));
    // calculate jacobian at integration points
    this_pip.push_back(new OpInvertMatrix<DIM>(jac_ptr, det_ptr, inv_jac_ptr));
    // calculate mass matrix to project derivatives
    this_pip.push_back(
        new OpBaseDerivativesMass<1>(base_mass, data_l2, base, L2));
    // calculate second derivative of base functions, i.e. hessian
    this_pip.push_back(new OpBaseDerivativesNext<1>(
        BaseDerivatives::SecondDerivative, base_mass, data_l2, base, space));
 
    switch (space) {
    case H1:
      // push first base derivatives tp physical element shape
      this_pip.push_back(
          new OpSetHOInvJacToScalarBases<DIM, 1>(space, inv_jac_ptr));
      // push second base derivatives tp physical element shape
      this_pip.push_back(
          new OpSetHOInvJacToScalarBases<DIM, 2>(space, inv_jac_ptr));
      break;
    default:
      MOFEM_LOG("WORLD", Sev::error)
          << "Unsupported space: " << space << " for field: " << field_name;
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Unsupported space");
    }
 
    auto m_grad =
        physical_equations_ptr->adolcDataPtr->getCommonDataPtr("grad");
    this_pip.push_back(
        new OpCalculateVectorFieldGradient<DIM, DIM>(field_name, m_grad));
 
    return op_this;
  }
};
 
 
 
} // namespace AdolCOps
} // namespace EshelbianPlasticity
 
#endif // ADOLC_OPS_STRETCH_HPP