diff --git a/tests/MeshDataBaseForTests.cpp b/tests/MeshDataBaseForTests.cpp
index dff010102d6470df3e75c6d0d303963bb602bf2d..7afec8326a85a2b2e5f2e5f5ce46c2aae1f2b953 100644
--- a/tests/MeshDataBaseForTests.cpp
+++ b/tests/MeshDataBaseForTests.cpp
@@ -1,8 +1,13 @@
 #include <MeshDataBaseForTests.hpp>
 #include <mesh/CartesianMeshBuilder.hpp>
 #include <mesh/Connectivity.hpp>
+#include <mesh/GmshReader.hpp>
+#include <utils/Messenger.hpp>
 #include <utils/PugsAssert.hpp>
 
+#include <filesystem>
+#include <fstream>
+
 const MeshDataBaseForTests* MeshDataBaseForTests::m_instance = nullptr;
 
 MeshDataBaseForTests::MeshDataBaseForTests()
@@ -15,6 +20,8 @@ MeshDataBaseForTests::MeshDataBaseForTests()
 
   m_cartesian_3d_mesh = std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(
     CartesianMeshBuilder{TinyVector<3>{0, 1, 0}, TinyVector<3>{2, -1, 3}, TinyVector<3, size_t>{6, 7, 4}}.mesh());
+
+  m_hybrid_2d_mesh = _buildHybrid2dMesh();
 }
 
 const MeshDataBaseForTests&
@@ -55,3 +62,181 @@ MeshDataBaseForTests::cartesian3DMesh() const
 {
   return m_cartesian_3d_mesh;
 }
