From 3dff9889b0cf111e9fc1baf4574941103cfb6a0d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Wed, 1 Sep 2021 17:24:14 +0200
Subject: [PATCH] Add tests for DiscreteFunctionInterpoler

---
 src/scheme/DiscreteFunctionInterpoler.cpp |  40 +-
 src/scheme/DiscreteFunctionInterpoler.hpp |   2 +-
 tests/CMakeLists.txt                      |   1 +
 tests/test_DiscreteFunctionInterpoler.cpp | 916 ++++++++++++++++++++++
 4 files changed, 949 insertions(+), 10 deletions(-)
 create mode 100644 tests/test_DiscreteFunctionInterpoler.cpp

diff --git a/src/scheme/DiscreteFunctionInterpoler.cpp b/src/scheme/DiscreteFunctionInterpoler.cpp
index 11b19835e..c39b6ee9d 100644
--- a/src/scheme/DiscreteFunctionInterpoler.cpp
+++ b/src/scheme/DiscreteFunctionInterpoler.cpp
@@ -4,7 +4,7 @@
 #include <scheme/DiscreteFunctionP0.hpp>
 #include <utils/Exceptions.hpp>
 
-template <size_t Dimension, typename DataType>
+template <size_t Dimension, typename DataType, typename ValueType>
 std::shared_ptr<IDiscreteFunction>
 DiscreteFunctionInterpoler::_interpolate() const
 {
@@ -12,10 +12,25 @@ DiscreteFunctionInterpoler::_interpolate() const
   using MeshDataType      = MeshData<Dimension>;
   MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
 
-  return std::make_shared<
-    DiscreteFunctionP0<Dimension, DataType>>(mesh,
-                                             InterpolateItemValue<DataType(TinyVector<Dimension>)>::
-                                               template interpolate<ItemType::cell>(m_function_id, mesh_data.xj()));
+  if constexpr (std::is_same_v<DataType, ValueType>) {
+    return std::make_shared<
+      DiscreteFunctionP0<Dimension, ValueType>>(mesh,
+                                                InterpolateItemValue<DataType(TinyVector<Dimension>)>::
+                                                  template interpolate<ItemType::cell>(m_function_id, mesh_data.xj()));
+  } else {
+    static_assert(std::is_convertible_v<DataType, ValueType>);
+
+    CellValue<DataType> cell_data =
+      InterpolateItemValue<DataType(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(m_function_id,
+                                                                                                  mesh_data.xj());
+
+    CellValue<ValueType> cell_value{mesh->connectivity()};
+
+    parallel_for(
+      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);
+  }
 }
 
 template <size_t Dimension>
@@ -29,13 +44,13 @@ DiscreteFunctionInterpoler::_interpolate() const
 
   switch (data_type) {
   case ASTNodeDataType::bool_t: {
-    return this->_interpolate<Dimension, bool>();
+    return this->_interpolate<Dimension, bool, double>();
   }
   case ASTNodeDataType::unsigned_int_t: {
-    return this->_interpolate<Dimension, uint64_t>();
+    return this->_interpolate<Dimension, uint64_t, double>();
   }
   case ASTNodeDataType::int_t: {
-    return this->_interpolate<Dimension, int64_t>();
+    return this->_interpolate<Dimension, int64_t, double>();
   }
   case ASTNodeDataType::double_t: {
     return this->_interpolate<Dimension, double>();
@@ -51,12 +66,14 @@ DiscreteFunctionInterpoler::_interpolate() const
     case 3: {
       return this->_interpolate<Dimension, TinyVector<3>>();
     }
+      // LCOV_EXCL_START
     default: {
       std::ostringstream os;
       os << "invalid vector dimension " << rang::fgB::red << data_type.dimension() << rang::style::reset;
 
       throw UnexpectedError(os.str());
     }
+      // LCOV_EXCL_STOP
     }
   }
   case ASTNodeDataType::matrix_t: {
@@ -71,20 +88,24 @@ DiscreteFunctionInterpoler::_interpolate() const
     case 3: {
       return this->_interpolate<Dimension, TinyMatrix<3>>();
     }
+      // LCOV_EXCL_START
     default: {
       std::ostringstream os;
       os << "invalid vector dimension " << rang::fgB::red << data_type.dimension() << rang::style::reset;
 
       throw UnexpectedError(os.str());
     }
+      // LCOV_EXCL_STOP
     }
   }
+    // LCOV_EXCL_START
   default: {
     std::ostringstream os;
     os << "invalid interpolation value type: " << rang::fgB::red << dataTypeName(data_type) << rang::style::reset;
 
     throw UnexpectedError(os.str());
   }
+    // LCOV_EXCL_STOP
   }
 }
 
@@ -102,9 +123,10 @@ DiscreteFunctionInterpoler::interpolate() const
   case 3: {
     return this->_interpolate<3>();
   }
+    // LCOV_EXCL_START
   default: {
     throw UnexpectedError("invalid dimension");
   }
+    // LCOV_EXCL_STOP
   }
-  return nullptr;
 }
