HODataOperators.hpp

Go to the documentation of this file.

/** \file HODataOperators.hpp
 * \brief High-order geometry operators for finite element computations
 *
 * This file provides operators for managing high-order (curved) geometry in
 * finite element analysis. These operators handle geometric transformations,
 * coordinate mappings, and basis function modifications required when elements
 * have curved boundaries or use high-order shape functions.
 *
 * ## Main Functionality:
 *
 * **Jacobian Calculations**: Compute transformation matrices between reference
 * and physical coordinates for curved elements (OpCalculateHOJacForVolume,
 * OpCalculateHOJacForFace, OpGetHOTangentsOnEdge, OpGetHONormalsOnFace)
 *
 * **Coordinate Transformations**: Calculate physical coordinates at integration
 * points using high-order approximations (OpCalculateHOCoords)
 *
 * **Basis Function Transformations**: Apply Piola transforms and inverse Jacobian
 * operations to basis functions for H1, Hdiv, and Hcurl spaces
 * (OpSetHOContravariantPiolaTransform, OpSetHOCovariantPiolaTransform)
 *
 * **Integration Weight Corrections**: Adjust integration weights to account for
 * geometric distortions in curved elements (OpSetHOWeights, OpSetHOWeightsOnFace)
 *
 * **Pipeline Management**: AddHOOps template specializations automatically add
 * appropriate HO operators to finite element pipelines based on element
 * dimension and problem type
 *
 * ## Usage Pattern:
 * 1. Use AddHOOps<FE_DIM, PROBLEM_DIM, SPACE_DIM>::add() to automatically
 *    add required HO operators to your finite element pipeline
 * 2. Individual operators can be used directly for specialized applications
 * 3. Operators handle both internal computations and external data access
 *    via shared pointers to Jacobian matrices and determinants
 *
 * ## Key Benefits:
 * - Accurate integration on curved domains
 * - Proper handling of geometric nonlinearity
 * - Seamless integration with MoFEM pipeline system
 * - Support for mixed finite element spaces (H1, Hdiv, Hcurl)
 *
 * @note Essential for problems involving curved boundaries, contact mechanics,
 *       fluid-structure interaction, and any application requiring geometric
 *       accuracy beyond linear approximation.
 */
 
#ifndef __HO_DATA_OPERATORS_HPP__
#define __HO_DATA_OPERATORS_HPP__
 
namespace MoFEM {
 
/**
 * @brief Calculate Jacobian on Hex or other volume elements which are not simplex
 *
 * This structure calculates the Jacobian matrix for high-order volume elements
 * using geometrical coordinates of the element rather than field data. It's
 * specifically designed for hexahedral and other non-simplex volume elements.
 *
 * @note Use OpCalculateVectorFieldGradient to calculate Jacobian from field data.
 */
struct OpCalculateHOJacForVolume
    : public VolumeElementForcesAndSourcesCore::UserDataOperator {
 
  /**
   * @brief Constructor for HO Jacobian calculation on volumes
   *
   * @param jac_ptr Shared pointer to matrix for storing calculated Jacobian values
   */
  OpCalculateHOJacForVolume(boost::shared_ptr<MatrixDouble> jac_ptr);
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
 
private:
  boost::shared_ptr<MatrixDouble> jacPtr;
};
 
/** \deprecated use OpCalculateHOJacForVolume
 */
using OpCalculateHOJacVolume = OpCalculateHOJacForVolume;
 
/**
 * @brief Calculate high-order coordinates at Gauss points
 *
 * This template structure computes high-order coordinate values at integration
 * (Gauss) points using field approximations. It supports different field
 * dimensions and is commonly used in high-order finite element methods.
 *
 * @tparam FIELD_DIM Dimension of the coordinate field (default: 3)
 */
template <int FIELD_DIM = 3>
struct OpCalculateHOCoords : public ForcesAndSourcesCore::UserDataOperator {
 
  /**
   * @brief Constructor for HO coordinates calculation
   *
   * @param field_name Name of the field containing coordinate information
   */
  OpCalculateHOCoords(const std::string field_name)
      : ForcesAndSourcesCore::UserDataOperator(field_name, OPROW) {}
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
};
 
/**
 * @deprecated This class should be DIM = 3 specialization when default
 * parameter is removed
 *
 */
struct OpSetHOInvJacToScalarBasesImpl : public OpSetInvJacToScalarBasesBasic {
 
