Add EdgeIntegrator (and associated tests)

src/scheme/EdgeIntegrator.hpp

+#ifndef EDGE_INTEGRATOR_HPP
+#define EDGE_INTEGRATOR_HPP
+
+#include <analysis/IQuadratureDescriptor.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <geometry/LineTransformation.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemId.hpp>
+#include <utils/Array.hpp>
+
+#include <functional>
+
+class EdgeIntegrator
+{
+ private:
+  template <size_t Dimension>
+  class AllEdges
+  {
+   private:
+    const Connectivity<Dimension>& m_connectivity;
+
+   public:
+    using index_type = EdgeId;
+
+    PUGS_INLINE
+    EdgeId
+    edgeId(const EdgeId edge_id) const
+    {
+      return edge_id;
+    }
+
+    PUGS_INLINE
+    size_t
+    size() const
+    {
+      return m_connectivity.numberOfEdges();
+    }
+
+    PUGS_INLINE
+    AllEdges(const Connectivity<Dimension>& connectivity) : m_connectivity{connectivity} {}
+  };
+
+  template <template <typename T> typename ArrayT>
+  class EdgeList
+  {
+   private:
+    const ArrayT<EdgeId>& m_edge_list;
+
+   public:
+    using index_type = typename ArrayT<EdgeId>::index_type;
+
+    PUGS_INLINE
+    EdgeId
+    edgeId(const index_type index) const
+    {
+      return m_edge_list[index];
+    }
+
+    PUGS_INLINE
+    size_t
+    size() const
+    {
+      return m_edge_list.size();
+    }
+
+    PUGS_INLINE
+    EdgeList(const ArrayT<EdgeId>& edge_list) : m_edge_list{edge_list} {}
+  };
+
+  template <typename MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
+  static PUGS_INLINE void
+  _integrateTo(std::function<OutputType(InputType)> function,
+               const IQuadratureDescriptor& quadrature_descriptor,
+               const MeshType& mesh,
+               const ListTypeT& edge_list,
+               OutputArrayT& value)
+  {
+    constexpr size_t Dimension = MeshType::Dimension;
+    static_assert(Dimension == 3);
+
+    static_assert(std::is_same_v<TinyVector<Dimension>, std::decay_t<InputType>>, "invalid input data type");
+    static_assert(std::is_same_v<std::remove_const_t<typename OutputArrayT::data_type>, OutputType>,
+                  "invalid output data type");
+
+    using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
+    Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
+
+    if constexpr (std::is_arithmetic_v<OutputType>) {
+      value.fill(0);
+    } else if constexpr (is_tiny_vector_v<OutputType> or is_tiny_matrix_v<OutputType>) {
+      value.fill(zero);
+    } else {
+      static_assert(std::is_same_v<OutputType, double>, "unexpected output type");
+    }
+
+    const auto& connectivity = mesh.connectivity();
+
+    const auto edge_to_node_matrix = connectivity.edgeToNodeMatrix();
+    const auto xr                  = mesh.xr();
+
+    const auto qf = QuadratureManager::instance().getLineFormula(quadrature_descriptor);
+
+    parallel_for(edge_list.size(), [=, &qf, &edge_list](typename ListTypeT::index_type index) {
+      const EdgeId edge_id = edge_list.edgeId(index);
+
+      const auto edge_to_node = edge_to_node_matrix[edge_id];
+
+      Assert(edge_to_node.size() == 2);
+      const LineTransformation<3> T(xr[edge_to_node[0]], xr[edge_to_node[1]]);
+
+      for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+        const auto xi = qf.point(i_point);
+        value[index] += qf.weight(i_point) * T.velocityNorm() * function(T(xi));
+      }
+    });
+  }
+
+ public:
+  template <typename MeshType, typename OutputArrayT, typename OutputType, typename InputType>
+  static PUGS_INLINE void
+  integrateTo(const std::function<OutputType(InputType)>& function,
+              const IQuadratureDescriptor& quadrature_descriptor,
+              const MeshType& mesh,
+              OutputArrayT& value)
+  {
+    constexpr size_t Dimension = MeshType::Dimension;
+
+    _integrateTo<MeshType, OutputArrayT>(function, quadrature_descriptor, mesh,
+                                         AllEdges<Dimension>{mesh.connectivity()}, value);
+  }
+
+  template <typename MeshType, typename OutputArrayT, typename FunctionType>
+  static PUGS_INLINE void
+  integrateTo(const FunctionType& function,
+              const IQuadratureDescriptor& quadrature_descriptor,
+              const MeshType& mesh,
+              OutputArrayT& value)
+  {
+    integrateTo(std::function(function), quadrature_descriptor, mesh, value);
+  }
+
+  template <typename MeshType, typename OutputType, typename InputType, template <typename DataType> typename ArrayT>
+  static PUGS_INLINE ArrayT<OutputType>
+  integrate(const std::function<OutputType(InputType)>& function,
+            const IQuadratureDescriptor& quadrature_descriptor,
+            const MeshType& mesh,
+            const ArrayT<EdgeId>& edge_list)
+  {
+    ArrayT<OutputType> value(size(edge_list));
+    _integrateTo<MeshType, ArrayT<OutputType>>(function, quadrature_descriptor, mesh, EdgeList{edge_list}, value);
+
+    return value;
+  }
+
+  template <typename MeshType, typename FunctionType, template <typename DataType> typename ArrayT>
+  static PUGS_INLINE auto
+  integrate(const FunctionType& function,
+            const IQuadratureDescriptor& quadrature_descriptor,
+            const MeshType& mesh,
+            const ArrayT<EdgeId>& edge_list)
+  {
+    return integrate(std::function(function), quadrature_descriptor, mesh, edge_list);
+  }
+};
+
+#endif   // EDGE_INTEGRATOR_HPP

