FormsIntegrators.hpp

Go to the documentation of this file.

/** \file FormsIntegrators.hpp
  * \brief Forms integrators
  * \ingroup mofem_form
  
  This file provides a higher-level abstraction layer over MoFEM's generic 
  UserDataOperator system for finite element computations. It implements
  structured approaches to linear and bilinear form integration within
  finite element operator pipelines.
  
  \section forms_abstraction Forms as Operator Abstractions
  
  Forms integrators represent one level of abstraction above generic user data
  operators, providing structured interfaces for common finite element operations:
  
  - **LinearForm**: Operators that integrate over single fields to produce vectors
    (load vectors, source terms, boundary conditions)
  - **BiLinearForm**: Operators that integrate over field pairs to produce matrices  
    (stiffness matrices, mass matrices, coupling terms)
  - **TriLinearForm**: Operators for nonlinear terms involving multiple fields
  
  These forms encapsulate the integration logic while maintaining compatibility
  with MoFEM's finite element pipeline architecture.
  
  \section forms_pipeline Integration with Finite Element Pipelines
  
  Forms integrators work within MoFEM's operator pipeline system:
  
  1. **Finite Element Context**: Forms operate within specific finite elements
     that define geometry, connectivity, and field definitions
  2. **Operator Sequencing**: Forms are added to finite element operator pipelines
     alongside preprocessing, integration, and postprocessing operators
  3. **Data Flow**: Forms access field data, geometry, and integration points
     through the standard ForcesAndSourcesCore infrastructure
  4. **Assembly**: Computed local contributions are assembled into global
     algebraic structures (vectors, matrices)
  
  \section assembly_types Assembly Type Abstraction
  
  Assembly types control how local finite element contributions are assembled
  into global algebraic structures:
  
  - **PETSC**: Standard assembly into PETSc Vec/Mat objects
  - **SCHUR**: Assembly for Schur complement methods  
  - **BLOCK_MAT**: Assembly into block-structured matrices
  - **BLOCK_SCHUR**: Block-based Schur complement assembly
  - **BLOCK_PRECONDITIONER_SCHUR**: Specialized preconditioner assembly
  - **USER_ASSEMBLE**: Custom assembly controlled by user functions
  
  Assembly types are compile-time template parameters, enabling optimized
  code generation for different algebraic backends and solution strategies.
  
  \section integration_methods Integration Method Selection
  
  Integration methods determine how numerical integration is performed:
  
  - **GAUSS**: Gaussian quadrature integration using standard rules
  - **USER_INTEGRATION**: Custom integration schemes defined by user functions
  
  Integration methods are template parameters that work with assembly types
  to provide flexible, efficient integration strategies for different
  physics and discretization requirements.
  
  Note: Specialized methods can be introduced for specific cases, such as
  fast integration exploiting Bernstein–Bézier basis with Duffy transforms
  on simplex elements, or dedicated GPU backends for element-wise integration
  and assembly. These can be expressed as additional integration tags and
  paired with corresponding assembly strategies.
  
  \section forms_workflow Typical Workflow
  
  1. Select assembly type based on algebraic solution strategy
  2. Choose integration method appropriate for physics and accuracy requirements  
  3. Instantiate specific form type (Linear/BiLinear/TriLinear)
  4. Configure field operations and boundary conditions
  5. Add form operators to finite element pipeline
  6. Execute pipeline to assemble global algebraic system
  
  This abstraction enables systematic construction of finite element methods
  while maintaining performance and flexibility for advanced solution techniques.
 
*/
 
#ifndef __FORMS_INTEGRATORS_HPP__
#define __FORMS_INTEGRATORS_HPP__
 
