diff --git a/tests/test_InterpolateItemArray.cpp b/tests/test_InterpolateItemArray.cpp
index ccee6e77877eadc359cc2c11313a7114081cc8c0..878503f1156e167796ebd8d3dec4019bc76e59f7 100644
--- a/tests/test_InterpolateItemArray.cpp
+++ b/tests/test_InterpolateItemArray.cpp
@@ -48,69 +48,131 @@ TEST_CASE("InterpolateItemArray", "[language]")
 
       std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
 
-      for (const auto& named_mesh : mesh_list) {
-        SECTION(named_mesh.name())
-        {
-          auto mesh_1d = named_mesh.mesh();
+      SECTION("from list of -> R")
+      {
+        for (const auto& named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_1d = named_mesh.mesh();
 
-          auto xj = MeshDataManager::instance().getMeshData(*mesh_1d).xj();
+            auto xj = MeshDataManager::instance().getMeshData(*mesh_1d).xj();
 
-          std::string_view data = R"(
+            std::string_view data = R"(
 import math;
 let scalar_affine_1d: R^1 -> R, x -> 2*x[0] + 2;
 let scalar_non_linear_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
 )";
-          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
 
-          std::vector<FunctionSymbolId> function_symbol_id_list;
+            std::vector<FunctionSymbolId> function_symbol_id_list;
 
-          {
-            auto [i_symbol, found] = symbol_table->find("scalar_affine_1d", position);
-            REQUIRE(found);
-            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_affine_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              function_symbol_id_list.push_back(
+                FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            }
+
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_non_linear_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
 
-            function_symbol_id_list.push_back(
-              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+              function_symbol_id_list.push_back(
+                FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            }
+
+            CellArray<double> cell_array{mesh_1d->connectivity(), 2};
+            parallel_for(
+              cell_array.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+                const TinyVector<Dimension>& x = xj[cell_id];
+                cell_array[cell_id][0]         = 2 * x[0] + 2;
+                cell_array[cell_id][1]         = 2 * exp(x[0]) + 3;
+              });
+
+            CellArray<const double> interpolate_array =
+              InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj);
+
+            REQUIRE(same_cell_array(cell_array, interpolate_array));
           }
+        }
+      }
 
+      SECTION("from -> (R)")
+      {
+        for (const auto& named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
           {
-            auto [i_symbol, found] = symbol_table->find("scalar_non_linear_1d", position);
-            REQUIRE(found);
-            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+            auto mesh_1d = named_mesh.mesh();
 
-            function_symbol_id_list.push_back(
-              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
-          }
+            auto xj = MeshDataManager::instance().getMeshData(*mesh_1d).xj();
+
+            std::string_view data = R"(
+import math;
+let f_1d: R^1 -> (R), x -> (2*x[0] + 2, 2 * exp(x[0]) + 3);
+)";
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
 
-          CellArray<double> cell_array{mesh_1d->connectivity(), 2};
-          parallel_for(
-            cell_array.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-              const TinyVector<Dimension>& x = xj[cell_id];
-              cell_array[cell_id][0]         = 2 * x[0] + 2;
-              cell_array[cell_id][1]         = 2 * exp(x[0]) + 3;
-            });
+            std::vector<FunctionSymbolId> function_symbol_id_list;
 
-          CellArray<const double> interpolate_array =
-            InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj);
+            {
+              auto [i_symbol, found] = symbol_table->find("f_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
 
-          REQUIRE(same_cell_array(cell_array, interpolate_array));
+              function_symbol_id_list.push_back(
+                FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            }
+
+            CellArray<double> cell_array{mesh_1d->connectivity(), 2};
+            parallel_for(
+              cell_array.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+                const TinyVector<Dimension>& x = xj[cell_id];
+                cell_array[cell_id][0]         = 2 * x[0] + 2;
+                cell_array[cell_id][1]         = 2 * exp(x[0]) + 3;
+              });
+
+            CellArray<const double> interpolate_array =
+              InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj);
+
+            REQUIRE(same_cell_array(cell_array, interpolate_array));
+          }
         }
       }
     }
@@ -121,69 +183,131 @@ let scalar_non_linear_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
 
       std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
 
