From 483ec22967b06f6a7f34d2764de6bd144dcee5ee Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Fri, 4 Feb 2022 18:56:20 +0100
Subject: [PATCH] Add missing tests for CellIntegrator

---
 tests/test_CellIntegrator.cpp | 3211 ++++++++++++++++++++++++++-------
 1 file changed, 2531 insertions(+), 680 deletions(-)

diff --git a/tests/test_CellIntegrator.cpp b/tests/test_CellIntegrator.cpp
index 49f750afb..ebd322d14 100644
--- a/tests/test_CellIntegrator.cpp
+++ b/tests/test_CellIntegrator.cpp
@@ -15,915 +15,2766 @@
 
 TEST_CASE("CellIntegrator", "[scheme]")
 {
-  SECTION("1D")
+  SECTION("scalar")
   {
-    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("1D")
     {
-      SECTION("all cells")
+      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("CellValue")
+        SECTION("all cells")
         {
-          CellValue<double> values(mesh->connectivity());
-          CellIntegrator::integrateTo([=](const R1 x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values);
+          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>;
+
+      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", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh)));
+
+      for (auto mesh_info : mesh_list) {
+        auto mesh_name = mesh_info.first;
+        auto mesh      = mesh_info.second;
+
+        SECTION(mesh_name)
+        {
+          SECTION("direct formula")
+          {
+            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;
+                  }
+                  case CellType::Diamond: {
+                    if (cell_node_list.size() == 5) {
+                      auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor(4));
+                      {   // top tetrahedron
+                        TetrahedronTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]],
+                                                     xr[cell_node_list[3]], xr[cell_node_list[4]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T0.jacobianDeterminant() * f(T0(xi));
+                        }
+                      }
+                      {   // bottom tetrahedron
+                        TetrahedronTransformation T1(xr[cell_node_list[3]], xr[cell_node_list[2]],
+                                                     xr[cell_node_list[1]], xr[cell_node_list[0]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T1.jacobianDeterminant() * f(T1(xi));
+                        }
+                      }
+                    } else if (cell_node_list.size() == 6) {
+                      auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor(4));
+                      {   // top pyramid
+                        PyramidTransformation T0(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]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi));
+                        }
+                      }
+                      {   // bottom pyramid
+                        PyramidTransformation T1(xr[cell_node_list[4]], xr[cell_node_list[3]], xr[cell_node_list[2]],
+                                                 xr[cell_node_list[1]], xr[cell_node_list[0]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi));
+                        }
+                      }
+                    } else {
+                      INFO("Diamond cells with more than 6 vertices are not tested");
+                      REQUIRE(false);
+                    }
+                    break;
+                  }
+                  default: {
+                    INFO("Diamond cells not tested yet");
+                    REQUIRE(cell_type[cell_id] != CellType::Diamond);
+                  }
+                  }
+                  int_f[cell_id] = integral;
+                });
+
+              return int_f;
+            }();
+
+            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 == 0);
+              }
+
+              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 == 0);
+              }
+
+              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 == 0);
+              }
+            }
+
+            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 == 0);
+              }
+
+              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 == 0);
+              }
+            }
+          }
+
+          SECTION("tensorial formula")
+          {
+            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();
+
+              auto qf = QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(4));
+
+              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: {
+                    CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+                                         xr[cell_node_list[2]], xr[cell_node_list[3]], xr[cell_node_list[3]],
+                                         xr[cell_node_list[3]], xr[cell_node_list[3]]);
+
+                    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::Pyramid: {
+                    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[4]],
+                                         xr[cell_node_list[4]], xr[cell_node_list[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: {
+                    CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+                                         xr[cell_node_list[2]], xr[cell_node_list[3]], xr[cell_node_list[4]],
+                                         xr[cell_node_list[5]], xr[cell_node_list[5]]);
+                    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]]);
+
+                    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::Diamond: {
+                    if (cell_node_list.size() == 5) {
+                      {   // top tetrahedron
+                        CubeTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]], xr[cell_node_list[3]],
+                                              xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[4]],
+                                              xr[cell_node_list[4]], xr[cell_node_list[4]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi));
+                        }
+                      }
+                      {   // bottom tetrahedron
+                        CubeTransformation T1(xr[cell_node_list[3]], xr[cell_node_list[2]], xr[cell_node_list[1]],
+                                              xr[cell_node_list[1]], xr[cell_node_list[0]], xr[cell_node_list[0]],
+                                              xr[cell_node_list[0]], xr[cell_node_list[0]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi));
+                        }
+                      }
+                    } else if (cell_node_list.size() == 6) {
+                      {   // top pyramid
+                        CubeTransformation T0(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[5]],
+                                              xr[cell_node_list[5]], xr[cell_node_list[5]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi));
+                        }
+                      }
+                      {   // bottom pyramid
+                        CubeTransformation T1(xr[cell_node_list[4]], xr[cell_node_list[3]], xr[cell_node_list[2]],
+                                              xr[cell_node_list[1]], xr[cell_node_list[0]], xr[cell_node_list[0]],
+                                              xr[cell_node_list[0]], xr[cell_node_list[0]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi));
+                        }
+                      }
+                    } else {
+                      INFO("Diamond cells with more than 6 vertices are not tested");
+                      REQUIRE(false);
+                    }
+                    break;
+                  }
+                  default: {
+                    INFO("Diamond cells not tested yet");
+                    REQUIRE(cell_type[cell_id] != CellType::Diamond);
+                  }
+                  }
+                  int_f[cell_id] = integral;
+                });
+
+              return int_f;
+            }();
+
+            SECTION("all cells")
+            {
+              SECTION("CellValue")
+              {
+                CellValue<double> values(mesh->connectivity());
+                CellIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10),
+                                            *mesh, values);
+
+                auto cell_type = mesh->connectivity().cellType();
+                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));
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+
+  SECTION("vector")
+  {
+    using R2 = TinyVector<2>;
+
+    SECTION("1D")
+    {
+      using R1 = TinyVector<1>;
+
+      const auto mesh = MeshDataBaseForTests::get().unordered1DMesh();
+      auto f          = [](const R1& x) -> R2 { return R2{x[0] * x[0] + 1, 2 * x[0]}; };
+
+      Array<const R2> int_f_per_cell = [=] {
+        Array<R2> 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] = R2{1. / 3 * (x_right * x_right * x_right - x_left * x_left * x_left) + x_right - x_left,
+                                x_right * x_right - x_left * x_left};
+          });
+
+        return int_f;
+      }();
+
+      SECTION("direct formula")
+      {
+        SECTION("all cells")
+        {
+          SECTION("CellValue")
+          {
+            CellValue<R2> 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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
+
+            REQUIRE(error == Catch::Approx(0).margin(1E-10));
+          }
+
+          SECTION("Array")
+          {
+            Array<R2> 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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
+
+            REQUIRE(error == Catch::Approx(0).margin(1E-10));
+          }
+
+          SECTION("SmallArray")
+          {
+            SmallArray<R2> 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 += l2Norm(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<R2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list);
+
+            double error = 0;
+            for (size_t i = 0; i < cell_list.size(); ++i) {
+              error += l2Norm(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<R2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list);
+
+            double error = 0;
+            for (size_t i = 0; i < cell_list.size(); ++i) {
+              error += l2Norm(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<R2> 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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
+
+            REQUIRE(error == Catch::Approx(0).margin(1E-10));
+          }
+
+          SECTION("Array")
+          {
+            Array<R2> 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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
+
+            REQUIRE(error == Catch::Approx(0).margin(1E-10));
+          }
+
+          SECTION("SmallArray")
+          {
+            SmallArray<R2> 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 += l2Norm(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<R2> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list);
+
+            double error = 0;
+            for (size_t i = 0; i < cell_list.size(); ++i) {
+              error += l2Norm(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<R2> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list);
+
+            double error = 0;
+            for (size_t i = 0; i < cell_list.size(); ++i) {
+              error += l2Norm(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) -> R2 {
+        const double x = X[0];
+        const double y = X[1];
+        return R2{x * x + 2 * x * y + 3 * y * y + 2, 2 * x + 3 * y * y};
+      };
+
+      Array<const R2> int_f_per_cell = [=] {
+        Array<R2> 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();
+            R2 integral         = zero;
+
+            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<R2> 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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
+
+            REQUIRE(error == Catch::Approx(0).margin(1E-10));
+          }
+
+          SECTION("Array")
+          {
+            Array<R2> 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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
+
+            REQUIRE(error == Catch::Approx(0).margin(1E-10));
+          }
+
+          SECTION("SmallArray")
+          {
+            SmallArray<R2> 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 += l2Norm(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<R2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list);
+
+            double error = 0;
+            for (size_t i = 0; i < cell_list.size(); ++i) {
+              error += l2Norm(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<R2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list);
+
+            double error = 0;
+            for (size_t i = 0; i < cell_list.size(); ++i) {
+              error += l2Norm(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<R2> 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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
+
+            REQUIRE(error == Catch::Approx(0).margin(1E-10));
+          }
+
+          SECTION("Array")
+          {
+            Array<R2> 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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
+
+            REQUIRE(error == Catch::Approx(0).margin(1E-10));
+          }
+
+          SECTION("SmallArray")
+          {
+            SmallArray<R2> 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 += l2Norm(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<R2> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list);
+
+            double error = 0;
+            for (size_t i = 0; i < cell_list.size(); ++i) {
+              error += l2Norm(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<R2> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list);
+
+            double error = 0;
+            for (size_t i = 0; i < cell_list.size(); ++i) {
+              error += l2Norm(int_f_per_cell[cell_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) -> R2 {
+        const double x = X[0];
+        const double y = X[1];
+        const double z = X[2];
+        return R2{x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1, 2 * x * y - 3 * y * z + 2 * x};
+      };
+
+      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", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh)));
+
+      for (auto mesh_info : mesh_list) {
+        auto mesh_name = mesh_info.first;
+        auto mesh      = mesh_info.second;
+
+        SECTION(mesh_name)
+        {
+          SECTION("direct formula")
+          {
+            Array<const R2> int_f_per_cell = [=] {
+              Array<R2> 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();
+                  R2 integral         = zero;
+
+                  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;
+                  }
+                  case CellType::Diamond: {
+                    if (cell_node_list.size() == 5) {
+                      auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor(4));
+                      {   // top tetrahedron
+                        TetrahedronTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]],
+                                                     xr[cell_node_list[3]], xr[cell_node_list[4]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T0.jacobianDeterminant() * f(T0(xi));
+                        }
+                      }
+                      {   // bottom tetrahedron
+                        TetrahedronTransformation T1(xr[cell_node_list[3]], xr[cell_node_list[2]],
+                                                     xr[cell_node_list[1]], xr[cell_node_list[0]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T1.jacobianDeterminant() * f(T1(xi));
+                        }
+                      }
+                    } else if (cell_node_list.size() == 6) {
+                      auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor(4));
+                      {   // top pyramid
+                        PyramidTransformation T0(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]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi));
+                        }
+                      }
+                      {   // bottom pyramid
+                        PyramidTransformation T1(xr[cell_node_list[4]], xr[cell_node_list[3]], xr[cell_node_list[2]],
+                                                 xr[cell_node_list[1]], xr[cell_node_list[0]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi));
+                        }
+                      }
+                    } else {
+                      INFO("Diamond cells with more than 6 vertices are not tested");
+                      REQUIRE(false);
+                    }
+                    break;
+                  }
+                  default: {
+                    INFO("Diamond cells not tested yet");
+                    REQUIRE(cell_type[cell_id] != CellType::Diamond);
+                  }
+                  }
+                  int_f[cell_id] = integral;
+                });
+
+              return int_f;
+            }();
+
+            SECTION("all cells")
+            {
+              SECTION("CellValue")
+              {
+                CellValue<R2> 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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+                }
+
+                REQUIRE(error == 0);
+              }
+
+              SECTION("Array")
+              {
+                Array<R2> 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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+                }
+
+                REQUIRE(error == 0);
+              }
+
+              SECTION("SmallArray")
+              {
+                SmallArray<R2> 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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+                }
+
+                REQUIRE(error == 0);
+              }
+            }
+
+            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<R2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, cell_list);
+
+                double error = 0;
+                for (size_t i = 0; i < cell_list.size(); ++i) {
+                  error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]);
+                }
+
+                REQUIRE(error == 0);
+              }
+
+              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<R2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, cell_list);
+
+                double error = 0;
+                for (size_t i = 0; i < cell_list.size(); ++i) {
+                  error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]);
+                }
+
+                REQUIRE(error == 0);
+              }
+            }
+          }
+
+          SECTION("tensorial formula")
+          {
+            Array<const R2> int_f_per_cell = [=] {
+              Array<R2> int_f(mesh->numberOfCells());
+              auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+              auto cell_type           = mesh->connectivity().cellType();
+
+              auto qf = QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(4));
+
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+                  auto cell_node_list = cell_to_node_matrix[cell_id];
+                  auto xr             = mesh->xr();
+                  R2 integral         = zero;
+
+                  switch (cell_type[cell_id]) {
+                  case CellType::Tetrahedron: {
+                    CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+                                         xr[cell_node_list[2]], xr[cell_node_list[3]], xr[cell_node_list[3]],
+                                         xr[cell_node_list[3]], xr[cell_node_list[3]]);
+
+                    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::Pyramid: {
+                    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[4]],
+                                         xr[cell_node_list[4]], xr[cell_node_list[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: {
+                    CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+                                         xr[cell_node_list[2]], xr[cell_node_list[3]], xr[cell_node_list[4]],
+                                         xr[cell_node_list[5]], xr[cell_node_list[5]]);
+                    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]]);
+
+                    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::Diamond: {
+                    if (cell_node_list.size() == 5) {
+                      {   // top tetrahedron
+                        CubeTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]], xr[cell_node_list[3]],
+                                              xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[4]],
+                                              xr[cell_node_list[4]], xr[cell_node_list[4]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi));
+                        }
+                      }
+                      {   // bottom tetrahedron
+                        CubeTransformation T1(xr[cell_node_list[3]], xr[cell_node_list[2]], xr[cell_node_list[1]],
+                                              xr[cell_node_list[1]], xr[cell_node_list[0]], xr[cell_node_list[0]],
+                                              xr[cell_node_list[0]], xr[cell_node_list[0]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi));
+                        }
+                      }
+                    } else if (cell_node_list.size() == 6) {
+                      {   // top pyramid
+                        CubeTransformation T0(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[5]],
+                                              xr[cell_node_list[5]], xr[cell_node_list[5]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi));
+                        }
+                      }
+                      {   // bottom pyramid
+                        CubeTransformation T1(xr[cell_node_list[4]], xr[cell_node_list[3]], xr[cell_node_list[2]],
+                                              xr[cell_node_list[1]], xr[cell_node_list[0]], xr[cell_node_list[0]],
+                                              xr[cell_node_list[0]], xr[cell_node_list[0]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi));
+                        }
+                      }
+                    } else {
+                      INFO("Diamond cells with more than 6 vertices are not tested");
+                      REQUIRE(false);
+                    }
+                    break;
+                  }
+                  default: {
+                    INFO("Diamond cells not tested yet");
+                    REQUIRE(cell_type[cell_id] != CellType::Diamond);
+                  }
+                  }
+                  int_f[cell_id] = integral;
+                });
+
+              return int_f;
+            }();
+
+            SECTION("all cells")
+            {
+              SECTION("CellValue")
+              {
+                CellValue<R2> values(mesh->connectivity());
+                CellIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10),
+                                            *mesh, values);
+
+                auto cell_type = mesh->connectivity().cellType();
+                double error   = 0;
+                for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+                  error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+                }
+
+                REQUIRE(error == Catch::Approx(0).margin(1E-10));
+              }
+
+              SECTION("Array")
+              {
+                Array<R2> 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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+                }
+
+                REQUIRE(error == Catch::Approx(0).margin(1E-10));
+              }
+
+              SECTION("SmallArray")
+              {
+                SmallArray<R2> 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 += l2Norm(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<R2> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, cell_list);
+
+                double error = 0;
+                for (size_t i = 0; i < cell_list.size(); ++i) {
+                  error += l2Norm(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<R2> values =
+                  CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, cell_list);
+
+                double error = 0;
+                for (size_t i = 0; i < cell_list.size(); ++i) {
+                  error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]);
+                }
+
+                REQUIRE(error == Catch::Approx(0).margin(1E-10));
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+
+  SECTION("matrix")
+  {
+    using R2x2  = TinyMatrix<2>;
+    auto l2Norm = [](const R2x2& A) -> double {
+      return std::sqrt(A(0, 0) * A(0, 0) + A(1, 0) * A(1, 0) + A(0, 1) * A(0, 1) + A(1, 1) * A(1, 1));
+    };
+
+    SECTION("1D")
+    {
+      using R1 = TinyVector<1>;
+
+      const auto mesh = MeshDataBaseForTests::get().unordered1DMesh();
+      auto f          = [](const R1& x) -> R2x2 { return R2x2{x[0] * x[0] + 1, 2 * x[0], 3 - x[0], 3 * x[0] - 1}; };
+
+      Array<const R2x2> int_f_per_cell = [=] {
+        Array<R2x2> 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] = R2x2{1. / 3 * (x_right * x_right * x_right - x_left * x_left * x_left) + x_right - x_left,
+                                  x_right * x_right - x_left * x_left,
+                                  -0.5 * (x_right * x_right - x_left * x_left) + 3 * (x_right - x_left),
+                                  3. / 2 * (x_right * x_right - x_left * x_left) - (x_right - x_left)};
+          });
+
+        return int_f;
+      }();
+
+      SECTION("direct formula")
+      {
+        SECTION("all cells")
+        {
+          SECTION("CellValue")
+          {
+            CellValue<R2x2> 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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
 
-          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));
           }
 
