Add finite element space

.gitignore

@@ -4,6 +4,12 @@ install-*
 .vscode
 
 
+#vscode
+.vscode
+
+#mac directories
+.DS_Store
+
 # Prerequisites
 *.d
 

src/discretizations/finiteElementMethod/FiniteElementSpace.hpp

@@ -0,0 +1,30 @@
+#pragma once
+
+
+
+template< typename PARENT_ELEMENT >
+class FiniteElementSpace
+{
+public:
+
+    /// Alias for the parent element type
+    using ParentElementType = PARENT_ELEMENT;
+
+    /// Alias for the floating point type
+    using RealType = typename ParentElementType::RealType;
+
+    /// Alias for the type that represents a coordinate
+    using CoordType = typename ParentElementType::CoordType;
+
+    /// Alias for the type that represents a coordinate
+    using IndexType = typename ParentElementType::IndexType;
+
+    /// Alias for the type that represents a coordinate
+    using JacobianType = typename ParentElementType::JacobianType;
+
+    /// Alias for the type that represents a coordinate
+    using BasisFunctionType = typename ParentElementType::BasisFunctionType;
+
+
+
+};

src/discretizations/unitTests/CMakeLists.txt

@@ -3,6 +3,7 @@
 #
 
 set( unit_tests_sources
+     testFiniteElementSpace.cpp
      testParentElement.cpp
    )
 

src/discretizations/unitTests/testFiniteElementSpace.cpp