-      for (const auto& named_mesh : mesh_list) {
-        SECTION(named_mesh.name())
-        {
-          auto mesh_2d = named_mesh.mesh();
+      SECTION("from list of -> R")
+      {
+        for (const auto& named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_2d = named_mesh.mesh();
 
-          auto xj = MeshDataManager::instance().getMeshData(*mesh_2d).xj();
+            auto xj = MeshDataManager::instance().getMeshData(*mesh_2d).xj();
 
-          std::string_view data = R"(
+            std::string_view data = R"(
 import math;
 let scalar_affine_2d: R^2 -> R, x -> 2*x[0] + 3*x[1] + 2;
 let scalar_non_linear_2d: R^2 -> R, x -> 2*exp(x[0])*sin(x[1])+3;
 )";
-          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
 
-          std::vector<FunctionSymbolId> function_symbol_id_list;
+            std::vector<FunctionSymbolId> function_symbol_id_list;
 
-          {
-            auto [i_symbol, found] = symbol_table->find("scalar_affine_2d", position);
-            REQUIRE(found);
-            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_affine_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              function_symbol_id_list.push_back(
+                FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            }
 
-            function_symbol_id_list.push_back(
-              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_non_linear_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              function_symbol_id_list.push_back(
+                FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            }
+
+            CellArray<double> cell_array{mesh_2d->connectivity(), 2};
+            parallel_for(
+              cell_array.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+                const TinyVector<Dimension>& x = xj[cell_id];
+                cell_array[cell_id][0]         = 2 * x[0] + 3 * x[1] + 2;
+                cell_array[cell_id][1]         = 2 * exp(x[0]) * sin(x[1]) + 3;
+              });
+
+            CellArray<const double> interpolate_array =
+              InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj);
+
+            REQUIRE(same_cell_array(cell_array, interpolate_array));
           }
+        }
+      }
 
+      SECTION("from -> (R)")
+      {
+        for (const auto& named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
           {
-            auto [i_symbol, found] = symbol_table->find("scalar_non_linear_2d", position);
-            REQUIRE(found);
-            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+            auto mesh_2d = named_mesh.mesh();
 
-            function_symbol_id_list.push_back(
-              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
-          }
+            auto xj = MeshDataManager::instance().getMeshData(*mesh_2d).xj();
+
+            std::string_view data = R"(
+import math;
+let f_2d: R^2 -> (R), x -> (2*x[0] + 3*x[1] + 2, 2*exp(x[0])*sin(x[1])+3);
+)";
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
 
-          CellArray<double> cell_array{mesh_2d->connectivity(), 2};
-          parallel_for(
-            cell_array.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-              const TinyVector<Dimension>& x = xj[cell_id];
-              cell_array[cell_id][0]         = 2 * x[0] + 3 * x[1] + 2;
-              cell_array[cell_id][1]         = 2 * exp(x[0]) * sin(x[1]) + 3;
-            });
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
 
-          CellArray<const double> interpolate_array =
-            InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj);
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
 
-          REQUIRE(same_cell_array(cell_array, interpolate_array));
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            std::vector<FunctionSymbolId> function_symbol_id_list;
+
+            {
+              auto [i_symbol, found] = symbol_table->find("f_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              function_symbol_id_list.push_back(
+                FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            }
+
+            CellArray<double> cell_array{mesh_2d->connectivity(), 2};
+            parallel_for(
+              cell_array.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+                const TinyVector<Dimension>& x = xj[cell_id];
+                cell_array[cell_id][0]         = 2 * x[0] + 3 * x[1] + 2;
+                cell_array[cell_id][1]         = 2 * exp(x[0]) * sin(x[1]) + 3;
+              });
+
+            CellArray<const double> interpolate_array =
+              InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj);
+
+            REQUIRE(same_cell_array(cell_array, interpolate_array));
+          }
         }
       }
     }
@@ -194,69 +318,131 @@ let scalar_non_linear_2d: R^2 -> R, x -> 2*exp(x[0])*sin(x[1])+3;
 
       std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
 