-          REQUIRE(error == Catch::Approx(0).margin(1E-10));
-        }
+          SECTION("Array")
+          {
+            Array<R2x2> values(mesh->numberOfCells());
 
-        SECTION("Array")
-        {
-          Array<double> values(mesh->numberOfCells());
+            CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
 
-          CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
+            double error = 0;
+            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+              error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
 
-          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));
           }
 
-          REQUIRE(error == Catch::Approx(0).margin(1E-10));
-        }
+          SECTION("SmallArray")
+          {
+            SmallArray<R2x2> values(mesh->numberOfCells());
+            CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
 
-        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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
 
-          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));
           }
-
-          REQUIRE(error == Catch::Approx(0).margin(1E-10));
         }
-      }
 
-      SECTION("cell list")
-      {
-        SECTION("Array")
+        SECTION("cell list")
         {
-          Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
-
+          SECTION("Array")
           {
-            size_t k = 0;
-            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
-              cell_list[k] = cell_id;
+            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());
             }
 
-            REQUIRE(k == cell_list.size());
-          }
+            Array<R2x2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list);
 
-          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 += l2Norm(int_f_per_cell[cell_list[i]] - values[i]);
+            }
 
-          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));
           }
 
-          REQUIRE(error == Catch::Approx(0).margin(1E-10));
-        }
+          SECTION("SmallArray")
+          {
+            SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
 
-        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;
+              }
 
