[NL] made interpolation easier to use

output parameters are not to passed as array of addresses anymore.

NumLib/Function/Interpolation.h

@@ -16,38 +16,75 @@
 namespace NumLib
 {
 
-/**
+namespace detail
+{
+
+//! \see ::NumLib::shapeFunctionInterpolate()
+template<unsigned DOFOffset, typename NodalValues, typename ShapeMatrix>
+void shapeFunctionInterpolate(
+        const NodalValues &/*nodal_values*/,
+        const ShapeMatrix &/*shape_matrix_N*/)
+{}
+
+//! \see ::NumLib::shapeFunctionInterpolate()
+template<unsigned DOFOffset, typename NodalValues, typename ShapeMatrix, typename... ScalarTypes>
+void shapeFunctionInterpolate(
+        const NodalValues &nodal_values,
+        const ShapeMatrix &shape_matrix_N,
+        double& interpolated_value,
+        ScalarTypes&... interpolated_values)
+{
+    auto const num_nodes = shape_matrix_N.size();
+    double iv = 0.0;
+
+    for (auto n=decltype(num_nodes){0}; n<num_nodes; ++n) {
+        iv += nodal_values[DOFOffset*num_nodes+n] * shape_matrix_N[n];
+    }
+
+    interpolated_value = iv;
+
+    shapeFunctionInterpolate<DOFOffset+1>(nodal_values, shape_matrix_N, interpolated_values...);
+}
+
+} // namespace detail
+
+/*!
  * Interpolates variables given at element nodes according to the given shape matrix.
  *
  * This function simply does the usual finite-element interpolation, i.e. multiplication
  * of nodal values with the shape function.
  *
- * @param nodal_values   vector of nodal values, ordered by component
- * @param shape_matrix_N shape matrix of the point to which will be interpolated
- * @param interpolated_values array of addresses to which the interpolated values will be written
+ * \param nodal_values vector of nodal values, ordered by component
+ * \param shape_matrix_N shape matrix of the point to which will be interpolated
+ * \param interpolated_value interpolated value of the first d.o.f. (output parameter)
+ * \param interpolated_values interpolated value of further d.o.f. (output parameter)
+ *
+ * \tparam NodalValues  type of the container where nodal values are stored
+ * \tparam ShapeMatrix  type of the shape matrix \f$N\f$.
+ * \tparam ScalarValues all of the types in this pack must currently be \c double.
+ *
+ * \note
+ * \c nodal_values have to be ordered by component and it is assumed that all passed d.o.f. are
+ * single-component and are interpolated using the same shape function.
  */
-template<typename NodalValues, typename ShapeMatrix, std::size_t NodalDOF>
+template<typename NodalValues, typename ShapeMatrix, typename... ScalarTypes>
 void shapeFunctionInterpolate(
         const NodalValues& nodal_values,
         const ShapeMatrix& shape_matrix_N,
-        std::array<double*, NodalDOF> interpolated_values
+        double& interpolated_value,
+        ScalarTypes&... interpolated_values
         )
 {
+    auto const num_nodal_dof = sizeof...(interpolated_values) + 1;
     auto const num_nodes = shape_matrix_N.size();
 
-    assert(num_nodes*NodalDOF == nodal_values.size());
-
-    for (auto d=decltype(NodalDOF){0}; d<NodalDOF; ++d)
-    {
-        *interpolated_values[d] = 0.0;
+    assert(num_nodes*num_nodal_dof == nodal_values.size());
+    (void) num_nodal_dof; // no warnings when not in debug build
 
-        for (auto n=decltype(num_nodes){0}; n<num_nodes; ++n)
-        {
-            *interpolated_values[d] += nodal_values[d*num_nodes+n] * shape_matrix_N[n];
-        }
-    }
+    detail::shapeFunctionInterpolate<0>(nodal_values, shape_matrix_N, interpolated_value,
+                                        interpolated_values...);
 }
 
-}
+} // namespace NumLib
 
 #endif // NUMLIB_INTERPOLATION_H

Tests/NumLib/TestFunctionInterpolation.cpp

@@ -26,7 +26,6 @@ TEST(NumLibFunctionInterpolationTest, TwoVariablesTwoNodes)
 {
     double variable1 = 0.0;
     double variable2 = 0.0;
-    std::array<double*, 2> interpolated_values = {{&variable1, &variable2}};
 
     const std::array<double, 4> nodal_values = {{
         0.0, 1.0,  // for variable1
@@ -36,7 +35,7 @@ TEST(NumLibFunctionInterpolationTest, TwoVariablesTwoNodes)
     const std::array<double, 2> shape_matrix = {{0.25, 0.75}};
 
     NumLib::shapeFunctionInterpolate(nodal_values, shape_matrix,
-                                     interpolated_values);
+                                     variable1, variable2);
 
     ASSERT_EQ(0.75, variable1);
     ASSERT_EQ(0.5,  variable2);
@@ -78,7 +77,6 @@ TEST(NumLibFunctionInterpolationTest, Linear1DElement)
     // actual test
     double variable1 = 0.0;
     double variable2 = 0.0;
-    std::array<double*, 2> interpolated_values = {{&variable1, &variable2}};
 
     const std::array<double, 4> nodal_values = {{
         0.0, 1.0,  // for variable1
@@ -86,7 +84,7 @@ TEST(NumLibFunctionInterpolationTest, Linear1DElement)
     }};
 
     NumLib::shapeFunctionInterpolate(nodal_values, shape_matrix.N,
-                                     interpolated_values);
+                                     variable1, variable2);
 
     const double n0 = shape_matrix.N[0];
     const double n1 = shape_matrix.N[1];