LoopMethods.hpp

Go to the documentation of this file.

/** \file LoopMethods.hpp
 * \brief MoFEM loop methods interface for finite element computations
 * \author Anonymous author(s) committing under MIT license 
 * 
 * This file defines data structures and interfaces for making loops over finite 
 * elements, entities, and degrees of freedom (DOFs) in MoFEM problems and database.
 *
 * \section loop_overview Loop Methods Overview
 * 
 * The core loop methods enable efficient iteration over:
 * - **Finite Elements**: Using FEMethod with Interface::loop_finite_elements()
 * - **Entities**: Using EntityMethod with Interface::loop_entities()  
 * - **DOFs**: Using DofMethod with Interface::loop_dofs()
 *
 * \section femethod_usage FEMethod Usage in Core Interface
 * 
 * FEMethod is the primary interface for finite element computations, designed for
 * high-performance element-level operations. It integrates with MoFEM's
 * core looping mechanisms:
 *
 * **Loop Integration:**
 * - Called via Interface::loop_finite_elements() for each element in a mesh
 * - User overloads FEMethod::operator() to define element-specific computations
 * - Provides access to element geometry, DOF data, and field values
 * - Supports assembly of global matrices/vectors from element contributions
 *
 * **PETSc Solver Integration:**
 * - **KSP (Linear)**: FEMethod assembles stiffness matrices and load vectors
 * - **SNES (Nonlinear)**: Computes residuals and Jacobians for Newton methods
 * - **TS (Time-stepping)**: Handles time-dependent mass matrices and residuals
 * - **TAO (Optimization)**: Evaluates objective functions and gradients
 *
 * **ForcesAndSources Integration:**
 * - FEMethod works with ForcesAndSources user data operators for element computations
 * - ForcesAndSources provides high-level element assembly and integration routines
 * - FEMethod offers direct low-level access for performance-critical applications
 *
 * Each solver context provides appropriate data structures (matrices, vectors,
 * contexts) accessible through the inheritance hierarchy (BasicMethod -> FEMethod).
 *
 * \section performance_design Performance Design
 * 
 * The loop methods are optimized for computational efficiency:
 * - Direct memory access to DOF data and field values
 * - Minimal overhead iteration over large meshes
 * - Cache-friendly data access patterns
 * - Integration with PETSc's high-performance linear algebra
 *
 * \see Interface::loop_finite_elements, Interface::loop_entities, Interface::loop_dofs
 *
 */
 
#ifndef __LOOPMETHODS_HPP__
#define __LOOPMETHODS_HPP__
 
namespace MoFEM {
 
/**
 * \brief Base data structure for PETSc-related contexts
 * \ingroup mofem_loops
 * 
 * This structure provides a foundation for managing PETSc data objects and contexts
 * used throughout MoFEM finite element computations. It handles vectors, matrices,
 * and context information needed for various PETSc solvers (KSP, SNES, TS, TAO).
 */
struct PetscData : public UnknownInterface {
 
  /**
   * \brief Query interface for type casting
   * \param type_index Type information for interface querying
   * \param iface Pointer to interface object
   * \return Error code
   */
  MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
                                 UnknownInterface **iface) const;
 
  /**
   * \brief Default constructor
   */
  PetscData();
 
  /**
   * \brief Virtual destructor
   */
  virtual ~PetscData() = default;
 
  /**
   * \brief Enumeration for data context flags
   * 
   * These flags indicate which PETSc data structures are currently set
   * and available for use in computations. Multiple flags can be combined
   * using bitwise operations.
   */
  enum DataContext {
    CTX_SET_NONE = 0,        ///< No data set
    CTX_SET_F = 1 << 0,      ///< Residual vector F is set
    CTX_SET_A = 1 << 1,      ///< Jacobian matrix A is set
    CTX_SET_B = 1 << 2,      ///< Preconditioner matrix B is set
    CTX_SET_X = 1 << 3,      ///< Solution vector X is set
    CTX_SET_DX = 1 << 4,     ///< Solution increment DX is set
    CTX_SET_X_T = 1 << 5,    ///< Time derivative X_t is set
    CTX_SET_X_TT = 1 << 6,   ///< Second time derivative X_tt is set
    CTX_SET_TIME = 1 << 7    ///< Time value is set
  };
 
