From 370fdad4d5b2fe979a45b2cdd4a17acd0571564f Mon Sep 17 00:00:00 2001 From: Randolph Settgast Date: Tue, 7 Nov 2023 19:06:52 -0800 Subject: [PATCH] move stuff into BasisProduct and InterpolatedShape classes from ParentElement --- src/discretizations/CMakeLists.txt | 2 +- .../parentElements/ParentElement.hpp | 147 ++------------ .../unitTests/testParentElement.cpp | 189 ------------------ src/functions/bases/BasisProduct.hpp | 170 ++++++++++++++++ src/functions/bases/LagrangeBasis.hpp | 11 +- src/functions/unitTests/testLagrangeBasis.cpp | 2 +- src/geometry/CMakeLists.txt | 1 + .../mapping/unitTests/testLinearTransform.cpp | 7 +- src/geometry/shapes/InterpolatedShape.hpp | 82 ++++++++ src/geometry/shapes/unitTests/CMakeLists.txt | 1 + .../unitTests/testInterpolatedShape.cpp | 188 +++++++++++++++++ .../testInterpolatedShapeSolutions.hpp} | 14 +- 12 files changed, 479 insertions(+), 335 deletions(-) delete mode 100644 src/discretizations/unitTests/testParentElement.cpp create mode 100644 src/functions/bases/BasisProduct.hpp create mode 100644 src/geometry/shapes/InterpolatedShape.hpp create mode 100644 src/geometry/shapes/unitTests/testInterpolatedShape.cpp rename src/{discretizations/unitTests/testParentElementSolutions.hpp => geometry/shapes/unitTests/testInterpolatedShapeSolutions.hpp} (94%) diff --git a/src/discretizations/CMakeLists.txt b/src/discretizations/CMakeLists.txt index a4c6741..96e06d2 100644 --- a/src/discretizations/CMakeLists.txt +++ b/src/discretizations/CMakeLists.txt @@ -27,4 +27,4 @@ install( TARGETS discretizations RUNTIME DESTINATION lib ) -add_subdirectory( unitTests ) +#add_subdirectory( unitTests ) diff --git a/src/discretizations/finiteElementMethod/parentElements/ParentElement.hpp b/src/discretizations/finiteElementMethod/parentElements/ParentElement.hpp index 16ed5ba..2555447 100644 --- a/src/discretizations/finiteElementMethod/parentElements/ParentElement.hpp +++ b/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 ) ); + // } diff --git a/src/discretizations/unitTests/testParentElement.cpp b/src/discretizations/unitTests/testParentElement.cpp deleted file mode 100644 index 3b453b2..0000000 --- a/src/discretizations/unitTests/testParentElement.cpp +++ /dev/null @@ -1,189 +0,0 @@ - -#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::discretizations::finiteElementMethod::basis; -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< double >, - 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< double >, - 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< double >, - 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/functions/bases/BasisProduct.hpp b/src/functions/bases/BasisProduct.hpp new file mode 100644 index 0000000..088524a --- /dev/null +++ b/src/functions/bases/BasisProduct.hpp @@ -0,0 +1,170 @@ +#pragma once + +#include "common/SequenceUtilities.hpp" + +namespace shiva +{ +namespace functions +{ + +/** + * @class BasisProduct + * @brief Defines a class that provides static functions to calculate quantities + * deriving from the product of basis functions. + * @tparam REAL_TYPE The floating point type to use + * @tparam BASIS_TYPES 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 ... BASIS_TYPES > +struct BasisProduct +{ + /// The number of dimensions/basis in the product + static inline constexpr int numDims = sizeof...(BASIS_TYPES); + + /// Alias for the floating point type + using RealType = REAL_TYPE; + + /// Alias for the type that represents a coordinate + using CoordType = REAL_TYPE[numDims]; + + + /** + * @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] + * + */ + 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_TYPES::template value< BASIS_FUNCTION_INDICES >( parentCoord[decltype(PRODUCT_TERM_INDEX)::value] ) * ... ); + } ); + +#endif + } + + /** + * @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] + */ + 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_TYPES, + 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_TYPES, + 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 + } + + +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] ) ); + } + } +}; + +} // namespace functions +} // namespace shiva \ No newline at end of file diff --git a/src/functions/bases/LagrangeBasis.hpp b/src/functions/bases/LagrangeBasis.hpp index fe7d10e..f99c06e 100644 --- a/src/functions/bases/LagrangeBasis.hpp +++ b/src/functions/bases/LagrangeBasis.hpp @@ -13,13 +13,10 @@ namespace shiva { -namespace discretizations -{ -namespace finiteElementMethod -{ -namespace basis +namespace functions { + /** * @class LagrangeBasis * @brief Defines a class to calculate quantities defined by a Lagrange @@ -219,7 +216,5 @@ class LagrangeBasis : public SPACING_TYPE< REAL_TYPE, ORDER + 1 > }; // class LagrangeBasis -} // namespace basis -} // namespace finiteElementMethod -} // namespace discretizations +} // namespace functions } // namespace shiva diff --git a/src/functions/unitTests/testLagrangeBasis.cpp b/src/functions/unitTests/testLagrangeBasis.cpp index 800a981..f65ef95 100644 --- a/src/functions/unitTests/testLagrangeBasis.cpp +++ b/src/functions/unitTests/testLagrangeBasis.cpp @@ -9,7 +9,7 @@ #include using namespace shiva; -using namespace shiva::discretizations::finiteElementMethod::basis; +using namespace shiva::functions; template< typename ... T > struct TestBasisHelper; diff --git a/src/geometry/CMakeLists.txt b/src/geometry/CMakeLists.txt index 27551a3..7ffbe74 100644 --- a/src/geometry/CMakeLists.txt +++ b/src/geometry/CMakeLists.txt @@ -4,6 +4,7 @@ set( geometry_headers mapping/LinearTransform.hpp mapping/Scaling.hpp mapping/UniformScaling.hpp + shapes/InterpolatedShape.hpp shapes/nCube.hpp ) diff --git a/src/geometry/mapping/unitTests/testLinearTransform.cpp b/src/geometry/mapping/unitTests/testLinearTransform.cpp index d85c66b..26f4431 100644 --- a/src/geometry/mapping/unitTests/testLinearTransform.cpp +++ b/src/geometry/mapping/unitTests/testLinearTransform.cpp @@ -4,7 +4,7 @@ #include "../JacobianTransforms.hpp" #include "functions/bases/LagrangeBasis.hpp" #include "functions/spacing/Spacing.hpp" -#include "discretizations/finiteElementMethod/parentElements/ParentElement.hpp" +#include "geometry/shapes/InterpolatedShape.hpp" #include "common/pmpl.hpp" @@ -14,8 +14,7 @@ using namespace shiva; using namespace shiva::geometry; using namespace shiva::geometry::utilities; -using namespace shiva::discretizations::finiteElementMethod; -using namespace shiva::discretizations::finiteElementMethod::basis; +using namespace shiva::functions; constexpr SHIVA_DEVICE double @@ -91,7 +90,7 @@ template< typename REAL_TYPE > SHIVA_HOST_DEVICE auto makeLinearTransform( REAL_TYPE const (&X)[8][3] ) { LinearTransform< REAL_TYPE, - ParentElement< double, + InterpolatedShape< double, Cube< double >, LagrangeBasis< double, 1, EqualSpacing >, LagrangeBasis< double, 1, EqualSpacing >, diff --git a/src/geometry/shapes/InterpolatedShape.hpp b/src/geometry/shapes/InterpolatedShape.hpp new file mode 100644 index 0000000..8639f71 --- /dev/null +++ b/src/geometry/shapes/InterpolatedShape.hpp @@ -0,0 +1,82 @@ +#pragma once + +#include "common/SequenceUtilities.hpp" +#include "common/ShivaMacros.hpp" +#include "common/types.hpp" + +#include "functions/bases/BasisProduct.hpp" + + +#include + +namespace shiva +{ +namespace geometry +{ + +/** + * @class InterpolatedShape + * @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_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_SHAPE, typename ... BASIS_TYPE > +class InterpolatedShape +{ + +public: + + /// The type used to represent the cell/geometry + using BaseShape = BASE_SHAPE; +// using FunctionalSpaceType = FUNCTIONAL_SPACE_TYPE; +// using IndexType = typename BaseShape::IndexType; + + /// Alias for the floating point type + using RealType = REAL_TYPE; + + /// Alias for the type that represents a coordinate + using CoordType = typename BaseShape::CoordType; + + /// The type used to represent the product of basis functions + using BASIS_PRODUCT_TYPE = functions::BasisProduct< REAL_TYPE, BASIS_TYPE... >; + + /// The number of dimensions on the InterpolatedShape + static inline constexpr int numDims = sizeof...(BASIS_TYPE); + + /// The number of vertices on the InterpolatedShape + static inline constexpr int numVertices = BaseShape::numVertices(); + + + static_assert( numDims == BaseShape::numDims(), "numDims mismatch between cell and number of basis specified" ); + + + /** + * @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 ( BASIS_PRODUCT_TYPE::template value< BASIS_FUNCTION_INDICES... >( parentCoord ) ); + } + + /** + * @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" ); + return ( BASIS_PRODUCT_TYPE::template gradient< BASIS_FUNCTION_INDICES... >( parentCoord ) ); + } +}; + + +} // namespace geometry +} // namespace shiva diff --git a/src/geometry/shapes/unitTests/CMakeLists.txt b/src/geometry/shapes/unitTests/CMakeLists.txt index 455d2d5..046c81c 100644 --- a/src/geometry/shapes/unitTests/CMakeLists.txt +++ b/src/geometry/shapes/unitTests/CMakeLists.txt @@ -4,6 +4,7 @@ set( unit_tests_sources testNCube.cpp + testInterpolatedShape.cpp ) set( dependencyList ${parallelDeps} gtest geometry ) diff --git a/src/geometry/shapes/unitTests/testInterpolatedShape.cpp b/src/geometry/shapes/unitTests/testInterpolatedShape.cpp new file mode 100644 index 0000000..24add7a --- /dev/null +++ b/src/geometry/shapes/unitTests/testInterpolatedShape.cpp @@ -0,0 +1,188 @@ + +#include "../InterpolatedShape.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::geometry; +using namespace shiva::functions; + +#include "testInterpolatedShapeSolutions.hpp" + + + +template< typename TEST_INTERPOLATED_SHAPE_HELPER > +SHIVA_GLOBAL void compileTimeKernel() +{ + using InterpolatedShapeType = typename TEST_INTERPOLATED_SHAPE_HELPER::InterpolatedShapeType; + constexpr int order = TEST_INTERPOLATED_SHAPE_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_INTERPOLATED_SHAPE_HELPER::testCoords[0], + TEST_INTERPOLATED_SHAPE_HELPER::testCoords[1], + TEST_INTERPOLATED_SHAPE_HELPER::testCoords[2] }; + + constexpr double value = InterpolatedShapeType::template value< a, b, c >( coord ); + constexpr CArray1d< double, 3 > gradient = InterpolatedShapeType::template gradient< a, b, c >( coord ); + constexpr double tolerance = 1.0e-12; + + static_assert( pmpl::check( value, TEST_INTERPOLATED_SHAPE_HELPER::referenceValues[a][b][c], tolerance ) ); + static_assert( pmpl::check( gradient[0], TEST_INTERPOLATED_SHAPE_HELPER::referenceGradients[a][b][c][0], tolerance ) ); + static_assert( pmpl::check( gradient[1], TEST_INTERPOLATED_SHAPE_HELPER::referenceGradients[a][b][c][1], tolerance ) ); + static_assert( pmpl::check( gradient[2], TEST_INTERPOLATED_SHAPE_HELPER::referenceGradients[a][b][c][2], tolerance ) ); + } ); + } ); + } ); +} + +template< typename TEST_INTERPOLATED_SHAPE_HELPER > +void testInterpolatedShapeAtCompileTime() +{ +#if defined(SHIVA_USE_DEVICE) + compileTimeKernel< TEST_INTERPOLATED_SHAPE_HELPER ><< < 1, 1 >> > (); +#else + compileTimeKernel< TEST_INTERPOLATED_SHAPE_HELPER >(); +#endif +} + + +template< typename TEST_INTERPOLATED_SHAPE_HELPER > +SHIVA_GLOBAL void runTimeKernel( double * const values, + double * const gradients ) +{ + using InterpolatedShapeType = typename TEST_INTERPOLATED_SHAPE_HELPER::InterpolatedShapeType; + constexpr int order = TEST_INTERPOLATED_SHAPE_HELPER::order; + constexpr int N = order + 1; + + + double const coord[3] = { TEST_INTERPOLATED_SHAPE_HELPER::testCoords[0], + TEST_INTERPOLATED_SHAPE_HELPER::testCoords[1], + TEST_INTERPOLATED_SHAPE_HELPER::testCoords[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 = InterpolatedShapeType::template value< a, b, c >( coord ); + CArray1d< double, 3 > const gradient = InterpolatedShapeType::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_INTERPOLATED_SHAPE_HELPER > +void testInterpolatedShapeAtRunTime() +{ + constexpr int order = TEST_INTERPOLATED_SHAPE_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_INTERPOLATED_SHAPE_HELPER ><< < 1, 1 >> > ( values, gradients ); + deviceDeviceSynchronize(); +#else + double values[N * N * N]; + double gradients[N * N * N * 3]; + runTimeKernel< TEST_INTERPOLATED_SHAPE_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_INTERPOLATED_SHAPE_HELPER::referenceValues[a][b][c], fabs( TEST_INTERPOLATED_SHAPE_HELPER::referenceValues[a][b][c] * tolerance ) ); + EXPECT_NEAR( gradients[ 3 * (a * N * N + b * N + c) + 0 ], TEST_INTERPOLATED_SHAPE_HELPER::referenceGradients[a][b][c][0], + fabs( TEST_INTERPOLATED_SHAPE_HELPER::referenceGradients[a][b][c][0] * tolerance ) ); + EXPECT_NEAR( gradients[ 3 * (a * N * N + b * N + c) + 1 ], TEST_INTERPOLATED_SHAPE_HELPER::referenceGradients[a][b][c][1], + fabs( TEST_INTERPOLATED_SHAPE_HELPER::referenceGradients[a][b][c][1] * tolerance ) ); + EXPECT_NEAR( gradients[ 3 * (a * N * N + b * N + c) + 2 ], TEST_INTERPOLATED_SHAPE_HELPER::referenceGradients[a][b][c][2], + fabs( TEST_INTERPOLATED_SHAPE_HELPER::referenceGradients[a][b][c][2] * tolerance ) ); + } + } + } + +} + + + +TEST( testInterpolatedShape, testBasicUsage ) +{ + + InterpolatedShape< double, + Cube< double >, + LagrangeBasis< double, 1, GaussLobattoSpacing >, + LagrangeBasis< double, 1, GaussLobattoSpacing >, + LagrangeBasis< double, 1, GaussLobattoSpacing > + >::template value< 1, 1, 1 >( {0.5, 0, 1} ); +} + +TEST( testInterpolatedShape, testCubeLagrangeBasisGaussLobatto_O1 ) +{ + using InterpolatedShapeType = InterpolatedShape< double, + Cube< double >, + LagrangeBasis< double, 1, GaussLobattoSpacing >, + LagrangeBasis< double, 1, GaussLobattoSpacing >, + LagrangeBasis< double, 1, GaussLobattoSpacing > + >; + using TEST_INTERPOLATED_SHAPE_HELPER = TestInterpolatedShapeHelper< InterpolatedShapeType >; + + testInterpolatedShapeAtCompileTime< TEST_INTERPOLATED_SHAPE_HELPER >(); + testInterpolatedShapeAtRunTime< TEST_INTERPOLATED_SHAPE_HELPER >(); +} + +TEST( testInterpolatedShape, testCubeLagrangeBasisGaussLobatto_O3 ) +{ + using InterpolatedShapeType = InterpolatedShape< double, + Cube< double >, + LagrangeBasis< double, 3, GaussLobattoSpacing >, + LagrangeBasis< double, 3, GaussLobattoSpacing >, + LagrangeBasis< double, 3, GaussLobattoSpacing > + >; + using TEST_INTERPOLATED_SHAPE_HELPER = TestInterpolatedShapeHelper< InterpolatedShapeType >; + + testInterpolatedShapeAtCompileTime< TEST_INTERPOLATED_SHAPE_HELPER >(); + testInterpolatedShapeAtRunTime< TEST_INTERPOLATED_SHAPE_HELPER >(); +} + + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + int const result = RUN_ALL_TESTS(); + return result; +} diff --git a/src/discretizations/unitTests/testParentElementSolutions.hpp b/src/geometry/shapes/unitTests/testInterpolatedShapeSolutions.hpp similarity index 94% rename from src/discretizations/unitTests/testParentElementSolutions.hpp rename to src/geometry/shapes/unitTests/testInterpolatedShapeSolutions.hpp index 9fdfc6a..a97def0 100644 --- a/src/discretizations/unitTests/testParentElementSolutions.hpp +++ b/src/geometry/shapes/unitTests/testInterpolatedShapeSolutions.hpp @@ -1,17 +1,17 @@ template< typename ... T > -struct TestParentElementHelper; +struct TestInterpolatedShapeHelper; template<> -struct TestParentElementHelper< ParentElement< double, +struct TestInterpolatedShapeHelper< InterpolatedShape< double, Cube< double >, LagrangeBasis< double, 1, GaussLobattoSpacing >, LagrangeBasis< double, 1, GaussLobattoSpacing >, LagrangeBasis< double, 1, GaussLobattoSpacing > > > { - using ParentElementType = ParentElement< double, + using InterpolatedShapeType = InterpolatedShape< double, Cube< double >, LagrangeBasis< double, 1, GaussLobattoSpacing >, LagrangeBasis< double, 1, GaussLobattoSpacing >, @@ -19,7 +19,7 @@ struct TestParentElementHelper< ParentElement< double, >; static constexpr int order = 1; - static constexpr double testParentCoords[3] = { 0.31415, -0.161803, 0.69314 }; + static constexpr double testCoords[3] = { 0.31415, -0.161803, 0.69314 }; static constexpr double referenceValues[2][2][2] = { { {0.030564122401949, 0.16864152448555}, {0.022050860348051, 0.12166849276445} }, { {0.058563594743051, 0.32313225836945}, {0.042251422506949, 0.23312772438055} } }; static constexpr double referenceGradients[2][2][2][3] = { { { {-0.0445638585725, -0.026307491375, -0.09960282344375}, {-0.2458868914275, -0.145155008625, 0.09960282344375} }, @@ -31,7 +31,7 @@ struct TestParentElementHelper< ParentElement< double, template<> -struct TestParentElementHelper< ParentElement< double, +struct TestInterpolatedShapeHelper< InterpolatedShape< double, Cube< double >, LagrangeBasis< double, 3, GaussLobattoSpacing >, LagrangeBasis< double, 3, GaussLobattoSpacing >, @@ -39,7 +39,7 @@ struct TestParentElementHelper< ParentElement< double, > > { - using ParentElementType = ParentElement< double, + using InterpolatedShapeType = InterpolatedShape< double, Cube< double >, LagrangeBasis< double, 3, GaussLobattoSpacing >, LagrangeBasis< double, 3, GaussLobattoSpacing >, @@ -47,7 +47,7 @@ struct TestParentElementHelper< ParentElement< double, >; static constexpr int order = 3; - static constexpr double testParentCoords[3] = { 0.31415, -0.161803, 0.69314}; + static constexpr double testCoords[3] = { 0.31415, -0.161803, 0.69314}; static constexpr double referenceValues[4][4][4] = { { {0.00029480662998955, -0.00097875857989107, 0.0045384751099593, 0.0016266339617432}, {-0.0019359683760398, 0.0064274187405959, -0.029803754035777, -0.01068195755787}, {-0.00090727558355123, 0.0030121566864284, -0.013967283076256, -0.0050060111501464}, {0.00021269185295386, -0.00070613736183238, 0.0032743384392556, 0.0011735549889536} },