-          {
-            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());
             }
 
-            REQUIRE(k == cell_list.size());
-          }
+            SmallArray<R2x2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list);
 
-          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 += l2Norm(int_f_per_cell[cell_list[i]] - values[i]);
+            }
 
-          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));
           }
-
-          REQUIRE(error == Catch::Approx(0).margin(1E-10));
         }
       }
-    }
 
-    SECTION("tensorial formula")
-    {
-      SECTION("all cells")
+      SECTION("tensorial formula")
       {
-        SECTION("CellValue")
+        SECTION("all cells")
         {
-          CellValue<double> values(mesh->connectivity());
-          CellIntegrator::integrateTo([=](const R1 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh,
-                                      values);
+          SECTION("CellValue")
+          {
+            CellValue<R2x2> 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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
 
-          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));
           }
 
-          REQUIRE(error == Catch::Approx(0).margin(1E-10));
-        }
+          SECTION("Array")
+          {
+            Array<R2x2> values(mesh->numberOfCells());
 
-        SECTION("Array")
-        {
-          Array<double> values(mesh->numberOfCells());
+            CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
 
-          CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
+            double error = 0;
+            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+              error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
 
-          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));
           }
 
-          REQUIRE(error == Catch::Approx(0).margin(1E-10));
-        }
+          SECTION("SmallArray")
+          {
+            SmallArray<R2x2> values(mesh->numberOfCells());
+            CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
 
-        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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
 
-          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));
           }
