diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index d9c175fec20..7682c470520 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -139,6 +139,11 @@ blt_add_executable( NAME geosx
                                ${extraComponentsLinkList}
                                ${externalComponentsLinkList} )
 
+if(APPLE)
+  target_link_options(geosx PRIVATE "-Wl,-export_dynamic")
+endif()
+
+
 # Seems to be required on some CMake versions (e.g. 3.16) to get enforce device linking
 if( ${ENABLE_HYPRE_DEVICE} STREQUAL "CUDA"  )
   set_target_properties( geosx PROPERTIES CUDA_RESOLVE_DEVICE_SYMBOLS TRUE )
diff --git a/src/coreComponents/finiteElement/CMakeLists.txt b/src/coreComponents/finiteElement/CMakeLists.txt
index 53ef7404158..7a12528dc0a 100644
--- a/src/coreComponents/finiteElement/CMakeLists.txt
+++ b/src/coreComponents/finiteElement/CMakeLists.txt
@@ -29,6 +29,7 @@ set( finiteElement_headers
      LinearFormUtilities.hpp
      PDEUtilities.hpp
      elementFormulations/FiniteElementBase.hpp
+     elementFormulations/FiniteElementOperators.hpp
      elementFormulations/H1_Hexahedron_Lagrange1_GaussLegendre2.hpp
      elementFormulations/H1_Pyramid_Lagrange1_Gauss5.hpp
      elementFormulations/H1_QuadrilateralFace_Lagrange1_GaussLegendre2.hpp
diff --git a/src/coreComponents/finiteElement/FiniteElementDiscretization.hpp b/src/coreComponents/finiteElement/FiniteElementDiscretization.hpp
index 21d7e9f71d5..d6f53c26b2e 100644
--- a/src/coreComponents/finiteElement/FiniteElementDiscretization.hpp
+++ b/src/coreComponents/finiteElement/FiniteElementDiscretization.hpp
@@ -153,7 +153,7 @@ FiniteElementDiscretization::
     typename FE_TYPE::StackVariables feStack;
     finiteElement.template setup< FE_TYPE >( k, meshData, feStack );
     real64 xLocal[numNodesPerElem][3];
-    localIndex numSupportPoints = finiteElement.template numSupportPoints< FE_TYPE >( feStack );
+    localIndex numSupportPoints = finiteElement.getNumSupportPoints( feStack );
     for( localIndex a=0; a< numSupportPoints; ++a )
     {
       localIndex const nodeIndex = elemsToNodes[ k][ a ];
diff --git a/src/coreComponents/finiteElement/elementFormulations/FiniteElementBase.hpp b/src/coreComponents/finiteElement/elementFormulations/FiniteElementBase.hpp
index b834a019513..4b206420e40 100644
--- a/src/coreComponents/finiteElement/elementFormulations/FiniteElementBase.hpp
+++ b/src/coreComponents/finiteElement/elementFormulations/FiniteElementBase.hpp
@@ -17,11 +17,6 @@
  * @file FiniteElementBase.hpp
  */
 
-#if defined(GEOS_USE_DEVICE)
-#define CALC_FEM_SHAPE_IN_KERNEL
-#endif
-
-
 
 #ifndef GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_FINITEELEMENTBASE_HPP_
 #define GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_FINITEELEMENTBASE_HPP_
@@ -39,6 +34,7 @@ namespace geos
 namespace finiteElement
 {
 
+
 /**
  * @brief Base class for FEM element implementations.
  */
@@ -57,18 +53,7 @@ class FiniteElementBase
    * @brief Copy Constructor
    * @param source The object to copy.
    */
-  GEOS_HOST_DEVICE
-  FiniteElementBase( FiniteElementBase const & source ):
-#ifdef CALC_FEM_SHAPE_IN_KERNEL
-    m_viewGradN(),
-    m_viewDetJ()
-#else
-    m_viewGradN( source.m_viewGradN ),
-    m_viewDetJ( source.m_viewDetJ )
-#endif
-  {
-    GEOS_UNUSED_VAR( source ); // suppress warning when CALC_FEM_SHAPE_IN_KERNEL is defined
-  }
+  FiniteElementBase( FiniteElementBase const & source ) = default;
 
   /// Default Move constructor
   FiniteElementBase( FiniteElementBase && ) = default;
@@ -228,19 +213,6 @@ class FiniteElementBase
   GEOS_HOST_DEVICE
   constexpr static PDEUtilities::FunctionSpace getFunctionSpace();
 
-  /**
-   * @brief Getter for the number of support points per element.
-   * @tparam LEAF Type of the derived finite element implementation.
-   * @param stack Stack variables created by a call to @ref setup.
-   * @return The number of support points per element.
-   */
-  template< typename LEAF >
-  GEOS_HOST_DEVICE
-  localIndex numSupportPoints( typename LEAF::StackVariables const & stack ) const
-  {
-    return LEAF::getNumSupportPoints( stack );
-  }
-
   /**
    * @brief Get the maximum number of support points for this element.
    * @details This should be used to know the size of pre-allocated objects whose size depend on the
@@ -250,80 +222,6 @@ class FiniteElementBase
   GEOS_HOST_DEVICE
   virtual localIndex getMaxSupportPoints() const = 0;
 
-  /**
-   * @brief Get the shape function gradients.
-   * @tparam LEAF Type of the derived finite element implementation.
-   * @param k The element index.
-   * @param q The quadrature point index.
-   * @param X Array of coordinates as the reference for the gradients.
-   * @param gradN Return array of the shape function gradients.
-   * @return The determinant of the Jacobian transformation matrix.
-   *
-   * This function calls the function to calculate shape function gradients.
-   */
-  template< typename LEAF >
-  GEOS_HOST_DEVICE
-  real64 getGradN( localIndex const k,
-                   localIndex const q,
-                   real64 const (&X)[LEAF::maxSupportPoints][3],
-                   real64 ( &gradN )[LEAF::maxSupportPoints][3] ) const;
-
-  /**
-   * @brief Get the shape function gradients.
-   * @tparam LEAF Type of the derived finite element implementation.
-   * @param k The element index.
-   * @param q The quadrature point index.
-   * @param X Array of coordinates as the reference for the gradients.
-   * @param stack Stack variables relative to the element @p k created by a call to @ref setup.
-   * @param gradN Return array of the shape function gradients.
-   * @return The determinant of the Jacobian transformation matrix.
-   *
-   * This function calls the function to calculate shape function gradients.
-   */
-  template< typename LEAF >
-  GEOS_HOST_DEVICE
-  real64 getGradN( localIndex const k,
-                   localIndex const q,
-                   real64 const (&X)[LEAF::maxSupportPoints][3],
-                   typename LEAF::StackVariables const & stack,
-                   real64 ( &gradN )[LEAF::maxSupportPoints][3] ) const;
-
-  /**
-   * @brief Get the shape function gradients.
-   * @tparam LEAF Type of the derived finite element implementation.
-   * @param k The element index.
-   * @param q The quadrature point index.
-   * @param X dummy variable.
-   * @param gradN Return array of the shape function gradients.
-   * @return The determinant of the Jacobian transformation matrix.
-   *
-   * This function returns pre-calculated shape function gradients.
-   */
-  template< typename LEAF >
-  GEOS_HOST_DEVICE
-  real64 getGradN( localIndex const k,
-                   localIndex const q,
-                   int const X,
-                   real64 ( &gradN )[LEAF::maxSupportPoints][3] ) const;
-  /**
-   * @brief Get the shape function gradients.
-   * @tparam LEAF Type of the derived finite element implementation.
-   * @param k The element index.
-   * @param q The quadrature point index.
-   * @param X dummy variable.
-   * @param stack Stack variables relative to the element @p k created by a call to @ref setup.
-   * @param gradN Return array of the shape function gradients.
-   * @return The determinant of the Jacobian transformation matrix.
-   *
-   * This function returns pre-calculated shape function gradients.
-   */
-  template< typename LEAF >
-  GEOS_HOST_DEVICE
-  real64 getGradN( localIndex const k,
-                   localIndex const q,
-                   int const X,
-                   typename LEAF::StackVariables const & stack,
-                   real64 ( &gradN )[LEAF::maxSupportPoints][3] ) const;
 
 
   /**
@@ -433,329 +331,10 @@ class FiniteElementBase
                                                                                               scaleFactor );
   }
 
-  /**
-   * @name Value Operator Functions
-   */
-  ///@{
-
-  /**
-   * @brief Compute the interpolated value of a variable.
-   * @tparam NUM_SUPPORT_POINTS The number of support points for the element.
-   * @param N Array (for each support point) of shape function values at the
-   *   coordinate the variable is to be interpolated.
-   * @param var Array of variable values for each support point.
-   * @param value The interpolated value of @p var.
-   *
-   * This is the standard finite element interpolation operator of a discrete
-   * variable defined at the support points.
-   * The operator is expressed as:
-   * \f[
-   * value  = \sum_a^{numSupport} \left ( N_a var_a \right ),
-   * \f]
-
-   * @note The shape function values @p N must be evaluated prior to calling this
-   * function.
-   *
-   */
-  template< int NUM_SUPPORT_POINTS >
-  GEOS_HOST_DEVICE
-  static
-  void value( real64 const (&N)[NUM_SUPPORT_POINTS],
-              real64 const (&var)[NUM_SUPPORT_POINTS],
-              real64 & value );
-
-  /**
-   * @brief Compute the interpolated value of a vector variable.
-   * @tparam NUM_COMPONENTS Number of components for the vector variable.
-   * @copydoc value
-   */
-  template< int NUM_SUPPORT_POINTS,
-            int NUM_COMPONENTS >
-  GEOS_HOST_DEVICE
-  static
-  void value( real64 const (&N)[NUM_SUPPORT_POINTS],
-              real64 const (&var)[NUM_SUPPORT_POINTS][NUM_COMPONENTS],
-              real64 ( &value )[NUM_COMPONENTS] );
-
-  ///@}
-
-  /**
-   * @name Gradient Operator Functions
-   */
-  ///@{
-
-  /**
-   * @brief Calculate the symmetric gradient of a vector valued support field
-   *   at a point using the stored basis function gradients for all support
-   *   points.
-   * @tparam NUM_SUPPORT_POINTS The number of support points for the element.
-   * @tparam GRADIENT_TYPE The type of the array object holding the shape
-   * @param gradN The basis function gradients at a point in the element.
-   * @param var The vector valued support field that the gradient operator will
-   *  be applied to.
-   * @param gradVar The symmetric gradient in Voigt notation.
-   *
-   * More precisely, the operator is defined as:
-   * \f[
-   * grad^s_{ij}  = \frac{1}{2} \sum_a^{nSupport} \left ( \frac{\partial N_a}{\partial X_j} var_{ai} + \frac{\partial N_a}{\partial X_i}
-   * var_{aj}\right ),
-   * \f]
-   *
-   */
-  template< int NUM_SUPPORT_POINTS,
-            typename GRADIENT_TYPE >
-  GEOS_HOST_DEVICE
-  static void symmetricGradient( GRADIENT_TYPE const & gradN,
-                                 real64 const (&var)[NUM_SUPPORT_POINTS][3],
-                                 real64 ( &gradVar )[6] );
-
-  /**
-   * @brief Calculate the trace of the symmetric gradient of a vector valued support
-   *   field (i.e. the volumetric strain for the displacement field) at a point using
-   *   the stored basis function gradients for all support points.
-   * @tparam NUM_SUPPORT_POINTS The number of support points for the element.
-   * @tparam GRADIENT_TYPE The type of the array object holding the shape
-   * @param gradN The basis function gradients at a point in the element.
-   * @param var The vector valued support field that the gradient operator will
-   *  be applied to.
-   * @return The trace of the symetric gradient tensor.
-   *
-   */
-  template< int NUM_SUPPORT_POINTS,
-            typename GRADIENT_TYPE >
-  GEOS_HOST_DEVICE
-  static real64 symmetricGradientTrace( GRADIENT_TYPE const & gradN,
-                                        real64 const (&var)[NUM_SUPPORT_POINTS][3] );
-
-  /**
-   * @brief Calculate the gradient of a scalar valued support field at a point
-   *   using the input basis function gradients.
-   * @tparam NUM_SUPPORT_POINTS The number of support points for the element.
-   * @tparam GRADIENT_TYPE The type of the array object holding the shape
-   *   function gradients.
-   * @param gradN The basis function gradients at a point in the element.
-   * @param var The vector valued support field that the gradient operator will
-   *  be applied to.
-   * @param gradVar The  gradient.
-   *
-   * More precisely, the operator is defined as:
-   * \f[
-   * grad_{j}  = \sum_a^{nSupport} \left ( \frac{\partial N_a}{\partial X_j} var_{a}\right ),
-   * \f]
-   */
-  template< int NUM_SUPPORT_POINTS,
-            typename GRADIENT_TYPE >
-  GEOS_HOST_DEVICE
-  static void gradient( GRADIENT_TYPE const & gradN,
-                        real64 const (&var)[NUM_SUPPORT_POINTS],
-                        real64 ( &gradVar )[3] );
-
-  /**
-   * @brief Calculate the gradient of a vector valued support field at a point
-   *   using the input basis function gradients.
-   * @copydoc gradient
-   *
-   * More precisely, the operator is defined as:
-   * \f[
-   * grad_{ij}  = \sum_a^{nSupport} \left ( \frac{\partial N_a}{\partial X_j} var_{ai}\right ),
-   * \f]
-   */
-  template< int NUM_SUPPORT_POINTS,
-            typename GRADIENT_TYPE >
-  GEOS_HOST_DEVICE
-  static void gradient( GRADIENT_TYPE const & gradN,
-                        real64 const (&var)[NUM_SUPPORT_POINTS][3],
-                        real64 ( &gradVar )[3][3] );
-  ///@}
-
-  /**
-   * @name Multi-Operator Functions
-   */
-  ///@{
-
-  /**
-   * @brief Calculate the value and gradient of a scalar valued support field
-   *   at a point using the input basis function gradients.
-   * @tparam NUM_SUPPORT_POINTS The number of support points for the element.
-   * @tparam GRADIENT_TYPE The type of the array object holding the shape
-   * @param N Array (for each support point) of shape function values at the
-   *   coordinate the variable is to be interpolated.
-   * @param gradN The basis function gradients at a point in the element.
-   * @param var The vector valued support field that the gradient operator will
-   *  be applied to.
-   * @param value The value at the point for which N was specified.
-   * @param gradVar The gradient at the point for which gradN was specified.
-   */
-  template< int NUM_SUPPORT_POINTS,
-            typename GRADIENT_TYPE >
-  GEOS_HOST_DEVICE
-  static void valueAndGradient( real64 const (&N)[NUM_SUPPORT_POINTS],
-                                GRADIENT_TYPE const & gradN,
-                                real64 const (&var)[NUM_SUPPORT_POINTS],
-                                real64 & value,
-                                real64 ( &gradVar )[3] );
-
-  ///@}
-
-  /**
-   * @name Scattering Operator Functions
-   *
-   * These functions take quadrature data and map it to the support points
-   * through some operator.
-   */
-  ///@{
-
-  /**
-   * @brief Inner product of each basis function gradient with a rank-2
-   *   symmetric tensor.
-   * @tparam NUM_SUPPORT_POINTS The number of support points for the element.
-   * @tparam GRADIENT_TYPE The type of the array object holding the shape
-   *   function gradients.
-   * @param gradN The basis function gradients at a point in the element.
-   * @param var_detJxW The rank-2 tensor at @p q scaled by J*W.
-   * @param R The vector at each support point which will hold the result from
-   *   the tensor contraction.
-   *
-   * More precisely, the operator is defined as:
-   *
-   * \f[
-   * R_i = \sum_a^{nSupport} \left( \frac{\partial N_a}{\partial X_j} var_{ij}\right),
-   * \f]
-   *
-   * where \f$\frac{\partial N_a}{\partial X_j}\f$ is the basis function gradient,
-   *   \f$var_{ij}\f$ is the rank-2 symmetric tensor.
-   */
-  template< int NUM_SUPPORT_POINTS,
-            typename GRADIENT_TYPE >
-  GEOS_HOST_DEVICE
-  static void plusGradNajAij( GRADIENT_TYPE const & gradN,
-                              real64 const (&var_detJxW)[6],
-                              real64 ( &R )[NUM_SUPPORT_POINTS][3] );
-
-  /**
-   * @copydoc plusGradNajAij
-   * @brief Inner product of each basis function gradient with a rank-2
-   *   tensor.
-   */
-  template< int NUM_SUPPORT_POINTS,
-            typename GRADIENT_TYPE >
-  GEOS_HOST_DEVICE
-  static void plusGradNajAij( GRADIENT_TYPE const & gradN,
-                              real64 const (&var_detJxW)[3][3],
-                              real64 ( &R )[NUM_SUPPORT_POINTS][3] );
-
-  /**
-   * @brief Product of each shape function with a vector forcing term.
-   * @tparam NUM_SUPPORT_POINTS The number of support points for the element.
-   * @param N The shape function value at a predetermined coordinate in the element.
-   * @param forcingTerm_detJxW A vector scaled by detJxW
-   * @param R The vector at each support point which will hold the result from
-   *   the tensor contraction.
-   */
-  template< int NUM_SUPPORT_POINTS >
-  GEOS_HOST_DEVICE
-  static void plusNaFi( real64 const (&N)[NUM_SUPPORT_POINTS],
-                        real64 const (&forcingTerm_detJxW)[3],
-                        real64 ( &R )[NUM_SUPPORT_POINTS][3] );
-
-  /**
-   * @brief Inner product of each basis function gradient with a rank-2
-   *   symmetric tensor added to the product each shape function with a vector.
-   * @tparam NUM_SUPPORT_POINTS The number of support points for the element.
-   * @tparam GRADIENT_TYPE The type of the array object holding the shape
-   *   function gradients.
-   * @param gradN The basis function gradients at a point in the element.
-   * @param var_detJxW The rank-2 symmetric tensor at @p q scaled by J*W.
-   * @param N The shape function value at a predetermined coordinate in the element.
-   * @param forcingTerm_detJxW A vector scaled by detJxW
-   * @param R The vector at each support point which will hold the result from
-   *   the tensor contraction.
-   *
-   * \f[
-   * R_i = \sum_a^{nSupport} \left ( \frac{\partial N_a}{\partial X_j} var_{ij} + N_a f_i \right ),
-   * \f]
-   */
-  template< int NUM_SUPPORT_POINTS,
-            typename GRADIENT_TYPE >
-  GEOS_HOST_DEVICE
-  static void plusGradNajAijPlusNaFi( GRADIENT_TYPE const & gradN,
-                                      real64 const (&var_detJxW)[3][3],
-                                      real64 const (&N)[NUM_SUPPORT_POINTS],
-                                      real64 const (&forcingTerm_detJxW)[3],
-                                      real64 ( &R )[NUM_SUPPORT_POINTS][3] );
-
-  /**
-   * @brief Inner product of each basis function gradient with a rank-2
-   *   tensor added to the product each shape function with a vector.
-   * @copydoc plusGradNajAijPlusNaFi
-   */
-  template< int NUM_SUPPORT_POINTS,
-            typename GRADIENT_TYPE >
-  GEOS_HOST_DEVICE
-  static void plusGradNajAijPlusNaFi( GRADIENT_TYPE const & gradN,
-                                      real64 const (&var_detJxW)[6],
-                                      real64 const (&N)[NUM_SUPPORT_POINTS],
-                                      real64 const (&forcingTerm_detJxW)[3],
-                                      real64 ( &R )[NUM_SUPPORT_POINTS][3] );
-
-
-  /**
-   * @brief Sets m_viewGradN equal to an input view.
-   * @param source The view to assign to m_viewGradN.
-   */
-  void setGradNView( arrayView4d< real64 const > const & source )
-  {
-    GEOS_ERROR_IF_NE_MSG( source.size( 1 ),
-                          getNumQuadraturePoints(),
-                          "2nd-dimension of gradN array does not match number of quadrature points" );
-    GEOS_ERROR_IF_NE_MSG( source.size( 2 ),
-                          getMaxSupportPoints(),
-                          "3rd-dimension of gradN array does not match number of support points" );
-    GEOS_ERROR_IF_NE_MSG( source.size( 3 ),
-                          3,
-                          "4th-dimension of gradN array does not match 3" );
-
-    m_viewGradN = source;
-  }
-
-  /**
-   * @brief Sets m_viewDetJ equal to an input view.
-   * @param source The view to assign to m_viewDetJ.
-   */
-  void setDetJView( arrayView2d< real64 const > const & source )
-  {
-    GEOS_ERROR_IF_NE_MSG( source.size( 1 ),
-                          getNumQuadraturePoints(),
-                          "2nd-dimension of gradN array does not match number of quadrature points" );
-    m_viewDetJ = source;
-  }
 
-  /**
-   * @brief Getter for m_viewGradN
-   * @return A new arrayView copy of m_viewGradN.
-   */
-  arrayView4d< real64 const > getGradNView() const
-  {
-    return m_viewGradN;
-  }
-
-  /**
-   * @brief Getter for m_viewDetJ
-   * @return A new arrayView copy of m_viewDetJ.
-   */
-  arrayView2d< real64 const > getDetJView() const
-  {
-    return m_viewDetJ;
-  }
 
 
-protected:
-  /// View to potentially hold pre-calculated shape function gradients.
-  arrayView4d< real64 const > m_viewGradN;
 
-  /// View to potentially hold pre-calculated weighted jacobian transformation
-  /// determinants.
-  arrayView2d< real64 const > m_viewDetJ;
 };
 
 /// @cond Doxygen_Suppress
@@ -791,282 +370,11 @@ constexpr PDEUtilities::FunctionSpace FiniteElementBase::getFunctionSpace()
   return FunctionSpaceHelper< N >::getFunctionSpace();
 }
 
-template< typename LEAF >
-GEOS_HOST_DEVICE
-GEOS_FORCE_INLINE
-real64 FiniteElementBase::getGradN( localIndex const k,
-                                    localIndex const q,
-                                    real64 const (&X)[LEAF::maxSupportPoints][3],
-                                    real64 (& gradN)[LEAF::maxSupportPoints][3] ) const
-{
-  GEOS_UNUSED_VAR( k );
-  return LEAF::calcGradN( q, X, gradN );
-}
-
-template< typename LEAF >
-GEOS_HOST_DEVICE
-GEOS_FORCE_INLINE
-real64 FiniteElementBase::getGradN( localIndex const k,
-                                    localIndex const q,
-                                    real64 const (&X)[LEAF::maxSupportPoints][3],
-                                    typename LEAF::StackVariables const & stack,
-                                    real64 ( & gradN )[LEAF::maxSupportPoints][3] ) const
-{
-  GEOS_UNUSED_VAR( k );
-  return LEAF::calcGradN( q, X, stack, gradN );
-}
-
-template< typename LEAF >
-GEOS_HOST_DEVICE
-GEOS_FORCE_INLINE
-real64 FiniteElementBase::getGradN( localIndex const k,
-                                    localIndex const q,
-                                    int const X,
-                                    real64 (& gradN)[LEAF::maxSupportPoints][3] ) const
-{
-  GEOS_UNUSED_VAR( X );
-
-  LvArray::tensorOps::copy< LEAF::maxSupportPoints, 3 >( gradN, m_viewGradN[ k ][ q ] );
-
-  return m_viewDetJ( k, q );
-}
-
-template< typename LEAF >
-GEOS_HOST_DEVICE
-GEOS_FORCE_INLINE
-real64 FiniteElementBase::getGradN( localIndex const k,
-                                    localIndex const q,
-                                    int const X,
-                                    typename LEAF::StackVariables const & stack,
-                                    real64 (& gradN)[LEAF::maxSupportPoints][3] ) const
-{
-  GEOS_UNUSED_VAR( X );
-  GEOS_UNUSED_VAR( stack );
-
-  LvArray::tensorOps::copy< LEAF::maxSupportPoints, 3 >( gradN, m_viewGradN[ k ][ q ] );
-
-  return m_viewDetJ( k, q );
-}
-
-//*************************************************************************************************
-//***** Interpolated Value Functions **************************************************************
-//*************************************************************************************************
-
-template< int NUM_SUPPORT_POINTS >
-GEOS_HOST_DEVICE
-GEOS_FORCE_INLINE
-void FiniteElementBase::value( real64 const (&N)[NUM_SUPPORT_POINTS],
-                               real64 const (&var)[NUM_SUPPORT_POINTS],
-                               real64 & value )
-{
-  value = LvArray::tensorOps::AiBi< NUM_SUPPORT_POINTS >( N, var );
-}
-
-template< int NUM_SUPPORT_POINTS,
-          int NUM_COMPONENTS >
-GEOS_HOST_DEVICE
-GEOS_FORCE_INLINE
-void FiniteElementBase::value( real64 const (&N)[NUM_SUPPORT_POINTS],
-                               real64 const (&var)[NUM_SUPPORT_POINTS][NUM_COMPONENTS],
-                               real64 (& value)[NUM_COMPONENTS] )
-{
-
-  LvArray::tensorOps::Ri_eq_AjiBj< 3, NUM_SUPPORT_POINTS >( value, var, N );
-}
-
-
-//*************************************************************************************************
-//***** Variable Gradient Functions ***************************************************************
-//*************************************************************************************************
-
-template< int NUM_SUPPORT_POINTS,
-          typename GRADIENT_TYPE >
-GEOS_HOST_DEVICE
-GEOS_FORCE_INLINE
-void FiniteElementBase::symmetricGradient( GRADIENT_TYPE const & gradN,
-                                           real64 const (&var)[NUM_SUPPORT_POINTS][3],
-                                           real64 (& gradVar)[6] )
-{
-  gradVar[0] = gradN[0][0] * var[0][0];
-  gradVar[1] = gradN[0][1] * var[0][1];
-  gradVar[2] = gradN[0][2] * var[0][2];
-  gradVar[3] = gradN[0][2] * var[0][1] + gradN[0][1] * var[0][2];
-  gradVar[4] = gradN[0][2] * var[0][0] + gradN[0][0] * var[0][2];
-  gradVar[5] = gradN[0][1] * var[0][0] + gradN[0][0] * var[0][1];
-
-  for( int a=1; a<NUM_SUPPORT_POINTS; ++a )
-  {
-    gradVar[0] = gradVar[0] + gradN[a][0] * var[ a ][0];
-    gradVar[1] = gradVar[1] + gradN[a][1] * var[ a ][1];
-    gradVar[2] = gradVar[2] + gradN[a][2] * var[ a ][2];
-    gradVar[3] = gradVar[3] + gradN[a][2] * var[ a ][1] + gradN[a][1] * var[ a ][2];
-    gradVar[4] = gradVar[4] + gradN[a][2] * var[ a ][0] + gradN[a][0] * var[ a ][2];
-    gradVar[5] = gradVar[5] + gradN[a][1] * var[ a ][0] + gradN[a][0] * var[ a ][1];
-  }
-}
-
-template< int NUM_SUPPORT_POINTS,
-          typename GRADIENT_TYPE >
-GEOS_HOST_DEVICE
-GEOS_FORCE_INLINE
-real64 FiniteElementBase::symmetricGradientTrace( GRADIENT_TYPE const & gradN,
-                                                  real64 const (&var)[NUM_SUPPORT_POINTS][3] )
-{
-  real64 result = gradN[0][0] * var[0][0] + gradN[0][1] * var[0][1] + gradN[0][2] * var[0][2];
-
-  for( int a=1; a<NUM_SUPPORT_POINTS; ++a )
-  {
-    result = result + gradN[a][0] * var[a][0] + gradN[a][1] * var[a][1] + gradN[a][2] * var[a][2];
-  }
-  return result;
-}
-
-template< int NUM_SUPPORT_POINTS,
-          typename GRADIENT_TYPE >
-GEOS_HOST_DEVICE
-GEOS_FORCE_INLINE
-void FiniteElementBase::gradient( GRADIENT_TYPE const & gradN,
-                                  real64 const (&var)[NUM_SUPPORT_POINTS],
-                                  real64 (& gradVar)[3] )
-{
-  LvArray::tensorOps::Ri_eq_AjiBj< 3, NUM_SUPPORT_POINTS >( gradVar, gradN, var );
-}
-
-template< int NUM_SUPPORT_POINTS,
-          typename GRADIENT_TYPE >
-GEOS_HOST_DEVICE
-GEOS_FORCE_INLINE
-void FiniteElementBase::gradient( GRADIENT_TYPE const & gradN,
-                                  real64 const (&var)[NUM_SUPPORT_POINTS][3],
-                                  real64 (& gradVar)[3][3] )
-{
-  LvArray::tensorOps::Rij_eq_AkiBkj< 3, 3, NUM_SUPPORT_POINTS >( gradVar, var, gradN );
-}
-
-
-
-template< int NUM_SUPPORT_POINTS,
-          typename GRADIENT_TYPE >
-GEOS_HOST_DEVICE
-GEOS_FORCE_INLINE
-void FiniteElementBase::valueAndGradient( real64 const (&N)[NUM_SUPPORT_POINTS],
-                                          GRADIENT_TYPE const & gradN,
-                                          real64 const (&var)[NUM_SUPPORT_POINTS],
-                                          real64 & value,
-                                          real64 (& gradVar)[3] )
-{
-  value = N[0] * var[0];
-  for( int i = 0; i < 3; ++i )
-  {
-    gradVar[i] = var[0] * gradN[0][i];
-  }
-
-  for( int a=1; a<NUM_SUPPORT_POINTS; ++a )
-  {
-    value = value + N[a] * var[a];
-    for( int i = 0; i < 3; ++i )
-    {
-      gradVar[i] = gradVar[i] + var[ a ] * gradN[a][i];
-    }
-  }
-}
-
-
 
-template< int NUM_SUPPORT_POINTS,
-          typename GRADIENT_TYPE >
-GEOS_HOST_DEVICE
-GEOS_FORCE_INLINE
-void FiniteElementBase::plusGradNajAij( GRADIENT_TYPE const & gradN,
-                                        real64 const (&var_detJxW)[6],
-                                        real64 (& R)[NUM_SUPPORT_POINTS][3] )
-{
-  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
-  {
-    R[a][0] = R[a][0] + var_detJxW[0] * gradN[a][0] + var_detJxW[5] * gradN[a][1] + var_detJxW[4] * gradN[a][2];
-    R[a][1] = R[a][1] + var_detJxW[5] * gradN[a][0] + var_detJxW[1] * gradN[a][1] + var_detJxW[3] * gradN[a][2];
-    R[a][2] = R[a][2] + var_detJxW[4] * gradN[a][0] + var_detJxW[3] * gradN[a][1] + var_detJxW[2] * gradN[a][2];
-  }
-}
-
-
-template< int NUM_SUPPORT_POINTS,
-          typename GRADIENT_TYPE >
-GEOS_HOST_DEVICE
-GEOS_FORCE_INLINE
-void FiniteElementBase::plusGradNajAij( GRADIENT_TYPE const & gradN,
-                                        real64 const (&var_detJxW)[3][3],
-                                        real64 (& R)[NUM_SUPPORT_POINTS][3] )
-{
-  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
-  {
-    LvArray::tensorOps::Ri_add_AijBj< 3, 3 >( R[a], var_detJxW, gradN[a] );
-  }
-}
-
-template< int NUM_SUPPORT_POINTS >
-GEOS_HOST_DEVICE
-GEOS_FORCE_INLINE
-void FiniteElementBase::plusNaFi( real64 const (&N)[NUM_SUPPORT_POINTS],
-                                  real64 const (&var_detJxW)[3],
-                                  real64 ( & R )[NUM_SUPPORT_POINTS][3] )
-{
-  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
-  {
-    LvArray::tensorOps::scaledAdd< 3 >( R[a], var_detJxW, N[a] );
-  }
-}
-
-
-template< int NUM_SUPPORT_POINTS,
-          typename GRADIENT_TYPE >
-GEOS_HOST_DEVICE
-GEOS_FORCE_INLINE
-void FiniteElementBase::plusGradNajAijPlusNaFi( GRADIENT_TYPE const & gradN,
-                                                real64 const (&var_detJxW)[6],
-                                                real64 const (&N)[NUM_SUPPORT_POINTS],
-                                                real64 const (&forcingTerm_detJxW)[3],
-                                                real64 (& R)[NUM_SUPPORT_POINTS][3] )
-{
-  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
-  {
-    R[a][0] = R[a][0] + var_detJxW[0] * gradN[a][0] + var_detJxW[5] * gradN[a][1] + var_detJxW[4] * gradN[a][2] + forcingTerm_detJxW[0] * N[a];
-    R[a][1] = R[a][1] + var_detJxW[5] * gradN[a][0] + var_detJxW[1] * gradN[a][1] + var_detJxW[3] * gradN[a][2] + forcingTerm_detJxW[1] * N[a];
-    R[a][2] = R[a][2] + var_detJxW[4] * gradN[a][0] + var_detJxW[3] * gradN[a][1] + var_detJxW[2] * gradN[a][2] + forcingTerm_detJxW[2] * N[a];
-  }
-}
-
-template< int NUM_SUPPORT_POINTS,
-          typename GRADIENT_TYPE >
-GEOS_HOST_DEVICE
-GEOS_FORCE_INLINE
-void FiniteElementBase::plusGradNajAijPlusNaFi( GRADIENT_TYPE const & gradN,
-                                                real64 const (&var_detJxW)[3][3],
-                                                real64 const (&N)[NUM_SUPPORT_POINTS],
-                                                real64 const (&forcingTerm_detJxW)[3],
-                                                real64 (& R)[NUM_SUPPORT_POINTS][3] )
-{
-  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
-  {
-    R[a][0] = R[a][0] + var_detJxW[0][0] * gradN[a][0] + var_detJxW[0][1] * gradN[a][1] + var_detJxW[0][2] * gradN[a][2] + forcingTerm_detJxW[0] * N[a];
-    R[a][1] = R[a][1] + var_detJxW[1][0] * gradN[a][0] + var_detJxW[1][1] * gradN[a][1] + var_detJxW[1][2] * gradN[a][2] + forcingTerm_detJxW[1] * N[a];
-    R[a][2] = R[a][2] + var_detJxW[2][0] * gradN[a][0] + var_detJxW[2][1] * gradN[a][1] + var_detJxW[2][2] * gradN[a][2] + forcingTerm_detJxW[2] * N[a];
-  }
-}
 /// @endcond
 
-}
-}
-
 
