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

+#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

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

tests/test_CellIntegrator.cpp

+#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));
+        }
+      }
+    }
+  }
+}