diff --git a/src/language/utils/IntegrateCellValue.hpp b/src/language/utils/IntegrateCellValue.hpp
index c4eafc151712219101a42d28ba2e8eb774ab0060..cd8bacc072bef0f5c68b1702fa9fef33c6d94f7d 100644
--- a/src/language/utils/IntegrateCellValue.hpp
+++ b/src/language/utils/IntegrateCellValue.hpp
@@ -34,7 +34,7 @@ class IntegrateCellValue<OutputType(InputType)>
             const Array<const CellId>& list_of_cells)
   {
     return IntegrateOnCells<OutputType(const InputType)>::integrate(function_symbol_id, quadrature_descriptor, mesh,
-                                                                    list_of_cells);
+                                                                    Array<const CellId>{list_of_cells});
   }
 
   template <typename MeshType>
@@ -44,8 +44,7 @@ class IntegrateCellValue<OutputType(InputType)>
             const MeshType& mesh,
             const Array<CellId>& list_of_cells)
   {
-    return IntegrateOnCells<OutputType(const InputType)>::integrate(function_symbol_id, quadrature_descriptor, mesh,
-                                                                    Array<const CellId>{list_of_cells});
+    return integrate(function_symbol_id, quadrature_descriptor, mesh, Array<const CellId>{list_of_cells});
   }
 };
 
diff --git a/src/language/utils/IntegrateOnCells.hpp b/src/language/utils/IntegrateOnCells.hpp
index 8b6008d73672a24ee5473369250f00ce05342d7e..fa47bb7856ec695d87d7702d9c561a2ebfe6ca6b 100644
--- a/src/language/utils/IntegrateOnCells.hpp
+++ b/src/language/utils/IntegrateOnCells.hpp
@@ -56,10 +56,10 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
   class CellList
   {
    private:
-    const Array<CellId>& m_cell_list;
+    const Array<const CellId>& m_cell_list;
 
    public:
-    using index_type = Array<CellId>::index_type;
+    using index_type = Array<const CellId>::index_type;
 
     PUGS_INLINE
     CellId
@@ -76,7 +76,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
     }
 
     PUGS_INLINE
-    CellList(const Array<CellId>& cell_list) : m_cell_list{cell_list} {}
+    CellList(const Array<const CellId>& cell_list) : m_cell_list{cell_list} {}
   };
 
   template <typename MeshType, typename OutputArrayT, typename ListTypeT>
@@ -150,10 +150,12 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
           }
           break;
         }
+          // LCOV_EXCL_START
         default: {
           invalid_cell = std::make_pair(true, cell_id);
           break;
         }
+          // LCOV_EXCL_STOP
         }
       } else if constexpr (Dimension == 2) {
         switch (cell_type[cell_id]) {
@@ -181,10 +183,12 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
           }
           break;
         }
+          // LCOV_EXCL_START
         default: {
           invalid_cell = std::make_pair(true, cell_id);
           break;
         }
+          // LCOV_EXCL_STOP
         }
       } else {
         static_assert(Dimension == 3);
@@ -218,8 +222,10 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
               value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
             }
           } else {
+            // LCOV_EXCL_START
             throw NotImplementedError("integration on pyramid with non-quadrangular base (number of nodes " +
                                       std::to_string(cell_to_node.size()) + ")");
+            // LCOV_EXCL_STOP
           }
           break;
         }
@@ -279,8 +285,10 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
               }
             }
           } else {
+            // LCOV_EXCL_START
             throw NotImplementedError("integration on diamond with non-quadrangular base (" +
                                       std::to_string(cell_to_node.size()) + ")");
+            // LCOV_EXCL_STOP
           }
           break;
         }
@@ -311,21 +319,25 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
           }
           break;
         }
+          // LCOV_EXCL_START
         default: {
           invalid_cell = std::make_pair(true, cell_id);
           break;
         }
+          // LCOV_EXCL_STOP
         }
       }
 
       tokens.release(t);
     });
 
+    // LCOV_EXCL_START
     if (invalid_cell.first) {
       std::ostringstream os;
       os << "invalid cell type for integration: " << name(cell_type[invalid_cell.second]);
       throw UnexpectedError(os.str());
     }
+    // LCOV_EXCL_STOP
   }
 
   template <typename MeshType, typename OutputArrayT, typename ListTypeT>
@@ -391,10 +403,12 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
           }
           break;
         }
+          // LCOV_EXCL_START
         default: {
           invalid_cell = std::make_pair(true, cell_id);
           break;
         }
+          // LCOV_EXCL_STOP
         }
       } else if constexpr (Dimension == 2) {
         switch (cell_type[cell_id]) {
@@ -423,10 +437,12 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
           }
           break;
         }
+          // LCOV_EXCL_START
         default: {
           invalid_cell = std::make_pair(true, cell_id);
           break;
         }
+          // LCOV_EXCL_STOP
         }
       } else {
         static_assert(Dimension == 3);
@@ -458,7 +474,9 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
               value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
             }
           } else {
+            // LCOV_EXCL_START
             throw NotImplementedError("integration on pyramid with non-quadrangular base");
+            // LCOV_EXCL_STOP
           }
           break;
         }
@@ -512,8 +530,10 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
               }
             }
           } else {
+            // LCOV_EXCL_START
             throw NotImplementedError("integration on diamond with non-quadrangular base (" +
                                       std::to_string(cell_to_node.size()) + ")");
+            // LCOV_EXCL_STOP
           }
           break;
         }
@@ -544,21 +564,25 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
           }
           break;
         }
+          // LCOV_EXCL_START
         default: {
           invalid_cell = std::make_pair(true, cell_id);
           break;
         }
+          // LCOV_EXCL_STOP
         }
       }
 
       tokens.release(t);
     });
 
+    // LCOV_EXCL_START
     if (invalid_cell.first) {
       std::ostringstream os;
       os << "invalid cell type for integration: " << name(cell_type[invalid_cell.second]);
       throw UnexpectedError(os.str());
     }
+    // LCOV_EXCL_STOP
   }
 
  public:
@@ -585,7 +609,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
   integrate(const FunctionSymbolId& function_symbol_id,
             const IQuadratureDescriptor& quadrature_descriptor,
             const MeshType& mesh,
-            const ArrayT<CellId>& cell_list)
+            const ArrayT<const CellId>& cell_list)
   {
     ArrayT<OutputType> value(size(cell_list));
     if (quadrature_descriptor.isTensorial()) {
@@ -598,6 +622,16 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
 
     return value;
   }
+
+  template <typename MeshType, template <typename DataType> typename ArrayT>
+  static PUGS_INLINE ArrayT<OutputType>
+  integrate(const FunctionSymbolId& function_symbol_id,
+            const IQuadratureDescriptor& quadrature_descriptor,
+            const MeshType& mesh,
+            const ArrayT<CellId>& cell_list)
+  {
+    return integrate(function_symbol_id, quadrature_descriptor, mesh, ArrayT<const CellId>{cell_list});
+  }
 };
 
 #endif   // INTEGRATE_ON_CELLS_HPP
diff --git a/src/mesh/ItemArray.hpp b/src/mesh/ItemArray.hpp
index 5a16861a37b5642ba10977d8c9e2d5ecd25fd194..3bd7de633659e81eb8fc178f275ee76a449048cf 100644
--- a/src/mesh/ItemArray.hpp
+++ b/src/mesh/ItemArray.hpp
@@ -36,6 +36,12 @@ class ItemArray
   // Allow const std:weak_ptr version to access our data
   friend ItemArray<std::add_const_t<DataType>, item_type, ConnectivityWeakPtr>;
 
+  // Allow non-const std:shared_ptr version to access our data
+  friend ItemArray<std::remove_const_t<DataType>, item_type, ConnectivitySharedPtr>;
+
+  // Allow non-const std:weak_ptr version to access our data
+  friend ItemArray<std::remove_const_t<DataType>, item_type, ConnectivityWeakPtr>;
+
  public:
   friend PUGS_INLINE ItemArray<std::remove_const_t<DataType>, item_type, ConnectivityPtr>
   copy(const ItemArray<DataType, item_type, ConnectivityPtr>& source)
diff --git a/src/mesh/ItemValue.hpp b/src/mesh/ItemValue.hpp
index d48468f3931a2dc0dc546cb2a2dfee630df9e87e..981ebcc196d4d637d9e775231fe14edc8a70a2f2 100644
--- a/src/mesh/ItemValue.hpp
+++ b/src/mesh/ItemValue.hpp
@@ -36,6 +36,12 @@ class ItemValue
   // Allow const std:weak_ptr version to access our data
   friend ItemValue<std::add_const_t<DataType>, item_type, ConnectivityWeakPtr>;
 
+  // Allow non-const std:shared_ptr version to access our data
+  friend ItemValue<std::remove_const_t<DataType>, item_type, ConnectivitySharedPtr>;
+
+  // Allow non-const std:weak_ptr version to access our data
+  friend ItemValue<std::remove_const_t<DataType>, item_type, ConnectivityWeakPtr>;
+
  public:
   friend PUGS_INLINE ItemValue<std::remove_const_t<DataType>, item_type, ConnectivityPtr>
   copy(const ItemValue<DataType, item_type, ConnectivityPtr>& source)
diff --git a/src/mesh/SubItemArrayPerItem.hpp b/src/mesh/SubItemArrayPerItem.hpp
index 9fd488bc5e6ebe3ea189368de8eeb67073e590ec..f2103789e9a49c92b921fe2f8b24b3bdaf8f69ec 100644
--- a/src/mesh/SubItemArrayPerItem.hpp
+++ b/src/mesh/SubItemArrayPerItem.hpp
@@ -42,15 +42,13 @@ class SubItemArrayPerItem
   // Allow const std:weak_ptr version to access our data
   friend SubItemArrayPerItem<std::add_const_t<DataType>, ItemOfItem, ConnectivityWeakPtr>;
 
-  // Allow const std:shared_ptr version to access our data
+  // Allow non-const std:shared_ptr version to access our data
   friend SubItemArrayPerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivitySharedPtr>;
 
-  // Allow const std:weak_ptr version to access our data
+  // Allow non-const std:weak_ptr version to access our data
   friend SubItemArrayPerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityWeakPtr>;
 
  public:
-  using ToShared = SubItemArrayPerItem<DataType, ItemOfItem, ConnectivitySharedPtr>;
-
   friend PUGS_INLINE SubItemArrayPerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityPtr>
   copy(SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& source)
   {
@@ -142,7 +140,7 @@ class SubItemArrayPerItem
   void
   fill(const DataType& data) const noexcept
   {
-    static_assert(not std::is_const_v<DataType>, "Cannot modify ItemValue of const");
+    static_assert(not std::is_const_v<DataType>, "Cannot modify SubItemArrayPerItem of const");
     m_values.fill(data);
   }
 
diff --git a/src/mesh/SubItemValuePerItem.hpp b/src/mesh/SubItemValuePerItem.hpp
index b86b075ecbb44ff0f95a1a28b1f8f2468d98b2e3..defdc8fc9711c6be265e022e0bb0c710d078226d 100644
--- a/src/mesh/SubItemValuePerItem.hpp
+++ b/src/mesh/SubItemValuePerItem.hpp
@@ -41,15 +41,13 @@ class SubItemValuePerItem
   // Allow const std:weak_ptr version to access our data
   friend SubItemValuePerItem<std::add_const_t<DataType>, ItemOfItem, ConnectivityWeakPtr>;
 
-  // Allow const std:shared_ptr version to access our data
+  // Allow non-const std:shared_ptr version to access our data
   friend SubItemValuePerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivitySharedPtr>;
 
-  // Allow const std:weak_ptr version to access our data
+  // Allow non-const std:weak_ptr version to access our data
   friend SubItemValuePerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityWeakPtr>;
 
  public:
-  using ToShared = SubItemValuePerItem<DataType, ItemOfItem, ConnectivitySharedPtr>;
-
   friend PUGS_INLINE SubItemValuePerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityPtr>
   copy(SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& source)
   {
diff --git a/src/mesh/Synchronizer.hpp b/src/mesh/Synchronizer.hpp
index adf127f238fc201a11ba5b503aed4b926e772a4d..dc0773ace208cc83e863df0ffada54e1db33f936 100644
--- a/src/mesh/Synchronizer.hpp
+++ b/src/mesh/Synchronizer.hpp
@@ -7,9 +7,13 @@
 #include <utils/Exceptions.hpp>
 #include <utils/Messenger.hpp>
 
+#include <utils/pugs_config.hpp>
+
 #include <iostream>
 #include <map>
 
+#ifdef PUGS_HAS_MPI
+
 class Synchronizer
 {
   template <ItemType item_type>
@@ -289,4 +293,28 @@ class Synchronizer
   }
 };
 
+#else   // PUGS_HAS_MPI
+
+class Synchronizer
+{
+ public:
+  template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+  PUGS_INLINE void
+  synchronize(ItemValue<DataType, item_type, ConnectivityPtr>&)
+  {}
+
+  template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+  PUGS_INLINE void
+  synchronize(ItemArray<DataType, item_type, ConnectivityPtr>&)
+  {}
+
+  PUGS_INLINE
+  Synchronizer()
+  {
+    ;
+  }
+};
+
+#endif   // PUGS_HAS_MPI
+
 #endif   // SYNCHRONIZER_HPP
diff --git a/src/mesh/SynchronizerManager.cpp b/src/mesh/SynchronizerManager.cpp
index b1301bd88a341cae40c348e1239f673636f5b818..a72191fd3bfc001c38e577384115bffaa441633c 100644
--- a/src/mesh/SynchronizerManager.cpp
+++ b/src/mesh/SynchronizerManager.cpp
@@ -18,6 +18,7 @@ SynchronizerManager::destroy()
   Assert(m_instance != nullptr, "SynchronizerManager was not created!");
 
   if (m_instance->m_connectivity_synchronizer_map.size() > 0) {
+    // LCOV_EXCL_START
     std::stringstream error;
     error << ": some connectivities are still registered\n";
     for (const auto& i_connectivity_synchronizer : m_instance->m_connectivity_synchronizer_map) {
@@ -25,6 +26,7 @@ SynchronizerManager::destroy()
             << '\n';
     }
     throw UnexpectedError(error.str());
+    // LCOV_EXCL_STOP
   }
   delete m_instance;
   m_instance = nullptr;
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 18d5656925cae0cb5d0cba949bebb2ce74ee4092..ddfa977a6ef8ca76df9f27bbdc2bda86395cc67f 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -98,6 +98,9 @@ add_executable (unit_tests
   test_GaussQuadratureDescriptor.cpp
   test_IfProcessor.cpp
   test_IncDecExpressionProcessor.cpp
+  test_IntegrateCellArray.cpp
+  test_IntegrateCellValue.cpp
+  test_IntegrateOnCells.cpp
   test_INodeProcessor.cpp
   test_ItemId.cpp
   test_ItemType.cpp
@@ -171,6 +174,7 @@ add_executable (mpi_unit_tests
   test_RandomEngine.cpp
   test_SubItemValuePerItem.cpp
   test_SubItemArrayPerItem.cpp
+  test_Synchronizer.cpp
   )
 
 add_library(test_Pugs_MeshDataBase
diff --git a/tests/test_CellType.cpp b/tests/test_CellType.cpp
index 0e4a42331cf6606c9d8c0edd60f372676e6bb4c0..17696985874a0a8da3a939a65028cd5af2c9ecbe 100644
--- a/tests/test_CellType.cpp
+++ b/tests/test_CellType.cpp
@@ -13,6 +13,7 @@ TEST_CASE("CellType", "[mesh]")
 
   REQUIRE(name(CellType::Triangle) == "triangle");
   REQUIRE(name(CellType::Quadrangle) == "quadrangle");
+  REQUIRE(name(CellType::Polygon) == "polygon");
 
   REQUIRE(name(CellType::Diamond) == "diamond");
   REQUIRE(name(CellType::Hexahedron) == "hexahedron");
diff --git a/tests/test_DualMeshManager.cpp b/tests/test_DualMeshManager.cpp
index e5c8dbcf83d5296390be98725e3fbca887464f5f..1d05b996ae569ea30e995e5b1c8422b1edb9aeef 100644
--- a/tests/test_DualMeshManager.cpp
+++ b/tests/test_DualMeshManager.cpp
@@ -11,89 +11,149 @@
 
 TEST_CASE("DualMeshManager", "[mesh]")
 {
-  using ConnectivityType = Connectivity<2>;
-  using MeshType         = Mesh<ConnectivityType>;
+  SECTION("same 1D dual connectivities ")
+  {
+    using ConnectivityType = Connectivity<1>;
+    using MeshType         = Mesh<ConnectivityType>;
 
-  std::shared_ptr<const MeshType> mesh = MeshDataBaseForTests::get().hybrid2DMesh();
+    std::shared_ptr<const MeshType> mesh = MeshDataBaseForTests::get().unordered1DMesh();
 
-  SECTION("diamond dual mesh access")
-  {
     std::shared_ptr p_diamond_dual_mesh = DualMeshManager::instance().getDiamondDualMesh(*mesh);
+    std::shared_ptr p_median_dual_mesh  = DualMeshManager::instance().getMedianDualMesh(*mesh);
+    std::shared_ptr p_dual1d_mesh       = DualMeshManager::instance().getDual1DMesh(*mesh);
+
+    // In 1d all these dual meshes are the same
+    REQUIRE(p_dual1d_mesh.get() == p_diamond_dual_mesh.get());
+    REQUIRE(p_dual1d_mesh.get() == p_median_dual_mesh.get());
+  }
+
+  SECTION("2D")
+  {
+    using ConnectivityType = Connectivity<2>;
+    using MeshType         = Mesh<ConnectivityType>;
 
-    const auto ref_counter = p_diamond_dual_mesh.use_count();
+    std::shared_ptr<const MeshType> mesh = MeshDataBaseForTests::get().hybrid2DMesh();
 
+    SECTION("diamond dual mesh access")
     {
-      std::shared_ptr p_diamond_dual_mesh2 = DualMeshManager::instance().getDiamondDualMesh(*mesh);
+      std::shared_ptr p_diamond_dual_mesh = DualMeshManager::instance().getDiamondDualMesh(*mesh);
 
-      REQUIRE(p_diamond_dual_mesh == p_diamond_dual_mesh2);
-      REQUIRE(p_diamond_dual_mesh.use_count() == ref_counter + 1);
-    }
+      const auto ref_counter = p_diamond_dual_mesh.use_count();
 
-    REQUIRE(p_diamond_dual_mesh.use_count() == ref_counter);
+      {
+        std::shared_ptr p_diamond_dual_mesh2 = DualMeshManager::instance().getDiamondDualMesh(*mesh);
 
-    DualMeshManager::instance().deleteMesh(mesh.get());
-    REQUIRE(p_diamond_dual_mesh.use_count() == ref_counter - 1);
+        REQUIRE(p_diamond_dual_mesh == p_diamond_dual_mesh2);
+        REQUIRE(p_diamond_dual_mesh.use_count() == ref_counter + 1);
+      }
 
-    // Can delete mesh from the list again. This means that no
-    // dual mesh associated with it is managed.
-    REQUIRE_NOTHROW(DualMeshManager::instance().deleteMesh(mesh.get()));
-    REQUIRE(p_diamond_dual_mesh.use_count() == ref_counter - 1);
+      REQUIRE(p_diamond_dual_mesh.use_count() == ref_counter);
 
-    // A new dual mesh is built
-    std::shared_ptr p_diamond_dual_mesh_rebuilt = DualMeshManager::instance().getDiamondDualMesh(*mesh);
-    REQUIRE(p_diamond_dual_mesh != p_diamond_dual_mesh_rebuilt);
-    REQUIRE(p_diamond_dual_mesh.get() != p_diamond_dual_mesh_rebuilt.get());
+      DualMeshManager::instance().deleteMesh(mesh.get());
+      REQUIRE(p_diamond_dual_mesh.use_count() == ref_counter - 1);
 
-    // Exactly two references to the dual mesh. One here and
-    // one in the manager.
-    REQUIRE(p_diamond_dual_mesh_rebuilt.use_count() == 2);
-  }
+      // Can delete mesh from the list again. This means that no
+      // dual mesh associated with it is managed.
+      REQUIRE_NOTHROW(DualMeshManager::instance().deleteMesh(mesh.get()));
+      REQUIRE(p_diamond_dual_mesh.use_count() == ref_counter - 1);
 
-  SECTION("median dual mesh access")
-  {
-    std::shared_ptr p_median_dual_mesh = DualMeshManager::instance().getMedianDualMesh(*mesh);
+      // A new dual mesh is built
+      std::shared_ptr p_diamond_dual_mesh_rebuilt = DualMeshManager::instance().getDiamondDualMesh(*mesh);
+      REQUIRE(p_diamond_dual_mesh != p_diamond_dual_mesh_rebuilt);
+      REQUIRE(p_diamond_dual_mesh.get() != p_diamond_dual_mesh_rebuilt.get());
 
-    const auto ref_counter = p_median_dual_mesh.use_count();
+      // Exactly two references to the dual mesh. One here and
+      // one in the manager.
+      REQUIRE(p_diamond_dual_mesh_rebuilt.use_count() == 2);
+    }
 
+    SECTION("median dual mesh access")
     {
-      std::shared_ptr p_median_dual_mesh2 = DualMeshManager::instance().getMedianDualMesh(*mesh);
+      std::shared_ptr p_median_dual_mesh = DualMeshManager::instance().getMedianDualMesh(*mesh);
 
-      REQUIRE(p_median_dual_mesh == p_median_dual_mesh2);
-      REQUIRE(p_median_dual_mesh.use_count() == ref_counter + 1);
-    }
+      const auto ref_counter = p_median_dual_mesh.use_count();
+
+      {
+        std::shared_ptr p_median_dual_mesh2 = DualMeshManager::instance().getMedianDualMesh(*mesh);
+
+        REQUIRE(p_median_dual_mesh == p_median_dual_mesh2);
+        REQUIRE(p_median_dual_mesh.use_count() == ref_counter + 1);
+      }
+
+      REQUIRE(p_median_dual_mesh.use_count() == ref_counter);
+
+      DualMeshManager::instance().deleteMesh(mesh.get());
+      REQUIRE(p_median_dual_mesh.use_count() == ref_counter - 1);
 
-    REQUIRE(p_median_dual_mesh.use_count() == ref_counter);
+      // Can delete mesh from the list again. This means that no
+      // dual mesh associated with it is managed.
+      REQUIRE_NOTHROW(DualMeshManager::instance().deleteMesh(mesh.get()));
+      REQUIRE(p_median_dual_mesh.use_count() == ref_counter - 1);
 
-    DualMeshManager::instance().deleteMesh(mesh.get());
-    REQUIRE(p_median_dual_mesh.use_count() == ref_counter - 1);
+      // A new dual mesh is built
+      std::shared_ptr p_median_dual_mesh_rebuilt = DualMeshManager::instance().getMedianDualMesh(*mesh);
+      REQUIRE(p_median_dual_mesh != p_median_dual_mesh_rebuilt);
+      REQUIRE(p_median_dual_mesh.get() != p_median_dual_mesh_rebuilt.get());
 
-    // Can delete mesh from the list again. This means that no
-    // dual mesh associated with it is managed.
-    REQUIRE_NOTHROW(DualMeshManager::instance().deleteMesh(mesh.get()));
-    REQUIRE(p_median_dual_mesh.use_count() == ref_counter - 1);
+      // Exactly two references to the dual mesh. One here and
+      // one in the manager.
+      REQUIRE(p_median_dual_mesh_rebuilt.use_count() == 2);
+    }
+
+    SECTION("check multiple dual mesh using/freeing")
+    {
+      std::shared_ptr p_median_dual_mesh  = DualMeshManager::instance().getMedianDualMesh(*mesh);
+      std::shared_ptr p_diamond_dual_mesh = DualMeshManager::instance().getDiamondDualMesh(*mesh);
+
+      const auto median_ref_counter  = p_median_dual_mesh.use_count();
+      const auto diamond_ref_counter = p_diamond_dual_mesh.use_count();
 
-    // A new dual mesh is built
-    std::shared_ptr p_median_dual_mesh_rebuilt = DualMeshManager::instance().getMedianDualMesh(*mesh);
-    REQUIRE(p_median_dual_mesh != p_median_dual_mesh_rebuilt);
-    REQUIRE(p_median_dual_mesh.get() != p_median_dual_mesh_rebuilt.get());
+      REQUIRE(p_median_dual_mesh != p_diamond_dual_mesh);
 
-    // Exactly two references to the dual mesh. One here and
-    // one in the manager.
-    REQUIRE(p_median_dual_mesh_rebuilt.use_count() == 2);
+      REQUIRE_NOTHROW(DualMeshManager::instance().deleteMesh(mesh.get()));
+      REQUIRE(p_median_dual_mesh.use_count() == median_ref_counter - 1);
+      REQUIRE(p_diamond_dual_mesh.use_count() == diamond_ref_counter - 1);
+    }
   }
 
-  SECTION("check multiple dual mesh using/freeing")
+  SECTION("3D")
   {
-    std::shared_ptr p_median_dual_mesh  = DualMeshManager::instance().getMedianDualMesh(*mesh);
-    std::shared_ptr p_diamond_dual_mesh = DualMeshManager::instance().getDiamondDualMesh(*mesh);
+    using ConnectivityType = Connectivity<3>;
+    using MeshType         = Mesh<ConnectivityType>;
+
+    std::shared_ptr<const MeshType> mesh = MeshDataBaseForTests::get().hybrid3DMesh();
 
-    const auto median_ref_counter  = p_median_dual_mesh.use_count();
-    const auto diamond_ref_counter = p_diamond_dual_mesh.use_count();
+    SECTION("diamond dual mesh access")
+    {
+      std::shared_ptr p_diamond_dual_mesh = DualMeshManager::instance().getDiamondDualMesh(*mesh);
+
+      const auto ref_counter = p_diamond_dual_mesh.use_count();
+
+      {
+        std::shared_ptr p_diamond_dual_mesh2 = DualMeshManager::instance().getDiamondDualMesh(*mesh);
+
+        REQUIRE(p_diamond_dual_mesh == p_diamond_dual_mesh2);
+        REQUIRE(p_diamond_dual_mesh.use_count() == ref_counter + 1);
+      }
 
-    REQUIRE(p_median_dual_mesh != p_diamond_dual_mesh);
+      REQUIRE(p_diamond_dual_mesh.use_count() == ref_counter);
 
-    REQUIRE_NOTHROW(DualMeshManager::instance().deleteMesh(mesh.get()));
-    REQUIRE(p_median_dual_mesh.use_count() == median_ref_counter - 1);
-    REQUIRE(p_diamond_dual_mesh.use_count() == diamond_ref_counter - 1);
+      DualMeshManager::instance().deleteMesh(mesh.get());
+      REQUIRE(p_diamond_dual_mesh.use_count() == ref_counter - 1);
+
+      // Can delete mesh from the list again. This means that no
+      // dual mesh associated with it is managed.
+      REQUIRE_NOTHROW(DualMeshManager::instance().deleteMesh(mesh.get()));
+      REQUIRE(p_diamond_dual_mesh.use_count() == ref_counter - 1);
+
+      // A new dual mesh is built
+      std::shared_ptr p_diamond_dual_mesh_rebuilt = DualMeshManager::instance().getDiamondDualMesh(*mesh);
+      REQUIRE(p_diamond_dual_mesh != p_diamond_dual_mesh_rebuilt);
+      REQUIRE(p_diamond_dual_mesh.get() != p_diamond_dual_mesh_rebuilt.get());
+
+      // Exactly two references to the dual mesh. One here and
+      // one in the manager.
+      REQUIRE(p_diamond_dual_mesh_rebuilt.use_count() == 2);
+    }
   }
 }
diff --git a/tests/test_DualMeshType.cpp b/tests/test_DualMeshType.cpp
index 531b1ae803edc7ca59c4aef5b6fe2be6e500741f..0c5b932bacdad06e67500cfd7a0045c47f96cb20 100644
--- a/tests/test_DualMeshType.cpp
+++ b/tests/test_DualMeshType.cpp
@@ -8,6 +8,7 @@
 TEST_CASE("DualMeshType", "[mesh]")
 {
   REQUIRE(name(DualMeshType::Diamond) == "diamond");
+  REQUIRE(name(DualMeshType::Dual1D) == "dual 1d");
   REQUIRE(name(DualMeshType::Median) == "median");
   REQUIRE_THROWS_WITH(name(DualMeshType{-1}), "unexpected error: unexpected dual mesh type");
 }
diff --git a/tests/test_IntegrateCellArray.cpp b/tests/test_IntegrateCellArray.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b6b8db8c8ca7bc0aa895f8be6b7d5f8dbccd89ff
--- /dev/null
+++ b/tests/test_IntegrateCellArray.cpp
@@ -0,0 +1,626 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <language/ast/ASTBuilder.hpp>
+#include <language/ast/ASTModulesImporter.hpp>
+#include <language/ast/ASTNodeDataTypeBuilder.hpp>
+#include <language/ast/ASTNodeExpressionBuilder.hpp>
+#include <language/ast/ASTNodeFunctionEvaluationExpressionBuilder.hpp>
+#include <language/ast/ASTNodeFunctionExpressionBuilder.hpp>
+#include <language/ast/ASTNodeTypeCleaner.hpp>
+#include <language/ast/ASTSymbolTableBuilder.hpp>
+#include <language/utils/PugsFunctionAdapter.hpp>
+#include <language/utils/SymbolTable.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/DualMeshManager.hpp>
+#include <mesh/Mesh.hpp>
+#include <scheme/CellIntegrator.hpp>
+
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/GaussLobattoQuadratureDescriptor.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
+
+#include <language/utils/IntegrateCellArray.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("IntegrateCellArray", "[language]")
+{
+  SECTION("integrate on all cells")
+  {
+    auto same_item_integral = [](auto f, auto g) -> bool {
+      using ItemIdType = typename decltype(f)::index_type;
+      for (ItemIdType item_id = 0; item_id < f.numberOfItems(); ++item_id) {
+        for (size_t i = 0; i < f.sizeOfArrays(); ++i) {
+          if (f[item_id][i] != g[item_id][i]) {
+            return false;
+          }
+        }
+      }
+
+      return true;
+    };
+
+    SECTION("1D")
+    {
+      constexpr size_t Dimension = 1;
+      auto quadrature_descriptor = GaussQuadratureDescriptor(3);
+
+      std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_1d = named_mesh.mesh();
+
+          std::string_view data = R"(
+import math;
+import math;
+let f: R^1 -> R, x -> 2*x[0] + 2;
+let g: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+)";
+          TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+          auto ast = ASTBuilder::build(input);
+
+          ASTModulesImporter{*ast};
+          ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+          ASTSymbolTableBuilder{*ast};
+          ASTNodeDataTypeBuilder{*ast};
+
+          ASTNodeTypeCleaner<language::var_declaration>{*ast};
+          ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+          ASTNodeExpressionBuilder{*ast};
+
+          std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+          TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+          position.byte = data.size();   // ensure that variables are declared at this point
+
+          std::vector<FunctionSymbolId> function_symbol_id_list;
+
+          {
+            auto [i_symbol, found] = symbol_table->find("f", position);
+            REQUIRE(found);
+            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+            function_symbol_id_list.push_back(
+              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+          }
+
+          {
+            auto [i_symbol, found] = symbol_table->find("g", position);
+            REQUIRE(found);
+            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+            function_symbol_id_list.push_back(
+              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+          }
+
+          CellArray<double> cell_integral_array{mesh_1d->connectivity(), 2};
+
+          {
+            CellValue<double> cell_f_integral{mesh_1d->connectivity()};
+            auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * x[0] + 2; };
+            CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_1d, cell_f_integral);
+
+            parallel_for(
+              mesh_1d->numberOfCells(),
+              PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][0] = cell_f_integral[cell_id]; });
+
+            CellValue<double> cell_g_integral{mesh_1d->connectivity()};
+            auto g = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) + 3; };
+            CellIntegrator::integrateTo(g, quadrature_descriptor, *mesh_1d, cell_g_integral);
+
+            parallel_for(
+              mesh_1d->numberOfCells(),
+              PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][1] = cell_g_integral[cell_id]; });
+          }
+
+          CellArray<double> integrate_array =
+            IntegrateCellArray<double(TinyVector<Dimension>)>::integrate(function_symbol_id_list, quadrature_descriptor,
+                                                                         *mesh_1d);
+
+          REQUIRE(same_item_integral(cell_integral_array, integrate_array));
+        }
+      }
+    }
+
+    SECTION("2D")
+    {
+      constexpr size_t Dimension = 2;
+      auto quadrature_descriptor = GaussLobattoQuadratureDescriptor(3);
+
+      std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_2d = named_mesh.mesh();
+
+          std::string_view data = R"(
+import math;
+let f: R^2 -> R, x -> 2*x[0] + 3*x[1] + 2;
+let g: R^2 -> R, x -> 2*exp(x[0])*sin(x[1])+3;
+)";
+
+          TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+          auto ast = ASTBuilder::build(input);
+
+          ASTModulesImporter{*ast};
+          ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+          ASTSymbolTableBuilder{*ast};
+          ASTNodeDataTypeBuilder{*ast};
+
+          ASTNodeTypeCleaner<language::var_declaration>{*ast};
+          ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+          ASTNodeExpressionBuilder{*ast};
+
+          std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+          TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+          position.byte = data.size();   // ensure that variables are declared at this point
+
+          std::vector<FunctionSymbolId> function_symbol_id_list;
+
+          {
+            auto [i_symbol, found] = symbol_table->find("f", position);
+            REQUIRE(found);
+            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+            function_symbol_id_list.push_back(
+              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+          }
+
+          {
+            auto [i_symbol, found] = symbol_table->find("g", position);
+            REQUIRE(found);
+            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+            function_symbol_id_list.push_back(
+              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+          }
+
+          CellArray<double> cell_integral_array{mesh_2d->connectivity(), 2};
+
+          {
+            CellValue<double> cell_f_integral{mesh_2d->connectivity()};
+            auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * x[0] + 3 * x[1] + 2; };
+            CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_2d, cell_f_integral);
+
+            parallel_for(
+              mesh_2d->numberOfCells(),
+              PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][0] = cell_f_integral[cell_id]; });
+
+            CellValue<double> cell_g_integral{mesh_2d->connectivity()};
+            auto g = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) + 3; };
+            CellIntegrator::integrateTo(g, quadrature_descriptor, *mesh_2d, cell_g_integral);
+
+            parallel_for(
+              mesh_2d->numberOfCells(),
+              PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][1] = cell_g_integral[cell_id]; });
+          }
+
+          CellArray<double> integrate_array =
+            IntegrateCellArray<double(TinyVector<Dimension>)>::integrate(function_symbol_id_list, quadrature_descriptor,
+                                                                         *mesh_2d);
+
+          REQUIRE(same_item_integral(cell_integral_array, integrate_array));
+        }
+      }
+    }
+
+    SECTION("3D")
+    {
+      constexpr size_t Dimension = 3;
+      auto quadrature_descriptor = GaussLegendreQuadratureDescriptor(3);
+
+      using NamedMesh = MeshDataBaseForTests::NamedMesh<Dimension>;
+
+      std::vector<NamedMesh> mesh_list = [] {
+        std::vector<NamedMesh> extended_mesh_list;
+        std::array mesh_array = MeshDataBaseForTests::get().all3DMeshes();
+        for (size_t i = 0; i < mesh_array.size(); ++i) {
+          extended_mesh_list.push_back(MeshDataBaseForTests::get().all3DMeshes()[i]);
+        }
+        extended_mesh_list.push_back(NamedMesh("diamond dual", DualMeshManager::instance().getDiamondDualMesh(
+                                                                 *MeshDataBaseForTests::get().hybrid3DMesh())));
+        return extended_mesh_list;
+      }();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_3d = named_mesh.mesh();
+
+          std::string_view data = R"(
+import math;
+let f: R^3 -> R, x -> 2 * x[0] + 3 * x[1] + 2 * x[2] - 1;
+let g: R^3 -> R, x -> 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
+)";
+
+          TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+          auto ast = ASTBuilder::build(input);
+
+          ASTModulesImporter{*ast};
+          ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+          ASTSymbolTableBuilder{*ast};
+          ASTNodeDataTypeBuilder{*ast};
+
+          ASTNodeTypeCleaner<language::var_declaration>{*ast};
+          ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+          ASTNodeExpressionBuilder{*ast};
+
+          std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+          TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+          position.byte = data.size();   // ensure that variables are declared at this point
+
+          std::vector<FunctionSymbolId> function_symbol_id_list;
+
+          {
+            auto [i_symbol, found] = symbol_table->find("f", position);
+            REQUIRE(found);
+            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+            function_symbol_id_list.push_back(
+              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+          }
+
+          {
+            auto [i_symbol, found] = symbol_table->find("g", position);
+            REQUIRE(found);
+            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+            function_symbol_id_list.push_back(
+              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+          }
+
+          CellArray<double> cell_integral_array{mesh_3d->connectivity(), 2};
+
+          {
+            CellValue<double> cell_f_integral{mesh_3d->connectivity()};
+            auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * x[0] + 3 * x[1] + 2 * x[2] - 1; };
+            CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_3d, cell_f_integral);
+
+            parallel_for(
+              mesh_3d->numberOfCells(),
+              PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][0] = cell_f_integral[cell_id]; });
+
+            CellValue<double> cell_g_integral{mesh_3d->connectivity()};
+            auto g = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) * x[2] + 3; };
+            CellIntegrator::integrateTo(g, quadrature_descriptor, *mesh_3d, cell_g_integral);
+
+            parallel_for(
+              mesh_3d->numberOfCells(),
+              PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][1] = cell_g_integral[cell_id]; });
+          }
+
+          CellArray<double> integrate_array =
+            IntegrateCellArray<double(TinyVector<Dimension>)>::integrate(function_symbol_id_list, quadrature_descriptor,
+                                                                         *mesh_3d);
+
+          REQUIRE(same_item_integral(cell_integral_array, integrate_array));
+        }
+      }
+    }
+  }
+
+  SECTION("integrate on cell list")
+  {
+    auto same_item_integral = [](auto f, auto g) -> bool {
+      using ItemIdType = typename decltype(f)::index_type;
+      for (ItemIdType item_id = 0; item_id < f.numberOfRows(); ++item_id) {
+        for (size_t i = 0; i < f.numberOfColumns(); ++i) {
+          if (f[item_id][i] != g[item_id][i]) {
+            return false;
+          }
+        }
+      }
+
+      return true;
+    };
+
+    SECTION("1D")
+    {
+      constexpr size_t Dimension = 1;
+      auto quadrature_descriptor = GaussLegendreQuadratureDescriptor(3);
+
+      std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_1d = named_mesh.mesh();
+          Array<CellId> cell_list{mesh_1d->numberOfCells() / 2 + mesh_1d->numberOfCells() % 2};
+
+          {
+            size_t k = 0;
+            for (CellId cell_id = 0; cell_id < mesh_1d->numberOfCells(); ++(++cell_id), ++k) {
+              cell_list[k] = cell_id;
+            }
+
+            REQUIRE(k == cell_list.size());
+          }
+
+          std::string_view data = R"(
+import math;
+let f: R^1 -> R, x -> 2*x[0] + 2;
+let g: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+)";
+          TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+          auto ast = ASTBuilder::build(input);
+
+          ASTModulesImporter{*ast};
+          ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+          ASTSymbolTableBuilder{*ast};
+          ASTNodeDataTypeBuilder{*ast};
+
+          ASTNodeTypeCleaner<language::var_declaration>{*ast};
+          ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+          ASTNodeExpressionBuilder{*ast};
+
+          std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+          TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+          position.byte = data.size();   // ensure that variables are declared at this point
+
+          std::vector<FunctionSymbolId> function_symbol_id_list;
+
+          {
+            auto [i_symbol, found] = symbol_table->find("f", position);
+            REQUIRE(found);
+            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+            function_symbol_id_list.push_back(
+              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+          }
+
+          {
+            auto [i_symbol, found] = symbol_table->find("g", position);
+            REQUIRE(found);
+            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+            function_symbol_id_list.push_back(
+              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+          }
+
+          Table<double> cell_integral_array{cell_list.size(), 2};
+
+          {
+            auto f                        = [](const TinyVector<Dimension>& x) -> double { return 2 * x[0] + 2; };
+            Array<double> cell_f_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_1d, cell_list);
+
+            parallel_for(
+              cell_integral_array.numberOfRows(),
+              PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][0] = cell_f_integral[cell_id]; });
+
+            auto g                        = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) + 3; };
+            Array<double> cell_g_integral = CellIntegrator::integrate(g, quadrature_descriptor, *mesh_1d, cell_list);
+
+            parallel_for(
+              cell_integral_array.numberOfRows(),
+              PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][1] = cell_g_integral[cell_id]; });
+          }
+
+          Table<const double> integrate_value =
+            IntegrateCellArray<double(TinyVector<Dimension>)>::integrate(function_symbol_id_list, quadrature_descriptor,
+                                                                         *mesh_1d, cell_list);
+
+          REQUIRE(same_item_integral(cell_integral_array, integrate_value));
+        }
+      }
+    }
+
+    SECTION("2D")
+    {
+      constexpr size_t Dimension = 2;
+      auto quadrature_descriptor = GaussLegendreQuadratureDescriptor(3);
+
+      std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_2d = named_mesh.mesh();
+
+          Array<CellId> cell_list{mesh_2d->numberOfCells() / 2 + mesh_2d->numberOfCells() % 2};
+
+          {
+            size_t k = 0;
+            for (CellId cell_id = 0; cell_id < mesh_2d->numberOfCells(); ++(++cell_id), ++k) {
+              cell_list[k] = cell_id;
+            }
+
+            REQUIRE(k == cell_list.size());
+          }
+
+          std::string_view data = R"(
+import math;
+let f: R^2 -> R, x -> 2*x[0] + 3*x[1] + 2;
+let g: R^2 -> R, x -> 2*exp(x[0])*sin(x[1])+3;
+)";
+
+          TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+          auto ast = ASTBuilder::build(input);
+
+          ASTModulesImporter{*ast};
+          ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+          ASTSymbolTableBuilder{*ast};
+          ASTNodeDataTypeBuilder{*ast};
+
+          ASTNodeTypeCleaner<language::var_declaration>{*ast};
+          ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+          ASTNodeExpressionBuilder{*ast};
+
+          std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+          TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+          position.byte = data.size();   // ensure that variables are declared at this point
+
+          std::vector<FunctionSymbolId> function_symbol_id_list;
+
+          {
+            auto [i_symbol, found] = symbol_table->find("f", position);
+            REQUIRE(found);
+            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+            function_symbol_id_list.push_back(
+              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+          }
+
+          {
+            auto [i_symbol, found] = symbol_table->find("g", position);
+            REQUIRE(found);
+            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+            function_symbol_id_list.push_back(
+              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+          }
+
+          Table<double> cell_integral_array{cell_list.size(), 2};
+
+          {
+            auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * x[0] + 3 * x[1] + 2; };
+            Array<double> cell_f_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_2d, cell_list);
+
+            parallel_for(
+              cell_integral_array.numberOfRows(),
+              PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][0] = cell_f_integral[cell_id]; });
+
+            auto g = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) + 3; };
+            Array<double> cell_g_integral = CellIntegrator::integrate(g, quadrature_descriptor, *mesh_2d, cell_list);
+
+            parallel_for(
+              cell_integral_array.numberOfRows(),
+              PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][1] = cell_g_integral[cell_id]; });
+          }
+
+          Table<const double> integrate_value =
+            IntegrateCellArray<double(TinyVector<Dimension>)>::integrate(function_symbol_id_list, quadrature_descriptor,
+                                                                         *mesh_2d, cell_list);
+
+          REQUIRE(same_item_integral(cell_integral_array, integrate_value));
+        }
+      }
+    }
+
+    SECTION("3D")
+    {
+      constexpr size_t Dimension = 3;
+      auto quadrature_descriptor = GaussQuadratureDescriptor(3);
+
+      using NamedMesh = MeshDataBaseForTests::NamedMesh<Dimension>;
+
+      std::vector<NamedMesh> mesh_list = [] {
+        std::vector<NamedMesh> extended_mesh_list;
+        std::array mesh_array = MeshDataBaseForTests::get().all3DMeshes();
+        for (size_t i = 0; i < mesh_array.size(); ++i) {
+          extended_mesh_list.push_back(MeshDataBaseForTests::get().all3DMeshes()[i]);
+        }
+        extended_mesh_list.push_back(NamedMesh("diamond dual", DualMeshManager::instance().getDiamondDualMesh(
+                                                                 *MeshDataBaseForTests::get().hybrid3DMesh())));
+        return extended_mesh_list;
+      }();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_3d = named_mesh.mesh();
+
+          Array<CellId> cell_list{mesh_3d->numberOfCells() / 2 + mesh_3d->numberOfCells() % 2};
+
+          {
+            size_t k = 0;
+            for (CellId cell_id = 0; cell_id < mesh_3d->numberOfCells(); ++(++cell_id), ++k) {
+              cell_list[k] = cell_id;
+            }
+
+            REQUIRE(k == cell_list.size());
+          }
+
+          std::string_view data = R"(
+import math;
+let f: R^3 -> R, x -> 2 * x[0] + 3 * x[1] + 2 * x[2] - 1;
+let g: R^3 -> R, x -> 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
+)";
+
+          TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+          auto ast = ASTBuilder::build(input);
+
+          ASTModulesImporter{*ast};
+          ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+          ASTSymbolTableBuilder{*ast};
+          ASTNodeDataTypeBuilder{*ast};
+
+          ASTNodeTypeCleaner<language::var_declaration>{*ast};
+          ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+          ASTNodeExpressionBuilder{*ast};
+
+          std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+          TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+          position.byte = data.size();   // ensure that variables are declared at this point
+
+          std::vector<FunctionSymbolId> function_symbol_id_list;
+
+          {
+            auto [i_symbol, found] = symbol_table->find("f", position);
+            REQUIRE(found);
+            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+            function_symbol_id_list.push_back(
+              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+          }
+
+          {
+            auto [i_symbol, found] = symbol_table->find("g", position);
+            REQUIRE(found);
+            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+            function_symbol_id_list.push_back(
+              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+          }
+
+          Table<double> cell_integral_array{cell_list.size(), 2};
+
+          {
+            auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * x[0] + 3 * x[1] + 2 * x[2] - 1; };
+            Array<double> cell_f_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_3d, cell_list);
+
+            parallel_for(
+              cell_integral_array.numberOfRows(),
+              PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][0] = cell_f_integral[cell_id]; });
+
+            auto g = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) * x[2] + 3; };
+            Array<double> cell_g_integral = CellIntegrator::integrate(g, quadrature_descriptor, *mesh_3d, cell_list);
+
+            parallel_for(
+              cell_integral_array.numberOfRows(),
+              PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][1] = cell_g_integral[cell_id]; });
+          }
+
+          Table<const double> integrate_value =
+            IntegrateCellArray<double(TinyVector<Dimension>)>::integrate(function_symbol_id_list, quadrature_descriptor,
+                                                                         *mesh_3d, cell_list);
+
+          REQUIRE(same_item_integral(cell_integral_array, integrate_value));
+        }
+      }
+    }
+  }
+}
diff --git a/tests/test_IntegrateCellValue.cpp b/tests/test_IntegrateCellValue.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..df35a215bea28b236cd2bcd7b243e2f759c4f40b
--- /dev/null
+++ b/tests/test_IntegrateCellValue.cpp
@@ -0,0 +1,451 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <language/ast/ASTBuilder.hpp>
+#include <language/ast/ASTModulesImporter.hpp>
+#include <language/ast/ASTNodeDataTypeBuilder.hpp>
+#include <language/ast/ASTNodeExpressionBuilder.hpp>
+#include <language/ast/ASTNodeFunctionEvaluationExpressionBuilder.hpp>
+#include <language/ast/ASTNodeFunctionExpressionBuilder.hpp>
+#include <language/ast/ASTNodeTypeCleaner.hpp>
+#include <language/ast/ASTSymbolTableBuilder.hpp>
+#include <language/utils/PugsFunctionAdapter.hpp>
+#include <language/utils/SymbolTable.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/DualMeshManager.hpp>
+#include <mesh/Mesh.hpp>
+#include <scheme/CellIntegrator.hpp>
+
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/GaussLobattoQuadratureDescriptor.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
+
+#include <language/utils/IntegrateCellValue.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("IntegrateCellValue", "[language]")
+{
+  SECTION("integrate on all cells")
+  {
+    auto same_item_integral = [](auto f, auto g) -> bool {
+      using ItemIdType = typename decltype(f)::index_type;
+      for (ItemIdType item_id = 0; item_id < f.numberOfItems(); ++item_id) {
+        if (f[item_id] != g[item_id]) {
+          return false;
+        }
+      }
+
+      return true;
+    };
+
+    SECTION("1D")
+    {
+      constexpr size_t Dimension = 1;
+      auto quadrature_descriptor = GaussQuadratureDescriptor(3);
+
+      std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_1d = named_mesh.mesh();
+
+          std::string_view data = R"(
+import math;
+let R2x2_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0]);
+)";
+          TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+          auto ast = ASTBuilder::build(input);
+
+          ASTModulesImporter{*ast};
+          ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+          ASTSymbolTableBuilder{*ast};
+          ASTNodeDataTypeBuilder{*ast};
+
+          ASTNodeTypeCleaner<language::var_declaration>{*ast};
+          ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+          ASTNodeExpressionBuilder{*ast};
+
+          std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+          TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+          position.byte = data.size();   // ensure that variables are declared at this point
+
+          using R2x2             = TinyMatrix<2>;
+          auto [i_symbol, found] = symbol_table->find("R2x2_1d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<R2x2> cell_integral{mesh_1d->connectivity()};
+          auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+            return R2x2{2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0]};
+          };
+          CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_1d, cell_integral);
+
+          CellValue<R2x2> integrate_value =
+            IntegrateCellValue<R2x2(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_1d);
+
+          REQUIRE(same_item_integral(cell_integral, integrate_value));
+        }
+      }
+    }
+
+    SECTION("2D")
+    {
+      constexpr size_t Dimension = 2;
+      auto quadrature_descriptor = GaussLobattoQuadratureDescriptor(3);
+
+      std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_2d = named_mesh.mesh();
+
+          std::string_view data = R"(
+import math;
+let R3_2d: R^2 -> R^3, x -> (2*exp(x[0])*sin(x[1])+3, x[0]-2*x[1], 3);
+)";
+
+          TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+          auto ast = ASTBuilder::build(input);
+
+          ASTModulesImporter{*ast};
+          ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+          ASTSymbolTableBuilder{*ast};
+          ASTNodeDataTypeBuilder{*ast};
+
+          ASTNodeTypeCleaner<language::var_declaration>{*ast};
+          ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+          ASTNodeExpressionBuilder{*ast};
+
+          std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+          TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+          position.byte = data.size();   // ensure that variables are declared at this point
+
+          using R3               = TinyVector<3>;
+          auto [i_symbol, found] = symbol_table->find("R3_2d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<R3> cell_integral{mesh_2d->connectivity()};
+          auto f = [](const TinyVector<Dimension>& x) -> R3 {
+            return R3{2 * exp(x[0]) * sin(x[1]) + 3, x[0] - 2 * x[1], 3};
+          };
+          CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_2d, cell_integral);
+
+          CellValue<R3> integrate_value =
+            IntegrateCellValue<R3(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                     *mesh_2d);
+
+          REQUIRE(same_item_integral(cell_integral, integrate_value));
+        }
+      }
+    }
+
+    SECTION("3D")
+    {
+      constexpr size_t Dimension = 3;
+      auto quadrature_descriptor = GaussLegendreQuadratureDescriptor(3);
+
+      using NamedMesh = MeshDataBaseForTests::NamedMesh<Dimension>;
+
+      std::vector<NamedMesh> mesh_list = [] {
+        std::vector<NamedMesh> extended_mesh_list;
+        std::array mesh_array = MeshDataBaseForTests::get().all3DMeshes();
+        for (size_t i = 0; i < mesh_array.size(); ++i) {
+          extended_mesh_list.push_back(MeshDataBaseForTests::get().all3DMeshes()[i]);
+        }
+        extended_mesh_list.push_back(NamedMesh("diamond dual", DualMeshManager::instance().getDiamondDualMesh(
+                                                                 *MeshDataBaseForTests::get().hybrid3DMesh())));
+        return extended_mesh_list;
+      }();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_3d = named_mesh.mesh();
+
+          std::string_view data = R"(
+import math;
+let scalar_3d: R^3 -> R, x -> 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
+)";
+
+          TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+          auto ast = ASTBuilder::build(input);
+
+          ASTModulesImporter{*ast};
+          ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+          ASTSymbolTableBuilder{*ast};
+          ASTNodeDataTypeBuilder{*ast};
+
+          ASTNodeTypeCleaner<language::var_declaration>{*ast};
+          ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+          ASTNodeExpressionBuilder{*ast};
+
+          std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+          TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+          position.byte = data.size();   // ensure that variables are declared at this point
+
+          auto [i_symbol, found] = symbol_table->find("scalar_3d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<double> cell_integral{mesh_3d->connectivity()};
+          auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) * x[2] + 3; };
+          CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_3d, cell_integral);
+
+          CellValue<double> integrate_value =
+            IntegrateCellValue<double(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_3d);
+
+          REQUIRE(same_item_integral(cell_integral, integrate_value));
+        }
+      }
+    }
+  }
+
+  SECTION("integrate on cell list")
+  {
+    auto same_item_integral = [](auto f, auto g) -> bool {
+      using ItemIdType = typename decltype(g)::index_type;
+      for (ItemIdType item_id = 0; item_id < f.size(); ++item_id) {
+        if (f[item_id] != g[item_id]) {
+          return false;
+        }
+      }
+
+      return true;
+    };
+
+    SECTION("1D")
+    {
+      constexpr size_t Dimension = 1;
+      auto quadrature_descriptor = GaussLegendreQuadratureDescriptor(3);
+
+      std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_1d = named_mesh.mesh();
+          Array<CellId> cell_list{mesh_1d->numberOfCells() / 2 + mesh_1d->numberOfCells() % 2};
+
+          {
+            size_t k = 0;
+            for (CellId cell_id = 0; cell_id < mesh_1d->numberOfCells(); ++(++cell_id), ++k) {
+              cell_list[k] = cell_id;
+            }
+
+            REQUIRE(k == cell_list.size());
+          }
+
+          std::string_view data = R"(
+import math;
+let scalar_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+)";
+          TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+          auto ast = ASTBuilder::build(input);
+
+          ASTModulesImporter{*ast};
+          ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+          ASTSymbolTableBuilder{*ast};
+          ASTNodeDataTypeBuilder{*ast};
+
+          ASTNodeTypeCleaner<language::var_declaration>{*ast};
+          ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+          ASTNodeExpressionBuilder{*ast};
+
+          std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+          TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+          position.byte = data.size();   // ensure that variables are declared at this point
+
+          auto [i_symbol, found] = symbol_table->find("scalar_1d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * std::exp(x[0]) + 3; };
+
+          Array<const double> cell_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_1d, cell_list);
+          Array<const double> integrate_value =
+            IntegrateCellValue<double(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_1d, cell_list);
+
+          REQUIRE(same_item_integral(cell_integral, integrate_value));
+        }
+      }
+    }
+
+    SECTION("2D")
+    {
+      constexpr size_t Dimension = 2;
+      auto quadrature_descriptor = GaussLegendreQuadratureDescriptor(3);
+
+      std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_2d = named_mesh.mesh();
+
+          Array<CellId> cell_list{mesh_2d->numberOfCells() / 2 + mesh_2d->numberOfCells() % 2};
+
+          {
+            size_t k = 0;
+            for (CellId cell_id = 0; cell_id < mesh_2d->numberOfCells(); ++(++cell_id), ++k) {
+              cell_list[k] = cell_id;
+            }
+
+            REQUIRE(k == cell_list.size());
+          }
+
+          std::string_view data = R"(
+import math;
+let R3_2d: R^2 -> R^3, x -> (2*exp(x[0])*sin(x[1])+3, x[0]-2*x[1], 3);
+)";
+
+          TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+          auto ast = ASTBuilder::build(input);
+
+          ASTModulesImporter{*ast};
+          ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+          ASTSymbolTableBuilder{*ast};
+          ASTNodeDataTypeBuilder{*ast};
+
+          ASTNodeTypeCleaner<language::var_declaration>{*ast};
+          ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+          ASTNodeExpressionBuilder{*ast};
+
+          std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+          TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+          position.byte = data.size();   // ensure that variables are declared at this point
+
+          using R3               = TinyVector<3>;
+          auto [i_symbol, found] = symbol_table->find("R3_2d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          auto f = [](const TinyVector<Dimension>& x) -> R3 {
+            return R3{2 * exp(x[0]) * sin(x[1]) + 3, x[0] - 2 * x[1], 3};
+          };
+
+          Array<const R3> cell_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_2d, cell_list);
+          Array<const R3> integrate_value =
+            IntegrateCellValue<R3(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                     *mesh_2d, cell_list);
+
+          REQUIRE(same_item_integral(cell_integral, integrate_value));
+        }
+      }
+    }
+
+    SECTION("3D")
+    {
+      constexpr size_t Dimension = 3;
+      auto quadrature_descriptor = GaussQuadratureDescriptor(3);
+
+      using NamedMesh = MeshDataBaseForTests::NamedMesh<Dimension>;
+
+      std::vector<NamedMesh> mesh_list = [] {
+        std::vector<NamedMesh> extended_mesh_list;
+        std::array mesh_array = MeshDataBaseForTests::get().all3DMeshes();
+        for (size_t i = 0; i < mesh_array.size(); ++i) {
+          extended_mesh_list.push_back(MeshDataBaseForTests::get().all3DMeshes()[i]);
+        }
+        extended_mesh_list.push_back(NamedMesh("diamond dual", DualMeshManager::instance().getDiamondDualMesh(
+                                                                 *MeshDataBaseForTests::get().hybrid3DMesh())));
+        return extended_mesh_list;
+      }();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_3d = named_mesh.mesh();
+
+          Array<CellId> cell_list{mesh_3d->numberOfCells() / 2 + mesh_3d->numberOfCells() % 2};
+
+          {
+            size_t k = 0;
+            for (CellId cell_id = 0; cell_id < mesh_3d->numberOfCells(); ++(++cell_id), ++k) {
+              cell_list[k] = cell_id;
+            }
+
+            REQUIRE(k == cell_list.size());
+          }
+
+          std::string_view data = R"(
+import math;
+let R2x2_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3, x[0] * x[1] * x[2]);
+)";
+
+          TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+          auto ast = ASTBuilder::build(input);
+
+          ASTModulesImporter{*ast};
+          ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+          ASTSymbolTableBuilder{*ast};
+          ASTNodeDataTypeBuilder{*ast};
+
+          ASTNodeTypeCleaner<language::var_declaration>{*ast};
+          ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+          ASTNodeExpressionBuilder{*ast};
+
+          std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+          TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+          position.byte = data.size();   // ensure that variables are declared at this point
+
+          using R2x2             = TinyMatrix<2>;
+          auto [i_symbol, found] = symbol_table->find("R2x2_3d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+            return R2x2{2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3, x[0] * x[1] * x[2]};
+          };
+
+          Array<const R2x2> cell_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_3d, cell_list);
+
+          Array<R2x2> integrate_value =
+            IntegrateCellValue<R2x2(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_3d, cell_list);
+
+          REQUIRE(same_item_integral(cell_integral, integrate_value));
+        }
+      }
+    }
+  }
+}
diff --git a/tests/test_IntegrateCellValue.hpp b/tests/test_IntegrateCellValue.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f9a9e32a9a0bb69e61626eab69334c3248ee2b74
--- /dev/null
+++ b/tests/test_IntegrateCellValue.hpp
@@ -0,0 +1 @@
+i_f_symboli_f_symboli_f_symbol
diff --git a/tests/test_IntegrateOnCells.cpp b/tests/test_IntegrateOnCells.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cf2d770265fd3dbf1812b14606aa7fadda972d05
--- /dev/null
+++ b/tests/test_IntegrateOnCells.cpp
@@ -0,0 +1,2118 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <language/ast/ASTBuilder.hpp>
+#include <language/ast/ASTModulesImporter.hpp>
+#include <language/ast/ASTNodeDataTypeBuilder.hpp>
+#include <language/ast/ASTNodeExpressionBuilder.hpp>
+#include <language/ast/ASTNodeFunctionEvaluationExpressionBuilder.hpp>
+#include <language/ast/ASTNodeFunctionExpressionBuilder.hpp>
+#include <language/ast/ASTNodeTypeCleaner.hpp>
+#include <language/ast/ASTSymbolTableBuilder.hpp>
+#include <language/utils/PugsFunctionAdapter.hpp>
+#include <language/utils/SymbolTable.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/DualMeshManager.hpp>
+#include <mesh/Mesh.hpp>
+#include <scheme/CellIntegrator.hpp>
+
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/GaussLobattoQuadratureDescriptor.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <language/utils/IntegrateOnCells.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("IntegrateOnCells", "[language]")
+{
+  SECTION("Gauss quadrature")
+  {
+    auto quadrature_descriptor = GaussQuadratureDescriptor(3);
+
+    SECTION("integrate on all cells")
+    {
+      auto same_item_integral = [](auto f, auto g) -> bool {
+        using ItemIdType = typename decltype(f)::index_type;
+        for (ItemIdType item_id = 0; item_id < f.numberOfItems(); ++item_id) {
+          if (f[item_id] != g[item_id]) {
+            return false;
+          }
+        }
+
+        return true;
+      };
+
+      SECTION("1D")
+      {
+        constexpr size_t Dimension = 1;
+
+        std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+        for (auto named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_1d = named_mesh.mesh();
+
+            std::string_view data = R"(
+import math;
+let scalar_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+let R3_1d: R^1 -> R^3, x -> (2 * exp(x[0]) + 3, x[0] - 2, 3);
+let R2x2_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0]);
+)";
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            SECTION("scalar 1d")
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<double> cell_integral{mesh_1d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * std::exp(x[0]) + 3; };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_1d, cell_integral);
+
+              Array<double> integrate_value(mesh_1d->numberOfCells());
+              IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_1d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("vector 1d")
+            {
+              using R3               = TinyVector<3>;
+              auto [i_symbol, found] = symbol_table->find("R3_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<R3> cell_integral{mesh_1d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> R3 { return R3{2 * exp(x[0]) + 3, x[0] - 2, 3}; };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_1d, cell_integral);
+
+              Array<R3> integrate_value(mesh_1d->numberOfCells());
+              IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_1d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("matrix 1d")
+            {
+              using R2x2             = TinyMatrix<2>;
+              auto [i_symbol, found] = symbol_table->find("R2x2_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<R2x2> cell_integral{mesh_1d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+                return R2x2{2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0]};
+              };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_1d, cell_integral);
+
+              Array<R2x2> integrate_value(mesh_1d->numberOfCells());
+              IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_1d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+          }
+        }
+      }
+
+      SECTION("2D")
+      {
+        constexpr size_t Dimension = 2;
+
+        std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+        for (auto named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_2d = named_mesh.mesh();
+
+            std::string_view data = R"(
+import math;
+let scalar_2d: R^2 -> R, x -> 2*exp(x[0])*sin(x[1])+3;
+let R3_2d: R^2 -> R^3, x -> (2*exp(x[0])*sin(x[1])+3, x[0]-2*x[1], 3);
+let R2x2_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*x[1]), 3, x[0]*x[1]);
+)";
+
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            SECTION("scalar 2d")
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<double> cell_integral{mesh_2d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) + 3; };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_2d, cell_integral);
+
+              Array<double> integrate_value(mesh_2d->numberOfCells());
+              IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_2d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("vector 2d")
+            {
+              using R3               = TinyVector<3>;
+              auto [i_symbol, found] = symbol_table->find("R3_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<R3> cell_integral{mesh_2d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> R3 {
+                return R3{2 * exp(x[0]) * sin(x[1]) + 3, x[0] - 2 * x[1], 3};
+              };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_2d, cell_integral);
+
+              Array<R3> integrate_value(mesh_2d->numberOfCells());
+              IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_2d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("matrix 2d")
+            {
+              using R2x2             = TinyMatrix<2>;
+              auto [i_symbol, found] = symbol_table->find("R2x2_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<R2x2> cell_integral{mesh_2d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+                return R2x2{2 * exp(x[0]) * sin(x[1]) + 3, sin(x[0] - 2 * x[1]), 3, x[0] * x[1]};
+              };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_2d, cell_integral);
+
+              Array<R2x2> integrate_value(mesh_2d->numberOfCells());
+              IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_2d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+          }
+        }
+      }
+
+      SECTION("3D")
+      {
+        constexpr size_t Dimension = 3;
+        using NamedMesh            = MeshDataBaseForTests::NamedMesh<Dimension>;
+
+        std::vector<NamedMesh> mesh_list = [] {
+          std::vector<NamedMesh> extended_mesh_list;
+          std::array mesh_array = MeshDataBaseForTests::get().all3DMeshes();
+          for (size_t i = 0; i < mesh_array.size(); ++i) {
+            extended_mesh_list.push_back(MeshDataBaseForTests::get().all3DMeshes()[i]);
+          }
+          extended_mesh_list.push_back(NamedMesh("diamond dual", DualMeshManager::instance().getDiamondDualMesh(
+                                                                   *MeshDataBaseForTests::get().hybrid3DMesh())));
+          return extended_mesh_list;
+        }();
+
+        for (auto named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_3d = named_mesh.mesh();
+
+            std::string_view data = R"(
+import math;
+let scalar_3d: R^3 -> R, x -> 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
+let R3_3d: R^3 -> R^3, x -> (2 * exp(x[0]) * sin(x[1]) + x[2] + 3, x[0] * x[2] - 2 * x[1], 3);
+let R2x2_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3, x[0] * x[1] * x[2]);
+)";
+
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            SECTION("scalar 3d")
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<double> cell_integral{mesh_3d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) * x[2] + 3; };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_3d, cell_integral);
+
+              Array<double> integrate_value(mesh_3d->numberOfCells());
+              IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_3d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("vector 3d")
+            {
+              using R3               = TinyVector<3>;
+              auto [i_symbol, found] = symbol_table->find("R3_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<R3> cell_integral{mesh_3d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> R3 {
+                return R3{2 * exp(x[0]) * sin(x[1]) + x[2] + 3, x[0] * x[2] - 2 * x[1], 3};
+              };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_3d, cell_integral);
+
+              Array<R3> integrate_value(mesh_3d->numberOfCells());
+              IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_3d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("matrix 3d")
+            {
+              using R2x2             = TinyMatrix<2>;
+              auto [i_symbol, found] = symbol_table->find("R2x2_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<R2x2> cell_integral{mesh_3d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+                return R2x2{2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3,
+                            x[0] * x[1] * x[2]};
+              };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_3d, cell_integral);
+
+              Array<R2x2> integrate_value(mesh_3d->numberOfCells());
+              IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_3d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+          }
+        }
+      }
+    }
+
+    SECTION("integrate on cell list")
+    {
+      auto same_item_integral = [](auto f, auto g) -> bool {
+        using ItemIdType = typename decltype(f)::index_type;
+        for (ItemIdType item_id = 0; item_id < f.size(); ++item_id) {
+          if (f[item_id] != g[item_id]) {
+            return false;
+          }
+        }
+
+        return true;
+      };
+
+      SECTION("1D")
+      {
+        constexpr size_t Dimension = 1;
+
+        std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+        for (auto named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_1d = named_mesh.mesh();
+            Array<CellId> cell_list{mesh_1d->numberOfCells() / 2 + mesh_1d->numberOfCells() % 2};
+
+            {
+              size_t k = 0;
+              for (CellId cell_id = 0; cell_id < mesh_1d->numberOfCells(); ++(++cell_id), ++k) {
+                cell_list[k] = cell_id;
+              }
+
+              REQUIRE(k == cell_list.size());
+            }
+
+            std::string_view data = R"(
+import math;
+let scalar_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+let R3_1d: R^1 -> R^3, x -> (2 * exp(x[0]) + 3, x[0] - 2, 3);
+let R2x2_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0]);
+)";
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            SECTION("scalar 1d")
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * std::exp(x[0]) + 3; };
+
+              Array<const double> cell_integral =
+                CellIntegrator::integrate(f, quadrature_descriptor, *mesh_1d, cell_list);
+              Array<const double> integrate_value =
+                IntegrateOnCells<double(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_1d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("vector 1d")
+            {
+              using R3               = TinyVector<3>;
+              auto [i_symbol, found] = symbol_table->find("R3_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> R3 { return R3{2 * exp(x[0]) + 3, x[0] - 2, 3}; };
+
+              Array<const R3> cell_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_1d, cell_list);
+              Array<const R3> integrate_value =
+                IntegrateOnCells<R3(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_1d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("matrix 1d")
+            {
+              using R2x2             = TinyMatrix<2>;
+              auto [i_symbol, found] = symbol_table->find("R2x2_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+                return R2x2{2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0]};
+              };
+
+              Array<const R2x2> cell_integral =
+                CellIntegrator::integrate(f, quadrature_descriptor, *mesh_1d, cell_list);
+              Array<const R2x2> integrate_value =
+                IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_1d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+          }
+        }
+      }
+
+      SECTION("2D")
+      {
+        constexpr size_t Dimension = 2;
+
+        std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+        for (auto named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_2d = named_mesh.mesh();
+
+            Array<CellId> cell_list{mesh_2d->numberOfCells() / 2 + mesh_2d->numberOfCells() % 2};
+
+            {
+              size_t k = 0;
+              for (CellId cell_id = 0; cell_id < mesh_2d->numberOfCells(); ++(++cell_id), ++k) {
+                cell_list[k] = cell_id;
+              }
+
+              REQUIRE(k == cell_list.size());
+            }
+
+            std::string_view data = R"(
+import math;
+let scalar_2d: R^2 -> R, x -> 2*exp(x[0])*sin(x[1])+3;
+let R3_2d: R^2 -> R^3, x -> (2*exp(x[0])*sin(x[1])+3, x[0]-2*x[1], 3);
+let R2x2_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*x[1]), 3, x[0]*x[1]);
+)";
+
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            SECTION("scalar 2d")
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) + 3; };
+
+              Array<const double> cell_integral =
+                CellIntegrator::integrate(f, quadrature_descriptor, *mesh_2d, cell_list);
+              Array<const double> integrate_value =
+                IntegrateOnCells<double(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_2d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("vector 2d")
+            {
+              using R3               = TinyVector<3>;
+              auto [i_symbol, found] = symbol_table->find("R3_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> R3 {
+                return R3{2 * exp(x[0]) * sin(x[1]) + 3, x[0] - 2 * x[1], 3};
+              };
+
+              Array<const R3> cell_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_2d, cell_list);
+              Array<const R3> integrate_value =
+                IntegrateOnCells<R3(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_2d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("matrix 2d")
+            {
+              using R2x2             = TinyMatrix<2>;
+              auto [i_symbol, found] = symbol_table->find("R2x2_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+                return R2x2{2 * exp(x[0]) * sin(x[1]) + 3, sin(x[0] - 2 * x[1]), 3, x[0] * x[1]};
+              };
+
+              Array<const R2x2> cell_integral =
+                CellIntegrator::integrate(f, quadrature_descriptor, *mesh_2d, cell_list);
+              Array<const R2x2> integrate_value =
+                IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_2d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+          }
+        }
+      }
+
+      SECTION("3D")
+      {
+        constexpr size_t Dimension = 3;
+        using NamedMesh            = MeshDataBaseForTests::NamedMesh<Dimension>;
+
+        std::vector<NamedMesh> mesh_list = [] {
+          std::vector<NamedMesh> extended_mesh_list;
+          std::array mesh_array = MeshDataBaseForTests::get().all3DMeshes();
+          for (size_t i = 0; i < mesh_array.size(); ++i) {
+            extended_mesh_list.push_back(MeshDataBaseForTests::get().all3DMeshes()[i]);
+          }
+          extended_mesh_list.push_back(NamedMesh("diamond dual", DualMeshManager::instance().getDiamondDualMesh(
+                                                                   *MeshDataBaseForTests::get().hybrid3DMesh())));
+          return extended_mesh_list;
+        }();
+
+        for (auto named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_3d = named_mesh.mesh();
+
+            Array<CellId> cell_list{mesh_3d->numberOfCells() / 2 + mesh_3d->numberOfCells() % 2};
+
+            {
+              size_t k = 0;
+              for (CellId cell_id = 0; cell_id < mesh_3d->numberOfCells(); ++(++cell_id), ++k) {
+                cell_list[k] = cell_id;
+              }
+
+              REQUIRE(k == cell_list.size());
+            }
+
+            std::string_view data = R"(
+import math;
+let scalar_3d: R^3 -> R, x -> 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
+let R3_3d: R^3 -> R^3, x -> (2 * exp(x[0]) * sin(x[1]) + x[2] + 3, x[0] * x[2] - 2 * x[1], 3);
+let R2x2_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3, x[0] * x[1] * x[2]);
+)";
+
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            SECTION("scalar 3d")
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) * x[2] + 3; };
+
+              Array<const double> cell_integral =
+                CellIntegrator::integrate(f, quadrature_descriptor, *mesh_3d, cell_list);
+              Array<const double> integrate_value =
+                IntegrateOnCells<double(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_3d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("vector 3d")
+            {
+              using R3               = TinyVector<3>;
+              auto [i_symbol, found] = symbol_table->find("R3_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> R3 {
+                return R3{2 * exp(x[0]) * sin(x[1]) + x[2] + 3, x[0] * x[2] - 2 * x[1], 3};
+              };
+
+              Array<const R3> cell_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_3d, cell_list);
+              Array<const R3> integrate_value =
+                IntegrateOnCells<R3(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_3d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("matrix 3d")
+            {
+              using R2x2             = TinyMatrix<2>;
+              auto [i_symbol, found] = symbol_table->find("R2x2_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+                return R2x2{2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3,
+                            x[0] * x[1] * x[2]};
+              };
+
+              Array<const R2x2> cell_integral =
+                CellIntegrator::integrate(f, quadrature_descriptor, *mesh_3d, cell_list);
+              Array<const R2x2> integrate_value =
+                IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_3d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+          }
+        }
+      }
+    }
+  }
+
+  SECTION("Gauss-Legendre quadrature")
+  {
+    auto quadrature_descriptor = GaussLegendreQuadratureDescriptor(3);
+
+    SECTION("integrate on all cells")
+    {
+      auto same_item_integral = [](auto f, auto g) -> bool {
+        using ItemIdType = typename decltype(f)::index_type;
+        for (ItemIdType item_id = 0; item_id < f.numberOfItems(); ++item_id) {
+          if (f[item_id] != g[item_id]) {
+            return false;
+          }
+        }
+
+        return true;
+      };
+
+      SECTION("1D")
+      {
+        constexpr size_t Dimension = 1;
+
+        std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+        for (auto named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_1d = named_mesh.mesh();
+
+            std::string_view data = R"(
+import math;
+let scalar_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+let R3_1d: R^1 -> R^3, x -> (2 * exp(x[0]) + 3, x[0] - 2, 3);
+let R2x2_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0]);
+)";
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            SECTION("scalar 1d")
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<double> cell_integral{mesh_1d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * std::exp(x[0]) + 3; };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_1d, cell_integral);
+
+              Array<double> integrate_value(mesh_1d->numberOfCells());
+              IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_1d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("vector 1d")
+            {
+              using R3               = TinyVector<3>;
+              auto [i_symbol, found] = symbol_table->find("R3_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<R3> cell_integral{mesh_1d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> R3 { return R3{2 * exp(x[0]) + 3, x[0] - 2, 3}; };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_1d, cell_integral);
+
+              Array<R3> integrate_value(mesh_1d->numberOfCells());
+              IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_1d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("matrix 1d")
+            {
+              using R2x2             = TinyMatrix<2>;
+              auto [i_symbol, found] = symbol_table->find("R2x2_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<R2x2> cell_integral{mesh_1d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+                return R2x2{2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0]};
+              };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_1d, cell_integral);
+
+              Array<R2x2> integrate_value(mesh_1d->numberOfCells());
+              IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_1d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+          }
+        }
+      }
+
+      SECTION("2D")
+      {
+        constexpr size_t Dimension = 2;
+
+        std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+        for (auto named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_2d = named_mesh.mesh();
+
+            std::string_view data = R"(
+import math;
+let scalar_2d: R^2 -> R, x -> 2*exp(x[0])*sin(x[1])+3;
+let R3_2d: R^2 -> R^3, x -> (2*exp(x[0])*sin(x[1])+3, x[0]-2*x[1], 3);
+let R2x2_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*x[1]), 3, x[0]*x[1]);
+)";
+
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            SECTION("scalar 2d")
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<double> cell_integral{mesh_2d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) + 3; };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_2d, cell_integral);
+
+              Array<double> integrate_value(mesh_2d->numberOfCells());
+              IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_2d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("vector 2d")
+            {
+              using R3               = TinyVector<3>;
+              auto [i_symbol, found] = symbol_table->find("R3_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<R3> cell_integral{mesh_2d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> R3 {
+                return R3{2 * exp(x[0]) * sin(x[1]) + 3, x[0] - 2 * x[1], 3};
+              };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_2d, cell_integral);
+
+              Array<R3> integrate_value(mesh_2d->numberOfCells());
+              IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_2d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("matrix 2d")
+            {
+              using R2x2             = TinyMatrix<2>;
+              auto [i_symbol, found] = symbol_table->find("R2x2_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<R2x2> cell_integral{mesh_2d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+                return R2x2{2 * exp(x[0]) * sin(x[1]) + 3, sin(x[0] - 2 * x[1]), 3, x[0] * x[1]};
+              };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_2d, cell_integral);
+
+              Array<R2x2> integrate_value(mesh_2d->numberOfCells());
+              IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_2d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+          }
+        }
+      }
+
+      SECTION("3D")
+      {
+        constexpr size_t Dimension = 3;
+        using NamedMesh            = MeshDataBaseForTests::NamedMesh<Dimension>;
+
+        std::vector<NamedMesh> mesh_list = [] {
+          std::vector<NamedMesh> extended_mesh_list;
+          std::array mesh_array = MeshDataBaseForTests::get().all3DMeshes();
+          for (size_t i = 0; i < mesh_array.size(); ++i) {
+            extended_mesh_list.push_back(MeshDataBaseForTests::get().all3DMeshes()[i]);
+          }
+          extended_mesh_list.push_back(NamedMesh("diamond dual", DualMeshManager::instance().getDiamondDualMesh(
+                                                                   *MeshDataBaseForTests::get().hybrid3DMesh())));
+          return extended_mesh_list;
+        }();
+
+        for (auto named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_3d = named_mesh.mesh();
+
+            std::string_view data = R"(
+import math;
+let scalar_3d: R^3 -> R, x -> 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
+let R3_3d: R^3 -> R^3, x -> (2 * exp(x[0]) * sin(x[1]) + x[2] + 3, x[0] * x[2] - 2 * x[1], 3);
+let R2x2_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3, x[0] * x[1] * x[2]);
+)";
+
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            SECTION("scalar 3d")
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<double> cell_integral{mesh_3d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) * x[2] + 3; };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_3d, cell_integral);
+
+              Array<double> integrate_value(mesh_3d->numberOfCells());
+              IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_3d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("vector 3d")
+            {
+              using R3               = TinyVector<3>;
+              auto [i_symbol, found] = symbol_table->find("R3_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<R3> cell_integral{mesh_3d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> R3 {
+                return R3{2 * exp(x[0]) * sin(x[1]) + x[2] + 3, x[0] * x[2] - 2 * x[1], 3};
+              };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_3d, cell_integral);
+
+              Array<R3> integrate_value(mesh_3d->numberOfCells());
+              IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_3d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("matrix 3d")
+            {
+              using R2x2             = TinyMatrix<2>;
+              auto [i_symbol, found] = symbol_table->find("R2x2_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<R2x2> cell_integral{mesh_3d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+                return R2x2{2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3,
+                            x[0] * x[1] * x[2]};
+              };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_3d, cell_integral);
+
+              Array<R2x2> integrate_value(mesh_3d->numberOfCells());
+              IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_3d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+          }
+        }
+      }
+    }
+
+    SECTION("integrate on cell list")
+    {
+      auto same_item_integral = [](auto f, auto g) -> bool {
+        using ItemIdType = typename decltype(f)::index_type;
+        for (ItemIdType item_id = 0; item_id < f.size(); ++item_id) {
+          if (f[item_id] != g[item_id]) {
+            return false;
+          }
+        }
+
+        return true;
+      };
+
+      SECTION("1D")
+      {
+        constexpr size_t Dimension = 1;
+
+        std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+        for (auto named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_1d = named_mesh.mesh();
+            Array<CellId> cell_list{mesh_1d->numberOfCells() / 2 + mesh_1d->numberOfCells() % 2};
+
+            {
+              size_t k = 0;
+              for (CellId cell_id = 0; cell_id < mesh_1d->numberOfCells(); ++(++cell_id), ++k) {
+                cell_list[k] = cell_id;
+              }
+
+              REQUIRE(k == cell_list.size());
+            }
+
+            std::string_view data = R"(
+import math;
+let scalar_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+let R3_1d: R^1 -> R^3, x -> (2 * exp(x[0]) + 3, x[0] - 2, 3);
+let R2x2_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0]);
+)";
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            SECTION("scalar 1d")
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * std::exp(x[0]) + 3; };
+
+              Array<const double> cell_integral =
+                CellIntegrator::integrate(f, quadrature_descriptor, *mesh_1d, cell_list);
+              Array<const double> integrate_value =
+                IntegrateOnCells<double(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_1d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("vector 1d")
+            {
+              using R3               = TinyVector<3>;
+              auto [i_symbol, found] = symbol_table->find("R3_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> R3 { return R3{2 * exp(x[0]) + 3, x[0] - 2, 3}; };
+
+              Array<const R3> cell_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_1d, cell_list);
+              Array<const R3> integrate_value =
+                IntegrateOnCells<R3(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_1d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("matrix 1d")
+            {
+              using R2x2             = TinyMatrix<2>;
+              auto [i_symbol, found] = symbol_table->find("R2x2_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+                return R2x2{2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0]};
+              };
+
+              Array<const R2x2> cell_integral =
+                CellIntegrator::integrate(f, quadrature_descriptor, *mesh_1d, cell_list);
+              Array<const R2x2> integrate_value =
+                IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_1d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+          }
+        }
+      }
+
+      SECTION("2D")
+      {
+        constexpr size_t Dimension = 2;
+
+        std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+        for (auto named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_2d = named_mesh.mesh();
+
+            Array<CellId> cell_list{mesh_2d->numberOfCells() / 2 + mesh_2d->numberOfCells() % 2};
+
+            {
+              size_t k = 0;
+              for (CellId cell_id = 0; cell_id < mesh_2d->numberOfCells(); ++(++cell_id), ++k) {
+                cell_list[k] = cell_id;
+              }
+
+              REQUIRE(k == cell_list.size());
+            }
+
+            std::string_view data = R"(
+import math;
+let scalar_2d: R^2 -> R, x -> 2*exp(x[0])*sin(x[1])+3;
+let R3_2d: R^2 -> R^3, x -> (2*exp(x[0])*sin(x[1])+3, x[0]-2*x[1], 3);
+let R2x2_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*x[1]), 3, x[0]*x[1]);
+)";
+
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            SECTION("scalar 2d")
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) + 3; };
+
+              Array<const double> cell_integral =
+                CellIntegrator::integrate(f, quadrature_descriptor, *mesh_2d, cell_list);
+              Array<const double> integrate_value =
+                IntegrateOnCells<double(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_2d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("vector 2d")
+            {
+              using R3               = TinyVector<3>;
+              auto [i_symbol, found] = symbol_table->find("R3_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> R3 {
+                return R3{2 * exp(x[0]) * sin(x[1]) + 3, x[0] - 2 * x[1], 3};
+              };
+
+              Array<const R3> cell_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_2d, cell_list);
+              Array<const R3> integrate_value =
+                IntegrateOnCells<R3(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_2d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("matrix 2d")
+            {
+              using R2x2             = TinyMatrix<2>;
+              auto [i_symbol, found] = symbol_table->find("R2x2_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+                return R2x2{2 * exp(x[0]) * sin(x[1]) + 3, sin(x[0] - 2 * x[1]), 3, x[0] * x[1]};
+              };
+
+              Array<const R2x2> cell_integral =
+                CellIntegrator::integrate(f, quadrature_descriptor, *mesh_2d, cell_list);
+              Array<const R2x2> integrate_value =
+                IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_2d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+          }
+        }
+      }
+
+      SECTION("3D")
+      {
+        constexpr size_t Dimension = 3;
+        using NamedMesh            = MeshDataBaseForTests::NamedMesh<Dimension>;
+
+        std::vector<NamedMesh> mesh_list = [] {
+          std::vector<NamedMesh> extended_mesh_list;
+          std::array mesh_array = MeshDataBaseForTests::get().all3DMeshes();
+          for (size_t i = 0; i < mesh_array.size(); ++i) {
+            extended_mesh_list.push_back(MeshDataBaseForTests::get().all3DMeshes()[i]);
+          }
+          extended_mesh_list.push_back(NamedMesh("diamond dual", DualMeshManager::instance().getDiamondDualMesh(
+                                                                   *MeshDataBaseForTests::get().hybrid3DMesh())));
+          return extended_mesh_list;
+        }();
+
+        for (auto named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_3d = named_mesh.mesh();
+
+            Array<CellId> cell_list{mesh_3d->numberOfCells() / 2 + mesh_3d->numberOfCells() % 2};
+
+            {
+              size_t k = 0;
+              for (CellId cell_id = 0; cell_id < mesh_3d->numberOfCells(); ++(++cell_id), ++k) {
+                cell_list[k] = cell_id;
+              }
+
+              REQUIRE(k == cell_list.size());
+            }
+
+            std::string_view data = R"(
+import math;
+let scalar_3d: R^3 -> R, x -> 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
+let R3_3d: R^3 -> R^3, x -> (2 * exp(x[0]) * sin(x[1]) + x[2] + 3, x[0] * x[2] - 2 * x[1], 3);
+let R2x2_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3, x[0] * x[1] * x[2]);
+)";
+
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            SECTION("scalar 3d")
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) * x[2] + 3; };
+
+              Array<const double> cell_integral =
+                CellIntegrator::integrate(f, quadrature_descriptor, *mesh_3d, cell_list);
+              Array<const double> integrate_value =
+                IntegrateOnCells<double(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_3d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("vector 3d")
+            {
+              using R3               = TinyVector<3>;
+              auto [i_symbol, found] = symbol_table->find("R3_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> R3 {
+                return R3{2 * exp(x[0]) * sin(x[1]) + x[2] + 3, x[0] * x[2] - 2 * x[1], 3};
+              };
+
+              Array<const R3> cell_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_3d, cell_list);
+              Array<const R3> integrate_value =
+                IntegrateOnCells<R3(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_3d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("matrix 3d")
+            {
+              using R2x2             = TinyMatrix<2>;
+              auto [i_symbol, found] = symbol_table->find("R2x2_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+                return R2x2{2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3,
+                            x[0] * x[1] * x[2]};
+              };
+
+              Array<const R2x2> cell_integral =
+                CellIntegrator::integrate(f, quadrature_descriptor, *mesh_3d, cell_list);
+              Array<const R2x2> integrate_value =
+                IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_3d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+          }
+        }
+      }
+    }
+  }
+
+  SECTION("Gauss-Lobatto quadrature")
+  {
+    auto quadrature_descriptor = GaussLobattoQuadratureDescriptor(3);
+
+    SECTION("integrate on all cells")
+    {
+      auto same_item_integral = [](auto f, auto g) -> bool {
+        using ItemIdType = typename decltype(f)::index_type;
+        for (ItemIdType item_id = 0; item_id < f.numberOfItems(); ++item_id) {
+          if (f[item_id] != g[item_id]) {
+            return false;
+          }
+        }
+
+        return true;
+      };
+
+      SECTION("1D")
+      {
+        constexpr size_t Dimension = 1;
+
+        std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+        for (auto named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_1d = named_mesh.mesh();
+
+            std::string_view data = R"(
+import math;
+let scalar_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+let R3_1d: R^1 -> R^3, x -> (2 * exp(x[0]) + 3, x[0] - 2, 3);
+let R2x2_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0]);
+)";
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            SECTION("scalar 1d")
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<double> cell_integral{mesh_1d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * std::exp(x[0]) + 3; };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_1d, cell_integral);
+
+              Array<double> integrate_value(mesh_1d->numberOfCells());
+              IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_1d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("vector 1d")
+            {
+              using R3               = TinyVector<3>;
+              auto [i_symbol, found] = symbol_table->find("R3_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<R3> cell_integral{mesh_1d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> R3 { return R3{2 * exp(x[0]) + 3, x[0] - 2, 3}; };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_1d, cell_integral);
+
+              Array<R3> integrate_value(mesh_1d->numberOfCells());
+              IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_1d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("matrix 1d")
+            {
+              using R2x2             = TinyMatrix<2>;
+              auto [i_symbol, found] = symbol_table->find("R2x2_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<R2x2> cell_integral{mesh_1d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+                return R2x2{2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0]};
+              };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_1d, cell_integral);
+
+              Array<R2x2> integrate_value(mesh_1d->numberOfCells());
+              IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_1d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+          }
+        }
+      }
+
+      SECTION("2D")
+      {
+        constexpr size_t Dimension = 2;
+
+        std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+        for (auto named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_2d = named_mesh.mesh();
+
+            std::string_view data = R"(
+import math;
+let scalar_2d: R^2 -> R, x -> 2*exp(x[0])*sin(x[1])+3;
+let R3_2d: R^2 -> R^3, x -> (2*exp(x[0])*sin(x[1])+3, x[0]-2*x[1], 3);
+let R2x2_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*x[1]), 3, x[0]*x[1]);
+)";
+
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            SECTION("scalar 2d")
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<double> cell_integral{mesh_2d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) + 3; };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_2d, cell_integral);
+
+              Array<double> integrate_value(mesh_2d->numberOfCells());
+              IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_2d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("vector 2d")
+            {
+              using R3               = TinyVector<3>;
+              auto [i_symbol, found] = symbol_table->find("R3_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<R3> cell_integral{mesh_2d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> R3 {
+                return R3{2 * exp(x[0]) * sin(x[1]) + 3, x[0] - 2 * x[1], 3};
+              };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_2d, cell_integral);
+
+              Array<R3> integrate_value(mesh_2d->numberOfCells());
+              IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_2d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("matrix 2d")
+            {
+              using R2x2             = TinyMatrix<2>;
+              auto [i_symbol, found] = symbol_table->find("R2x2_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<R2x2> cell_integral{mesh_2d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+                return R2x2{2 * exp(x[0]) * sin(x[1]) + 3, sin(x[0] - 2 * x[1]), 3, x[0] * x[1]};
+              };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_2d, cell_integral);
+
+              Array<R2x2> integrate_value(mesh_2d->numberOfCells());
+              IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_2d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+          }
+        }
+      }
+
+      SECTION("3D")
+      {
+        constexpr size_t Dimension = 3;
+        using NamedMesh            = MeshDataBaseForTests::NamedMesh<Dimension>;
+
+        std::vector<NamedMesh> mesh_list = [] {
+          std::vector<NamedMesh> extended_mesh_list;
+          std::array mesh_array = MeshDataBaseForTests::get().all3DMeshes();
+          for (size_t i = 0; i < mesh_array.size(); ++i) {
+            extended_mesh_list.push_back(MeshDataBaseForTests::get().all3DMeshes()[i]);
+          }
+          extended_mesh_list.push_back(NamedMesh("diamond dual", DualMeshManager::instance().getDiamondDualMesh(
+                                                                   *MeshDataBaseForTests::get().hybrid3DMesh())));
+          return extended_mesh_list;
+        }();
+
+        for (auto named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_3d = named_mesh.mesh();
+
+            std::string_view data = R"(
+import math;
+let scalar_3d: R^3 -> R, x -> 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
+let R3_3d: R^3 -> R^3, x -> (2 * exp(x[0]) * sin(x[1]) + x[2] + 3, x[0] * x[2] - 2 * x[1], 3);
+let R2x2_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3, x[0] * x[1] * x[2]);
+)";
+
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            SECTION("scalar 3d")
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<double> cell_integral{mesh_3d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) * x[2] + 3; };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_3d, cell_integral);
+
+              Array<double> integrate_value(mesh_3d->numberOfCells());
+              IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_3d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("vector 3d")
+            {
+              using R3               = TinyVector<3>;
+              auto [i_symbol, found] = symbol_table->find("R3_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<R3> cell_integral{mesh_3d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> R3 {
+                return R3{2 * exp(x[0]) * sin(x[1]) + x[2] + 3, x[0] * x[2] - 2 * x[1], 3};
+              };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_3d, cell_integral);
+
+              Array<R3> integrate_value(mesh_3d->numberOfCells());
+              IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_3d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("matrix 3d")
+            {
+              using R2x2             = TinyMatrix<2>;
+              auto [i_symbol, found] = symbol_table->find("R2x2_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              CellValue<R2x2> cell_integral{mesh_3d->connectivity()};
+              auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+                return R2x2{2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3,
+                            x[0] * x[1] * x[2]};
+              };
+              CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_3d, cell_integral);
+
+              Array<R2x2> integrate_value(mesh_3d->numberOfCells());
+              IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_3d, integrate_value);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+          }
+        }
+      }
+    }
+
+    SECTION("integrate on cell list")
+    {
+      auto same_item_integral = [](auto f, auto g) -> bool {
+        using ItemIdType = typename decltype(f)::index_type;
+        for (ItemIdType item_id = 0; item_id < f.size(); ++item_id) {
+          if (f[item_id] != g[item_id]) {
+            return false;
+          }
+        }
+
+        return true;
+      };
+
+      SECTION("1D")
+      {
+        constexpr size_t Dimension = 1;
+
+        std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+        for (auto named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_1d = named_mesh.mesh();
+            Array<CellId> cell_list{mesh_1d->numberOfCells() / 2 + mesh_1d->numberOfCells() % 2};
+
+            {
+              size_t k = 0;
+              for (CellId cell_id = 0; cell_id < mesh_1d->numberOfCells(); ++(++cell_id), ++k) {
+                cell_list[k] = cell_id;
+              }
+
+              REQUIRE(k == cell_list.size());
+            }
+
+            std::string_view data = R"(
+import math;
+let scalar_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+let R3_1d: R^1 -> R^3, x -> (2 * exp(x[0]) + 3, x[0] - 2, 3);
+let R2x2_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0]);
+)";
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            SECTION("scalar 1d")
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * std::exp(x[0]) + 3; };
+
+              Array<const double> cell_integral =
+                CellIntegrator::integrate(f, quadrature_descriptor, *mesh_1d, cell_list);
+              Array<const double> integrate_value =
+                IntegrateOnCells<double(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_1d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("vector 1d")
+            {
+              using R3               = TinyVector<3>;
+              auto [i_symbol, found] = symbol_table->find("R3_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> R3 { return R3{2 * exp(x[0]) + 3, x[0] - 2, 3}; };
+
+              Array<const R3> cell_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_1d, cell_list);
+              Array<const R3> integrate_value =
+                IntegrateOnCells<R3(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_1d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("matrix 1d")
+            {
+              using R2x2             = TinyMatrix<2>;
+              auto [i_symbol, found] = symbol_table->find("R2x2_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+                return R2x2{2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0]};
+              };
+
+              Array<const R2x2> cell_integral =
+                CellIntegrator::integrate(f, quadrature_descriptor, *mesh_1d, cell_list);
+              Array<const R2x2> integrate_value =
+                IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_1d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+          }
+        }
+      }
+
+      SECTION("2D")
+      {
+        constexpr size_t Dimension = 2;
+
+        std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+        for (auto named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_2d = named_mesh.mesh();
+
+            Array<CellId> cell_list{mesh_2d->numberOfCells() / 2 + mesh_2d->numberOfCells() % 2};
+
+            {
+              size_t k = 0;
+              for (CellId cell_id = 0; cell_id < mesh_2d->numberOfCells(); ++(++cell_id), ++k) {
+                cell_list[k] = cell_id;
+              }
+
+              REQUIRE(k == cell_list.size());
+            }
+
+            std::string_view data = R"(
+import math;
+let scalar_2d: R^2 -> R, x -> 2*exp(x[0])*sin(x[1])+3;
+let R3_2d: R^2 -> R^3, x -> (2*exp(x[0])*sin(x[1])+3, x[0]-2*x[1], 3);
+let R2x2_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*x[1]), 3, x[0]*x[1]);
+)";
+
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            SECTION("scalar 2d")
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) + 3; };
+
+              Array<const double> cell_integral =
+                CellIntegrator::integrate(f, quadrature_descriptor, *mesh_2d, cell_list);
+              Array<const double> integrate_value =
+                IntegrateOnCells<double(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_2d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("vector 2d")
+            {
+              using R3               = TinyVector<3>;
+              auto [i_symbol, found] = symbol_table->find("R3_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> R3 {
+                return R3{2 * exp(x[0]) * sin(x[1]) + 3, x[0] - 2 * x[1], 3};
+              };
+
+              Array<const R3> cell_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_2d, cell_list);
+              Array<const R3> integrate_value =
+                IntegrateOnCells<R3(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_2d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("matrix 2d")
+            {
+              using R2x2             = TinyMatrix<2>;
+              auto [i_symbol, found] = symbol_table->find("R2x2_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+                return R2x2{2 * exp(x[0]) * sin(x[1]) + 3, sin(x[0] - 2 * x[1]), 3, x[0] * x[1]};
+              };
+
+              Array<const R2x2> cell_integral =
+                CellIntegrator::integrate(f, quadrature_descriptor, *mesh_2d, cell_list);
+              Array<const R2x2> integrate_value =
+                IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_2d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+          }
+        }
+      }
+
+      SECTION("3D")
+      {
+        constexpr size_t Dimension = 3;
+        using NamedMesh            = MeshDataBaseForTests::NamedMesh<Dimension>;
+
+        std::vector<NamedMesh> mesh_list = [] {
+          std::vector<NamedMesh> extended_mesh_list;
+          std::array mesh_array = MeshDataBaseForTests::get().all3DMeshes();
+          for (size_t i = 0; i < mesh_array.size(); ++i) {
+            extended_mesh_list.push_back(MeshDataBaseForTests::get().all3DMeshes()[i]);
+          }
+          extended_mesh_list.push_back(NamedMesh("diamond dual", DualMeshManager::instance().getDiamondDualMesh(
+                                                                   *MeshDataBaseForTests::get().hybrid3DMesh())));
+          return extended_mesh_list;
+        }();
+
+        for (auto named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_3d = named_mesh.mesh();
+
+            Array<CellId> cell_list{mesh_3d->numberOfCells() / 2 + mesh_3d->numberOfCells() % 2};
+
+            {
+              size_t k = 0;
+              for (CellId cell_id = 0; cell_id < mesh_3d->numberOfCells(); ++(++cell_id), ++k) {
+                cell_list[k] = cell_id;
+              }
+
+              REQUIRE(k == cell_list.size());
+            }
+
+            std::string_view data = R"(
+import math;
+let scalar_3d: R^3 -> R, x -> 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
+let R3_3d: R^3 -> R^3, x -> (2 * exp(x[0]) * sin(x[1]) + x[2] + 3, x[0] * x[2] - 2 * x[1], 3);
+let R2x2_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3, x[0] * x[1] * x[2]);
+)";
+
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            SECTION("scalar 3d")
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) * x[2] + 3; };
+
+              Array<const double> cell_integral =
+                CellIntegrator::integrate(f, quadrature_descriptor, *mesh_3d, cell_list);
+              Array<const double> integrate_value =
+                IntegrateOnCells<double(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_3d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("vector 3d")
+            {
+              using R3               = TinyVector<3>;
+              auto [i_symbol, found] = symbol_table->find("R3_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> R3 {
+                return R3{2 * exp(x[0]) * sin(x[1]) + x[2] + 3, x[0] * x[2] - 2 * x[1], 3};
+              };
+
+              Array<const R3> cell_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_3d, cell_list);
+              Array<const R3> integrate_value =
+                IntegrateOnCells<R3(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_3d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+
+            SECTION("matrix 3d")
+            {
+              using R2x2             = TinyMatrix<2>;
+              auto [i_symbol, found] = symbol_table->find("R2x2_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+              auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+                return R2x2{2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3,
+                            x[0] * x[1] * x[2]};
+              };
+
+              Array<const R2x2> cell_integral =
+                CellIntegrator::integrate(f, quadrature_descriptor, *mesh_3d, cell_list);
+              Array<const R2x2> integrate_value =
+                IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_3d, cell_list);
+
+              REQUIRE(same_item_integral(cell_integral, integrate_value));
+            }
+          }
+        }
+      }
+    }
+  }
+}
diff --git a/tests/test_InterpolateItemArray.cpp b/tests/test_InterpolateItemArray.cpp
index 2974e01d8732d35702b32493f0cb810932de2df2..7fd968be9aaeed6cd7774353df69b2153eb68d51 100644
--- a/tests/test_InterpolateItemArray.cpp
+++ b/tests/test_InterpolateItemArray.cpp
@@ -297,10 +297,10 @@ let scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
           }();
 
           std::string_view data = R"(