-
-          REQUIRE(error == Catch::Approx(0).margin(1E-10));
         }
-      }
 
-      SECTION("cell list")
-      {
-        SECTION("Array")
+        SECTION("cell list")
         {
-          Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
-
+          SECTION("Array")
           {
-            size_t k = 0;
-            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
-              cell_list[k] = cell_id;
+            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());
             }
 
-            REQUIRE(k == cell_list.size());
-          }
+            Array<R2x2> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list);
 
-          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 += l2Norm(int_f_per_cell[cell_list[i]] - values[i]);
+            }
 
-          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));
           }
 
-          REQUIRE(error == Catch::Approx(0).margin(1E-10));
-        }
+          SECTION("SmallArray")
+          {
+            SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
 
-        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;
+              }
 
-          {
-            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());
             }
 
-            REQUIRE(k == cell_list.size());
-          }
+            SmallArray<R2x2> values =
+              CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list);
 
-          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 += l2Norm(int_f_per_cell[cell_list[i]] - values[i]);
+            }
 
-          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));
           }
-
-          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));
+    SECTION("2D")
+    {
+      using R2 = TinyVector<2>;
+
+      const auto mesh = MeshDataBaseForTests::get().hybrid2DMesh();
+
+      auto f = [](const R2& X) -> R2x2 {
+        const double x = X[0];
+        const double y = X[1];
+        return R2x2{x * x + 2 * x * y + 3 * y * y + 2, 2 * x + 3 * y * y, 2 * x - 3 * y, 3 * x * x - 2 * y};
+      };
+
+      Array<const R2x2> int_f_per_cell = [=] {
+        Array<R2x2> 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();
+            R2x2 integral       = zero;
+
+            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;
             }
-            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));
+            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;
             }