  using OpSetInvJacToScalarBasesBasic::OpSetInvJacToScalarBasesBasic;
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
};
 
/**
 * \brief Set inverse jacobian to base functions
 *
 * \deprecated  Version with default variant DIM = 3 is deprecated, keep for
 * back compatibility. Should be removed in future. Use of default variant make
 * code implicit, what can be source of some hidden error. Explict interface is
 * better, when user see and control outcome, and is aware of existing variants.
 *
 */
template <int DIM = 3>
struct OpSetHOInvJacToScalarBases : public OpSetHOInvJacToScalarBasesImpl {
  using OpSetHOInvJacToScalarBasesImpl::OpSetHOInvJacToScalarBasesImpl;
};
 
template <>
struct OpSetHOInvJacToScalarBases<2>
    : public OpSetInvJacSpaceForFaceImpl<2, 1> {
 
  using OpSetInvJacSpaceForFaceImpl<2, 1>::OpSetInvJacSpaceForFaceImpl;
};
 
/**
 * @brief Transform local reference derivatives of shape functions to global derivatives
 *
 * This structure transforms derivatives of shape functions from local reference
 * coordinates to global coordinates when higher-order geometry is used. It's
 * essential for vector-based finite element spaces requiring coordinate transformations.
 *
 * @ingroup mofem_forces_and_sources
 */
struct OpSetHOInvJacVectorBase : public ForcesAndSourcesCore::UserDataOperator {
 
  /**
   * @brief Constructor for HO inverse Jacobian transformation for vector bases
   *
   * @param space Field space specification for the transformation
   * @param inv_jac_ptr Shared pointer to matrix containing inverse Jacobian data
   */
  OpSetHOInvJacVectorBase(const FieldSpace space,
                          boost::shared_ptr<MatrixDouble> inv_jac_ptr)
      : ForcesAndSourcesCore::UserDataOperator(space), invJacPtr(inv_jac_ptr) {    if (!invJacPtr)
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Pointer for invJacPtr not allocated");
 
    doVertices = false;
    if (space == HDIV)
      doEdges = false;
  }
 
    
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
 
private:
  boost::shared_ptr<MatrixDouble> invJacPtr;
  MatrixDouble diffHdivInvJac;
  MatrixDouble diffNinvJac;
};
 
/**
 * @brief Modify integration weights on face to take into account higher-order geometry
 *
 * This structure adjusts integration weights for face elements when higher-order
 * geometry is present. It accounts for geometric distortions introduced by
 * curved boundaries and high-order shape functions.
 *
 * @ingroup mofem_forces_and_sources_tri_element
 */
struct OpSetHOWeightsOnFace
    : public FaceElementForcesAndSourcesCore::UserDataOperator {
 
  /**
   * @brief Default constructor for HO weights on face
   */
  OpSetHOWeightsOnFace()
      : FaceElementForcesAndSourcesCore::UserDataOperator(NOSPACE) {}
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
};
 
/**
 * @brief Modify integration weights on edge to take into account higher-order geometry
 *
 * This structure adjusts integration weights for edge elements when higher-order
 * geometry is present. It accounts for geometric distortions introduced by
 * curved edges and high-order shape functions along the edge.
 *
 * @ingroup mofem_forces_and_sources_tri_element
 */
struct OpSetHOWeightsOnEdge
    : public EdgeElementForcesAndSourcesCore::UserDataOperator {
 
  /**
   * @brief Default constructor for HO weights on edge
   */
  OpSetHOWeightsOnEdge()
      : EdgeElementForcesAndSourcesCore::UserDataOperator(NOSPACE) {}
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
};
 
template <int SPACE_DIM>
struct OpSetHOWeightsOnSubDim;
 
template <> struct OpSetHOWeightsOnSubDim<2> : public OpSetHOWeightsOnEdge {
  using OpSetHOWeightsOnEdge::OpSetHOWeightsOnEdge;
};
 
template <> struct OpSetHOWeightsOnSubDim<3> : public OpSetHOWeightsOnFace {
  using OpSetHOWeightsOnFace::OpSetHOWeightsOnFace;
};
 
/**
 * @brief Set integration weights for higher-order elements
 *
 * This structure sets integration weights by computing the determinant of the
 * Jacobian matrix for higher-order elements. It ensures proper integration
 * when dealing with curved geometries and high-order shape functions.
 */
struct OpSetHOWeights : public ForcesAndSourcesCore::UserDataOperator {
 
  /**
   * @brief Constructor for HO weights calculation
   *
   * @param det_ptr Shared pointer to vector for storing determinant values
   */
  OpSetHOWeights(boost::shared_ptr<VectorDouble> det_ptr)
      : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPSPACE),
        detPtr(det_ptr) {
    if (!detPtr)
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Pointer for detPtr not allocated");
  }
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
 
private:
  boost::shared_ptr<VectorDouble> detPtr;
};
 
/**
 * @brief Apply contravariant (Piola) transformation to Hdiv space for HO geometry
 *
 * This structure applies the contravariant Piola transformation to basis functions
 * in Hdiv spaces when higher-order geometry is present. This transformation is
 * essential for preserving the divergence properties of vector fields across
 * element boundaries in curved geometries.
 *
 * @ingroup mofem_forces_and_sources
 */