-  import math;
-  let scalar_affine_1d: R^1 -> R, x -> 2*x[0] + 2;
-  let scalar_non_linear_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
-  )";
+import math;
+let scalar_affine_1d: R^1 -> R, x -> 2*x[0] + 2;
+let scalar_non_linear_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+)";
           TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
 
           auto ast = ASTBuilder::build(input);
diff --git a/tests/test_InterpolateItemValue.cpp b/tests/test_InterpolateItemValue.cpp
index fabbda2560ac969e6305586041f531d0b348a0f8..1d7032c660ee001d652339abc0d2daee8a825fff 100644
--- a/tests/test_InterpolateItemValue.cpp
+++ b/tests/test_InterpolateItemValue.cpp
@@ -28,7 +28,7 @@ TEST_CASE("InterpolateItemValue", "[language]")
 {
   SECTION("interpolate on all items")
   {
-    auto same_cell_value = [](auto f, auto g) -> bool {
+    auto same_item_value = [](auto f, auto g) -> bool {
       using ItemIdType = typename decltype(f)::index_type;
       for (ItemIdType item_id = 0; item_id < f.numberOfItems(); ++item_id) {
         if (f[item_id] != g[item_id]) {
@@ -98,7 +98,7 @@ let R2x2_non_linear_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x
             CellValue<const double> interpolate_value =
               InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("scalar_non_linear_1d")
@@ -119,7 +119,7 @@ let R2x2_non_linear_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x
             CellValue<const double> interpolate_value =
               InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R3_affine_1d")
@@ -140,7 +140,7 @@ let R2x2_non_linear_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x
             CellValue<const TinyVector<3>> interpolate_value =
               InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R3_non_linear_1d")
@@ -161,7 +161,7 @@ let R2x2_non_linear_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x
             CellValue<const TinyVector<3>> interpolate_value =
               InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R2x2_affine_1d")
@@ -182,7 +182,7 @@ let R2x2_non_linear_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x
             CellValue<const TinyMatrix<2>> interpolate_value =
               InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R2x2_non_linear_1d")
@@ -204,7 +204,7 @@ let R2x2_non_linear_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x
             CellValue<const TinyMatrix<2>> interpolate_value =
               InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
         }
       }
@@ -268,7 +268,7 @@ let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*
             CellValue<const double> interpolate_value =
               InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("scalar_non_linear_2d")
@@ -288,7 +288,7 @@ let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*
             CellValue<const double> interpolate_value =
               InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R3_affine_2d")
@@ -308,7 +308,7 @@ let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*
             CellValue<const TinyVector<3>> interpolate_value =
               InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R3_non_linear_2d")
@@ -328,7 +328,7 @@ let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*
             CellValue<const TinyVector<3>> interpolate_value =
               InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R2x2_affine_2d")
@@ -348,7 +348,7 @@ let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*
             CellValue<const TinyMatrix<2>> interpolate_value =
               InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R2x2_non_linear_2d")
@@ -369,7 +369,7 @@ let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*
             CellValue<const TinyMatrix<2>> interpolate_value =
               InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
         }
       }
@@ -433,7 +433,7 @@ let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(
             CellValue<const double> interpolate_value =
               InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("scalar_non_linear_3d")
@@ -453,7 +453,7 @@ let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(
             CellValue<const double> interpolate_value =
               InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R3_affine_3d")
@@ -473,7 +473,7 @@ let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(
             CellValue<const TinyVector<3>> interpolate_value =
               InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R3_non_linear_3d")
@@ -493,7 +493,7 @@ let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(
             CellValue<const TinyVector<3>> interpolate_value =
               InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R2x2_affine_3d")
@@ -514,7 +514,7 @@ let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(
             CellValue<const TinyMatrix<2>> interpolate_value =
               InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R2x2_non_linear_3d")
@@ -535,16 +535,16 @@ let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(
             CellValue<const TinyMatrix<2>> interpolate_value =
               InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
         }
       }
     }
   }
 
-  SECTION("interpolate on items list")
+  SECTION("interpolate on item list")
   {
-    auto same_cell_value = [](auto interpolated, auto reference) -> bool {
+    auto same_item_value = [](auto interpolated, auto reference) -> bool {
       for (size_t i = 0; i < interpolated.size(); ++i) {
         if (interpolated[i] != reference[i]) {
           return false;
@@ -621,7 +621,7 @@ let R2x2_non_linear_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x
             Array<const double> interpolate_value =
               InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("scalar_non_linear_1d")
@@ -642,7 +642,7 @@ let R2x2_non_linear_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x
             Array<const double> interpolate_value =
               InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R3_affine_1d")
@@ -664,7 +664,7 @@ let R2x2_non_linear_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x
               InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj,
                                                                                       cell_id_list);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R3_non_linear_1d")
@@ -686,7 +686,7 @@ let R2x2_non_linear_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x
               InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj,
                                                                                       cell_id_list);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R2x2_affine_1d")
@@ -708,7 +708,7 @@ let R2x2_non_linear_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x
               InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj,
                                                                                       cell_id_list);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R2x2_non_linear_1d")
@@ -730,7 +730,7 @@ let R2x2_non_linear_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x
               InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj,
                                                                                       cell_id_list);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
         }
       }
@@ -800,7 +800,7 @@ let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*
             Array<const double> interpolate_value =
               InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("scalar_non_linear_2d")
@@ -820,7 +820,7 @@ let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*
             Array<const double> interpolate_value =
               InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R3_affine_2d")
@@ -841,7 +841,7 @@ let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*
               InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj,
                                                                                       cell_id_list);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R3_non_linear_2d")
@@ -862,7 +862,7 @@ let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*
               InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj,
                                                                                       cell_id_list);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R2x2_affine_2d")
@@ -883,7 +883,7 @@ let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*
               InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj,
                                                                                       cell_id_list);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R2x2_non_linear_2d")
