MatHuHu.hpp

Go to the documentation of this file.

/**
 * @file MatHuHu.hpp
 * @author your name (you@domain.com)
 * @brief 
 * @version 0.1
 * @date 2026-03-29
 * 
 * @copyright Copyright (c) 2026
 * 
 */
 
#ifndef MAT_HUHU_HPP
#define MAT_HUHU_HPP
 
namespace MatOps {
 
struct HUHU {
  HUHU() = delete;
};
 
template <int FIELD_DIM, int SPACE_DIM, AssemblyType A, IntegrationType I,
          typename OpBase>
struct OpRhsHu;
 
template <int FIELD_DIM, int SPACE_DIM, AssemblyType A, IntegrationType I,
          typename OpBase>
struct OpLhsHuHu;
 
template <int FIELD_DIM, int SPACE_DIM, AssemblyType A, IntegrationType I,
          typename OpBase>
struct OpLhsHuGrad;
 
using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
 
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, false));
    auto m_k = physical_equations_ptr->matOpsDataPtr->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, false));
    auto m_k = physical_equations_ptr->matOpsDataPtr->getCommonDataPtr("k");
    auto m_diff_k =
        physical_equations_ptr->matOpsDataPtr->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->matOpsDataPtr->getCommonDataPtr("grad");
    this_pip.push_back(
        new OpCalculateVectorFieldGradient<DIM, DIM>(field_name, m_grad));
 
    return op_this;
  }
};
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<HUHU, MODEL_3D>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag);
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<HUHU, MODEL_2D_PLANE_STRAIN>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag);
 
template <int FIELD_DIM, int SPACE_DIM, AssemblyType A, typename OpBase>
struct OpRhsHu<FIELD_DIM, SPACE_DIM, A, GAUSS, OpBase>
    : public FormsIntegrators<OpBase>::template Assembly<A>::OpBase {
 
  using OP = typename FormsIntegrators<OpBase>::template Assembly<A>::OpBase;
 
  OpRhsHu(const std::string field_name,
          boost::shared_ptr<MatrixDouble> mat_vals,
          boost::shared_ptr<MatrixDouble> mat_K,
          boost::shared_ptr<Range> ents_ptr = nullptr)
      : OP(field_name, field_name, OP::OPROW), matVals(mat_vals), matK(mat_K) {}
 
protected:
  boost::shared_ptr<MatrixDouble> matVals;
  boost::shared_ptr<MatrixDouble> matK;
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
};
 
template <int FIELD_DIM, int SPACE_DIM, AssemblyType A, typename OpBase>
struct OpLhsHuHu<FIELD_DIM, SPACE_DIM, A, GAUSS, OpBase>
    : public FormsIntegrators<OpBase>::template Assembly<A>::OpBase {
 
  using OP = typename FormsIntegrators<OpBase>::template Assembly<A>::OpBase;
  OpLhsHuHu(const std::string field_name, boost::shared_ptr<MatrixDouble> mat_K,
            boost::shared_ptr<Range> ents_ptr = nullptr)
      : OP(field_name, field_name, OP::OPROWCOL), matK(mat_K) {}
 
protected:
  boost::shared_ptr<MatrixDouble> matK;
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data);
};
 
template <int FIELD_DIM, int SPACE_DIM, AssemblyType A, typename OpBase>
struct OpLhsHuGrad<FIELD_DIM, SPACE_DIM, A, GAUSS, OpBase>
    : public FormsIntegrators<OpBase>::template Assembly<A>::OpBase {
  using OP = typename FormsIntegrators<OpBase>::template Assembly<A>::OpBase;
  OpLhsHuGrad(const std::string field_name,
              boost::shared_ptr<MatrixDouble> mat_vals,
              boost::shared_ptr<MatrixDouble> mat_diff_K,
              boost::shared_ptr<Range> ents_ptr = nullptr)
      : OP(field_name, field_name, OP::OPROWCOL), matVals(mat_vals),
        matDiffK(mat_diff_K) {
    OP::sYmm = false;
  }
 
protected:
  boost::shared_ptr<MatrixDouble> matVals;
  boost::shared_ptr<MatrixDouble> matDiffK;
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data);
};
 