struct OpSetHOContravariantPiolaTransform
    : public ForcesAndSourcesCore::UserDataOperator {
 
  /**
   * @brief Constructor for contravariant Piola transformation
   *
   * @param space Field space specification (typically HDIV)
   * @param det_ptr Shared pointer to determinant values for transformation
   * @param jac_ptr Shared pointer to Jacobian matrix for transformation
   */
  OpSetHOContravariantPiolaTransform(const FieldSpace space,
                                     boost::shared_ptr<VectorDouble> det_ptr,
                                     boost::shared_ptr<MatrixDouble> jac_ptr)
      : ForcesAndSourcesCore::UserDataOperator(space, OPSPACE), detPtr(det_ptr),
        jacPtr(jac_ptr) {
    doVertices = false;
    if (space == HDIV)
      doEdges = false;
  }
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
 
private:
  boost::shared_ptr<VectorDouble> detPtr;
  boost::shared_ptr<MatrixDouble> jacPtr;
 
  MatrixDouble piolaN;
  MatrixDouble piolaDiffN;
};
 
/**
 * @brief Apply covariant (Piola) transformation to Hcurl space for HO geometry
 *
 * This structure applies the covariant Piola transformation to basis functions
 * in Hcurl spaces when higher-order geometry is present. This transformation
 * preserves the curl properties of vector fields and maintains tangential
 * continuity across element boundaries in curved geometries.
 */
struct OpSetHOCovariantPiolaTransform
    : public ForcesAndSourcesCore::UserDataOperator {
 
  /**
   * @brief Constructor for covariant Piola transformation
   *
   * @param space Field space specification (typically HCURL)
   * @param jac_inv_ptr Shared pointer to inverse Jacobian matrix for transformation
   */
  OpSetHOCovariantPiolaTransform(const FieldSpace space,
                                 boost::shared_ptr<MatrixDouble> jac_inv_ptr)
      : ForcesAndSourcesCore::UserDataOperator(space, OPSPACE),
        jacInvPtr(jac_inv_ptr) {
    doVertices = false;
    if (space == HDIV)
      doEdges = false;
    if (!jacInvPtr)
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Pointer for jacPtr not allocated");
  }
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
 
private:
  boost::shared_ptr<MatrixDouble> jacInvPtr;
 
  MatrixDouble piolaN;
  MatrixDouble piolaDiffN;
};
 
/** \brief Calculate jacobian for face element
 
  It is assumed that face element is XY plane. Applied
  only for 2d problems and 2d problems embedded in 3d space
 
  \note If you push operators for HO normal befor this operator, HO geometry is
  taken into account when you calculate jacobian.
 
*/
template <int DIM> struct OpCalculateHOJacForFaceImpl;
 
template <>
struct OpCalculateHOJacForFaceImpl<2>
    : public FaceElementForcesAndSourcesCore::UserDataOperator {
 
  OpCalculateHOJacForFaceImpl(boost::shared_ptr<MatrixDouble> jac_ptr);
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
 
protected:
  boost::shared_ptr<MatrixDouble> jacPtr;
};
 
template <>
struct OpCalculateHOJacForFaceImpl<3> : public OpCalculateHOJacForFaceImpl<2> {
 
  using OpCalculateHOJacForFaceImpl<2>::OpCalculateHOJacForFaceImpl;
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
};
 
using OpCalculateHOJacForFace = OpCalculateHOJacForFaceImpl<2>;
using OpCalculateHOJacForFaceEmbeddedIn3DSpace = OpCalculateHOJacForFaceImpl<3>;
 
template <int DIM> struct OpCalculateHOJac;
 
template <> struct OpCalculateHOJac<3> : public OpCalculateHOJacForVolume {
  using OpCalculateHOJacForVolume::OpCalculateHOJacForVolume;
};
 
template <> struct OpCalculateHOJac<2> : public OpCalculateHOJacForFaceImpl<2> {
  using OpCalculateHOJacForFaceImpl<2>::OpCalculateHOJacForFaceImpl;
};
 
/**
 * @brief Calculate normal vectors at Gauss points of face elements
 *
 * This template structure computes normal vectors at integration points for
 * face elements using high-order geometry approximation. It's essential for
 * boundary conditions, surface integrals, and flux calculations in finite
 * element methods.
 *
 * @tparam FIELD_DIM Dimension of the field for normal calculation (default: 3)
 *
 * @ingroup mofem_forces_and_sources
 */