@@ -0,0 +1,187 @@
+
+#include "../finiteElementMethod/FiniteElementSpace.hpp"
+#include "../finiteElementMethod/parentElements/ParentElement.hpp"
+#include "functions/bases/LagrangeBasis.hpp"
+#include "functions/spacing/Spacing.hpp"
+#include "geometry/shapes/NCube.hpp"
+#include "common/ShivaMacros.hpp"
+#include "common/pmpl.hpp"
+
+
+
+#include <gtest/gtest.h>
+#include <cmath>
+
+using namespace shiva;
+using namespace shiva::discretizations::finiteElementMethod;
+using namespace shiva::functions;
+using namespace shiva::geometry;
+
+#include "testParentElementSolutions.hpp"
+
+
+
+template< typename TEST_PARENT_ELEMENT_HELPER >
+SHIVA_GLOBAL void compileTimeKernel()
+{
+  using ParentElementType = typename TEST_PARENT_ELEMENT_HELPER::ParentElementType;
+  constexpr int order = TEST_PARENT_ELEMENT_HELPER::order;
+
+
+  forSequence< order + 1 >( [] ( auto ica ) constexpr
+  {
+    constexpr int a = decltype(ica)::value;
+    forSequence< order + 1 >( [] ( auto icb ) constexpr
+    {
+      constexpr int b = decltype(icb)::value;
+      forSequence< order + 1 >( [] ( auto icc ) constexpr
+      {
+        constexpr int c = decltype(icc)::value;
+
+        constexpr double coord[3] = { TEST_PARENT_ELEMENT_HELPER::testParentCoords[0],
+                                      TEST_PARENT_ELEMENT_HELPER::testParentCoords[1],
+                                      TEST_PARENT_ELEMENT_HELPER::testParentCoords[2] };
+
+        constexpr double value = ParentElementType::template value< a, b, c >( coord );
+        constexpr CArray1d< double, 3 > gradient = ParentElementType::template gradient< a, b, c >( coord );
+        constexpr double tolerance = 1.0e-12;
+
+        static_assert( pmpl::check( value, TEST_PARENT_ELEMENT_HELPER::referenceValues[a][b][c], tolerance ) );
+        static_assert( pmpl::check( gradient[0], TEST_PARENT_ELEMENT_HELPER::referenceGradients[a][b][c][0], tolerance ) );
+        static_assert( pmpl::check( gradient[1], TEST_PARENT_ELEMENT_HELPER::referenceGradients[a][b][c][1], tolerance ) );
+        static_assert( pmpl::check( gradient[2], TEST_PARENT_ELEMENT_HELPER::referenceGradients[a][b][c][2], tolerance ) );
+      } );
+    } );
+  } );
+}
+
+template< typename TEST_PARENT_ELEMENT_HELPER >
+void testParentElementAtCompileTime()
+{
+#if defined(SHIVA_USE_DEVICE)
+  compileTimeKernel<TEST_PARENT_ELEMENT_HELPER><<<1,1>>>();
+#else
+  compileTimeKernel<TEST_PARENT_ELEMENT_HELPER>();
+#endif
+}
+
+
+template< typename TEST_PARENT_ELEMENT_HELPER >
+SHIVA_GLOBAL void runTimeKernel( double * const values, 
+                                 double * const gradients )
+{
+  using ParentElementType = typename TEST_PARENT_ELEMENT_HELPER::ParentElementType;
+  constexpr int order = TEST_PARENT_ELEMENT_HELPER::order;
+  constexpr int N = order + 1;
+
+
+  double const coord[3] = { TEST_PARENT_ELEMENT_HELPER::testParentCoords[0],
+                            TEST_PARENT_ELEMENT_HELPER::testParentCoords[1],
+                            TEST_PARENT_ELEMENT_HELPER::testParentCoords[2] };
+
+  forSequence< N >( [&] ( auto const ica ) constexpr
+  {
+    constexpr int a = decltype(ica)::value;
+    forSequence< N >( [&] ( auto const icb ) constexpr
+    {
+    constexpr int b = decltype(icb)::value;
+      forSequence< N >( [&] ( auto const icc ) constexpr
+      {
+        constexpr int c = decltype(icc)::value;
+        double const value = ParentElementType::template value< a, b, c >( coord );
+        CArray1d< double, 3 > const gradient = ParentElementType::template gradient< a, b, c >( coord );
+
+        values[ a*N*N + b*N + c ] = value;
+        gradients[ 3*(a*N*N + b*N + c) + 0 ] = gradient[0];
+        gradients[ 3*(a*N*N + b*N + c) + 1 ] = gradient[1];
+        gradients[ 3*(a*N*N + b*N + c) + 2 ] = gradient[2];
+      } );
+    } );
+  } );
+}
+
+template< typename TEST_PARENT_ELEMENT_HELPER >
+void testParentElementAtRunTime()
+{
+  constexpr int order = TEST_PARENT_ELEMENT_HELPER::order;
+  constexpr int N = order + 1;
+
+#if defined(SHIVA_USE_DEVICE)
+  constexpr int bytes = N*N*N*sizeof(double);
+  double * values;
+  double * gradients;
+  deviceMallocManaged( &values, bytes );
+  deviceMallocManaged( &gradients, 3*bytes );
+  runTimeKernel<TEST_PARENT_ELEMENT_HELPER><<<1,1>>>( values, gradients );
+  deviceDeviceSynchronize();
+#else
+  double values[N*N*N];
+  double gradients[N*N*N*3];
+  runTimeKernel<TEST_PARENT_ELEMENT_HELPER>( values, gradients );
+#endif
+
+  constexpr double tolerance = 1.0e-12;
+  for( int a=0; a<N; ++a )
+  {
+    for( int b=0; b<N; ++b )
+    {
+      for( int c=0; c<N; ++c)
+      {
+        EXPECT_NEAR( values[ a*N*N + b*N + c ], TEST_PARENT_ELEMENT_HELPER::referenceValues[a][b][c], fabs( TEST_PARENT_ELEMENT_HELPER::referenceValues[a][b][c] * tolerance ) );
+        EXPECT_NEAR( gradients[ 3*(a*N*N + b*N + c) + 0 ], TEST_PARENT_ELEMENT_HELPER::referenceGradients[a][b][c][0], fabs( TEST_PARENT_ELEMENT_HELPER::referenceGradients[a][b][c][0] * tolerance ) );
+        EXPECT_NEAR( gradients[ 3*(a*N*N + b*N + c) + 1 ], TEST_PARENT_ELEMENT_HELPER::referenceGradients[a][b][c][1], fabs( TEST_PARENT_ELEMENT_HELPER::referenceGradients[a][b][c][1] * tolerance ) );
+        EXPECT_NEAR( gradients[ 3*(a*N*N + b*N + c) + 2 ], TEST_PARENT_ELEMENT_HELPER::referenceGradients[a][b][c][2], fabs( TEST_PARENT_ELEMENT_HELPER::referenceGradients[a][b][c][2] * tolerance ) );
+      }
+    }
+  }
+
+}
+
+
+
+TEST( testParentElement, testBasicUsage )
+{
+
+  ParentElement< double,
+                 Cube,
+                 LagrangeBasis< double, 1, GaussLobattoSpacing >,
+                 LagrangeBasis< double, 1, GaussLobattoSpacing >,
+                 LagrangeBasis< double, 1, GaussLobattoSpacing >
+                 >::template value< 1, 1, 1 >( {0.5, 0, 1} );
+}
+
+TEST( testParentElement, testCubeLagrangeBasisGaussLobatto_O1 )
+{
+  using ParentElementType = ParentElement< double,
+                                           Cube,
+                                           LagrangeBasis< double, 1, GaussLobattoSpacing >,
+                                           LagrangeBasis< double, 1, GaussLobattoSpacing >,
+                                           LagrangeBasis< double, 1, GaussLobattoSpacing >
+                                           >;
+  using TEST_PARENT_ELEMENT_HELPER = TestParentElementHelper< ParentElementType >;
+
+  testParentElementAtCompileTime< TEST_PARENT_ELEMENT_HELPER >();
+  testParentElementAtRunTime< TEST_PARENT_ELEMENT_HELPER >();
+}
+
+TEST( testParentElement, testCubeLagrangeBasisGaussLobatto_O3 )
+{
+  using ParentElementType = ParentElement< double,
+                                           Cube,
+                                           LagrangeBasis< double, 3, GaussLobattoSpacing >,
+                                           LagrangeBasis< double, 3, GaussLobattoSpacing >,
+                                           LagrangeBasis< double, 3, GaussLobattoSpacing >
+                                           >;
+  using TEST_PARENT_ELEMENT_HELPER = TestParentElementHelper< ParentElementType >;
+
+ testParentElementAtCompileTime< TEST_PARENT_ELEMENT_HELPER >();
+ testParentElementAtRunTime< TEST_PARENT_ELEMENT_HELPER >();
+}
+
+
+int main( int argc, char * * argv )
+{
+  ::testing::InitGoogleTest( &argc, argv );
+  int const result = RUN_ALL_TESTS();
+  return result;
+}

src/geometry/shapes/InterpolatedShape.hpp

@@ -70,6 +70,7 @@ class InterpolatedShape
   /// Instantiation of the integer sequence that holds the number of support points in each basis.
   static inline constexpr numSupportPointsSequence basisSupportCounts{};
 
+  static inline constexpr int numVerticesInBasis[numDims] = {BASIS_TYPE::numSupportPoints...};
 
 
   static_assert( numDims == StandardGeom::numDims(), "numDims mismatch between cell and number of basis specified" );