From 83e51bee3090cd469955b2c498a8cfc431ead873 Mon Sep 17 00:00:00 2001 From: Randolph Settgast Date: Tue, 7 Oct 2025 10:06:46 -0700 Subject: [PATCH 1/2] add jacobian function for Scaling class to match LinearTransform interface --- include/shiva/geometry/mapping/Scaling.hpp | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/include/shiva/geometry/mapping/Scaling.hpp b/include/shiva/geometry/mapping/Scaling.hpp index 0cbccfe..0b1f7b0 100644 --- a/include/shiva/geometry/mapping/Scaling.hpp +++ b/include/shiva/geometry/mapping/Scaling.hpp @@ -67,6 +67,8 @@ class Scaling */ constexpr SHIVA_HOST_DEVICE SHIVA_FORCE_INLINE DataType const & getData() const { return m_length; } + 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. @@ -113,6 +115,19 @@ jacobian( Scaling< REAL_TYPE > const & cell, } +template< typename QUADRATURE, + int ... QA, + typename REAL_TYPE > +SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE void +jacobian( Scaling< REAL_TYPE > const & cell, + typename Scaling< REAL_TYPE >::JacobianType & J ) +{ + typename Scaling< REAL_TYPE >::DataType const & h = cell.getData(); + J( 0 ) = 0.5 * h[0]; + J( 1 ) = 0.5 * h[1]; + J( 2 ) = 0.5 * h[2]; +} + /** * @brief Calculates the inverse Jacobian transformation and detJ for a * rectangular cuboid. From 81c21e462da878d5dee78aa8379d8d00574c6f60 Mon Sep 17 00:00:00 2001 From: Randolph Settgast Date: Tue, 8 Sep 2026 18:08:18 -0700 Subject: [PATCH 2/2] Fix Scaling Jacobian overload compatibility --- include/shiva/geometry/mapping/Scaling.hpp | 34 +++++-- tests/geometry/mapping/testScaling.cpp | 106 +++++++++++++++++++++ 2 files changed, 132 insertions(+), 8 deletions(-) diff --git a/include/shiva/geometry/mapping/Scaling.hpp b/include/shiva/geometry/mapping/Scaling.hpp index 0b1f7b0..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,14 @@ 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; } /** @@ -100,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 @@ -115,17 +121,29 @@ 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 > + 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 ) { - typename Scaling< REAL_TYPE >::DataType const & h = cell.getData(); - J( 0 ) = 0.5 * h[0]; - J( 1 ) = 0.5 * h[1]; - J( 2 ) = 0.5 * h[2]; + jacobian( cell, J ); } /** 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 };