namespace MoFEM {
 
//! [Storage and set boundary conditions]
 
struct EssentialBcStorage : public EntityStorage {
  EssentialBcStorage(VectorInt &indices) : entityIndices(indices) {}
  VectorInt entityIndices;
  /**
   * @brief Store modified indices by field
   *
   * Hash map, key is field name, value is storage.
   *
   */
  using HashVectorStorage =
      map<std::string, std::vector<boost::shared_ptr<EssentialBcStorage>>>;
  static HashVectorStorage feStorage;
};
 
/**
 * @brief Set indices on entities on finite element
 * @ingroup mofem_forms
 *
 * If index is marked, its value is set to -1. DOF with such index is not
 * assembled into the system.
 *
 * Indices are stored on the entity.
 *
 */
struct OpSetBc : public ForcesAndSourcesCore::UserDataOperator {
  OpSetBc(std::string field_name, bool yes_set,
          boost::shared_ptr<std::vector<unsigned char>> boundary_marker);
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
 
public:
  bool yesSet;
  boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
};
 
/**
 * @brief Unset indices on entities on finite element
 * @ingroup mofem_forms
 *
 * Removes boundary condition markings from DOF indices, allowing them to be
 * assembled into the system. This operator reverses the effect of OpSetBc.
 *
 */
struct OpUnSetBc : public ForcesAndSourcesCore::UserDataOperator {
  /**
   * @brief Constructor
   * @param field_name Name of field to unset boundary conditions on
   */
  OpUnSetBc(std::string field_name);
  
  /**
   * @brief Perform unset operation on entity data
   * @param side Entity side number
   * @param type Entity type
   * @param data Entity field data
   * @return Error code
   */
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
};
 
/**
 * @brief Set values to vector in operator
 * @ingroup mofem_forms
 *
 * MoFEM::FieldEntity provides MoFEM::FieldEntity::getWeakStoragePtr() storage
 * function which allows to transfer data between FEs or operators processing
 * the same entities.
 *
 * When MoFEM::OpSetBc is pushed in weak storage indices taking in account
 * indices which are skip to take boundary conditions are stored. Those entities
 * are used by VecSetValues.
 *
 * @param V
 * @param data
 * @param ptr
 * @param iora
 * @return MoFEMErrorCode
 */
template <>
MoFEMErrorCode
VecSetValues<EssentialBcStorage>(Vec V, const EntitiesFieldData::EntData &data,
                                 const double *ptr, InsertMode iora);
 
/**
 * @brief Set values to matrix in operator
 *
 * See MoFEM::VecSetValues<EssentialBcStorage> for explanation.
 *
 * @param M
 * @param row_data
 * @param col_data
 * @param ptr
 * @param iora
 * @return MoFEMErrorCode
 */
template <>
MoFEMErrorCode
MatSetValues<EssentialBcStorage>(Mat M,
                                 const EntitiesFieldData::EntData &row_data,
                                 const EntitiesFieldData::EntData &col_data,
                                 const double *ptr, InsertMode iora);
 
//! [Storage and set boundary conditions]
 
/**
 * @brief Form integrator assembly types
 * @ingroup mofem_forms
 *
 */
enum AssemblyType {
  PETSC,                       ///< Standard PETSc assembly
  SCHUR,                       ///< Schur complement assembly
  BLOCK_MAT,                   ///< Block matrix assembly
  BLOCK_SCHUR,                 ///< Block Schur assembly
  BLOCK_PRECONDITIONER_SCHUR,  ///< Block preconditioner Schur assembly
  USER_ASSEMBLE,               ///< User-defined assembly
  LAST_ASSEMBLE                ///< Sentinel value
};
 
/**
 * @brief Template selector for assembly type
 * @ingroup mofem_forms
 * @tparam A Assembly type value
 */
template <int A> struct AssemblyTypeSelector {};
 
template <>
inline MoFEMErrorCode MatSetValues<AssemblyTypeSelector<PETSC>>(
    Mat M, const EntitiesFieldData::EntData &row_data,
    const EntitiesFieldData::EntData &col_data, const double *ptr,
    InsertMode iora) {
  return MatSetValues<EssentialBcStorage>(M, row_data, col_data, ptr, iora);
}
 
template <>
inline MoFEMErrorCode VecSetValues<AssemblyTypeSelector<PETSC>>(
    Vec V, const EntitiesFieldData::EntData &data, const double *ptr,
    InsertMode iora) {
  return VecSetValues<EssentialBcStorage>(V, data, ptr, iora);
}
 
/**
 * @brief Form integrator integration types
 * @ingroup mofem_forms
 *
 */
enum IntegrationType { 
  GAUSS,            ///< Gaussian quadrature integration
  USER_INTEGRATION, ///< User-defined integration rules
  LAST_INTEGRATION  ///< Sentinel value
};
 
/**
 * @brief Scalar function type
 * @ingroup mofem_forms
 *
 */
using ScalarFun =
    boost::function<double(const double, const double, const double)>;
 
/**
 * @brief Default scalar function returning 1
 * @ingroup mofem_forms
 *
 * @param x coordinate
 * @param y coordinate  
 * @param z coordinate
 * @return 1.0
 */
inline double scalar_fun_one(const double, const double, const double) {
  return 1;
}
 
/**
 * @brief Lambda function used to scale with time
 *
 */
using TimeFun = boost::function<double(double)>;
 
/**
 * @brief Lambda function used to scale with time
 *
 */
using FEFun = boost::function<double(const FEMethod *fe_ptr)>;
 
/**
 * @brief Constant function type
 *
 */
using ConstantFun = boost::function<double()>;
 
/**
 * @brief Vector function type
 * @ingroup mofem_forms
 *
 * @tparam DIM dimension of the return
 */
template <int DIM>
using VectorFun = boost::function<FTensor::Tensor1<double, DIM>(
    const double, const double, const double)>;
 
template <AssemblyType A, typename EleOp> struct OpBaseImpl : public EleOp {
  using OpType = typename EleOp::OpType;
  using EntData = EntitiesFieldData::EntData;
 
  /**
   * @brief Constructor for base operator implementation
   * @param row_field_name Name of field for row entities
   * @param col_field_name Name of field for column entities
   * @param type Operator type
   * @param ents_ptr Optional pointer to range of entities to operate on
   */
  OpBaseImpl(const std::string row_field_name, const std::string col_field_name,
             const OpType type, boost::shared_ptr<Range> ents_ptr = nullptr)
      : EleOp(row_field_name, col_field_name, type, false),
        assembleTranspose(false), onlyTranspose(false), entsPtr(ents_ptr) {}
 
  /**
   * \brief Do calculations for the left hand side
   * @param  row_side row side number (local number) of entity on element
   * @param  col_side column side number (local number) of entity on element
   * @param  row_type type of row entity
   * @param  col_type type of column entity
   * @param  row_data data for row
   * @param  col_data data for column
   * @return          error code
   */
  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                        EntityType col_type, EntData &row_data,
                        EntData &col_data);
 