-      for (const auto& named_mesh : mesh_list) {
-        SECTION(named_mesh.name())
-        {
-          auto mesh_3d = named_mesh.mesh();
+      SECTION("from list of -> R")
+      {
+        for (const auto& named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_3d = named_mesh.mesh();
 
-          auto xj = MeshDataManager::instance().getMeshData(*mesh_3d).xj();
+            auto xj = MeshDataManager::instance().getMeshData(*mesh_3d).xj();
 
-          std::string_view data = R"(
+            std::string_view data = R"(
 import math;
 let scalar_affine_3d: R^3 -> R, x -> 2 * x[0] + 3 * x[1] + 2 * x[2] - 1;
 let scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
 )";
-          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
 
-          std::vector<FunctionSymbolId> function_symbol_id_list;
+            std::vector<FunctionSymbolId> function_symbol_id_list;
 
-          {
-            auto [i_symbol, found] = symbol_table->find("scalar_affine_3d", position);
-            REQUIRE(found);
-            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_affine_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              function_symbol_id_list.push_back(
+                FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            }
+
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_non_linear_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              function_symbol_id_list.push_back(
+                FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            }
+
+            CellArray<double> cell_array{mesh_3d->connectivity(), 2};
+            parallel_for(
+              cell_array.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+                const TinyVector<Dimension>& x = xj[cell_id];
+                cell_array[cell_id][0]         = 2 * x[0] + 3 * x[1] + 2 * x[2] - 1;
+                cell_array[cell_id][1]         = 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
+              });
+
+            CellArray<const double> interpolate_array =
+              InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj);
 
-            function_symbol_id_list.push_back(
-              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            REQUIRE(same_cell_array(cell_array, interpolate_array));
           }
+        }
+      }
 
+      SECTION("from -> (R)")
+      {
+        for (const auto& named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
           {
-            auto [i_symbol, found] = symbol_table->find("scalar_non_linear_3d", position);
-            REQUIRE(found);
-            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+            auto mesh_3d = named_mesh.mesh();
 
-            function_symbol_id_list.push_back(
-              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
-          }
+            auto xj = MeshDataManager::instance().getMeshData(*mesh_3d).xj();
+
+            std::string_view data = R"(
+import math;
+let f_3d: R^3 -> (R), x -> (2 * x[0] + 3 * x[1] + 2 * x[2] - 1, 2 * exp(x[0]) * sin(x[1]) * x[2] + 3);
+)";
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
 
-          CellArray<double> cell_array{mesh_3d->connectivity(), 2};
-          parallel_for(
-            cell_array.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-              const TinyVector<Dimension>& x = xj[cell_id];
-              cell_array[cell_id][0]         = 2 * x[0] + 3 * x[1] + 2 * x[2] - 1;
-              cell_array[cell_id][1]         = 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
-            });
+            std::vector<FunctionSymbolId> function_symbol_id_list;
 
-          CellArray<const double> interpolate_array =
-            InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj);
+            {
+              auto [i_symbol, found] = symbol_table->find("f_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
 
-          REQUIRE(same_cell_array(cell_array, interpolate_array));
+              function_symbol_id_list.push_back(
+                FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            }
+
+            CellArray<double> cell_array{mesh_3d->connectivity(), 2};
+            parallel_for(
+              cell_array.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+                const TinyVector<Dimension>& x = xj[cell_id];
+                cell_array[cell_id][0]         = 2 * x[0] + 3 * x[1] + 2 * x[2] - 1;
+                cell_array[cell_id][1]         = 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
+              });
+
+            CellArray<const double> interpolate_array =
+              InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj);
+
+            REQUIRE(same_cell_array(cell_array, interpolate_array));
+          }
         }
       }
     }
@@ -281,77 +467,149 @@ let scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
 
       std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
 
