diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index 3d097d30761fda305ea6490d74da0fd19438a972..b0f1a2313a6db1a5164ea46862d4fd84d0ddd6a1 100644
--- a/src/language/modules/MeshModule.cpp
+++ b/src/language/modules/MeshModule.cpp
@@ -2,8 +2,11 @@
 
 #include <algebra/TinyVector.hpp>
 #include <language/node_processor/ExecutionPolicy.hpp>
+#include <language/utils/BinaryOperatorProcessorBuilder.hpp>
 #include <language/utils/BuiltinFunctionEmbedder.hpp>
 #include <language/utils/FunctionTable.hpp>
+#include <language/utils/OStream.hpp>
+#include <language/utils/OperatorRepository.hpp>
 #include <language/utils/SymbolTable.hpp>
 #include <language/utils/TypeDescriptor.hpp>
 #include <mesh/CartesianMeshBuilder.hpp>
@@ -250,4 +253,10 @@ MeshModule::MeshModule()
 
 void
 MeshModule::registerOperators() const
-{}
+{
+  OperatorRepository& repository = OperatorRepository::instance();
+
+  repository.addBinaryOperator<language::shift_left_op>(
+    std::make_shared<BinaryOperatorProcessorBuilder<language::shift_left_op, std::shared_ptr<const OStream>,
+                                                    std::shared_ptr<const OStream>, std::shared_ptr<const IMesh>>>());
+}
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index c724e1674415b1a2be0c8848f0a9cf923bd24363..5d397223996c592687707711150705056e57879c 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -11,6 +11,7 @@
 #include <language/utils/TypeDescriptor.hpp>
 #include <mesh/Connectivity.hpp>
 #include <mesh/IBoundaryDescriptor.hpp>
+#include <mesh/IZoneDescriptor.hpp>
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
@@ -90,6 +91,27 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this->_addBuiltinFunction("integrate",
+                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+                              const IDiscreteFunction>(std::shared_ptr<const IMesh>,
+                                                       const std::vector<std::shared_ptr<const IZoneDescriptor>>&,
+                                                       std::shared_ptr<const IQuadratureDescriptor>,
+                                                       std::shared_ptr<const IDiscreteFunctionDescriptor>,
+                                                       const std::vector<FunctionSymbolId>&)>>(
+                              [](std::shared_ptr<const IMesh> mesh,
+                                 const std::vector<std::shared_ptr<const IZoneDescriptor>>& integration_zone_list,
+                                 std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor,
+                                 std::shared_ptr<const IDiscreteFunctionDescriptor> discrete_function_descriptor,
+                                 const std::vector<FunctionSymbolId>& function_id_list)
+                                -> std::shared_ptr<const IDiscreteFunction> {
+                                return DiscreteFunctionVectorIntegrator{mesh, integration_zone_list,
+                                                                        quadrature_descriptor,
+                                                                        discrete_function_descriptor, function_id_list}
+                                  .integrate();
+                              }
+
+                              ));
+
   this
     ->_addBuiltinFunction("integrate",
                           std::make_shared<BuiltinFunctionEmbedder<
@@ -109,6 +131,20 @@ SchemeModule::SchemeModule()
 
                             ));
 
+  this->_addBuiltinFunction(
+    "integrate",
+    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+      const IDiscreteFunction>(std::shared_ptr<const IMesh>, const std::vector<std::shared_ptr<const IZoneDescriptor>>&,
+                               std::shared_ptr<const IQuadratureDescriptor>, const FunctionSymbolId&)>>(
+      [](std::shared_ptr<const IMesh> mesh,
+         const std::vector<std::shared_ptr<const IZoneDescriptor>>& integration_zone_list,
+         std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor,
+         const FunctionSymbolId& function_id) -> std::shared_ptr<const IDiscreteFunction> {
+        return DiscreteFunctionIntegrator{mesh, integration_zone_list, quadrature_descriptor, function_id}.integrate();
+      }
+
+      ));
+
   this->_addBuiltinFunction("integrate",
                             std::make_shared<BuiltinFunctionEmbedder<
                               std::shared_ptr<const IDiscreteFunction>(std::shared_ptr<const IMesh>,
@@ -135,6 +171,24 @@ SchemeModule::SchemeModule()
 
       ));
 
+  this->_addBuiltinFunction("interpolate",
+                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+                              const IDiscreteFunction>(std::shared_ptr<const IMesh>,
+                                                       const std::vector<std::shared_ptr<const IZoneDescriptor>>&,
+                                                       std::shared_ptr<const IDiscreteFunctionDescriptor>,
+                                                       const std::vector<FunctionSymbolId>&)>>(
+                              [](std::shared_ptr<const IMesh> mesh,
+                                 const std::vector<std::shared_ptr<const IZoneDescriptor>>& interpolation_zone_list,
+                                 std::shared_ptr<const IDiscreteFunctionDescriptor> discrete_function_descriptor,
+                                 const std::vector<FunctionSymbolId>& function_id_list)
+                                -> std::shared_ptr<const IDiscreteFunction> {
+                                return DiscreteFunctionVectorInterpoler{mesh, interpolation_zone_list,
+                                                                        discrete_function_descriptor, function_id_list}
+                                  .interpolate();
+                              }
+
+                              ));
+
   this->_addBuiltinFunction(
     "interpolate",
     std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
@@ -158,6 +212,71 @@ SchemeModule::SchemeModule()
 
       ));
 
+  this->_addBuiltinFunction("interpolate",
+                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+                              const IDiscreteFunction>(std::shared_ptr<const IMesh>,
+                                                       const std::vector<std::shared_ptr<const IZoneDescriptor>>&,
+                                                       std::shared_ptr<const IDiscreteFunctionDescriptor>,
+                                                       const FunctionSymbolId&)>>(
+                              [](std::shared_ptr<const IMesh> mesh,
+                                 const std::vector<std::shared_ptr<const IZoneDescriptor>>& interpolation_zone_list,
+                                 std::shared_ptr<const IDiscreteFunctionDescriptor> discrete_function_descriptor,
+                                 const FunctionSymbolId& function_id) -> std::shared_ptr<const IDiscreteFunction> {
+                                switch (discrete_function_descriptor->type()) {
+                                case DiscreteFunctionType::P0: {
+                                  return DiscreteFunctionInterpoler{mesh, interpolation_zone_list,
+                                                                    discrete_function_descriptor, function_id}
+                                    .interpolate();
+                                }
+                                case DiscreteFunctionType::P0Vector: {
+                                  return DiscreteFunctionVectorInterpoler{mesh,
+                                                                          interpolation_zone_list,
+                                                                          discrete_function_descriptor,
+                                                                          {function_id}}
+                                    .interpolate();
+                                }
+                                default: {
+                                  throw NormalError("invalid function descriptor type");
+                                }
+                                }
+                              }
+
+                              ));
+
+  this
+    ->_addBuiltinFunction("interpolate",
+                          std::make_shared<BuiltinFunctionEmbedder<
+                            std::shared_ptr<const IDiscreteFunction>(std::shared_ptr<const IMesh>,
+                                                                     std::shared_ptr<const IZoneDescriptor>,
+                                                                     std::shared_ptr<const IDiscreteFunctionDescriptor>,
+                                                                     const FunctionSymbolId&)>>(
+                            [](std::shared_ptr<const IMesh> mesh,
+                               std::shared_ptr<const IZoneDescriptor> interpolation_zone,
+                               std::shared_ptr<const IDiscreteFunctionDescriptor> discrete_function_descriptor,
+                               const FunctionSymbolId& function_id) -> std::shared_ptr<const IDiscreteFunction> {
+                              switch (discrete_function_descriptor->type()) {
+                              case DiscreteFunctionType::P0: {
+                                return DiscreteFunctionInterpoler{mesh,
+                                                                  {interpolation_zone},
+                                                                  discrete_function_descriptor,
+                                                                  function_id}
+                                  .interpolate();
+                              }
+                              case DiscreteFunctionType::P0Vector: {
+                                return DiscreteFunctionVectorInterpoler{mesh,
+                                                                        {interpolation_zone},
+                                                                        discrete_function_descriptor,
+                                                                        {function_id}}
+                                  .interpolate();
+                              }
+                              default: {
+                                throw NormalError("invalid function descriptor type");
+                              }
+                              }
+                            }
+
+                            ));
+
   this->_addBuiltinFunction("randomizeMesh",
                             std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
                               const IMesh>(std::shared_ptr<const IMesh>,
diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index a819af495a1dbb358a2705231e3ff30f2bcbed02..7b81eeb0326e7cd5bc5df2b9c6d45b642dcb5e56 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -241,6 +241,63 @@ Connectivity<Dimension>::_buildIsBoundaryNode() const
   }
 }
 
+template <ItemType item_type, size_t Dimension>
+inline void
+_printReference(std::ostream& os, const Connectivity<Dimension>& connectivity)
+{
+  auto count_all_items = [](const auto& item_is_owned) -> size_t {
+    using ItemId  = typename std::decay_t<decltype(item_is_owned)>::index_type;
+    size_t number = 0;
+    for (ItemId item_id = 0; item_id < item_is_owned.numberOfItems(); ++item_id) {
+      number += item_is_owned[item_id];
+    }
+    return parallel::allReduceSum(number);
+  };
+
+  auto count_zone_items = [](const auto& item_is_owned, const auto& item_list) -> size_t {
+    size_t number = 0;
+    for (size_t i_item = 0; i_item < item_list.size(); ++i_item) {
+      number += item_is_owned[item_list[i_item]];
+    }
+    return parallel::allReduceSum(number);
+  };
+
+  os << "- number of " << itemName(item_type) << "s: " << rang::fgB::yellow
+     << count_all_items(connectivity.template isOwned<item_type>()) << rang::style::reset << '\n';
+  os << "  " << rang::fgB::yellow << connectivity.template numberOfRefItemList<item_type>() << rang::style::reset
+     << " references\n";
+  if (connectivity.template numberOfRefItemList<item_type>() > 0) {
+    for (size_t i_ref_item = 0; i_ref_item < connectivity.template numberOfRefItemList<item_type>(); ++i_ref_item) {
+      auto ref_item_list = connectivity.template refItemList<item_type>(i_ref_item);
+      os << "  - " << rang::fgB::green << ref_item_list.refId().tagName() << rang::style::reset << " ("
+         << rang::fgB::green << ref_item_list.refId().tagNumber() << rang::style::reset << ") number "
+         << rang::fgB::yellow << count_zone_items(connectivity.template isOwned<item_type>(), ref_item_list.list())
+         << rang::style::reset << '\n';
+    }
+  }
+}
+
+template <size_t Dimension>
+std::ostream&
+Connectivity<Dimension>::_write(std::ostream& os) const
+{
+  os << "connectivity of dimension " << Dimension << '\n';
+  _printReference<ItemType::cell>(os, *this);
+  if constexpr (Dimension > 1) {
+    _printReference<ItemType::face>(os, *this);
+  }
+  if constexpr (Dimension > 2) {
+    _printReference<ItemType::edge>(os, *this);
+  }
+  _printReference<ItemType::node>(os, *this);
+
+  return os;
+}
+
+template std::ostream& Connectivity<1>::_write(std::ostream&) const;
+template std::ostream& Connectivity<2>::_write(std::ostream&) const;
+template std::ostream& Connectivity<3>::_write(std::ostream&) const;
+
 template void Connectivity<1>::_buildIsBoundaryFace() const;
 template void Connectivity<2>::_buildIsBoundaryFace() const;
 template void Connectivity<3>::_buildIsBoundaryFace() const;
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index 86b6606c7f24c2271569f07cbb5fe9d10a8b8fe3..b0279507ee41e6f6d3433bda2b5dd7beaf406265 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -29,6 +29,9 @@ class ConnectivityDescriptor;
 template <size_t Dim>
 class Connectivity final : public IConnectivity
 {
+ private:
+  std::ostream& _write(std::ostream&) const final;
+
  public:
   PUGS_INLINE
   static std::shared_ptr<Connectivity<Dim>> build(const ConnectivityDescriptor&);
diff --git a/src/mesh/IConnectivity.hpp b/src/mesh/IConnectivity.hpp
index def3475c386a3e3b321ea640cdf378456f39fac6..3ea34ab6e4fcf6935b456ed0c9cec1b65b79d44f 100644
--- a/src/mesh/IConnectivity.hpp
+++ b/src/mesh/IConnectivity.hpp
@@ -9,6 +9,9 @@
 
 class IConnectivity : public std::enable_shared_from_this<IConnectivity>
 {
+ protected:
+  virtual std::ostream& _write(std::ostream&) const = 0;
+
  public:
   template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
   friend class SubItemValuePerItem;
@@ -16,6 +19,12 @@ class IConnectivity : public std::enable_shared_from_this<IConnectivity>
   template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
   friend class SubItemArrayPerItem;
 
+  friend std::ostream&
+  operator<<(std::ostream& os, const IConnectivity& connectivity)
+  {
+    return connectivity._write(os);
+  }
+
   virtual const ConnectivityMatrix& getMatrix(const ItemType& item_type_0, const ItemType& item_type_1) const = 0;
 
   virtual size_t dimension() const = 0;
diff --git a/src/mesh/IMesh.hpp b/src/mesh/IMesh.hpp
index cb524849aaf0a7f994d8fb34fe2e97508c5e4f8a..350f0f2e660df2ae7070e3435b25eb53d19f61b5 100644
--- a/src/mesh/IMesh.hpp
+++ b/src/mesh/IMesh.hpp
@@ -2,12 +2,22 @@
 #define I_MESH_HPP
 
 #include <cstddef>
+#include <iostream>
 
 class IMesh
 {
+ protected:
+  virtual std::ostream& _write(std::ostream&) const = 0;
+
  public:
   virtual size_t dimension() const = 0;
 
+  friend std::ostream&
+  operator<<(std::ostream& os, const IMesh& mesh)
+  {
+    return mesh._write(os);
+  }
+
   IMesh(const IMesh&) = delete;
   IMesh(IMesh&&)      = delete;
 
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index e5817a2a9394b96bea04e78d061e166d3f0fbc02..52cdd8bbfb834c98aaf6242cefa6b0d726e06ef2 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -20,6 +20,12 @@ class Mesh final : public IMesh
   const std::shared_ptr<const Connectivity> m_connectivity;
   const NodeValue<const Rd> m_xr;
 
+  std::ostream&
+  _write(std::ostream& os) const final
+  {
+    return os << *m_connectivity;
+  }
+
  public:
   PUGS_INLINE
   size_t
diff --git a/src/mesh/MeshCellZone.cpp b/src/mesh/MeshCellZone.cpp
index 36418e2b510dae6a15d4f30229c5330ff09af2c9..240d8ea4f1be27e879d0455d312f12f67545e7fb 100644
--- a/src/mesh/MeshCellZone.cpp
+++ b/src/mesh/MeshCellZone.cpp
@@ -26,7 +26,17 @@ getMeshCellZone(const Mesh<Connectivity<Dimension>>& mesh, const IZoneDescriptor
   }
 
   std::ostringstream ost;
-  ost << "cannot find zone with name " << rang::fgB::red << zone_descriptor << rang::style::reset;
+  ost << "cannot find zone with name " << rang::fgB::red << zone_descriptor << rang::style::reset << '\n';
+  ost << "The mesh contains " << mesh.connectivity().template numberOfRefItemList<ItemType::cell>() << " zones: ";
+  for (size_t i_ref_cell_list = 0; i_ref_cell_list < mesh.connectivity().template numberOfRefItemList<ItemType::cell>();
+       ++i_ref_cell_list) {
+    const auto& ref_cell_list = mesh.connectivity().template refItemList<ItemType::cell>(i_ref_cell_list);
+    const RefId& ref          = ref_cell_list.refId();
+    if (i_ref_cell_list > 0) {
+      ost << ", ";
+    }
+    ost << rang::fgB::yellow << ref << rang::style::reset;
+  }
 
   throw NormalError(ost.str());
 }
diff --git a/src/scheme/DiscreteFunctionIntegrator.cpp b/src/scheme/DiscreteFunctionIntegrator.cpp
index 91cf921df2764ad13f4470687fc4cc5f4d6bced1..93c836461ce14ae05b51b879cf90525f73a793eb 100644
--- a/src/scheme/DiscreteFunctionIntegrator.cpp
+++ b/src/scheme/DiscreteFunctionIntegrator.cpp
@@ -1,35 +1,94 @@
 #include <scheme/DiscreteFunctionIntegrator.hpp>
 
 #include <language/utils/IntegrateCellValue.hpp>
+#include <mesh/MeshCellZone.hpp>
 #include <scheme/DiscreteFunctionP0.hpp>
 #include <utils/Exceptions.hpp>
 
 template <size_t Dimension, typename DataType, typename ValueType>
 std::shared_ptr<IDiscreteFunction>
-DiscreteFunctionIntegrator::_integrate() const
+DiscreteFunctionIntegrator::_integrateOnZoneList() const
 {
-  using MeshType       = Mesh<Connectivity<Dimension>>;
-  std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(m_mesh);
+  static_assert(std::is_convertible_v<DataType, ValueType>);
+  Assert(m_zone_list.size() > 0, "no zone list provided");
+
+  using MeshType = Mesh<Connectivity<Dimension>>;
 
-  if constexpr (std::is_same_v<DataType, ValueType>) {
-    return std::make_shared<
-      DiscreteFunctionP0<Dimension, ValueType>>(mesh,
-                                                IntegrateCellValue<DataType(TinyVector<Dimension>)>::template integrate<
-                                                  MeshType>(m_function_id, *m_quadrature_descriptor, *mesh));
+  std::shared_ptr p_mesh = std::dynamic_pointer_cast<const MeshType>(m_mesh);
+
+  CellValue<bool> is_in_zone{p_mesh->connectivity()};
+  is_in_zone.fill(false);
+
+  size_t number_of_cells = 0;
+  for (const auto& zone : m_zone_list) {
+    auto mesh_cell_zone   = getMeshCellZone(*p_mesh, *zone);
+    const auto& cell_list = mesh_cell_zone.cellList();
+    for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+      const CellId cell_id = cell_list[i_cell];
+      if (is_in_zone[cell_id]) {
+        std::ostringstream os;
+        os << "cell " << cell_id << " (number " << p_mesh->connectivity().cellNumber()[cell_id]
+           << ") is present multiple times in zone list";
+        throw NormalError(os.str());
+      }
+      ++number_of_cells;
+      is_in_zone[cell_id] = true;
+    }
+  }
+
+  Array<CellId> zone_cell_list{number_of_cells};
+  {
+    size_t i_cell = 0;
+    for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+      if (is_in_zone[cell_id]) {
+        zone_cell_list[i_cell++] = cell_id;
+      }
+    }
+  }
+
+  Array<ValueType> array =
+    IntegrateCellValue<ValueType(TinyVector<Dimension>)>::template integrate<MeshType>(m_function_id,
+                                                                                       *m_quadrature_descriptor,
+                                                                                       *p_mesh, zone_cell_list);
+
+  CellValue<ValueType> cell_value{p_mesh->connectivity()};
+  if constexpr (is_tiny_vector_v<ValueType> or is_tiny_matrix_v<ValueType>) {
+    cell_value.fill(zero);
   } else {
-    static_assert(std::is_convertible_v<DataType, ValueType>);
+    cell_value.fill(0);
+  }
 
-    CellValue<DataType> cell_data =
-      IntegrateCellValue<DataType(TinyVector<Dimension>)>::template integrate<MeshType>(m_function_id,
-                                                                                        *m_quadrature_descriptor,
-                                                                                        *mesh);
+  parallel_for(
+    zone_cell_list.size(), PUGS_LAMBDA(const size_t i_cell) { cell_value[zone_cell_list[i_cell]] = array[i_cell]; });
 
-    CellValue<ValueType> cell_value{mesh->connectivity()};
+  return std::make_shared<DiscreteFunctionP0<Dimension, ValueType>>(p_mesh, cell_value);
+}
 
-    parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { cell_value[cell_id] = cell_data[cell_id]; });
+template <size_t Dimension, typename DataType, typename ValueType>
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionIntegrator::_integrateGlobally() const
+{
+  Assert(m_zone_list.size() == 0, "invalid call when zones are defined");
+
+  using MeshType       = Mesh<Connectivity<Dimension>>;
+  std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(m_mesh);
+
+  static_assert(std::is_convertible_v<DataType, ValueType>);
 
-    return std::make_shared<DiscreteFunctionP0<Dimension, ValueType>>(mesh, cell_value);
+  return std::make_shared<
+    DiscreteFunctionP0<Dimension, ValueType>>(mesh,
+                                              IntegrateCellValue<ValueType(TinyVector<Dimension>)>::template integrate<
+                                                MeshType>(m_function_id, *m_quadrature_descriptor, *mesh));
+}
+
+template <size_t Dimension, typename DataType, typename ValueType>
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionIntegrator::_integrate() const
+{
+  if (m_zone_list.size() == 0) {
+    return this->_integrateGlobally<Dimension, DataType, ValueType>();
+  } else {
+    return this->_integrateOnZoneList<Dimension, DataType, ValueType>();
   }
 }
 