  using Switches = std::bitset<8>; ///< Bitset type for context switches
 
  static constexpr Switches CtxSetNone = PetscData::Switches(CTX_SET_NONE);     ///< No data switch
  static constexpr Switches CtxSetF = PetscData::Switches(CTX_SET_F);           ///< Residual vector switch
  static constexpr Switches CtxSetA = PetscData::Switches(CTX_SET_A);           ///< Jacobian matrix switch
  static constexpr Switches CtxSetB = PetscData::Switches(CTX_SET_B);           ///< Preconditioner matrix switch
  static constexpr Switches CtxSetX = PetscData::Switches(CTX_SET_X);           ///< Solution vector switch
  static constexpr Switches CtxSetDX = PetscData::Switches(CTX_SET_DX);         ///< Solution increment switch
  static constexpr Switches CtxSetX_T = PetscData::Switches(CTX_SET_X_T);       ///< First time derivative switch
  static constexpr Switches CtxSetX_TT = PetscData::Switches(CTX_SET_X_TT);     ///< Second time derivative switch
  static constexpr Switches CtxSetTime = PetscData::Switches(CTX_SET_TIME);     ///< Time value switch
 
  Switches data_ctx; ///< Current data context switches
 
  /**
   * \brief Copy PETSc data from another instance
   * \param petsc_data Source PETSc data structure
   * \return Error code
   */
  MoFEMErrorCode copyPetscData(const PetscData &petsc_data);
 
  Vec f;    ///< PETSc residual vector
  Mat A;    ///< PETSc Jacobian matrix
  Mat B;    ///< PETSc preconditioner matrix
  Vec x;    ///< PETSc solution vector
  Vec dx;   ///< PETSc solution increment vector
  Vec x_t;  ///< PETSc first time derivative vector
  Vec x_tt; ///< PETSc second time derivative vector
};
 
/**
 * \brief Data structure for KSP (linear solver) context
 * \ingroup mofem_loops
 *
 * This structure extends PetscData to provide context and data management
 * specifically for PETSc KSP (Krylov Subspace) linear solvers. It stores
 * the KSP solver instance and maintains context information about the
 * current computation phase.
 */
struct KspMethod : virtual public PetscData {
 
  /**
   * \brief Query interface for type casting
   * \param type_index Type information for interface querying
   * \param iface Pointer to interface object
   * \return Error code
   */
  MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
                                 UnknownInterface **iface) const;
 
  /**
   * \brief Context enumeration for KSP solver phases
   * 
   * Indicates the current phase of KSP computation, which determines
   * which data structures are being used and what operations are valid.
   */
  enum KSPContext { 
    CTX_SETFUNCTION, ///< Setting up the linear system function
    CTX_OPERATORS,   ///< Setting up linear operators
    CTX_KSPNONE      ///< No specific KSP context
  };
 
  /**
   * \brief Default constructor
   */
  KspMethod();
 
  /**
   * \brief Virtual destructor
   */
  virtual ~KspMethod() = default;
 
  /**
   * \brief Copy data from another KSP method
   * \param ksp Source KSP method
   * \return Error code
   */
  MoFEMErrorCode copyKsp(const KspMethod &ksp);
 
  KSPContext ksp_ctx; ///< Current KSP computation context
  KSP ksp;            ///< PETSc KSP linear solver object
 
  Vec &ksp_f; ///< Reference to residual vector in KSP context
  Mat &ksp_A; ///< Reference to system matrix in KSP context
  Mat &ksp_B; ///< Reference to preconditioner matrix in KSP context
};
 
/**
 * \brief Data structure for SNES (nonlinear solver) context
 * \ingroup mofem_loops
 *
 * This structure extends PetscData to provide context and data management
 * specifically for PETSc SNES (Scalable Nonlinear Equation Solvers). It stores
 * the SNES solver instance and maintains context information about nonlinear
 * solution phases including function evaluation and Jacobian computation.
 */