-      for (const auto& named_mesh : mesh_list) {
-        SECTION(named_mesh.name())
-        {
-          auto mesh_1d = named_mesh.mesh();
+      SECTION("from list of -> R")
+      {
+        for (const auto& named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_1d = named_mesh.mesh();
 
-          auto xj = MeshDataManager::instance().getMeshData(*mesh_1d).xj();
+            auto xj = MeshDataManager::instance().getMeshData(*mesh_1d).xj();
 
-          Array<const CellId> cell_id_list = [&] {
-            Array<CellId> cell_ids{mesh_1d->numberOfCells() / 2};
-            for (size_t i_cell = 0; i_cell < cell_ids.size(); ++i_cell) {
-              cell_ids[i_cell] = static_cast<CellId>(2 * i_cell);
-            }
-            return cell_ids;
-          }();
+            Array<const CellId> cell_id_list = [&] {
+              Array<CellId> cell_ids{mesh_1d->numberOfCells() / 2};
+              for (size_t i_cell = 0; i_cell < cell_ids.size(); ++i_cell) {
+                cell_ids[i_cell] = static_cast<CellId>(2 * i_cell);
+              }
+              return cell_ids;
+            }();
 
-          std::string_view data = R"(
+            std::string_view data = R"(
 import math;
 let scalar_affine_1d: R^1 -> R, x -> 2*x[0] + 2;
 let scalar_non_linear_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
 )";
-          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
 
-          std::vector<FunctionSymbolId> function_symbol_id_list;
+            std::vector<FunctionSymbolId> function_symbol_id_list;
 
-          {
-            auto [i_symbol, found] = symbol_table->find("scalar_affine_1d", position);
-            REQUIRE(found);
-            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_affine_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              function_symbol_id_list.push_back(
+                FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            }
 
-            function_symbol_id_list.push_back(
-              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_non_linear_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              function_symbol_id_list.push_back(
+                FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            }
+
+            Table<double> cell_array{cell_id_list.size(), 2};
+            parallel_for(
+              cell_id_list.size(), PUGS_LAMBDA(const size_t i) {
+                const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+                cell_array[i][0]               = 2 * x[0] + 2;
+                cell_array[i][1]               = 2 * exp(x[0]) + 3;
+              });
+
+            Table<const double> interpolate_array =
+              InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj,
+                                                                               cell_id_list);
+
+            REQUIRE(same_cell_value(cell_array, interpolate_array));
           }
+        }
+      }
 
+      SECTION("from -> (R)")
+      {
+        for (const auto& named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
           {
-            auto [i_symbol, found] = symbol_table->find("scalar_non_linear_1d", position);
-            REQUIRE(found);
-            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+            auto mesh_1d = named_mesh.mesh();
 
-            function_symbol_id_list.push_back(
-              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
-          }
+            auto xj = MeshDataManager::instance().getMeshData(*mesh_1d).xj();
+
+            Array<const CellId> cell_id_list = [&] {
+              Array<CellId> cell_ids{mesh_1d->numberOfCells() / 2};
+              for (size_t i_cell = 0; i_cell < cell_ids.size(); ++i_cell) {
+                cell_ids[i_cell] = static_cast<CellId>(2 * i_cell);
+              }
+              return cell_ids;
+            }();
+
+            std::string_view data = R"(
+import math;
+let f_1d: R^1 -> (R), x -> (2*x[0] + 2, 2 * exp(x[0]) + 3);
+)";
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            std::vector<FunctionSymbolId> function_symbol_id_list;
+
+            {
+              auto [i_symbol, found] = symbol_table->find("f_1d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              function_symbol_id_list.push_back(
+                FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            }
 
-          Table<double> cell_array{cell_id_list.size(), 2};
-          parallel_for(
-            cell_id_list.size(), PUGS_LAMBDA(const size_t i) {
-              const TinyVector<Dimension>& x = xj[cell_id_list[i]];
-              cell_array[i][0]               = 2 * x[0] + 2;
-              cell_array[i][1]               = 2 * exp(x[0]) + 3;
-            });
+            Table<double> cell_array{cell_id_list.size(), 2};
+            parallel_for(
+              cell_id_list.size(), PUGS_LAMBDA(const size_t i) {
+                const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+                cell_array[i][0]               = 2 * x[0] + 2;
+                cell_array[i][1]               = 2 * exp(x[0]) + 3;
+              });
 
-          Table<const double> interpolate_array =
-            InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj, cell_id_list);
+            Table<const double> interpolate_array =
+              InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj,
+                                                                               cell_id_list);
 
-          REQUIRE(same_cell_value(cell_array, interpolate_array));
+            REQUIRE(same_cell_value(cell_array, interpolate_array));
+          }
         }
       }
     }
@@ -362,74 +620,143 @@ let scalar_non_linear_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
 
       std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
 