diff --git a/src/scheme/DiscreteFunctionIntegrator.hpp b/src/scheme/DiscreteFunctionIntegrator.hpp
index 1784f891eaa6fd7471b17006b3c9f277c31e51dd..21a461d54561d67f1f745d3257ec76ee514d264a 100644
--- a/src/scheme/DiscreteFunctionIntegrator.hpp
+++ b/src/scheme/DiscreteFunctionIntegrator.hpp
@@ -4,6 +4,7 @@
 #include <analysis/IQuadratureDescriptor.hpp>
 #include <language/utils/FunctionSymbolId.hpp>
 #include <mesh/IMesh.hpp>
+#include <mesh/IZoneDescriptor.hpp>
 #include <scheme/IDiscreteFunction.hpp>
 
 #include <memory>
@@ -12,9 +13,16 @@ class DiscreteFunctionIntegrator
 {
  private:
   std::shared_ptr<const IMesh> m_mesh;
+  const std::vector<std::shared_ptr<const IZoneDescriptor>> m_zone_list;
   std::shared_ptr<const IQuadratureDescriptor> m_quadrature_descriptor;
   const FunctionSymbolId m_function_id;
 
+  template <size_t Dimension, typename DataType, typename ValueType = DataType>
+  std::shared_ptr<IDiscreteFunction> _integrateOnZoneList() const;
+
+  template <size_t Dimension, typename DataType, typename ValueType = DataType>
+  std::shared_ptr<IDiscreteFunction> _integrateGlobally() const;
+
   template <size_t Dimension, typename DataType, typename ValueType = DataType>
   std::shared_ptr<IDiscreteFunction> _integrate() const;
 
@@ -30,6 +38,13 @@ class DiscreteFunctionIntegrator
     : m_mesh{mesh}, m_quadrature_descriptor{quadrature_descriptor}, m_function_id{function_id}
   {}
 
+  DiscreteFunctionIntegrator(const std::shared_ptr<const IMesh>& mesh,
+                             const std::vector<std::shared_ptr<const IZoneDescriptor>>& zone_list,
+                             const std::shared_ptr<const IQuadratureDescriptor>& quadrature_descriptor,
+                             const FunctionSymbolId& function_id)
+    : m_mesh{mesh}, m_zone_list{zone_list}, m_quadrature_descriptor{quadrature_descriptor}, m_function_id{function_id}
+  {}
+
   DiscreteFunctionIntegrator(const DiscreteFunctionIntegrator&) = delete;
   DiscreteFunctionIntegrator(DiscreteFunctionIntegrator&&)      = delete;
 
diff --git a/src/scheme/DiscreteFunctionInterpoler.cpp b/src/scheme/DiscreteFunctionInterpoler.cpp
index c39b6ee9d7020b1ffee3582ee36bd461a0c03435..daf7646ee266c7c874c17b6b7403388fc33c3ce6 100644
--- a/src/scheme/DiscreteFunctionInterpoler.cpp
+++ b/src/scheme/DiscreteFunctionInterpoler.cpp
@@ -1,20 +1,82 @@
 #include <scheme/DiscreteFunctionInterpoler.hpp>
 
 #include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/MeshCellZone.hpp>
 #include <scheme/DiscreteFunctionP0.hpp>
 #include <utils/Exceptions.hpp>
 
 template <size_t Dimension, typename DataType, typename ValueType>
 std::shared_ptr<IDiscreteFunction>
-DiscreteFunctionInterpoler::_interpolate() const
+DiscreteFunctionInterpoler::_interpolateOnZoneList() const
+{
+  static_assert(std::is_convertible_v<DataType, ValueType>);
+  Assert(m_zone_list.size() > 0, "no zone list provided");
+
+  std::shared_ptr p_mesh  = std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(m_mesh);
+  using MeshDataType      = MeshData<Dimension>;
+  MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);
+
+  CellValue<bool> is_in_zone{p_mesh->connectivity()};
+  is_in_zone.fill(false);
+
+  size_t number_of_cells = 0;
+  for (const auto& zone : m_zone_list) {
+    auto mesh_cell_zone   = getMeshCellZone(*p_mesh, *zone);
+    const auto& cell_list = mesh_cell_zone.cellList();
+    for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+      const CellId cell_id = cell_list[i_cell];
+      if (is_in_zone[cell_id]) {
+        std::ostringstream os;
+        os << "cell " << cell_id << " (number " << p_mesh->connectivity().cellNumber()[cell_id]
+           << ") is present multiple times in zone list";
+        throw NormalError(os.str());
+      }
+      ++number_of_cells;
+      is_in_zone[cell_id] = true;
+    }
+  }
+
+  Array<CellId> zone_cell_list{number_of_cells};
+  {
+    size_t i_cell = 0;
+    for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+      if (is_in_zone[cell_id]) {
+        zone_cell_list[i_cell++] = cell_id;
+      }
+    }
+  }
+
+  Array<const DataType> array =
+    InterpolateItemValue<DataType(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(m_function_id,
+                                                                                                mesh_data.xj(),
+                                                                                                zone_cell_list);
+
+  CellValue<ValueType> cell_value{p_mesh->connectivity()};
+  if constexpr (is_tiny_vector_v<ValueType> or is_tiny_matrix_v<ValueType>) {
+    cell_value.fill(zero);
+  } else {
+    cell_value.fill(0);
+  }
+
+  parallel_for(
+    zone_cell_list.size(), PUGS_LAMBDA(const size_t i_cell) { cell_value[zone_cell_list[i_cell]] = array[i_cell]; });
+
+  return std::make_shared<DiscreteFunctionP0<Dimension, ValueType>>(p_mesh, cell_value);
+}
+
+template <size_t Dimension, typename DataType, typename ValueType>
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionInterpoler::_interpolateGlobally() const
 {
-  std::shared_ptr mesh    = std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(m_mesh);
+  Assert(m_zone_list.size() == 0, "invalid call when zones are defined");
+
+  std::shared_ptr p_mesh  = std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(m_mesh);
   using MeshDataType      = MeshData<Dimension>;
-  MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+  MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);
 
   if constexpr (std::is_same_v<DataType, ValueType>) {
     return std::make_shared<
-      DiscreteFunctionP0<Dimension, ValueType>>(mesh,
+      DiscreteFunctionP0<Dimension, ValueType>>(p_mesh,
                                                 InterpolateItemValue<DataType(TinyVector<Dimension>)>::
                                                   template interpolate<ItemType::cell>(m_function_id, mesh_data.xj()));
   } else {
@@ -24,12 +86,23 @@ DiscreteFunctionInterpoler::_interpolate() const
       InterpolateItemValue<DataType(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(m_function_id,
                                                                                                   mesh_data.xj());
 
-    CellValue<ValueType> cell_value{mesh->connectivity()};
+    CellValue<ValueType> cell_value{p_mesh->connectivity()};
 
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { cell_value[cell_id] = cell_data[cell_id]; });
+      p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { cell_value[cell_id] = cell_data[cell_id]; });
 
-    return std::make_shared<DiscreteFunctionP0<Dimension, ValueType>>(mesh, cell_value);
+    return std::make_shared<DiscreteFunctionP0<Dimension, ValueType>>(p_mesh, cell_value);
+  }
+}
+
+template <size_t Dimension, typename DataType, typename ValueType>
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionInterpoler::_interpolate() const
+{
+  if (m_zone_list.size() == 0) {
+    return this->_interpolateGlobally<Dimension, DataType, ValueType>();
+  } else {
+    return this->_interpolateOnZoneList<Dimension, DataType, ValueType>();
   }
 }
 
diff --git a/src/scheme/DiscreteFunctionInterpoler.hpp b/src/scheme/DiscreteFunctionInterpoler.hpp
index e5724e3fec3ea75098772827537a1267f68328b7..295c2fd9a540269bf484f485716131b1346964e1 100644
--- a/src/scheme/DiscreteFunctionInterpoler.hpp
+++ b/src/scheme/DiscreteFunctionInterpoler.hpp
@@ -3,6 +3,7 @@
 
 #include <language/utils/FunctionSymbolId.hpp>
 #include <mesh/IMesh.hpp>
+#include <mesh/IZoneDescriptor.hpp>
 #include <scheme/IDiscreteFunction.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 