template <int FIELD_DIM = 3>
struct OpGetHONormalsOnFace
    : public FaceElementForcesAndSourcesCore::UserDataOperator {
 
  using OP = FaceElementForcesAndSourcesCore::UserDataOperator;
 
  /**
   * @brief Constructor for HO normal vectors calculation on faces
   *
   * @param field_name Field name where normal vectors will be stored
   * @param add_field If true, add to existing field values; if false, overwrite (default: false)
   */
  OpGetHONormalsOnFace(std::string field_name, bool add_field = false)
      : OP(field_name, OPROW), addField(add_field) {}
 
  OpGetHONormalsOnFace(std::string field_name,
                       boost::shared_ptr<MatrixDouble> tangent1_at_pts,
                       boost::shared_ptr<MatrixDouble> tangent2_at_pts,
                       boost::shared_ptr<MatrixDouble> normals_at_pts,
                       bool add_field = false)
      : OP(field_name, OPROW), tangent1AtPts(tangent1_at_pts),
        tangent2AtPts(tangent2_at_pts), normalsAtPts(normals_at_pts),
        addField(add_field) {}
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
 
private:
  MatrixDouble &getTangent1AtGaussPts() {
    return tangent1AtPts ? *tangent1AtPts : OP::getTangent1AtGaussPts();
  }
 
  MatrixDouble &getTangent2AtGaussPts() {
    return tangent2AtPts ? *tangent2AtPts : OP::getTangent2AtGaussPts();
  }
 
  MatrixDouble &getNormalsAtGaussPts() {
    return normalsAtPts ? *normalsAtPts : OP::getNormalsAtGaussPts();
  }
 
  boost::shared_ptr<MatrixDouble> tangent1AtPts;
  boost::shared_ptr<MatrixDouble> tangent2AtPts;
  boost::shared_ptr<MatrixDouble> normalsAtPts;
 
  bool addField = false; ///< if true add field values, do not overwrite
};
 
/**
 * @brief Transform Hdiv base fluxes from reference element to physical face in 3D
 *
 * This structure applies contravariant Piola transformation to Hdiv basis
 * functions on 3D face elements. It transforms flux values from reference
 * coordinates to physical coordinates while preserving divergence properties.
 *
 * @note Matrix which keeps normal vectors is assumed to have three columns,
 *       and number of rows should be equal to number of integration points.
 *
 * @ingroup mofem_forces_and_sources
 */
struct OpHOSetContravariantPiolaTransformOnFace3D
    : public FaceElementForcesAndSourcesCore::UserDataOperator {
 
  /**
   * @brief Constructor for contravariant Piola transformation on 3D faces
   *
   * @param space Field space specification (typically HDIV)
   * @param normals_at_gauss_pts Shared pointer to matrix containing normal vectors at integration points (optional)
   */
  OpHOSetContravariantPiolaTransformOnFace3D(
      const FieldSpace space,
      boost::shared_ptr<MatrixDouble> normals_at_gauss_pts = nullptr)
      : FaceElementForcesAndSourcesCore::UserDataOperator(space, OPSPACE),
        normalsAtGaussPts(normals_at_gauss_pts) {}
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
 
private:
  boost::shared_ptr<MatrixDouble> normalsAtGaussPts;
};
 
/**
 * @brief Transform Hcurl base functions from reference element to physical edge in 3D
 *
 * This structure applies contravariant Piola transformation to Hcurl basis
 * functions on 3D edge elements. It transforms curl-conforming basis functions
 * from reference coordinates to physical coordinates while preserving tangential
 * continuity across element boundaries.
 *
 * @ingroup mofem_forces_and_sources
 */
struct OpHOSetContravariantPiolaTransformOnEdge3D
    : public EdgeElementForcesAndSourcesCore::UserDataOperator {
 
  /**
   * @brief Constructor for contravariant Piola transformation on 3D edges
   *
   * @param space Field space specification (default: HCURL)
   * @param tangent_at_pts Shared pointer to matrix containing tangent vectors at integration points (optional)
   */
  OpHOSetContravariantPiolaTransformOnEdge3D(
      const FieldSpace space = HCURL,
      boost::shared_ptr<MatrixDouble> tangent_at_pts = nullptr)
      : EdgeElementForcesAndSourcesCore::UserDataOperator(space),
        tangentAtGaussPts(tangent_at_pts) {}
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
 
private:
  boost::shared_ptr<MatrixDouble> tangentAtGaussPts;
};
 
/**
 * @brief Transform Hcurl base functions from reference element to physical face in 3D
 *
 * This structure applies covariant Piola transformation to Hcurl basis functions
 * on 3D face elements. It transforms curl-conforming basis functions from
 * reference coordinates to physical coordinates using normal and tangent vectors,
 * preserving the curl properties of the field.
 *
 * @ingroup mofem_forces_and_sources
 */