@@ -904,7 +904,7 @@ let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*
               InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj,
                                                                                       cell_id_list);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
         }
       }
@@ -973,7 +973,7 @@ let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(
             Array<const double> interpolate_value =
               InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("scalar_non_linear_3d")
@@ -993,7 +993,7 @@ let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(
             Array<const double> interpolate_value =
               InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R3_affine_3d")
@@ -1014,7 +1014,7 @@ let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(
               InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj,
                                                                                       cell_id_list);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R3_non_linear_3d")
@@ -1035,7 +1035,7 @@ let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(
               InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj,
                                                                                       cell_id_list);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R2x2_affine_3d")
@@ -1057,7 +1057,7 @@ let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(
               InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj,
                                                                                       cell_id_list);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
 
           SECTION("R2x2_non_linear_3d")
@@ -1079,7 +1079,7 @@ let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(
               InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj,
                                                                                       cell_id_list);
 
-            REQUIRE(same_cell_value(cell_value, interpolate_value));
+            REQUIRE(same_item_value(cell_value, interpolate_value));
           }
         }
       }
diff --git a/tests/test_ItemArray.cpp b/tests/test_ItemArray.cpp
index 3422d1da9c8e9627abf75fae743a826809ca5505..9b3b95fded76400eb5953d4a3d70f50481dadde8 100644
--- a/tests/test_ItemArray.cpp
+++ b/tests/test_ItemArray.cpp
@@ -232,12 +232,41 @@ TEST_CASE("ItemArray", "[mesh]")
         const Connectivity<2>& connectivity = mesh_2d->connectivity();
 
         WeakFaceArray<int> weak_face_array{connectivity, 5};