-/// Macro to simplify name resolution in derived classes.
-#define USING_FINITEELEMENTBASE                       \
-  using FiniteElementBase::value;                     \
-  using FiniteElementBase::symmetricGradient;         \
-  using FiniteElementBase::gradient;                  \
-  using FiniteElementBase::valueAndGradient;          \
-  using FiniteElementBase::plusGradNajAij;           \
-  using FiniteElementBase::plusNaFi;                 \
-  using FiniteElementBase::plusGradNajAijPlusNaFi;
+} // namespace geos::finiteElement
+} // namespace geos
 
 #endif //GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_FINITEELEMENTBASE_HPP_
diff --git a/src/coreComponents/finiteElement/elementFormulations/FiniteElementOperators.hpp b/src/coreComponents/finiteElement/elementFormulations/FiniteElementOperators.hpp
new file mode 100644
index 00000000000..5b84b5d827b
--- /dev/null
+++ b/src/coreComponents/finiteElement/elementFormulations/FiniteElementOperators.hpp
@@ -0,0 +1,383 @@
+#pragma once
+
+#include "common/DataTypes.hpp"
+#include "LvArray/src/tensorOps.hpp"
+
+namespace geos
+{
+namespace finiteElement
+{
+namespace feOps
+{
+
+/**
+ * @name Value Operator Functions
+ */
+///@{
+
+/**
+ * @brief Compute the interpolated value of a variable.
+ * @tparam NUM_SUPPORT_POINTS The number of support points for the element.
+ * @param N Array (for each support point) of shape function values at the
+ *   coordinate the variable is to be interpolated.
+ * @param var Array of variable values for each support point.
+ * @param value The interpolated value of @p var.
+ *
+ * This is the standard finite element interpolation operator of a discrete
+ * variable defined at the support points.
+ * The operator is expressed as:
+ * \f[
+ * value  = \sum_a^{numSupport} \left ( N_a var_a \right ),
+ * \f]
+ * @note The shape function values @p N must be evaluated prior to calling this
+ * function.
+ *
+ */
+template< int NUM_SUPPORT_POINTS >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+void value( real64 const (&N)[NUM_SUPPORT_POINTS],
+                               real64 const (&var)[NUM_SUPPORT_POINTS],
+                               real64 & value )
+{
+  value = LvArray::tensorOps::AiBi< NUM_SUPPORT_POINTS >( N, var );
+}
+
+/**
+ * @brief Compute the interpolated value of a vector variable.
+ * @tparam NUM_COMPONENTS Number of components for the vector variable.
+ * @copydoc value
+ */
+template< int NUM_SUPPORT_POINTS,
+          int NUM_COMPONENTS >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+void value( real64 const (&N)[NUM_SUPPORT_POINTS],
+                               real64 const (&var)[NUM_SUPPORT_POINTS][NUM_COMPONENTS],
+                               real64 (& value)[NUM_COMPONENTS] )
+{
+  LvArray::tensorOps::Ri_eq_AjiBj< 3, NUM_SUPPORT_POINTS >( value, var, N );
+}
+
+///@}
+
+/**
+ * @name Gradient Operator Functions
+ */
+///@{
+
+/**
+ * @brief Calculate the symmetric gradient of a vector valued support field
+ *   at a point using the stored basis function gradients for all support
+ *   points.
+ * @tparam NUM_SUPPORT_POINTS The number of support points for the element.
+ * @tparam GRADIENT_TYPE The type of the array object holding the shape
+ * @param gradN The basis function gradients at a point in the element.
+ * @param var The vector valued support field that the gradient operator will
+ *  be applied to.
+ * @param gradVar The symmetric gradient in Voigt notation.
+ *
+ * More precisely, the operator is defined as:
+ * \f[
+ * grad^s_{ij}  = \frac{1}{2} \sum_a^{nSupport} \left ( \frac{\partial N_a}{\partial X_j} var_{ai} + \frac{\partial N_a}{\partial X_i}
+ * var_{aj}\right ),
+ * \f]
+ *
+ */
+template< int NUM_SUPPORT_POINTS,
+          typename GRADIENT_TYPE >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+void symmetricGradient( GRADIENT_TYPE const & gradN,
+                                           real64 const (&var)[NUM_SUPPORT_POINTS][3],
+                                           real64 (& gradVar)[6] )
+{
+  gradVar[0] = gradN[0][0] * var[0][0];
+  gradVar[1] = gradN[0][1] * var[0][1];
+  gradVar[2] = gradN[0][2] * var[0][2];
+  gradVar[3] = gradN[0][2] * var[0][1] + gradN[0][1] * var[0][2];
+  gradVar[4] = gradN[0][2] * var[0][0] + gradN[0][0] * var[0][2];
+  gradVar[5] = gradN[0][1] * var[0][0] + gradN[0][0] * var[0][1];
+
+  for( int a=1; a<NUM_SUPPORT_POINTS; ++a )
+  {
+    gradVar[0] = gradVar[0] + gradN[a][0] * var[ a ][0];
+    gradVar[1] = gradVar[1] + gradN[a][1] * var[ a ][1];
+    gradVar[2] = gradVar[2] + gradN[a][2] * var[ a ][2];
+    gradVar[3] = gradVar[3] + gradN[a][2] * var[ a ][1] + gradN[a][1] * var[ a ][2];
+    gradVar[4] = gradVar[4] + gradN[a][2] * var[ a ][0] + gradN[a][0] * var[ a ][2];
+    gradVar[5] = gradVar[5] + gradN[a][1] * var[ a ][0] + gradN[a][0] * var[ a ][1];
+  }
+}
+
+/**
+ * @brief Calculate the trace of the symmetric gradient of a vector valued support
+ *   field (i.e. the volumetric strain for the displacement field) at a point using
+ *   the stored basis function gradients for all support points.
+ * @tparam NUM_SUPPORT_POINTS The number of support points for the element.
+ * @tparam GRADIENT_TYPE The type of the array object holding the shape
+ * @param gradN The basis function gradients at a point in the element.
+ * @param var The vector valued support field that the gradient operator will
+ *  be applied to.
+ * @return The trace of the symetric gradient tensor.
+ *
+ */
+template< int NUM_SUPPORT_POINTS,
+          typename GRADIENT_TYPE >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+real64 symmetricGradientTrace( GRADIENT_TYPE const & gradN,
+                                                  real64 const (&var)[NUM_SUPPORT_POINTS][3] )
+{
+  real64 result = gradN[0][0] * var[0][0] + gradN[0][1] * var[0][1] + gradN[0][2] * var[0][2];
+
+  for( int a=1; a<NUM_SUPPORT_POINTS; ++a )
+  {
+    result = result + gradN[a][0] * var[a][0] + gradN[a][1] * var[a][1] + gradN[a][2] * var[a][2];
+  }
+  return result;
+}
+
+/**
+ * @brief Calculate the gradient of a scalar valued support field at a point
+ *   using the input basis function gradients.
+ * @tparam NUM_SUPPORT_POINTS The number of support points for the element.
+ * @tparam GRADIENT_TYPE The type of the array object holding the shape
+ *   function gradients.
+ * @param gradN The basis function gradients at a point in the element.
+ * @param var The vector valued support field that the gradient operator will
+ *  be applied to.
+ * @param gradVar The  gradient.
+ *
+ * More precisely, the operator is defined as:
+ * \f[
+ * grad_{j}  = \sum_a^{nSupport} \left ( \frac{\partial N_a}{\partial X_j} var_{a}\right ),
+ * \f]
+ */
+template< int NUM_SUPPORT_POINTS,
+          typename GRADIENT_TYPE >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+void gradient( GRADIENT_TYPE const & gradN,
+                                  real64 const (&var)[NUM_SUPPORT_POINTS],
+                                  real64 (& gradVar)[3] )
+{
+  LvArray::tensorOps::Ri_eq_AjiBj< 3, NUM_SUPPORT_POINTS >( gradVar, gradN, var );
+}
+
+/**
+ * @brief Calculate the gradient of a vector valued support field at a point
+ *   using the input basis function gradients.
+ * @copydoc gradient
+ *
+ * More precisely, the operator is defined as:
+ * \f[
+ * grad_{ij}  = \sum_a^{nSupport} \left ( \frac{\partial N_a}{\partial X_j} var_{ai}\right ),
+ * \f]
+ */
+template< int NUM_SUPPORT_POINTS,
+          typename GRADIENT_TYPE >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+void gradient( GRADIENT_TYPE const & gradN,
+                                  real64 const (&var)[NUM_SUPPORT_POINTS][3],
+                                  real64 (& gradVar)[3][3] )
+{
+  LvArray::tensorOps::Rij_eq_AkiBkj< 3, 3, NUM_SUPPORT_POINTS >( gradVar, var, gradN );
+}
+
+///@}
+
+/**
+ * @name Multi-Operator Functions
+ */
+///@{
+
+/**
+ * @brief Calculate the value and gradient of a scalar valued support field
+ *   at a point using the input basis function gradients.
+ * @tparam NUM_SUPPORT_POINTS The number of support points for the element.
+ * @tparam GRADIENT_TYPE The type of the array object holding the shape
+ * @param N Array (for each support point) of shape function values at the
+ *   coordinate the variable is to be interpolated.
+ * @param gradN The basis function gradients at a point in the element.
+ * @param var The vector valued support field that the gradient operator will
+ *  be applied to.
+ * @param value The value at the point for which N was specified.
+ * @param gradVar The gradient at the point for which gradN was specified.
+ */
+template< int NUM_SUPPORT_POINTS,
+          typename GRADIENT_TYPE >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+void valueAndGradient( real64 const (&N)[NUM_SUPPORT_POINTS],
+                                          GRADIENT_TYPE const & gradN,
+                                          real64 const (&var)[NUM_SUPPORT_POINTS],
+                                          real64 & value,
+                                          real64 (& gradVar)[3] )
+{
+  value = N[0] * var[0];
+  for( int i = 0; i < 3; ++i )
+  {
+    gradVar[i] = var[0] * gradN[0][i];
+  }
+
+  for( int a=1; a<NUM_SUPPORT_POINTS; ++a )
+  {
+    value = value + N[a] * var[a];
+    for( int i = 0; i < 3; ++i )
+    {
+      gradVar[i] = gradVar[i] + var[ a ] * gradN[a][i];
+    }
+  }
+}
+
+///@}
+
+/**
+ * @name Scattering Operator Functions
+ *
+ * These functions take quadrature data and map it to the support points
+ * through some operator.
+ */
+///@{
+
+/**
+ * @brief Inner product of each basis function gradient with a rank-2
+ *   symmetric tensor.
+ * @tparam NUM_SUPPORT_POINTS The number of support points for the element.
+ * @tparam GRADIENT_TYPE The type of the array object holding the shape
+ *   function gradients.
+ * @param gradN The basis function gradients at a point in the element.
+ * @param var_detJxW The rank-2 tensor at @p q scaled by J*W.
+ * @param R The vector at each support point which will hold the result from
+ *   the tensor contraction.
+ *
+ * More precisely, the operator is defined as:
+ *
+ * \f[
+ * R_i = \sum_a^{nSupport} \left( \frac{\partial N_a}{\partial X_j} var_{ij}\right),
+ * \f]
+ *
+ * where \f$\frac{\partial N_a}{\partial X_j}\f$ is the basis function gradient,
+ *   \f$var_{ij}\f$ is the rank-2 symmetric tensor.
+ */
+template< int NUM_SUPPORT_POINTS,
+          typename GRADIENT_TYPE >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+void plusGradNajAij( GRADIENT_TYPE const & gradN,
+                                        real64 const (&var_detJxW)[6],
+                                        real64 (& R)[NUM_SUPPORT_POINTS][3] )
+{
+  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
+  {
+    R[a][0] = R[a][0] + var_detJxW[0] * gradN[a][0] + var_detJxW[5] * gradN[a][1] + var_detJxW[4] * gradN[a][2];
+    R[a][1] = R[a][1] + var_detJxW[5] * gradN[a][0] + var_detJxW[1] * gradN[a][1] + var_detJxW[3] * gradN[a][2];
+    R[a][2] = R[a][2] + var_detJxW[4] * gradN[a][0] + var_detJxW[3] * gradN[a][1] + var_detJxW[2] * gradN[a][2];
+  }
+}
+
+/**
+ * @copydoc plusGradNajAij
+ * @brief Inner product of each basis function gradient with a rank-2
+ *   tensor.
+ */
+template< int NUM_SUPPORT_POINTS,
+          typename GRADIENT_TYPE >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+void plusGradNajAij( GRADIENT_TYPE const & gradN,
+                                        real64 const (&var_detJxW)[3][3],
+                                        real64 (& R)[NUM_SUPPORT_POINTS][3] )
+{
+  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
+  {
+    LvArray::tensorOps::Ri_add_AijBj< 3, 3 >( R[a], var_detJxW, gradN[a] );
+  }
+}
+
+/**
+ * @brief Product of each shape function with a vector forcing term.
+ * @tparam NUM_SUPPORT_POINTS The number of support points for the element.
+ * @param N The shape function value at a predetermined coordinate in the element.
+ * @param forcingTerm_detJxW A vector scaled by detJxW
+ * @param R The vector at each support point which will hold the result from
+ *   the tensor contraction.
+ */
+template< int NUM_SUPPORT_POINTS >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+void plusNaFi( real64 const (&N)[NUM_SUPPORT_POINTS],
+                                  real64 const (&var_detJxW)[3],
+                                  real64 ( & R )[NUM_SUPPORT_POINTS][3] )
+{
+  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
+  {
+    LvArray::tensorOps::scaledAdd< 3 >( R[a], var_detJxW, N[a] );
+  }
+}
+
+/**
+ * @brief Inner product of each basis function gradient with a rank-2
+ *   symmetric tensor added to the product each shape function with a vector.
+ * @tparam NUM_SUPPORT_POINTS The number of support points for the element.
+ * @tparam GRADIENT_TYPE The type of the array object holding the shape
+ *   function gradients.
+ * @param gradN The basis function gradients at a point in the element.
+ * @param var_detJxW The rank-2 symmetric tensor at @p q scaled by J*W.
+ * @param N The shape function value at a predetermined coordinate in the element.
+ * @param forcingTerm_detJxW A vector scaled by detJxW
+ * @param R The vector at each support point which will hold the result from
+ *   the tensor contraction.
+ *
+ * \f[
+ * R_i = \sum_a^{nSupport} \left ( \frac{\partial N_a}{\partial X_j} var_{ij} + N_a f_i \right ),
+ * \f]
+ */
+template< int NUM_SUPPORT_POINTS,
+          typename GRADIENT_TYPE >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+void plusGradNajAijPlusNaFi( GRADIENT_TYPE const & gradN,
+                                                real64 const (&var_detJxW)[6],
+                                                real64 const (&N)[NUM_SUPPORT_POINTS],
+                                                real64 const (&forcingTerm_detJxW)[3],
+                                                real64 (& R)[NUM_SUPPORT_POINTS][3] )
+{
+  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
+  {
+    R[a][0] = R[a][0] + var_detJxW[0] * gradN[a][0] + var_detJxW[5] * gradN[a][1] + var_detJxW[4] * gradN[a][2] + forcingTerm_detJxW[0] * N[a];
+    R[a][1] = R[a][1] + var_detJxW[5] * gradN[a][0] + var_detJxW[1] * gradN[a][1] + var_detJxW[3] * gradN[a][2] + forcingTerm_detJxW[1] * N[a];
+    R[a][2] = R[a][2] + var_detJxW[4] * gradN[a][0] + var_detJxW[3] * gradN[a][1] + var_detJxW[2] * gradN[a][2] + forcingTerm_detJxW[2] * N[a];
+  }
+}
+
+/**
+ * @brief Inner product of each basis function gradient with a rank-2
+ *   tensor added to the product each shape function with a vector.
+ * @copydoc plusGradNajAijPlusNaFi
+ */
+template< int NUM_SUPPORT_POINTS,
+          typename GRADIENT_TYPE >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+void plusGradNajAijPlusNaFi( GRADIENT_TYPE const & gradN,
+                                                real64 const (&var_detJxW)[3][3],
+                                                real64 const (&N)[NUM_SUPPORT_POINTS],
+                                                real64 const (&forcingTerm_detJxW)[3],
+                                                real64 (& R)[NUM_SUPPORT_POINTS][3] )
+{
+  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
+  {
+    R[a][0] = R[a][0] + var_detJxW[0][0] * gradN[a][0] + var_detJxW[0][1] * gradN[a][1] + var_detJxW[0][2] * gradN[a][2] + forcingTerm_detJxW[0] * N[a];
+    R[a][1] = R[a][1] + var_detJxW[1][0] * gradN[a][0] + var_detJxW[1][1] * gradN[a][1] + var_detJxW[1][2] * gradN[a][2] + forcingTerm_detJxW[1] * N[a];
+    R[a][2] = R[a][2] + var_detJxW[2][0] * gradN[a][0] + var_detJxW[2][1] * gradN[a][1] + var_detJxW[2][2] * gradN[a][2] + forcingTerm_detJxW[2] * N[a];
+  }
+}
+
+///@}
+
+} // namespace finiteElement::finiteElement::feOps
+} // namespace geos::finiteElement
+} // namespace geos
diff --git a/src/coreComponents/finiteElement/elementFormulations/H1_Hexahedron_Lagrange1_GaussLegendre2.hpp b/src/coreComponents/finiteElement/elementFormulations/H1_Hexahedron_Lagrange1_GaussLegendre2.hpp
index 0b164a174c2..49e56182214 100644
--- a/src/coreComponents/finiteElement/elementFormulations/H1_Hexahedron_Lagrange1_GaussLegendre2.hpp
+++ b/src/coreComponents/finiteElement/elementFormulations/H1_Hexahedron_Lagrange1_GaussLegendre2.hpp
@@ -89,10 +89,6 @@ class H1_Hexahedron_Lagrange1_GaussLegendre2 final : public FiniteElementBase
   /// The number of sampling points per element
   constexpr static int numSamplingPoints = numSamplingPointsPerDirection * numSamplingPointsPerDirection * numSamplingPointsPerDirection;
 
