diff --git a/include/shiva/geometry/mapping/Scaling.hpp b/include/shiva/geometry/mapping/Scaling.hpp index 0cbccfe..30d2e4b 100644 --- a/include/shiva/geometry/mapping/Scaling.hpp +++ b/include/shiva/geometry/mapping/Scaling.hpp @@ -20,6 +20,8 @@ #include "shiva/common/ShivaMacros.hpp" #include "shiva/common/types.hpp" #include "shiva/common/CArray.hpp" + +#include namespace shiva { @@ -63,10 +65,16 @@ class Scaling /** * @brief Returns the length dimensions of the rectangular cuboid. - * @return The length dimensions of the rectangular cuboid. + * @return A const reference to the stored length dimensions. */ constexpr SHIVA_HOST_DEVICE SHIVA_FORCE_INLINE DataType const & getData() const { return m_length; } + /** + * @brief Provides mutable access to the length dimensions of the rectangular cuboid. + * @return A reference to the stored length dimensions. + */ + constexpr SHIVA_HOST_DEVICE SHIVA_FORCE_INLINE DataType & getData() { return m_length; } + /** * @brief provides a reference to the member data. * @return a mutable reference to the member data. @@ -98,8 +106,8 @@ namespace utilities /** * @brief Calculates the Jacobian transformation for a rectangular cuboid. * @tparam REAL_TYPE The floating point type. - * @param cell The rectangular cuboid for which the Jacobian is calculated. - * @param J The Jacobian transformation operator. + * @param[in] cell The rectangular cuboid for which the Jacobian is calculated. + * @param[out] J The diagonal Jacobian entries, overwritten with half the current lengths. */ template< typename REAL_TYPE > SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE void @@ -113,6 +121,31 @@ jacobian( Scaling< REAL_TYPE > const & cell, } +/** + * @brief Calculates the constant Jacobian, independent of quadrature arguments. + * @tparam QUADRATURE The quadrature type, ignored because the Jacobian is constant. + * @tparam QA The quadrature indices, ignored for any index pack length, including zero. + * @tparam REAL_TYPE The floating point type. + * @param[in] cell The rectangular cuboid for which the Jacobian is calculated. + * @param[out] J The diagonal Jacobian entries, overwritten with half the current lengths. + * + * The constraint excludes QUADRATURE equal to REAL_TYPE so explicit scalar calls + * select the non-quadrature overload: + * @code + * jacobian< REAL_TYPE >( cell, J ); + * @endcode + */ +template< typename QUADRATURE, + int ... QA, + typename REAL_TYPE, + std::enable_if_t< !std::is_same_v< QUADRATURE, REAL_TYPE >, int > = 0 > +SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE void +jacobian( Scaling< REAL_TYPE > const & cell, + typename Scaling< REAL_TYPE >::JacobianType & J ) +{ + jacobian( cell, J ); +} + /** * @brief Calculates the inverse Jacobian transformation and detJ for a * rectangular cuboid. diff --git a/tests/geometry/mapping/testScaling.cpp b/tests/geometry/mapping/testScaling.cpp index 4f14fc9..16fd29b 100644 --- a/tests/geometry/mapping/testScaling.cpp +++ b/tests/geometry/mapping/testScaling.cpp @@ -15,8 +15,10 @@ #include "shiva/geometry/mapping/Scaling.hpp" #include "shiva/geometry/mapping/JacobianTransforms.hpp" #include "shiva/common/pmpl.hpp" +#include "shiva/functions/quadrature/Quadrature.hpp" #include +#include using namespace shiva; using namespace shiva::geometry; @@ -83,6 +85,110 @@ TEST( testScaling, testJacobianFunctionModifyLvalueRefArg ) testJacobianFunctionModifyLvalueRefArgHelper(); } +template< typename REAL_TYPE > +SHIVA_CONSTEXPR_HOSTDEVICE_FORCEINLINE bool matchesExpectedValue( REAL_TYPE const actual, int const expected ) +{ + REAL_TYPE const difference = actual - expected; + return difference > -1e-6 && difference < 1e-6; +} + +// Use void for deduction, REAL_TYPE for an explicit scalar argument, or a quadrature type. +// Initialize all entries to sentinels to verify that each overload overwrites the output. +template< typename REAL_TYPE, typename QUADRATURE = void, int ... QA > +SHIVA_CONSTEXPR_HOSTDEVICE_FORCEINLINE bool checkJacobianOverload( Scaling< REAL_TYPE > const & cell ) +{ + typename Scaling< REAL_TYPE >::JacobianType J{}; + J( 0 ) = -1; + J( 1 ) = -2; + J( 2 ) = -3; + if constexpr ( std::is_same_v< QUADRATURE, void > ) + { + jacobian( cell, J ); + } + else + { + jacobian< QUADRATURE, QA ... >( cell, J ); + } + return matchesExpectedValue( J( 0 ), 5 ) && matchesExpectedValue( J( 1 ), 10 ) && matchesExpectedValue( J( 2 ), 15 ); +} + +// Share the same interface checks between constant evaluation and host/device execution. +template< typename REAL_TYPE > +SHIVA_CONSTEXPR_HOSTDEVICE_FORCEINLINE void checkScalingInterface( bool * const checks ) +{ + using Cell = Scaling< REAL_TYPE >; + using Quadrature = QuadratureGaussLegendre< REAL_TYPE, 2 >; + constexpr REAL_TYPE h[3] = { 10, 20, 30 }; + auto cell = makeScaling( h ); + Cell const & constCell = cell; + static_assert( std::is_same_v< decltype(cell.getData()), typename Cell::DataType & > ); + static_assert( std::is_same_v< decltype(constCell.getData()), typename Cell::DataType const & > ); + + // Deduced and explicit scalar calls, then empty, partial, and full quadrature index packs. + checks[0] = checkJacobianOverload< REAL_TYPE >( cell ); + checks[1] = checkJacobianOverload< REAL_TYPE, REAL_TYPE >( cell ); + checks[2] = checkJacobianOverload< REAL_TYPE, Quadrature >( cell ); + checks[3] = checkJacobianOverload< REAL_TYPE, Quadrature, 0 >( cell ); + checks[4] = checkJacobianOverload< REAL_TYPE, Quadrature, 0, 1 >( cell ); + checks[5] = checkJacobianOverload< REAL_TYPE, Quadrature, 0, 0, 0 >( cell ); + checks[6] = checkJacobianOverload< REAL_TYPE, Quadrature, 1, 1, 1 >( cell ); + checks[7] = checkJacobianOverload< REAL_TYPE, Quadrature, 0, 1, 0 >( cell ); + + // Both const access and Jacobian evaluation must observe mutations through the getter. + auto & lengths = cell.getData(); + lengths[0] = 12; + lengths[1] = 24; + lengths[2] = 36; + checks[8] = matchesExpectedValue( constCell.getData()[0], 12 ) && + matchesExpectedValue( constCell.getData()[1], 24 ) && + matchesExpectedValue( constCell.getData()[2], 36 ); + typename Cell::JacobianType J{}; + jacobian< Quadrature, 1, 0, 1 >( constCell, J ); + checks[9] = matchesExpectedValue( J( 0 ), 6 ) && matchesExpectedValue( J( 1 ), 12 ) && matchesExpectedValue( J( 2 ), 18 ); +} + +template< typename REAL_TYPE > +constexpr bool checkScalingInterfaceConstexpr() +{ + bool checks[10]{}; + checkScalingInterface< REAL_TYPE >( checks ); + for( bool const check : checks ) + { + if( !check ) + { + return false; + } + } + return true; +} + +static_assert( checkScalingInterfaceConstexpr< float >() ); +static_assert( checkScalingInterfaceConstexpr< double >() ); + +template< typename REAL_TYPE > +void testScalingInterfaceHelper() +{ + bool checks[10]{}; + pmpl::genericKernelWrapper( 10, checks, [ = ] SHIVA_HOST_DEVICE ( bool * const kdata ) + { + checkScalingInterface< REAL_TYPE >( kdata ); + } ); + for( int i = 0; i < 10; ++i ) + { + EXPECT_TRUE( checks[i] ) << "Scaling interface check " << i; + } +} + +TEST( testScaling, testScalingInterfaceFloat ) +{ + testScalingInterfaceHelper< float >(); +} + +TEST( testScaling, testScalingInterfaceDouble ) +{ + testScalingInterfaceHelper< double >(); +} + void testJacobianFunctionReturnByValueHelper() { constexpr double h[3] = { 10, 20, 30 };