Add spline quadrature and norm calculations.

Simplify the code for quadrature_coeffs_nd. Add a quadrature which is equivalent to calculating a spline approximation of a function and integrating it. Add functions which take a quadrature instance and the values of a function at the grid points and return the norm of that function. Add documentation to SplineBuilder2D and add accessor functions.

src/quadrature/compute_norms.hpp

@@ -0,0 +1,106 @@
+/**
+ * @file compute_norms.hpp
+ * File providing the L1 and the L2 norms.
+ */
+
+#pragma once
+
+#include <ddc/ddc.hpp>
+
+#include "geometry.hpp"
+#include "quadrature.hpp"
+
+
+/**
+ * @brief Compute L1 norm of a function with a given quadrature.
+ *
+ * @f$ \int_{\Omega} |f(X)|  dX @f$
+ *
+ * @param[in] quadrature
+ *      The quadrature used to compute the integral.
+ * @param[in] function
+ *      A ChunkSpan to the value of the function on the quadrature grid.
+ *
+ * @return A double containing the L1 norm of the function.
+ */
+template <class... IDim>
+double compute_L1_norm(
+        Quadrature<IDim...>& quadrature,
+        ddc::ChunkSpan<double, ddc::DiscreteDomain<IDim...>> function)
+{
+    auto grid = ddc::get_domain<IDim...>(function);
+    ddc::Chunk<double, ddc::DiscreteDomain<IDim...>> abs_function(grid);
+    ddc::for_each(grid, [&](ddc::DiscreteElement<IDim...> const idx) {
+        abs_function(idx) = fabs(function(idx));
+    });
+
+    return quadrature(abs_function);
+}
+
+
+
+/**
+ * @brief Compute L2 norm of a function with a given quadrature.
+ *
+ * @f$ \sqrt{\int_{\Omega} |f(X)|^2  dX} @f$
+ *
+ * @param[in] quadrature
+ *      The quadrature used to compute the integral.
+ * @param[in] function
+ *      A ChunkSpan to the value of the function on the quadrature grid.
+ *
+ * @return A double containing the L2 norm of the function.
+ */
+template <class... IDim>
+double compute_L2_norm(
+        Quadrature<IDim...>& quadrature,
+        ddc::ChunkSpan<double, ddc::DiscreteDomain<IDim...>> function)
+{
+    auto grid = ddc::get_domain<IDim...>(function);
+    ddc::Chunk<double, ddc::DiscreteDomain<IDim...>> square_function(grid);
+    ddc::for_each(grid, [&](ddc::DiscreteElement<IDim...> const idx) {
+        square_function(idx) = function(idx) * function(idx);
+    });
+
+    return std::sqrt(quadrature(square_function));
+}
+
+
+
+/**
+ * @brief Add the Jacobian determinant to the coefficients.
+ *
+ * For polar integration, we can add the Jacobian determinant to the
+ * quadrature coefficients.
+ *
+ * - 2D example: (but it is implemented for ND)
+ *
+ * @f$ \int_{\Omega_{\text{cart}}} f(x,y) dxdy
+ *  = \int_{\Omega} f(r,\theta) |det(J(r,\theta))| drd\theta
+ *  \simeq \sum_{i,j} [q_{i,j}| det(J(r_i,\theta_j))|] f_{ij}@f$
+ *
+ * This function uses rvalues. It means that coefficients is a temporary input
+ * parameter and it returns a temporary coefficient object. The Quadrature
+ * object can only be instantiate with rvalues.
+ *
+ * @param[in] mapping
+ *      The mapping function from the logical domain @f$ (r,\theta) @f$
+ *      to the physical domain @f$ (x, y) @f$.
+ * @param[in, out] coefficients
+ *      The quadrature coefficients @f$\{q_{ij}\}_{ij} @f$.
+ *
+ * @return A rvalue Chunk to the modified coefficients  @f$\{q_{ij}| det(J(r_i,\theta_j))|\}_{ij} @f$.
+ */
+template <class Mapping, class... IDim>
+ddc::Chunk<double, ddc::DiscreteDomain<IDim...>>&& compute_coeffs_on_mapping(
+        Mapping& mapping,
+        ddc::Chunk<double, ddc::DiscreteDomain<IDim...>>&& coefficients)
+{
+    ddc::DiscreteDomain<IDim...> grid = ddc::get_domain<IDim...>(coefficients);
+    ddc::for_each(grid, [&](ddc::DiscreteElement<IDim...> const idx) {
+        ddc::Coordinate<typename Mapping::curvilinear_tag_r, typename Mapping::curvilinear_tag_p>
+                coord(ddc::coordinate(idx));
+        coefficients(idx) *= mapping.jacobian(coord);
+    });
+    return std::move(coefficients);
+}