-            break;
-          }
-          default: {
-            throw UnexpectedError("invalid cell type in 2d");
-          }
-          }
-          int_f[cell_id] = integral;
-        });
+            default: {
+              throw UnexpectedError("invalid cell type in 2d");
+            }
+            }
+            int_f[cell_id] = integral;
+          });
 
-      return int_f;
-    }();
+        return int_f;
+      }();
 
-    SECTION("direct formula")
-    {
-      SECTION("all cells")
+      SECTION("direct formula")
       {
-        SECTION("CellValue")
+        SECTION("all cells")
         {
-          CellValue<double> values(mesh->connectivity());
-          CellIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values);
+          SECTION("CellValue")
+          {
+            CellValue<R2x2> 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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
 
-          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));
           }
 
-          REQUIRE(error == Catch::Approx(0).margin(1E-10));
-        }
+          SECTION("Array")
+          {
+            Array<R2x2> values(mesh->numberOfCells());
 
-        SECTION("Array")
-        {
-          Array<double> values(mesh->numberOfCells());
+            CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
 
-          CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
+            double error = 0;
+            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+              error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
 
-          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));
           }
 
-          REQUIRE(error == Catch::Approx(0).margin(1E-10));
-        }
+          SECTION("SmallArray")
+          {
+            SmallArray<R2x2> values(mesh->numberOfCells());
+            CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
 
-        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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
 
-          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));
           }
-
-          REQUIRE(error == Catch::Approx(0).margin(1E-10));
         }
-      }
 
-      SECTION("cell list")
-      {
-        SECTION("Array")
+        SECTION("cell list")
         {
-          Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
-
+          SECTION("Array")
           {
-            size_t k = 0;
-            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
-              cell_list[k] = cell_id;
+            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());
             }
 
-            REQUIRE(k == cell_list.size());
-          }
+            Array<R2x2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list);
 
-          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 += l2Norm(int_f_per_cell[cell_list[i]] - values[i]);
+            }
 
-          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));
           }
 
-          REQUIRE(error == Catch::Approx(0).margin(1E-10));
-        }
+          SECTION("SmallArray")
+          {
+            SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
 
-        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;
+              }
 
-          {
-            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());
             }
 
-            REQUIRE(k == cell_list.size());
-          }
+            SmallArray<R2x2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list);
 
-          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 += l2Norm(int_f_per_cell[cell_list[i]] - values[i]);
+            }
 
-          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));
           }
-
-          REQUIRE(error == Catch::Approx(0).margin(1E-10));
         }
       }
-    }
 
-    SECTION("tensorial formula")
-    {
-      SECTION("all cells")
+      SECTION("tensorial formula")
       {
-        SECTION("CellValue")
+        SECTION("all cells")
         {
-          CellValue<double> values(mesh->connectivity());
-          CellIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh,
-                                      values);
+          SECTION("CellValue")
+          {
+            CellValue<R2x2> 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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
 
-          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));
           }
 
-          REQUIRE(error == Catch::Approx(0).margin(1E-10));
-        }
+          SECTION("Array")
+          {
+            Array<R2x2> values(mesh->numberOfCells());
 
-        SECTION("Array")
-        {
-          Array<double> values(mesh->numberOfCells());
+            CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
 
-          CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
+            double error = 0;
+            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+              error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
 
-          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));
           }
 
-          REQUIRE(error == Catch::Approx(0).margin(1E-10));
-        }
+          SECTION("SmallArray")
+          {
+            SmallArray<R2x2> values(mesh->numberOfCells());
+            CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
 
-        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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+            }
 
-          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));
           }
-
-          REQUIRE(error == Catch::Approx(0).margin(1E-10));
         }
-      }
 
-      SECTION("cell list")
-      {
-        SECTION("Array")
+        SECTION("cell list")
         {
-          Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
-
+          SECTION("Array")
           {
-            size_t k = 0;
-            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
-              cell_list[k] = cell_id;
+            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());
             }
 
-            REQUIRE(k == cell_list.size());
-          }
+            Array<R2x2> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list);
 