-  /** @cond Doxygen_Suppress */
-  USING_FINITEELEMENTBASE
-  /** @endcond Doxygen_Suppress */
-
   GEOS_HOST_DEVICE
   virtual ~H1_Hexahedron_Lagrange1_GaussLegendre2() override
   {}
diff --git a/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp b/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp
index 25a3f6c7783..8f4796aa83c 100644
--- a/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp
+++ b/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp
@@ -118,10 +118,6 @@ class Qk_Hexahedron_Lagrange_GaussLobatto final : public FiniteElementBase
                              ( num1dNodes - 1 ) * ( k / 2 ) );
   }
 
-  /** @cond Doxygen_Suppress */
-  USING_FINITEELEMENTBASE
-  /** @endcond Doxygen_Suppress */
-
   ~Qk_Hexahedron_Lagrange_GaussLobatto() = default;
 
   GEOS_HOST_DEVICE
diff --git a/src/coreComponents/finiteElement/kernelInterface/ImplicitKernelBase.hpp b/src/coreComponents/finiteElement/kernelInterface/ImplicitKernelBase.hpp
index 8e86d076e68..cb6a2369484 100644
--- a/src/coreComponents/finiteElement/kernelInterface/ImplicitKernelBase.hpp
+++ b/src/coreComponents/finiteElement/kernelInterface/ImplicitKernelBase.hpp
@@ -184,7 +184,7 @@ class ImplicitKernelBase : public KernelBase< SUBREGION_TYPE,
               StackVariables & stack ) const
   {
     m_finiteElementSpace.template setup< FE_TYPE >( k, m_meshData, stack.feStack );
-    localIndex numTestSupportPoints = m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
+    localIndex numTestSupportPoints = m_finiteElementSpace.getNumSupportPoints( stack.feStack );
     localIndex numTrialSupportPoints = numTestSupportPoints;
     stack.numRows = numTestSupportPoints * numDofPerTestSupportPoint;
     stack.numCols = numTrialSupportPoints * numDofPerTrialSupportPoint;
diff --git a/src/coreComponents/finiteElement/unitTests/testFiniteElementBase.cpp b/src/coreComponents/finiteElement/unitTests/testFiniteElementBase.cpp
index b6e3c1f4610..9d66c3be10d 100644
--- a/src/coreComponents/finiteElement/unitTests/testFiniteElementBase.cpp
+++ b/src/coreComponents/finiteElement/unitTests/testFiniteElementBase.cpp
@@ -15,6 +15,7 @@
 
 
 #include "finiteElement/elementFormulations/FiniteElementBase.hpp"
+#include "finiteElement/elementFormulations/FiniteElementOperators.hpp"
 #include "gtest/gtest.h"
 #include "testFiniteElementHelpers.hpp"
 #include "common/GEOS_RAJA_Interface.hpp"
@@ -51,123 +52,7 @@ class TestFiniteElementBase final : public FiniteElementBase
   {}
 };
 