@@ -12,9 +13,16 @@ class DiscreteFunctionInterpoler
 {
  private:
   std::shared_ptr<const IMesh> m_mesh;
+  const std::vector<std::shared_ptr<const IZoneDescriptor>> m_zone_list;
   std::shared_ptr<const IDiscreteFunctionDescriptor> m_discrete_function_descriptor;
   const FunctionSymbolId m_function_id;
 
+  template <size_t Dimension, typename DataType, typename ValueType = DataType>
+  std::shared_ptr<IDiscreteFunction> _interpolateOnZoneList() const;
+
+  template <size_t Dimension, typename DataType, typename ValueType = DataType>
+  std::shared_ptr<IDiscreteFunction> _interpolateGlobally() const;
+
   template <size_t Dimension, typename DataType, typename ValueType = DataType>
   std::shared_ptr<IDiscreteFunction> _interpolate() const;
 
@@ -30,6 +38,16 @@ class DiscreteFunctionInterpoler
     : m_mesh{mesh}, m_discrete_function_descriptor{discrete_function_descriptor}, m_function_id{function_id}
   {}
 
+  DiscreteFunctionInterpoler(const std::shared_ptr<const IMesh>& mesh,
+                             const std::vector<std::shared_ptr<const IZoneDescriptor>>& zone_list,
+                             const std::shared_ptr<const IDiscreteFunctionDescriptor>& discrete_function_descriptor,
+                             const FunctionSymbolId& function_id)
+    : m_mesh{mesh},
+      m_zone_list{zone_list},
+      m_discrete_function_descriptor{discrete_function_descriptor},
+      m_function_id{function_id}
+  {}
+
   DiscreteFunctionInterpoler(const DiscreteFunctionInterpoler&) = delete;
   DiscreteFunctionInterpoler(DiscreteFunctionInterpoler&&)      = delete;
 
diff --git a/src/scheme/DiscreteFunctionVectorIntegrator.cpp b/src/scheme/DiscreteFunctionVectorIntegrator.cpp
index 6d8937f84aeba6ae96ac5fc6e5ba523e5781f239..cfa65b1787090a06c6de0ac6bf774b3771a4a1eb 100644
--- a/src/scheme/DiscreteFunctionVectorIntegrator.cpp
+++ b/src/scheme/DiscreteFunctionVectorIntegrator.cpp
@@ -1,14 +1,76 @@
 #include <scheme/DiscreteFunctionVectorIntegrator.hpp>
 
 #include <language/utils/IntegrateCellArray.hpp>
-
+#include <mesh/MeshCellZone.hpp>
 #include <scheme/DiscreteFunctionP0Vector.hpp>
 #include <utils/Exceptions.hpp>
 
 template <size_t Dimension, typename DataType>
 std::shared_ptr<IDiscreteFunction>
-DiscreteFunctionVectorIntegrator::_integrate() const
+DiscreteFunctionVectorIntegrator::_integrateOnZoneList() const
+{
+  Assert(m_zone_list.size() > 0, "no zone list provided");
+
+  std::shared_ptr p_mesh = std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(m_mesh);
+
+  CellValue<bool> is_in_zone{p_mesh->connectivity()};
+  is_in_zone.fill(false);
+
+  size_t number_of_cells = 0;
+  for (const auto& zone : m_zone_list) {
+    auto mesh_cell_zone   = getMeshCellZone(*p_mesh, *zone);
+    const auto& cell_list = mesh_cell_zone.cellList();
+    for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+      const CellId cell_id = cell_list[i_cell];
+      if (is_in_zone[cell_id]) {
+        std::ostringstream os;
+        os << "cell " << cell_id << " (number " << p_mesh->connectivity().cellNumber()[cell_id]
+           << ") is present multiple times in zone list";
+        throw NormalError(os.str());
+      }
+      ++number_of_cells;
+      is_in_zone[cell_id] = true;
+    }
+  }
+
+  Array<CellId> zone_cell_list{number_of_cells};
+  {
+    size_t i_cell = 0;
+    for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+      if (is_in_zone[cell_id]) {
+        zone_cell_list[i_cell++] = cell_id;
+      }
+    }
+  }
+
+  Table<const DataType> table =
+    IntegrateCellArray<DataType(TinyVector<Dimension>)>::template integrate(m_function_id_list,
+                                                                            *m_quadrature_descriptor, *p_mesh,
+                                                                            zone_cell_list);
+
+  CellArray<DataType> cell_array{p_mesh->connectivity(), m_function_id_list.size()};
+  if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) {
+    cell_array.fill(zero);
+  } else {
+    cell_array.fill(0);
+  }
+
+  parallel_for(
+    zone_cell_list.size(), PUGS_LAMBDA(const size_t i_cell) {
+      for (size_t i = 0; i < table.numberOfColumns(); ++i) {
+        cell_array[zone_cell_list[i_cell]][i] = table[i_cell][i];
+      }
+    });
+
+  return std::make_shared<DiscreteFunctionP0Vector<Dimension, DataType>>(p_mesh, cell_array);
+}
+
+template <size_t Dimension, typename DataType>
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionVectorIntegrator::_integrateGlobally() const
 {
+  Assert(m_zone_list.size() == 0, "invalid call when zones are defined");
+
   std::shared_ptr mesh = std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(m_mesh);
 
   return std::make_shared<
@@ -17,6 +79,17 @@ DiscreteFunctionVectorIntegrator::_integrate() const
                                                                               *m_quadrature_descriptor, *mesh));
 }
 
+template <size_t Dimension, typename DataType>
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionVectorIntegrator::_integrate() const
+{
+  if (m_zone_list.size() == 0) {
+    return this->_integrateGlobally<Dimension, DataType>();
+  } else {
+    return this->_integrateOnZoneList<Dimension, DataType>();
+  }
+}
+
 template <size_t Dimension>
 std::shared_ptr<IDiscreteFunction>
 DiscreteFunctionVectorIntegrator::_integrate() const
@@ -36,7 +109,7 @@ DiscreteFunctionVectorIntegrator::_integrate() const
     default: {
       std::ostringstream os;
       os << "vector functions require scalar value type.\n"
-         << "Invalid interpolation value type: " << rang::fgB::red << dataTypeName(data_type) << rang::style::reset;
+         << "Invalid value type: " << rang::fgB::red << dataTypeName(data_type) << rang::style::reset;
       throw NormalError(os.str());
     }
     }
@@ -48,7 +121,7 @@ std::shared_ptr<IDiscreteFunction>
 DiscreteFunctionVectorIntegrator::integrate() const
 {
   if (m_discrete_function_descriptor->type() != DiscreteFunctionType::P0Vector) {
-    throw NormalError("invalid discrete function type for vector interpolation");
+    throw NormalError("invalid discrete function type for vector integration");
   }
 
   switch (m_mesh->dimension()) {
diff --git a/src/scheme/DiscreteFunctionVectorIntegrator.hpp b/src/scheme/DiscreteFunctionVectorIntegrator.hpp
index 01e086e7efd9c2334d07dfdb3cf8dce1f33e49f1..55f5fe3572cbb237e28b9bf86ff7bda19db20a98 100644
--- a/src/scheme/DiscreteFunctionVectorIntegrator.hpp
+++ b/src/scheme/DiscreteFunctionVectorIntegrator.hpp
@@ -4,6 +4,7 @@
 #include <analysis/IQuadratureDescriptor.hpp>
 #include <language/utils/FunctionSymbolId.hpp>
 #include <mesh/IMesh.hpp>
+#include <mesh/IZoneDescriptor.hpp>
 #include <scheme/IDiscreteFunction.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 
@@ -14,10 +15,17 @@ class DiscreteFunctionVectorIntegrator
 {
  private:
   std::shared_ptr<const IMesh> m_mesh;
+  const std::vector<std::shared_ptr<const IZoneDescriptor>> m_zone_list;
   std::shared_ptr<const IQuadratureDescriptor> m_quadrature_descriptor;
   std::shared_ptr<const IDiscreteFunctionDescriptor> m_discrete_function_descriptor;
   const std::vector<FunctionSymbolId> m_function_id_list;
 
+  template <size_t Dimension, typename DataType>
+  std::shared_ptr<IDiscreteFunction> _integrateOnZoneList() const;
+
+  template <size_t Dimension, typename DataType>
+  std::shared_ptr<IDiscreteFunction> _integrateGlobally() const;
+
   template <size_t Dimension, typename DataType>
   std::shared_ptr<IDiscreteFunction> _integrate() const;
 
@@ -38,6 +46,19 @@ class DiscreteFunctionVectorIntegrator
       m_function_id_list{function_id_list}
   {}
 
+  DiscreteFunctionVectorIntegrator(
+    const std::shared_ptr<const IMesh>& mesh,
+    const std::vector<std::shared_ptr<const IZoneDescriptor>>& zone_list,
+    const std::shared_ptr<const IQuadratureDescriptor>& quadrature_descriptor,
+    const std::shared_ptr<const IDiscreteFunctionDescriptor>& discrete_function_descriptor,
+    const std::vector<FunctionSymbolId>& function_id_list)
+    : m_mesh{mesh},
+      m_zone_list{zone_list},
+      m_quadrature_descriptor(quadrature_descriptor),
+      m_discrete_function_descriptor{discrete_function_descriptor},
+      m_function_id_list{function_id_list}
+  {}
+
   DiscreteFunctionVectorIntegrator(const DiscreteFunctionVectorIntegrator&) = delete;
   DiscreteFunctionVectorIntegrator(DiscreteFunctionVectorIntegrator&&)      = delete;
 
diff --git a/src/scheme/DiscreteFunctionVectorInterpoler.cpp b/src/scheme/DiscreteFunctionVectorInterpoler.cpp
index 98d6513c9a896b81d8e6909e057a8300b66d6dc7..128c735e4ae8e158f3dc42bf3d5ed47b17ac3fb5 100644
--- a/src/scheme/DiscreteFunctionVectorInterpoler.cpp
+++ b/src/scheme/DiscreteFunctionVectorInterpoler.cpp
@@ -1,22 +1,94 @@
 #include <scheme/DiscreteFunctionVectorInterpoler.hpp>
 
 #include <language/utils/InterpolateItemArray.hpp>
+#include <mesh/MeshCellZone.hpp>
 #include <scheme/DiscreteFunctionP0Vector.hpp>
 #include <utils/Exceptions.hpp>
 
 template <size_t Dimension, typename DataType>
 std::shared_ptr<IDiscreteFunction>
-DiscreteFunctionVectorInterpoler::_interpolate() const
+DiscreteFunctionVectorInterpoler::_interpolateOnZoneList() const
 {
-  std::shared_ptr mesh = std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(m_mesh);
+  Assert(m_zone_list.size() > 0, "no zone list provided");
 
+  std::shared_ptr p_mesh  = std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(m_mesh);
   using MeshDataType      = MeshData<Dimension>;
-  MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+  MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);
+
+  CellValue<bool> is_in_zone{p_mesh->connectivity()};
+  is_in_zone.fill(false);
+
+  size_t number_of_cells = 0;
+  for (const auto& zone : m_zone_list) {
+    auto mesh_cell_zone   = getMeshCellZone(*p_mesh, *zone);
+    const auto& cell_list = mesh_cell_zone.cellList();
+    for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+      const CellId cell_id = cell_list[i_cell];
+      if (is_in_zone[cell_id]) {
+        std::ostringstream os;
+        os << "cell " << cell_id << " (number " << p_mesh->connectivity().cellNumber()[cell_id]
+           << ") is present multiple times in zone list";
+        throw NormalError(os.str());
+      }
+      ++number_of_cells;
+      is_in_zone[cell_id] = true;
+    }
+  }
+
+  Array<CellId> zone_cell_list{number_of_cells};
+  {
+    size_t i_cell = 0;
+    for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+      if (is_in_zone[cell_id]) {
+        zone_cell_list[i_cell++] = cell_id;
+      }
+    }
+  }
+
+  Table<const DataType> table =
+    InterpolateItemArray<DataType(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(m_function_id_list,
+                                                                                                mesh_data.xj(),
+                                                                                                zone_cell_list);
+
+  CellArray<DataType> cell_array{p_mesh->connectivity(), m_function_id_list.size()};
+  cell_array.fill(0);
+
+  parallel_for(
+    zone_cell_list.size(), PUGS_LAMBDA(const size_t i_cell) {
+      for (size_t i = 0; i < table.numberOfColumns(); ++i) {
+        cell_array[zone_cell_list[i_cell]][i] = table(i_cell, i);
+      }
+    });
+
+  return std::make_shared<DiscreteFunctionP0Vector<Dimension, DataType>>(p_mesh, cell_array);
+}
+
+template <size_t Dimension, typename DataType>
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionVectorInterpoler::_interpolateGlobally() const
+{
+  Assert(m_zone_list.size() == 0, "invalid call when zones are defined");
+
+  std::shared_ptr p_mesh = std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(m_mesh);
+
+  using MeshDataType      = MeshData<Dimension>;
+  MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);
 
   return std::make_shared<
-    DiscreteFunctionP0Vector<Dimension, DataType>>(mesh, InterpolateItemArray<DataType(TinyVector<Dimension>)>::
-                                                           template interpolate<ItemType::cell>(m_function_id_list,
-                                                                                                mesh_data.xj()));
+    DiscreteFunctionP0Vector<Dimension, DataType>>(p_mesh, InterpolateItemArray<DataType(TinyVector<Dimension>)>::
+                                                             template interpolate<ItemType::cell>(m_function_id_list,
+                                                                                                  mesh_data.xj()));
+}
+
+template <size_t Dimension, typename DataType>
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionVectorInterpoler::_interpolate() const
+{
+  if (m_zone_list.size() == 0) {
+    return this->_interpolateGlobally<Dimension, DataType>();
+  } else {
+    return this->_interpolateOnZoneList<Dimension, DataType>();
+  }
 }
 
 template <size_t Dimension>
diff --git a/src/scheme/DiscreteFunctionVectorInterpoler.hpp b/src/scheme/DiscreteFunctionVectorInterpoler.hpp
index 093f290c40b99dc7c3c8f6db1ab32cbae10dfcb2..4bd917264fcfdb8992d6096db824bc17951c2603 100644
--- a/src/scheme/DiscreteFunctionVectorInterpoler.hpp
+++ b/src/scheme/DiscreteFunctionVectorInterpoler.hpp
@@ -3,6 +3,7 @@
 
 #include <language/utils/FunctionSymbolId.hpp>
 #include <mesh/IMesh.hpp>
+#include <mesh/IZoneDescriptor.hpp>
 #include <scheme/IDiscreteFunction.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 