-          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 += l2Norm(int_f_per_cell[cell_list[i]] - values[i]);
+            }
 
-          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));
           }
 
-          REQUIRE(error == Catch::Approx(0).margin(1E-10));
-        }
+          SECTION("SmallArray")
+          {
+            SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
 
-        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;
+              }
 
-          {
-            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());
             }
 
-            REQUIRE(k == cell_list.size());
-          }
+            SmallArray<R2x2> values =
+              CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list);
 
-          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 += l2Norm(int_f_per_cell[cell_list[i]] - values[i]);
+            }
 
-          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));
           }
-
-          REQUIRE(error == Catch::Approx(0).margin(1E-10));
         }
       }
     }
-  }
 
-  SECTION("3D")
-  {
-    using R3 = TinyVector<3>;
+    SECTION("3D")
+    {
+      using R3 = TinyVector<3>;
 
-    auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh();
+      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;
-    };
+      auto f = [](const R3& X) -> R2x2 {
+        const double x = X[0];
+        const double y = X[1];
+        const double z = X[2];
+        return R2x2{x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1, 2 * x * y - 3 * y * z + 2 * x,
+                    x * x - 2 * x * z - 3 * y, 3 * z + x * y};
+      };
 
-    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", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh)));
+      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", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh)));
 
-    for (auto mesh_info : mesh_list) {
-      auto mesh_name = mesh_info.first;
-      auto mesh      = mesh_info.second;
+      for (auto mesh_info : mesh_list) {
+        auto mesh_name = mesh_info.first;
+        auto mesh      = mesh_info.second;
 
-      SECTION(mesh_name)
-      {
-        SECTION("direct formula")
+        SECTION(mesh_name)
         {
-          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));
+          SECTION("direct formula")
+          {
+            Array<const R2x2> int_f_per_cell = [=] {
+              Array<R2x2> 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();
+                  R2x2 integral       = zero;
+
+                  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));
+                    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;
                   }
-                  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));
+                  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));
+                    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;
                   }
-                  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));
+                  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;
                   }
-                  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));
+                  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;
                   }
-                  break;
-                }
-                case CellType::Diamond: {
-                  if (cell_node_list.size() == 5) {
-                    auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor(4));
-                    {   // top tetrahedron
-                      TetrahedronTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]], xr[cell_node_list[3]],
-                                                   xr[cell_node_list[4]]);
-
-                      for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
-                        const auto& xi = qf.point(i);
-                        integral += qf.weight(i) * T0.jacobianDeterminant() * f(T0(xi));
+                  case CellType::Diamond: {
+                    if (cell_node_list.size() == 5) {
+                      auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor(4));
+                      {   // top tetrahedron
+                        TetrahedronTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]],
+                                                     xr[cell_node_list[3]], xr[cell_node_list[4]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T0.jacobianDeterminant() * f(T0(xi));
+                        }
                       }
-                    }
-                    {   // bottom tetrahedron
-                      TetrahedronTransformation T1(xr[cell_node_list[3]], xr[cell_node_list[2]], xr[cell_node_list[1]],
-                                                   xr[cell_node_list[0]]);
-
-                      for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
-                        const auto& xi = qf.point(i);
-                        integral += qf.weight(i) * T1.jacobianDeterminant() * f(T1(xi));
+                      {   // bottom tetrahedron
+                        TetrahedronTransformation T1(xr[cell_node_list[3]], xr[cell_node_list[2]],
+                                                     xr[cell_node_list[1]], xr[cell_node_list[0]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T1.jacobianDeterminant() * f(T1(xi));
+                        }
                       }
-                    }
-                  } else if (cell_node_list.size() == 6) {
-                    auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor(4));
-                    {   // top pyramid
-                      PyramidTransformation T0(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]]);
-
-                      for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
-                        const auto& xi = qf.point(i);
-                        integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi));
+                    } else if (cell_node_list.size() == 6) {
+                      auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor(4));
+                      {   // top pyramid
+                        PyramidTransformation T0(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]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi));
+                        }
                       }
-                    }
-                    {   // bottom pyramid
-                      PyramidTransformation T1(xr[cell_node_list[4]], xr[cell_node_list[3]], xr[cell_node_list[2]],
-                                               xr[cell_node_list[1]], xr[cell_node_list[0]]);
-
-                      for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
-                        const auto& xi = qf.point(i);
-                        integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi));
+                      {   // bottom pyramid
+                        PyramidTransformation T1(xr[cell_node_list[4]], xr[cell_node_list[3]], xr[cell_node_list[2]],
+                                                 xr[cell_node_list[1]], xr[cell_node_list[0]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi));
+                        }
                       }
+                    } else {
+                      INFO("Diamond cells with more than 6 vertices are not tested");
+                      REQUIRE(false);
                     }
-                  } else {
-                    INFO("Diamond cells with more than 6 vertices are not tested");
-                    REQUIRE(false);
+                    break;
                   }
-                  break;
-                }
-                default: {
-                  INFO("Diamond cells not tested yet");
-                  REQUIRE(cell_type[cell_id] != CellType::Diamond);
-                }
-                }
-                int_f[cell_id] = integral;
-              });
+                  default: {
+                    INFO("Diamond cells not tested yet");
+                    REQUIRE(cell_type[cell_id] != CellType::Diamond);
+                  }
+                  }
+                  int_f[cell_id] = integral;
+                });
 
