diff --git a/src/language/utils/IntegrateOnCells.hpp b/src/language/utils/IntegrateOnCells.hpp
index 8b6008d73672a24ee5473369250f00ce05342d7e..82a4c3242aff86d2b719dde80035dd72cf8fef89 100644
--- a/src/language/utils/IntegrateOnCells.hpp
+++ b/src/language/utils/IntegrateOnCells.hpp
@@ -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:
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 18d5656925cae0cb5d0cba949bebb2ce74ee4092..cf294f9ecb8ea51c6c6ec140915d0fb6107fc197 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -98,6 +98,7 @@ add_executable (unit_tests
   test_GaussQuadratureDescriptor.cpp
   test_IfProcessor.cpp
   test_IncDecExpressionProcessor.cpp
+  test_IntegrateOnCells.cpp
   test_INodeProcessor.cpp
   test_ItemId.cpp
   test_ItemType.cpp
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));
+            }
+          }
+        }
+      }
+    }
+  }
+}