struct OpHOSetCovariantPiolaTransformOnFace3D
    : public FaceElementForcesAndSourcesCore::UserDataOperator {
 
  /**
   * @brief Constructor for covariant Piola transformation on 3D faces
   *
   * @param space Field space specification (typically HCURL)
   * @param normals_at_pts Shared pointer to matrix containing normal vectors at integration points (optional)
   * @param tangent1_at_pts Shared pointer to matrix containing first tangent vectors at integration points (optional)
   * @param tangent2_at_pts Shared pointer to matrix containing second tangent vectors at integration points (optional)
   */
  OpHOSetCovariantPiolaTransformOnFace3D(
      const FieldSpace space,
      boost::shared_ptr<MatrixDouble> normals_at_pts = nullptr,
      boost::shared_ptr<MatrixDouble> tangent1_at_pts = nullptr,
      boost::shared_ptr<MatrixDouble> tangent2_at_pts = nullptr)
      : FaceElementForcesAndSourcesCore::UserDataOperator(space, OPSPACE),
        normalsAtPts(normals_at_pts), tangent1AtPts(tangent1_at_pts),
        tangent2AtPts(tangent2_at_pts) {
    if (normals_at_pts || tangent1_at_pts || tangent2_at_pts)
      if (normals_at_pts && tangent1_at_pts && tangent2_at_pts)
        CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                          "All elements in constructor have to set pointer");
  }
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
 
private:
  boost::shared_ptr<MatrixDouble> normalsAtPts;
  boost::shared_ptr<MatrixDouble> tangent1AtPts;
  boost::shared_ptr<MatrixDouble> tangent2AtPts;
 
  MatrixDouble piolaN;
  MatrixDouble diffPiolaN;
};
 
/**
 * @brief Calculate tangent vector on edge from HO geometry approximation
 *
 * This template structure computes tangent vectors along edges using high-order
 * geometry approximation. It extracts tangent information at integration points
 * and stores the results for further processing in finite element calculations.
 *
 * @tparam FIELD_DIM Dimension of the field (default: 3)
 *
 * @ingroup mofem_forces_and_sources
 */
template <int FIELD_DIM = 3>
struct OpGetHOTangentsOnEdge
    : public EdgeElementForcesAndSourcesCore::UserDataOperator {
 
  /**
   * @brief Constructor for HO tangents calculation on edges
   *
   * @param field_name Name of the field used for tangent calculation
   * @param tangents_at_pts Shared pointer to matrix for storing calculated tangents (optional)
   */
  OpGetHOTangentsOnEdge(
      std::string field_name,
      boost::shared_ptr<MatrixDouble> tangents_at_pts = nullptr)
      : EdgeElementForcesAndSourcesCore::UserDataOperator(field_name, OPROW),
        tangentsAtPts(tangents_at_pts) {}
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
 
private:
  boost::shared_ptr<MatrixDouble> tangentsAtPts;
};
 
/**
 * @brief Scale base functions by inverses of measure of element
 *
 * @tparam OP
 */
struct OpScaleBaseBySpaceInverseOfMeasure
    : public ForcesAndSourcesCore::UserDataOperator {
 
  using OP = ForcesAndSourcesCore::UserDataOperator;
 
  OpScaleBaseBySpaceInverseOfMeasure(
      const FieldSpace space, const FieldApproximationBase base,
      boost::shared_ptr<VectorDouble> det_jac_ptr,
      boost::shared_ptr<double> scale_det_ptr = nullptr);
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
 
private:
  FieldSpace fieldSpace;
  FieldApproximationBase fieldBase;
  boost::shared_ptr<VectorDouble> detJacPtr;
  boost::shared_ptr<double> scaleDetPtr = nullptr;
};
 
/**
 * @brief Add operators pushing bases from local to physical configuration
 *
 * @tparam FE_DIM dimension of element
 * @tparam PROBLEM_DIM problem dimension
 * @tparam SPACE_DIM space dimension
 */
template <int FE_DIM, int PROBLEM_DIM, int SPACE_DIM> struct AddHOOps;
 
/**
 * @brief Specialization for 2D face elements in 2D problems with 2D space
 *
 * This specialization adds high-order operators for 2D face elements operating
 * in 2D problems within 2D coordinate space. Commonly used for planar face
 * elements in 2D finite element analysis.
 */