-            return int_f;
-          }();
+              return int_f;
+            }();
 
-          SECTION("all cells")
-          {
-            SECTION("CellValue")
+            SECTION("all cells")
             {
-              CellValue<double> values(mesh->connectivity());
-              CellIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh,
-                                          values);
+              SECTION("CellValue")
+              {
+                CellValue<R2x2> 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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+                }
 
-              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 == 0);
               }
 
-              REQUIRE(error == 0);
-            }
+              SECTION("Array")
+              {
+                Array<R2x2> values(mesh->numberOfCells());
 
-            SECTION("Array")
-            {
-              Array<double> values(mesh->numberOfCells());
+                CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
 
-              CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+                double error = 0;
+                for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+                  error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+                }
 
-              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 == 0);
               }
 
-              REQUIRE(error == 0);
-            }
+              SECTION("SmallArray")
+              {
+                SmallArray<R2x2> values(mesh->numberOfCells());
+                CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
 
-            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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+                }
 
-              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 == 0);
               }
-
-              REQUIRE(error == 0);
             }
-          }
 
-          SECTION("cell list")
-          {
-            SECTION("Array")
+            SECTION("cell list")
             {
-              Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
-
+              SECTION("Array")
               {
-                size_t k = 0;
-                for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
-                  cell_list[k] = cell_id;
+                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());
                 }
 
-                REQUIRE(k == cell_list.size());
-              }
+                Array<R2x2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, cell_list);
 
-              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 += l2Norm(int_f_per_cell[cell_list[i]] - values[i]);
+                }
 
-              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 == 0);
               }
 
-              REQUIRE(error == 0);
-            }
+              SECTION("SmallArray")
+              {
+                SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
 
-            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;
+                  }
 
-              {
-                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());
                 }
 
-                REQUIRE(k == cell_list.size());
-              }
+                SmallArray<R2x2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, cell_list);
 
-              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 += l2Norm(int_f_per_cell[cell_list[i]] - values[i]);
+                }
 
-              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 == 0);
               }
-
-              REQUIRE(error == 0);
             }
           }
-        }
 
-        SECTION("tensorial formula")
-        {
-          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();
-
-            auto qf = QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(4));
-
-            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: {
-                  CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
-                                       xr[cell_node_list[2]], xr[cell_node_list[3]], xr[cell_node_list[3]],
-                                       xr[cell_node_list[3]], xr[cell_node_list[3]]);
-
-                  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::Pyramid: {
-                  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[4]],
-                                       xr[cell_node_list[4]], xr[cell_node_list[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: {
-                  CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
-                                       xr[cell_node_list[2]], xr[cell_node_list[3]], xr[cell_node_list[4]],
-                                       xr[cell_node_list[5]], xr[cell_node_list[5]]);
-                  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]]);
-
-                  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::Diamond: {
-                  if (cell_node_list.size() == 5) {
-                    {   // top tetrahedron
-                      CubeTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]], xr[cell_node_list[3]],
-                                            xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[4]],
-                                            xr[cell_node_list[4]], xr[cell_node_list[4]]);
-
-                      for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
-                        const auto& xi = qf.point(i);
-                        integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi));
-                      }
+          SECTION("tensorial formula")
+          {
+            Array<const R2x2> int_f_per_cell = [=] {
+              Array<R2x2> int_f(mesh->numberOfCells());
+              auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+              auto cell_type           = mesh->connectivity().cellType();
+
+              auto qf = QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(4));
+
+              parallel_for(
+                mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+                  auto cell_node_list = cell_to_node_matrix[cell_id];
+                  auto xr             = mesh->xr();
+                  R2x2 integral       = zero;
+
+                  switch (cell_type[cell_id]) {
+                  case CellType::Tetrahedron: {
+                    CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+                                         xr[cell_node_list[2]], xr[cell_node_list[3]], xr[cell_node_list[3]],
+                                         xr[cell_node_list[3]], xr[cell_node_list[3]]);
+
+                    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));
                     }
-                    {   // bottom tetrahedron
-                      CubeTransformation T1(xr[cell_node_list[3]], xr[cell_node_list[2]], xr[cell_node_list[1]],
-                                            xr[cell_node_list[1]], xr[cell_node_list[0]], xr[cell_node_list[0]],
-                                            xr[cell_node_list[0]], xr[cell_node_list[0]]);
-
-                      for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
-                        const auto& xi = qf.point(i);
-                        integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi));
-                      }
+                    break;
+                  }
+                  case CellType::Pyramid: {
+                    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[4]],
+                                         xr[cell_node_list[4]], xr[cell_node_list[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));
                     }
-                  } else if (cell_node_list.size() == 6) {
-                    {   // top pyramid
-                      CubeTransformation T0(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[5]],
-                                            xr[cell_node_list[5]], xr[cell_node_list[5]]);
-
-                      for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
-                        const auto& xi = qf.point(i);
-                        integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi));
-                      }
+                    break;
+                  }
+                  case CellType::Prism: {
+                    CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+                                         xr[cell_node_list[2]], xr[cell_node_list[3]], xr[cell_node_list[4]],
+                                         xr[cell_node_list[5]], xr[cell_node_list[5]]);
+                    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]]);
+
+                    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));
                     }