template <int FIELD_DIM, int SPACE_DIM, AssemblyType A, typename OpBase>
MoFEMErrorCode OpRhsHu<FIELD_DIM, SPACE_DIM, A, GAUSS, OpBase>::iNtegrate(
    EntitiesFieldData::EntData &row_data) {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEX(FIELD_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, J);
  FTENSOR_INDEX(SPACE_DIM, K);
 
  auto get_val_grad_at_pts = MatrixSizeHelper<
      GetFTensor3FromMatType<FIELD_DIM, SPACE_DIM, SPACE_DIM, -1, DL>,
      DL>::get(*matVals, OP::nbIntegrationPts);
  auto get_K_at_pts =
      MatrixSizeHelper<GetFTensor1FromMatType<1, -1, DL>, DL>::get(
          *matK, OP::nbIntegrationPts);
 
  // get element volume
  const double vol = OP::getMeasure();
  // get integration weights
  auto t_w = OP::getFTensor0IntegrationWeight();
  // get base function gradient on rows
  auto t_row_grad = row_data.getFTensor2DiffN2<SPACE_DIM>();
  // get field gradient values
  auto t_val_grad_at_pts =
      get_val_grad_at_pts(); // tensor of second derivatives
  // get K values
  auto t_K_at_pts = get_K_at_pts(); // tensor of second derivatives
  // get coordinate at integration points
  auto t_coords = OP::getFTensor1CoordsAtGaussPts();
  // loop over integration points
  for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
    // take into account Jacobian
    const double alpha = t_w * vol * t_K_at_pts(0);
    // calculate rhs
    auto t_nf = getFTensor1FromArray<FIELD_DIM, FIELD_DIM>(OP::locF);
    // loop over rows base functions
    int rr = 0;
    for (; rr != OP::nbRows / FIELD_DIM; rr++) {
      // calculate element of local rhs vector
      t_nf(i) += alpha * (t_row_grad(J, K) * t_val_grad_at_pts(i, J, K));
      ++t_row_grad; // move to another element of gradient of base
                    // function on row
    }
    for (; rr < OP::nbRowBaseFunctions; ++rr)
      ++t_row_grad;
 
    ++t_coords;
    ++t_val_grad_at_pts;
    ++t_K_at_pts;
    ++t_w; // move to another integration weight
  }
  MoFEMFunctionReturn(0);
}
 
template <int FIELD_DIM, int SPACE_DIM, AssemblyType A, typename OpBase>
MoFEMErrorCode OpLhsHuHu<FIELD_DIM, SPACE_DIM, A, GAUSS, OpBase>::iNtegrate(
    EntitiesFieldData::EntData &row_data,
    EntitiesFieldData::EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEX(FIELD_DIM, i);
  FTENSOR_INDEX(FIELD_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, J);
  FTENSOR_INDEX(SPACE_DIM, K);
 
  auto get_K_at_pts =
      MatrixSizeHelper<GetFTensor1FromMatType<1, -1, DL>, DL>::get(
          *matK, OP::nbIntegrationPts);
 
  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
 
  // get element volume
  const double vol = OP::getMeasure();
  // get integration weights
  auto t_w = OP::getFTensor0IntegrationWeight();
  // get base function gradient on rows
  auto t_row_grad = row_data.getFTensor2DiffN2<SPACE_DIM>();
  // get k values
  auto t_K_at_pts = get_K_at_pts(); // tensor of second derivatives
  // get coordinate at integration points
  auto t_coords = OP::getFTensor1CoordsAtGaussPts();
  // loop over integration points
  for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
    // take into account Jacobian
    const double alpha = t_w * vol * t_K_at_pts(0);
    // loop over rows base functions
    int rr = 0;
    for (; rr != OP::nbRows / FIELD_DIM; rr++) {
 
      auto t_mat = getFTensor2FromArray<FIELD_DIM, FIELD_DIM, FIELD_DIM>(
          OP::locMat, rr * FIELD_DIM, 0);
      auto t_col_grad = col_data.getFTensor2DiffN2<SPACE_DIM>(gg, 0);
 
      // calculate element of local matrix
      for (int cc = 0; cc != OP::nbCols / FIELD_DIM; cc++) {
        t_mat(i, j) +=
            alpha * (t_row_grad(J, K) * t_col_grad(J, K)) * t_kd(i, j);
        ++t_col_grad; // move to another element of gradient of base
                      // function on column
        ++t_mat;      // move to another element of local matrix
      }
      ++t_row_grad; // move to another element of gradient of base
                    // function on row
    }
    for (; rr < OP::nbRowBaseFunctions; ++rr)
      ++t_row_grad;
 
    ++t_coords;
    ++t_K_at_pts;
    ++t_w; // move to another integration weight
  }
  MoFEMFunctionReturn(0);
}
 