  /**
   * @brief Do calculations for the right hand side
   *
   * @param row_side
   * @param row_type
   * @param row_data
   * @return MoFEMErrorCode
   */
  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
 
  TimeFun timeScalingFun;           ///< assumes that time variable is set
  FEFun feScalingFun;               ///< set by fe entity handle
  boost::shared_ptr<Range> entsPtr; ///< Entities on which element is run
 
  using MatSetValuesHook = boost::function<MoFEMErrorCode(
      ForcesAndSourcesCore::UserDataOperator *op_ptr,
      const EntitiesFieldData::EntData &row_data,
      const EntitiesFieldData::EntData &col_data, MatrixDouble &m)>;
 
  static MatSetValuesHook matSetValuesHook;
 
protected:
  using EleOp::EleOp;
 
  /**
   * @brief Get local vector tensor for assembly
   * @tparam DIM Dimension of the tensor
   * @return FTensor wrapper around local vector
   */
  template <int DIM>
  inline auto getNf() {
    return getFTensor1FromArray<DIM, DIM>(locF);
  }
 
  /**
   * @brief Get local matrix tensor for assembly
   * @tparam DIM Dimension of the tensor
   * @param rr Row offset in the local matrix
   * @return FTensor wrapper around local matrix
   */
  template <int DIM>
  inline FTensor::Tensor2<FTensor::PackPtr<double *, DIM>, DIM, DIM>
  getLocMat(const int rr) {
    return getFTensor2FromArray<DIM, DIM, DIM>(locMat, rr);
  }
 
  int nbRows;             ///< number of dofs on rows
  int nbCols;             ///< number if dof on column
  int nbIntegrationPts;   ///< number of integration points
  int nbRowBaseFunctions; ///< number or row base functions
 
  int rowSide;        ///< row side number
  int colSide;        ///< column side number
  EntityType rowType; ///< row type
  EntityType colType; ///< column type
 
  bool assembleTranspose;
  bool onlyTranspose;
 
  MatrixDouble locMat;          ///< local entity block matrix
  MatrixDouble locMatTranspose; ///< local entity block matrix
  VectorDouble locF;            ///< local entity vector
 