-
-        weak_face_array.fill(parallel::rank());
+        for (FaceId face_id = 0; face_id < mesh_2d->numberOfFaces(); ++face_id) {
+          for (size_t i = 0; i < weak_face_array.sizeOfArrays(); ++i) {
+            weak_face_array[face_id][i] = 2 * face_id + 3 * i;
+          }
+        }
 
         FaceArray<const int> face_array{weak_face_array};
 
         REQUIRE(face_array.connectivity_ptr() == weak_face_array.connectivity_ptr());
+
+        FaceArray<int> copied_face_array = copy(weak_face_array);
+        REQUIRE(copied_face_array.connectivity_ptr() == weak_face_array.connectivity_ptr());
+        REQUIRE(weak_face_array.sizeOfArrays() == copied_face_array.sizeOfArrays());
+
+        weak_face_array.fill(0);
+
+        {
+          bool is_same = true;
+          for (FaceId face_id = 0; face_id < mesh_2d->numberOfFaces(); ++face_id) {
+            for (size_t i = 0; i < face_array.sizeOfArrays(); ++i) {
+              is_same &= (face_array[face_id][i] == 0);
+            }
+          }
+          REQUIRE(is_same);
+        }
+
+        {
+          bool is_same = true;
+          for (FaceId face_id = 0; face_id < mesh_2d->numberOfFaces(); ++face_id) {
+            for (size_t i = 0; i < copied_face_array.sizeOfArrays(); ++i) {
+              is_same &= (copied_face_array[face_id][i] == static_cast<int>(2 * face_id + 3 * i));
+            }
+          }
+          REQUIRE(is_same);
+        }
       }
     }
   }