diff --git a/src/scheme/DiscreteFunctionInterpoler.hpp b/src/scheme/DiscreteFunctionInterpoler.hpp
index 758ff2454..e5724e3fe 100644
--- a/src/scheme/DiscreteFunctionInterpoler.hpp
+++ b/src/scheme/DiscreteFunctionInterpoler.hpp
@@ -15,7 +15,7 @@ class DiscreteFunctionInterpoler
   std::shared_ptr<const IDiscreteFunctionDescriptor> m_discrete_function_descriptor;
   const FunctionSymbolId m_function_id;
 
-  template <size_t Dimension, typename DataType>
+  template <size_t Dimension, typename DataType, typename ValueType = DataType>
   std::shared_ptr<IDiscreteFunction> _interpolate() const;
 
   template <size_t Dimension>
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index b107e3e02..b52aec0b9 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -108,6 +108,7 @@ add_executable (unit_tests
 
 add_executable (mpi_unit_tests
   mpi_test_main.cpp
+  test_DiscreteFunctionInterpoler.cpp
   test_DiscreteFunctionP0.cpp
   test_DiscreteFunctionP0Vector.cpp
   test_InterpolateItemArray.cpp
diff --git a/tests/test_DiscreteFunctionInterpoler.cpp b/tests/test_DiscreteFunctionInterpoler.cpp
new file mode 100644
index 000000000..d24bae5f0
--- /dev/null
+++ b/tests/test_DiscreteFunctionInterpoler.cpp
@@ -0,0 +1,916 @@
+#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 <scheme/DiscreteFunctionDescriptorP0.hpp>
+#include <scheme/DiscreteFunctionInterpoler.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+
+#include <pegtl/string_input.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("DiscreteFunctionInterpoler", "[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;
+
+    const auto& mesh_1d = MeshDataBaseForTests::get().cartesianMesh1D();
+    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;
+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) {
+          const TinyVector<Dimension>& x = xj[cell_id];
+          cell_value[cell_id]            = std::exp(2 * x[0]) + 3 > 4;
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, 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) {
+          const TinyVector<Dimension>& x = xj[cell_id];
+          cell_value[cell_id]            = std::floor(3 * x[0] * x[0] + 2);
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, 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) {
+          const TinyVector<Dimension>& x = xj[cell_id];
+          cell_value[cell_id]            = std::floor(std::exp(2 * x[0]) - 1);
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, 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) {
+          const TinyVector<Dimension>& x = xj[cell_id];
+          cell_value[cell_id]            = 2 * std::exp(x[0]) + 3;
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, 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) {
+          const TinyVector<Dimension>& x = xj[cell_id];
+          cell_value[cell_id]            = DataType{2 * std::exp(x[0])};
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, 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) {
+          const TinyVector<Dimension>& x = xj[cell_id];
+          cell_value[cell_id]            = DataType{2 * std::exp(x[0]), -3 * x[0]};
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, 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) {
+          const TinyVector<Dimension>& x = xj[cell_id];
+          cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + 3, x[0] - 2, 3};
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, 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) {
+          const TinyVector<Dimension>& x = xj[cell_id];
+          cell_value[cell_id]            = DataType{2 * std::exp(x[0]) * std::sin(x[0]) + 3};
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, 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) {
+          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]};
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, 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) {
+          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])};
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_1d, 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;
+
+    const auto& mesh_2d = MeshDataBaseForTests::get().cartesianMesh2D();
+    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) {
+          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)));
+    }
+  }
+
+  SECTION("3D")
+  {
+    constexpr size_t Dimension = 3;
+
+    const auto& mesh_3d = MeshDataBaseForTests::get().cartesianMesh3D();
+    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) {
+          const TinyVector<Dimension>& x = xj[cell_id];
+          cell_value[cell_id]            = std::exp(2 * x[0]) < 2 * x[1] + x[2];
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, 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) {
+          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]);
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, 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) {
+          const TinyVector<Dimension>& x = xj[cell_id];
+          cell_value[cell_id]            = std::floor(std::exp(2 * x[0]) - 3 * x[1] + x[2]);
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, 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) {
+          const TinyVector<Dimension>& x = xj[cell_id];
+          cell_value[cell_id]            = 2 * std::exp(x[0] + x[2]) + 3 * x[1];
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, 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) {
+          const TinyVector<Dimension>& x = xj[cell_id];
+          cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + std::sin(x[1] + x[2])};
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, 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) {
+          const TinyVector<Dimension>& x = xj[cell_id];
+          cell_value[cell_id]            = DataType{2 * std::exp(x[0]), -3 * x[1] * x[2]};
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, 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) {
+          const TinyVector<Dimension>& x = xj[cell_id];
+          cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + 3, x[1] - 2, 3 * x[2]};
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, 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) {
+          const TinyVector<Dimension>& x = xj[cell_id];
+          cell_value[cell_id]            = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3 * x[2]};
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, 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) {
+          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]};
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, 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) {
+          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])};
+        });
+
+      DiscreteFunctionInterpoler interpoler(mesh_3d, 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)));
+    }
+  }
+}
-- 
GitLab