move stuff into BasisProduct and InterpolatedShape classes from Paren…

…tElement

src/discretizations/CMakeLists.txt

@@ -27,4 +27,4 @@ install( TARGETS discretizations
          RUNTIME DESTINATION lib ) 
 
 
-add_subdirectory( unitTests )
+#add_subdirectory( unitTests )

src/discretizations/finiteElementMethod/parentElements/ParentElement.hpp

@@ -19,173 +19,70 @@ namespace finiteElementMethod
  * @brief Defines a class that provides static functions to calculate quantities
  * required from the parent element in a finite element method.
  * @tparam REAL_TYPE The floating point type to use
- * @tparam BASE_GEOMETRY The cell type/geometry
+ * @tparam SHAPE The cell type/geometry
  * @tparam FUNCTIONAL_SPACE_TYPE The functional space type
  * @tparam BASIS_TYPE Pack of basis types to apply to each direction of the
  * parent element. There should be a basis defined for each direction.
  */
-template< typename REAL_TYPE, typename BASE_GEOMETRY, typename ... BASIS_TYPE >
+template< typename REAL_TYPE, typename SHAPE, typename ... BASIS_TYPE >
 class ParentElement
 {
 
 public:
 
   /// The type used to represent the cell/geometry
-  using BaseGeometry = BASE_GEOMETRY;
+  using ShapeType = SHAPE;
 //  using FunctionalSpaceType = FUNCTIONAL_SPACE_TYPE;
-//  using IndexType = typename BaseGeometry::IndexType;
+//  using IndexType = typename ShapeType::IndexType;
 
   /// Alias for the floating point type
   using RealType = REAL_TYPE;
 
   /// Alias for the type that represents a coordinate
-  using CoordType = typename BaseGeometry::CoordType;
+  using CoordType = typename ShapeType::CoordType;
 
   /// The number of dimensions on the ParentElement
   static inline constexpr int numDims = sizeof...(BASIS_TYPE);
 
-  /// The number of vertices on the ParentElement
-  static inline constexpr int numVertices = BaseGeometry::numVertices();
+  /// The number of degrees of freedom on the ParentElement in each 
+  /// dimension/basis.
+  static inline constexpr int numDofs[numDims] = {BASIS_TYPE::numDofs()...};
 
 
-  static_assert( numDims == BaseGeometry::numDims(), "numDims mismatch between cell and number of basis specified" );
+  static_assert( numDims == ShapeType::numDims(), "numDims mismatch between cell and number of basis specified" );
 
   /**
-   * @brief Calculates the value of the basis function at the specified parent
-   * coordinate.
-   * @tparam BASIS_FUNCTION_INDICES Pack of indices of the basis function to
-   * evaluate in each dimension.
-   * @param parentCoord The parent coordinate at which to evaluate the basis
-   * function.
-   * @return The value of the basis function at the specified parent coordinate.
-   *
-   * The equation for the value of a basis is:
-   * \f[
-   * \Phi_{i_0 i_1 ... i_{(numDims-1)}}(\boldsymbol{\xi}) =
-   * \prod_{k=0}^{(numDims-1)} \phi_{i_k}(\xi_k)\text{, where } \\
-   * i_j \text{is index of the basis function in the jth dimension, and
-   * ranges from [0...(order+1)]}
-   * \f]
-   *
+   * @copydoc functions::BasisProduct::value
    */
   template< int ... BASIS_FUNCTION_INDICES >
   SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE RealType
   value( CoordType const & parentCoord )
   {
     static_assert( sizeof...(BASIS_FUNCTION_INDICES) == numDims, "Wrong number of basis function indicies specified" );
-
-    return
-#if __cplusplus >= 202002L
-      // expand pack over number of dimensions
-      executeSequence< numDims >( [&]< int ... PRODUCT_TERM_INDEX > () constexpr
-    {
-      return ( BASIS_TYPE::template value< BASIS_FUNCTION_INDICES >( parentCoord[PRODUCT_TERM_INDEX] ) * ... );
-    } );
-#else
-      executeSequence< numDims >( [&] ( auto ... PRODUCT_TERM_INDEX ) constexpr
-    {
-      // fold expression to multiply the value of each BASIS_TYPE in each
-      // dimension. In other words the fold expands on BASIS_TYPE...,
-      // BASIS_FUNCTION_INDICES..., and PRODUCT_TERM_INDEX... together.
-      return ( BASIS_TYPE::template value< BASIS_FUNCTION_INDICES >( parentCoord[decltype(PRODUCT_TERM_INDEX)::value] ) * ... );
-    } );
-
-#endif
+    return ( BASIS_PRODUCT_TYPE::template value< BASIS_FUNCTION_INDICES... >( parentCoord ) );
   }
 
+
   /**
-   * @brief Calculates the gradient of the basis function at the specified
-   * parent coordinate.
-   * @tparam BASIS_FUNCTION_INDICES Pack of indices of the basis function to
-   * evaluate in each dimension.
-   * @param parentCoord The parent coordinate at which to evaluate the basis
-   * function gradient.
-   * @return The gradient of the basis function at the specified parent
-   * coordinate.
-   *
-   * The equation for the gradient of a basis is:
-   * \f[
-   * \frac{ d\left( \Phi_{i_0 i_1 ... i_{(numDims-1)}}(\boldsymbol{\xi})
-   * \right) }{ d \xi_j} =
-   * \nabla \phi_{i_j}(\xi_j) \prod_{ {k=0}\atop {k\neq j} }^{(numDims-1)}
-   * \phi_{i_k}(\xi_k)
-   * \f]
+   * @copydoc functions::BasisProduct::gradient
    */
   template< int ... BASIS_FUNCTION_INDICES >
   SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE CArray1d< RealType, numDims >
   gradient( CoordType const & parentCoord )
   {
     static_assert( sizeof...(BASIS_FUNCTION_INDICES) == numDims, "Wrong number of basis function indicies specified" );
-
-#if __cplusplus >= 202002L
-    return executeSequence< numDims >( [&]< int ... i >() constexpr -> CArray1d< RealType, numDims >
-    {
-      auto gradientComponent = [&] ( auto const iGrad,
-                                     auto const  ... PRODUCT_TERM_INDICES ) constexpr
-      {
-        // Ca
-        return ( gradientComponentHelper< BASIS_TYPE,
-                                          decltype(iGrad)::value,
-                                          BASIS_FUNCTION_INDICES,
-                                          PRODUCT_TERM_INDICES >( parentCoord ) * ... );
-      };
-
-      return { (executeSequence< numDims >( gradientComponent, std::integral_constant< int, i >{} ) )...  };
-    } );
-#else
-    // Expand over the dimensions.
-    return executeSequence< numDims >( [&] ( auto ... a ) constexpr -> CArray1d< RealType, numDims >
-    {
-      // define a lambda that calculates the gradient of the basis function in
-      // a single dimension/direction.
-      auto gradientComponent = [&] ( auto GRADIENT_COMPONENT, auto ... PRODUCT_TERM_INDICES ) constexpr
-      {
-        // fold expression calling gradientComponentHelper using expanding on
-        // BASIS_TYPE, BASIS_FUNCTION_INDICES, and PRODUCT_TERM_INDICES.
-        return ( gradientComponentHelper< BASIS_TYPE,
-                                          decltype(GRADIENT_COMPONENT)::value,
-                                          BASIS_FUNCTION_INDICES,
-                                          decltype(PRODUCT_TERM_INDICES)::value >( parentCoord ) * ... );
-      };
-
-      // execute the gradientComponent lambda on each direction, expand the
-      // pack on "i" corresponding to each direction of the gradient.
-      return { (executeSequence< numDims >( gradientComponent, a ) )...  };
-    } );
-#endif
+    return ( BASIS_PRODUCT_TYPE::template gradient< BASIS_FUNCTION_INDICES... >( parentCoord ) );
   }
 
+  // template< typename VAR_DATA, typename VALUE_TYPE >
+  // SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE VALUE_TYPE
+  // value( CoordType const & parentCoord,
+  //        VAR_DATA const & var )
+  // {
+  //   static_assert( sizeof...(BASIS_FUNCTION_INDICES) == numDims, "Wrong number of basis function indicies specified" );
 
-private:
-  /**
-   * @brief Helper function to return the gradient of a basis function, or the
-   * value of the basis function depending on whether or not the index of the
-   * basis function (@p BASIS_FUNCTION) matches the index of the direction
-   * component index (@p GRADIENT_COMPONENT). This filter is illustrated in the
-   * documentation for the gradient function.
-   * @tparam BASIS_FUNCTION The basis function type
-   * @tparam GRADIENT_COMPONENT The dimension component of the gradient.
-   * @tparam BASIS_FUNCTION_INDEX The index of the basis function that is being
-   * evaluated.
-   * @tparam COORD_INDEX The dimension component of the coordinate.
-   * @param parentCoord The parent coordinate at which to evaluate the basis
-   * function gradient.
-   * @return The gradient component of the basis function at the specified
-   * parent coordinate.
-   */
-  template< typename BASIS_FUNCTION, int GRADIENT_COMPONENT, int BASIS_FUNCTION_INDEX, int COORD_INDEX >
-  SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE RealType
-  gradientComponentHelper( CoordType const & parentCoord )
-  {
-    if constexpr ( GRADIENT_COMPONENT == COORD_INDEX )
-    {
-      return BASIS_FUNCTION::template gradient< BASIS_FUNCTION_INDEX >( parentCoord[COORD_INDEX] );
-    }
-    else
-    {
-      return ( BASIS_FUNCTION::template value< BASIS_FUNCTION_INDEX >( parentCoord[COORD_INDEX] ) );
-    }
-  }
+  //   return ( BASIS_PRODUCT_TYPE::template value< BASIS_FUNCTION_INDICES... >( parentCoord ) );
+  // }