diff --git a/tests/test_ItemValue.cpp b/tests/test_ItemValue.cpp
index 84778227aa28e7a23267f9ec355822d257718c89..5ec3fdb892eb5474b80381a250d1f2dc487574de 100644
--- a/tests/test_ItemValue.cpp
+++ b/tests/test_ItemValue.cpp
@@ -130,8 +130,12 @@ TEST_CASE("ItemValue", "[mesh]")
 
         cell_value = values;
 
-        for (CellId i_cell = 0; i_cell < mesh_3d->numberOfCells(); ++i_cell) {
-          REQUIRE(cell_value[i_cell] == i_cell);
+        {
+          bool is_same = true;
+          for (CellId i_cell = 0; i_cell < mesh_3d->numberOfCells(); ++i_cell) {
+            is_same &= (cell_value[i_cell] == i_cell);
+          }
+          REQUIRE(is_same);
         }
       }
     }
@@ -156,12 +160,20 @@ TEST_CASE("ItemValue", "[mesh]")
 
         cell_value.fill(0);
 
-        for (CellId i_cell = 0; i_cell < mesh_3d->numberOfCells(); ++i_cell) {
-          REQUIRE(cell_value[i_cell] == 0);
+        {
+          bool is_same = true;
+          for (CellId i_cell = 0; i_cell < mesh_3d->numberOfCells(); ++i_cell) {
+            is_same &= (cell_value[i_cell] == 0);
+          }
+          REQUIRE(is_same);
         }
 
-        for (CellId i_cell = 0; i_cell < mesh_3d->numberOfCells(); ++i_cell) {
-          REQUIRE(const_cell_value[i_cell] == static_cast<std::int64_t>(parallel::rank()));
+        {
+          bool is_same = true;
+          for (CellId i_cell = 0; i_cell < mesh_3d->numberOfCells(); ++i_cell) {
+            is_same &= (const_cell_value[i_cell] == static_cast<std::int64_t>(parallel::rank()));
+          }
+          REQUIRE(is_same);
         }
       }
     }