struct SnesMethod : virtual protected PetscData {
 
  /**
   * \brief Query interface for type casting
   * \param type_index Type information for interface querying
   * \param iface Pointer to interface object
   * \return Error code
   */
  MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
                                 UnknownInterface **iface) const;
 
  /**
   * \brief Context enumeration for SNES solver phases
   * 
   * Indicates the current phase of SNES computation, determining which
   * operations are being performed and what data structures are active.
   */
  enum SNESContext { 
    CTX_SNESSETFUNCTION, ///< Setting up nonlinear function evaluation
    CTX_SNESSETJACOBIAN, ///< Setting up Jacobian matrix computation
    CTX_SNESNONE         ///< No specific SNES context
  };
 
  /**
   * \brief Default constructor
   */
  SnesMethod();
 
  /**
   * \brief Virtual destructor
   */
  virtual ~SnesMethod() = default;
 
  /**
   * \brief Copy SNES data from another instance
   * \param snes Source SNES method
   * \return Error code
   */
  MoFEMErrorCode copySnes(const SnesMethod &snes);
 
  SNESContext snes_ctx; ///< Current SNES computation context
 
  SNES snes;    ///< PETSc SNES nonlinear solver object
  Vec &snes_x;  ///< Reference to current solution state vector
  Vec &snes_dx; ///< Reference to solution update/increment vector
  Vec &snes_f;  ///< Reference to residual vector
  Mat &snes_A;  ///< Reference to Jacobian matrix
  Mat &snes_B;  ///< Reference to preconditioner of Jacobian matrix
};
 
/**
 * \brief Data structure for TS (time stepping) context
 * \ingroup mofem_loops
 *
 * This structure extends PetscData to provide context and data management
 * specifically for PETSc TS (Time Stepping) solvers. It manages time-dependent
 * computations including time derivatives, time step information, and temporal
 * integration contexts for transient finite element analyses.
 */
struct TSMethod : virtual protected PetscData {
 
  /**
   * \brief Query interface for type casting
   * \param type_index Type information for interface querying
   * \param iface Pointer to interface object
   * \return Error code
   */
  MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
                                 UnknownInterface **iface) const;
 
  /**
   * \brief Context enumeration for TS solver phases
   * 
   * Indicates the current phase of time stepping computation, determining
   * which temporal operations are being performed and what data is active.
   */
  enum TSContext {
    CTX_TSSETRHSFUNCTION, ///< Setting up right-hand side function
    CTX_TSSETRHSJACOBIAN, ///< Setting up right-hand side Jacobian
    CTX_TSSETIFUNCTION,   ///< Setting up implicit function
    CTX_TSSETIJACOBIAN,   ///< Setting up implicit Jacobian
    CTX_TSTSMONITORSET,   ///< Setting up time step monitoring
    CTX_TSNONE            ///< No specific TS context
  };
 
  /**
   * \brief Default constructor
   */
  TSMethod();
 
  /**
   * \brief Virtual destructor
   */
  virtual ~TSMethod() = default;
 
  /**
   * \brief Copy TS solver data from another instance
   * \param ts Source TS method
   * \return Error code
   */
  MoFEMErrorCode copyTs(const TSMethod &ts);
 
  TS ts; ///< PETSc time stepping solver object
 
  TSContext ts_ctx; ///< Current TS computation context
 
  PetscInt ts_step; ///< Current time step number
  PetscReal ts_a;   ///< Shift parameter for U_t (see PETSc Time Solver documentation)
  PetscReal ts_aa;  ///< Shift parameter for U_tt (second time derivative)
  PetscReal ts_t;   ///< Current time value
  PetscReal ts_dt;  ///< Current time step size
 
  Vec &ts_u;    ///< Reference to current state vector U(t)
  Vec &ts_u_t;  ///< Reference to first time derivative of state vector dU/dt
  Vec &ts_u_tt; ///< Reference to second time derivative of state vector d²U/dt²
  Vec &ts_F;    ///< Reference to residual vector F(t,U,U_t,U_tt)
 
  Mat &ts_A; ///< Reference to Jacobian matrix: dF/dU + a*dF/dU_t + aa*dF/dU_tt
  Mat &ts_B; ///< Reference to preconditioner matrix for ts_A
};
 
