From 91ae89cb0f9e97dbbc1d767a6171fbdac754898f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Mon, 6 Dec 2021 18:28:50 +0100
Subject: [PATCH] Add EdgeIntegrator (and associated tests)
---
src/scheme/EdgeIntegrator.hpp | 166 +++++++++++++++++++++
tests/CMakeLists.txt | 1 +
tests/test_EdgeIntegrator.cpp | 261 ++++++++++++++++++++++++++++++++++
3 files changed, 428 insertions(+)
create mode 100644 src/scheme/EdgeIntegrator.hpp
create mode 100644 tests/test_EdgeIntegrator.cpp
diff --git a/src/scheme/EdgeIntegrator.hpp b/src/scheme/EdgeIntegrator.hpp
new file mode 100644
index 000000000..2bf287c3d
--- /dev/null
+++ b/src/scheme/EdgeIntegrator.hpp
@@ -0,0 +1,166 @@
+#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
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index b10741748..f825a2cfd 100644
--- a/tests/CMakeLists.txt
+++ b/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
diff --git a/tests/test_EdgeIntegrator.cpp b/tests/test_EdgeIntegrator.cpp
new file mode 100644
index 000000000..0cd2feb9d
--- /dev/null
+++ b/tests/test_EdgeIntegrator.cpp
@@ -0,0 +1,261 @@
+#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));
+ }
+ }
+ }
+ }
+ }
+ }
+}
--
GitLab