@@ -179,12 +191,29 @@ TEST_CASE("ItemValue", "[mesh]")
         const Connectivity<2>& connectivity = mesh_2d->connectivity();
 
         WeakFaceValue<int> weak_face_value{connectivity};
-
-        weak_face_value.fill(parallel::rank());
+        for (FaceId face_id = 0; face_id < mesh_2d->numberOfFaces(); ++face_id) {
+          weak_face_value[face_id] = 2 * face_id + 1;
+        }
 
         FaceValue<const int> face_value{weak_face_value};
 
         REQUIRE(face_value.connectivity_ptr() == weak_face_value.connectivity_ptr());
+
+        FaceValue<int> copied_face_value = copy(weak_face_value);
+        REQUIRE(copied_face_value.connectivity_ptr() == weak_face_value.connectivity_ptr());
+
+        weak_face_value.fill(0);
+        for (FaceId face_id = 0; face_id < mesh_2d->numberOfFaces(); ++face_id) {
+          REQUIRE(face_value[face_id] == static_cast<int>(0));
+        }
+
+        {
+          bool is_same = true;
+          for (FaceId face_id = 0; face_id < mesh_2d->numberOfFaces(); ++face_id) {
+            is_same &= (copied_face_value[face_id] == static_cast<int>(2 * face_id + 1));
+          }
+          REQUIRE(is_same);
+        }
       }
     }
   }
diff --git a/tests/test_SubItemArrayPerItem.cpp b/tests/test_SubItemArrayPerItem.cpp
index c72f98b76de01f6562aac21c5cc008716606a0a5..71a226f01379cccd86ca48c27683ca0b467f80fa 100644
--- a/tests/test_SubItemArrayPerItem.cpp
+++ b/tests/test_SubItemArrayPerItem.cpp
@@ -835,7 +835,9 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
             }
           }
 
-          NodeArrayPerCell<size_t> copy_node_array_per_cell = copy(node_array_per_cell);
+          WeakNodeArrayPerCell<const size_t> node_const_array_per_cell = node_array_per_cell;
+
+          NodeArrayPerCell<size_t> copy_node_array_per_cell = copy(node_const_array_per_cell);
           {
             bool is_same = true;
             for (size_t i = 0; i < copy_node_array_per_cell.numberOfArrays(); ++i) {
@@ -874,6 +876,117 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
     }
   }
 
+  SECTION("fill")
+  {
+    std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+    for (auto named_mesh : mesh_list) {
+      SECTION(named_mesh.name())
+      {
+        auto mesh_3d = named_mesh.mesh();
+
+        const Connectivity<3>& connectivity = mesh_3d->connectivity();
+
+        SECTION("classic")
+        {
+          NodeArrayPerCell<size_t> node_array_per_cell{connectivity, 3};
+          {
+            size_t k = 0;
+            for (size_t i = 0; i < node_array_per_cell.numberOfArrays(); ++i) {
+              for (size_t j = 0; j < node_array_per_cell.sizeOfArrays(); ++j, ++k) {
+                node_array_per_cell[i][j] = k;
+              }
+            }
+          }
+          NodeArrayPerCell<size_t> copy_node_array_per_cell = copy(node_array_per_cell);
+          {
+            bool is_same = true;
+            for (size_t i = 0; i < copy_node_array_per_cell.numberOfArrays(); ++i) {
+              for (size_t j = 0; j < node_array_per_cell.sizeOfArrays(); ++j) {
+                is_same &= (copy_node_array_per_cell[i][j] == node_array_per_cell[i][j]);
+              }
+            }
+
+            REQUIRE(is_same);
+          }
+
+          node_array_per_cell.fill(11);
+
+          {
+            bool is_same = true;
+            for (size_t i = 0; i < copy_node_array_per_cell.numberOfArrays(); ++i) {
+              for (size_t j = 0; j < copy_node_array_per_cell.sizeOfArrays(); ++j) {
+                is_same &= (copy_node_array_per_cell[i][j] == node_array_per_cell[i][j]);
+              }
+            }
+
+            REQUIRE(not is_same);
+          }
+
+          {
+            bool is_filled_with_11 = true;
+            for (size_t i = 0; i < copy_node_array_per_cell.numberOfArrays(); ++i) {
+              for (size_t j = 0; j < copy_node_array_per_cell.sizeOfArrays(); ++j) {
+                is_filled_with_11 &= (node_array_per_cell[i][j] == 11);
+              }
+            }
+
+            REQUIRE(is_filled_with_11);
+          }
+        }
+
+        SECTION("from weak")
+        {
+          WeakNodeArrayPerCell<size_t> node_array_per_cell{connectivity, 3};
+          {
+            size_t k = 0;
+            for (size_t i = 0; i < node_array_per_cell.numberOfArrays(); ++i) {
+              for (size_t j = 0; j < node_array_per_cell.sizeOfArrays(); ++j, ++k) {
+                node_array_per_cell[i][j] = k;
+              }
+            }
+          }
+
+          NodeArrayPerCell<size_t> copy_node_array_per_cell = copy(node_array_per_cell);
+          {
+            bool is_same = true;
+            for (size_t i = 0; i < copy_node_array_per_cell.numberOfArrays(); ++i) {
+              for (size_t j = 0; j < node_array_per_cell.sizeOfArrays(); ++j) {
+                is_same &= (copy_node_array_per_cell[i][j] == node_array_per_cell[i][j]);
+              }
+            }
+
+            REQUIRE(is_same);
+          }
+
+          node_array_per_cell.fill(13);
+
+          {
+            bool is_same = true;
+            for (size_t i = 0; i < copy_node_array_per_cell.numberOfArrays(); ++i) {
+              for (size_t j = 0; j < node_array_per_cell.sizeOfArrays(); ++j) {
+                is_same &= (copy_node_array_per_cell[i][j] == node_array_per_cell[i][j]);
+              }
+            }
+
+            REQUIRE(not is_same);
+          }
+
+          {
+            bool is_filled_with_13 = true;
+            for (size_t i = 0; i < copy_node_array_per_cell.numberOfArrays(); ++i) {
+              for (size_t j = 0; j < copy_node_array_per_cell.sizeOfArrays(); ++j) {
+                is_filled_with_13 &= (node_array_per_cell[i][j] == 13);
+              }
+            }
+
+            REQUIRE(is_filled_with_13);
+          }
+        }
+      }
+    }
+  }
+
   SECTION("WeakSubItemArrayPerItem")
   {
     std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
diff --git a/tests/test_SubItemValuePerItem.cpp b/tests/test_SubItemValuePerItem.cpp
index 5eaadd94d88d07d2c802e9da867ab5ae2442aa5a..4063133a30f71e2304d865fb1250e975470c7ed0 100644
--- a/tests/test_SubItemValuePerItem.cpp
+++ b/tests/test_SubItemValuePerItem.cpp
@@ -720,7 +720,8 @@ TEST_CASE("SubItemValuePerItem", "[mesh]")
             }
           }
 