/**
 * \brief Data structure for TAO (optimization) context
 * \ingroup mofem_loops
 *
 * This structure extends PetscData to provide context and data management
 * specifically for PETSc TAO (Toolkit for Advanced Optimization) solvers.
 * It manages optimization computations including objective function evaluation,
 * gradient computation, and Hessian matrix assembly for constrained and
 * unconstrained optimization problems.
 */
struct TaoMethod : virtual protected PetscData {
 
  /**
   * \brief Query interface for type casting
   * \param type_index Type information for interface querying
   * \param iface Pointer to interface object
   * \return Error code
   */
  MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
                                 UnknownInterface **iface) const;
 
  /**
   * \brief Context enumeration for TAO solver phases
   * 
   * Indicates the current phase of optimization computation, determining
   * which optimization operations are being performed.
   */
  enum TAOContext {
    CTX_TAO_OBJECTIVE, ///< Evaluating objective function
    CTX_TAO_GRADIENT,  ///< Computing gradient vector
    CTX_TAO_HESSIAN,   ///< Assembling Hessian matrix
    CTX_TAO_NONE       ///< No specific TAO context
  };
 
  /**
   * \brief Default constructor
   */
  TaoMethod();
 
  /**
   * \brief Virtual destructor
   */
  virtual ~TaoMethod() = default;
 
  /**
   * \brief Copy TAO data from another instance
   * \param tao Source TAO method
   * \return Error code
   */
  MoFEMErrorCode copyTao(const TaoMethod &tao);
 
  TAOContext tao_ctx; ///< Current TAO computation context
 
  Tao tao;    ///< PETSc TAO optimization solver object
  Vec &tao_x; ///< Reference to optimization variables vector
  Vec &tao_f; ///< Reference to gradient vector
  Mat &tao_A; ///< Reference to Hessian matrix
  Mat &tao_B; ///< Reference to preconditioner matrix for Hessian
};
 
/**
 * \brief Data structure to exchange data between MoFEM and User Loop Methods
 * \ingroup mofem_loops
 *
 * This is a comprehensive base class that combines all PETSc solver contexts
 * (KSP, SNES, TS, TAO) and provides fundamental data exchange capabilities
 * between MoFEM and user-defined functions. It stores information about
 * multi-indices, problem databases, and provides access to finite element
 * data structures for loop-based computations.
 */
struct BasicMethod : public KspMethod, SnesMethod, TSMethod, TaoMethod {
 