-TEST( FiniteElementBase, test_setGradNView )
-{
-  TestFiniteElementBase feBase;
-  {
-    array4d< real64 > gradN( 2, 8, 8, 3 );
-    feBase.setGradNView( gradN.toViewConst() );
-
-    EXPECT_EQ( feBase.getGradNView().size( 0 ), gradN.size( 0 ) );
-    EXPECT_EQ( feBase.getGradNView().size( 1 ), gradN.size( 1 ) );
-    EXPECT_EQ( feBase.getGradNView().size( 2 ), gradN.size( 2 ) );
-    EXPECT_EQ( feBase.getGradNView().size( 3 ), gradN.size( 3 ) );
-  }
-
-  {
-    array4d< real64 > gradN( 2, 7, 8, 3 );
-    EXPECT_DEATH_IF_SUPPORTED( feBase.setGradNView( gradN.toViewConst() ), "" );
-  }
-
-  {
-    array4d< real64 > gradN( 2, 8, 7, 3 );
-    EXPECT_DEATH_IF_SUPPORTED( feBase.setGradNView( gradN.toViewConst() ), "" );
-  }
-
-  {
-    array4d< real64 > gradN( 2, 8, 8, 2 );
-    EXPECT_DEATH_IF_SUPPORTED( feBase.setGradNView( gradN.toViewConst() ), "" );
-  }
-}
-
-TEST( FiniteElementBase, test_setDetJView )
-{
-  TestFiniteElementBase feBase;
-  {
-    array2d< real64 > detJ( 4, 8 );
-    feBase.setDetJView( detJ.toViewConst() );
-    EXPECT_EQ( feBase.getDetJView().size( 0 ), detJ.size( 0 ) );
-    EXPECT_EQ( feBase.getDetJView().size( 1 ), detJ.size( 1 ) );
-  }
 
-  {
-    array2d< real64 > detJ( 4, 7 );
-    EXPECT_DEATH_IF_SUPPORTED( feBase.setDetJView( detJ.toViewConst() ), "" );
-  }
-}
-
-//***** TEST getGradN *****************************************************************************
-TEST( FiniteElementBase, test_capture )
-{
-  TestFiniteElementBase feBase;
-  array4d< real64 > gradN( 4, 8, 8, 3 );
-  array2d< real64 > detJ( 4, 8 );
-  feBase.setGradNView( gradN.toViewConst() );
-  feBase.setDetJView( detJ.toViewConst() );
-
-
-  array1d< localIndex > gradNDims( 4 );
-  array1d< localIndex > detJDims( 2 );
-  arrayView1d< localIndex > gradNDimsView = gradNDims.toView();
-  arrayView1d< localIndex > detJDimsView = detJDims.toView();
-
-#if defined(CALC_FEM_SHAPE_IN_KERNEL)
-
-  forAll< parallelDevicePolicy<> >( 1, [ feBase, gradNDimsView, detJDimsView ]( int const i )
-  {
-    gradNDimsView[0] = feBase.getGradNView().size( 0 );
-    gradNDimsView[1] = feBase.getGradNView().size( 1 );
-    gradNDimsView[2] = feBase.getGradNView().size( 2 );
-    gradNDimsView[3] = feBase.getGradNView().size( 3 );
-    detJDimsView[0] = feBase.getDetJView().size( 0 );
-    detJDimsView[1] = feBase.getDetJView().size( 1 );
-
-    printf( "gradNDimsView = { %ld, %ld, %ld, %ld }\n",
-            gradNDimsView[0],
-            gradNDimsView[1],
-            gradNDimsView[2],
-            gradNDimsView[3] );
-
-    printf( "detJDimsView = { %ld, %ld }\n",
-            detJDimsView[0],
-            detJDimsView[1] );
-  } );
-
-  forAll< serialPolicy >( 1, [ feBase, gradNDimsView, detJDimsView ]( int const i )
-  {} );
-
-  EXPECT_EQ( gradNDimsView[0], 0 );
-  EXPECT_EQ( gradNDimsView[1], 0 );
-  EXPECT_EQ( gradNDimsView[2], 0 );
-  EXPECT_EQ( gradNDimsView[3], 0 );
-
-  EXPECT_EQ( detJDimsView[0], 0 );
-  EXPECT_EQ( detJDimsView[1], 0 );
-#endif
-
-
-  forAll< serialPolicy >( 1, [ feBase, gradNDimsView, detJDimsView ]( int const )
-  {
-    gradNDimsView[0] = feBase.getGradNView().size( 0 );
-    gradNDimsView[1] = feBase.getGradNView().size( 1 );
-    gradNDimsView[2] = feBase.getGradNView().size( 2 );
-    gradNDimsView[3] = feBase.getGradNView().size( 3 );
-    detJDimsView[0] = feBase.getDetJView().size( 0 );
-    detJDimsView[1] = feBase.getDetJView().size( 1 );
-  } );
-
-
-  EXPECT_EQ( gradNDimsView[0], gradN.size( 0 ) );
-  EXPECT_EQ( gradNDimsView[1], gradN.size( 1 ) );
-  EXPECT_EQ( gradNDimsView[2], gradN.size( 2 ) );
-  EXPECT_EQ( gradNDimsView[3], gradN.size( 3 ) );
-
-  EXPECT_EQ( detJDimsView[0], detJ.size( 0 ) );
-  EXPECT_EQ( detJDimsView[1], detJ.size( 1 ) );
-
-
-
-}
 
 //***** TEST value() ******************************************************************************
 
@@ -219,7 +104,7 @@ TEST( FiniteElementBase, test_value )
     real64 referenceScalarValue = -1.0;
     real64 feBaseScalarValue = -1.0;
     value( N, scalar, referenceScalarValue );
-    FiniteElementBase::value( N, scalar, feBaseScalarValue );
+    finiteElement::feOps::value( N, scalar, feBaseScalarValue );
 
     EXPECT_FLOAT_EQ( feBaseScalarValue, referenceScalarValue );
 
@@ -227,7 +112,7 @@ TEST( FiniteElementBase, test_value )
     real64 referenceVectorValue[3] = {-1, -1, -1};
     real64 feBaseVectorValue[3] = {-1, -1, -1};
     value( N, vector, referenceVectorValue );
-    FiniteElementBase::value( N, vector, feBaseVectorValue );
+    finiteElement::feOps::value( N, vector, feBaseVectorValue );
 
     EXPECT_FLOAT_EQ( feBaseVectorValue[0], referenceVectorValue[0] );
     EXPECT_FLOAT_EQ( feBaseVectorValue[1], referenceVectorValue[1] );
@@ -274,7 +159,7 @@ TEST( FiniteElementBase, test_symmetricGradient )
     real64 referenceVectorGradient[6] = {-1, -1, -1, -1, -1, -1};
     real64 feBaseVectorGradient[6] = {-1, -1, -1, -1, -1, -1};
     symmetricGradient( gradN, vector, referenceVectorGradient );
-    FiniteElementBase::symmetricGradient( gradN, vector, feBaseVectorGradient );
+    finiteElement::feOps::symmetricGradient( gradN, vector, feBaseVectorGradient );
 
     for( int i=0; i<6; ++i )
     {
@@ -344,13 +229,13 @@ TEST( FiniteElementBase, test_gradient )
     real64 referenceScalarGradient[3] = {-1, -1, -1};
     real64 feBaseScalarGradient[3] = {-1, -1, -1};
     gradient( gradN, scalar, referenceScalarGradient );
-    FiniteElementBase::gradient( gradN, scalar, feBaseScalarGradient );
+    finiteElement::feOps::gradient( gradN, scalar, feBaseScalarGradient );
 
 
     real64 referenceVectorGradient[3][3] = {{-1, -1, -1}, {-1, -1, -1}, {-1, -1, -1}};
     real64 feBaseVectorGradient[3][3] = {{-1, -1, -1}, {-1, -1, -1}, {-1, -1, -1}};
     gradient( gradN, vector, referenceVectorGradient );
-    FiniteElementBase::gradient( gradN, vector, feBaseVectorGradient );
+    finiteElement::feOps::gradient( gradN, vector, feBaseVectorGradient );
 
     for( int i=0; i<3; ++i )
     {
@@ -411,7 +296,7 @@ TEST( FiniteElementBase, test_valueAndGradient )
     real64 referenceScalarGradient[3] = {-1, -1, -1};
     real64 feBaseScalarGradient[3] = {-1, -1, -1};
     valueAndGradient( N, gradN, scalar, referenceScalarValue, referenceScalarGradient );
-    FiniteElementBase::valueAndGradient( N, gradN, scalar, feBaseScalarValue, feBaseScalarGradient );
+    finiteElement::feOps::valueAndGradient( N, gradN, scalar, feBaseScalarValue, feBaseScalarGradient );
 
     EXPECT_FLOAT_EQ( feBaseScalarValue, referenceScalarValue );
     for( int i=0; i<3; ++i )
@@ -425,7 +310,7 @@ TEST( FiniteElementBase, test_valueAndGradient )
 //    real64 referenceVectorGradient[3][3] = {{-1,-1,-1},{-1,-1,-1},{-1,-1,-1}};
 //    real64 feBaseVectorGradient[3][3] = {{-1,-1,-1},{-1,-1,-1},{-1,-1,-1}};
 //    valueAndGradient( N, gradN, vector, referenceVectorValue, referenceVectorGradient );
-//    FiniteElementBase::valueAndGradient( N, gradN, vector, feBaseVectorValue, feBaseVectorGradient );
+//    finiteElement::feOps::valueAndGradient( N, gradN, vector, feBaseVectorValue, feBaseVectorGradient );
 //
 //
 //    for( int i=0; i<3; ++i )
@@ -495,10 +380,10 @@ TEST( FiniteElementBase, test_plusGradNajAij )
     randomVar( r2SymmTensor );
 
     plusGradNajAij( gradN, r2Tensor, baselineResult );
-    FiniteElementBase::plusGradNajAij( gradN, r2Tensor, feResult );
+    finiteElement::feOps::plusGradNajAij( gradN, r2Tensor, feResult );
 
     plusGradNajAij( gradN, r2SymmTensor, baselineResultSym );
-    FiniteElementBase::plusGradNajAij( gradN, r2SymmTensor, feResultSym );
+    finiteElement::feOps::plusGradNajAij( gradN, r2SymmTensor, feResultSym );
   }
 
   for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
@@ -542,7 +427,7 @@ TEST( FiniteElementBase, test_plusNaFi )
     randomVar( f );
 
     plusNaFi( N, f, baselineResult );
-    FiniteElementBase::plusNaFi( N, f, feResult );
+    finiteElement::feOps::plusNaFi( N, f, feResult );
   }
 
   for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