-                    {   // bottom pyramid
-                      CubeTransformation T1(xr[cell_node_list[4]], xr[cell_node_list[3]], xr[cell_node_list[2]],
-                                            xr[cell_node_list[1]], xr[cell_node_list[0]], xr[cell_node_list[0]],
-                                            xr[cell_node_list[0]], xr[cell_node_list[0]]);
-
-                      for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
-                        const auto& xi = qf.point(i);
-                        integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi));
+                    break;
+                  }
+                  case CellType::Diamond: {
+                    if (cell_node_list.size() == 5) {
+                      {   // top tetrahedron
+                        CubeTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]], xr[cell_node_list[3]],
+                                              xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[4]],
+                                              xr[cell_node_list[4]], xr[cell_node_list[4]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi));
+                        }
                       }
+                      {   // bottom tetrahedron
+                        CubeTransformation T1(xr[cell_node_list[3]], xr[cell_node_list[2]], xr[cell_node_list[1]],
+                                              xr[cell_node_list[1]], xr[cell_node_list[0]], xr[cell_node_list[0]],
+                                              xr[cell_node_list[0]], xr[cell_node_list[0]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi));
+                        }
+                      }
+                    } else if (cell_node_list.size() == 6) {
+                      {   // top pyramid
+                        CubeTransformation T0(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[5]],
+                                              xr[cell_node_list[5]], xr[cell_node_list[5]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi));
+                        }
+                      }
+                      {   // bottom pyramid
+                        CubeTransformation T1(xr[cell_node_list[4]], xr[cell_node_list[3]], xr[cell_node_list[2]],
+                                              xr[cell_node_list[1]], xr[cell_node_list[0]], xr[cell_node_list[0]],
+                                              xr[cell_node_list[0]], xr[cell_node_list[0]]);
+
+                        for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+                          const auto& xi = qf.point(i);
+                          integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi));
+                        }
+                      }
+                    } else {
+                      INFO("Diamond cells with more than 6 vertices are not tested");
+                      REQUIRE(false);
                     }
-                  } else {
-                    INFO("Diamond cells with more than 6 vertices are not tested");
-                    REQUIRE(false);
+                    break;
                   }
-                  break;
-                }
-                default: {
-                  INFO("Diamond cells not tested yet");
-                  REQUIRE(cell_type[cell_id] != CellType::Diamond);
-                }
-                }
-                int_f[cell_id] = integral;
-              });
+                  default: {
+                    INFO("Diamond cells not tested yet");
+                    REQUIRE(cell_type[cell_id] != CellType::Diamond);
+                  }
+                  }
+                  int_f[cell_id] = integral;
+                });
 
-            return int_f;
-          }();
+              return int_f;
+            }();
 
-          SECTION("all cells")
-          {
-            SECTION("CellValue")
+            SECTION("all cells")
             {
-              CellValue<double> values(mesh->connectivity());
-              CellIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10),
-                                          *mesh, values);
+              SECTION("CellValue")
+              {
+                CellValue<R2x2> values(mesh->connectivity());
+                CellIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10),
+                                            *mesh, values);
+
+                auto cell_type = mesh->connectivity().cellType();
+                double error   = 0;
+                for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+                  error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+                }
 
-              auto cell_type = mesh->connectivity().cellType();
-              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));
               }
 
-              REQUIRE(error == Catch::Approx(0).margin(1E-10));
-            }
+              SECTION("Array")
+              {
+                Array<R2x2> values(mesh->numberOfCells());
 
-            SECTION("Array")
-            {
-              Array<double> values(mesh->numberOfCells());
+                CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
 
-              CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+                double error = 0;
+                for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+                  error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+                }
 
-              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));
               }
 
-              REQUIRE(error == Catch::Approx(0).margin(1E-10));
-            }
+              SECTION("SmallArray")
+              {
+                SmallArray<R2x2> values(mesh->numberOfCells());
+                CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
 
-            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 += l2Norm(int_f_per_cell[cell_id] - values[cell_id]);
+                }
 
-              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));
               }
-
-              REQUIRE(error == Catch::Approx(0).margin(1E-10));
             }
-          }
 
-          SECTION("cell list")
-          {
-            SECTION("Array")
+            SECTION("cell list")
             {
-              Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
-
+              SECTION("Array")
               {
-                size_t k = 0;
-                for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
-                  cell_list[k] = cell_id;
+                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());
                 }
 
-                REQUIRE(k == cell_list.size());
-              }
+                Array<R2x2> values =
+                  CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, cell_list);
 
-              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 += l2Norm(int_f_per_cell[cell_list[i]] - values[i]);
+                }
 
-              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));
               }
 
-              REQUIRE(error == Catch::Approx(0).margin(1E-10));
-            }
+              SECTION("SmallArray")
+              {
+                SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
 
-            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;
+                  }
 
-              {
-                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());
                 }
 
-                REQUIRE(k == cell_list.size());
-              }
+                SmallArray<R2x2> values =
+                  CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, cell_list);
 
-              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 += l2Norm(int_f_per_cell[cell_list[i]] - values[i]);
+                }
 
-              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));
               }
-
-              REQUIRE(error == Catch::Approx(0).margin(1E-10));
             }
           }
         }
-- 
GitLab