  /**
   * \brief Query interface for type casting
   * \param type_index Type information for interface querying
   * \param iface Pointer to interface object
   * \return Error code
   */
  MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
                                 UnknownInterface **iface) const {
    MoFEMFunctionBeginHot;
    *iface = const_cast<BasicMethod *>(this);
    MoFEMFunctionReturnHot(0);
  }
 
  /**
   * \brief Default constructor
   */
  BasicMethod();
  
  /**
   * \brief Virtual destructor
   */
  virtual ~BasicMethod() = default;
 
  /**
   * @brief Current index of processed method in the loop
   * 
   * This counter tracks which method/element is currently being processed
   * in the loop iteration, useful for debugging and progress monitoring.
   */
  int nInTheLoop;
 
  /**
   * @brief Total number of methods to process in the loop
   * 
   * This value represents the total count of items (elements, entities, etc.)
   * that will be processed in the current loop execution.
   */
  int loopSize;
 
  /**
   * \brief Get current loop iteration index
   * \return Current position in the loop (0-based)
   */
  inline int getNinTheLoop() const { return nInTheLoop; }
 
  /**
   * \brief Get total loop size
   * \return Total number of items to process in the loop
   */
  inline int getLoopSize() const { return loopSize; }
 
  /**
   * @brief Processor rank range for distributed finite element iteration
   * 
   * This pair stores the low and high processor ranks that define the range
   * of processors involved in the current finite element iteration. Used for
   * parallel processing and load balancing across MPI processes.
   */
  std::pair<int, int> loHiFERank;
 
  /**
   * @brief Get processor rank range for finite element iteration
   * 
   * Returns the pair containing the lowest and highest processor ranks
   * involved in the current finite element loop iteration.
   * 
   * @return Pair<low_rank, high_rank> of processor ranges
   */
  inline auto getLoHiFERank() const { return loHiFERank; }
 
  /**
   * @brief Get lower processor rank for finite element iteration
   * 
   * Returns the lowest processor rank in the range of processors
   * that will process finite elements in the current loop.
   * 
   * @return Lowest processor rank in the iteration range
   */
  inline auto getLoFERank() const { return loHiFERank.first; }
 
  /**
   * @brief Get upper processor rank for finite element iteration
   * 
   * Returns the highest processor rank in the range of processors
   * that will process finite elements in the current loop.
   * 
   * @return Highest processor rank in the iteration range
   */
  inline auto getHiFERank() const { return loHiFERank.second; }
 
  int rAnk; ///< Current processor rank in MPI communicator
 
  int sIze; ///< Total number of processors in MPI communicator
 
  const RefEntity_multiIndex
      *refinedEntitiesPtr; ///< Pointer to container of refined MoFEM DOF entities
 
  const RefElement_multiIndex
      *refinedFiniteElementsPtr; ///< Pointer to container of refined finite element entities
 
  const Problem *problemPtr; ///< Raw pointer to current MoFEM problem instance
 
  const Field_multiIndex *fieldsPtr; ///< Raw pointer to fields multi-index container
 
  const FieldEntity_multiIndex
      *entitiesPtr; ///< Raw pointer to container of field entities
 
  const DofEntity_multiIndex *dofsPtr; ///< Raw pointer to container of degree of freedom entities
 
  const FiniteElement_multiIndex
      *finiteElementsPtr; ///< Raw pointer to container of finite elements
 
  const EntFiniteElement_multiIndex
      *finiteElementsEntitiesPtr; ///< Raw pointer to container of finite element entities
 
  const FieldEntityEntFiniteElementAdjacencyMap_multiIndex
      *adjacenciesPtr; ///< Raw pointer to container of adjacencies between DOFs and finite elements
 
  /**
   * \brief Get bit number for a specific field by name
   * 
   * This function looks up a field by name in the fields container and returns
   * its bit number, which is used for field identification and bit operations
   * in MoFEM's field management system.
   * 
   * \param field_name Name of the field to look up
   * \return Bit number of the field, or BITFEID_SIZE if field not found
   * \throws MoFEM exception if fieldsPtr is not set
   */
  inline unsigned int getFieldBitNumber(std::string field_name) const {
    if (fieldsPtr) {
      auto field_it = fieldsPtr->get<FieldName_mi_tag>().find(field_name);
      if (field_it != fieldsPtr->get<FieldName_mi_tag>().end())
        return (*field_it)->getBitNumber();
      else
        return BITFEID_SIZE;
    } else {
      THROW_MESSAGE("Pointer to fields multi-index is not set");
      return BITFEID_SIZE;
    }
  }
 
  /**
   * @brief Copy data from another BasicMethod instance
   *
   * This function copies all relevant data from another BasicMethod instance
   * including solver contexts, database pointers, and configuration settings.
   *
   * @param basic Source BasicMethod to copy from
   * @return Error code indicating success or failure
   */
  MoFEMErrorCode copyBasicMethod(const BasicMethod &basic);
 
  /**
   * @brief Hook function for pre-processing operations
   * 
   * User-defined function that is executed before the main loop begins.
   * Can be used for initialization, setup operations, or custom preprocessing.
   */
  boost::function<MoFEMErrorCode()> preProcessHook;
 
  /**
   * @brief Hook function for post-processing operations
   * 
   * User-defined function that is executed after the main loop completes.
   * Can be used for cleanup, finalization, or custom post-processing.
   */
  boost::function<MoFEMErrorCode()> postProcessHook;
 
  /**
   * @brief Hook function for main operator execution
   * 
   * User-defined function that can be executed during the main operator
   * call, providing additional customization points in the computation.
   */
  boost::function<MoFEMErrorCode()> operatorHook;
 
  /**
   * \brief Pre-processing function executed at loop initialization
   *
   * This virtual function is called once at the beginning of the loop.
   * It is typically used for:
   * - Zeroing matrices and vectors
   * - Computing shape functions on reference elements
   * - Preprocessing boundary conditions
   * - Setting up initial computation state
   * 
   * \return Error code indicating success or failure
   */
  virtual MoFEMErrorCode preProcess();
 
  /**
   * \brief Main operator function executed for each loop iteration
   *
   * This virtual function is called for every item (finite element, entity, etc.)
   * in the loop. It is the core computation function typically used for:
   * - Calculating element local matrices and vectors
   * - Assembling contributions to global system
   * - Element-level post-processing operations
   * 
   * \return Error code indicating success or failure
   */
  virtual MoFEMErrorCode operator()();
 
  /**
   * \brief Post-processing function executed at loop completion
   *
   * This virtual function is called once at the end of the loop.
   * It is typically used for:
   * - Assembling matrices and vectors to global system
   * - Computing global variables (e.g., total internal energy)
   * - Finalizing computation results
   * - Cleanup operations
   *
   * Example of iterating over DOFs:
   * \code
   * for(_IT_GET_FEROW_BY_NAME_DOFS_FOR_LOOP_(this,"DISPLACEMENT",it)) { 
   *   // Process each DOF
   * }
   * \endcode
   * 
   * \return Error code indicating success or failure
   */
  virtual MoFEMErrorCode postProcess();
 
  /**
   * @brief Get the cache weak pointer object
   *
   * Returns a weak pointer to the cache tuple that stores problem information
   * on entities about DOFs. Each problem stores different information.
   * 
   * \warning When iterating over finite elements in a TS (time-stepping) solver
   * preprocessor, use the TS cache in the loop. Using wrong cache will cause
   * undefined behavior or segmentation faults. This design is a necessary
   * compromise between bug resilience and memory/performance optimization.
   *
   * @return Weak pointer to the cache tuple
   */
  inline boost::weak_ptr<CacheTuple> getCacheWeakPtr() const {
    return cacheWeakPtr;
  }
 
  boost::movelib::unique_ptr<bool> vecAssembleSwitch; ///< Switch for vector assembly operations
  boost::movelib::unique_ptr<bool> matAssembleSwitch; ///< Switch for matrix assembly operations
 
  boost::weak_ptr<CacheTuple> cacheWeakPtr; ///< Weak pointer to cached entity data
};
 