@@ -614,10 +499,10 @@ TEST( FiniteElementBase, test_plusGradNajAijPlusNaFi )
     randomVar( f );
 
     plusGradNajAijPlusNaFi( gradN, r2Tensor, N, f, baselineResult );
-    FiniteElementBase::plusGradNajAijPlusNaFi( gradN, r2Tensor, N, f, feResult );
+    finiteElement::feOps::plusGradNajAijPlusNaFi( gradN, r2Tensor, N, f, feResult );
 
     plusGradNajAijPlusNaFi( gradN, r2SymmTensor, N, f, baselineResultSym );
-    FiniteElementBase::plusGradNajAijPlusNaFi( gradN, r2SymmTensor, N, f, feResultSym );
+    finiteElement::feOps::plusGradNajAijPlusNaFi( gradN, r2SymmTensor, N, f, feResultSym );
   }
 
   for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
diff --git a/src/coreComponents/mainInterface/ProblemManager.cpp b/src/coreComponents/mainInterface/ProblemManager.cpp
index e979d815d71..120308b4d02 100644
--- a/src/coreComponents/mainInterface/ProblemManager.cpp
+++ b/src/coreComponents/mainInterface/ProblemManager.cpp
@@ -935,7 +935,6 @@ map< std::tuple< string, string, string, string >, localIndex > ProblemManager::
         ElementRegionManager & elemManager = meshLevel.getElemManager();
         FaceManager const & faceManager = meshLevel.getFaceManager();
         EdgeManager const & edgeManager = meshLevel.getEdgeManager();