  /**
   * \brief Integrate grad-grad operator
   * @param  row_data row data (consist base functions on row entity)
   * @param  col_data column data (consist base functions on column entity)
   * @return          error code
   */
  virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data) {
    return MOFEM_NOT_IMPLEMENTED;
  }
 
  /**
   * @brief Assemble local matrix into global matrix
   * @param row_data Row entity data
   * @param col_data Column entity data  
   * @param trans Whether to assemble transpose
   * @return Error code
   */
  virtual MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data,
                                  const bool trans);
 
  /**
   * \brief Class dedicated to integrate operator
   * @param  data entity data on element row
   * @return      error code
   */
  virtual MoFEMErrorCode iNtegrate(EntData &data) {
    return MOFEM_NOT_IMPLEMENTED;
  }
 
  /**
   * @brief Assemble local vector into global vector
   * @param data Entity data
   * @return Error code
   */
  virtual MoFEMErrorCode aSsemble(EntData &data);
 
  /** \brief Get number of base functions
   *
   * @param data
   * @return number of base functions
   */
  virtual size_t getNbOfBaseFunctions(EntitiesFieldData::EntData &data);
};
 
template <AssemblyType A, typename EleOp>
typename OpBaseImpl<A, EleOp>::MatSetValuesHook
    OpBaseImpl<A, EleOp>::matSetValuesHook =
        [](ForcesAndSourcesCore::UserDataOperator *op_ptr,
           const EntitiesFieldData::EntData &row_data,
           const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
          return MatSetValues<AssemblyTypeSelector<A>>(
              op_ptr->getKSPB(), row_data, col_data, m, ADD_VALUES);
        };
 
template <typename OpBase>
struct OpBrokenBaseImpl; //< This is base operator for broken spaces.
                         //Implementation is in
                         //FormsBrokenSpaceConstraintImpl.hpp
 
/**
 * @brief Integrator forms
 * @ingroup mofem_forms
 *
 * @tparam EleOp
 */
template <typename EleOp> struct FormsIntegrators {
 
  using EntData = EntitiesFieldData::EntData;
  using OpType = typename EleOp::OpType;
 
  /**
   * @brief Assembly methods
   * @ingroup mofem_forms
   *
   * @tparam A
   */
  template <AssemblyType A> struct Assembly {
 
    using OpBase = OpBaseImpl<A, EleOp>;
    using OpBrokenBase = OpBrokenBaseImpl<OpBase>;
 
    /**
     * @brief Linear form
     * @ingroup mofem_forms
     *
     * @tparam I
     */
    template <IntegrationType I> struct LinearForm;
 
    /**
     * @brief Bi linear form
     * @ingroup mofem_forms
     *
     * @tparam I
     */
    template <IntegrationType I> struct BiLinearForm;
 
    /**
     * @brief Tri linear form
     * @ingroup mofem_forms
     *
     * @tparam I
     */
    template <IntegrationType I> struct TriLinearForm;
 
  }; // Assembly
};   // namespace MoFEM
 
template <AssemblyType A, typename EleOp>
size_t
OpBaseImpl<A, EleOp>::getNbOfBaseFunctions(EntitiesFieldData::EntData &data) {
  if (!data.getNSharedPtr(data.getBase())) {
  #ifndef NDEBUG
    MOFEM_LOG("SELF", Sev::warning)
        << "Ptr to base " << ApproximationBaseNames[data.getBase()]
        << " functions is null, returning 0";
  #endif // NDEBUG
    return 0;
  }
  auto nb_base_functions = data.getN().size2();
  if (data.getBase() != USER_BASE) {
    switch (data.getSpace()) {
    case NOSPACE:
      break;
    case NOFIELD:
      break;
    case H1:
      break;
    case HCURL:
    case HDIV:
      nb_base_functions /= 3;
#ifndef NDEBUG
      if (data.getN().size2() % 3) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Number of base functions is not divisible by 3");
      }
#endif
      break;
    case L2:
      break;
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
               "Space %s not implemented", FieldSpaceNames[data.getSpace()]);
    }
  }
  return nb_base_functions;
}
 