/**
 * \brief Structure for user loop methods on finite elements
 * \ingroup mofem_loops
 *
 * This class provides a high-performance interface for finite element computations.
 * It can be used to calculate stiffness matrices, residual vectors, load vectors,
 * and other element-level quantities. While it is a low-level class, users seeking
 * maximum speed and efficiency can use it directly.
 *
 * This class is designed to work with Interface::loop_finite_elements, where the
 * user-overloaded operator FEMethod::operator() is executed for each element in
 * the problem. The class provides additional virtual methods (preProcess and
 * postProcess) that can be overloaded by users for custom preprocessing and
 * postprocessing operations.
 *
 * Key features:
 * - Direct access to finite element data structures
 * - High-performance element iteration
 * - Customizable preprocessing and postprocessing
 * - Integration with MoFEM's data management system
 */
struct FEMethod : public BasicMethod {
 
  /**
   * \brief Query interface for type casting
   * \param type_index Type information for interface querying  
   * \param iface Pointer to interface object
   * \return Error code
   */
  MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
                                 UnknownInterface **iface) const {
    MoFEMFunctionBeginHot;
    *iface = const_cast<FEMethod *>(this);
    MoFEMFunctionReturnHot(0);
  }
 
  /**
   * \brief Default constructor
   */
  FEMethod() = default;
 
  std::string feName; ///< Name of the finite element being processed
 
  boost::shared_ptr<const NumeredEntFiniteElement>
      numeredEntFiniteElementPtr; ///< Shared pointer to finite element database structure
 
  /**
   * @brief Get the name of the current finite element
   * 
   * Returns the name string identifier of the finite element currently
   * being processed in the loop.
   *
   * @return String containing the finite element name
   */
  inline auto getFEName() const { return feName; }
 
  /**
   * @brief Test function to determine if element should be skipped
   *
   * This optional hook function allows users to implement custom logic
   * to determine whether a particular finite element should be skipped
   * during the loop execution. If this function is set and returns false,
   * the element will be skipped in MoFEM::Core::loop_finite_elements.
   * 
   * \note This functionality is particularly useful for running elements
   * only on specific bit levels or under certain conditions.
   *
   * @param fe_method_ptr Pointer to the current FEMethod instance
   * @return true if element should be processed, false if it should be skipped
   */
  boost::function<bool(FEMethod *fe_method_ptr)> exeTestHook;
 
  /**
   * \brief Get pointer to DOF data for the current finite element
   * \return Pointer to DOF data container
   */
  inline auto getDataDofsPtr() const {
    return numeredEntFiniteElementPtr->getDataDofsPtr();
  };
 
  /**
   * \brief Get pointer to vector DOF data for the current finite element
   * \return Pointer to vector DOF data container
   */
  inline auto getDataVectorDofsPtr() const {
    return numeredEntFiniteElementPtr->getDataVectorDofsPtr();
  };
 
  /**
   * \brief Get reference to data field entities for the current finite element
   * \return Constant reference to field entities view
   */
  inline const FieldEntity_vector_view &getDataFieldEnts() const {
    return numeredEntFiniteElementPtr->getDataFieldEnts();
  }
 
  /**
   * \brief Get shared pointer to data field entities for the current finite element
   * \return Shared pointer to field entities view (mutable access)
   */
  inline boost::shared_ptr<FieldEntity_vector_view> &
  getDataFieldEntsPtr() const {
    return const_cast<NumeredEntFiniteElement *>(
               numeredEntFiniteElementPtr.get())
        ->getDataFieldEntsPtr();
  }
 
  /**
   * \brief Get reference to row field entities for the current finite element
   * \return Reference to row field entities container
   */
  inline auto &getRowFieldEnts() const {
    return numeredEntFiniteElementPtr->getRowFieldEnts();
  };
 
  /**
   * \brief Get shared pointer to row field entities for the current finite element
   * \return Shared pointer to row field entities container
   */
  inline auto &getRowFieldEntsPtr() const {
    return numeredEntFiniteElementPtr->getRowFieldEntsPtr();
  };
 
  /**
   * \brief Get reference to column field entities for the current finite element
   * \return Reference to column field entities container
   */
  inline auto &getColFieldEnts() const {
    return numeredEntFiniteElementPtr->getColFieldEnts();
  };
 
  /**
   * \brief Get shared pointer to column field entities for the current finite element
   * \return Shared pointer to column field entities container
   */
  inline auto &getColFieldEntsPtr() const {
    return numeredEntFiniteElementPtr->getColFieldEntsPtr();
  };
 
  /**
   * \brief Get pointer to row DOFs for the current finite element
   * 
   * Returns a pointer to the container of row degrees of freedom
   * associated with the current finite element. Row DOFs are used
   * in matrix assembly for the row indices.
   * 
   * \return Pointer to row DOFs container
   */
  inline auto getRowDofsPtr() const {
    return numeredEntFiniteElementPtr->getRowDofsPtr();
  };
 
  /**
   * \brief Get pointer to column DOFs for the current finite element
   * 
   * Returns a pointer to the container of column degrees of freedom
   * associated with the current finite element. Column DOFs are used
   * in matrix assembly for the column indices.
   * 
   * \return Pointer to column DOFs container
   */
  inline auto getColDofsPtr() const {
    return numeredEntFiniteElementPtr->getColDofsPtr();
  };
 
  /**
   * \brief Get the number of nodes in the current finite element
   * 
   * Returns the number of vertex nodes for the current finite element
   * based on its entity type (e.g., 4 for tetrahedron, 8 for hexahedron).
   * 
   * \return Number of nodes in the finite element
   */
  inline auto getNumberOfNodes() const;
 
  /**
   * \brief Get the entity handle of the current finite element
   * 
   * Returns the MOAB entity handle that uniquely identifies the
   * current finite element in the mesh database.
   * 
   * \return MOAB entity handle of the finite element
   */
  inline EntityHandle getFEEntityHandle() const;
 
  /**
   * \brief Get nodal data for a specific field
   * 
   * Extracts nodal values for the specified field from the current
   * finite element and stores them in the provided data vector.
   * 
   * \param field_name Name of the field to extract data for
   * \param data Vector to store the extracted nodal data
   * \param reset_dofs Flag to reset DOF values (default: true)
   * \return Error code indicating success or failure
   */
  MoFEMErrorCode getNodeData(const std::string field_name, VectorDouble &data,
                             const bool reset_dofs = true);
 
};
 