-        arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & X = nodeManager.referencePosition();
 
         for( auto const & regionName : regionNames )
         {
@@ -987,9 +986,6 @@ map< std::tuple< string, string, string, string >, localIndex > ProblemManager::
 
                   localIndex const numQuadraturePoints = FE_TYPE::numQuadraturePoints;
 
-//#if ! defined( CALC_FEM_SHAPE_IN_KERNEL )
-                  feDiscretization->calculateShapeFunctionGradients< SUBREGION_TYPE, FE_TYPE >( X, &subRegion, meshData, finiteElement );
-//#endif
 
                   localIndex & numQuadraturePointsInList = regionQuadrature[ std::make_tuple( meshBodyName,
                                                                                               meshLevel.getName(),
diff --git a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.hpp
index 359f428a3c3..75de5e22b49 100644
--- a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.hpp
+++ b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.hpp
@@ -144,7 +144,7 @@ struct DamageAndDamageGradientInterpolationKernel
 
         real64 qDamage = 0.0;
         real64 qDamageGrad[3] = {0, 0, 0};
-        FE_TYPE::valueAndGradient( N, dNdX, nodalDamageLocal, qDamage, qDamageGrad );
+        finiteElement::feOps::valueAndGradient( N, dNdX, nodalDamageLocal, qDamage, qDamageGrad );
 
         damageFieldOnMaterial( k, q ) = qDamage;
 
diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics_impl.hpp
index b5f2621b165..ee4a1888fce 100644
--- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics_impl.hpp
+++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics_impl.hpp
@@ -23,6 +23,7 @@
 #include "constitutive/fluid/multifluid/MultiFluidBase.hpp"
 #include "finiteElement/BilinearFormUtilities.hpp"
 #include "finiteElement/LinearFormUtilities.hpp"
+#include "finiteElement/elementFormulations/FiniteElementOperators.hpp"
 #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp"
 #include "physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp"
 #include "physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics.hpp"
@@ -617,12 +618,12 @@ quadraturePointKernel( localIndex const k,
   real64 N[numNodesPerElem]{};
   real64 dNdX[numNodesPerElem][3]{};
   FE_TYPE::calcN( q, stack.feStack, N );
-  real64 const detJxW = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal,
+  real64 const detJxW = FE_TYPE::calcGradN( q, stack.xLocal,
                                                                            stack.feStack, dNdX );
 
   // Step 2: compute strain increment
   LvArray::tensorOps::fill< 6 >( stack.strainIncrement, 0.0 );
-  FE_TYPE::symmetricGradient( dNdX, stack.uhat_local, stack.strainIncrement );
+  finiteElement::feOps::symmetricGradient( dNdX, stack.uhat_local, stack.strainIncrement );
 
   // Step 3: compute 1) the total stress, 2) the body force terms, and 3) the fluidMassIncrement
   // using quantities returned by the PorousSolid constitutive model.
@@ -656,7 +657,7 @@ complete( localIndex const k,
 
   real64 maxForce = 0;
   localIndex const numSupportPoints =
-    m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
+    m_finiteElementSpace.getNumSupportPoints( stack.feStack );
   integer numDisplacementDofs = numSupportPoints * numDofPerTestSupportPoint;
   constexpr integer maxNumDisplacementDofs = FE_TYPE::maxSupportPoints * numDofPerTestSupportPoint;
 
diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/PoromechanicsBase.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/PoromechanicsBase.hpp
index a71ed7ec34f..25b3e749bbb 100644
--- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/PoromechanicsBase.hpp
+++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/PoromechanicsBase.hpp
@@ -150,13 +150,8 @@ class PoromechanicsBase :
       localPressureDofIndex{ 0 }
     {}
 
-#if !defined(CALC_FEM_SHAPE_IN_KERNEL)
-    /// Dummy
-    int xLocal;
-#else
     /// C-array stack storage for element local the nodal positions.
     real64 xLocal[numNodesPerElem][3];
-#endif
 
     // Storage for displacements
 
@@ -224,7 +219,7 @@ class PoromechanicsBase :
 
     m_finiteElementSpace.template setup< FE_TYPE >( k, m_meshData, stack.feStack );
     localIndex const numSupportPoints =
-      m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
+      m_finiteElementSpace.getNumSupportPoints( stack.feStack );
 
     for( localIndex a=0; a<numSupportPoints; ++a )
     {
@@ -232,9 +227,7 @@ class PoromechanicsBase :
 
       for( integer i = 0; i < numDims; ++i )
       {
-#if defined(CALC_FEM_SHAPE_IN_KERNEL)
         stack.xLocal[a][i] = m_X[localNodeIndex][i];
-#endif
         stack.u_local[a][i] = m_disp[localNodeIndex][i];
         stack.uhat_local[a][i] = m_uhat[localNodeIndex][i];
         stack.localRowDofIndex[a*numDims+i] = m_dofNumber[localNodeIndex]+i;
diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsDamage_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsDamage_impl.hpp
index 67250f7b36d..c5da60ac797 100644
--- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsDamage_impl.hpp
+++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsDamage_impl.hpp
@@ -260,12 +260,12 @@ quadraturePointKernel( localIndex const k,
   real64 N[numNodesPerElem]{};
   real64 dNdX[numNodesPerElem][3]{};
   FE_TYPE::calcN( q, stack.feStack, N );
-  real64 const detJxW = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal,
+  real64 const detJxW = FE_TYPE::calcGradN( q, stack.xLocal,
                                                                            stack.feStack, dNdX );
 
   // Step 2: compute strain increment
   LvArray::tensorOps::fill< 6 >( stack.strainIncrement, 0.0 );
-  FE_TYPE::symmetricGradient( dNdX, stack.uhat_local, stack.strainIncrement );
+  finiteElement::feOps::symmetricGradient( dNdX, stack.uhat_local, stack.strainIncrement );
 
   // Step 3: compute 1) the total stress, 2) the body force terms, and 3) the fluidMassIncrement
   // using quantities returned by the PorousSolid constitutive model.
diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp
index 06887b4ce6d..9a2981a2636 100644
--- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp
+++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp
@@ -199,7 +199,7 @@ quadraturePointKernel( localIndex const k,
   // derivatives (dNdX), and (iii) determinant of the Jacobian transformation
   // matrix times the quadrature weight (detJxW)
   real64 dNdX[numNodesPerElem][3];
-  real64 const detJ = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal, dNdX );
+  real64 const detJ = FE_TYPE::calcGradN( q, stack.xLocal, dNdX );
 
   // EFEM part starts here
   constexpr int nUdof = numNodesPerElem*3;
diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics_impl.hpp
index 158372d2e3f..5de66273486 100644
--- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics_impl.hpp
+++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics_impl.hpp
@@ -22,6 +22,7 @@
 
 #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp"
 #include "finiteElement/BilinearFormUtilities.hpp"
+#include "finiteElement/elementFormulations/FiniteElementOperators.hpp"
 #include "finiteElement/LinearFormUtilities.hpp"
 #include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp"
 
@@ -365,12 +366,12 @@ quadraturePointKernel( localIndex const k,
   real64 N[numNodesPerElem]{};
   real64 dNdX[numNodesPerElem][3]{};
   FE_TYPE::calcN( q, stack.feStack, N );
-  real64 const detJxW = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal,
+  real64 const detJxW = FE_TYPE::calcGradN( q, stack.xLocal,
                                                                            stack.feStack, dNdX );
 
   // Step 2: compute strain increment
   LvArray::tensorOps::fill< 6 >( stack.strainIncrement, 0.0 );
-  FE_TYPE::symmetricGradient( dNdX, stack.uhat_local, stack.strainIncrement );
+  finiteElement::feOps::symmetricGradient( dNdX, stack.uhat_local, stack.strainIncrement );
 
   // Step 3: compute 1) the total stress, 2) the body force terms, and 3) the fluidMassIncrement
   // using quantities returned by the PorousSolid constitutive model.
@@ -398,7 +399,7 @@ complete( localIndex const k,
   GEOS_UNUSED_VAR( k );
   real64 maxForce = 0;
   localIndex const numSupportPoints =
-    m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
+    m_finiteElementSpace.getNumSupportPoints( stack.feStack );
   integer numDisplacementDofs = numSupportPoints * numDofPerTestSupportPoint;
 
   for( int localNode = 0; localNode < numSupportPoints; ++localNode )
diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics_impl.hpp
index 98e5a37075f..2be4e17296b 100644
--- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics_impl.hpp
+++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics_impl.hpp
@@ -493,12 +493,12 @@ quadraturePointKernel( localIndex const k,
   real64 N[numNodesPerElem]{};
   real64 dNdX[numNodesPerElem][3]{};
   FE_TYPE::calcN( q, stack.feStack, N );
-  real64 const detJxW = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal,
+  real64 const detJxW = FE_TYPE::calcGradN( q, stack.xLocal,
                                                                            stack.feStack, dNdX );
 
   // Step 2: compute strain increment
   LvArray::tensorOps::fill< 6 >( stack.strainIncrement, 0.0 );
-  FE_TYPE::symmetricGradient( dNdX, stack.uhat_local, stack.strainIncrement );
+  finiteElement::feOps::symmetricGradient( dNdX, stack.uhat_local, stack.strainIncrement );
 
   // Step 3: compute 1) the total stress, 2) the body force terms, and 3) the fluidMassIncrement
   // using quantities returned by the PorousSolid constitutive model.
@@ -528,7 +528,7 @@ complete( localIndex const k,
   real64 const maxForce = Base::complete( k, stack );
 
   localIndex const numSupportPoints =
-    m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
+    m_finiteElementSpace.getNumSupportPoints( stack.feStack );
   integer numDisplacementDofs = numSupportPoints * numDofPerTestSupportPoint;
 
   if( m_useTotalMassEquation > 0 )
diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics_impl.hpp
index eede62e0ca2..2baa8d3e1e4 100644
--- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics_impl.hpp
+++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics_impl.hpp
@@ -362,12 +362,12 @@ quadraturePointKernel( localIndex const k,
   real64 N[numNodesPerElem]{};
   real64 dNdX[numNodesPerElem][3]{};
   FE_TYPE::calcN( q, stack.feStack, N );
-  real64 const detJxW = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal,
+  real64 const detJxW = FE_TYPE::calcGradN( q, stack.xLocal,
                                                                            stack.feStack, dNdX );
 
   // Step 2: compute strain increment
   LvArray::tensorOps::fill< 6 >( stack.strainIncrement, 0.0 );
-  FE_TYPE::symmetricGradient( dNdX, stack.uhat_local, stack.strainIncrement );
+  finiteElement::feOps::symmetricGradient( dNdX, stack.uhat_local, stack.strainIncrement );
 
   // Step 3: compute 1) the total stress, 2) the body force terms, and 3) the fluidMassIncrement
   // using quantities returned by the PorousSolid constitutive model.
@@ -395,7 +395,7 @@ complete( localIndex const k,
   real64 const maxForce = Base::complete( k, stack );
 
   localIndex const numSupportPoints =
-    m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
+    m_finiteElementSpace.getNumSupportPoints( stack.feStack );
   integer numDisplacementDofs = numSupportPoints * numDofPerTestSupportPoint;
 
   // Step 1: assemble the derivatives of linear momentum balance wrt temperature into the global matrix
diff --git a/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEMKernels.hpp b/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEMKernels.hpp
index 6d992ba0134..3f8db629455 100644
--- a/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEMKernels.hpp
+++ b/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEMKernels.hpp
@@ -138,13 +138,8 @@ class LaplaceFEMKernel :
             primaryField_local{ 0.0 }
     {}
 
-#if !defined(CALC_FEM_SHAPE_IN_KERNEL)
-    /// Dummy
-    int xLocal;
-#else
     /// C-array stack storage for element local the nodal positions.
     real64 xLocal[ maxNumTestSupportPointsPerElem ][ 3 ];
-#endif
 
     /// C-array storage for the element local primary field variable.
     real64 primaryField_local[ maxNumTestSupportPointsPerElem ];
@@ -165,18 +160,16 @@ class LaplaceFEMKernel :
               StackVariables & stack ) const
   {
     m_finiteElementSpace.template setup< FE_TYPE >( k, m_meshData, stack.feStack );
-    stack.numRows = m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
+    stack.numRows = m_finiteElementSpace.getNumSupportPoints( stack.feStack );
     stack.numCols = stack.numRows;
     for( localIndex a = 0; a < stack.numRows; ++a )
     {
       localIndex const localNodeIndex = m_elemsToNodes( k, a );
 
-#if defined(CALC_FEM_SHAPE_IN_KERNEL)
       for( int i=0; i<3; ++i )
       {
         stack.xLocal[ a ][ i ] = m_X[ localNodeIndex ][ i ];
       }
-#endif
 
       stack.primaryField_local[ a ] = m_primaryField[ localNodeIndex ];
       stack.localRowDofIndex[a] = m_dofNumber[localNodeIndex];
@@ -196,9 +189,9 @@ class LaplaceFEMKernel :
                               localIndex const q,
                               StackVariables & stack ) const
   {
+    GEOS_UNUSED_VAR( k );
     real64 dNdX[ maxNumTestSupportPointsPerElem ][ 3 ];
-    real64 const detJ = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal,
-                                                                           stack.feStack, dNdX );
+    real64 const detJ = FE_TYPE::calcGradN( q, stack.xLocal, stack.feStack, dNdX );
     for( localIndex a = 0; a < stack.numRows; ++a )
     {
       for( localIndex b = 0; b < stack.numCols; ++b )
diff --git a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEMKernels.hpp b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEMKernels.hpp
index 32768c1074d..8bca1266398 100644
--- a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEMKernels.hpp
+++ b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEMKernels.hpp
@@ -21,6 +21,7 @@
 #define GEOS_PHYSICSSOLVERS_SIMPLEPDE_PHASEFIELDDAMAGEKERNELS_HPP_
 
 #include "finiteElement/kernelInterface/ImplicitKernelBase.hpp"
+#include "finiteElement/elementFormulations/FiniteElementOperators.hpp"
 
 namespace geos
 {
@@ -152,13 +153,8 @@ class PhaseFieldDamageKernel :
             nodalDamageLocal{ 0.0 }
     {}
 
-#if !defined(CALC_FEM_SHAPE_IN_KERNEL)
-    /// Dummy
-    int xLocal;
-#else
     /// C-array stack storage for element local the nodal positions.
     real64 xLocal[ numNodesPerElem ][ 3 ];
-#endif
 
     /// C-array storage for the element local primary field variable.
     real64 nodalDamageLocal[numNodesPerElem];
@@ -182,9 +178,7 @@ class PhaseFieldDamageKernel :
     {
       localIndex const localNodeIndex = m_elemsToNodes( k, a );
 
-#if defined(CALC_FEM_SHAPE_IN_KERNEL)
       LvArray::tensorOps::copy< 3 >( stack.xLocal[ a ], m_X[ localNodeIndex ] );
-#endif
 
       stack.nodalDamageLocal[ a ] = m_nodalDamage[ localNodeIndex ];
       stack.localRowDofIndex[a] = m_dofNumber[localNodeIndex];
@@ -210,12 +204,12 @@ class PhaseFieldDamageKernel :
     //Interpolate d and grad_d
     real64 N[ numNodesPerElem ];
     real64 dNdX[ numNodesPerElem ][ 3 ];
-    real64 const detJ = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal, dNdX );
+    real64 const detJ = FE_TYPE::calcGradN( q, stack.xLocal, dNdX );
     FE_TYPE::calcN( q, N );
 
     real64 qp_damage = 0.0;
     real64 qp_grad_damage[3] = {0, 0, 0};
-    FE_TYPE::valueAndGradient( N, dNdX, stack.nodalDamageLocal, qp_damage, qp_grad_damage );
+    finiteElement::feOps::valueAndGradient( N, dNdX, stack.nodalDamageLocal, qp_damage, qp_grad_damage );
 
     real64 D = 0;                                                                   //max between threshold and
                                                                                     // Elastic energy
diff --git a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldPressurizedDamageFEMKernels.hpp b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldPressurizedDamageFEMKernels.hpp
index 44747dd0546..0623667d56c 100644
--- a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldPressurizedDamageFEMKernels.hpp
+++ b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldPressurizedDamageFEMKernels.hpp
@@ -172,14 +172,14 @@ class PhaseFieldPressurizedDamageKernel :
     //Interpolate d and grad_d
     real64 N[ numNodesPerElem ];
     real64 dNdX[ numNodesPerElem ][ 3 ];
-    real64 const detJ = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal, dNdX );
+    real64 const detJ = FE_TYPE::calcGradN( q, stack.xLocal, dNdX );
     FE_TYPE::calcN( q, N );
 
     real64 qp_damage = 0.0;
-    FE_TYPE::value( N, stack.nodalDamageLocal, qp_damage );
+    finiteElement::feOps::value( N, stack.nodalDamageLocal, qp_damage );
 
     real64 qp_disp[3] = {0, 0, 0};
-    FE_TYPE::value( N, stack.u_local, qp_disp );
+    finiteElement::feOps::value( N, stack.u_local, qp_disp );
 
     real64 elemPresGradient[3] = {0, 0, 0};
     for( integer i=0; i<3; ++i )
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp
index 74dad430581..3b4ef52d2fe 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp
@@ -400,9 +400,8 @@ void SolidMechanicsLagrangianFEM::initializePostInitialConditionsPreSubGroups()
             typename FE_TYPE::StackVariables feStack;
             finiteElement.template setup< FE_TYPE >( k, meshData, feStack );
             localIndex const numSupportPoints =
-              finiteElement.template numSupportPoints< FE_TYPE >( feStack );
+              finiteElement.getNumSupportPoints( feStack );
 
-//#if ! defined( CALC_FEM_SHAPE_IN_KERNEL ) // we don't calculate detJ in this case
             for( localIndex q=0; q<numQuadraturePointsPerElem; ++q )
             {
               FE_TYPE::calcN( q, feStack, N );
@@ -412,7 +411,6 @@ void SolidMechanicsLagrangianFEM::initializePostInitialConditionsPreSubGroups()
                 mass[elemsToNodes[k][a]] += rho[k][q] * detJ[k][q] * N[a];
               }
             }
-//#endif
 
             bool isAttachedToGhostNode = false;
             for( localIndex a=0; a<elementSubRegion.numNodesPerElement(); ++a )
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsEFEMKernelsBase.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsEFEMKernelsBase.hpp
index eef5c878ab1..f05817185e7 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsEFEMKernelsBase.hpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsEFEMKernelsBase.hpp
@@ -247,7 +247,7 @@ class EFEMKernelsBase :
     localIndex const embSurfIndex = m_cellsToEmbeddedSurfaces[k][0];
 
     real64 dNdX[ numNodesPerElem ][ 3 ];
-    real64 const detJ = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.X, dNdX );
+    real64 const detJ = FE_TYPE::calcGradN( q, stack.X, dNdX );
 
     constexpr int nUdof = numNodesPerElem*3;
 
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitFiniteStrain_impl.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitFiniteStrain_impl.hpp
index 8454837b438..eee43c6a04a 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitFiniteStrain_impl.hpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitFiniteStrain_impl.hpp
@@ -67,9 +67,7 @@ void ExplicitFiniteStrain< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >::setup(
     localIndex const nodeIndex = m_elemsToNodes( k, a );
     for( int i=0; i<numDofPerTrialSupportPoint; ++i )
     {
-#if defined(CALC_FEM_SHAPE_IN_KERNEL)
       stack.xLocal[ a ][ i ] = m_X[ nodeIndex ][ i ];
-#endif
       stack.uLocal[ a ][ i ] = m_u[ nodeIndex ][ i ];
       stack.varLocal[ a ][ i ] = m_vel[ nodeIndex ][ i ];
     }
@@ -86,7 +84,7 @@ void ExplicitFiniteStrain< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >::quadrat
                                                                                                 StackVariables & stack ) const
 {
   real64 dNdX[ numNodesPerElem ][ 3 ];
-  real64 const detJ = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal, dNdX );
+  real64 const detJ = FE_TYPE::calcGradN( q, stack.xLocal, dNdX );
 
   real64 dUhatdX[3][3] = { {0} };
   real64 dUdX[3][3] = { {0} };
@@ -94,8 +92,8 @@ void ExplicitFiniteStrain< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >::quadrat
   real64 Ldt[3][3] = { {0} };
   real64 fInv[3][3] = { {0} };
 
-  FE_TYPE::gradient( dNdX, stack.varLocal, dUhatdX );
-  FE_TYPE::gradient( dNdX, stack.uLocal, dUdX );
+  finiteElement::feOps::gradient( dNdX, stack.varLocal, dUhatdX );
+  finiteElement::feOps::gradient( dNdX, stack.uLocal, dUdX );
 
   LvArray::tensorOps::scale< 3, 3 >( dUhatdX, m_dt );
 
@@ -127,7 +125,7 @@ void ExplicitFiniteStrain< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >::quadrat
   LvArray::tensorOps::Rij_eq_symAikBjk< 3 >( P, stress, fInv );
   LvArray::tensorOps::scale< 3, 3 >( P, -detJ * detF );
 
-  FE_TYPE::plusGradNajAij( dNdX, P, stack.fLocal );
+  finiteElement::feOps::plusGradNajAij( dNdX, P, stack.fLocal );
 }
 
 
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitSmallStrain.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitSmallStrain.hpp
index b4518ca4cbd..aa590299f28 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitSmallStrain.hpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitSmallStrain.hpp
@@ -131,13 +131,8 @@ class ExplicitSmallStrain : public finiteElement::KernelBase< SUBREGION_TYPE,
     /// C-array stack storage for element local primary variable values.
     real64 varLocal[ numNodesPerElem ][ numDofPerTestSupportPoint ];
 
-#if !defined(CALC_FEM_SHAPE_IN_KERNEL)
-    /// Dummy
-    int xLocal;
-#else
     /// C-array stack storage for element local the nodal positions.
     real64 xLocal[ numNodesPerElem ][ 3 ];
-#endif
   };
   //***************************************************************************
 
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitSmallStrain_impl.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitSmallStrain_impl.hpp
index 30910aac4b0..de1de659b3b 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitSmallStrain_impl.hpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitSmallStrain_impl.hpp
@@ -23,6 +23,7 @@
 //#include "ExplicitFiniteStrain.hpp"
 #include "ExplicitSmallStrain.hpp"
 
+#include "finiteElement/elementFormulations/FiniteElementOperators.hpp"
 
 namespace geos
 {
@@ -73,10 +74,7 @@ void ExplicitSmallStrain< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >::setup( l
     localIndex const nodeIndex = m_elemsToNodes( k, a );
     for( int i=0; i<numDofPerTrialSupportPoint; ++i )
     {
-#if defined(CALC_FEM_SHAPE_IN_KERNEL)
       stack.xLocal[ a ][ i ] = m_X[ nodeIndex ][ i ];
-#endif
-
 #if UPDATE_STRESS==2
       stack.varLocal[ a ][ i ] = m_vel[ nodeIndex ][ i ] * m_dt;
 #else
@@ -98,11 +96,11 @@ void ExplicitSmallStrain< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >::quadratu
 //#define USE_JACOBIAN
 #if !defined( USE_JACOBIAN )
   real64 dNdX[ numNodesPerElem ][ 3 ];
-  real64 const detJ = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal, dNdX );
+  real64 const detJ = FE_TYPE::calcGradN( q, stack.xLocal, dNdX );
   /// Macro to substitute in the shape function derivatives.
   real64 strain[6] = {0};
   //real64 timeIncrement = 0.0;
-  FE_TYPE::symmetricGradient( dNdX, stack.varLocal, strain );
+  finiteElement::feOps::symmetricGradient( dNdX, stack.varLocal, strain );
 
   real64 stressLocal[ 6 ] = {0};
 #if UPDATE_STRESS == 2
@@ -124,14 +122,14 @@ void ExplicitSmallStrain< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >::quadratu
 #endif
   }
 
-  FE_TYPE::plusGradNajAij( dNdX, stressLocal, stack.fLocal );
+  finiteElement::feOps::plusGradNajAij( dNdX, stressLocal, stack.fLocal );
 
 #else
   real64 invJ[3][3];
   real64 const detJ = FE_TYPE::inverseJacobianTransformation( q, stack.xLocal, invJ );
 
   real64 strain[6] = {0};
-  FE_TYPE::symmetricGradient( q, invJ, stack.varLocal, strain );
+  finiteElement::feOps::symmetricGradient( q, invJ, stack.varLocal, strain );
 
   real64 stressLocal[ 6 ] = {0};
 #if UPDATE_STRESS == 2
@@ -153,7 +151,7 @@ void ExplicitSmallStrain< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >::quadratu
 #endif
   }
 
-  FE_TYPE::plusGradNajAij( q, invJ, stressLocal, stack.fLocal );
+  finiteElement::feOps::plusGradNajAij( q, invJ, stressLocal, stack.fLocal );
 #endif
 }
 
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/FixedStressThermoPoromechanics.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/FixedStressThermoPoromechanics.hpp
index 4b67faf912a..a2a751c7ae2 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/FixedStressThermoPoromechanics.hpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/FixedStressThermoPoromechanics.hpp
@@ -117,13 +117,8 @@ class FixedStressThermoPoromechanics :
                                        constitutiveStiffness()
     {}
 
-#if !defined(CALC_FEM_SHAPE_IN_KERNEL)
-    /// Dummy
-    int xLocal;
-#else
     /// C-array stack storage for element local the nodal positions.
     real64 xLocal[ numNodesPerElem ][ 3 ];
-#endif
 
     /// Stack storage for the element local nodal displacement
     real64 u_local[numNodesPerElem][numDofPerTrialSupportPoint];
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/FixedStressThermoPoromechanics_impl.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/FixedStressThermoPoromechanics_impl.hpp
index 36d70558513..0ecea9672db 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/FixedStressThermoPoromechanics_impl.hpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/FixedStressThermoPoromechanics_impl.hpp
@@ -21,6 +21,7 @@
 #define GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_FIXEDSTRESSTHERMOPOROMECHANICS_IMPL_HPP_
 
 #include "FixedStressThermoPoromechanics.hpp"
+#include "finiteElement/elementFormulations/FiniteElementOperators.hpp"
 #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp"
 #include "physicsSolvers/multiphysics/PoromechanicsFields.hpp"
 #include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp"
@@ -84,7 +85,7 @@ setup( localIndex const k,
 {
   m_finiteElementSpace.template setup< FE_TYPE >( k, m_meshData, stack.feStack );
   localIndex const numSupportPoints =
-    m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
+    m_finiteElementSpace.getNumSupportPoints( stack.feStack );
   stack.numRows = 3 * numSupportPoints;
   stack.numCols = stack.numRows;
 
@@ -94,9 +95,7 @@ setup( localIndex const k,
 
     for( int i = 0; i < 3; ++i )
     {
-#if defined(CALC_FEM_SHAPE_IN_KERNEL)
       stack.xLocal[ a ][ i ] = m_X[ localNodeIndex ][ i ];
-#endif
       stack.u_local[ a ][i] = m_disp[ localNodeIndex ][i];
       stack.uhat_local[ a ][i] = m_uhat[ localNodeIndex ][i];
       stack.localRowDofIndex[a*3+i] = m_dofNumber[localNodeIndex]+i;
@@ -124,7 +123,7 @@ quadraturePointKernel( localIndex const k,
                        StackVariables & stack ) const
 {
   real64 dNdX[ numNodesPerElem ][ 3 ];
-  real64 const detJxW = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal,
+  real64 const detJxW = FE_TYPE::calcGradN( q, stack.xLocal,
                                                                            stack.feStack, dNdX );
 
   real64 strainInc[6] = {0};
@@ -132,7 +131,7 @@ quadraturePointKernel( localIndex const k,
 
   typename CONSTITUTIVE_TYPE::KernelWrapper::DiscretizationOps stiffness;
 
-  FE_TYPE::symmetricGradient( dNdX, stack.uhat_local, strainInc );
+  finiteElement::feOps::symmetricGradient( dNdX, stack.uhat_local, strainInc );
 
   // Evaluate total stress and its derivatives
   // TODO: allow for a customization of the kernel to pass the average pressure to the small strain update (to account for cap pressure
@@ -160,7 +159,7 @@ quadraturePointKernel( localIndex const k,
 
   real64 N[numNodesPerElem];
   FE_TYPE::calcN( q, stack.feStack, N );
-  FE_TYPE::plusGradNajAijPlusNaFi( dNdX,
+  finiteElement::feOps::plusGradNajAijPlusNaFi( dNdX,
                                    totalStress,
                                    N,
                                    gravityForce,
@@ -190,7 +189,7 @@ complete( localIndex const k,
   // TODO: Does this work if BTDB is non-symmetric?
   CONSTITUTIVE_TYPE::KernelWrapper::DiscretizationOps::template fillLowerBTDB< numNodesPerElem >( stack.localJacobian );
   localIndex const numSupportPoints =
-    m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
+    m_finiteElementSpace.getNumSupportPoints( stack.feStack );
   for( int localNode = 0; localNode < numSupportPoints; ++localNode )
   {
     for( int dim = 0; dim < numDofPerTestSupportPoint; ++dim )
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainNewmark.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainNewmark.hpp
index 5b97f0022da..2399ad03ead 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainNewmark.hpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainNewmark.hpp
@@ -120,7 +120,7 @@ class ImplicitSmallStrainNewmark : public ImplicitSmallStrainQuasiStatic< SUBREG
     GEOS_HOST_DEVICE
     StackVariables():
       Base::StackVariables(),
-            dRdU_InertiaMassDamping{ {0.0} },
+      dRdU_InertiaMassDamping{ {0.0} },
       vtilde_local(),
       uhattilde_local()
     {}
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic.hpp
index 0e9bcf87b18..b5f1d5a05d6 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic.hpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic.hpp
@@ -128,13 +128,8 @@ class ImplicitSmallStrainQuasiStatic :
                                        constitutiveStiffness()
     {}
 
-#if !defined(CALC_FEM_SHAPE_IN_KERNEL)
-    /// Dummy
-    int xLocal;
-#else
     /// C-array stack storage for element local the nodal positions.
     real64 xLocal[ numNodesPerElem ][ 3 ];
-#endif
 
     /// Stack storage for the element local nodal displacement
     real64 u_local[numNodesPerElem][numDofPerTrialSupportPoint];
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic_impl.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic_impl.hpp
index 12753ab393f..56ee2ec9940 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic_impl.hpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic_impl.hpp
@@ -77,7 +77,7 @@ setup( localIndex const k,
 {
   m_finiteElementSpace.template setup< FE_TYPE >( k, m_meshData, stack.feStack );
 
-  localIndex const numSupportPoints = m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
+  localIndex const numSupportPoints = m_finiteElementSpace.getNumSupportPoints( stack.feStack );
 
   stack.numRows =  3 * numSupportPoints;
   stack.numCols = stack.numRows;
@@ -90,9 +90,7 @@ setup( localIndex const k,
     // #pragma unroll
     for( int i = 0; i < numDofPerTestSupportPoint; ++i )
     {
-#if defined(CALC_FEM_SHAPE_IN_KERNEL)
       stack.xLocal[ a ][ i ] = m_X[ localNodeIndex ][ i ];
-#endif
       stack.u_local[ a ][i] = m_disp[ localNodeIndex ][i];
       stack.uhat_local[ a ][i] = m_uhat[ localNodeIndex ][i];
       stack.localRowDofIndex[a*3+i] = m_dofNumber[localNodeIndex]+i;
@@ -120,14 +118,14 @@ void ImplicitSmallStrainQuasiStatic< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE
                                                                                                           STRESS_MODIFIER && stressModifier ) const
 {
   real64 dNdX[ numNodesPerElem ][ 3 ];
-  real64 const detJxW = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal, stack.feStack, dNdX );
+  real64 const detJxW = FE_TYPE::calcGradN( q, stack.xLocal, stack.feStack, dNdX );
 
   real64 strainInc[6] = {0};
   real64 stress[6] = {0};
 
   typename CONSTITUTIVE_TYPE::KernelWrapper::DiscretizationOps stiffness;
 
-  FE_TYPE::symmetricGradient( dNdX, stack.uhat_local, strainInc );
+  finiteElement::feOps::symmetricGradient( dNdX, stack.uhat_local, strainInc );
 
   m_constitutiveUpdate.smallStrainUpdate( k, q, m_dt, strainInc, stress, stiffness );
 
@@ -144,7 +142,7 @@ void ImplicitSmallStrainQuasiStatic< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE
 
   real64 N[numNodesPerElem];
   FE_TYPE::calcN( q, stack.feStack, N );
-  FE_TYPE::plusGradNajAijPlusNaFi( dNdX,
+  finiteElement::feOps::plusGradNajAijPlusNaFi( dNdX,
                                    stress,
                                    N,
                                    gravityForce,
@@ -177,7 +175,7 @@ real64 ImplicitSmallStrainQuasiStatic< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYP
   // TODO: Does this work if BTDB is non-symmetric?
   CONSTITUTIVE_TYPE::KernelWrapper::DiscretizationOps::template fillLowerBTDB< numNodesPerElem >( stack.localJacobian );
 #endif
-  localIndex const numSupportPoints = m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
+  localIndex const numSupportPoints = m_finiteElementSpace.getNumSupportPoints( stack.feStack );
 
   // #pragma unroll
   for( int localNode = 0; localNode < numSupportPoints; ++localNode )
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp
index faf51f89351..79e86129f47 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp
@@ -23,6 +23,7 @@
 #include "common/DataTypes.hpp"
 #include "common/GEOS_RAJA_Interface.hpp"
 #include "finiteElement/FiniteElementDispatch.hpp"
+#include "finiteElement/elementFormulations/FiniteElementOperators.hpp"
 #include "constitutive/ConstitutivePassThru.hpp"
 #include "mesh/CellElementSubRegion.hpp"
 #include "mesh/utilities/AverageOverQuadraturePointsKernel.hpp"
@@ -133,11 +134,11 @@ class AverageStrainOverQuadraturePoints :
     //real64 const weight = FE_TYPE::transformedQuadratureWeight( q, stack.xLocal, stack.feStack ) / m_elementVolume[k];
 
     real64 dNdX[ FE_TYPE::maxSupportPoints ][3];
-    real64 const detJxW = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal, stack.feStack, dNdX );
+    real64 const detJxW = FE_TYPE::calcGradN( q, stack.xLocal, stack.feStack, dNdX );
     real64 strain[6] = {0.0};
     real64 strainInc[6] = {0.0};
-    FE_TYPE::symmetricGradient( dNdX, stack.uLocal, strain );
-    FE_TYPE::symmetricGradient( dNdX, stack.uHatLocal, strainInc );
+    finiteElement::feOps::symmetricGradient( dNdX, stack.uLocal, strain );
+    finiteElement::feOps::symmetricGradient( dNdX, stack.uHatLocal, strainInc );
 
     real64 elasticStrainInc[6] = {0.0};
     m_solidUpdate.getElasticStrainInc( k, q, elasticStrainInc );
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticPMLSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticPMLSEMKernel.hpp
index 8f1b96077ba..322046c293c 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticPMLSEMKernel.hpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticPMLSEMKernel.hpp
@@ -20,6 +20,8 @@
 #ifndef GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICPMLSEMKERNEL_HPP_
 #define GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICPMLSEMKERNEL_HPP_
 
+#include "finiteElement/elementFormulations/FiniteElementOperators.hpp"
+
 namespace geos
 {
 
@@ -190,15 +192,15 @@ struct AcousticPMLSEM
         {
 
           /// compute the shape functions gradients
-          real32 const detJ = m_finiteElement.template getGradN< FE_TYPE >( k, i, xLocal, gradN );
+          real32 const detJ = FE_TYPE::calcGradN( i, xLocal, gradN );
           GEOS_UNUSED_VAR ( detJ );
 
           /// compute the gradient of the pressure and the PML auxiliary variables at the node
-          m_finiteElement.template gradient< numNodesPerElem, GRADIENT_TYPE >( gradN, pressure, pressureGrad );
-          m_finiteElement.template gradient< numNodesPerElem, GRADIENT_TYPE >( gradN, auxU, auxUGrad );
+          finiteElement::feOps::gradient< numNodesPerElem, GRADIENT_TYPE >( gradN, pressure, pressureGrad );
+          finiteElement::feOps::gradient< numNodesPerElem, GRADIENT_TYPE >( gradN, auxU, auxUGrad );
           for( int j=0; j<3; ++j )
           {
-            m_finiteElement.template gradient< numNodesPerElem, GRADIENT_TYPE >( gradN, auxV[j], auxVGrad[j] );
+            finiteElement::feOps::gradient< numNodesPerElem, GRADIENT_TYPE >( gradN, auxV[j], auxVGrad[j] );
           }
 
           /// compute the PML damping profile
diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd
index aca2fe256fa..3f4ccf0728d 100644
--- a/src/coreComponents/schema/schema.xsd
+++ b/src/coreComponents/schema/schema.xsd
@@ -3506,7 +3506,7 @@ Information output from lower logLevels is added with the desired log level
 		<!--maxNumResolves => Value to indicate how many resolves may be executed to perform surface generation after the execution of flow and mechanics solver. -->
 		<xsd:attribute name="maxNumResolves" type="integer" default="10" />
 		<!--newFractureInitializationType => Type of new fracture element initialization. Can be Pressure or Displacement. -->
-		<xsd:attribute name="newFractureInitializationType" type="geos_HydrofractureSolver_lt_geos_SinglePhasePoromechanics_lt_geos_SinglePhaseBase_cm_-geos_SolidMechanicsLagrangianFEM_gt_-_gt__InitializationType" default="Pressure" />
+		<xsd:attribute name="newFractureInitializationType" type="geos_HydrofractureSolver_lt_geos_SinglePhasePoromechanics_lt_geos_SinglePhaseBase_cm_-geos_SolidMechanicsLagrangianFEM_gt__gt__InitializationType" default="Pressure" />
 		<!--solidSolverName => Name of the solid solver used by the coupled solver-->
 		<xsd:attribute name="solidSolverName" type="groupNameRef" use="required" />
 		<!--stabilizationMultiplier => Constant multiplier of stabilization strength-->
@@ -3529,7 +3529,7 @@ Local- Add jump stabilization on interior of macro elements-->
 		<!--name => A name is required for any non-unique nodes-->
 		<xsd:attribute name="name" type="groupName" use="required" />
 	</xsd:complexType>
-	<xsd:simpleType name="geos_HydrofractureSolver_lt_geos_SinglePhasePoromechanics_lt_geos_SinglePhaseBase_cm_-geos_SolidMechanicsLagrangianFEM_gt_-_gt__InitializationType">
+	<xsd:simpleType name="geos_HydrofractureSolver_lt_geos_SinglePhasePoromechanics_lt_geos_SinglePhaseBase_cm_-geos_SolidMechanicsLagrangianFEM_gt__gt__InitializationType">
 		<xsd:restriction base="xsd:string">
 			<xsd:pattern value=".*[\[\]`$].*|Pressure|Displacement" />
 		</xsd:restriction>
diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other
index 90bd4aacb07..d37384973b0 100644
--- a/src/coreComponents/schema/schema.xsd.other
+++ b/src/coreComponents/schema/schema.xsd.other
@@ -486,7 +486,7 @@
 		<!--fractureStencil => (no description available)-->
 		<xsd:attribute name="fractureStencil" type="geos_SurfaceElementStencil" />
 		<!--targetRegions => List of regions to build the stencil for-->
-		<xsd:attribute name="targetRegions" type="geos_mapBase_lt_string_cm_-geos_internal_StdVectorWrapper_lt_string_cm_-std_allocator_lt_string-_gt__cm_-true_gt__cm_-std_integral_constant_lt_bool_cm_-true_gt_-_gt_" />
+		<xsd:attribute name="targetRegions" type="geos_mapBase_lt_std___1_basic_string_lt_char_cm_-std___1_char_traits_lt_char_gt__cm_-std___1_allocator_lt_char_gt__gt__cm_-std___1_vector_lt_std___1_basic_string_lt_char_cm_-std___1_char_traits_lt_char_gt__cm_-std___1_allocator_lt_char_gt__gt__cm_-std___1_allocator_lt_std___1_basic_string_lt_char_cm_-std___1_char_traits_lt_char_gt__cm_-std___1_allocator_lt_char_gt__gt__gt__gt__cm_-std___1_integral_constant_lt_bool_cm_-true_gt__gt_" />
 	</xsd:complexType>
 	<xsd:complexType name="OutputsType">
 		<xsd:choice minOccurs="0" maxOccurs="unbounded">
@@ -1343,7 +1343,7 @@
 			<xsd:element name="MeshBodies" type="MeshBodiesType" />
 		</xsd:choice>
 		<!--Neighbors => (no description available)-->
-		<xsd:attribute name="Neighbors" type="geos_internal_StdVectorWrapper_lt_geos_NeighborCommunicator_cm_-std_allocator_lt_geos_NeighborCommunicator_gt__cm_-true_gt_" />
+		<xsd:attribute name="Neighbors" type="std___1_vector_lt_geos_NeighborCommunicator_cm_-std___1_allocator_lt_geos_NeighborCommunicator_gt__gt_" />
 		<!--partitionManager => (no description available)-->
 		<xsd:attribute name="partitionManager" type="geos_PartitionBase" />
 	</xsd:complexType>
@@ -3079,7 +3079,7 @@
 		<!--ghostRank => (no description available)-->
 		<xsd:attribute name="ghostRank" type="integer_array" />
 		<!--globalToLocalMap => (no description available)-->
-		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std_integral_constant_lt_bool_cm_-false_gt_-_gt_" />
+		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std___1_integral_constant_lt_bool_cm_-false_gt__gt_" />
 		<!--isExternal => (no description available)-->
 		<xsd:attribute name="isExternal" type="integer_array" />
 		<!--localMaxGlobalIndex => (no description available)-->
@@ -3107,7 +3107,7 @@
 		<!--ghostRank => (no description available)-->
 		<xsd:attribute name="ghostRank" type="integer_array" />
 		<!--globalToLocalMap => (no description available)-->
-		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std_integral_constant_lt_bool_cm_-false_gt_-_gt_" />
+		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std___1_integral_constant_lt_bool_cm_-false_gt__gt_" />
 		<!--isExternal => (no description available)-->
 		<xsd:attribute name="isExternal" type="integer_array" />
 		<!--localMaxGlobalIndex => (no description available)-->
@@ -3126,11 +3126,11 @@
 		<!--domainBoundaryIndicator => (no description available)-->
 		<xsd:attribute name="domainBoundaryIndicator" type="integer_array" />
 		<!--faceList => (no description available)-->
-		<xsd:attribute name="faceList" type="geos_InterObjectRelation_lt_LvArray_ArrayOfSets_lt_int_cm_-int_cm_-LvArray_ChaiBuffer_gt_-_gt_" />
+		<xsd:attribute name="faceList" type="geos_InterObjectRelation_lt_LvArray_ArrayOfSets_lt_int_cm_-int_cm_-LvArray_ChaiBuffer_gt__gt_" />
 		<!--ghostRank => (no description available)-->
 		<xsd:attribute name="ghostRank" type="integer_array" />
 		<!--globalToLocalMap => (no description available)-->
-		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std_integral_constant_lt_bool_cm_-false_gt_-_gt_" />
+		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std___1_integral_constant_lt_bool_cm_-false_gt__gt_" />
 		<!--isExternal => (no description available)-->
 		<xsd:attribute name="isExternal" type="integer_array" />
 		<!--localMaxGlobalIndex => (no description available)-->
@@ -3140,7 +3140,7 @@
 		<!--maxGlobalIndex => (no description available)-->
 		<xsd:attribute name="maxGlobalIndex" type="globalIndex" />
 		<!--nodeList => (no description available)-->
-		<xsd:attribute name="nodeList" type="geos_InterObjectRelation_lt_LvArray_Array_lt_int_cm_-2_cm_-camp_int_seq_lt_long_cm_-0l_cm_-1l_gt__cm_-int_cm_-LvArray_ChaiBuffer_gt_-_gt_" />
+		<xsd:attribute name="nodeList" type="geos_InterObjectRelation_lt_LvArray_Array_lt_int_cm_-2_cm_-camp_int_seq_lt_long_cm_-0l_cm_-1l_gt__cm_-int_cm_-LvArray_ChaiBuffer_gt__gt_" />
 	</xsd:complexType>
 	<xsd:complexType name="embeddedSurfacesEdgeManagerType">
 		<xsd:choice minOccurs="0" maxOccurs="unbounded">
@@ -3150,11 +3150,11 @@
 		<!--domainBoundaryIndicator => (no description available)-->
 		<xsd:attribute name="domainBoundaryIndicator" type="integer_array" />
 		<!--faceList => (no description available)-->
-		<xsd:attribute name="faceList" type="geos_InterObjectRelation_lt_LvArray_ArrayOfSets_lt_int_cm_-int_cm_-LvArray_ChaiBuffer_gt_-_gt_" />
+		<xsd:attribute name="faceList" type="geos_InterObjectRelation_lt_LvArray_ArrayOfSets_lt_int_cm_-int_cm_-LvArray_ChaiBuffer_gt__gt_" />
 		<!--ghostRank => (no description available)-->
 		<xsd:attribute name="ghostRank" type="integer_array" />
 		<!--globalToLocalMap => (no description available)-->
-		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std_integral_constant_lt_bool_cm_-false_gt_-_gt_" />
+		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std___1_integral_constant_lt_bool_cm_-false_gt__gt_" />
 		<!--isExternal => (no description available)-->
 		<xsd:attribute name="isExternal" type="integer_array" />
 		<!--localMaxGlobalIndex => (no description available)-->
@@ -3164,7 +3164,7 @@
 		<!--maxGlobalIndex => (no description available)-->
 		<xsd:attribute name="maxGlobalIndex" type="globalIndex" />
 		<!--nodeList => (no description available)-->
-		<xsd:attribute name="nodeList" type="geos_InterObjectRelation_lt_LvArray_Array_lt_int_cm_-2_cm_-camp_int_seq_lt_long_cm_-0l_cm_-1l_gt__cm_-int_cm_-LvArray_ChaiBuffer_gt_-_gt_" />
+		<xsd:attribute name="nodeList" type="geos_InterObjectRelation_lt_LvArray_Array_lt_int_cm_-2_cm_-camp_int_seq_lt_long_cm_-0l_cm_-1l_gt__cm_-int_cm_-LvArray_ChaiBuffer_gt__gt_" />
 	</xsd:complexType>
 	<xsd:complexType name="embeddedSurfacesNodeManagerType">
 		<xsd:choice minOccurs="0" maxOccurs="unbounded">
@@ -3174,7 +3174,7 @@
 		<!--domainBoundaryIndicator => (no description available)-->
 		<xsd:attribute name="domainBoundaryIndicator" type="integer_array" />
 		<!--edgeList => (no description available)-->
-		<xsd:attribute name="edgeList" type="geos_InterObjectRelation_lt_LvArray_ArrayOfSets_lt_int_cm_-int_cm_-LvArray_ChaiBuffer_gt_-_gt_" />
+		<xsd:attribute name="edgeList" type="geos_InterObjectRelation_lt_LvArray_ArrayOfSets_lt_int_cm_-int_cm_-LvArray_ChaiBuffer_gt__gt_" />
 		<!--elemList => (no description available)-->
 		<xsd:attribute name="elemList" type="LvArray_ArrayOfArrays_lt_int_cm_-int_cm_-LvArray_ChaiBuffer_gt_" />
 		<!--elemRegionList => (no description available)-->
@@ -3184,7 +3184,7 @@
 		<!--ghostRank => (no description available)-->
 		<xsd:attribute name="ghostRank" type="integer_array" />
 		<!--globalToLocalMap => (no description available)-->
-		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std_integral_constant_lt_bool_cm_-false_gt_-_gt_" />
+		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std___1_integral_constant_lt_bool_cm_-false_gt__gt_" />
 		<!--isExternal => (no description available)-->
 		<xsd:attribute name="isExternal" type="integer_array" />
 		<!--localMaxGlobalIndex => (no description available)-->
@@ -3208,7 +3208,7 @@
 		<!--domainBoundaryIndicator => (no description available)-->
 		<xsd:attribute name="domainBoundaryIndicator" type="integer_array" />
 		<!--edgeList => (no description available)-->
-		<xsd:attribute name="edgeList" type="geos_InterObjectRelation_lt_LvArray_ArrayOfArrays_lt_int_cm_-int_cm_-LvArray_ChaiBuffer_gt_-_gt_" />
+		<xsd:attribute name="edgeList" type="geos_InterObjectRelation_lt_LvArray_ArrayOfArrays_lt_int_cm_-int_cm_-LvArray_ChaiBuffer_gt__gt_" />
 		<!--elemList => (no description available)-->
 		<xsd:attribute name="elemList" type="integer_array2d" />
 		<!--elemRegionList => (no description available)-->
@@ -3226,7 +3226,7 @@
 		<!--ghostRank => (no description available)-->
 		<xsd:attribute name="ghostRank" type="integer_array" />
 		<!--globalToLocalMap => (no description available)-->
-		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std_integral_constant_lt_bool_cm_-false_gt_-_gt_" />
+		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std___1_integral_constant_lt_bool_cm_-false_gt__gt_" />
 		<!--isExternal => (no description available)-->
 		<xsd:attribute name="isExternal" type="integer_array" />
 		<!--localMaxGlobalIndex => (no description available)-->
@@ -3238,7 +3238,7 @@
 		<!--mimGravityCoefficient => Mimetic gravity coefficient => CompositionalMultiphaseHybridFVM-->
 		<xsd:attribute name="mimGravityCoefficient" type="real64_array" />
 		<!--nodeList => (no description available)-->
-		<xsd:attribute name="nodeList" type="geos_InterObjectRelation_lt_LvArray_ArrayOfArrays_lt_int_cm_-int_cm_-LvArray_ChaiBuffer_gt_-_gt_" />
+		<xsd:attribute name="nodeList" type="geos_InterObjectRelation_lt_LvArray_ArrayOfArrays_lt_int_cm_-int_cm_-LvArray_ChaiBuffer_gt__gt_" />
 	</xsd:complexType>
 	<xsd:complexType name="nodeManagerType">
 		<xsd:choice minOccurs="0" maxOccurs="unbounded">
@@ -3250,7 +3250,7 @@
 		<!--domainBoundaryIndicator => (no description available)-->
 		<xsd:attribute name="domainBoundaryIndicator" type="integer_array" />
 		<!--edgeList => (no description available)-->
-		<xsd:attribute name="edgeList" type="geos_InterObjectRelation_lt_LvArray_ArrayOfSets_lt_int_cm_-int_cm_-LvArray_ChaiBuffer_gt_-_gt_" />
+		<xsd:attribute name="edgeList" type="geos_InterObjectRelation_lt_LvArray_ArrayOfSets_lt_int_cm_-int_cm_-LvArray_ChaiBuffer_gt__gt_" />
 		<!--elemList => (no description available)-->
 		<xsd:attribute name="elemList" type="LvArray_ArrayOfArrays_lt_int_cm_-int_cm_-LvArray_ChaiBuffer_gt_" />
 		<!--elemRegionList => (no description available)-->
@@ -3258,11 +3258,11 @@
 		<!--elemSubRegionList => (no description available)-->
 		<xsd:attribute name="elemSubRegionList" type="LvArray_ArrayOfArrays_lt_int_cm_-int_cm_-LvArray_ChaiBuffer_gt_" />
 		<!--faceList => (no description available)-->
-		<xsd:attribute name="faceList" type="geos_InterObjectRelation_lt_LvArray_ArrayOfSets_lt_int_cm_-int_cm_-LvArray_ChaiBuffer_gt_-_gt_" />
+		<xsd:attribute name="faceList" type="geos_InterObjectRelation_lt_LvArray_ArrayOfSets_lt_int_cm_-int_cm_-LvArray_ChaiBuffer_gt__gt_" />
 		<!--ghostRank => (no description available)-->
 		<xsd:attribute name="ghostRank" type="integer_array" />
 		<!--globalToLocalMap => (no description available)-->
-		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std_integral_constant_lt_bool_cm_-false_gt_-_gt_" />
+		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std___1_integral_constant_lt_bool_cm_-false_gt__gt_" />
 		<!--isExternal => (no description available)-->
 		<xsd:attribute name="isExternal" type="integer_array" />
 		<!--localMaxGlobalIndex => (no description available)-->
@@ -3285,7 +3285,7 @@
 		<!--ghostRank => (no description available)-->
 		<xsd:attribute name="ghostRank" type="integer_array" />
 		<!--globalToLocalMap => (no description available)-->
-		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std_integral_constant_lt_bool_cm_-false_gt_-_gt_" />
+		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std___1_integral_constant_lt_bool_cm_-false_gt__gt_" />
 		<!--isExternal => (no description available)-->
 		<xsd:attribute name="isExternal" type="integer_array" />
 		<!--localMaxGlobalIndex => (no description available)-->
@@ -3311,7 +3311,7 @@
 		<!--ghostRank => (no description available)-->
 		<xsd:attribute name="ghostRank" type="integer_array" />
 		<!--globalToLocalMap => (no description available)-->
-		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std_integral_constant_lt_bool_cm_-false_gt_-_gt_" />
+		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std___1_integral_constant_lt_bool_cm_-false_gt__gt_" />
 		<!--isExternal => (no description available)-->
 		<xsd:attribute name="isExternal" type="integer_array" />
 		<!--localMaxGlobalIndex => (no description available)-->
@@ -3332,7 +3332,7 @@
 		<!--ghostRank => (no description available)-->
 		<xsd:attribute name="ghostRank" type="integer_array" />
 		<!--globalToLocalMap => (no description available)-->
-		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std_integral_constant_lt_bool_cm_-false_gt_-_gt_" />
+		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std___1_integral_constant_lt_bool_cm_-false_gt__gt_" />
 		<!--isExternal => (no description available)-->
 		<xsd:attribute name="isExternal" type="integer_array" />
 		<!--localMaxGlobalIndex => (no description available)-->
@@ -3362,7 +3362,7 @@
 		<!--ghostRank => (no description available)-->
 		<xsd:attribute name="ghostRank" type="integer_array" />
 		<!--globalToLocalMap => (no description available)-->
-		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std_integral_constant_lt_bool_cm_-false_gt_-_gt_" />
+		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std___1_integral_constant_lt_bool_cm_-false_gt__gt_" />
 		<!--isExternal => (no description available)-->
 		<xsd:attribute name="isExternal" type="integer_array" />
 		<!--localMaxGlobalIndex => (no description available)-->
@@ -3376,7 +3376,7 @@
 		<!--nextWellElementIndexGlobal => (no description available)-->
 		<xsd:attribute name="nextWellElementIndexGlobal" type="integer_array" />
 		<!--nodeList => (no description available)-->
-		<xsd:attribute name="nodeList" type="geos_InterObjectRelation_lt_LvArray_Array_lt_int_cm_-2_cm_-camp_int_seq_lt_long_cm_-0l_cm_-1l_gt__cm_-int_cm_-LvArray_ChaiBuffer_gt_-_gt_" />
+		<xsd:attribute name="nodeList" type="geos_InterObjectRelation_lt_LvArray_Array_lt_int_cm_-2_cm_-camp_int_seq_lt_long_cm_-0l_cm_-1l_gt__cm_-int_cm_-LvArray_ChaiBuffer_gt__gt_" />
 		<!--numEdgesPerElement => (no description available)-->
 		<xsd:attribute name="numEdgesPerElement" type="integer" />
 		<!--numFacesPerElement => (no description available)-->
@@ -3403,7 +3403,7 @@
 		<!--ghostRank => (no description available)-->
 		<xsd:attribute name="ghostRank" type="integer_array" />
 		<!--globalToLocalMap => (no description available)-->
-		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std_integral_constant_lt_bool_cm_-false_gt_-_gt_" />
+		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std___1_integral_constant_lt_bool_cm_-false_gt__gt_" />
 		<!--isExternal => (no description available)-->
 		<xsd:attribute name="isExternal" type="integer_array" />
 		<!--localMaxGlobalIndex => (no description available)-->
@@ -3442,7 +3442,7 @@
 		<!--ghostRank => (no description available)-->
 		<xsd:attribute name="ghostRank" type="integer_array" />
 		<!--globalToLocalMap => (no description available)-->
-		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std_integral_constant_lt_bool_cm_-false_gt_-_gt_" />
+		<xsd:attribute name="globalToLocalMap" type="geos_mapBase_lt_long-long_cm_-int_cm_-std___1_integral_constant_lt_bool_cm_-false_gt__gt_" />
 		<!--isExternal => (no description available)-->
 		<xsd:attribute name="isExternal" type="integer_array" />
 		<!--localMaxGlobalIndex => (no description available)-->
diff --git a/src/coreComponents/unitTests/virtualElementTests/testConformingVirtualElementOrder1.cpp b/src/coreComponents/unitTests/virtualElementTests/testConformingVirtualElementOrder1.cpp
index 266ce88d436..cf1c9a23344 100644
--- a/src/coreComponents/unitTests/virtualElementTests/testConformingVirtualElementOrder1.cpp
+++ b/src/coreComponents/unitTests/virtualElementTests/testConformingVirtualElementOrder1.cpp
@@ -35,8 +35,7 @@ constexpr real64 relTol = geos::testing::DEFAULT_REL_TOL*10;
 
 template< typename VEM >
 GEOS_HOST_DEVICE
-static void checkIntegralMeanConsistency( FiniteElementBase const & feBase,
-                                          typename VEM::StackVariables const & stack,
+static void checkIntegralMeanConsistency( typename VEM::StackVariables const & stack,
                                           real64 & sumBasisFunctions )
 {
   static constexpr localIndex
@@ -44,8 +43,7 @@ static void checkIntegralMeanConsistency( FiniteElementBase const & feBase,
   real64 basisFunctionsIntegralMean[maxSupportPoints];
   VEM::calcN( 0, stack, basisFunctionsIntegralMean );
   sumBasisFunctions = 0;
-  for( localIndex iBasisFun = 0;
-       iBasisFun < feBase.template numSupportPoints< VEM >( stack ); ++iBasisFun )
+  for( localIndex iBasisFun = 0; iBasisFun < VEM::getNumSupportPoints( stack ); ++iBasisFun )
   {
     sumBasisFunctions += basisFunctionsIntegralMean[iBasisFun];
   }
@@ -54,8 +52,7 @@ static void checkIntegralMeanConsistency( FiniteElementBase const & feBase,
 template< typename VEM >
 GEOS_HOST_DEVICE
 static void
-checkIntegralMeanDerivativesConsistency( FiniteElementBase const & feBase,
-                                         typename VEM::StackVariables const & stack,
+checkIntegralMeanDerivativesConsistency( typename VEM::StackVariables const & stack,
                                          real64 & sumXDerivatives,
                                          real64 & sumYDerivatives,
                                          real64 & sumZDerivatives )
@@ -63,15 +60,12 @@ checkIntegralMeanDerivativesConsistency( FiniteElementBase const & feBase,
   static constexpr localIndex
     maxSupportPoints = VEM::maxSupportPoints;
   real64 const dummy[VEM::numNodes][3] { { 0.0 } };
-  localIndex const k = 0;
   for( localIndex q = 0; q < VEM::numQuadraturePoints; ++q )
   {
     real64 basisDerivativesIntegralMean[maxSupportPoints][3]{};
-    feBase.template getGradN< VEM >( k, q, dummy, stack, basisDerivativesIntegralMean );
+    VEM::calcGradN( q, dummy, stack, basisDerivativesIntegralMean );
     sumXDerivatives = 0; sumYDerivatives = 0; sumZDerivatives = 0;
-    for( localIndex iBasisFun = 0;
-         iBasisFun < feBase.template numSupportPoints< VEM >( stack );
-         ++iBasisFun )
+    for( localIndex iBasisFun = 0; iBasisFun < VEM::getNumSupportPoints( stack ); ++iBasisFun )
     {
       sumXDerivatives += basisDerivativesIntegralMean[iBasisFun][0];
       sumYDerivatives += basisDerivativesIntegralMean[iBasisFun][1];
@@ -221,9 +215,9 @@ static void testCellsInMeshLevel( MeshLevel const & mesh )
     VEM virtualElement;
     virtualElement.template setup< VEM >( cellIndex, meshData, stack );
 
-    checkIntegralMeanConsistency< VEM >( virtualElement, stack,
+    checkIntegralMeanConsistency< VEM >( stack, 
                                          sumBasisFunctionsView( cellIndex ) );
-    checkIntegralMeanDerivativesConsistency< VEM >( virtualElement, stack,
+    checkIntegralMeanDerivativesConsistency< VEM >( stack,
                                                     sumXDerivativesView( cellIndex ),
                                                     sumYDerivativesView( cellIndex ),
                                                     sumZDerivativesView( cellIndex )
diff --git a/src/docs/doxygen/Doxyfile.in b/src/docs/doxygen/Doxyfile.in
index f3869f78f70..70ff31e11d8 100644
--- a/src/docs/doxygen/Doxyfile.in
+++ b/src/docs/doxygen/Doxyfile.in
@@ -2034,7 +2034,6 @@ EXPAND_AS_DEFINED      = LVARRAY_HOST_DEVICE \
                          DISABLE_HD_WARNING \
                          DECLARE_FIELD \
                          GEOS_FORCE_INLINE \
-                         USING_FINITEELEMENTBASE \ 
                          MAKE_PRINT_MATVIEW
 
 # If the SKIP_FUNCTION_MACROS tag is set to YES then doxygen's preprocessor will