template <> struct AddHOOps<2, 2, 2> {
  AddHOOps() = delete;
  
  /**
   * @brief Add high-order operators to the pipeline for 2D face elements
   *
   * @param pipeline Pipeline to add operators to
   * @param spaces Vector of field spaces to be handled
   * @param geom_field_name Name of the geometry field (default: empty string)
   * @param jac Shared pointer to Jacobian matrix. If provided, values will be set by this function.
   *            If null, Jacobian will be handled internally.
   * @param det Shared pointer to determinant vector. If provided, values will be set by this function.
   *            If null, determinant will be handled internally.
   * @param inv_jac Shared pointer to inverse Jacobian matrix. If provided, values will be set by this function.
   *                If null, inverse Jacobian will be handled internally.
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  static MoFEMErrorCode
  add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
      std::vector<FieldSpace> spaces, std::string geom_field_name = "",
      boost::shared_ptr<MatrixDouble> jac = nullptr,
      boost::shared_ptr<VectorDouble> det = nullptr,
      boost::shared_ptr<MatrixDouble> inv_jac = nullptr);
};
 
/**
 * @brief Specialization for 2D face elements in 2D problems embedded in 3D space
 *
 * This specialization adds high-order operators for 2D face elements operating
 * in 2D problems but embedded within 3D coordinate space. Typically used for
 * surface elements or membranes in 3D environments.
 */
template <> struct AddHOOps<2, 2, 3> {
  AddHOOps() = delete;
  
  /**
   * @brief Add high-order operators to the pipeline for 2D face elements embedded in 3D
   *
   * @param pipeline Pipeline to add operators to
   * @param spaces Vector of field spaces to be handled
   * @param geom_field_name Name of the geometry field (default: empty string)
   * @param jac Shared pointer to Jacobian matrix. If provided, values will be set by this function.
   *            If null, Jacobian will be handled internally.
   * @param det Shared pointer to determinant vector. If provided, values will be set by this function.
   *            If null, determinant will be handled internally.
   * @param inv_jac Shared pointer to inverse Jacobian matrix. If provided, values will be set by this function.
   *                If null, inverse Jacobian will be handled internally.
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  static MoFEMErrorCode
  add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
      std::vector<FieldSpace> spaces, std::string geom_field_name = "",
      boost::shared_ptr<MatrixDouble> jac = nullptr,
      boost::shared_ptr<VectorDouble> det = nullptr,
      boost::shared_ptr<MatrixDouble> inv_jac = nullptr);
};
 
/**
 * @brief Specialization for 1D edge elements in 2D problems with 2D space
 *
 * This specialization adds high-order operators for 1D edge elements operating
 * in 2D problems within 2D coordinate space. Used for edge-based computations
 * such as boundary conditions or line integrals in 2D finite element analysis.
 */
template <> struct AddHOOps<1, 2, 2> {
  AddHOOps() = delete;
  
  /**
   * @brief Add high-order operators to the pipeline for 1D edge elements in 2D
   *
   * @param pipeline Pipeline to add operators to
   * @param spaces Vector of field spaces to be handled
   * @param geom_field_name Name of the geometry field (default: empty string)
   * @return MoFEMErrorCode Error code indicating success or failure
   * 
   * @note This specialization handles Jacobian, determinant, and inverse Jacobian internally.
   */
  static MoFEMErrorCode
  add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
      std::vector<FieldSpace> spaces, std::string geom_field_name = "");
};
 
/**
 * @brief Specialization for 1D edge elements in 3D problems with 3D space
 *
 * This specialization adds high-order operators for 1D edge elements operating
 * in 3D problems within 3D coordinate space. Used for edge-based computations
 * such as line integrals, boundary conditions, or wire-like structures in 3D
 * finite element analysis.
 */
template <> struct AddHOOps<1, 3, 3> {
  AddHOOps() = delete;
  
  /**
   * @brief Add high-order operators to the pipeline for 1D edge elements in 3D
   *
   * @param pipeline Pipeline to add operators to
   * @param space Vector of field spaces to be handled
   * @param geom_field_name Name of the geometry field (default: empty string)
   * @return MoFEMErrorCode Error code indicating success or failure
   * 
   * @note This specialization handles Jacobian, determinant, and inverse Jacobian internally.
   */
  static MoFEMErrorCode
  add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
      std::vector<FieldSpace> space, std::string geom_field_name = "");
};
 
/**
 * @brief Specialization for 2D face elements in 3D problems with 3D space
 *
 * This specialization adds high-order operators for 2D face elements operating
 * in 3D problems within 3D coordinate space. Commonly used for surface elements,
 * boundary conditions, membranes, or shell structures in 3D finite element
 * analysis.
 */
template <> struct AddHOOps<2, 3, 3> {
  AddHOOps() = delete;
  
  /**
   * @brief Add high-order operators to the pipeline for 2D face elements in 3D
   *
   * @param pipeline Pipeline to add operators to
   * @param space Vector of field spaces to be handled
   * @param geom_field_name Name of the geometry field (default: empty string)
   * @return MoFEMErrorCode Error code indicating success or failure
   * 
   * @note This specialization handles Jacobian, determinant, and inverse Jacobian internally.
   */
  static MoFEMErrorCode
  add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
      std::vector<FieldSpace> space, std::string geom_field_name = "");
};
 
/**
 * @brief Specialization for 3D volume elements in 3D problems with 3D space
 *
 * This specialization adds high-order operators for 3D volume elements operating
 * in 3D problems within 3D coordinate space. This is the most common case for
 * solid mechanics, heat transfer, and other 3D finite element applications
 * involving volumetric elements like tetrahedra, hexahedra, and prisms.
 */