inline auto FEMethod::getNumberOfNodes() const {
  return moab::CN::VerticesPerEntity(numeredEntFiniteElementPtr->getEntType());
};
 
inline EntityHandle FEMethod::getFEEntityHandle() const {
  return numeredEntFiniteElementPtr->getEnt();
}
 
/**
 * \brief Data structure for user loop methods on entities
 * \ingroup mofem_loops
 *
 * This structure extends BasicMethod to provide entity-level loop operations.
 * It allows exchange of data between MoFEM and user functions when iterating
 * over mesh entities (vertices, edges, faces, volumes). It stores pointers
 * to field and entity information for entity-based computations.
 * 
 * Key features:
 * - Entity-level iteration capabilities
 * - Field data access for entities
 * - Integration with MoFEM's entity management system
 */
struct EntityMethod : public BasicMethod {
 
  /**
   * \brief Query interface for type casting
   * \param type_index Type information for interface querying
   * \param iface Pointer to interface object
   * \return Error code
   */
  MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
                                 UnknownInterface **iface) const {
    MoFEMFunctionBegin;
    *iface = const_cast<EntityMethod *>(this);
    MoFEMFunctionReturn(0);
  }
 
  /**
   * \brief Default constructor
   */
  EntityMethod() = default;
 
  boost::shared_ptr<Field> fieldPtr;       ///< Shared pointer to field information
  boost::shared_ptr<FieldEntity> entPtr;   ///< Shared pointer to field entity data
};
 