+
+std::shared_ptr<const Mesh<Connectivity<2>>>
+MeshDataBaseForTests::hybrid2DMesh() const
+{
+  return m_hybrid_2d_mesh;
+}
+
+std::shared_ptr<const Mesh<Connectivity<2>>>
+MeshDataBaseForTests::_buildHybrid2dMesh()
+{
+  const std::string filename = std::filesystem::path{PUGS_BINARY_DIR}.append("tests").append("hybrid-2d.msh");
+  if (parallel::rank() == 0) {
+    std::ofstream fout(filename);
+    fout << R"($MeshFormat
+2.2 0 8
+$EndMeshFormat
+$PhysicalNames
+10
+0 7 "XMINYMIN"
+0 8 "XMINYMAX"
+0 9 "XMAXYMIN"
+0 10 "XMAXYMAX"
+1 1 "XMIN"
+1 2 "XMAX"
+1 3 "YMAX"
+1 4 "YMIN"
+2 5 "LEFT"
+2 6 "RIGHT"
+$EndPhysicalNames
+$Nodes
+53
+1 0 0 0
+2 1 0 0
+3 1 1 0
+4 0 1 0
+5 2 0 0
+6 2 1 0
+7 0 0.7500000000003477 0
+8 0 0.5000000000020616 0
+9 0 0.2500000000010419 0
+10 1 0.2499999999994109 0
+11 1 0.4999999999986918 0
+12 1 0.7499999999993401 0
+13 2 0.2499999999994109 0
+14 2 0.4999999999986918 0
+15 2 0.7499999999993401 0
+16 0.7500000000003477 1 0
+17 0.5000000000020616 1 0
+18 0.2500000000010419 1 0
+19 1.749999999999999 1 0
+20 1.499999999999998 1 0
+21 1.249999999999999 1 0
+22 0.2499999999994109 0 0
+23 0.4999999999986918 0 0
+24 0.7499999999993401 0 0
+25 1.250000000000001 0 0
+26 1.500000000000002 0 0
+27 1.750000000000001 0 0
+28 0.5645539210988633 0.7822119881665017 0
+29 0.8034002623653069 0.616177001724073 0
+30 0.6429040982169911 0.5281266166868597 0
+31 0.8070137458488089 0.4050273873671746 0
+32 0.193468206450602 0.7967828567280011 0
+33 0.192644796000149 0.6419508693748935 0
+34 0.3924123452244094 0.5739639975328588 0
+35 0.3555609898784376 0.7915857937795373 0
+36 0.1932938945983771 0.4291816994842411 0
+37 0.3837729988924632 0.3582605556259788 0
+38 0.5797710532625071 0.2819566416178859 0
+39 0.7822395465040892 0.7780567773069014 0
+40 0.355038784333282 0.2031847484281184 0
+41 0.2059222030287144 0.2167595959294507 0
+42 0.7859869002303115 0.2399319972346836 0
+43 1.212848849190124 0.2732369527435796 0
+44 1.307250451008527 0.5758879308566747 0
+45 1.493594137086004 0.410605759034483 0
+46 1.439097356236337 0.1761649793157971 0
+47 1.574782862067863 0.7190230111887566 0
+48 1.740502300977129 0.4934083027495074 0
+49 1.853320217111459 0.2440179674328906 0
+50 1.629419001691624 0.2206995018497341 0
+51 1.844214193840832 0.743738552322716 0
+52 1.151907331010239 0.7762099486172758 0
+53 1.354193535052909 0.8313717608489495 0
+$EndNodes
+$Elements
+86
+1 15 2 7 1 1
+2 15 2 8 4 4
+3 15 2 9 5 5
+4 15 2 10 6 6
+5 1 2 1 1 4 7
+6 1 2 1 1 7 8
+7 1 2 1 1 8 9
+8 1 2 1 1 9 1
+9 1 2 2 3 5 13
+10 1 2 2 3 13 14
+11 1 2 2 3 14 15
+12 1 2 2 3 15 6
+13 1 2 3 4 3 16
+14 1 2 3 4 16 17
+15 1 2 3 4 17 18
+16 1 2 3 4 18 4
+17 1 2 3 5 6 19
+18 1 2 3 5 19 20
+19 1 2 3 5 20 21
+20 1 2 3 5 21 3
+21 1 2 4 6 1 22
+22 1 2 4 6 22 23
+23 1 2 4 6 23 24
+24 1 2 4 6 24 2
+25 1 2 4 7 2 25
+26 1 2 4 7 25 26
+27 1 2 4 7 26 27
+28 1 2 4 7 27 5
+29 2 2 6 2 43 45 44
+30 2 2 6 2 46 45 43
+31 2 2 6 2 11 43 44
+32 2 2 6 2 45 48 47
+33 2 2 6 2 46 43 25
+34 2 2 6 2 44 45 47
+35 2 2 6 2 43 11 10
+36 2 2 6 2 48 50 49
+37 2 2 6 2 14 51 48
+38 2 2 6 2 45 50 48
+39 2 2 6 2 27 49 50
+40 2 2 6 2 52 11 44
+41 2 2 6 2 49 27 5
+42 2 2 6 2 47 19 20
+43 2 2 6 2 12 52 3
+44 2 2 6 2 52 21 3
+45 2 2 6 2 47 53 44
+46 2 2 6 2 27 50 26
+47 2 2 6 2 43 10 2
+48 2 2 6 2 25 43 2
+49 2 2 6 2 20 53 47
+50 2 2 6 2 21 53 20
+51 2 2 6 2 46 50 45
+52 2 2 6 2 26 50 46
+53 2 2 6 2 44 53 52
+54 2 2 6 2 46 25 26
+55 2 2 6 2 11 52 12
+56 2 2 6 2 47 48 51
+57 2 2 6 2 51 14 15
+58 2 2 6 2 13 49 5
+59 2 2 6 2 51 19 47
+60 2 2 6 2 21 52 53
+61 2 2 6 2 19 51 6
+62 2 2 6 2 51 15 6
+63 2 2 6 2 48 49 14
+64 2 2 6 2 49 13 14
+65 3 2 5 1 37 34 33 36
+66 3 2 5 1 30 31 11 29
+67 3 2 5 1 35 34 30 28
+68 3 2 5 1 30 29 39 28
+69 3 2 5 1 34 37 38 30
+70 3 2 5 1 37 40 23 38
+71 3 2 5 1 41 40 37 36
+72 3 2 5 1 28 17 18 35
+73 3 2 5 1 30 38 42 31
+74 3 2 5 1 11 31 42 10
+75 3 2 5 1 39 16 17 28
+76 3 2 5 1 33 32 4 7
+77 3 2 5 1 11 12 39 29
+78 3 2 5 1 35 18 4 32
+79 3 2 5 1 23 40 41 22
+80 3 2 5 1 16 39 12 3
+81 3 2 5 1 41 36 8 9
+82 3 2 5 1 9 1 22 41
+83 3 2 5 1 10 42 24 2
+84 3 2 5 1 23 24 42 38
+85 3 2 5 1 7 8 36 33
+86 3 2 5 1 35 32 33 34
+$EndElements
+)";
+  }
+  return std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(GmshReader{filename}.mesh());
+}
diff --git a/tests/MeshDataBaseForTests.hpp b/tests/MeshDataBaseForTests.hpp
index 09ab816f52aef55c207ae1bea52898c14fdf2e8f..4c34019134f36ce6dbaaffe4528e20c808942926 100644
--- a/tests/MeshDataBaseForTests.hpp
+++ b/tests/MeshDataBaseForTests.hpp
@@ -22,11 +22,16 @@ class MeshDataBaseForTests
   std::shared_ptr<const Mesh<Connectivity<2>>> m_cartesian_2d_mesh;
   std::shared_ptr<const Mesh<Connectivity<3>>> m_cartesian_3d_mesh;
 