@@ -13,9 +14,16 @@ class DiscreteFunctionVectorInterpoler
 {
  private:
   std::shared_ptr<const IMesh> m_mesh;
+  const std::vector<std::shared_ptr<const IZoneDescriptor>> m_zone_list;
   std::shared_ptr<const IDiscreteFunctionDescriptor> m_discrete_function_descriptor;
   const std::vector<FunctionSymbolId> m_function_id_list;
 
+  template <size_t Dimension, typename DataType>
+  std::shared_ptr<IDiscreteFunction> _interpolateOnZoneList() const;
+
+  template <size_t Dimension, typename DataType>
+  std::shared_ptr<IDiscreteFunction> _interpolateGlobally() const;
+
   template <size_t Dimension, typename DataType>
   std::shared_ptr<IDiscreteFunction> _interpolate() const;
 
@@ -32,6 +40,17 @@ class DiscreteFunctionVectorInterpoler
     : m_mesh{mesh}, m_discrete_function_descriptor{discrete_function_descriptor}, m_function_id_list{function_id_list}
   {}
 
+  DiscreteFunctionVectorInterpoler(
+    const std::shared_ptr<const IMesh>& mesh,
+    const std::vector<std::shared_ptr<const IZoneDescriptor>>& zone_list,
+    const std::shared_ptr<const IDiscreteFunctionDescriptor>& discrete_function_descriptor,
+    const std::vector<FunctionSymbolId>& function_id_list)
+    : m_mesh{mesh},
+      m_zone_list{zone_list},
+      m_discrete_function_descriptor{discrete_function_descriptor},
+      m_function_id_list{function_id_list}
+  {}
+
   DiscreteFunctionVectorInterpoler(const DiscreteFunctionVectorInterpoler&) = delete;
   DiscreteFunctionVectorInterpoler(DiscreteFunctionVectorInterpoler&&)      = delete;
 
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 37a87d841e2ad88946edcdfe46abb3fe8ca4ea71..52a41bcb69fc8181fce73b92059b49cf1cccf314 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -157,10 +157,16 @@ add_executable (unit_tests
 add_executable (mpi_unit_tests
   mpi_test_main.cpp
   test_Connectivity.cpp
+  test_DiscreteFunctionIntegrator.cpp
+  test_DiscreteFunctionIntegratorByZone.cpp
   test_DiscreteFunctionInterpoler.cpp
+  test_DiscreteFunctionInterpolerByZone.cpp
   test_DiscreteFunctionP0.cpp
   test_DiscreteFunctionP0Vector.cpp
+  test_DiscreteFunctionVectorIntegrator.cpp
+  test_DiscreteFunctionVectorIntegratorByZone.cpp
   test_DiscreteFunctionVectorInterpoler.cpp
+  test_DiscreteFunctionVectorInterpolerByZone.cpp
   test_EmbeddedIDiscreteFunctionMathFunctions.cpp
   test_EmbeddedIDiscreteFunctionOperators.cpp
   test_InterpolateItemArray.cpp
diff --git a/tests/test_DiscreteFunctionIntegrator.cpp b/tests/test_DiscreteFunctionIntegrator.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..465288f67d4efcd9a7a78eff53ab9649d381e572
--- /dev/null
+++ b/tests/test_DiscreteFunctionIntegrator.cpp
@@ -0,0 +1,781 @@
+#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/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <language/utils/IntegrateCellValue.hpp>
+#include <scheme/DiscreteFunctionDescriptorP0.hpp>
+#include <scheme/DiscreteFunctionIntegrator.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+
+#include <pegtl/string_input.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("DiscreteFunctionIntegrator", "[scheme]")
+{
+  auto same_cell_value = [](auto f, auto g) -> bool {
+    using ItemIdType = typename decltype(f)::index_type;
+    for (ItemIdType item_id = 0; item_id < f.numberOfItems(); ++item_id) {
+      if (f[item_id] != g[item_id]) {
+        return false;
+      }
+    }
+
+    return true;
+  };
+
+  SECTION("1D")
+  {
+    constexpr size_t Dimension = 1;
+
+    std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+    std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+    for (auto named_mesh : mesh_list) {
+      SECTION(named_mesh.name())
+      {
+        auto mesh_1d = named_mesh.mesh();
+
+        std::string_view data = R"(
+import math;
+let B_scalar_non_linear_1d: R^1 -> B, x -> (exp(2 * x[0]) + 3 > 4);
+let N_scalar_non_linear_1d: R^1 -> N, x -> floor(3 * x[0] * x[0] + 2);
+let Z_scalar_non_linear_1d: R^1 -> Z, x -> floor(exp(2 * x[0]) - 1);
+let R_scalar_non_linear_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+let R1_non_linear_1d: R^1 -> R^1, x -> 2 * exp(x[0]);
+let R2_non_linear_1d: R^1 -> R^2, x -> (2 * exp(x[0]), -3*x[0]);
+let R3_non_linear_1d: R^1 -> R^3, x -> (2 * exp(x[0]) + 3, x[0] - 2, 3);
+let R1x1_non_linear_1d: R^1 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[0]) + 3);
+let R2x2_non_linear_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]);
+let R3x3_non_linear_1d: R^1 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0], -4*x[0], 2*x[0]+1, 3, -6*x[0], exp(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("B_scalar_non_linear_1d")
+        {
+          auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_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_value =
+            IntegrateCellValue<double(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d);
+
+          DiscreteFunctionIntegrator integrator(mesh_1d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("N_scalar_non_linear_1d")
+        {
+          auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_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_value =
+            IntegrateCellValue<double(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d);
+
+          DiscreteFunctionIntegrator integrator(mesh_1d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("Z_scalar_non_linear_1d")
+        {
+          auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_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_value =
+            IntegrateCellValue<double(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d);
+
+          DiscreteFunctionIntegrator integrator(mesh_1d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("R_scalar_non_linear_1d")
+        {
+          auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_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_value =
+            IntegrateCellValue<double(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d);
+
+          DiscreteFunctionIntegrator integrator(mesh_1d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("R1_non_linear_1d")
+        {
+          using DataType = TinyVector<1>;
+
+          auto [i_symbol, found] = symbol_table->find("R1_non_linear_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<DataType> cell_value =
+            IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor,
+                                                                   *mesh_1d);
+
+          DiscreteFunctionIntegrator integrator(mesh_1d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R2_non_linear_1d")
+        {
+          using DataType = TinyVector<2>;
+
+          auto [i_symbol, found] = symbol_table->find("R2_non_linear_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<DataType> cell_value =
+            IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor,
+                                                                   *mesh_1d);
+
+          DiscreteFunctionIntegrator integrator(mesh_1d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R3_non_linear_1d")
+        {
+          using DataType = TinyVector<3>;
+
+          auto [i_symbol, found] = symbol_table->find("R3_non_linear_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<DataType> cell_value =
+            IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor,
+                                                                   *mesh_1d);
+
+          DiscreteFunctionIntegrator integrator(mesh_1d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R1x1_non_linear_1d")
+        {
+          using DataType = TinyMatrix<1>;
+
+          auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_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<DataType> cell_value =
+            IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor,
+                                                                   *mesh_1d);
+
+          DiscreteFunctionIntegrator integrator(mesh_1d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R2x2_non_linear_1d")
+        {
+          using DataType = TinyMatrix<2>;
+
+          auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_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<DataType> cell_value =
+            IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor,
+                                                                   *mesh_1d);
+
+          DiscreteFunctionIntegrator integrator(mesh_1d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R3x3_non_linear_1d")
+        {
+          using DataType = TinyMatrix<3>;
+
+          auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_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<DataType> cell_value =
+            IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor,
+                                                                   *mesh_1d);
+
+          DiscreteFunctionIntegrator integrator(mesh_1d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+      }
+    }
+  }
+
+  SECTION("2D")
+  {
+    constexpr size_t Dimension = 2;
+
+    std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+    std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+    for (auto named_mesh : mesh_list) {
+      SECTION(named_mesh.name())
+      {
+        auto mesh_2d = named_mesh.mesh();
+
+        std::string_view data = R"(
+import math;
+let B_scalar_non_linear_2d: R^2 -> B, x -> (exp(2 * x[0])< 2*x[1]);
+let N_scalar_non_linear_2d: R^2 -> N, x -> floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + 2);
+let Z_scalar_non_linear_2d: R^2 -> Z, x -> floor(exp(2 * x[0]) - 3 * x[1]);
+let R_scalar_non_linear_2d: R^2 -> R, x -> 2 * exp(x[0]) + 3 * x[1];
+let R1_non_linear_2d: R^2 -> R^1, x -> 2 * exp(x[0]);
+let R2_non_linear_2d: R^2 -> R^2, x -> (2 * exp(x[0]), -3*x[1]);
+let R3_non_linear_2d: R^2 -> R^3, x -> (2 * exp(x[0]) + 3, x[1] - 2, 3);
+let R1x1_non_linear_2d: R^2 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[1]) + 3);
+let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[0]), 3, x[1] * x[0]);
+let R3x3_non_linear_2d: R^2 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[0]), 3, x[1] * x[0], -4*x[1], 2*x[0]+1, 3, -6*x[0], exp(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("B_scalar_non_linear_2d")
+        {
+          auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_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_value =
+            IntegrateCellValue<double(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d);
+
+          DiscreteFunctionIntegrator integrator(mesh_2d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("N_scalar_non_linear_2d")
+        {
+          auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_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_value =
+            IntegrateCellValue<double(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d);
+
+          DiscreteFunctionIntegrator integrator(mesh_2d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("Z_scalar_non_linear_2d")
+        {
+          auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_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_value =
+            IntegrateCellValue<double(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d);
+
+          DiscreteFunctionIntegrator integrator(mesh_2d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("R_scalar_non_linear_2d")
+        {
+          auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_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_value =
+            IntegrateCellValue<double(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d);
+
+          DiscreteFunctionIntegrator integrator(mesh_2d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("R1_non_linear_2d")
+        {
+          using DataType = TinyVector<1>;
+
+          auto [i_symbol, found] = symbol_table->find("R1_non_linear_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<DataType> cell_value =
+            IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor,
+                                                                   *mesh_2d);
+
+          DiscreteFunctionIntegrator integrator(mesh_2d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R2_non_linear_2d")
+        {
+          using DataType = TinyVector<2>;
+
+          auto [i_symbol, found] = symbol_table->find("R2_non_linear_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<DataType> cell_value =
+            IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor,
+                                                                   *mesh_2d);
+
+          DiscreteFunctionIntegrator integrator(mesh_2d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R3_non_linear_2d")
+        {
+          using DataType = TinyVector<3>;
+
+          auto [i_symbol, found] = symbol_table->find("R3_non_linear_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<DataType> cell_value =
+            IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor,
+                                                                   *mesh_2d);
+
+          DiscreteFunctionIntegrator integrator(mesh_2d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R1x1_non_linear_2d")
+        {
+          using DataType = TinyMatrix<1>;
+
+          auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_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<DataType> cell_value =
+            IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor,
+                                                                   *mesh_2d);
+
+          DiscreteFunctionIntegrator integrator(mesh_2d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R2x2_non_linear_2d")
+        {
+          using DataType = TinyMatrix<2>;
+
+          auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_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<DataType> cell_value =
+            IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor,
+                                                                   *mesh_2d);
+
+          DiscreteFunctionIntegrator integrator(mesh_2d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R3x3_non_linear_2d")
+        {
+          using DataType = TinyMatrix<3>;
+
+          auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_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<DataType> cell_value =
+            IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor,
+                                                                   *mesh_2d);
+
+          DiscreteFunctionIntegrator integrator(mesh_2d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+      }
+    }
+  }
+
+  SECTION("3D")
+  {
+    constexpr size_t Dimension = 3;
+
+    std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+    std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+    for (auto named_mesh : mesh_list) {
+      SECTION(named_mesh.name())
+      {
+        auto mesh_3d = named_mesh.mesh();
+
+        std::string_view data = R"(
+import math;
+let B_scalar_non_linear_3d: R^3 -> B, x -> (exp(2 * x[0])< 2*x[1]+x[2]);
+let N_scalar_non_linear_3d: R^3 -> N, x -> floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + x[2] * x[2]);
+let Z_scalar_non_linear_3d: R^3 -> Z, x -> floor(exp(2 * x[0]) - 3 * x[1] + x[2]);
+let R_scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0]+x[2]) + 3 * x[1];
+let R1_non_linear_3d: R^3 -> R^1, x -> 2 * exp(x[0])+sin(x[1] + x[2]);
+let R2_non_linear_3d: R^3 -> R^2, x -> (2 * exp(x[0]), -3*x[1] * x[2]);
+let R3_non_linear_3d: R^3 -> R^3, x -> (2 * exp(x[0]) + 3, x[1] - 2, 3 * x[2]);
+let R1x1_non_linear_3d: R^3 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * x[2]);
+let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[2] - 2 * x[0]), 3, x[1] * x[0] - x[2]);
+let R3x3_non_linear_3d: R^3 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[2]), 3, x[1] * x[2], -4*x[1], 2*x[2]+1, 3, -6*x[2], exp(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("B_scalar_non_linear_3d")
+        {
+          auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_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_value =
+            IntegrateCellValue<double(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d);
+
+          DiscreteFunctionIntegrator integrator(mesh_3d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("N_scalar_non_linear_3d")
+        {
+          auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_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_value =
+            IntegrateCellValue<double(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d);
+
+          DiscreteFunctionIntegrator integrator(mesh_3d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("Z_scalar_non_linear_3d")
+        {
+          auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_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_value =
+            IntegrateCellValue<double(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d);
+
+          DiscreteFunctionIntegrator integrator(mesh_3d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("R_scalar_non_linear_3d")
+        {
+          auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_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_value =
+            IntegrateCellValue<double(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d);
+
+          DiscreteFunctionIntegrator integrator(mesh_3d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("R1_non_linear_3d")
+        {
+          using DataType = TinyVector<1>;
+
+          auto [i_symbol, found] = symbol_table->find("R1_non_linear_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<DataType> cell_value =
+            IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor,
+                                                                   *mesh_3d);
+
+          DiscreteFunctionIntegrator integrator(mesh_3d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R2_non_linear_3d")
+        {
+          using DataType = TinyVector<2>;
+
+          auto [i_symbol, found] = symbol_table->find("R2_non_linear_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<DataType> cell_value =
+            IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor,
+                                                                   *mesh_3d);
+
+          DiscreteFunctionIntegrator integrator(mesh_3d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R3_non_linear_3d")
+        {
+          using DataType = TinyVector<3>;
+
+          auto [i_symbol, found] = symbol_table->find("R3_non_linear_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<DataType> cell_value =
+            IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor,
+                                                                   *mesh_3d);
+
+          DiscreteFunctionIntegrator integrator(mesh_3d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R1x1_non_linear_3d")
+        {
+          using DataType = TinyMatrix<1>;
+
+          auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_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<DataType> cell_value =
+            IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor,
+                                                                   *mesh_3d);
+
+          DiscreteFunctionIntegrator integrator(mesh_3d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R2x2_non_linear_3d")
+        {
+          using DataType = TinyMatrix<2>;
+
+          auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_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<DataType> cell_value =
+            IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor,
+                                                                   *mesh_3d);
+
+          DiscreteFunctionIntegrator integrator(mesh_3d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R3x3_non_linear_3d")
+        {
+          using DataType = TinyMatrix<3>;
+
+          auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_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<DataType> cell_value =
+            IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor,
+                                                                   *mesh_3d);
+
+          DiscreteFunctionIntegrator integrator(mesh_3d, quadrature_descriptor, function_symbol_id);
+          std::shared_ptr discrete_function = integrator.integrate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+      }
+    }
+  }
+}
diff --git a/tests/test_DiscreteFunctionIntegratorByZone.cpp b/tests/test_DiscreteFunctionIntegratorByZone.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..071f8cddb7e9f71e0b102bf361afb8715a7afe87
--- /dev/null
+++ b/tests/test_DiscreteFunctionIntegratorByZone.cpp
@@ -0,0 +1,1062 @@
+#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/Mesh.hpp>
+#include <mesh/MeshCellZone.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/NamedZoneDescriptor.hpp>
+
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <language/utils/IntegrateCellValue.hpp>
+#include <scheme/DiscreteFunctionDescriptorP0.hpp>
+#include <scheme/DiscreteFunctionIntegrator.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+
+#include <pegtl/string_input.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("DiscreteFunctionIntegratorByZone", "[scheme]")
+{
+  auto same_cell_value = [](auto f, auto g) -> bool {
+    using ItemIdType = typename decltype(f)::index_type;
+    for (ItemIdType item_id = 0; item_id < f.numberOfItems(); ++item_id) {
+      if (f[item_id] != g[item_id]) {
+        return false;
+      }
+    }
+
+    return true;
+  };
+
+  SECTION("1D")
+  {
+    constexpr size_t Dimension = 1;
+
+    std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+    auto mesh_1d = MeshDataBaseForTests::get().unordered1DMesh();
+
+    std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+    zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+    auto mesh_cell_zone = getMeshCellZone(*mesh_1d, *zone_list[0]);
+    auto zone_cell_list = mesh_cell_zone.cellList();
+
+    std::string_view data = R"(
+import math;
+let B_scalar_non_linear_1d: R^1 -> B, x -> (exp(2 * x[0]) + 3 < 3.3);
+let N_scalar_non_linear_1d: R^1 -> N, x -> floor(3 * x[0] * x[0] + 2);
+let Z_scalar_non_linear_1d: R^1 -> Z, x -> floor(exp(2 * x[0]) - 1);
+let R_scalar_non_linear_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+let R1_non_linear_1d: R^1 -> R^1, x -> 2 * exp(x[0]);
+let R2_non_linear_1d: R^1 -> R^2, x -> (2 * exp(x[0]), -3*x[0]);
+let R3_non_linear_1d: R^1 -> R^3, x -> (2 * exp(x[0]) + 3, x[0] - 2, 3);
+let R1x1_non_linear_1d: R^1 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[0]) + 3);
+let R2x2_non_linear_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]);
+let R3x3_non_linear_1d: R^1 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0], -4*x[0], 2*x[0]+1, 3, -6*x[0], exp(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("B_scalar_non_linear_1d")
+    {
+      auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_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);
+
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_1d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("N_scalar_non_linear_1d")
+    {
+      auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_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);
+
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_1d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("Z_scalar_non_linear_1d")
+    {
+      auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_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);
+
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_1d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("R_scalar_non_linear_1d")
+    {
+      auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_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);
+
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_1d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("R1_non_linear_1d")
+    {
+      using DataType = TinyVector<1>;
+
+      auto [i_symbol, found] = symbol_table->find("R1_non_linear_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);
+
+      Array<DataType> array =
+        IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d,
+                                                               zone_cell_list);
+
+      CellValue<DataType> cell_value{mesh_1d->connectivity()};
+      cell_value.fill(zero);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R2_non_linear_1d")
+    {
+      using DataType = TinyVector<2>;
+
+      auto [i_symbol, found] = symbol_table->find("R2_non_linear_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);
+
+      Array<DataType> array =
+        IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d,
+                                                               zone_cell_list);
+
+      CellValue<DataType> cell_value{mesh_1d->connectivity()};
+      cell_value.fill(zero);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R3_non_linear_1d")
+    {
+      using DataType = TinyVector<3>;
+
+      auto [i_symbol, found] = symbol_table->find("R3_non_linear_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);
+
+      Array<DataType> array =
+        IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d,
+                                                               zone_cell_list);
+
+      CellValue<DataType> cell_value{mesh_1d->connectivity()};
+      cell_value.fill(zero);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R1x1_non_linear_1d")
+    {
+      using DataType = TinyMatrix<1>;
+
+      auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_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);
+
+      Array<DataType> array =
+        IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d,
+                                                               zone_cell_list);
+
+      CellValue<DataType> cell_value{mesh_1d->connectivity()};
+      cell_value.fill(zero);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R2x2_non_linear_1d")
+    {
+      using DataType = TinyMatrix<2>;
+
+      auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_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);
+
+      Array<DataType> array =
+        IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d,
+                                                               zone_cell_list);
+
+      CellValue<DataType> cell_value{mesh_1d->connectivity()};
+      cell_value.fill(zero);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R3x3_non_linear_1d")
+    {
+      using DataType = TinyMatrix<3>;
+
+      auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_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);
+
+      Array<DataType> array =
+        IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d,
+                                                               zone_cell_list);
+
+      CellValue<DataType> cell_value{mesh_1d->connectivity()};
+      cell_value.fill(zero);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+  }
+
+  SECTION("2D")
+  {
+    constexpr size_t Dimension = 2;
+
+    std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+    auto mesh_2d = MeshDataBaseForTests::get().hybrid2DMesh();
+
+    std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+    zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+    auto mesh_cell_zone = getMeshCellZone(*mesh_2d, *zone_list[0]);
+    auto zone_cell_list = mesh_cell_zone.cellList();
+
+    std::string_view data = R"(
+import math;
+let B_scalar_non_linear_2d: R^2 -> B, x -> (exp(2 * x[0])< 2*x[1]);
+let N_scalar_non_linear_2d: R^2 -> N, x -> floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + 2);
+let Z_scalar_non_linear_2d: R^2 -> Z, x -> floor(exp(2 * x[0]) - 3 * x[1]);
+let R_scalar_non_linear_2d: R^2 -> R, x -> 2 * exp(x[0]) + 3 * x[1];
+let R1_non_linear_2d: R^2 -> R^1, x -> 2 * exp(x[0]);
+let R2_non_linear_2d: R^2 -> R^2, x -> (2 * exp(x[0]), -3*x[1]);
+let R3_non_linear_2d: R^2 -> R^3, x -> (2 * exp(x[0]) + 3, x[1] - 2, 3);
+let R1x1_non_linear_2d: R^2 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[1]) + 3);
+let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[0]), 3, x[1] * x[0]);
+let R3x3_non_linear_2d: R^2 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[0]), 3, x[1] * x[0], -4*x[1], 2*x[0]+1, 3, -6*x[0], exp(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("B_scalar_non_linear_2d")
+    {
+      auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_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);
+
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_2d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("N_scalar_non_linear_2d")
+    {
+      auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_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);
+
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_2d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("Z_scalar_non_linear_2d")
+    {
+      auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_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);
+
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_2d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("R_scalar_non_linear_2d")
+    {
+      auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_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);
+
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_2d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("R1_non_linear_2d")
+    {
+      using DataType = TinyVector<1>;
+
+      auto [i_symbol, found] = symbol_table->find("R1_non_linear_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);
+
+      Array<DataType> array =
+        IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d,
+                                                               zone_cell_list);
+
+      CellValue<DataType> cell_value{mesh_2d->connectivity()};
+      cell_value.fill(zero);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R2_non_linear_2d")
+    {
+      using DataType = TinyVector<2>;
+
+      auto [i_symbol, found] = symbol_table->find("R2_non_linear_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);
+
+      Array<DataType> array =
+        IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d,
+                                                               zone_cell_list);
+
+      CellValue<DataType> cell_value{mesh_2d->connectivity()};
+      cell_value.fill(zero);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R3_non_linear_2d")
+    {
+      using DataType = TinyVector<3>;
+
+      auto [i_symbol, found] = symbol_table->find("R3_non_linear_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);
+
+      Array<DataType> array =
+        IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d,
+                                                               zone_cell_list);
+
+      CellValue<DataType> cell_value{mesh_2d->connectivity()};
+      cell_value.fill(zero);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R1x1_non_linear_2d")
+    {
+      using DataType = TinyMatrix<1>;
+
+      auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_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);
+
+      Array<DataType> array =
+        IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d,
+                                                               zone_cell_list);
+
+      CellValue<DataType> cell_value{mesh_2d->connectivity()};
+      cell_value.fill(zero);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R2x2_non_linear_2d")
+    {
+      using DataType = TinyMatrix<2>;
+
+      auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_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);
+
+      Array<DataType> array =
+        IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d,
+                                                               zone_cell_list);
+
+      CellValue<DataType> cell_value{mesh_2d->connectivity()};
+      cell_value.fill(zero);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R3x3_non_linear_2d")
+    {
+      using DataType = TinyMatrix<3>;
+
+      auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_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);
+
+      Array<DataType> array =
+        IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d,
+                                                               zone_cell_list);
+
+      CellValue<DataType> cell_value{mesh_2d->connectivity()};
+      cell_value.fill(zero);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+  }
+
+  SECTION("3D")
+  {
+    constexpr size_t Dimension = 3;
+
+    std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+    auto mesh_3d = MeshDataBaseForTests::get().hybrid3DMesh();
+
+    std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+    zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+    auto mesh_cell_zone = getMeshCellZone(*mesh_3d, *zone_list[0]);
+    auto zone_cell_list = mesh_cell_zone.cellList();
+
+    std::string_view data = R"(
+import math;
+let B_scalar_non_linear_3d: R^3 -> B, x -> (exp(2 * x[0])< 2*x[1]+x[2]);
+let N_scalar_non_linear_3d: R^3 -> N, x -> floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + x[2] * x[2]);
+let Z_scalar_non_linear_3d: R^3 -> Z, x -> floor(exp(2 * x[0]) - 3 * x[1] + x[2]);
+let R_scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0]+x[2]) + 3 * x[1];
+let R1_non_linear_3d: R^3 -> R^1, x -> 2 * exp(x[0])+sin(x[1] + x[2]);
+let R2_non_linear_3d: R^3 -> R^2, x -> (2 * exp(x[0]), -3*x[1] * x[2]);
+let R3_non_linear_3d: R^3 -> R^3, x -> (2 * exp(x[0]) + 3, x[1] - 2, 3 * x[2]);
+let R1x1_non_linear_3d: R^3 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * x[2]);
+let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[2] - 2 * x[0]), 3, x[1] * x[0] - x[2]);
+let R3x3_non_linear_3d: R^3 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[2]), 3, x[1] * x[2], -4*x[1], 2*x[2]+1, 3, -6*x[2], exp(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("B_scalar_non_linear_3d")
+    {
+      auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_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);
+
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_3d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("N_scalar_non_linear_3d")
+    {
+      auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_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);
+
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_3d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("Z_scalar_non_linear_3d")
+    {
+      auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_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);
+
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_3d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("R_scalar_non_linear_3d")
+    {
+      auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_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);
+
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_3d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("R1_non_linear_3d")
+    {
+      using DataType = TinyVector<1>;
+
+      auto [i_symbol, found] = symbol_table->find("R1_non_linear_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);
+
+      Array<DataType> array =
+        IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d,
+                                                               zone_cell_list);
+
+      CellValue<DataType> cell_value{mesh_3d->connectivity()};
+      cell_value.fill(zero);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R2_non_linear_3d")
+    {
+      using DataType = TinyVector<2>;
+
+      auto [i_symbol, found] = symbol_table->find("R2_non_linear_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);
+
+      Array<DataType> array =
+        IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d,
+                                                               zone_cell_list);
+
+      CellValue<DataType> cell_value{mesh_3d->connectivity()};
+      cell_value.fill(zero);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R3_non_linear_3d")
+    {
+      using DataType = TinyVector<3>;
+
+      auto [i_symbol, found] = symbol_table->find("R3_non_linear_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);
+
+      Array<DataType> array =
+        IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d,
+                                                               zone_cell_list);
+
+      CellValue<DataType> cell_value{mesh_3d->connectivity()};
+      cell_value.fill(zero);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R1x1_non_linear_3d")
+    {
+      using DataType = TinyMatrix<1>;
+
+      auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_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);
+
+      Array<DataType> array =
+        IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d,
+                                                               zone_cell_list);
+
+      CellValue<DataType> cell_value{mesh_3d->connectivity()};
+      cell_value.fill(zero);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R2x2_non_linear_3d")
+    {
+      using DataType = TinyMatrix<2>;
+
+      auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_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);
+
+      Array<DataType> array =
+        IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d,
+                                                               zone_cell_list);
+
+      CellValue<DataType> cell_value{mesh_3d->connectivity()};
+      cell_value.fill(zero);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R3x3_non_linear_3d")
+    {
+      using DataType = TinyMatrix<3>;
+
+      auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_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);
+
+      Array<DataType> array =
+        IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d,
+                                                               zone_cell_list);
+
+      CellValue<DataType> cell_value{mesh_3d->connectivity()};
+      cell_value.fill(zero);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t i) {
+          const CellId cell_id = zone_cell_list[i];
+          cell_value[cell_id]  = array[i];
+        });
+
+      DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id);
+      std::shared_ptr discrete_function = integrator.integrate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+  }
+}
diff --git a/tests/test_DiscreteFunctionInterpolerByZone.cpp b/tests/test_DiscreteFunctionInterpolerByZone.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4ad5f464f5fe276bf489df459cd2cfc90e42a519
--- /dev/null
+++ b/tests/test_DiscreteFunctionInterpolerByZone.cpp
@@ -0,0 +1,1074 @@
+#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/Mesh.hpp>
+#include <mesh/MeshCellZone.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/NamedZoneDescriptor.hpp>
+
+#include <scheme/DiscreteFunctionDescriptorP0.hpp>
+#include <scheme/DiscreteFunctionInterpoler.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+
+#include <pegtl/string_input.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("DiscreteFunctionInterpolerByZone", "[scheme]")
+{
+  auto same_cell_value = [](auto f, auto g) -> bool {
+    using ItemIdType = typename decltype(f)::index_type;
+    for (ItemIdType item_id = 0; item_id < f.numberOfItems(); ++item_id) {
+      if (f[item_id] != g[item_id]) {
+        return false;
+      }
+    }
+
+    return true;
+  };
+
+  SECTION("1D")
+  {
+    constexpr size_t Dimension = 1;
+
+    auto mesh_1d = MeshDataBaseForTests::get().unordered1DMesh();
+
+    std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+    zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+    auto mesh_cell_zone = getMeshCellZone(*mesh_1d, *zone_list[0]);
+    CellValue<bool> is_cell_in_zone{mesh_1d->connectivity()};
+    is_cell_in_zone.fill(false);
+    auto zone_cell_list = mesh_cell_zone.cellList();
+    for (size_t i_cell = 0; i_cell < zone_cell_list.size(); ++i_cell) {
+      is_cell_in_zone[zone_cell_list[i_cell]] = true;
+    }
+
+    auto xj = MeshDataManager::instance().getMeshData(*mesh_1d).xj();
+
+    std::string_view data = R"(
+import math;
+let B_scalar_non_linear_1d: R^1 -> B, x -> (exp(2 * x[0]) + 3 < 3.3);
+let N_scalar_non_linear_1d: R^1 -> N, x -> floor(3 * x[0] * x[0] + 2);
+let Z_scalar_non_linear_1d: R^1 -> Z, x -> floor(exp(2 * x[0]) - 1);
+let R_scalar_non_linear_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+let R1_non_linear_1d: R^1 -> R^1, x -> 2 * exp(x[0]);
+let R2_non_linear_1d: R^1 -> R^2, x -> (2 * exp(x[0]), -3*x[0]);
+let R3_non_linear_1d: R^1 -> R^3, x -> (2 * exp(x[0]) + 3, x[0] - 2, 3);
+let R1x1_non_linear_1d: R^1 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[0]) + 3);
+let R2x2_non_linear_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]);
+let R3x3_non_linear_1d: R^1 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0], -4*x[0], 2*x[0]+1, 3, -6*x[0], exp(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("B_scalar_non_linear_1d")
+    {
+      auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_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_value{mesh_1d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::exp(2 * x[0]) + 3 < 3.3;
+          } else {
+            cell_value[cell_id] = false;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("N_scalar_non_linear_1d")
+    {
+      auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_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_value{mesh_1d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::floor(3 * x[0] * x[0] + 2);
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("Z_scalar_non_linear_1d")
+    {
+      auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_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_value{mesh_1d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::floor(std::exp(2 * x[0]) - 1);
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("R_scalar_non_linear_1d")
+    {
+      auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_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_value{mesh_1d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = 2 * std::exp(x[0]) + 3;
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("R1_non_linear_1d")
+    {
+      using DataType = TinyVector<1>;
+
+      auto [i_symbol, found] = symbol_table->find("R1_non_linear_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<DataType> cell_value{mesh_1d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0])};
+          } else {
+            cell_value[cell_id] = zero;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R2_non_linear_1d")
+    {
+      using DataType = TinyVector<2>;
+
+      auto [i_symbol, found] = symbol_table->find("R2_non_linear_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<DataType> cell_value{mesh_1d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0]), -3 * x[0]};
+          } else {
+            cell_value[cell_id] = zero;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R3_non_linear_1d")
+    {
+      using DataType = TinyVector<3>;
+
+      auto [i_symbol, found] = symbol_table->find("R3_non_linear_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<DataType> cell_value{mesh_1d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + 3, x[0] - 2, 3};
+          } else {
+            cell_value[cell_id] = zero;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R1x1_non_linear_1d")
+    {
+      using DataType = TinyMatrix<1>;
+
+      auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_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<DataType> cell_value{mesh_1d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0]) * std::sin(x[0]) + 3};
+          } else {
+            cell_value[cell_id] = zero;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R2x2_non_linear_1d")
+    {
+      using DataType = TinyMatrix<2>;
+
+      auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_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<DataType> cell_value{mesh_1d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id] =
+              DataType{2 * std::exp(x[0]) * std::sin(x[0]) + 3, std::sin(x[0] - 2 * x[0]), 3, x[0] * x[0]};
+          } else {
+            cell_value[cell_id] = zero;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R3x3_non_linear_1d")
+    {
+      using DataType = TinyMatrix<3>;
+
+      auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_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<DataType> cell_value{mesh_1d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * exp(x[0]) * std::sin(x[0]) + 3,
+                                           std::sin(x[0] - 2 * x[0]),
+                                           3,
+                                           x[0] * x[0],
+                                           -4 * x[0],
+                                           2 * x[0] + 1,
+                                           3,
+                                           -6 * x[0],
+                                           std::exp(x[0])};
+          } else {
+            cell_value[cell_id] = zero;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+  }
+
+  SECTION("2D")
+  {
+    constexpr size_t Dimension = 2;
+
+    auto mesh_2d = MeshDataBaseForTests::get().hybrid2DMesh();
+
+    std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+    zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+    auto mesh_cell_zone = getMeshCellZone(*mesh_2d, *zone_list[0]);
+    CellValue<bool> is_cell_in_zone{mesh_2d->connectivity()};
+    is_cell_in_zone.fill(false);
+    auto zone_cell_list = mesh_cell_zone.cellList();
+    for (size_t i_cell = 0; i_cell < zone_cell_list.size(); ++i_cell) {
+      is_cell_in_zone[zone_cell_list[i_cell]] = true;
+    }
+
+    auto xj = MeshDataManager::instance().getMeshData(*mesh_2d).xj();
+
+    std::string_view data = R"(
+import math;
+let B_scalar_non_linear_2d: R^2 -> B, x -> (exp(2 * x[0])< 2*x[1]);
+let N_scalar_non_linear_2d: R^2 -> N, x -> floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + 2);
+let Z_scalar_non_linear_2d: R^2 -> Z, x -> floor(exp(2 * x[0]) - 3 * x[1]);
+let R_scalar_non_linear_2d: R^2 -> R, x -> 2 * exp(x[0]) + 3 * x[1];
+let R1_non_linear_2d: R^2 -> R^1, x -> 2 * exp(x[0]);
+let R2_non_linear_2d: R^2 -> R^2, x -> (2 * exp(x[0]), -3*x[1]);
+let R3_non_linear_2d: R^2 -> R^3, x -> (2 * exp(x[0]) + 3, x[1] - 2, 3);
+let R1x1_non_linear_2d: R^2 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[1]) + 3);
+let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[0]), 3, x[1] * x[0]);
+let R3x3_non_linear_2d: R^2 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[0]), 3, x[1] * x[0], -4*x[1], 2*x[0]+1, 3, -6*x[0], exp(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("B_scalar_non_linear_2d")
+    {
+      auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_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_value{mesh_2d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::exp(2 * x[0]) < 2 * x[1];
+          } else {
+            cell_value[cell_id] = false;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("N_scalar_non_linear_2d")
+    {
+      auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_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_value{mesh_2d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + 2);
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("Z_scalar_non_linear_2d")
+    {
+      auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_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_value{mesh_2d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::floor(std::exp(2 * x[0]) - 3 * x[1]);
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("R_scalar_non_linear_2d")
+    {
+      auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_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_value{mesh_2d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = 2 * std::exp(x[0]) + 3 * x[1];
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("R1_non_linear_2d")
+    {
+      using DataType = TinyVector<1>;
+
+      auto [i_symbol, found] = symbol_table->find("R1_non_linear_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<DataType> cell_value{mesh_2d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0])};
+          } else {
+            cell_value[cell_id] = zero;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R2_non_linear_2d")
+    {
+      using DataType = TinyVector<2>;
+
+      auto [i_symbol, found] = symbol_table->find("R2_non_linear_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<DataType> cell_value{mesh_2d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0]), -3 * x[1]};
+          } else {
+            cell_value[cell_id] = zero;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R3_non_linear_2d")
+    {
+      using DataType = TinyVector<3>;
+
+      auto [i_symbol, found] = symbol_table->find("R3_non_linear_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<DataType> cell_value{mesh_2d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + 3, x[1] - 2, 3};
+          } else {
+            cell_value[cell_id] = zero;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R1x1_non_linear_2d")
+    {
+      using DataType = TinyMatrix<1>;
+
+      auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_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<DataType> cell_value{mesh_2d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3};
+          } else {
+            cell_value[cell_id] = zero;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R2x2_non_linear_2d")
+    {
+      using DataType = TinyMatrix<2>;
+
+      auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_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<DataType> cell_value{mesh_2d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id] =
+              DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3, std::sin(x[1] - 2 * x[0]), 3, x[1] * x[0]};
+          } else {
+            cell_value[cell_id] = zero;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R3x3_non_linear_2d")
+    {
+      using DataType = TinyMatrix<3>;
+
+      auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_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<DataType> cell_value{mesh_2d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+
+            cell_value[cell_id] = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3,
+                                           std::sin(x[1] - 2 * x[0]),
+                                           3,
+                                           x[1] * x[0],
+                                           -4 * x[1],
+                                           2 * x[0] + 1,
+                                           3,
+                                           -6 * x[0],
+                                           std::exp(x[1])};
+          } else {
+            cell_value[cell_id] = zero;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+  }
+
+  SECTION("3D")
+  {
+    constexpr size_t Dimension = 3;
+
+    auto mesh_3d = MeshDataBaseForTests::get().hybrid3DMesh();
+
+    std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+    zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+    auto mesh_cell_zone = getMeshCellZone(*mesh_3d, *zone_list[0]);
+    CellValue<bool> is_cell_in_zone{mesh_3d->connectivity()};
+    is_cell_in_zone.fill(false);
+    auto zone_cell_list = mesh_cell_zone.cellList();
+    for (size_t i_cell = 0; i_cell < zone_cell_list.size(); ++i_cell) {
+      is_cell_in_zone[zone_cell_list[i_cell]] = true;
+    }
+
+    auto xj = MeshDataManager::instance().getMeshData(*mesh_3d).xj();
+
+    std::string_view data = R"(
+import math;
+let B_scalar_non_linear_3d: R^3 -> B, x -> (exp(2 * x[0])< 2*x[1]+x[2]);
+let N_scalar_non_linear_3d: R^3 -> N, x -> floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + x[2] * x[2]);
+let Z_scalar_non_linear_3d: R^3 -> Z, x -> floor(exp(2 * x[0]) - 3 * x[1] + x[2]);
+let R_scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0]+x[2]) + 3 * x[1];
+let R1_non_linear_3d: R^3 -> R^1, x -> 2 * exp(x[0])+sin(x[1] + x[2]);
+let R2_non_linear_3d: R^3 -> R^2, x -> (2 * exp(x[0]), -3*x[1] * x[2]);
+let R3_non_linear_3d: R^3 -> R^3, x -> (2 * exp(x[0]) + 3, x[1] - 2, 3 * x[2]);
+let R1x1_non_linear_3d: R^3 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * x[2]);
+let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[2] - 2 * x[0]), 3, x[1] * x[0] - x[2]);
+let R3x3_non_linear_3d: R^3 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[2]), 3, x[1] * x[2], -4*x[1], 2*x[2]+1, 3, -6*x[2], exp(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("B_scalar_non_linear_3d")
+    {
+      auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_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_value{mesh_3d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::exp(2 * x[0]) < 2 * x[1] + x[2];
+          } else {
+            cell_value[cell_id] = false;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("N_scalar_non_linear_3d")
+    {
+      auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_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_value{mesh_3d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + x[2] * x[2]);
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("Z_scalar_non_linear_3d")
+    {
+      auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_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_value{mesh_3d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::floor(std::exp(2 * x[0]) - 3 * x[1] + x[2]);
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("R_scalar_non_linear_3d")
+    {
+      auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_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_value{mesh_3d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = 2 * std::exp(x[0] + x[2]) + 3 * x[1];
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+    }
+
+    SECTION("R1_non_linear_3d")
+    {
+      using DataType = TinyVector<1>;
+
+      auto [i_symbol, found] = symbol_table->find("R1_non_linear_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<DataType> cell_value{mesh_3d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + std::sin(x[1] + x[2])};
+          } else {
+            cell_value[cell_id] = zero;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R2_non_linear_3d")
+    {
+      using DataType = TinyVector<2>;
+
+      auto [i_symbol, found] = symbol_table->find("R2_non_linear_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<DataType> cell_value{mesh_3d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0]), -3 * x[1] * x[2]};
+          } else {
+            cell_value[cell_id] = zero;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R3_non_linear_3d")
+    {
+      using DataType = TinyVector<3>;
+
+      auto [i_symbol, found] = symbol_table->find("R3_non_linear_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<DataType> cell_value{mesh_3d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + 3, x[1] - 2, 3 * x[2]};
+          } else {
+            cell_value[cell_id] = zero;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R1x1_non_linear_3d")
+    {
+      using DataType = TinyMatrix<1>;
+
+      auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_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<DataType> cell_value{mesh_3d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3 * x[2]};
+          } else {
+            cell_value[cell_id] = zero;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R2x2_non_linear_3d")
+    {
+      using DataType = TinyMatrix<2>;
+
+      auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_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<DataType> cell_value{mesh_3d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id] =
+              DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3, std::sin(x[2] - 2 * x[0]), 3, x[1] * x[0] - x[2]};
+          } else {
+            cell_value[cell_id] = zero;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+
+    SECTION("R3x3_non_linear_3d")
+    {
+      using DataType = TinyMatrix<3>;
+
+      auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_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<DataType> cell_value{mesh_3d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+
+            cell_value[cell_id] = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3,
+                                           std::sin(x[1] - 2 * x[2]),
+                                           3,
+                                           x[1] * x[2],
+                                           -4 * x[1],
+                                           2 * x[2] + 1,
+                                           3,
+                                           -6 * x[2],
+                                           std::exp(x[1] + x[2])};
+          } else {
+            cell_value[cell_id] = zero;
+          }
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                            function_symbol_id);
+      std::shared_ptr discrete_function = interpoler.interpolate();
+
+      REQUIRE(
+        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+    }
+  }
+}
diff --git a/tests/test_DiscreteFunctionVectorIntegrator.cpp b/tests/test_DiscreteFunctionVectorIntegrator.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d0a9d6a78fd260c0b259aa630d9e79f85271ae89
--- /dev/null
+++ b/tests/test_DiscreteFunctionVectorIntegrator.cpp
@@ -0,0 +1,410 @@
+#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/IntegrateCellValue.hpp>
+#include <language/utils/PugsFunctionAdapter.hpp>
+#include <language/utils/SymbolTable.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <analysis/IQuadratureDescriptor.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+
+#include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
+#include <scheme/DiscreteFunctionVectorIntegrator.hpp>
+
+#include <pegtl/string_input.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("DiscreteFunctionVectorIntegrator", "[scheme]")
+{
+  auto same_cell_value = [](const CellValue<const double>& fi, const size_t i, const auto& f) -> bool {
+    for (CellId cell_id = 0; cell_id < fi.numberOfItems(); ++cell_id) {
+      if (fi[cell_id] != f[cell_id][i]) {
+        return false;
+      }
+    }
+
+    return true;
+  };
+
+  auto register_function = [](const TAO_PEGTL_NAMESPACE::position& position,
+                              const std::shared_ptr<SymbolTable>& symbol_table, const std::string& name,
+                              std::vector<FunctionSymbolId>& function_id_list) {
+    auto [i_symbol, found] = symbol_table->find(name, 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);
+    function_id_list.push_back(function_symbol_id);
+  };
+
+  SECTION("1D")
+  {
+    constexpr size_t Dimension = 1;
+
+    std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+    std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+    for (auto named_mesh : mesh_list) {
+      SECTION(named_mesh.name())
+      {
+        auto mesh_1d = named_mesh.mesh();
+
+        std::string_view data = R"(
+import math;
+let B_scalar_non_linear_1d: R^1 -> B, x -> (exp(2 * x[0]) + 3 > 4);
+let N_scalar_non_linear_1d: R^1 -> N, x -> floor(3 * x[0] * x[0] + 2);
+let Z_scalar_non_linear_1d: R^1 -> Z, x -> floor(exp(2 * x[0]) - 1);
+let R_scalar_non_linear_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+)";
+        TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+        auto ast = ASTBuilder::build(input);
+
+        ASTModulesImporter{*ast};
+        ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+        ASTSymbolTableBuilder{*ast};
+        ASTNodeDataTypeBuilder{*ast};
+
+        ASTNodeTypeCleaner<language::var_declaration>{*ast};
+        ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+        ASTNodeExpressionBuilder{*ast};
+
+        TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+        position.byte = data.size();   // ensure that variables are declared at this point
+
+        std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+        std::vector<FunctionSymbolId> function_id_list;
+        register_function(position, symbol_table, "B_scalar_non_linear_1d", function_id_list);
+        register_function(position, symbol_table, "N_scalar_non_linear_1d", function_id_list);
+        register_function(position, symbol_table, "Z_scalar_non_linear_1d", function_id_list);
+        register_function(position, symbol_table, "R_scalar_non_linear_1d", function_id_list);
+
+        DiscreteFunctionVectorIntegrator integrator(mesh_1d, quadrature_descriptor,
+                                                    std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+                                                    function_id_list);
+        std::shared_ptr discrete_function = integrator.integrate();
+
+        size_t i = 0;
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<1>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_1d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<1>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_1d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<1>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_1d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<1>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_1d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        REQUIRE(i == function_id_list.size());
+      }
+    }
+  }
+
+  SECTION("2D")
+  {
+    constexpr size_t Dimension = 2;
+
+    std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+    std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+    for (auto named_mesh : mesh_list) {
+      SECTION(named_mesh.name())
+      {
+        auto mesh_2d = named_mesh.mesh();
+
+        std::string_view data = R"(
+import math;
+let B_scalar_non_linear_2d: R^2 -> B, x -> (exp(2 * x[0]) + 3 > 4);
+let N_scalar_non_linear_2d: R^2 -> N, x -> floor(3 * (x[0] * x[1]) * (x[0] * x[1]) + 2);
+let Z_scalar_non_linear_2d: R^2 -> Z, x -> floor(exp(2 * x[1]) - 1);
+let R_scalar_non_linear_2d: R^2 -> R, x -> 2 * exp(x[0] + x[1]) + 3;
+)";
+        TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+        auto ast = ASTBuilder::build(input);
+
+        ASTModulesImporter{*ast};
+        ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+        ASTSymbolTableBuilder{*ast};
+        ASTNodeDataTypeBuilder{*ast};
+
+        ASTNodeTypeCleaner<language::var_declaration>{*ast};
+        ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+        ASTNodeExpressionBuilder{*ast};
+
+        TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+        position.byte = data.size();   // ensure that variables are declared at this point
+
+        std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+        std::vector<FunctionSymbolId> function_id_list;
+        register_function(position, symbol_table, "B_scalar_non_linear_2d", function_id_list);
+        register_function(position, symbol_table, "N_scalar_non_linear_2d", function_id_list);
+        register_function(position, symbol_table, "Z_scalar_non_linear_2d", function_id_list);
+        register_function(position, symbol_table, "R_scalar_non_linear_2d", function_id_list);
+
+        DiscreteFunctionVectorIntegrator integrator(mesh_2d, quadrature_descriptor,
+                                                    std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+                                                    function_id_list);
+        std::shared_ptr discrete_function = integrator.integrate();
+
+        size_t i = 0;
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<2>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_2d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<2>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_2d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<2>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_2d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<2>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_2d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        REQUIRE(i == function_id_list.size());
+      }
+    }
+  }
+
+  SECTION("3D")
+  {
+    constexpr size_t Dimension = 3;
+
+    std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+    std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+    for (auto named_mesh : mesh_list) {
+      SECTION(named_mesh.name())
+      {
+        auto mesh_3d = named_mesh.mesh();
+
+        std::string_view data = R"(
+import math;
+let B_scalar_non_linear_3d: R^3 -> B, x -> (exp(2 * x[0] + x[2]) + 3 > 4);
+let N_scalar_non_linear_3d: R^3 -> N, x -> floor(3 * (x[0] * x[1]) * (x[0] * x[1]) + 2);
+let Z_scalar_non_linear_3d: R^3 -> Z, x -> floor(exp(2 * x[1]) - x[2]);
+let R_scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0] + x[1]) + 3 * 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};
+
+        TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+        position.byte = data.size();   // ensure that variables are declared at this point
+
+        std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+        std::vector<FunctionSymbolId> function_id_list;
+        register_function(position, symbol_table, "B_scalar_non_linear_3d", function_id_list);
+        register_function(position, symbol_table, "N_scalar_non_linear_3d", function_id_list);
+        register_function(position, symbol_table, "Z_scalar_non_linear_3d", function_id_list);
+        register_function(position, symbol_table, "R_scalar_non_linear_3d", function_id_list);
+
+        DiscreteFunctionVectorIntegrator integrator(mesh_3d, quadrature_descriptor,
+                                                    std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+                                                    function_id_list);
+        std::shared_ptr discrete_function = integrator.integrate();
+
+        size_t i = 0;
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<3>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_3d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<3>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_3d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<3>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_3d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<3>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_3d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        REQUIRE(i == function_id_list.size());
+      }
+    }
+  }
+
+  SECTION("errors")
+  {
+    std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+    std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+    for (auto named_mesh : mesh_list) {
+      SECTION(named_mesh.name())
+      {
+        auto mesh_3d = named_mesh.mesh();
+
+        std::string_view data = R"(
+import math;
+let B_scalar_non_linear_2d: R^2 -> B, x -> (exp(2 * x[0] + x[1]) + 3 > 4);
+let N_scalar_non_linear_1d: R^1 -> N, x -> floor(3 * x[0] * x[0] + 2);
+let Z_scalar_non_linear_3d: R^3 -> Z, x -> floor(exp(2 * x[1]) - x[2]);
+let R_scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0] + x[1]) + 3 * x[2];
+let R2_scalar_non_linear_3d: R^3 -> R^2, x -> (2 * exp(x[0] + x[1]) + 3 * x[2], 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};
+
+        TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+        position.byte = data.size();   // ensure that variables are declared at this point
+
+        std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+        SECTION("invalid function type")
+        {
+          std::vector<FunctionSymbolId> function_id_list;
+          register_function(position, symbol_table, "B_scalar_non_linear_2d", function_id_list);
+          register_function(position, symbol_table, "N_scalar_non_linear_1d", function_id_list);
+          register_function(position, symbol_table, "Z_scalar_non_linear_3d", function_id_list);
+          register_function(position, symbol_table, "R_scalar_non_linear_3d", function_id_list);
+
+          DiscreteFunctionVectorIntegrator integrator(mesh_3d, quadrature_descriptor,
+                                                      std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+                                                      function_id_list);
+
+          const std::string error_msg = R"(error: invalid function type
+note: expecting R^3 -> R
+note: provided function B_scalar_non_linear_2d: R^2 -> B)";
+
+          REQUIRE_THROWS_WITH(integrator.integrate(), error_msg);
+        }
+
+        SECTION("invalid value type")
+        {
+          std::vector<FunctionSymbolId> function_id_list;
+          register_function(position, symbol_table, "Z_scalar_non_linear_3d", function_id_list);
+          register_function(position, symbol_table, "R_scalar_non_linear_3d", function_id_list);
+          register_function(position, symbol_table, "R2_scalar_non_linear_3d", function_id_list);
+
+          DiscreteFunctionVectorIntegrator integrator(mesh_3d, quadrature_descriptor,
+                                                      std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+                                                      function_id_list);
+
+          const std::string error_msg = R"(error: vector functions require scalar value type.
+Invalid value type: R^2)";
+
+          REQUIRE_THROWS_WITH(integrator.integrate(), error_msg);
+        }
+      }
+    }
+  }
+}
diff --git a/tests/test_DiscreteFunctionVectorIntegratorByZone.cpp b/tests/test_DiscreteFunctionVectorIntegratorByZone.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..98111d4c6b7f1dd4788745db480e72c06c083945
--- /dev/null
+++ b/tests/test_DiscreteFunctionVectorIntegratorByZone.cpp
@@ -0,0 +1,440 @@
+#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/IntegrateCellValue.hpp>
+#include <language/utils/PugsFunctionAdapter.hpp>
+#include <language/utils/SymbolTable.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <analysis/IQuadratureDescriptor.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshCellZone.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/NamedZoneDescriptor.hpp>
+
+#include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
+#include <scheme/DiscreteFunctionVectorIntegrator.hpp>
+
+#include <pegtl/string_input.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("DiscreteFunctionVectorIntegratorByZone", "[scheme]")
+{
+  auto same_cell_value = [](const CellValue<const double>& fi, const size_t i, const auto& f) -> bool {
+    for (CellId cell_id = 0; cell_id < fi.numberOfItems(); ++cell_id) {
+      if (fi[cell_id] != f[cell_id][i]) {
+        return false;
+      }
+    }
+
+    return true;
+  };
+
+  auto register_function = [](const TAO_PEGTL_NAMESPACE::position& position,
+                              const std::shared_ptr<SymbolTable>& symbol_table, const std::string& name,
+                              std::vector<FunctionSymbolId>& function_id_list) {
+    auto [i_symbol, found] = symbol_table->find(name, 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);
+    function_id_list.push_back(function_symbol_id);
+  };
+
+  SECTION("1D")
+  {
+    constexpr size_t Dimension = 1;
+
+    std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+    auto mesh_1d = MeshDataBaseForTests::get().unordered1DMesh();
+
+    std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+    zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+    auto mesh_cell_zone = getMeshCellZone(*mesh_1d, *zone_list[0]);
+    auto zone_cell_list = mesh_cell_zone.cellList();
+
+    std::string_view data = R"(
+import math;
+let B_scalar_non_linear_1d: R^1 -> B, x -> (exp(2 * x[0]) + 3 < 3.3);
+let N_scalar_non_linear_1d: R^1 -> N, x -> floor(3 * x[0] * x[0] + 2);
+let Z_scalar_non_linear_1d: R^1 -> Z, x -> floor(exp(2 * x[0]) - 1);
+let R_scalar_non_linear_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+)";
+    TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+    auto ast = ASTBuilder::build(input);
+
+    ASTModulesImporter{*ast};
+    ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+    ASTSymbolTableBuilder{*ast};
+    ASTNodeDataTypeBuilder{*ast};
+
+    ASTNodeTypeCleaner<language::var_declaration>{*ast};
+    ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+    ASTNodeExpressionBuilder{*ast};
+
+    TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+    position.byte = data.size();   // ensure that variables are declared at this point
+
+    std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+    std::vector<FunctionSymbolId> function_id_list;
+    register_function(position, symbol_table, "B_scalar_non_linear_1d", function_id_list);
+    register_function(position, symbol_table, "N_scalar_non_linear_1d", function_id_list);
+    register_function(position, symbol_table, "Z_scalar_non_linear_1d", function_id_list);
+    register_function(position, symbol_table, "R_scalar_non_linear_1d", function_id_list);
+
+    DiscreteFunctionVectorIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor,
+                                                std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+                                                function_id_list);
+    std::shared_ptr discrete_function = integrator.integrate();
+
+    size_t i = 0;
+
+    {
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<1>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_1d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_1d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+          const CellId cell_id = zone_cell_list[j];
+          cell_value[cell_id]  = array[j];
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    {
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<1>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_1d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_1d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+          const CellId cell_id = zone_cell_list[j];
+          cell_value[cell_id]  = array[j];
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    {
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<1>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_1d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_1d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+          const CellId cell_id = zone_cell_list[j];
+          cell_value[cell_id]  = array[j];
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    {
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<1>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_1d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_1d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+          const CellId cell_id = zone_cell_list[j];
+          cell_value[cell_id]  = array[j];
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    REQUIRE(i == function_id_list.size());
+  }
+
+  SECTION("2D")
+  {
+    constexpr size_t Dimension = 2;
+
+    std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+    auto mesh_2d = MeshDataBaseForTests::get().hybrid2DMesh();
+
+    std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+    zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+    auto mesh_cell_zone = getMeshCellZone(*mesh_2d, *zone_list[0]);
+    auto zone_cell_list = mesh_cell_zone.cellList();
+
+    std::string_view data = R"(
+import math;
+let B_scalar_non_linear_2d: R^2 -> B, x -> (exp(2 * x[0]) + 3 > 3.3);
+let N_scalar_non_linear_2d: R^2 -> N, x -> floor(3 * (x[0] * x[1]) * (x[0] * x[1]) + 2);
+let Z_scalar_non_linear_2d: R^2 -> Z, x -> floor(exp(2 * x[1]) - 1);
+let R_scalar_non_linear_2d: R^2 -> R, x -> 2 * exp(x[0] + x[1]) + 3;
+)";
+    TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+    auto ast = ASTBuilder::build(input);
+
+    ASTModulesImporter{*ast};
+    ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+    ASTSymbolTableBuilder{*ast};
+    ASTNodeDataTypeBuilder{*ast};
+
+    ASTNodeTypeCleaner<language::var_declaration>{*ast};
+    ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+    ASTNodeExpressionBuilder{*ast};
+
+    TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+    position.byte = data.size();   // ensure that variables are declared at this point
+
+    std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+    std::vector<FunctionSymbolId> function_id_list;
+    register_function(position, symbol_table, "B_scalar_non_linear_2d", function_id_list);
+    register_function(position, symbol_table, "N_scalar_non_linear_2d", function_id_list);
+    register_function(position, symbol_table, "Z_scalar_non_linear_2d", function_id_list);
+    register_function(position, symbol_table, "R_scalar_non_linear_2d", function_id_list);
+
+    DiscreteFunctionVectorIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor,
+                                                std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+                                                function_id_list);
+    std::shared_ptr discrete_function = integrator.integrate();
+
+    size_t i = 0;
+
+    {
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<2>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_2d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_2d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+          const CellId cell_id = zone_cell_list[j];
+          cell_value[cell_id]  = array[j];
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    {
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<2>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_2d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_2d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+          const CellId cell_id = zone_cell_list[j];
+          cell_value[cell_id]  = array[j];
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    {
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<2>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_2d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_2d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+          const CellId cell_id = zone_cell_list[j];
+          cell_value[cell_id]  = array[j];
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    {
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<2>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_2d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_2d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+          const CellId cell_id = zone_cell_list[j];
+          cell_value[cell_id]  = array[j];
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    REQUIRE(i == function_id_list.size());
+  }
+
+  SECTION("3D")
+  {
+    constexpr size_t Dimension = 3;
+
+    std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+    auto mesh_3d = MeshDataBaseForTests::get().hybrid3DMesh();
+
+    std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+    zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+    auto mesh_cell_zone = getMeshCellZone(*mesh_3d, *zone_list[0]);
+    auto zone_cell_list = mesh_cell_zone.cellList();
+
+    std::string_view data = R"(
+import math;
+let B_scalar_non_linear_3d: R^3 -> B, x -> (exp(2 * x[0] + x[2]) + 3 > 3.4);
+let N_scalar_non_linear_3d: R^3 -> N, x -> floor(3 * (x[0] * x[1]) * (x[0] * x[1]) + 2);
+let Z_scalar_non_linear_3d: R^3 -> Z, x -> floor(exp(2 * x[1]) - x[2]);
+let R_scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0] + x[1]) + 3 * 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};
+
+    TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+    position.byte = data.size();   // ensure that variables are declared at this point
+
+    std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+    std::vector<FunctionSymbolId> function_id_list;
+    register_function(position, symbol_table, "B_scalar_non_linear_3d", function_id_list);
+    register_function(position, symbol_table, "N_scalar_non_linear_3d", function_id_list);
+    register_function(position, symbol_table, "Z_scalar_non_linear_3d", function_id_list);
+    register_function(position, symbol_table, "R_scalar_non_linear_3d", function_id_list);
+
+    DiscreteFunctionVectorIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor,
+                                                std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+                                                function_id_list);
+    std::shared_ptr discrete_function = integrator.integrate();
+
+    size_t i = 0;
+
+    {
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<3>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_3d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_3d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+          const CellId cell_id = zone_cell_list[j];
+          cell_value[cell_id]  = array[j];
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    {
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<3>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_3d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_3d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+          const CellId cell_id = zone_cell_list[j];
+          cell_value[cell_id]  = array[j];
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    {
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<3>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_3d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_3d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+          const CellId cell_id = zone_cell_list[j];
+          cell_value[cell_id]  = array[j];
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    {
+      Array<double> array =
+        IntegrateCellValue<double(TinyVector<3>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_3d,
+                                                             zone_cell_list);
+
+      CellValue<double> cell_value{mesh_3d->connectivity()};
+      cell_value.fill(0);
+
+      parallel_for(
+        zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+          const CellId cell_id = zone_cell_list[j];
+          cell_value[cell_id]  = array[j];
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    REQUIRE(i == function_id_list.size());
+  }
+}
diff --git a/tests/test_DiscreteFunctionVectorInterpolerByZone.cpp b/tests/test_DiscreteFunctionVectorInterpolerByZone.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..172c68d4f59f318d0f36c315aaf4962e71195dbe
--- /dev/null
+++ b/tests/test_DiscreteFunctionVectorInterpolerByZone.cpp
@@ -0,0 +1,511 @@
+#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/Mesh.hpp>
+#include <mesh/MeshCellZone.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/NamedZoneDescriptor.hpp>
+
+#include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
+#include <scheme/DiscreteFunctionVectorInterpoler.hpp>
+
+#include <pegtl/string_input.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("DiscreteFunctionVectorInterpolerByZone", "[scheme]")
+{
+  auto same_cell_value = [](const CellValue<const double>& fi, const size_t i, const auto& f) -> bool {
+    for (CellId cell_id = 0; cell_id < fi.numberOfItems(); ++cell_id) {
+      if (fi[cell_id] != f[cell_id][i]) {
+        return false;
+      }
+    }
+
+    return true;
+  };
+
+  auto register_function = [](const TAO_PEGTL_NAMESPACE::position& position,
+                              const std::shared_ptr<SymbolTable>& symbol_table, const std::string& name,
+                              std::vector<FunctionSymbolId>& function_id_list) {
+    auto [i_symbol, found] = symbol_table->find(name, 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);
+    function_id_list.push_back(function_symbol_id);
+  };
+
+  SECTION("1D")
+  {
+    constexpr size_t Dimension = 1;
+
+    auto mesh_1d = MeshDataBaseForTests::get().unordered1DMesh();
+
+    std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+    zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+    auto mesh_cell_zone = getMeshCellZone(*mesh_1d, *zone_list[0]);
+    CellValue<bool> is_cell_in_zone{mesh_1d->connectivity()};
+    is_cell_in_zone.fill(false);
+    auto zone_cell_list = mesh_cell_zone.cellList();
+    for (size_t i_cell = 0; i_cell < zone_cell_list.size(); ++i_cell) {
+      is_cell_in_zone[zone_cell_list[i_cell]] = true;
+    }
+
+    auto xj = MeshDataManager::instance().getMeshData(*mesh_1d).xj();
+
+    std::string_view data = R"(
+import math;
+let B_scalar_non_linear_1d: R^1 -> B, x -> (exp(2 * x[0]) + 3 > 4);
+let N_scalar_non_linear_1d: R^1 -> N, x -> floor(3 * x[0] * x[0] + 2);
+let Z_scalar_non_linear_1d: R^1 -> Z, x -> floor(exp(2 * x[0]) - 1);
+let R_scalar_non_linear_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+)";
+    TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+    auto ast = ASTBuilder::build(input);
+
+    ASTModulesImporter{*ast};
+    ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+    ASTSymbolTableBuilder{*ast};
+    ASTNodeDataTypeBuilder{*ast};
+
+    ASTNodeTypeCleaner<language::var_declaration>{*ast};
+    ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+    ASTNodeExpressionBuilder{*ast};
+
+    TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+    position.byte = data.size();   // ensure that variables are declared at this point
+
+    std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+    std::vector<FunctionSymbolId> function_id_list;
+    register_function(position, symbol_table, "B_scalar_non_linear_1d", function_id_list);
+    register_function(position, symbol_table, "N_scalar_non_linear_1d", function_id_list);
+    register_function(position, symbol_table, "Z_scalar_non_linear_1d", function_id_list);
+    register_function(position, symbol_table, "R_scalar_non_linear_1d", function_id_list);
+
+    DiscreteFunctionVectorInterpoler interpoler(mesh_1d, zone_list,
+                                                std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+                                                function_id_list);
+    std::shared_ptr discrete_function = interpoler.interpolate();
+
+    size_t i = 0;
+
+    {
+      CellValue<double> cell_value{mesh_1d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::exp(2 * x[0]) + 3 > 4;
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    {
+      CellValue<double> cell_value{mesh_1d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::floor(3 * x[0] * x[0] + 2);
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    {
+      CellValue<double> cell_value{mesh_1d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::floor(std::exp(2 * x[0]) - 1);
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    {
+      CellValue<double> cell_value{mesh_1d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = 2 * std::exp(x[0]) + 3;
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    REQUIRE(i == function_id_list.size());
+  }
+
+  SECTION("2D")
+  {
+    constexpr size_t Dimension = 2;
+
+    auto mesh_2d = MeshDataBaseForTests::get().hybrid2DMesh();
+
+    std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+    zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+    auto mesh_cell_zone = getMeshCellZone(*mesh_2d, *zone_list[0]);
+    CellValue<bool> is_cell_in_zone{mesh_2d->connectivity()};
+    is_cell_in_zone.fill(false);
+    auto zone_cell_list = mesh_cell_zone.cellList();
+    for (size_t i_cell = 0; i_cell < zone_cell_list.size(); ++i_cell) {
+      is_cell_in_zone[zone_cell_list[i_cell]] = true;
+    }
+
+    auto xj = MeshDataManager::instance().getMeshData(*mesh_2d).xj();
+
+    std::string_view data = R"(
+import math;
+let B_scalar_non_linear_2d: R^2 -> B, x -> (exp(2 * x[0]) + 3 > 4);
+let N_scalar_non_linear_2d: R^2 -> N, x -> floor(3 * (x[0] * x[1]) * (x[0] * x[1]) + 2);
+let Z_scalar_non_linear_2d: R^2 -> Z, x -> floor(exp(2 * x[1]) - 1);
+let R_scalar_non_linear_2d: R^2 -> R, x -> 2 * exp(x[0] + x[1]) + 3;
+)";
+    TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+    auto ast = ASTBuilder::build(input);
+
+    ASTModulesImporter{*ast};
+    ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+    ASTSymbolTableBuilder{*ast};
+    ASTNodeDataTypeBuilder{*ast};
+
+    ASTNodeTypeCleaner<language::var_declaration>{*ast};
+    ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+    ASTNodeExpressionBuilder{*ast};
+
+    TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+    position.byte = data.size();   // ensure that variables are declared at this point
+
+    std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+    std::vector<FunctionSymbolId> function_id_list;
+    register_function(position, symbol_table, "B_scalar_non_linear_2d", function_id_list);
+    register_function(position, symbol_table, "N_scalar_non_linear_2d", function_id_list);
+    register_function(position, symbol_table, "Z_scalar_non_linear_2d", function_id_list);
+    register_function(position, symbol_table, "R_scalar_non_linear_2d", function_id_list);
+
+    DiscreteFunctionVectorInterpoler interpoler(mesh_2d, zone_list,
+                                                std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+                                                function_id_list);
+    std::shared_ptr discrete_function = interpoler.interpolate();
+
+    size_t i = 0;
+
+    {
+      CellValue<double> cell_value{mesh_2d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::exp(2 * x[0]) + 3 > 4;
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    {
+      CellValue<double> cell_value{mesh_2d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::floor(3 * (x[0] * x[1]) * (x[0] * x[1]) + 2);
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    {
+      CellValue<double> cell_value{mesh_2d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::floor(std::exp(2 * x[1]) - 1);
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    {
+      CellValue<double> cell_value{mesh_2d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = 2 * std::exp(x[0] + x[1]) + 3;
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    REQUIRE(i == function_id_list.size());
+  }
+
+  SECTION("3D")
+  {
+    constexpr size_t Dimension = 3;
+
+    auto mesh_3d = MeshDataBaseForTests::get().hybrid3DMesh();
+
+    std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+    zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+    auto mesh_cell_zone = getMeshCellZone(*mesh_3d, *zone_list[0]);
+    CellValue<bool> is_cell_in_zone{mesh_3d->connectivity()};
+    is_cell_in_zone.fill(false);
+    auto zone_cell_list = mesh_cell_zone.cellList();
+    for (size_t i_cell = 0; i_cell < zone_cell_list.size(); ++i_cell) {
+      is_cell_in_zone[zone_cell_list[i_cell]] = true;
+    }
+
+    auto xj = MeshDataManager::instance().getMeshData(*mesh_3d).xj();
+
+    std::string_view data = R"(
+import math;
+let B_scalar_non_linear_3d: R^3 -> B, x -> (exp(2 * x[0] + x[2]) + 3 > 4);
+let N_scalar_non_linear_3d: R^3 -> N, x -> floor(3 * (x[0] * x[1]) * (x[0] * x[1]) + 2);
+let Z_scalar_non_linear_3d: R^3 -> Z, x -> floor(exp(2 * x[1]) - x[2]);
+let R_scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0] + x[1]) + 3 * 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};
+
+    TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+    position.byte = data.size();   // ensure that variables are declared at this point
+
+    std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+    std::vector<FunctionSymbolId> function_id_list;
+    register_function(position, symbol_table, "B_scalar_non_linear_3d", function_id_list);
+    register_function(position, symbol_table, "N_scalar_non_linear_3d", function_id_list);
+    register_function(position, symbol_table, "Z_scalar_non_linear_3d", function_id_list);
+    register_function(position, symbol_table, "R_scalar_non_linear_3d", function_id_list);
+
+    DiscreteFunctionVectorInterpoler interpoler(mesh_3d, zone_list,
+                                                std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+                                                function_id_list);
+    std::shared_ptr discrete_function = interpoler.interpolate();
+
+    size_t i = 0;
+
+    {
+      CellValue<double> cell_value{mesh_3d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::exp(2 * x[0] + x[2]) + 3 > 4;
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    {
+      CellValue<double> cell_value{mesh_3d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::floor(3 * (x[0] * x[1]) * (x[0] * x[1]) + 2);
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    {
+      CellValue<double> cell_value{mesh_3d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::floor(std::exp(2 * x[1]) - x[2]);
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    {
+      CellValue<double> cell_value{mesh_3d->connectivity()};
+      parallel_for(
+        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+          if (is_cell_in_zone[cell_id]) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = 2 * std::exp(x[0] + x[1]) + 3 * x[2];
+          } else {
+            cell_value[cell_id] = 0;
+          }
+        });
+
+      REQUIRE(same_cell_value(cell_value, i++,
+                              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+    }
+
+    REQUIRE(i == function_id_list.size());
+  }
+
+  SECTION("errors")
+  {
+    std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+    for (auto named_mesh : mesh_list) {
+      SECTION(named_mesh.name())
+      {
+        auto mesh_3d = named_mesh.mesh();
+
+        auto xj = MeshDataManager::instance().getMeshData(*mesh_3d).xj();
+
+        std::string_view data = R"(
+import math;
+let B_scalar_non_linear_2d: R^2 -> B, x -> (exp(2 * x[0] + x[1]) + 3 > 4);
+let N_scalar_non_linear_1d: R^1 -> N, x -> floor(3 * x[0] * x[0] + 2);
+let Z_scalar_non_linear_3d: R^3 -> Z, x -> floor(exp(2 * x[1]) - x[2]);
+let R_scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0] + x[1]) + 3 * x[2];
+let R2_scalar_non_linear_3d: R^3 -> R^2, x -> (2 * exp(x[0] + x[1]) + 3 * x[2], 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};
+
+        TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+        position.byte = data.size();   // ensure that variables are declared at this point
+
+        std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+        SECTION("invalid function type")
+        {
+          std::vector<FunctionSymbolId> function_id_list;
+          register_function(position, symbol_table, "B_scalar_non_linear_2d", function_id_list);
+          register_function(position, symbol_table, "N_scalar_non_linear_1d", function_id_list);
+          register_function(position, symbol_table, "Z_scalar_non_linear_3d", function_id_list);
+          register_function(position, symbol_table, "R_scalar_non_linear_3d", function_id_list);
+
+          DiscreteFunctionVectorInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+                                                      function_id_list);
+
+          const std::string error_msg = R"(error: invalid function type
+note: expecting R^3 -> R
+note: provided function B_scalar_non_linear_2d: R^2 -> B)";
+
+          REQUIRE_THROWS_WITH(interpoler.interpolate(), error_msg);
+        }
+
+        SECTION("invalid value type")
+        {
+          std::vector<FunctionSymbolId> function_id_list;
+          register_function(position, symbol_table, "Z_scalar_non_linear_3d", function_id_list);
+          register_function(position, symbol_table, "R_scalar_non_linear_3d", function_id_list);
+          register_function(position, symbol_table, "R2_scalar_non_linear_3d", function_id_list);
+
+          DiscreteFunctionVectorInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+                                                      function_id_list);
+
+          const std::string error_msg = R"(error: vector functions require scalar value type.
+Invalid interpolation value type: R^2)";
+
+          REQUIRE_THROWS_WITH(interpoler.interpolate(), error_msg);
+        }
+
+        SECTION("invalid discrete function type")
+        {
+          const std::string error_msg = "error: invalid discrete function type for vector interpolation";
+
+          DiscreteFunctionVectorInterpoler interpoler{mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(), {}};
+          REQUIRE_THROWS_WITH(interpoler.interpolate(), error_msg);
+        }
+      }
+    }
+  }
+}
diff --git a/tests/test_SubItemValuePerItem.cpp b/tests/test_SubItemValuePerItem.cpp
index 08d421fc1333ae8344402f9806446ab22c8de3b4..65741d70089785bb58b34230569afcd33a961c62 100644
--- a/tests/test_SubItemValuePerItem.cpp
+++ b/tests/test_SubItemValuePerItem.cpp
@@ -63,7 +63,8 @@ TEST_CASE("SubItemValuePerItem", "[mesh]")
     };
 
     auto check_same_memory = [](const auto& sub_item_value_per_item, auto array) {
-      REQUIRE(sub_item_value_per_item.numberOfValues() == array.size());
+      const bool same_size = (sub_item_value_per_item.numberOfValues() == array.size());
+      REQUIRE(same_size);
       bool is_same = true;
       for (size_t i = 0; i < sub_item_value_per_item.numberOfValues(); ++i) {
         is_same &= (&(sub_item_value_per_item[i]) == &(array[i]));