From be75689e462e78705da032a5f501b2960e6c4314 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Wed, 1 Dec 2021 18:58:36 +0100
Subject: [PATCH] Add CellIntegrator
This allows to compute numerical quadrature of C++ functions on all
the mesh cells or on a given set of cells.
---
src/scheme/CellIntegrator.hpp | 554 +++++++++++++++++++++++++
tests/CMakeLists.txt | 1 +
tests/test_CellIntegrator.cpp | 751 ++++++++++++++++++++++++++++++++++
3 files changed, 1306 insertions(+)
create mode 100644 src/scheme/CellIntegrator.hpp
create mode 100644 tests/test_CellIntegrator.cpp
diff --git a/src/scheme/CellIntegrator.hpp b/src/scheme/CellIntegrator.hpp
new file mode 100644
index 000000000..51758560a
--- /dev/null
+++ b/src/scheme/CellIntegrator.hpp
@@ -0,0 +1,554 @@
+#ifndef CELL_INTEGRATOR_HPP
+#define CELL_INTEGRATOR_HPP
+
+#include <analysis/IQuadratureDescriptor.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <geometry/CubeTransformation.hpp>
+#include <geometry/LineTransformation.hpp>
+#include <geometry/PrismTransformation.hpp>
+#include <geometry/PyramidTransformation.hpp>
+#include <geometry/SquareTransformation.hpp>
+#include <geometry/TetrahedronTransformation.hpp>
+#include <geometry/TriangleTransformation.hpp>
+#include <mesh/CellType.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemId.hpp>
+#include <utils/Array.hpp>
+
+#include <functional>
+
+class CellIntegrator
+{
+ private:
+ template <size_t Dimension>
+ class AllCells
+ {
+ private:
+ const Connectivity<Dimension>& m_connectivity;
+
+ public:
+ using index_type = CellId;
+
+ PUGS_INLINE
+ CellId
+ cellId(const CellId cell_id) const
+ {
+ return cell_id;
+ }
+
+ PUGS_INLINE
+ size_t
+ size() const
+ {
+ return m_connectivity.numberOfCells();
+ }
+
+ PUGS_INLINE
+ AllCells(const Connectivity<Dimension>& connectivity) : m_connectivity{connectivity} {}
+ };
+
+ template <template <typename T> typename ArrayT>
+ class CellList
+ {
+ private:
+ const ArrayT<CellId>& m_cell_list;
+
+ public:
+ using index_type = typename ArrayT<CellId>::index_type;
+
+ PUGS_INLINE
+ CellId
+ cellId(const index_type index) const
+ {
+ return m_cell_list[index];
+ }
+
+ PUGS_INLINE
+ size_t
+ size() const
+ {
+ return m_cell_list.size();
+ }
+
+ PUGS_INLINE
+ CellList(const ArrayT<CellId>& cell_list) : m_cell_list{cell_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& cell_list,
+ OutputArrayT& value)
+ {
+ constexpr size_t Dimension = MeshType::Dimension;
+
+ 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 cell_to_node_matrix = connectivity.cellToNodeMatrix();
+ const auto cell_type = connectivity.cellType();
+ const auto xr = mesh.xr();
+
+ auto invalid_cell = std::make_pair(false, CellId{});
+
+ const auto qf = [&] {
+ if constexpr (Dimension == 1) {
+ return QuadratureManager::instance().getLineFormula(quadrature_descriptor);
+ } else if constexpr (Dimension == 2) {
+ return QuadratureManager::instance().getSquareFormula(quadrature_descriptor);
+ } else {
+ static_assert(Dimension == 3);
+ return QuadratureManager::instance().getCubeFormula(quadrature_descriptor);
+ }
+ }();
+
+ parallel_for(cell_list.size(), [=, &invalid_cell, &qf, &cell_list](typename ListTypeT::index_type index) {
+ const CellId cell_id = cell_list.cellId(index);
+
+ const auto cell_to_node = cell_to_node_matrix[cell_id];
+
+ if constexpr (Dimension == 1) {
+ switch (cell_type[cell_id]) {
+ case CellType::Line: {
+ const LineTransformation<1> T(xr[cell_to_node[0]], xr[cell_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.jacobianDeterminant() * function(T(xi));
+ }
+ break;
+ }
+ default: {
+ invalid_cell = std::make_pair(true, cell_id);
+ break;
+ }
+ }
+ } else if constexpr (Dimension == 2) {
+ switch (cell_type[cell_id]) {
+ case CellType::Triangle: {
+ const SquareTransformation<2> T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+ xr[cell_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.jacobianDeterminant(xi) * function(T(xi));
+ }
+ break;
+ }
+ case CellType::Quadrangle: {
+ const SquareTransformation<2> T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+ xr[cell_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.jacobianDeterminant(xi) * function(T(xi));
+ }
+ break;
+ }
+ default: {
+ invalid_cell = std::make_pair(true, cell_id);
+ break;
+ }
+ }
+ } else {
+ static_assert(Dimension == 3);
+
+ switch (cell_type[cell_id]) {
+ case CellType::Tetrahedron: {
+ const CubeTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+ xr[cell_to_node[2]], //
+ xr[cell_to_node[3]], xr[cell_to_node[3]], xr[cell_to_node[3]],
+ xr[cell_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.jacobianDeterminant(xi) * function(T(xi));
+ }
+ break;
+ }
+ case CellType::Pyramid: {
+ if (cell_to_node.size() == 5) {
+ const CubeTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]], //
+ xr[cell_to_node[3]], //
+ xr[cell_to_node[4]], xr[cell_to_node[4]], xr[cell_to_node[4]],
+ xr[cell_to_node[4]]);
+
+ 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.jacobianDeterminant(xi) * function(T(xi));
+ }
+ } else {
+ throw NotImplementedError("integration on pyramid with non-quadrangular base (number of nodes " +
+ std::to_string(cell_to_node.size()) + ")");
+ }
+ break;
+ }
+ case CellType::Diamond: {
+ if (cell_to_node.size() == 5) {
+ { // top tetrahedron
+ const CubeTransformation T0(xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
+ xr[cell_to_node[3]], //
+ xr[cell_to_node[4]], xr[cell_to_node[4]], xr[cell_to_node[4]],
+ xr[cell_to_node[4]]);
+
+ 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) * T0.jacobianDeterminant(xi) * function(T0(xi));
+ }
+ }
+ { // bottom tetrahedron
+ const CubeTransformation T1(xr[cell_to_node[3]], xr[cell_to_node[3]], xr[cell_to_node[2]],
+ xr[cell_to_node[1]], //
+ xr[cell_to_node[0]], xr[cell_to_node[0]], xr[cell_to_node[0]],
+ xr[cell_to_node[0]]);
+
+ 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) * T1.jacobianDeterminant(xi) * function(T1(xi));
+ }
+ }
+ } else if (cell_to_node.size() == 6) {
+ { // top pyramid
+ const CubeTransformation T0(xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
+ xr[cell_to_node[4]], //
+ xr[cell_to_node[5]], xr[cell_to_node[5]], xr[cell_to_node[5]],
+ xr[cell_to_node[5]]);
+
+ 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) * T0.jacobianDeterminant(xi) * function(T0(xi));
+ }
+ }
+ { // bottom pyramid
+ const CubeTransformation T1(xr[cell_to_node[4]], xr[cell_to_node[3]], xr[cell_to_node[2]],
+ xr[cell_to_node[1]], //
+ xr[cell_to_node[0]], xr[cell_to_node[0]], xr[cell_to_node[0]],
+ xr[cell_to_node[0]]);
+
+ 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) * T1.jacobianDeterminant(xi) * function(T1(xi));
+ }
+ }
+ } else {
+ throw NotImplementedError("integration on diamond with non-quadrangular base (" +
+ std::to_string(cell_to_node.size()) + ")");
+ }
+ break;
+ }
+ case CellType::Prism: {
+ const CubeTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], //
+ xr[cell_to_node[2]], xr[cell_to_node[2]], //
+ xr[cell_to_node[3]], xr[cell_to_node[4]], //
+ xr[cell_to_node[5]], xr[cell_to_node[5]]);
+
+ 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.jacobianDeterminant(xi) * function(T(xi));
+ }
+ break;
+ }
+ case CellType::Hexahedron: {
+ const CubeTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
+ xr[cell_to_node[4]], xr[cell_to_node[5]], xr[cell_to_node[6]],
+ xr[cell_to_node[7]]);
+
+ 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.jacobianDeterminant(xi) * function(T(xi));
+ }
+ break;
+ }
+ default: {
+ invalid_cell = std::make_pair(true, cell_id);
+ break;
+ }
+ }
+ }
+ });
+
+ if (invalid_cell.first) {
+ std::ostringstream os;
+ os << "invalid cell type for integration: " << name(cell_type[invalid_cell.second]);
+ throw UnexpectedError(os.str());
+ }
+ }
+
+ 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& cell_list,
+ OutputArrayT& value)
+ {
+ Assert(not quadrature_descriptor.isTensorial());
+
+ constexpr size_t Dimension = MeshType::Dimension;
+
+ 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 cell_to_node_matrix = connectivity.cellToNodeMatrix();
+ const auto cell_type = connectivity.cellType();
+ const auto xr = mesh.xr();
+
+ auto invalid_cell = std::make_pair(false, CellId{});
+
+ parallel_for(cell_list.size(), [=, &invalid_cell, &quadrature_descriptor,
+ &cell_list](const typename ListTypeT::index_type& index) {
+ const CellId cell_id = cell_list.cellId(index);
+
+ const auto cell_to_node = cell_to_node_matrix[cell_id];
+
+ if constexpr (Dimension == 1) {
+ switch (cell_type[cell_id]) {
+ case CellType::Line: {
+ const LineTransformation<1> T(xr[cell_to_node[0]], xr[cell_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.jacobianDeterminant() * function(T(xi));
+ }
+ break;
+ }
+ default: {
+ invalid_cell = std::make_pair(true, cell_id);
+ break;
+ }
+ }
+ } else if constexpr (Dimension == 2) {
+ switch (cell_type[cell_id]) {
+ case CellType::Triangle: {
+ const TriangleTransformation<2> T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_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.jacobianDeterminant() * function(T(xi));
+ }
+ break;
+ }
+ case CellType::Quadrangle: {
+ const SquareTransformation<2> T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+ xr[cell_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.jacobianDeterminant(xi) * function(T(xi));
+ }
+ break;
+ }
+ default: {
+ invalid_cell = std::make_pair(true, cell_id);
+ break;
+ }
+ }
+ } else {
+ static_assert(Dimension == 3);
+
+ switch (cell_type[cell_id]) {
+ case CellType::Tetrahedron: {
+ const TetrahedronTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], //
+ xr[cell_to_node[2]], xr[cell_to_node[3]]);
+ const auto qf = QuadratureManager::instance().getTetrahedronFormula(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.jacobianDeterminant() * function(T(xi));
+ }
+ break;
+ }
+ case CellType::Pyramid: {
+ if (cell_to_node.size() == 5) {
+ const PyramidTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]], //
+ xr[cell_to_node[3]], xr[cell_to_node[4]]);
+ const auto qf = QuadratureManager::instance().getPyramidFormula(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.jacobianDeterminant(xi) * function(T(xi));
+ }
+ } else {
+ throw NotImplementedError("integration on pyramid with non-quadrangular base");
+ }
+ break;
+ }
+ case CellType::Diamond: {
+ if (cell_to_node.size() == 5) {
+ const auto qf = QuadratureManager::instance().getTetrahedronFormula(quadrature_descriptor);
+ { // top tetrahedron
+ const TetrahedronTransformation T0(xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]], //
+ xr[cell_to_node[4]]);
+
+ 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) * T0.jacobianDeterminant() * function(T0(xi));
+ }
+ }
+ { // bottom tetrahedron
+ const TetrahedronTransformation T1(xr[cell_to_node[3]], xr[cell_to_node[2]], xr[cell_to_node[1]], //
+ xr[cell_to_node[0]]);
+
+ 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) * T1.jacobianDeterminant() * function(T1(xi));
+ }
+ }
+ } else if (cell_to_node.size() == 6) {
+ const auto qf = QuadratureManager::instance().getCubeFormula(quadrature_descriptor);
+ { // top pyramid
+ const CubeTransformation T0(xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
+ xr[cell_to_node[4]], //
+ xr[cell_to_node[5]], xr[cell_to_node[5]], xr[cell_to_node[5]],
+ xr[cell_to_node[5]]);
+
+ 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) * T0.jacobianDeterminant(xi) * function(T0(xi));
+ }
+ }
+ { // bottom pyramid
+ const CubeTransformation T1(xr[cell_to_node[4]], xr[cell_to_node[3]], xr[cell_to_node[2]],
+ xr[cell_to_node[1]], //
+ xr[cell_to_node[0]], xr[cell_to_node[0]], xr[cell_to_node[0]],
+ xr[cell_to_node[0]]);
+
+ 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) * T1.jacobianDeterminant(xi) * function(T1(xi));
+ }
+ }
+ } else {
+ throw NotImplementedError("integration on diamond with non-quadrangular base (" +
+ std::to_string(cell_to_node.size()) + ")");
+ }
+ break;
+ }
+ case CellType::Prism: {
+ const PrismTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+ xr[cell_to_node[3]], xr[cell_to_node[4]], xr[cell_to_node[5]]);
+ const auto qf = QuadratureManager::instance().getPrismFormula(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.jacobianDeterminant(xi) * function(T(xi));
+ }
+ break;
+ }
+ case CellType::Hexahedron: {
+ const CubeTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
+ xr[cell_to_node[4]], xr[cell_to_node[5]], xr[cell_to_node[6]],
+ xr[cell_to_node[7]]);
+ const auto qf = QuadratureManager::instance().getCubeFormula(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.jacobianDeterminant(xi) * function(T(xi));
+ }
+ break;
+ }
+ default: {
+ invalid_cell = std::make_pair(true, cell_id);
+ break;
+ }
+ }
+ }
+ });
+
+ if (invalid_cell.first) {
+ std::ostringstream os;
+ os << "invalid cell type for integration: " << name(cell_type[invalid_cell.second]);
+ throw UnexpectedError(os.str());
+ }
+ }
+
+ 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,
+ AllCells<Dimension>{mesh.connectivity()}, value);
+ } else {
+ _directIntegrateTo<MeshType, OutputArrayT>(function, quadrature_descriptor, mesh,
+ AllCells<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<CellId>& cell_list)
+ {
+ ArrayT<OutputType> value(size(cell_list));
+ if (quadrature_descriptor.isTensorial()) {
+ _tensorialIntegrateTo<MeshType, ArrayT<OutputType>>(function, quadrature_descriptor, mesh, CellList{cell_list},
+ value);
+ } else {
+ _directIntegrateTo<MeshType, ArrayT<OutputType>>(function, quadrature_descriptor, mesh, CellList{cell_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<CellId>& cell_list)
+ {
+ return integrate(std::function(function), quadrature_descriptor, mesh, cell_list);
+ }
+};
+
+#endif // CELL_INTEGRATOR_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 0994985c3..7f4c6a323 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -127,6 +127,7 @@ add_executable (unit_tests
add_executable (mpi_unit_tests
mpi_test_main.cpp
+ test_CellIntegrator.cpp
test_DiscreteFunctionInterpoler.cpp
test_DiscreteFunctionP0.cpp
test_DiscreteFunctionP0Vector.cpp
diff --git a/tests/test_CellIntegrator.cpp b/tests/test_CellIntegrator.cpp
new file mode 100644
index 000000000..ef0fca0b9
--- /dev/null
+++ b/tests/test_CellIntegrator.cpp
@@ -0,0 +1,751 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <analysis/GaussLobattoQuadratureDescriptor.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <mesh/ItemValue.hpp>
+#include <mesh/Mesh.hpp>
+#include <scheme/CellIntegrator.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("CellIntegrator", "[scheme]")
+{
+ SECTION("1D")
+ {
+ using R1 = TinyVector<1>;
+
+ const auto mesh = MeshDataBaseForTests::get().unordered1DMesh();
+ auto f = [](const R1& x) -> double { return x[0] * x[0] + 1; };
+
+ Array<const double> int_f_per_cell = [=] {
+ Array<double> int_f(mesh->numberOfCells());
+ auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto cell_node_list = cell_to_node_matrix[cell_id];
+ auto xr = mesh->xr();
+ const double x_left = xr[cell_node_list[0]][0];
+ const double x_right = xr[cell_node_list[1]][0];
+ int_f[cell_id] = 1. / 3 * (x_right * x_right * x_right - x_left * x_left * x_left) + x_right - x_left;
+ });
+
+ return int_f;
+ }();
+
+ SECTION("direct formula")
+ {
+ SECTION("all cells")
+ {
+ SECTION("CellValue")
+ {
+ CellValue<double> values(mesh->connectivity());
+ CellIntegrator::integrateTo([=](const R1 x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("Array")
+ {
+ Array<double> values(mesh->numberOfCells());
+
+ CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<double> values(mesh->numberOfCells());
+ CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+
+ SECTION("cell list")
+ {
+ SECTION("Array")
+ {
+ Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
+
+ {
+ size_t k = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
+ cell_list[k] = cell_id;
+ }
+
+ REQUIRE(k == cell_list.size());
+ }
+
+ Array<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list);
+
+ double error = 0;
+ for (size_t i = 0; i < cell_list.size(); ++i) {
+ error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
+
+ {
+ size_t k = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
+ cell_list[k] = cell_id;
+ }
+
+ REQUIRE(k == cell_list.size());
+ }
+
+ SmallArray<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list);
+
+ double error = 0;
+ for (size_t i = 0; i < cell_list.size(); ++i) {
+ error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+ }
+
+ SECTION("tensorial formula")
+ {
+ SECTION("all cells")
+ {
+ SECTION("CellValue")
+ {
+ CellValue<double> values(mesh->connectivity());
+ CellIntegrator::integrateTo([=](const R1 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh,
+ values);
+
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("Array")
+ {
+ Array<double> values(mesh->numberOfCells());
+
+ CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<double> values(mesh->numberOfCells());
+ CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+
+ SECTION("cell list")
+ {
+ SECTION("Array")
+ {
+ Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
+
+ {
+ size_t k = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
+ cell_list[k] = cell_id;
+ }
+
+ REQUIRE(k == cell_list.size());
+ }
+
+ Array<double> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list);
+
+ double error = 0;
+ for (size_t i = 0; i < cell_list.size(); ++i) {
+ error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
+
+ {
+ size_t k = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
+ cell_list[k] = cell_id;
+ }
+
+ REQUIRE(k == cell_list.size());
+ }
+
+ SmallArray<double> values =
+ CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list);
+
+ double error = 0;
+ for (size_t i = 0; i < cell_list.size(); ++i) {
+ error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+ }
+ }
+
+ 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_cell = [=] {
+ Array<double> int_f(mesh->numberOfCells());
+ auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+ auto cell_type = mesh->connectivity().cellType();
+
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto cell_node_list = cell_to_node_matrix[cell_id];
+ auto xr = mesh->xr();
+ double integral = 0;
+
+ switch (cell_type[cell_id]) {
+ case CellType::Triangle: {
+ TriangleTransformation<2> T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]]);
+ auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(4));
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.jacobianDeterminant() * f(T(xi));
+ }
+ break;
+ }
+ case CellType::Quadrangle: {
+ SquareTransformation<2> T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]]);
+ auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(4));
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi));
+ }
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid cell type in 2d");
+ }
+ }
+ int_f[cell_id] = integral;
+ });
+
+ return int_f;
+ }();
+
+ SECTION("direct formula")
+ {
+ SECTION("all cells")
+ {
+ SECTION("CellValue")
+ {
+ CellValue<double> values(mesh->connectivity());
+ CellIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("Array")
+ {
+ Array<double> values(mesh->numberOfCells());
+
+ CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<double> values(mesh->numberOfCells());
+ CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+
+ SECTION("cell list")
+ {
+ SECTION("Array")
+ {
+ Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
+
+ {
+ size_t k = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
+ cell_list[k] = cell_id;
+ }
+
+ REQUIRE(k == cell_list.size());
+ }
+
+ Array<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list);
+
+ double error = 0;
+ for (size_t i = 0; i < cell_list.size(); ++i) {
+ error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
+
+ {
+ size_t k = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
+ cell_list[k] = cell_id;
+ }
+
+ REQUIRE(k == cell_list.size());
+ }
+
+ SmallArray<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list);
+
+ double error = 0;
+ for (size_t i = 0; i < cell_list.size(); ++i) {
+ error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+ }
+
+ SECTION("tensorial formula")
+ {
+ SECTION("all cells")
+ {
+ SECTION("CellValue")
+ {
+ CellValue<double> values(mesh->connectivity());
+ CellIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh,
+ values);
+
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("Array")
+ {
+ Array<double> values(mesh->numberOfCells());
+
+ CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<double> values(mesh->numberOfCells());
+ CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+
+ SECTION("cell list")
+ {
+ SECTION("Array")
+ {
+ Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
+
+ {
+ size_t k = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
+ cell_list[k] = cell_id;
+ }
+
+ REQUIRE(k == cell_list.size());
+ }
+
+ Array<double> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list);
+
+ double error = 0;
+ for (size_t i = 0; i < cell_list.size(); ++i) {
+ error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
+
+ {
+ size_t k = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
+ cell_list[k] = cell_id;
+ }
+
+ REQUIRE(k == cell_list.size());
+ }
+
+ SmallArray<double> values =
+ CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list);
+
+ double error = 0;
+ for (size_t i = 0; i < cell_list.size(); ++i) {
+ error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+ }
+ }
+
+ SECTION("3D")
+ {
+ using R3 = TinyVector<3>;
+
+ const auto 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;
+ };
+
+ Array<const double> int_f_per_cell = [=] {
+ Array<double> int_f(mesh->numberOfCells());
+ auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+ auto cell_type = mesh->connectivity().cellType();
+
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto cell_node_list = cell_to_node_matrix[cell_id];
+ auto xr = mesh->xr();
+ double integral = 0;
+
+ switch (cell_type[cell_id]) {
+ case CellType::Tetrahedron: {
+ TetrahedronTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]]);
+ auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor(4));
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.jacobianDeterminant() * f(T(xi));
+ }
+ break;
+ }
+ case CellType::Pyramid: {
+ PyramidTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]], xr[cell_node_list[4]]);
+ auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor(4));
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi));
+ }
+ break;
+ }
+ case CellType::Prism: {
+ PrismTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[5]]);
+ auto qf = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor(4));
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi));
+ }
+ break;
+ }
+ case CellType::Hexahedron: {
+ CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[5]],
+ xr[cell_node_list[6]], xr[cell_node_list[7]]);
+ auto qf = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor(4));
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi));
+ }
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid cell type in 2d");
+ }
+ }
+ int_f[cell_id] = integral;
+ });
+
+ return int_f;
+ }();
+
+ SECTION("direct formula")
+ {
+ SECTION("all cells")
+ {
+ SECTION("CellValue")
+ {
+ CellValue<double> values(mesh->connectivity());
+ CellIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("Array")
+ {
+ Array<double> values(mesh->numberOfCells());
+
+ CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<double> values(mesh->numberOfCells());
+ CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+
+ SECTION("cell list")
+ {
+ SECTION("Array")
+ {
+ Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
+
+ {
+ size_t k = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
+ cell_list[k] = cell_id;
+ }
+
+ REQUIRE(k == cell_list.size());
+ }
+
+ Array<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, cell_list);
+
+ double error = 0;
+ for (size_t i = 0; i < cell_list.size(); ++i) {
+ error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
+
+ {
+ size_t k = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
+ cell_list[k] = cell_id;
+ }
+
+ REQUIRE(k == cell_list.size());
+ }
+
+ SmallArray<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, cell_list);
+
+ double error = 0;
+ for (size_t i = 0; i < cell_list.size(); ++i) {
+ error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+ }
+
+ SECTION("tensorial formula")
+ {
+ SECTION("all cells")
+ {
+ SECTION("CellValue")
+ {
+ CellValue<double> values(mesh->connectivity());
+ CellIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLobattoQuadratureDescriptor(4), *mesh,
+ values);
+
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("Array")
+ {
+ Array<double> values(mesh->numberOfCells());
+
+ CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<double> values(mesh->numberOfCells());
+ CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+
+ SECTION("cell list")
+ {
+ SECTION("Array")
+ {
+ Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
+
+ {
+ size_t k = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
+ cell_list[k] = cell_id;
+ }
+
+ REQUIRE(k == cell_list.size());
+ }
+
+ Array<double> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, cell_list);
+
+ double error = 0;
+ for (size_t i = 0; i < cell_list.size(); ++i) {
+ error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
+
+ {
+ size_t k = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
+ cell_list[k] = cell_id;
+ }
+
+ REQUIRE(k == cell_list.size());
+ }
+
+ SmallArray<double> values =
+ CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, cell_list);
+
+ double error = 0;
+ for (size_t i = 0; i < cell_list.size(); ++i) {
+ error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+ }
+ }
+}
--
GitLab