+  std::shared_ptr<const Mesh<Connectivity<2>>> m_hybrid_2d_mesh;
+
+  std::shared_ptr<const Mesh<Connectivity<2>>> _buildHybrid2dMesh();
+
  public:
   std::shared_ptr<const Mesh<Connectivity<1>>> cartesian1DMesh() const;
   std::shared_ptr<const Mesh<Connectivity<2>>> cartesian2DMesh() const;
   std::shared_ptr<const Mesh<Connectivity<3>>> cartesian3DMesh() const;
 
+  std::shared_ptr<const Mesh<Connectivity<2>>> hybrid2DMesh() const;
 
   static const MeshDataBaseForTests& get();
   static void create();
diff --git a/tests/test_DiscreteFunctionInterpoler.cpp b/tests/test_DiscreteFunctionInterpoler.cpp
index df0589d195aa0be51cadae6ce948d8206a46d177..ed7ae7c32dcab61467bb0bea2a975a0bd885baaf 100644
--- a/tests/test_DiscreteFunctionInterpoler.cpp
+++ b/tests/test_DiscreteFunctionInterpoler.cpp
@@ -334,10 +334,12 @@ let R3x3_non_linear_1d: R^1 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x
   {
     constexpr size_t Dimension = 2;
 
-    const auto& mesh_2d = MeshDataBaseForTests::get().cartesian2DMesh();
-    auto xj             = MeshDataManager::instance().getMeshData(*mesh_2d).xj();
+    SECTION("cartesian grid")
+    {
+      const auto& mesh_2d = MeshDataBaseForTests::get().cartesian2DMesh();
+      auto xj             = MeshDataManager::instance().getMeshData(*mesh_2d).xj();
 
-    std::string_view data = R"(
+      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);
@@ -350,275 +352,566 @@ 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"};
+      TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
 
-    auto ast = ASTBuilder::build(input);
+      auto ast = ASTBuilder::build(input);
 
-    ASTModulesImporter{*ast};
-    ASTNodeTypeCleaner<language::import_instruction>{*ast};
+      ASTModulesImporter{*ast};
+      ASTNodeTypeCleaner<language::import_instruction>{*ast};
 
-    ASTSymbolTableBuilder{*ast};
-    ASTNodeDataTypeBuilder{*ast};
+      ASTSymbolTableBuilder{*ast};
+      ASTNodeDataTypeBuilder{*ast};
 
-    ASTNodeTypeCleaner<language::var_declaration>{*ast};
-    ASTNodeTypeCleaner<language::fct_declaration>{*ast};
-    ASTNodeExpressionBuilder{*ast};
+      ASTNodeTypeCleaner<language::var_declaration>{*ast};
+      ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+      ASTNodeExpressionBuilder{*ast};
 
-    std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+      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
+      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);
+      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);
+        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) {
-          const TinyVector<Dimension>& x = xj[cell_id];
-          cell_value[cell_id]            = std::exp(2 * x[0]) < 2 * x[1];
-        });
+        CellValue<double> cell_value{mesh_2d->connectivity()};
+        parallel_for(
+          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::exp(2 * x[0]) < 2 * x[1];
+          });
 
-      DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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)));
-    }
+        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);
+      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);
+        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) {
-          const TinyVector<Dimension>& x = xj[cell_id];
-          cell_value[cell_id]            = std::floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + 2);
-        });
+        CellValue<double> cell_value{mesh_2d->connectivity()};
+        parallel_for(
+          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId 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);
+          });
 
-      DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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)));
-    }
+        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);
+      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);
+        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) {
-          const TinyVector<Dimension>& x = xj[cell_id];
-          cell_value[cell_id]            = std::floor(std::exp(2 * x[0]) - 3 * x[1]);
-        });
+        CellValue<double> cell_value{mesh_2d->connectivity()};
+        parallel_for(
+          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::floor(std::exp(2 * x[0]) - 3 * x[1]);
+          });
 
-      DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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)));
-    }
+        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);
+      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);
+        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) {
-          const TinyVector<Dimension>& x = xj[cell_id];
-          cell_value[cell_id]            = 2 * std::exp(x[0]) + 3 * x[1];
-        });
+        CellValue<double> cell_value{mesh_2d->connectivity()};
+        parallel_for(
+          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = 2 * std::exp(x[0]) + 3 * x[1];
+          });
 