/**
 * \brief Data structure for user loop methods on degrees of freedom (DOFs)
 * \ingroup mofem_loops
 *
 * This structure extends BasicMethod to provide DOF-level loop operations.
 * It allows exchange of data between MoFEM and user functions when iterating
 * over degrees of freedom. It stores pointers to field, DOF entity, and
 * numbered DOF entity information for DOF-based computations.
 * 
 * Key features:
 * - DOF-level iteration capabilities
 * - Access to both raw and numbered DOF data
 * - Field information associated with DOFs
 * - Integration with MoFEM's DOF management system
 */
struct DofMethod : public BasicMethod {
 
  /**
   * \brief Query interface for type casting
   * \param type_index Type information for interface querying
   * \param iface Pointer to interface object
   * \return Error code
   */
  MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
                                 UnknownInterface **iface) const {
    MoFEMFunctionBeginHot;
    *iface = const_cast<DofMethod *>(this);
    MoFEMFunctionReturnHot(0);
  }
 
  /**
   * \brief Default constructor
   */
  DofMethod() = default;
 
  boost::shared_ptr<Field> fieldPtr;              ///< Shared pointer to field information
  boost::shared_ptr<DofEntity> dofPtr;            ///< Shared pointer to DOF entity data
  boost::shared_ptr<NumeredDofEntity> dofNumeredPtr; ///< Shared pointer to numbered DOF entity data
};
 
/// \deprecated name changed use DofMethod instead EntMethod
DEPRECATED typedef DofMethod EntMethod;
 
} // namespace MoFEM
 
#endif // __LOOPMETHODS_HPP__
 
/**
 * \defgroup mofem_loops Loops
 * \ingroup mofem
 */