template <AssemblyType A, typename EleOp>
MoFEMErrorCode
OpBaseImpl<A, EleOp>::doWork(int row_side, int col_side, EntityType row_type,
                             EntityType col_type,
                             EntitiesFieldData::EntData &row_data,
                             EntitiesFieldData::EntData &col_data) {
  MoFEMFunctionBegin;
 
  if (entsPtr) {
    if (entsPtr->find(this->getFEEntityHandle()) == entsPtr->end())
      MoFEMFunctionReturnHot(0);
  }
 
  // get number of dofs on row
  nbRows = row_data.getIndices().size();
  // if no dofs on row, exit that work, nothing to do here
  if (!nbRows)
    MoFEMFunctionReturnHot(0);
  rowSide = row_side;
  rowType = row_type;
  // get number of dofs on column
  nbCols = col_data.getIndices().size();
  // if no dofs on column, exit nothing to do here
  if (!nbCols)
    MoFEMFunctionReturnHot(0);
  colSide = col_side;
  colType = col_type;
  // get number of integration points
  nbIntegrationPts = EleOp::getGaussPts().size2();
  // get row base functions
  nbRowBaseFunctions = getNbOfBaseFunctions(row_data);
 
  // set size of local entity bock
  locMat.resize(nbRows, nbCols, false);
  // clear matrix
  locMat.clear();
  // integrate local matrix for entity block
  CHKERR this->iNtegrate(row_data, col_data);
 
  // assemble local matrix
  auto check_if_assemble_transpose = [&] {
    if (this->sYmm) {
      if (row_side != col_side || row_type != col_type)
        return true;
      else
        return false;
    } else if (assembleTranspose) {
      return true;
    }
 
    return false;
  };
  CHKERR aSsemble(row_data, col_data, check_if_assemble_transpose());
  MoFEMFunctionReturn(0);
}
 
template <AssemblyType A, typename EleOp>
MoFEMErrorCode OpBaseImpl<A, EleOp>::doWork(int row_side, EntityType row_type,
                                            EntData &row_data) {
  MoFEMFunctionBegin;
 
  if (entsPtr) {
    if (entsPtr->find(this->getFEEntityHandle()) == entsPtr->end())
      MoFEMFunctionReturnHot(0);
  }
 
  // get number of dofs on row
  nbRows = row_data.getIndices().size();
  rowSide = row_side;
  rowType = row_type;
 
  if (!nbRows)
    MoFEMFunctionReturnHot(0);
  // get number of integration points
  nbIntegrationPts = EleOp::getGaussPts().size2();
  // get row base functions
  nbRowBaseFunctions = getNbOfBaseFunctions(row_data);
  // resize and clear the right hand side vector
  locF.resize(nbRows);
  locF.clear();
  // integrate local vector
  CHKERR this->iNtegrate(row_data);
  // assemble local vector
  CHKERR this->aSsemble(row_data);
  MoFEMFunctionReturn(0);
}
 
template <AssemblyType A, typename EleOp>
MoFEMErrorCode OpBaseImpl<A, EleOp>::aSsemble(EntData &row_data,
                                              EntData &col_data,
                                              const bool transpose) {
  MoFEMFunctionBegin;
 
  if (!this->timeScalingFun.empty())
    this->locMat *= this->timeScalingFun(this->getFEMethod()->ts_t);
  if (!this->feScalingFun.empty())
    this->locMat *= this->feScalingFun(this->getFEMethod());
 
  // Assemble transpose
  if (transpose) {
    this->locMatTranspose.resize(this->locMat.size2(), this->locMat.size1(),
                                 false);
    noalias(this->locMatTranspose) = trans(this->locMat);
    CHKERR matSetValuesHook(this, col_data, row_data, this->locMatTranspose);
  }
 
  if (!this->onlyTranspose) {
    // assemble local matrix
    CHKERR matSetValuesHook(this, row_data, col_data, this->locMat);
  }
 
  MoFEMFunctionReturn(0);
}
 
template <AssemblyType A, typename EleOp>
MoFEMErrorCode OpBaseImpl<A, EleOp>::aSsemble(EntData &data) {
  if (!this->timeScalingFun.empty())
    this->locF *= this->timeScalingFun(this->getFEMethod()->ts_t);
  if (!this->feScalingFun.empty())
    this->locF *= this->feScalingFun(this->getFEMethod());
 
  return VecSetValues<AssemblyTypeSelector<A>>(this->getKSPf(), data,
                                               this->locF, ADD_VALUES);
}
 
} // namespace MoFEM
 
/**
 * \defgroup mofem_forms Forms Integrators
 *
 * \brief Classes and functions used to evaluate fields at integration pts,
 *jacobians, etc..
 *
 * \ingroup mofem_forces_and_sources
 **/
 
#endif //__FORMS_INTEGRATORS_HPP__