-          NodeValuePerCell<size_t> copy_node_value_per_cell = copy(node_value_per_cell);
+          WeakNodeValuePerCell<const size_t> node_const_value_per_cell = node_value_per_cell;
+          NodeValuePerCell<size_t> copy_node_value_per_cell            = copy(node_const_value_per_cell);
 
           {
             bool is_same = true;
diff --git a/tests/test_Synchronizer.cpp b/tests/test_Synchronizer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6a4a9b52a6ce05a4a0b2758c868a49dc758d3e3d
--- /dev/null
+++ b/tests/test_Synchronizer.cpp
@@ -0,0 +1,450 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/Synchronizer.hpp>
+#include <utils/pugs_config.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("Synchronizer", "[mesh]")
+{
+  auto is_same_item_value = [](auto a, auto b) {
+    using IndexT = typename decltype(a)::index_type;
+    bool is_same = true;
+    for (IndexT i = 0; i < a.numberOfItems(); ++i) {
+      is_same &= (a[i] == b[i]);
+    }
+    return parallel::allReduceAnd(is_same);
+  };
+
+  auto is_same_item_array = [](auto a, auto b) {
+    using IndexT = typename decltype(a)::index_type;
+    bool is_same = true;
+    for (IndexT i = 0; i < a.numberOfItems(); ++i) {
+      for (size_t j = 0; j < a.sizeOfArrays(); ++j) {
+        is_same &= (a[i][j] == b[i][j]);
+      }
+    }
+    return parallel::allReduceAnd(is_same);
+  };
+
+  SECTION("1D")
+  {
+    constexpr size_t Dimension = 1;
+    using ConnectivityType     = Connectivity<Dimension>;
+
+    const ConnectivityType& connectivity = MeshDataBaseForTests::get().unordered1DMesh()->connectivity();
+
+    SECTION("synchonize NodeValue")
+    {
+      const auto node_owner  = connectivity.nodeOwner();
+      const auto node_number = connectivity.nodeNumber();
+
+      NodeValue<int> node_value_ref{connectivity};
+      parallel_for(
+        connectivity.numberOfNodes(),
+        PUGS_LAMBDA(const NodeId node_id) { node_value_ref[node_id] = node_owner[node_id] + node_number[node_id]; });
+
+      NodeValue<int> node_value{connectivity};
+      parallel_for(
+        connectivity.numberOfNodes(),
+        PUGS_LAMBDA(const NodeId node_id) { node_value[node_id] = parallel::rank() + node_number[node_id]; });
+
+      if (parallel::size() > 1) {
+        REQUIRE(not is_same_item_value(node_value, node_value_ref));
+      }
+
+      Synchronizer synchronizer;
+      synchronizer.synchronize(node_value);
+
+      REQUIRE(is_same_item_value(node_value, node_value_ref));
+    }
+
+    SECTION("synchonize EdgeValue")
+    {
+      const auto edge_owner  = connectivity.edgeOwner();
+      const auto edge_number = connectivity.edgeNumber();
+
+      EdgeValue<int> edge_value_ref{connectivity};
+      parallel_for(
+        connectivity.numberOfEdges(),
+        PUGS_LAMBDA(const EdgeId edge_id) { edge_value_ref[edge_id] = edge_owner[edge_id] + edge_number[edge_id]; });
+
+      EdgeValue<int> edge_value{connectivity};
+      parallel_for(
+        connectivity.numberOfEdges(),
+        PUGS_LAMBDA(const EdgeId edge_id) { edge_value[edge_id] = parallel::rank() + edge_number[edge_id]; });
+
+      if (parallel::size() > 1) {
+        REQUIRE(not is_same_item_value(edge_value, edge_value_ref));
+      }
+
+      Synchronizer synchronizer;
+      synchronizer.synchronize(edge_value);
+
+      REQUIRE(is_same_item_value(edge_value, edge_value_ref));
+    }
+
+    SECTION("synchonize FaceValue")
+    {
+      const auto face_owner  = connectivity.faceOwner();
+      const auto face_number = connectivity.faceNumber();
+
+      FaceValue<int> face_value_ref{connectivity};
+      parallel_for(
+        connectivity.numberOfFaces(),
+        PUGS_LAMBDA(const FaceId face_id) { face_value_ref[face_id] = face_owner[face_id] + face_number[face_id]; });
+
+      FaceValue<int> face_value{connectivity};
+      parallel_for(
+        connectivity.numberOfFaces(),
+        PUGS_LAMBDA(const FaceId face_id) { face_value[face_id] = parallel::rank() + face_number[face_id]; });
+
+      if (parallel::size() > 1) {
+        REQUIRE(not is_same_item_value(face_value, face_value_ref));
+      }
+
+      Synchronizer synchronizer;
+      synchronizer.synchronize(face_value);
+
+      REQUIRE(is_same_item_value(face_value, face_value_ref));
+    }
+
+    SECTION("synchonize CellValue")
+    {
+      const auto cell_owner  = connectivity.cellOwner();
+      const auto cell_number = connectivity.cellNumber();
+
+      CellValue<int> cell_value_ref{connectivity};
+      parallel_for(
+        connectivity.numberOfCells(),
+        PUGS_LAMBDA(const CellId cell_id) { cell_value_ref[cell_id] = cell_owner[cell_id] + cell_number[cell_id]; });
+
+      CellValue<int> cell_value{connectivity};
+      parallel_for(
+        connectivity.numberOfCells(),
+        PUGS_LAMBDA(const CellId cell_id) { cell_value[cell_id] = parallel::rank() + cell_number[cell_id]; });
+
+      if (parallel::size() > 1) {
+        REQUIRE(not is_same_item_value(cell_value, cell_value_ref));
+      }
+
+      Synchronizer synchronizer;
+      synchronizer.synchronize(cell_value);
+
+      REQUIRE(is_same_item_value(cell_value, cell_value_ref));
+    }
+
+    SECTION("synchonize CellArray")
+    {
+      const auto cell_owner  = connectivity.cellOwner();
+      const auto cell_number = connectivity.cellNumber();
+
+      CellArray<int> cell_array_ref{connectivity, 3};
+      parallel_for(
+        connectivity.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+          for (size_t i = 0; i < cell_array_ref.sizeOfArrays(); ++i) {
+            cell_array_ref[cell_id][i] = (i + 1) * cell_owner[cell_id] + i + cell_number[cell_id];
+          }
+        });
+
+      CellArray<int> cell_array{connectivity, 3};
+      parallel_for(
+        connectivity.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+          for (size_t i = 0; i < cell_array.sizeOfArrays(); ++i) {
+            cell_array[cell_id][i] = (i + 1) * parallel::rank() + i + cell_number[cell_id];
+          }
+        });
+
+      if (parallel::size() > 1) {
+        REQUIRE(not is_same_item_array(cell_array, cell_array_ref));
+      }
+
+      Synchronizer synchronizer;
+      synchronizer.synchronize(cell_array);
+
+      REQUIRE(is_same_item_array(cell_array, cell_array_ref));
+    }
+  }
+
+  SECTION("2D")
+  {
+    constexpr size_t Dimension = 2;
+    using ConnectivityType     = Connectivity<Dimension>;
+
+    const ConnectivityType& connectivity = MeshDataBaseForTests::get().hybrid2DMesh()->connectivity();
+
+    SECTION("synchonize NodeValue")
+    {
+      const auto node_owner  = connectivity.nodeOwner();
+      const auto node_number = connectivity.nodeNumber();
+
+      NodeValue<int> node_value_ref{connectivity};
+      parallel_for(
+        connectivity.numberOfNodes(),
+        PUGS_LAMBDA(const NodeId node_id) { node_value_ref[node_id] = node_owner[node_id] + node_number[node_id]; });
+
+      NodeValue<int> node_value{connectivity};
+      parallel_for(
+        connectivity.numberOfNodes(),
+        PUGS_LAMBDA(const NodeId node_id) { node_value[node_id] = parallel::rank() + node_number[node_id]; });
+
+      if (parallel::size() > 1) {
+        REQUIRE(not is_same_item_value(node_value, node_value_ref));
+      }
+
+      Synchronizer synchronizer;
+      synchronizer.synchronize(node_value);
+
+      REQUIRE(is_same_item_value(node_value, node_value_ref));
+    }
+
+    SECTION("synchonize EdgeValue")
+    {
+      const auto edge_owner  = connectivity.edgeOwner();
+      const auto edge_number = connectivity.edgeNumber();
+
+      EdgeValue<int> edge_value_ref{connectivity};
+      parallel_for(
+        connectivity.numberOfEdges(),
+        PUGS_LAMBDA(const EdgeId edge_id) { edge_value_ref[edge_id] = edge_owner[edge_id] + edge_number[edge_id]; });
+
+      EdgeValue<int> edge_value{connectivity};
+      parallel_for(
+        connectivity.numberOfEdges(),
+        PUGS_LAMBDA(const EdgeId edge_id) { edge_value[edge_id] = parallel::rank() + edge_number[edge_id]; });
+
+      if (parallel::size() > 1) {
+        REQUIRE(not is_same_item_value(edge_value, edge_value_ref));
+      }
+
+      Synchronizer synchronizer;
+      synchronizer.synchronize(edge_value);
+
+      REQUIRE(is_same_item_value(edge_value, edge_value_ref));
+    }
+
+    SECTION("synchonize FaceValue")
+    {
+      const auto face_owner  = connectivity.faceOwner();
+      const auto face_number = connectivity.faceNumber();
+
+      FaceValue<int> face_value_ref{connectivity};
+      parallel_for(
+        connectivity.numberOfFaces(),
+        PUGS_LAMBDA(const FaceId face_id) { face_value_ref[face_id] = face_owner[face_id] + face_number[face_id]; });
+
+      FaceValue<int> face_value{connectivity};
+      parallel_for(
+        connectivity.numberOfFaces(),
+        PUGS_LAMBDA(const FaceId face_id) { face_value[face_id] = parallel::rank() + face_number[face_id]; });
+
+      if (parallel::size() > 1) {
+        REQUIRE(not is_same_item_value(face_value, face_value_ref));
+      }
+
+      Synchronizer synchronizer;
+      synchronizer.synchronize(face_value);
+
+      REQUIRE(is_same_item_value(face_value, face_value_ref));
+    }
+
+    SECTION("synchonize CellValue")
+    {
+      const auto cell_owner  = connectivity.cellOwner();
+      const auto cell_number = connectivity.cellNumber();
+
+      CellValue<int> cell_value_ref{connectivity};
+      parallel_for(
+        connectivity.numberOfCells(),
+        PUGS_LAMBDA(const CellId cell_id) { cell_value_ref[cell_id] = cell_owner[cell_id] + cell_number[cell_id]; });
+
+      CellValue<int> cell_value{connectivity};
+      parallel_for(
+        connectivity.numberOfCells(),
+        PUGS_LAMBDA(const CellId cell_id) { cell_value[cell_id] = parallel::rank() + cell_number[cell_id]; });
+
+      if (parallel::size() > 1) {
+        REQUIRE(not is_same_item_value(cell_value, cell_value_ref));
+      }
+
+      Synchronizer synchronizer;
+      synchronizer.synchronize(cell_value);
+
+      REQUIRE(is_same_item_value(cell_value, cell_value_ref));
+    }
+
+    SECTION("synchonize CellArray")
+    {
+      const auto cell_owner  = connectivity.cellOwner();
+      const auto cell_number = connectivity.cellNumber();
+
+      CellArray<int> cell_array_ref{connectivity, 3};
+      parallel_for(
+        connectivity.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+          for (size_t i = 0; i < cell_array_ref.sizeOfArrays(); ++i) {
+            cell_array_ref[cell_id][i] = (i + 1) * cell_owner[cell_id] + i + cell_number[cell_id];
+          }
+        });
+
+      CellArray<int> cell_array{connectivity, 3};
+      parallel_for(
+        connectivity.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+          for (size_t i = 0; i < cell_array.sizeOfArrays(); ++i) {
+            cell_array[cell_id][i] = (i + 1) * parallel::rank() + i + cell_number[cell_id];
+          }
+        });
+
+      if (parallel::size() > 1) {
+        REQUIRE(not is_same_item_array(cell_array, cell_array_ref));
+      }
+
+      Synchronizer synchronizer;
+      synchronizer.synchronize(cell_array);
+
+      REQUIRE(is_same_item_array(cell_array, cell_array_ref));
+    }
+  }
+
+  SECTION("3D")
+  {
+    constexpr size_t Dimension = 3;
+    using ConnectivityType     = Connectivity<Dimension>;
+
+    const ConnectivityType& connectivity = MeshDataBaseForTests::get().hybrid3DMesh()->connectivity();
+
+    SECTION("synchonize NodeValue")
+    {
+      const auto node_owner  = connectivity.nodeOwner();
+      const auto node_number = connectivity.nodeNumber();
+
+      NodeValue<int> node_value_ref{connectivity};
+      parallel_for(
+        connectivity.numberOfNodes(),
+        PUGS_LAMBDA(const NodeId node_id) { node_value_ref[node_id] = node_owner[node_id] + node_number[node_id]; });
+
+      NodeValue<int> node_value{connectivity};
+      parallel_for(
+        connectivity.numberOfNodes(),
+        PUGS_LAMBDA(const NodeId node_id) { node_value[node_id] = parallel::rank() + node_number[node_id]; });
+
+      if (parallel::size() > 1) {
+        REQUIRE(not is_same_item_value(node_value, node_value_ref));
+      }
+
+      Synchronizer synchronizer;
+      synchronizer.synchronize(node_value);
+
+      REQUIRE(is_same_item_value(node_value, node_value_ref));
+    }
+
+    SECTION("synchonize EdgeValue")
+    {
+      const auto edge_owner  = connectivity.edgeOwner();
+      const auto edge_number = connectivity.edgeNumber();
+
+      EdgeValue<int> edge_value_ref{connectivity};
+      parallel_for(
+        connectivity.numberOfEdges(),
+        PUGS_LAMBDA(const EdgeId edge_id) { edge_value_ref[edge_id] = edge_owner[edge_id] + edge_number[edge_id]; });
+
+      EdgeValue<int> edge_value{connectivity};
+      parallel_for(
+        connectivity.numberOfEdges(),
+        PUGS_LAMBDA(const EdgeId edge_id) { edge_value[edge_id] = parallel::rank() + edge_number[edge_id]; });
+
+      if (parallel::size() > 1) {
+        REQUIRE(not is_same_item_value(edge_value, edge_value_ref));
+      }
+
+      Synchronizer synchronizer;
+      synchronizer.synchronize(edge_value);
+
+      REQUIRE(is_same_item_value(edge_value, edge_value_ref));
+    }
+
+    SECTION("synchonize FaceValue")
+    {
+      const auto face_owner  = connectivity.faceOwner();
+      const auto face_number = connectivity.faceNumber();
+
+      FaceValue<int> face_value_ref{connectivity};
+      parallel_for(
+        connectivity.numberOfFaces(),
+        PUGS_LAMBDA(const FaceId face_id) { face_value_ref[face_id] = face_owner[face_id] + face_number[face_id]; });
+
+      FaceValue<int> face_value{connectivity};
+      parallel_for(
+        connectivity.numberOfFaces(),
+        PUGS_LAMBDA(const FaceId face_id) { face_value[face_id] = parallel::rank() + face_number[face_id]; });
+
+      if (parallel::size() > 1) {
+        REQUIRE(not is_same_item_value(face_value, face_value_ref));
+      }
+
+      Synchronizer synchronizer;
+      synchronizer.synchronize(face_value);
+
+      REQUIRE(is_same_item_value(face_value, face_value_ref));
+    }
+
+    SECTION("synchonize CellValue")
+    {
+      const auto cell_owner  = connectivity.cellOwner();
+      const auto cell_number = connectivity.cellNumber();
+
+      CellValue<int> cell_value_ref{connectivity};
+      parallel_for(
+        connectivity.numberOfCells(),
+        PUGS_LAMBDA(const CellId cell_id) { cell_value_ref[cell_id] = cell_owner[cell_id] + cell_number[cell_id]; });
+
+      CellValue<int> cell_value{connectivity};
+      parallel_for(
+        connectivity.numberOfCells(),
+        PUGS_LAMBDA(const CellId cell_id) { cell_value[cell_id] = parallel::rank() + cell_number[cell_id]; });
+
+      if (parallel::size() > 1) {
+        REQUIRE(not is_same_item_value(cell_value, cell_value_ref));
+      }
+
+      Synchronizer synchronizer;
+      synchronizer.synchronize(cell_value);
+
+      REQUIRE(is_same_item_value(cell_value, cell_value_ref));
+    }
+
+    SECTION("synchonize CellArray")
+    {
+      const auto cell_owner  = connectivity.cellOwner();
+      const auto cell_number = connectivity.cellNumber();
+
+      CellArray<int> cell_array_ref{connectivity, 3};
+      parallel_for(
+        connectivity.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+          for (size_t i = 0; i < cell_array_ref.sizeOfArrays(); ++i) {
+            cell_array_ref[cell_id][i] = (i + 1) * cell_owner[cell_id] + i + cell_number[cell_id];
+          }
+        });
+
+      CellArray<int> cell_array{connectivity, 3};
+      parallel_for(
+        connectivity.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+          for (size_t i = 0; i < cell_array.sizeOfArrays(); ++i) {
+            cell_array[cell_id][i] = (i + 1) * parallel::rank() + i + cell_number[cell_id];
+          }
+        });
+
+      if (parallel::size() > 1) {
+        REQUIRE(not is_same_item_array(cell_array, cell_array_ref));
+      }
+
+      Synchronizer synchronizer;
+      synchronizer.synchronize(cell_array);
+
+      REQUIRE(is_same_item_array(cell_array, cell_array_ref));
+    }
+  }
+}