-      for (const auto& named_mesh : mesh_list) {
-        SECTION(named_mesh.name())
-        {
-          auto mesh_2d = named_mesh.mesh();
+      SECTION("from list of -> R")
+      {
+        for (const auto& named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_2d = named_mesh.mesh();
 
-          auto xj = MeshDataManager::instance().getMeshData(*mesh_2d).xj();
+            auto xj = MeshDataManager::instance().getMeshData(*mesh_2d).xj();
 
-          Array<CellId> cell_id_list{mesh_2d->numberOfCells() / 2};
-          for (size_t i_cell = 0; i_cell < cell_id_list.size(); ++i_cell) {
-            cell_id_list[i_cell] = static_cast<CellId>(2 * i_cell);
-          }
+            Array<CellId> cell_id_list{mesh_2d->numberOfCells() / 2};
+            for (size_t i_cell = 0; i_cell < cell_id_list.size(); ++i_cell) {
+              cell_id_list[i_cell] = static_cast<CellId>(2 * i_cell);
+            }
 
-          std::string_view data = R"(
+            std::string_view data = R"(
 import math;
 let scalar_affine_2d: R^2 -> R, x -> 2*x[0] + 3*x[1] + 2;
 let scalar_non_linear_2d: R^2 -> R, x -> 2*exp(x[0])*sin(x[1])+3;
 )";
-          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
 
-          std::vector<FunctionSymbolId> function_symbol_id_list;
+            std::vector<FunctionSymbolId> function_symbol_id_list;
 
-          {
-            auto [i_symbol, found] = symbol_table->find("scalar_affine_2d", position);
-            REQUIRE(found);
-            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_affine_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              function_symbol_id_list.push_back(
+                FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            }
+
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_non_linear_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              function_symbol_id_list.push_back(
+                FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            }
 
-            function_symbol_id_list.push_back(
-              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            Table<double> cell_array{cell_id_list.size(), 2};
+            parallel_for(
+              cell_id_list.size(), PUGS_LAMBDA(const size_t i) {
+                const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+                cell_array[i][0]               = 2 * x[0] + 3 * x[1] + 2;
+                cell_array[i][1]               = 2 * exp(x[0]) * sin(x[1]) + 3;
+              });
+
+            Table<const double> interpolate_array =
+              InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj,
+                                                                               cell_id_list);
+
+            REQUIRE(same_cell_value(cell_array, interpolate_array));
           }
+        }
+      }
 
+      SECTION("from -> (R)")
+      {
+        for (const auto& named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
           {
-            auto [i_symbol, found] = symbol_table->find("scalar_non_linear_2d", position);
-            REQUIRE(found);
-            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+            auto mesh_2d = named_mesh.mesh();
 
-            function_symbol_id_list.push_back(
-              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
-          }
+            auto xj = MeshDataManager::instance().getMeshData(*mesh_2d).xj();
+
+            Array<CellId> cell_id_list{mesh_2d->numberOfCells() / 2};
+            for (size_t i_cell = 0; i_cell < cell_id_list.size(); ++i_cell) {
+              cell_id_list[i_cell] = static_cast<CellId>(2 * i_cell);
+            }
+
+            std::string_view data = R"(
+import math;
+let f_2d: R^2 -> (R), x -> (2*x[0] + 3*x[1] + 2, 2*exp(x[0])*sin(x[1])+3);
+)";
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
 
-          Table<double> cell_array{cell_id_list.size(), 2};
-          parallel_for(
-            cell_id_list.size(), PUGS_LAMBDA(const size_t i) {
-              const TinyVector<Dimension>& x = xj[cell_id_list[i]];
-              cell_array[i][0]               = 2 * x[0] + 3 * x[1] + 2;
-              cell_array[i][1]               = 2 * exp(x[0]) * sin(x[1]) + 3;
-            });
+            auto ast = ASTBuilder::build(input);
 
-          Table<const double> interpolate_array =
-            InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj, cell_id_list);
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
 
-          REQUIRE(same_cell_value(cell_array, interpolate_array));
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
+
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            std::vector<FunctionSymbolId> function_symbol_id_list;
+
+            {
+              auto [i_symbol, found] = symbol_table->find("f_2d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              function_symbol_id_list.push_back(
+                FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            }
+
+            Table<double> cell_array{cell_id_list.size(), 2};
+            parallel_for(
+              cell_id_list.size(), PUGS_LAMBDA(const size_t i) {
+                const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+                cell_array[i][0]               = 2 * x[0] + 3 * x[1] + 2;
+                cell_array[i][1]               = 2 * exp(x[0]) * sin(x[1]) + 3;
+              });
+
+            Table<const double> interpolate_array =
+              InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj,
+                                                                               cell_id_list);
+
+            REQUIRE(same_cell_value(cell_array, interpolate_array));
+          }
         }
       }
     }
@@ -440,74 +767,143 @@ let scalar_non_linear_2d: R^2 -> R, x -> 2*exp(x[0])*sin(x[1])+3;
 
       std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
 