-      DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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)));
-    }
+        REQUIRE(
+          same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+      }
 
-    SECTION("R1_non_linear_2d")
-    {
-      using DataType = TinyVector<1>;
+      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);
+        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);
+        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) {
-          const TinyVector<Dimension>& x = xj[cell_id];
-          cell_value[cell_id]            = DataType{2 * std::exp(x[0])};
-        });
+        CellValue<DataType> cell_value{mesh_2d->connectivity()};
+        parallel_for(
+          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0])};
+          });
 
-      DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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)));
-    }
+        REQUIRE(same_cell_value(cell_value,
+                                dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+      }
 
-    SECTION("R2_non_linear_2d")
-    {
-      using DataType = TinyVector<2>;
+      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);
+        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);
+        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) {
-          const TinyVector<Dimension>& x = xj[cell_id];
-          cell_value[cell_id]            = DataType{2 * std::exp(x[0]), -3 * x[1]};
-        });
+        CellValue<DataType> cell_value{mesh_2d->connectivity()};
+        parallel_for(
+          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0]), -3 * x[1]};
+          });
 
-      DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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)));
-    }
+        REQUIRE(same_cell_value(cell_value,
+                                dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+      }
 
-    SECTION("R3_non_linear_2d")
-    {
-      using DataType = TinyVector<3>;
+      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);
+        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);
+        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) {
-          const TinyVector<Dimension>& x = xj[cell_id];
-          cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + 3, x[1] - 2, 3};
-        });
+        CellValue<DataType> cell_value{mesh_2d->connectivity()};
+        parallel_for(
+          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + 3, x[1] - 2, 3};
+          });
 
-      DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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)));
-    }
+        REQUIRE(same_cell_value(cell_value,
+                                dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+      }
 
-    SECTION("R1x1_non_linear_2d")
-    {
-      using DataType = TinyMatrix<1>;
+      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);
+        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);
+        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) {
-          const TinyVector<Dimension>& x = xj[cell_id];
-          cell_value[cell_id]            = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3};
-        });
+        CellValue<DataType> cell_value{mesh_2d->connectivity()};
+        parallel_for(
+          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3};
+          });
 
-      DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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)));
-    }
+        REQUIRE(same_cell_value(cell_value,
+                                dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+      }
 
-    SECTION("R2x2_non_linear_2d")
-    {
-      using DataType = TinyMatrix<2>;
+      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);
+        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);
+        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) {
-          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]};
-        });
+        CellValue<DataType> cell_value{mesh_2d->connectivity()};
+        parallel_for(
+          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId 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]};
+          });
 
-      DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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)));
+        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) {
+            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])};
+          });
+
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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")
+    SECTION("hybrid grid")
     {
-      using DataType = TinyMatrix<3>;
+      const auto& mesh_2d = MeshDataBaseForTests::get().hybrid2DMesh();
+      auto xj             = MeshDataManager::instance().getMeshData(*mesh_2d).xj();
 
-      auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_2d", position);
-      REQUIRE(found);
-      REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+      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"};
 
-      FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+      auto ast = ASTBuilder::build(input);
 
-      CellValue<DataType> cell_value{mesh_2d->connectivity()};
-      parallel_for(
-        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-          const TinyVector<Dimension>& x = xj[cell_id];
+      ASTModulesImporter{*ast};
+      ASTNodeTypeCleaner<language::import_instruction>{*ast};
 
-          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])};
-        });
+      ASTSymbolTableBuilder{*ast};
+      ASTNodeDataTypeBuilder{*ast};
 
-      DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
+      ASTNodeTypeCleaner<language::var_declaration>{*ast};
+      ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+      ASTNodeExpressionBuilder{*ast};
 
-      REQUIRE(
-        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+      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) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::exp(2 * x[0]) < 2 * x[1];
+          });
+
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + 2);
+          });
+
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = std::floor(std::exp(2 * x[0]) - 3 * x[1]);
+          });
+
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = 2 * std::exp(x[0]) + 3 * x[1];
+          });
+
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0])};
+          });
+
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0]), -3 * x[1]};
+          });
+
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + 3, x[1] - 2, 3};
+          });
+
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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) {
+            const TinyVector<Dimension>& x = xj[cell_id];
+            cell_value[cell_id]            = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3};
+          });
+
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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) {
+            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]};
+          });
+
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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) {
+            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])};
+          });
+
+        DiscreteFunctionInterpoler interpoler(mesh_2d, 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)));
+      }
     }
   }