diff --git a/src/scheme/FaceIntegrator.hpp b/src/scheme/FaceIntegrator.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..30ee39077ea12d898b4bdee00f12970be79e8753
--- /dev/null
+++ b/src/scheme/FaceIntegrator.hpp
@@ -0,0 +1,322 @@
+#ifndef FACE_INTEGRATOR_HPP
+#define FACE_INTEGRATOR_HPP
+
+#include <analysis/IQuadratureDescriptor.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <geometry/LineTransformation.hpp>
+#include <geometry/SquareTransformation.hpp>
+#include <geometry/TriangleTransformation.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemId.hpp>
+#include <utils/Array.hpp>
+
+#include <functional>
+
+class FaceIntegrator
+{
+ private:
+ template <size_t Dimension>
+ class AllFaces
+ {
+ private:
+ const Connectivity<Dimension>& m_connectivity;
+
+ public:
+ using index_type = FaceId;
+
+ PUGS_INLINE
+ FaceId
+ faceId(const FaceId face_id) const
+ {
+ return face_id;
+ }
+
+ PUGS_INLINE
+ size_t
+ size() const
+ {
+ return m_connectivity.numberOfFaces();
+ }
+
+ PUGS_INLINE
+ AllFaces(const Connectivity<Dimension>& connectivity) : m_connectivity{connectivity} {}
+ };
+
+ template <template <typename T> typename ArrayT>
+ class FaceList
+ {
+ private:
+ const ArrayT<FaceId>& m_face_list;
+
+ public:
+ using index_type = typename ArrayT<FaceId>::index_type;
+
+ PUGS_INLINE
+ FaceId
+ faceId(const index_type index) const
+ {
+ return m_face_list[index];
+ }
+
+ PUGS_INLINE
+ size_t
+ size() const
+ {
+ return m_face_list.size();
+ }
+
+ PUGS_INLINE
+ FaceList(const ArrayT<FaceId>& face_list) : m_face_list{face_list} {}
+ };
+
+ template <typename MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
+ static PUGS_INLINE void
+ _tensorialIntegrateTo(std::function<OutputType(InputType)> function,
+ const IQuadratureDescriptor& quadrature_descriptor,
+ const MeshType& mesh,
+ const ListTypeT& face_list,
+ OutputArrayT& value)
+ {
+ constexpr size_t Dimension = MeshType::Dimension;
+ static_assert(Dimension > 1);
+
+ 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 face_to_node_matrix = connectivity.faceToNodeMatrix();
+ const auto xr = mesh.xr();
+
+ if constexpr (Dimension == 2) {
+ const auto qf = QuadratureManager::instance().getLineFormula(quadrature_descriptor);
+
+ parallel_for(face_list.size(), [=, &qf, &face_list](typename ListTypeT::index_type index) {
+ const FaceId face_id = face_list.faceId(index);
+
+ const auto face_to_node = face_to_node_matrix[face_id];
+
+ Assert(face_to_node.size() == 2);
+ const LineTransformation<2> T(xr[face_to_node[0]], xr[face_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));
+ }
+ });
+
+ } else if constexpr (Dimension == 3) {
+ const auto qf = QuadratureManager::instance().getSquareFormula(quadrature_descriptor);
+
+ auto invalid_face = std::make_pair(false, FaceId{});
+
+ parallel_for(face_list.size(), [=, &invalid_face, &qf, &face_list](typename ListTypeT::index_type index) {
+ const FaceId face_id = face_list.faceId(index);
+
+ const auto face_to_node = face_to_node_matrix[face_id];
+
+ switch (face_to_node.size()) {
+ case 3: {
+ SquareTransformation<3> T(xr[face_to_node[0]], xr[face_to_node[1]], xr[face_to_node[2]], xr[face_to_node[2]]);
+ 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.areaVariationNorm(xi) * function(T(xi));
+ }
+ break;
+ }
+ case 4: {
+ SquareTransformation<3> T(xr[face_to_node[0]], xr[face_to_node[1]], xr[face_to_node[2]], xr[face_to_node[3]]);
+ 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.areaVariationNorm(xi) * function(T(xi));
+ }
+ break;
+ }
+ default: {
+ invalid_face = std::make_pair(true, face_id);
+ }
+ }
+ });
+
+ // LCOV_EXCL_START
+ if (invalid_face.first) {
+ std::ostringstream os;
+ os << "invalid face type for integration: " << face_to_node_matrix[invalid_face.second].size() << " points";
+ throw UnexpectedError(os.str());
+ }
+ // LCOV_EXCL_STOP
+
+ } else {
+ throw UnexpectedError("invalid dimension for face integration");
+ }
+ }
+
+ template <typename MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
+ static PUGS_INLINE void
+ _directIntegrateTo(std::function<OutputType(InputType)> function,
+ const IQuadratureDescriptor& quadrature_descriptor,
+ const MeshType& mesh,
+ const ListTypeT& face_list,
+ OutputArrayT& value)
+ {
+ Assert(not quadrature_descriptor.isTensorial());
+
+ constexpr size_t Dimension = MeshType::Dimension;
+ static_assert(Dimension > 1);
+
+ 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");
+
+ 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 face_to_node_matrix = connectivity.faceToNodeMatrix();
+ const auto xr = mesh.xr();
+
+ if constexpr (Dimension == 2) {
+ parallel_for(face_list.size(),
+ [=, &quadrature_descriptor, &face_list](const typename ListTypeT::index_type& index) {
+ const FaceId face_id = face_list.faceId(index);
+
+ const auto face_to_node = face_to_node_matrix[face_id];
+
+ Assert(face_to_node.size() == 2);
+ const LineTransformation<2> T(xr[face_to_node[0]], xr[face_to_node[1]]);
+ const auto qf = QuadratureManager::instance().getLineFormula(quadrature_descriptor);
+
+ 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));
+ }
+ });
+
+ } else if constexpr (Dimension == 3) {
+ auto invalid_face = std::make_pair(false, FaceId{});
+
+ parallel_for(face_list.size(),
+ [=, &invalid_face, &quadrature_descriptor, &face_list](const typename ListTypeT::index_type& index) {
+ const FaceId face_id = face_list.faceId(index);
+
+ const auto face_to_node = face_to_node_matrix[face_id];
+
+ switch (face_to_node.size()) {
+ case 3: {
+ const TriangleTransformation<3> T(xr[face_to_node[0]], xr[face_to_node[1]], xr[face_to_node[2]]);
+ const auto qf = QuadratureManager::instance().getTriangleFormula(quadrature_descriptor);
+
+ 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.areaVariationNorm() * function(T(xi));
+ }
+ break;
+ }
+ case 4: {
+ const SquareTransformation<3> T(xr[face_to_node[0]], xr[face_to_node[1]], xr[face_to_node[2]],
+ xr[face_to_node[3]]);
+ const auto qf = QuadratureManager::instance().getSquareFormula(quadrature_descriptor);
+
+ 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.areaVariationNorm(xi) * function(T(xi));
+ }
+ break;
+ }
+ // LCOV_EXCL_START
+ default: {
+ invalid_face = std::make_pair(true, face_id);
+ break;
+ }
+ // LCOV_EXCL_STOP
+ }
+ });
+
+ // LCOV_EXCL_START
+ if (invalid_face.first) {
+ std::ostringstream os;
+ os << "invalid face type for integration: " << face_to_node_matrix[invalid_face.second].size() << " points";
+ throw UnexpectedError(os.str());
+ }
+ // LCOV_EXCL_STOP
+ }
+ }
+
+ 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;
+
+ if (quadrature_descriptor.isTensorial()) {
+ _tensorialIntegrateTo<MeshType, OutputArrayT>(function, quadrature_descriptor, mesh,
+ AllFaces<Dimension>{mesh.connectivity()}, value);
+ } else {
+ _directIntegrateTo<MeshType, OutputArrayT>(function, quadrature_descriptor, mesh,
+ AllFaces<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<FaceId>& face_list)
+ {
+ ArrayT<OutputType> value(size(face_list));
+ if (quadrature_descriptor.isTensorial()) {
+ _tensorialIntegrateTo<MeshType, ArrayT<OutputType>>(function, quadrature_descriptor, mesh, FaceList{face_list},
+ value);
+ } else {
+ _directIntegrateTo<MeshType, ArrayT<OutputType>>(function, quadrature_descriptor, mesh, FaceList{face_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<FaceId>& face_list)
+ {
+ return integrate(std::function(function), quadrature_descriptor, mesh, face_list);
+ }
+};
+
+#endif // FACE_INTEGRATOR_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 8d7595e7ca2e9809829f9636e6d52de4c343910b..b10741748ad167b6d0f806ece0571cdbdb00eeab 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -76,6 +76,7 @@ add_executable (unit_tests
test_EscapedString.cpp
test_Exceptions.cpp
test_ExecutionPolicy.cpp
+ test_FaceIntegrator.cpp
test_FakeProcessor.cpp
test_ForProcessor.cpp
test_FunctionArgumentConverter.cpp
diff --git a/tests/test_FaceIntegrator.cpp b/tests/test_FaceIntegrator.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9dee9c032b48db4ad51c362b5f47681613039fb5
--- /dev/null
+++ b/tests/test_FaceIntegrator.cpp
@@ -0,0 +1,506 @@
+#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/FaceIntegrator.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("FaceIntegrator", "[scheme]")
+{
+ SECTION("2D")
+ {
+ using R2 = TinyVector<2>;
+
+ const auto mesh = MeshDataBaseForTests::get().hybrid2DMesh();
+
+ auto f = [](const R2& X) -> double {
+ const double x = X[0];
+ const double y = X[1];
+ return x * x + 2 * x * y + 3 * y * y + 2;
+ };
+
+ Array<const double> int_f_per_face = [=] {
+ Array<double> int_f(mesh->numberOfFaces());
+ auto face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+
+ parallel_for(
+ mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
+ auto face_node_list = face_to_node_matrix[face_id];
+ auto xr = mesh->xr();
+ double integral = 0;
+ LineTransformation<2> T(xr[face_node_list[0]], xr[face_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[face_id] = integral;
+ });
+
+ return int_f;
+ }();
+
+ SECTION("direct formula")
+ {
+ SECTION("all faces")
+ {
+ SECTION("FaceValue")
+ {
+ FaceValue<double> values(mesh->connectivity());
+ FaceIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("Array")
+ {
+ Array<double> values(mesh->numberOfFaces());
+
+ FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<double> values(mesh->numberOfFaces());
+ FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+
+ SECTION("face list")
+ {
+ SECTION("Array")
+ {
+ Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ Array<double> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ SmallArray<double> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+ }
+
+ SECTION("tensorial formula")
+ {
+ SECTION("all faces")
+ {
+ SECTION("FaceValue")
+ {
+ FaceValue<double> values(mesh->connectivity());
+ FaceIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh,
+ values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("Array")
+ {
+ Array<double> values(mesh->numberOfFaces());
+
+ FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<double> values(mesh->numberOfFaces());
+ FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+
+ SECTION("face list")
+ {
+ SECTION("Array")
+ {
+ Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ Array<double> values = FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ SmallArray<double> values =
+ FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+ }
+ }
+
+ 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_face = [=] {
+ Array<double> int_f(mesh->numberOfFaces());
+ auto face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+
+ parallel_for(
+ mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
+ auto face_node_list = face_to_node_matrix[face_id];
+ auto xr = mesh->xr();
+ double integral = 0;
+
+ switch (face_node_list.size()) {
+ case 3: {
+ TriangleTransformation<3> T(xr[face_node_list[0]], xr[face_node_list[1]], xr[face_node_list[2]]);
+ auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(2));
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.areaVariationNorm() * f(T(xi));
+ }
+ break;
+ }
+ case 4: {
+ SquareTransformation<3> T(xr[face_node_list[0]], xr[face_node_list[1]], xr[face_node_list[2]],
+ xr[face_node_list[3]]);
+ auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(2));
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.areaVariationNorm(xi) * f(T(xi));
+ }
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid face (node number must be 3 or 4)");
+ }
+ }
+ int_f[face_id] = integral;
+ });
+
+ return int_f;
+ }();
+
+ SECTION(mesh_name)
+ {
+ SECTION("direct formula")
+ {
+ SECTION("all faces")
+ {
+ SECTION("FaceValue")
+ {
+ FaceValue<double> values(mesh->connectivity());
+ FaceIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh,
+ values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("Array")
+ {
+ Array<double> values(mesh->numberOfFaces());
+
+ FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<double> values(mesh->numberOfFaces());
+ FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+
+ SECTION("face list")
+ {
+ SECTION("Array")
+ {
+ Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ Array<double> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ SmallArray<double> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+ }
+
+ SECTION("tensorial formula")
+ {
+ SECTION("all faces")
+ {
+ SECTION("FaceValue")
+ {
+ FaceValue<double> values(mesh->connectivity());
+ FaceIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10),
+ *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("Array")
+ {
+ Array<double> values(mesh->numberOfFaces());
+
+ FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<double> values(mesh->numberOfFaces());
+ FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+
+ SECTION("face list")
+ {
+ SECTION("Array")
+ {
+ Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ Array<double> values =
+ FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ SmallArray<double> values =
+ FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+ }
+ }
+ }
+ }
+}