-      for (const auto& named_mesh : mesh_list) {
-        SECTION(named_mesh.name())
-        {
-          auto mesh_3d = named_mesh.mesh();
+      SECTION("from list of -> R")
+      {
+        for (const auto& named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_3d = named_mesh.mesh();
 
-          auto xj = MeshDataManager::instance().getMeshData(*mesh_3d).xj();
+            auto xj = MeshDataManager::instance().getMeshData(*mesh_3d).xj();
 
-          Array<CellId> cell_id_list{mesh_3d->numberOfCells() / 2};
-          for (size_t i_cell = 0; i_cell < cell_id_list.size(); ++i_cell) {
-            cell_id_list[i_cell] = static_cast<CellId>(2 * i_cell);
-          }
+            Array<CellId> cell_id_list{mesh_3d->numberOfCells() / 2};
+            for (size_t i_cell = 0; i_cell < cell_id_list.size(); ++i_cell) {
+              cell_id_list[i_cell] = static_cast<CellId>(2 * i_cell);
+            }
 
-          std::string_view data = R"(
+            std::string_view data = R"(
 import math;
 let scalar_affine_3d: R^3 -> R, x -> 2 * x[0] + 3 * x[1] + 2 * x[2] - 1;
 let scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
 )";
-          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
 
-          std::vector<FunctionSymbolId> function_symbol_id_list;
+            std::vector<FunctionSymbolId> function_symbol_id_list;
 
-          {
-            auto [i_symbol, found] = symbol_table->find("scalar_affine_3d", position);
-            REQUIRE(found);
-            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_affine_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              function_symbol_id_list.push_back(
+                FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            }
+
+            {
+              auto [i_symbol, found] = symbol_table->find("scalar_non_linear_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              function_symbol_id_list.push_back(
+                FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            }
+
+            Table<double> cell_array{cell_id_list.size(), 2};
+            parallel_for(
+              cell_id_list.size(), PUGS_LAMBDA(const size_t i) {
+                const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+                cell_array[i][0]               = 2 * x[0] + 3 * x[1] + 2 * x[2] - 1;
+                cell_array[i][1]               = 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
+              });
 
-            function_symbol_id_list.push_back(
-              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            Table<const double> interpolate_array =
+              InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj,
+                                                                               cell_id_list);
+
+            REQUIRE(same_cell_value(cell_array, interpolate_array));
           }
+        }
+      }
 
+      SECTION("from -> (R)")
+      {
+        for (const auto& named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
           {
-            auto [i_symbol, found] = symbol_table->find("scalar_non_linear_3d", position);
-            REQUIRE(found);
-            REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+            auto mesh_3d = named_mesh.mesh();
 
-            function_symbol_id_list.push_back(
-              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
-          }
+            auto xj = MeshDataManager::instance().getMeshData(*mesh_3d).xj();
+
+            Array<CellId> cell_id_list{mesh_3d->numberOfCells() / 2};
+            for (size_t i_cell = 0; i_cell < cell_id_list.size(); ++i_cell) {
+              cell_id_list[i_cell] = static_cast<CellId>(2 * i_cell);
+            }
+
+            std::string_view data = R"(
+import math;
+let f_3d: R^3 -> (R), x -> (2 * x[0] + 3 * x[1] + 2 * x[2] - 1, 2 * exp(x[0]) * sin(x[1]) * x[2] + 3);
+)";
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+            auto ast = ASTBuilder::build(input);
 
-          Table<double> cell_array{cell_id_list.size(), 2};
-          parallel_for(
-            cell_id_list.size(), PUGS_LAMBDA(const size_t i) {
-              const TinyVector<Dimension>& x = xj[cell_id_list[i]];
-              cell_array[i][0]               = 2 * x[0] + 3 * x[1] + 2 * x[2] - 1;
-              cell_array[i][1]               = 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
-            });
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
 
-          Table<const double> interpolate_array =
-            InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj, cell_id_list);
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
 
-          REQUIRE(same_cell_value(cell_array, interpolate_array));
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
+
+            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+            TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+            position.byte = data.size();   // ensure that variables are declared at this point
+
+            std::vector<FunctionSymbolId> function_symbol_id_list;
+
+            {
+              auto [i_symbol, found] = symbol_table->find("f_3d", position);
+              REQUIRE(found);
+              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+              function_symbol_id_list.push_back(
+                FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            }
+
+            Table<double> cell_array{cell_id_list.size(), 2};
+            parallel_for(
+              cell_id_list.size(), PUGS_LAMBDA(const size_t i) {
+                const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+                cell_array[i][0]               = 2 * x[0] + 3 * x[1] + 2 * x[2] - 1;
+                cell_array[i][1]               = 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
+              });
+
+            Table<const double> interpolate_array =
+              InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj,
+                                                                               cell_id_list);
+
+            REQUIRE(same_cell_value(cell_array, interpolate_array));
+          }
         }
       }
     }