unify interface

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 <type_traits>
 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.

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 <gtest/gtest.h>
+#include <type_traits>
 
 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 };