src/quadrature/quadrature_coeffs_nd.hpp

@@ -5,68 +5,43 @@
  */
 #pragma once
 
+#include <ddc/ddc.hpp>
+
+
+namespace {
+template <class IDim>
+using CoefficientChunk1D = ddc::Chunk<double, ddc::DiscreteDomain<IDim>>;
+}
 
 /**
  * @brief Helper function which creates ND dimensions from N 1D quadrature coefficient functions.
  *
  * @param domain
- * 	The domain on which the coefficients will be defined.
- * @param func
- * 	The function which defines quadrature coefficients in the first dimension.
- * @param o_funcs
- * 	The functions which define quadrature coefficients in the subsequent dimensions.
+ *      The domain on which the coefficients will be defined.
+ * @param funcs
+ *      The functions which define quadrature coefficients in the different dimensions.
  *
  * @returns The coefficients which define the quadrature method in ND.
  */
-template <class IDim, class... ODims>
-ddc::Chunk<double, ddc::DiscreteDomain<IDim, ODims...>> quadrature_coeffs_nd(
-        ddc::DiscreteDomain<IDim, ODims...> const& domain,
-        std::function<ddc::Chunk<double, ddc::DiscreteDomain<IDim>>(ddc::DiscreteDomain<IDim>)>
-                func,
-        std::function<ddc::Chunk<double, ddc::DiscreteDomain<ODims>>(
-                ddc::DiscreteDomain<ODims>)>... o_funcs)
+template <class... DDims>
+ddc::Chunk<double, ddc::DiscreteDomain<DDims...>> quadrature_coeffs_nd(
+        ddc::DiscreteDomain<DDims...> const& domain,
+        std::function<ddc::Chunk<double, ddc::DiscreteDomain<DDims>>(
+                ddc::DiscreteDomain<DDims>)>... funcs)
 {
-    // Get domains
-    ddc::DiscreteDomain<IDim> const current_dim_domain = ddc::select<IDim>(domain);
-    ddc::DiscreteDomain<ODims...> const other_dims_domain = ddc::select<ODims...>(domain);
-
-    // Get coefficients in the first dimension
-    ddc::Chunk<double, ddc::DiscreteDomain<IDim>> current_dim_coeffs = func(current_dim_domain);
-    // Get coefficients over all other dimensions
-    ddc::Chunk<double, ddc::DiscreteDomain<ODims...>> other_dim_coeffs
-            = quadrature_coeffs_nd(other_dims_domain, o_funcs...);
-
-    using DElemC = ddc::DiscreteElement<IDim>;
-    using DElemO = ddc::DiscreteElement<ODims...>;
-
-    // Calculate the combined coefficients over all dimensions
-    ddc::Chunk<double, ddc::DiscreteDomain<IDim, ODims...>> coefficients(domain);
-
-    ddc::for_each(current_dim_domain, [&](DElemC const idim) {
-        ddc::for_each(other_dims_domain, [&](DElemO const odim) {
-            coefficients(idim, odim) = current_dim_coeffs(idim) * other_dim_coeffs(odim);
-        });
+    // Get coefficients for each dimension
+    std::tuple<CoefficientChunk1D<DDims>...> current_dim_coeffs(
+            funcs(ddc::select<DDims>(domain))...);
+
+    // Allocate ND coefficients
+    ddc::Chunk<double, ddc::DiscreteDomain<DDims...>> coefficients(domain);
+
+    ddc::for_each(domain, [&](ddc::DiscreteElement<DDims...> const idim) {
+        // multiply the 1D coefficients by one another
+        coefficients(idim)
+                = (std::get<CoefficientChunk1D<DDims>>(current_dim_coeffs)(ddc::select<DDims>(idim))
+                   * ... * 1);
     });
 
     return coefficients;
 }
-
-/**
- * @brief Specialised 1D version of quadrature_coeffs_nd<IDim, ODims...>.
- *
- * @param domain
- * 	The domain on which the coefficients will be defined.
- * @param func
- * 	The function which defines quadrature coefficients.
- *
- * @returns The coefficients which define the quadrature method in 1D.
- */
-#include <ddc/ddc.hpp>
-template <class IDim>
-ddc::Chunk<double, ddc::DiscreteDomain<IDim>> quadrature_coeffs_nd(
-        ddc::DiscreteDomain<IDim> const& domain,
-        std::function<ddc::Chunk<double, ddc::DiscreteDomain<IDim>>(ddc::DiscreteDomain<IDim>)>
-                func)
-{
-    return func(domain);
-}

src/quadrature/spline_quadrature.hpp

@@ -0,0 +1,86 @@
+// SPDX-License-Identifier: MIT
+
+/**
+ * @file spline_quadrature.hpp
+ * File providing quadrature coefficients via a spline quadrature.
+ */
+
+#pragma once
+
+#include <cassert>
+
+#include <ddc/ddc.hpp>
+
+#include <sll/bsplines_non_uniform.hpp>
+#include <sll/matrix.hpp>
+#include <sll/null_boundary_value.hpp>
+#include <sll/spline_builder.hpp>
+#include <sll/spline_evaluator.hpp>
+
+
+
+namespace {
+template <class IDim>
+using CoefficientChunk1D = ddc::Chunk<double, ddc::DiscreteDomain<IDim>>;
+}
+
+
+/**
+ * @brief Get the spline quadrature coefficients in 1D.
+ *
+ * This function mainly uses the SplineBuilder::quadrature_coefficients function.
+ *
+ * @param[in] domain
+ *      The domain on which the splines quadrature will be carried out.
+ * @param[in] builder
+ *      The spline builder used for the quadrature coefficients.
+ *
+ * @return The quadrature coefficients for the trapezoid method defined on the provided domain.
+ */
+template <class IDim, class SplineBuilder>
+ddc::Chunk<double, ddc::DiscreteDomain<IDim>> spline_quadrature_coefficients_1d(
+        ddc::DiscreteDomain<IDim> const& domain,
+        SplineBuilder&& builder)
+{
+    return builder.quadrature_coefficients(domain);
+}
+
+
+
+/**
+ * @brief Get the spline quadrature coefficients in ND from N 1D quadrature coefficient.
+ *
+ * Calculate the quadrature coefficients for the spline quadrature method defined on the provided domain.
+ *
+ * @param[in] domain
+ *      The domain on which the coefficients will be defined.
+ * @param[in] builders
+ *      The spline builder used for the quadrature coefficients in the different dimensions.
+ *
+ * @return The coefficients which define the spline quadrature method in ND.
+ */
+template <class... DDims, class... SplineBuilders>
+ddc::Chunk<double, ddc::DiscreteDomain<DDims...>> spline_quadrature_coefficients(
+        ddc::DiscreteDomain<DDims...> const& domain,
+        SplineBuilders const&... builders)
+{
+    assert((std::is_same_v<
+                    typename DDims::continuous_dimension_type,
+                    typename SplineBuilders::bsplines_type::tag_type> and ...));
+
+    // Get coefficients for each dimension
+    std::tuple<CoefficientChunk1D<DDims>...> current_dim_coeffs(
+            spline_quadrature_coefficients_1d(ddc::select<DDims>(domain), builders)...);
+
+    // Allocate ND coefficients
+    ddc::Chunk<double, ddc::DiscreteDomain<DDims...>> coefficients(domain);
+
+    ddc::for_each(domain, [&](ddc::DiscreteElement<DDims...> const idim) {
+        // multiply the 1D coefficients by one another
+        coefficients(idim)
+                = (std::get<CoefficientChunk1D<DDims>>(current_dim_coeffs)(ddc::select<DDims>(idim))
+                   * ... * 1);
+    });
+
+    return coefficients;
+}

tests/geometryRTheta/CMakeLists.txt

@@ -2,3 +2,4 @@
 
 add_subdirectory(polar_poisson)
 add_subdirectory(2d_spline_interpolator)
+add_subdirectory(quadrature)

tests/geometryRTheta/quadrature/CMakeLists.txt

@@ -0,0 +1,18 @@
+cmake_minimum_required(VERSION 3.15)
+
+
+add_executable(Li_norms_spline_quadrature_tests
+	../../main.cpp
+	tests_L1_and_L2_norms.cpp
+)
+target_compile_features(Li_norms_spline_quadrature_tests PUBLIC cxx_std_17)
+target_link_libraries(Li_norms_spline_quadrature_tests
+    PUBLIC
+        GTest::gtest
+        GTest::gmock
+        DDC::DDC
+        sll::splines
+        vcx::geometry_RTheta
+        vcx::quadrature
+)
+gtest_discover_tests(Li_norms_spline_quadrature_tests)