diff --git a/.gitignore b/.gitignore index 7dbe5bd..68641e9 100644 --- a/.gitignore +++ b/.gitignore @@ -4,6 +4,12 @@ install-* .vscode +#vscode +.vscode + +#mac directories +.DS_Store + # Prerequisites *.d diff --git a/src/discretizations/finiteElementMethod/FiniteElementSpace.hpp b/src/discretizations/finiteElementMethod/FiniteElementSpace.hpp new file mode 100644 index 0000000..fbe3dc9 --- /dev/null +++ b/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; + + + +}; \ No newline at end of file diff --git a/src/discretizations/unitTests/CMakeLists.txt b/src/discretizations/unitTests/CMakeLists.txt index 2644351..c4666d4 100644 --- a/src/discretizations/unitTests/CMakeLists.txt +++ b/src/discretizations/unitTests/CMakeLists.txt @@ -3,6 +3,7 @@ # set( unit_tests_sources + testFiniteElementSpace.cpp testParentElement.cpp ) diff --git a/src/discretizations/unitTests/testFiniteElementSpace.cpp b/src/discretizations/unitTests/testFiniteElementSpace.cpp new file mode 100644 index 0000000..451f7b1 --- /dev/null +++ b/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 +#include + +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<<<1,1>>>(); +#else + compileTimeKernel(); +#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<<<1,1>>>( values, gradients ); + deviceDeviceSynchronize(); +#else + double values[N*N*N]; + double gradients[N*N*N*3]; + runTimeKernel( values, gradients ); +#endif + + constexpr double tolerance = 1.0e-12; + for( int a=0; a, + 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; +} diff --git a/src/geometry/shapes/InterpolatedShape.hpp b/src/geometry/shapes/InterpolatedShape.hpp index 28afddb..39e8c7a 100644 --- a/src/geometry/shapes/InterpolatedShape.hpp +++ b/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" );