template <int FIELD_DIM, int SPACE_DIM, AssemblyType A, typename OpBase>
MoFEMErrorCode OpLhsHuGrad<FIELD_DIM, SPACE_DIM, A, GAUSS, OpBase>::iNtegrate(
    EntitiesFieldData::EntData &row_data,
    EntitiesFieldData::EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEX(FIELD_DIM, i);
  FTENSOR_INDEX(FIELD_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, J);
  FTENSOR_INDEX(SPACE_DIM, K);
  FTENSOR_INDEX(SPACE_DIM, M);
 
  auto get_val_grad_at_pts = MatrixSizeHelper<
      GetFTensor3FromMatType<FIELD_DIM, SPACE_DIM, SPACE_DIM, -1, DL>,
      DL>::get(*matVals, OP::nbIntegrationPts);
  auto get_diff_K_at_pts =
      MatrixSizeHelper<GetFTensor2FromMatType<FIELD_DIM, SPACE_DIM, -1, DL>,
                       DL>::get(*matDiffK, OP::nbIntegrationPts);
 
  // get element volume
  const double vol = OP::getMeasure();
  // get integration weights
  auto t_w = OP::getFTensor0IntegrationWeight();
  // get base function gradient on rows
  auto t_row_grad = row_data.getFTensor2DiffN2<SPACE_DIM>();
  // get field gradient values
  auto t_val_grad_at_pts =
      get_val_grad_at_pts(); // tensor of second derivatives
  // get K values
  auto t_diff_K_at_pts = get_diff_K_at_pts(); // tensor of second derivatives
  // get coordinate at integration points
  auto t_coords = OP::getFTensor1CoordsAtGaussPts();
  // loop over integration points
  for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
    // take into account Jacobian
    const double alpha = t_w * vol;
    // loop over rows base functions
    int rr = 0;
    for (; rr != OP::nbRows / FIELD_DIM; rr++) {
 
      auto t_mat = getFTensor2FromArray<FIELD_DIM, FIELD_DIM, FIELD_DIM>(
          OP::locMat, rr * FIELD_DIM, 0);
      auto t_col_grad = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
      for (int bb = 0; bb != OP::nbCols / FIELD_DIM; ++bb) {
        t_mat(i, j) += alpha * (t_row_grad(J, K) * t_val_grad_at_pts(i, J, K)) *
                       (t_diff_K_at_pts(j, M) * t_col_grad(M));
 
        ++t_mat;      // move to another element of local matrix
        ++t_col_grad; // move to another element of gradient of base
                      // function on column
      }
 
      ++t_row_grad; // move to another element of gradient of base
                    // function on row
    }
    for (; rr < OP::nbRowBaseFunctions; ++rr)
      ++t_row_grad;
 
    ++t_coords;
    ++t_val_grad_at_pts;
    ++t_diff_K_at_pts;
    ++t_w; // move to another integration weight
  }
  MoFEMFunctionReturn(0);
}
 
} // namespace MatOps
 
#endif // MAT_HUHU_HPP