tests/CMakeLists.txt

@@ -70,6 +70,7 @@ add_executable (unit_tests
   test_DiscreteFunctionType.cpp
   test_DiscreteFunctionUtils.cpp
   test_DoWhileProcessor.cpp
+  test_EdgeIntegrator.cpp
   test_EigenvalueSolver.cpp
   test_EmbeddedData.cpp
   test_EmbeddedIDiscreteFunctionUtils.cpp

tests/test_EdgeIntegrator.cpp

+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/GaussLobattoQuadratureDescriptor.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <mesh/DiamondDualMeshManager.hpp>
+#include <mesh/ItemValue.hpp>
+#include <mesh/Mesh.hpp>
+#include <scheme/EdgeIntegrator.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("EdgeIntegrator", "[scheme]")
+{
+  SECTION("3D")
+  {
+    using R3 = TinyVector<3>;
+
+    auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh();
+
+    auto f = [](const R3& X) -> double {
+      const double x = X[0];
+      const double y = X[1];
+      const double z = X[2];
+      return x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1;
+    };
+
+    std::vector<std::pair<std::string, decltype(hybrid_mesh)>> mesh_list;
+    mesh_list.push_back(std::make_pair("hybrid mesh", hybrid_mesh));
+    mesh_list.push_back(
+      std::make_pair("diamond mesh", DiamondDualMeshManager::instance().getDiamondDualMesh(hybrid_mesh)));
+
+    for (auto mesh_info : mesh_list) {
+      auto mesh_name = mesh_info.first;
+      auto mesh      = mesh_info.second;
+
+      Array<const double> int_f_per_edge = [=] {
+        Array<double> int_f(mesh->numberOfEdges());
+        auto edge_to_node_matrix = mesh->connectivity().edgeToNodeMatrix();
+
+        parallel_for(
+          mesh->numberOfEdges(), PUGS_LAMBDA(const EdgeId edge_id) {
+            auto edge_node_list = edge_to_node_matrix[edge_id];
+            auto xr             = mesh->xr();
+            double integral     = 0;
+
+            LineTransformation<3> T(xr[edge_node_list[0]], xr[edge_node_list[1]]);
+            auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2));
+            for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+              const auto& xi = qf.point(i);
+              integral += qf.weight(i) * T.velocityNorm() * f(T(xi));
+            }
+
+            int_f[edge_id] = integral;
+          });
+
+        return int_f;
+      }();
+
+      SECTION(mesh_name)
+      {
+        SECTION("direct formula")
+        {
+          SECTION("all edges")
+          {
+            SECTION("EdgeValue")
+            {
+              EdgeValue<double> values(mesh->connectivity());
+              EdgeIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh,
+                                          values);
+
+              double error = 0;
+              for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+                error += std::abs(int_f_per_edge[edge_id] - values[edge_id]);
+              }
+
+              REQUIRE(error == Catch::Approx(0).margin(1E-10));
+            }
+
+            SECTION("Array")
+            {
+              Array<double> values(mesh->numberOfEdges());
+
+              EdgeIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+
+              double error = 0;
+              for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+                error += std::abs(int_f_per_edge[edge_id] - values[edge_id]);
+              }
+
+              REQUIRE(error == Catch::Approx(0).margin(1E-10));
+            }
+
+            SECTION("SmallArray")
+            {
+              SmallArray<double> values(mesh->numberOfEdges());
+              EdgeIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+
+              double error = 0;
+              for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+                error += std::abs(int_f_per_edge[edge_id] - values[edge_id]);
+              }
+
+              REQUIRE(error == Catch::Approx(0).margin(1E-10));
+            }
+          }
+
+          SECTION("edge list")
+          {
+            SECTION("Array")
+            {
+              Array<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
+
+              {
+                size_t k = 0;
+                for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
+                  edge_list[k] = edge_id;
+                }
+
+                REQUIRE(k == edge_list.size());
+              }
+
+              Array<double> values = EdgeIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, edge_list);
+
+              double error = 0;
+              for (size_t i = 0; i < edge_list.size(); ++i) {
+                error += std::abs(int_f_per_edge[edge_list[i]] - values[i]);
+              }
+
+              REQUIRE(error == Catch::Approx(0).margin(1E-10));
+            }
+
+            SECTION("SmallArray")
+            {
+              SmallArray<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
+
+              {
+                size_t k = 0;
+                for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
+                  edge_list[k] = edge_id;
+                }
+
+                REQUIRE(k == edge_list.size());
+              }
+
+              SmallArray<double> values = EdgeIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, edge_list);
+
+              double error = 0;
+              for (size_t i = 0; i < edge_list.size(); ++i) {
+                error += std::abs(int_f_per_edge[edge_list[i]] - values[i]);
+              }
+
+              REQUIRE(error == Catch::Approx(0).margin(1E-10));
+            }
+          }
+        }
+
+        SECTION("tensorial formula")
+        {
+          SECTION("all edges")
+          {
+            SECTION("EdgeValue")
+            {
+              EdgeValue<double> values(mesh->connectivity());
+              EdgeIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10),
+                                          *mesh, values);
+
+              double error = 0;
+              for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+                error += std::abs(int_f_per_edge[edge_id] - values[edge_id]);
+              }
+
+              REQUIRE(error == Catch::Approx(0).margin(1E-10));
+            }
+
+            SECTION("Array")
+            {
+              Array<double> values(mesh->numberOfEdges());
+
+              EdgeIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+
+              double error = 0;
+              for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+                error += std::abs(int_f_per_edge[edge_id] - values[edge_id]);
+              }
+
+              REQUIRE(error == Catch::Approx(0).margin(1E-10));
+            }
+
+            SECTION("SmallArray")
+            {
+              SmallArray<double> values(mesh->numberOfEdges());
+              EdgeIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+
+              double error = 0;
+              for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+                error += std::abs(int_f_per_edge[edge_id] - values[edge_id]);
+              }
+
+              REQUIRE(error == Catch::Approx(0).margin(1E-10));
+            }
+          }
+
+          SECTION("edge list")
+          {
+            SECTION("Array")
+            {
+              Array<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
+
+              {
+                size_t k = 0;
+                for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
+                  edge_list[k] = edge_id;
+                }
+
+                REQUIRE(k == edge_list.size());
+              }
+
+              Array<double> values =
+                EdgeIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, edge_list);
+
+              double error = 0;
+              for (size_t i = 0; i < edge_list.size(); ++i) {
+                error += std::abs(int_f_per_edge[edge_list[i]] - values[i]);
+              }
+
+              REQUIRE(error == Catch::Approx(0).margin(1E-10));
+            }
+
+            SECTION("SmallArray")
+            {
+              SmallArray<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
+
+              {
+                size_t k = 0;
+                for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
+                  edge_list[k] = edge_id;
+                }
+
+                REQUIRE(k == edge_list.size());
+              }
+
+              SmallArray<double> values =
+                EdgeIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, edge_list);
+
+              double error = 0;
+              for (size_t i = 0; i < edge_list.size(); ++i) {
+                error += std::abs(int_f_per_edge[edge_list[i]] - values[i]);
+              }
+
+              REQUIRE(error == Catch::Approx(0).margin(1E-10));
+            }
+          }
+        }
+      }
+    }
+  }
+}