feat: split FiniteElementBase into nonvirtual and virtual interfaces

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 )

src/cmake/thirdparty/SetupGeosxThirdParty.cmake

@@ -887,7 +887,7 @@ if( ${CMAKE_VERSION} VERSION_LESS "3.19" )
     set( PYTHON_AND_VERSION Python3 )
     set( PYTHON_OPTIONAL_COMPONENTS)
 else()
-    set( PYTHON_AND_VERSION Python3 3.6.0...<4 )
+    set( PYTHON_AND_VERSION Python3 3.6.0...<4.0.0 )
     set( PYTHON_OPTIONAL_COMPONENTS OPTIONAL_COMPONENTS Development NumPy)
 endif()
 if(ENABLE_PYGEOSX)

src/coreComponents/constitutive/fluid/twophaseimmisciblefluid/TwoPhaseImmiscibleFluid.hpp

@@ -253,7 +253,7 @@ void TwoPhaseImmiscibleFluid::KernelWrapper::
 {
   using Deriv = constitutive::multifluid::DerivativeOffset;
 
-  LvArray::forValuesInSlice( phaseDensity.derivs, []( real64 & val ) { val = 0.0; } );
+  // LvArray::forValuesInSlice( phaseDensity.derivs, []( real64 & val ) { val = 0.0; } );
 
   for( integer iph = 0; iph < 2; ++iph )
   {
@@ -274,7 +274,7 @@ void TwoPhaseImmiscibleFluid::KernelWrapper::
 {
   using Deriv = constitutive::multifluid::DerivativeOffset;
 
-  LvArray::forValuesInSlice( phaseViscosity.derivs, []( real64 & val ) { val = 0.0; } );
+  // LvArray::forValuesInSlice( phaseViscosity.derivs, []( real64 & val ) { val = 0.0; } );
 
   for( integer iph = 0; iph < 2; ++iph )
   {

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

src/coreComponents/finiteElement/FiniteElementDiscretization.hpp

@@ -27,7 +27,6 @@
 #include "LvArray/src/tensorOps.hpp"
 #include "FiniteElementDispatch.hpp"
 
-
 namespace geos
 {
 
@@ -42,8 +41,6 @@ string const detJ = "detJ";
 }
 }
 
-
-
 class FiniteElementDiscretization : public dataRepository::Group
 {
 public:
@@ -157,7 +154,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 ];
@@ -183,7 +180,6 @@ FiniteElementDiscretization::
 
 }
 
-
 } /* namespace geos */
 
 #endif /* GEOS_FINITEELEMENT_FINITEELEMENTDISCRETIZATION_HPP_ */

src/coreComponents/finiteElement/FiniteElementDispatch.hpp

@@ -159,7 +159,7 @@ struct FiniteElementDispatchHandler< FE_TYPE, FE_TYPES... >
   dispatch3D( FiniteElementBase const & input,
               LAMBDA && lambda )
   {
-    if( auto const * const ptr = dynamic_cast< FE_TYPE const * >(&input) )
+    if( auto const * const ptr = dynamic_cast< FE_TYPE const * >(&input)->getImpl() )
     {
       lambda( *ptr );
     }
@@ -174,7 +174,7 @@ struct FiniteElementDispatchHandler< FE_TYPE, FE_TYPES... >
   dispatch3D( FiniteElementBase & input,
               LAMBDA && lambda )
   {
-    if( auto * const ptr = dynamic_cast< FE_TYPE * >(&input) )
+    if( auto * const ptr = dynamic_cast< FE_TYPE * >(&input)->getImpl() )
     {
       lambda( *ptr );
     }
@@ -189,7 +189,7 @@ struct FiniteElementDispatchHandler< FE_TYPE, FE_TYPES... >
   dispatch2D( FiniteElementBase const & input,
               LAMBDA && lambda )
   {
-    if( auto const * const ptr = dynamic_cast< FE_TYPE const * >(&input) )
+    if( auto const * const ptr = dynamic_cast< FE_TYPE const * >(&input)->getImpl() )
     {
       lambda( *ptr );
     }
@@ -207,19 +207,19 @@ void
 dispatchlowOrder3D( FiniteElementBase const & input,
                     LAMBDA && lambda )
 {
-  if( auto const * const ptr1 = dynamic_cast< H1_Hexahedron_Lagrange1_GaussLegendre2 const * >(&input) )
+  if( auto const * const ptr1 = dynamic_cast< H1_Hexahedron_Lagrange1_GaussLegendre2 const * >(&input)->getImpl() )
   {
     lambda( *ptr1 );
   }
-  else if( auto const * const ptr2 = dynamic_cast< H1_Wedge_Lagrange1_Gauss6 const * >(&input) )
+  else if( auto const * const ptr2 = dynamic_cast< H1_Wedge_Lagrange1_Gauss6 const * >(&input)->getImpl() )
   {
     lambda( *ptr2 );
   }
-  else if( auto const * const ptr3 = dynamic_cast< H1_Tetrahedron_Lagrange1_Gauss1 const * >(&input) )
+  else if( auto const * const ptr3 = dynamic_cast< H1_Tetrahedron_Lagrange1_Gauss1 const * >(&input)->getImpl() )
   {
     lambda( *ptr3 );
   }
-  else if( auto const * const ptr4 = dynamic_cast< H1_Pyramid_Lagrange1_Gauss5 const * >(&input) )
+  else if( auto const * const ptr4 = dynamic_cast< H1_Pyramid_Lagrange1_Gauss5 const * >(&input)->getImpl() )
   {
     lambda( *ptr4 );
   }
@@ -233,6 +233,4 @@ dispatchlowOrder3D( FiniteElementBase const & input,
 
 } // namespace geos
 
-
-
 #endif /* GEOS_FINITEELEMENT_FINITEELEMENTDISPATCH_HPP_ */

src/coreComponents/finiteElement/elementFormulations/BB_Tetrahedron.hpp

@@ -30,6 +30,12 @@ namespace geos
 namespace finiteElement
 {
 
+namespace
+{
+template< int ORDER >
+constexpr int BB_Tetrahedron_NumNodes = ( ORDER + 1 ) * ( ORDER + 2 ) * ( ORDER + 3 ) / 6;
+}
+
 /**
  * This class contains the kernel accessible functions specific to the
  * Bernstein-Bézier (BB) modal any-order tetrahedron finite element with
@@ -42,33 +48,37 @@ namespace finiteElement
  * by 0<=l1,l2,l3,l3<=1, l1+l2+l3+l4=1
  */
 template< int ORDER >
-class BB_Tetrahedron final : public FiniteElementBase
+class BB_Tetrahedron_impl : public FiniteElementBase_impl< BB_Tetrahedron_NumNodes< ORDER >,
+                                                           4,
+                                                           BB_Tetrahedron_NumNodes< ORDER > >
 {
 public:
+  /// Convenience alias for base class type.
+  using Base = FiniteElementBase_impl< BB_Tetrahedron_NumNodes< ORDER >,
+                                       4,
+                                       BB_Tetrahedron_NumNodes< ORDER > >;
 
-  /// The order of the finite element.
-  static constexpr int order = ORDER;
+  /// struct to hold stack variables.
+  struct StackVariables {};
 
   /// The number of shape functions per element.
-  constexpr static localIndex numNodes = ( ORDER + 1 ) * ( ORDER + 2 ) * ( ORDER + 3 ) / 6;
+  using Base::numNodes;
 
-  /// The number of shape functions per face
-  constexpr static localIndex numNodesPerFace = ( ORDER + 1 ) * ( ORDER + 2 ) / 2;
+  /// The number of quadrature points per element.
+  using Base::numQuadraturePoints;
 
   /// The maximum number of support points per element.
-  constexpr static localIndex maxSupportPoints = numNodes;
+  using Base::maxSupportPoints;
 
-  /// The number of quadrature points per element.
-  constexpr static localIndex numQuadraturePoints = numNodes;
-
-  /** @cond Doxygen_Suppress */
-  USING_FINITEELEMENTBASE
-  /** @endcond Doxygen_Suppress */
+  /// The order of the finite element.
+  static constexpr int order = ORDER;
 
-  virtual ~BB_Tetrahedron() = default;
+  /// The number of shape functions per face
+  constexpr static localIndex numNodesPerFace = ( ORDER + 1 ) * ( ORDER + 2 ) / 2;
 
+  /// @copydoc FiniteElementBase_impl::getNumQuadraturePoints()
   GEOS_HOST_DEVICE
-  virtual localIndex getNumQuadraturePoints() const override
+  static localIndex getNumQuadraturePoints()
   {
     return numQuadraturePoints;
   }
@@ -86,16 +96,18 @@ class BB_Tetrahedron final : public FiniteElementBase
     return numQuadraturePoints;
   }
 
+  /// @copydoc FiniteElementBase_impl::getNumSupportPoints()
   GEOS_HOST_DEVICE
   GEOS_FORCE_INLINE
-  virtual localIndex getNumSupportPoints() const override
+  static localIndex getNumSupportPoints()
   {
     return numNodes;
   }
 
+  /// @copydoc FiniteElementBase_impl::getMaxSupportPoints()
   GEOS_HOST_DEVICE
   GEOS_FORCE_INLINE
-  virtual localIndex getMaxSupportPoints() const override
+  static localIndex getMaxSupportPoints()
   {
     return maxSupportPoints;
   }
@@ -122,7 +134,6 @@ class BB_Tetrahedron final : public FiniteElementBase
   GEOS_HOST_DEVICE
   GEOS_FORCE_INLINE
   static real64 jacobianDeterminant( real64 const (&X)[4][3] )
-
   {
     real64 m[3][3] = {};
     for( int i = 0; i < 3; i++ )
@@ -1320,21 +1331,81 @@ class BB_Tetrahedron final : public FiniteElementBase
 
 };
 
-/**
- *  Tetrahedron element with Bernstein-Bézier basis functions of order 1.
- */
+/// @copydoc BB_Tetrahedron_impl
+template< int ORDER >
+class BB_Tetrahedron final : public BB_Tetrahedron_impl< ORDER >, public FiniteElementBase
+{
+public:
+
+  /// The Implementation type
+  using ImplType = BB_Tetrahedron_impl< ORDER >;
+
+  /// The number of nodes per element.
+  using ImplType::numNodes;
+
+  /// The max number of support points per element.
+  using ImplType::maxSupportPoints;
+
+  /// The number of quadrature points per element.
+  using ImplType::numQuadraturePoints;
+
+  BB_Tetrahedron():
+    FiniteElementBase( numNodes,
+                       maxSupportPoints,
+                       numQuadraturePoints )
+  { }
+
+#ifdef __CUDACC__
+  #pragma diag_push
+  #pragma nv_diag_suppress 20012
+#endif
+  GEOS_HOST_DEVICE
+  virtual ~BB_Tetrahedron() override final = default;
+#ifdef __CUDACC__
+  #pragma diag_pop
+#endif
+
+  /**
+   * @brief Get the device-compatible implementation type.
+   * @return A pointer to the device-compatible implementation type.
+   */
+  ImplType * getImpl()
+  {
+    return static_cast< ImplType * >(this);
+  }
+
+  /**
+   * @brief Get the device-compatible implementation type.
+   * @return A pointer to the device-compatible implementation type.
+   */
+  ImplType const * getImpl() const
+  {
+    return static_cast< ImplType const * >(this);
+  }
+
+};
+
+///  Tetrahedron element with Bernstein-Bézier basis functions of order 1.
 using BB1_Tetrahedron = BB_Tetrahedron< 1 >;
-/**
- *  Tetrahedron element with Bernstein-Bézier basis functions of order 2.
- */
+///  Tetrahedron element with Bernstein-Bézier basis functions of order 2.
 using BB2_Tetrahedron = BB_Tetrahedron< 2 >;
-/**
- *  Tetrahedron element with Bernstein-Bézier basis functions of order 3.
- */
+///  Tetrahedron element with Bernstein-Bézier basis functions of order 3.
 using BB3_Tetrahedron = BB_Tetrahedron< 3 >;
+
 //using BB4_Tetrahedron = BB_Tetrahedron< 4 >;
 //using BB5_Tetrahedron = BB_Tetrahedron< 5 >;
 
+
+///  Tetrahedron element with Bernstein-Bézier basis functions of order 1.
+using BB1_Tetrahedron_impl = BB_Tetrahedron_impl< 1 >;
+///  Tetrahedron element with Bernstein-Bézier basis functions of order 2.
+using BB2_Tetrahedron_impl = BB_Tetrahedron_impl< 2 >;
+///  Tetrahedron element with Bernstein-Bézier basis functions of order 3.
+using BB3_Tetrahedron_impl = BB_Tetrahedron_impl< 3 >;
+
+// using BB2_Tetrahedron_impl = BB_Tetrahedron_impl< 4 >;
+// using BB3_Tetrahedron_impl = BB_Tetrahedron_impl< 5 >;
+
 }
 }