template <> struct AddHOOps<3, 3, 3> {
  AddHOOps() = delete;
  
  /**
   * @brief Add high-order operators to the pipeline for 3D volume elements
   *
   * @param pipeline Pipeline to add operators to
   * @param space Vector of field spaces to be handled
   * @param geom_field_name Name of the geometry field (default: empty string)
   * @param jac Shared pointer to Jacobian matrix. If provided, values will be set by this function.
   *            If null, Jacobian will be handled internally.
   * @param det Shared pointer to determinant vector. If provided, values will be set by this function.
   *            If null, determinant will be handled internally.
   * @param inv_jac Shared pointer to inverse Jacobian matrix. If provided, values will be set by this function.
   *                If null, inverse Jacobian will be handled internally.
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  static MoFEMErrorCode
  add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
      std::vector<FieldSpace> space, std::string geom_field_name = "",
      boost::shared_ptr<MatrixDouble> jac = nullptr,
      boost::shared_ptr<VectorDouble> det = nullptr,
      boost::shared_ptr<MatrixDouble> inv_jac = nullptr);
};
 
template <int FIELD_DIM>
MoFEMErrorCode
OpCalculateHOCoords<FIELD_DIM>::doWork(int side, EntityType type,
                                       EntitiesFieldData::EntData &data) {
  FTensor::Index<'i', FIELD_DIM> i;
  MoFEMFunctionBegin;
  const auto nb_dofs = data.getFieldData().size() / FIELD_DIM;
  if (nb_dofs) {
    if (type == MBVERTEX)
      getCoordsAtGaussPts().clear();
    auto t_base = data.getFTensor0N();
    auto t_coords = getFTensor1CoordsAtGaussPts();
    const auto nb_integration_pts = data.getN().size1();
    const auto nb_base_functions = data.getN().size2();
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
      auto t_dof = data.getFTensor1FieldData<FIELD_DIM>();
      size_t bb = 0;
      for (; bb != nb_dofs; ++bb) {
        t_coords(i) += t_base * t_dof(i);
        ++t_dof;
        ++t_base;
      }
      for (; bb < nb_base_functions; ++bb)
        ++t_base;
      ++t_coords;
    }
  }
  MoFEMFunctionReturn(0);
}
 
template <int FIELD_DIM>
MoFEMErrorCode
OpGetHONormalsOnFace<FIELD_DIM>::doWork(int side, EntityType type,
                                        EntitiesFieldData::EntData &data) {
  MoFEMFunctionBeginHot;
 
  FTensor::Number<0> N0;
  FTensor::Number<1> N1;
 
  auto get_ftensor1 = [](MatrixDouble &m) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
        &m(0, 0), &m(0, 1), &m(0, 2));
  };
 
  unsigned int nb_dofs = data.getFieldData().size();
  if (nb_dofs != 0) {
 
    int nb_gauss_pts = data.getN().size1();
    auto &tangent1_at_gauss_pts = getTangent1AtGaussPts();
    auto &tangent2_at_gauss_pts = getTangent2AtGaussPts();
    tangent1_at_gauss_pts.resize(nb_gauss_pts, 3, false);
    tangent2_at_gauss_pts.resize(nb_gauss_pts, 3, false);
 
    switch (type) {
    case MBVERTEX: {
      if (!addField) {
        tangent1_at_gauss_pts.clear();
        tangent2_at_gauss_pts.clear();
      }
    }
    case MBEDGE:
    case MBTRI:
    case MBQUAD: {
 
#ifndef NDEBUG
      if (2 * data.getN().size2() != data.getDiffN().size2()) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "expected two derivative in local element coordinates");
      }
      if (nb_dofs % FIELD_DIM != 0) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "expected that number of dofs is multiplicative of field "
                "dimension");
      }
#endif
 
      if (nb_dofs > FIELD_DIM * data.getN().size2()) {
        unsigned int nn = 0;
        for (; nn != nb_dofs; nn++) {
          if (!data.getFieldDofs()[nn]->getActive())
            break;
        }
        if (nn > FIELD_DIM * data.getN().size2()) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                   "Data inconsistency for base %s",
                   ApproximationBaseNames[data.getBase()]);
        } else {
          nb_dofs = nn;
          if (!nb_dofs)
            MoFEMFunctionReturnHot(0);
        }
      }
      const int nb_base_functions = data.getN().size2();
      auto t_base = data.getFTensor0N();
      auto t_diff_base = data.getFTensor1DiffN<2>();
      auto t_t1 = get_ftensor1(tangent1_at_gauss_pts);
      auto t_t2 = get_ftensor1(tangent2_at_gauss_pts);
      for (int gg = 0; gg != nb_gauss_pts; ++gg) {
        auto t_data = data.getFTensor1FieldData<FIELD_DIM>();
        int bb = 0;
        for (; bb != nb_dofs / FIELD_DIM; ++bb) {
          FTensor::Index<'i', FIELD_DIM> i;
          t_t1(i) += t_data(i) * t_diff_base(N0);
          t_t2(i) += t_data(i) * t_diff_base(N1);
          ++t_data;
          ++t_base;
          ++t_diff_base;
        }
        for (; bb != nb_base_functions; ++bb) {
          ++t_base;
          ++t_diff_base;
        }
        ++t_t1;
        ++t_t2;
      }
    } break;
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
    }
  }
 
  if (moab::CN::Dimension(type) == 2) {
 
    auto &normal_at_gauss_pts = getNormalsAtGaussPts();
    auto &tangent1_at_gauss_pts = getTangent1AtGaussPts();
    auto &tangent2_at_gauss_pts = getTangent2AtGaussPts();
 
    const auto nb_gauss_pts = tangent1_at_gauss_pts.size1();
    normal_at_gauss_pts.resize(nb_gauss_pts, 3, false);
 
    auto t_normal = get_ftensor1(normal_at_gauss_pts);
    auto t_t1 = get_ftensor1(tangent1_at_gauss_pts);
    auto t_t2 = get_ftensor1(tangent2_at_gauss_pts);
    for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
      FTensor::Index<'i', 3> i;
      FTensor::Index<'j', 3> j;
      FTensor::Index<'k', 3> k;
      t_normal(j) = FTensor::levi_civita(i, j, k) * t_t1(k) * t_t2(i);
      ++t_normal;
      ++t_t1;
      ++t_t2;
    }
  }
 
  MoFEMFunctionReturnHot(0);
}
 
template <int FIELD_DIM>
MoFEMErrorCode
OpGetHOTangentsOnEdge<FIELD_DIM>::doWork(int side, EntityType type,
                                         EntitiesFieldData::EntData &data) {
  MoFEMFunctionBegin;
 
  auto get_tangent = [&]() -> MatrixDouble & {
    if (tangentsAtPts)
      return *tangentsAtPts;
    else
      return getTangentAtGaussPts();
  };
 
  auto &tangent = get_tangent();
  int nb_gauss_pts = getGaussPts().size2();
 
  if (type == MBVERTEX) {
    tangent.resize(nb_gauss_pts, 3, false);
    tangent.clear();
  }
 
  int nb_dofs = data.getFieldData().size();
  if (nb_dofs != 0) {
    const int nb_base_functions = data.getN().size2();
    double *diff_base_function = &*data.getDiffN().data().begin();
    auto tangent_at_gauss_pts =
        getFTensor1FromPtr<FIELD_DIM, 3>(&*tangent.data().begin());
 
    FTensor::Index<'i', FIELD_DIM> i;
    int size = nb_dofs / FIELD_DIM;
    if (nb_dofs % FIELD_DIM) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Inconsistent number of dofs and field dimension");
    }
    for (int gg = 0; gg != nb_gauss_pts; ++gg) {
      auto field_data = data.getFTensor1FieldData<FIELD_DIM>();
      int bb = 0;
      for (; bb < size; ++bb) {
        tangent_at_gauss_pts(i) += field_data(i) * (*diff_base_function);
        ++field_data;
        ++diff_base_function;
      }
      for (; bb != nb_base_functions; ++bb) {
        ++diff_base_function;
      }
      ++tangent_at_gauss_pts;
    }
  }
  MoFEMFunctionReturn(0);
}
 
/**
 * @deprecated do not use this function, instead use AddHOOps.
 *
 * @tparam E
 * @param field
 * @param e
 * @param h1
 * @param hcurl
 * @param hdiv
 * @param l2
 * @return MoFEMErrorCode
 */
template <typename E>
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1,
                                      bool hcurl, bool hdiv, bool l2) {
  std::vector<FieldSpace> spaces;
  if (h1)
    spaces.push_back(H1);
  if (hcurl)
    spaces.push_back(HCURL);
  if (hdiv)
    spaces.push_back(HDIV);
  if (l2)
    spaces.push_back(L2);
  return AddHOOps<3, 3, 3>::add(e.getOpPtrVector(), spaces, field);
}
 
/**
 * @deprecated do not use this function, instead use AddHOOps.
 *
 * @tparam E
 * @param field
 * @param e
 * @param hcurl
 * @param hdiv
 * @return MoFEMErrorCode
 */
template <typename E>
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e,
                                         bool hcurl, bool hdiv) {
  std::vector<FieldSpace> spaces;
  if (hcurl)
    spaces.push_back(HCURL);
  if (hdiv)
    spaces.push_back(HDIV);
  return AddHOOps<2, 3, 3>::add(e.getOpPtrVector(), spaces, field);
}
 
}; // namespace MoFEM
 
#endif