From abae271d222231b934d7e5a0191b8e58b10a0fe7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Thu, 22 Apr 2021 18:33:27 +0200
Subject: [PATCH] Add tests for InterpolateItemValue and proceed to a few
 clean-up

---
 src/language/utils/InterpolateItemValue.hpp |   13 +-
 tests/CMakeLists.txt                        |    1 +
 tests/test_InterpolateItemValue.cpp         | 1026 +++++++++++++++++++
 3 files changed, 1038 insertions(+), 2 deletions(-)
 create mode 100644 tests/test_InterpolateItemValue.cpp

diff --git a/src/language/utils/InterpolateItemValue.hpp b/src/language/utils/InterpolateItemValue.hpp
index 16a87b495..de12b563a 100644
--- a/src/language/utils/InterpolateItemValue.hpp
+++ b/src/language/utils/InterpolateItemValue.hpp
@@ -15,7 +15,7 @@ class InterpolateItemValue<OutputType(InputType)> : public PugsFunctionAdapter<O
 
  public:
   template <ItemType item_type>
-  static inline ItemValue<OutputType, item_type>
+  PUGS_INLINE static ItemValue<OutputType, item_type>
   interpolate(const FunctionSymbolId& function_symbol_id, const ItemValue<const InputType, item_type>& position)
   {
     auto& expression    = Adapter::getFunctionExpression(function_symbol_id);
@@ -46,7 +46,7 @@ class InterpolateItemValue<OutputType(InputType)> : public PugsFunctionAdapter<O
   }
 
   template <ItemType item_type>
-  static inline Array<OutputType>
+  PUGS_INLINE static Array<OutputType>
   interpolate(const FunctionSymbolId& function_symbol_id,
               const ItemValue<const InputType, item_type>& position,
               const Array<const ItemIdT<item_type>>& list_of_items)
@@ -77,6 +77,15 @@ class InterpolateItemValue<OutputType(InputType)> : public PugsFunctionAdapter<O
 
     return value;
   }
+
+  template <ItemType item_type>
+  PUGS_INLINE static Array<OutputType>
+  interpolate(const FunctionSymbolId& function_symbol_id,
+              const ItemValue<const InputType, item_type>& position,
+              const Array<ItemIdT<item_type>>& list_of_items)
+  {
+    return interpolate(function_symbol_id, position, Array<const ItemIdT<item_type>>{list_of_items});
+  }
 };
 
 #endif   // INTERPOLATE_ITEM_VALUE_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index b95035260..098e28a20 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -103,6 +103,7 @@ add_executable (mpi_unit_tests
   mpi_test_main.cpp
   test_DiscreteFunctionP0.cpp
   test_DiscreteFunctionP0Vector.cpp
+  test_InterpolateItemValue.cpp
   test_ItemArray.cpp
   test_ItemArrayUtils.cpp
   test_ItemValue.cpp
diff --git a/tests/test_InterpolateItemValue.cpp b/tests/test_InterpolateItemValue.cpp
new file mode 100644
index 000000000..fa786632c
--- /dev/null
+++ b/tests/test_InterpolateItemValue.cpp
@@ -0,0 +1,1026 @@
+#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 <language/utils/InterpolateItemValue.hpp>
+
+#include <pegtl/string_input.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("InterpolateItemValue", "[language]")
+{
+  SECTION("interpolate on all items")
+  {
+    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 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;
+let R3_affine_1d: R^1 -> R^3, x -> (2 * x[0] + 2, 3 * x[0], 2);
+let R3_non_linear_1d: R^1 -> R^3, x -> (2 * exp(x[0]) + 3, x[0] - 2, 3);
+let R2x2_affine_1d: R^1 -> R^2x2, x -> (2 * x[0] + 3 + 2, 3 * x[0], 2 * x[0], 2);
+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]);
+)";
+      TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+      auto ast = ASTBuilder::build(input);
+
+      ASTModulesImporter{*ast};
+      ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+      ASTSymbolTableBuilder{*ast};
+      ASTNodeDataTypeBuilder{*ast};
+
+      ASTNodeTypeCleaner<language::var_declaration>{*ast};
+      ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+      ASTNodeExpressionBuilder{*ast};
+
+      std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+      TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+      position.byte = data.size();   // ensure that variables are declared at this point
+
+      SECTION("scalar_affine_1d")
+      {
+        auto [i_symbol, found] = symbol_table->find("scalar_affine_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 * x[0] + 2;
+          });
+
+        CellValue<const double> interpolate_value =
+          InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("scalar_non_linear_1d")
+      {
+        auto [i_symbol, found] = symbol_table->find("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 * exp(x[0]) + 3;
+          });
+
+        CellValue<const double> interpolate_value =
+          InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R3_affine_1d")
+      {
+        auto [i_symbol, found] = symbol_table->find("R3_affine_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<TinyVector<3>> 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]            = TinyVector<3>{2 * x[0] + 2, 3 * x[0], 2};
+          });
+
+        CellValue<const TinyVector<3>> interpolate_value =
+          InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R3_non_linear_1d")
+      {
+        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<TinyVector<3>> 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]            = TinyVector<3>{2 * exp(x[0]) + 3, x[0] - 2, 3};
+          });
+
+        CellValue<const TinyVector<3>> interpolate_value =
+          InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R2x2_affine_1d")
+      {
+        auto [i_symbol, found] = symbol_table->find("R2x2_affine_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<TinyMatrix<2>> 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]            = TinyMatrix<2>{2 * x[0] + 3 + 2, 3 * x[0], 2 * x[0], 2};
+          });
+
+        CellValue<const TinyMatrix<2>> interpolate_value =
+          InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R2x2_non_linear_1d")
+      {
+        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<TinyMatrix<2>> 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] = TinyMatrix<2>{2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0]};
+          });
+
+        CellValue<const TinyMatrix<2>> interpolate_value =
+          InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+    }
+
+    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 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;
+let R3_affine_2d: R^2 -> R^3, x -> (2 * x[0] + 3 * x[1] + 2, 3 * x[0] + x[1], 2 * x[1]);
+let R3_non_linear_2d: R^2 -> R^3, x -> (2*exp(x[0])*sin(x[1])+3, x[0]-2*x[1], 3);
+let R2x2_affine_2d: R^2 -> R^2x2, x -> (2 * x[0] + 3 * x[1] + 2, 3 * x[0] + x[1], 2 * x[0] + x[1], 2);
+let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*x[1]), 3, x[0]*x[1]);
+)";
+      TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+      auto ast = ASTBuilder::build(input);
+
+      ASTModulesImporter{*ast};
+      ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+      ASTSymbolTableBuilder{*ast};
+      ASTNodeDataTypeBuilder{*ast};
+
+      ASTNodeTypeCleaner<language::var_declaration>{*ast};
+      ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+      ASTNodeExpressionBuilder{*ast};
+
+      std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+      TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+      position.byte = data.size();   // ensure that variables are declared at this point
+
+      SECTION("scalar_affine_2d")
+      {
+        auto [i_symbol, found] = symbol_table->find("scalar_affine_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 * x[0] + 3 * x[1] + 2;
+          });
+        CellValue<const double> interpolate_value =
+          InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("scalar_non_linear_2d")
+      {
+        auto [i_symbol, found] = symbol_table->find("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 * exp(x[0]) * sin(x[1]) + 3;
+          });
+        CellValue<const double> interpolate_value =
+          InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R3_affine_2d")
+      {
+        auto [i_symbol, found] = symbol_table->find("R3_affine_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<TinyVector<3>> 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]            = TinyVector<3>{2 * x[0] + 3 * x[1] + 2, 3 * x[0] + x[1], 2 * x[1]};
+          });
+        CellValue<const TinyVector<3>> interpolate_value =
+          InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R3_non_linear_2d")
+      {
+        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<TinyVector<3>> 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]            = TinyVector<3>{2 * exp(x[0]) * sin(x[1]) + 3, x[0] - 2 * x[1], 3};
+          });
+        CellValue<const TinyVector<3>> interpolate_value =
+          InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R2x2_affine_2d")
+      {
+        auto [i_symbol, found] = symbol_table->find("R2x2_affine_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<TinyMatrix<2>> 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] = TinyMatrix<2>{2 * x[0] + 3 * x[1] + 2, 3 * x[0] + x[1], 2 * x[0] + x[1], 2};
+          });
+        CellValue<const TinyMatrix<2>> interpolate_value =
+          InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R2x2_non_linear_2d")
+      {
+        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<TinyMatrix<2>> 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] = TinyMatrix<2>{2 * exp(x[0]) * sin(x[1]) + 3, sin(x[0] - 2 * x[1]), 3, x[0] * x[1]};
+          });
+        CellValue<const TinyMatrix<2>> interpolate_value =
+          InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+    }
+
+    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 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;
+let R3_affine_3d: R^3 -> R^3, x -> (2 * x[0] + 3 * x[1] + 2, 3 * x[0] + x[1] + 2 * x[2], 2 * x[2]);
+let R3_non_linear_3d: R^3 -> R^3, x -> (2 * exp(x[0]) * sin(x[1]) + x[2] + 3, x[0] * x[2] - 2 * x[1], 3);
+let R2x2_affine_3d: R^3 -> R^2x2, x -> (2 * x[0] + 3 * x[1] + 2 * x[2] + 1, 3 * x[0] + x[1] + 2 * x[2], 2 * x[0] + x[1] + x[2], 2);
+let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3, x[0] * x[1] * x[2]);
+)";
+      TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+      auto ast = ASTBuilder::build(input);
+
+      ASTModulesImporter{*ast};
+      ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+      ASTSymbolTableBuilder{*ast};
+      ASTNodeDataTypeBuilder{*ast};
+
+      ASTNodeTypeCleaner<language::var_declaration>{*ast};
+      ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+      ASTNodeExpressionBuilder{*ast};
+
+      std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+      TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+      position.byte = data.size();   // ensure that variables are declared at this point
+
+      SECTION("scalar_affine_3d")
+      {
+        auto [i_symbol, found] = symbol_table->find("scalar_affine_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 * x[0] + 3 * x[1] + 2 * x[2] - 1;
+          });
+        CellValue<const double> interpolate_value =
+          InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("scalar_non_linear_3d")
+      {
+        auto [i_symbol, found] = symbol_table->find("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 * exp(x[0]) * sin(x[1]) * x[2] + 3;
+          });
+        CellValue<const double> interpolate_value =
+          InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R3_affine_3d")
+      {
+        auto [i_symbol, found] = symbol_table->find("R3_affine_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<TinyVector<3>> 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] = TinyVector<3>{2 * x[0] + 3 * x[1] + 2, 3 * x[0] + x[1] + 2 * x[2], 2 * x[2]};
+          });
+        CellValue<const TinyVector<3>> interpolate_value =
+          InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R3_non_linear_3d")
+      {
+        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<TinyVector<3>> 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] = TinyVector<3>{2 * exp(x[0]) * sin(x[1]) + x[2] + 3, x[0] * x[2] - 2 * x[1], 3};
+          });
+        CellValue<const TinyVector<3>> interpolate_value =
+          InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R2x2_affine_3d")
+      {
+        auto [i_symbol, found] = symbol_table->find("R2x2_affine_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<TinyMatrix<2>> 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] =
+              TinyMatrix<2>{2 * x[0] + 3 * x[1] + 2 * x[2] + 1, 3 * x[0] + x[1] + 2 * x[2], 2 * x[0] + x[1] + x[2], 2};
+          });
+        CellValue<const TinyMatrix<2>> interpolate_value =
+          InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R2x2_non_linear_3d")
+      {
+        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<TinyMatrix<2>> 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] = TinyMatrix<2>{2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]),
+                                                3, x[0] * x[1] * x[2]};
+          });
+        CellValue<const TinyMatrix<2>> interpolate_value =
+          InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+    }
+  }
+
+  SECTION("interpolate on items list")
+  {
+    auto same_cell_value = [](auto interpolated, auto reference) -> bool {
+      for (size_t i = 0; i < interpolated.size(); ++i) {
+        if (interpolated[i] != reference[i]) {
+          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();
+
+      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 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;
+let R3_affine_1d: R^1 -> R^3, x -> (2 * x[0] + 2, 3 * x[0], 2);
+let R3_non_linear_1d: R^1 -> R^3, x -> (2 * exp(x[0]) + 3, x[0] - 2, 3);
+let R2x2_affine_1d: R^1 -> R^2x2, x -> (2 * x[0] + 3 + 2, 3 * x[0], 2 * x[0], 2);
+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]);
+)";
+      TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+      auto ast = ASTBuilder::build(input);
+
+      ASTModulesImporter{*ast};
+      ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+      ASTSymbolTableBuilder{*ast};
+      ASTNodeDataTypeBuilder{*ast};
+
+      ASTNodeTypeCleaner<language::var_declaration>{*ast};
+      ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+      ASTNodeExpressionBuilder{*ast};
+
+      std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+      TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+      position.byte = data.size();   // ensure that variables are declared at this point
+
+      SECTION("scalar_affine_1d")
+      {
+        auto [i_symbol, found] = symbol_table->find("scalar_affine_1d", position);
+        REQUIRE(found);
+        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+        Array<double> cell_value{cell_id_list.size()};
+        parallel_for(
+          cell_value.size(), PUGS_LAMBDA(const size_t i) {
+            const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+            cell_value[i]                  = 2 * x[0] + 2;
+          });
+
+        Array<const double> interpolate_value =
+          InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("scalar_non_linear_1d")
+      {
+        auto [i_symbol, found] = symbol_table->find("scalar_non_linear_1d", position);
+        REQUIRE(found);
+        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+        Array<double> cell_value{cell_id_list.size()};
+        parallel_for(
+          cell_value.size(), PUGS_LAMBDA(const size_t i) {
+            const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+            cell_value[i]                  = 2 * exp(x[0]) + 3;
+          });
+
+        Array<const double> interpolate_value =
+          InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R3_affine_1d")
+      {
+        auto [i_symbol, found] = symbol_table->find("R3_affine_1d", position);
+        REQUIRE(found);
+        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+        Array<TinyVector<3>> cell_value{cell_id_list.size()};
+        parallel_for(
+          cell_value.size(), PUGS_LAMBDA(const size_t i) {
+            const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+            cell_value[i]                  = TinyVector<3>{2 * x[0] + 2, 3 * x[0], 2};
+          });
+
+        Array<const TinyVector<3>> interpolate_value =
+          InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R3_non_linear_1d")
+      {
+        auto [i_symbol, found] = symbol_table->find("R3_non_linear_1d", position);
+        REQUIRE(found);
+        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+        Array<TinyVector<3>> cell_value{cell_id_list.size()};
+        parallel_for(
+          cell_value.size(), PUGS_LAMBDA(const size_t i) {
+            const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+            cell_value[i]                  = TinyVector<3>{2 * exp(x[0]) + 3, x[0] - 2, 3};
+          });
+
+        Array<const TinyVector<3>> interpolate_value =
+          InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R2x2_affine_1d")
+      {
+        auto [i_symbol, found] = symbol_table->find("R2x2_affine_1d", position);
+        REQUIRE(found);
+        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+        Array<TinyMatrix<2>> cell_value{cell_id_list.size()};
+        parallel_for(
+          cell_value.size(), PUGS_LAMBDA(const size_t i) {
+            const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+            cell_value[i]                  = TinyMatrix<2>{2 * x[0] + 3 + 2, 3 * x[0], 2 * x[0], 2};
+          });
+
+        Array<const TinyMatrix<2>> interpolate_value =
+          InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R2x2_non_linear_1d")
+      {
+        auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_1d", position);
+        REQUIRE(found);
+        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+        Array<TinyMatrix<2>> cell_value{cell_id_list.size()};
+        parallel_for(
+          cell_value.size(), PUGS_LAMBDA(const size_t i) {
+            const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+            cell_value[i] = TinyMatrix<2>{2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0]};
+          });
+
+        Array<const TinyMatrix<2>> interpolate_value =
+          InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+    }
+
+    SECTION("2D")
+    {
+      constexpr size_t Dimension = 2;
+
+      const auto& mesh_2d = MeshDataBaseForTests::get().cartesianMesh2D();
+      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 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;
+let R3_affine_2d: R^2 -> R^3, x -> (2 * x[0] + 3 * x[1] + 2, 3 * x[0] + x[1], 2 * x[1]);
+let R3_non_linear_2d: R^2 -> R^3, x -> (2*exp(x[0])*sin(x[1])+3, x[0]-2*x[1], 3);
+let R2x2_affine_2d: R^2 -> R^2x2, x -> (2 * x[0] + 3 * x[1] + 2, 3 * x[0] + x[1], 2 * x[0] + x[1], 2);
+let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2*exp(x[0])*sin(x[1])+3, sin(x[0]-2*x[1]), 3, x[0]*x[1]);
+)";
+      TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+      auto ast = ASTBuilder::build(input);
+
+      ASTModulesImporter{*ast};
+      ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+      ASTSymbolTableBuilder{*ast};
+      ASTNodeDataTypeBuilder{*ast};
+
+      ASTNodeTypeCleaner<language::var_declaration>{*ast};
+      ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+      ASTNodeExpressionBuilder{*ast};
+
+      std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+      TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+      position.byte = data.size();   // ensure that variables are declared at this point
+
+      SECTION("scalar_affine_2d")
+      {
+        auto [i_symbol, found] = symbol_table->find("scalar_affine_2d", position);
+        REQUIRE(found);
+        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+        Array<double> cell_value{cell_id_list.size()};
+        parallel_for(
+          cell_value.size(), PUGS_LAMBDA(const size_t i) {
+            const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+            cell_value[i]                  = 2 * x[0] + 3 * x[1] + 2;
+          });
+
+        Array<const double> interpolate_value =
+          InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("scalar_non_linear_2d")
+      {
+        auto [i_symbol, found] = symbol_table->find("scalar_non_linear_2d", position);
+        REQUIRE(found);
+        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+        Array<double> cell_value{cell_id_list.size()};
+        parallel_for(
+          cell_value.size(), PUGS_LAMBDA(const size_t i) {
+            const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+            cell_value[i]                  = 2 * exp(x[0]) * sin(x[1]) + 3;
+          });
+        Array<const double> interpolate_value =
+          InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R3_affine_2d")
+      {
+        auto [i_symbol, found] = symbol_table->find("R3_affine_2d", position);
+        REQUIRE(found);
+        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+        Array<TinyVector<3>> cell_value{cell_id_list.size()};
+        parallel_for(
+          cell_value.size(), PUGS_LAMBDA(const size_t i) {
+            const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+            cell_value[i]                  = TinyVector<3>{2 * x[0] + 3 * x[1] + 2, 3 * x[0] + x[1], 2 * x[1]};
+          });
+        Array<const TinyVector<3>> interpolate_value =
+          InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R3_non_linear_2d")
+      {
+        auto [i_symbol, found] = symbol_table->find("R3_non_linear_2d", position);
+        REQUIRE(found);
+        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+        Array<TinyVector<3>> cell_value{cell_id_list.size()};
+        parallel_for(
+          cell_value.size(), PUGS_LAMBDA(const size_t i) {
+            const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+            cell_value[i]                  = TinyVector<3>{2 * exp(x[0]) * sin(x[1]) + 3, x[0] - 2 * x[1], 3};
+          });
+        Array<const TinyVector<3>> interpolate_value =
+          InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R2x2_affine_2d")
+      {
+        auto [i_symbol, found] = symbol_table->find("R2x2_affine_2d", position);
+        REQUIRE(found);
+        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+        Array<TinyMatrix<2>> cell_value{cell_id_list.size()};
+        parallel_for(
+          cell_value.size(), PUGS_LAMBDA(const size_t i) {
+            const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+            cell_value[i] = TinyMatrix<2>{2 * x[0] + 3 * x[1] + 2, 3 * x[0] + x[1], 2 * x[0] + x[1], 2};
+          });
+        Array<const TinyMatrix<2>> interpolate_value =
+          InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R2x2_non_linear_2d")
+      {
+        auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_2d", position);
+        REQUIRE(found);
+        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+        Array<TinyMatrix<2>> cell_value{cell_id_list.size()};
+        parallel_for(
+          cell_value.size(), PUGS_LAMBDA(const size_t i) {
+            const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+            cell_value[i] = TinyMatrix<2>{2 * exp(x[0]) * sin(x[1]) + 3, sin(x[0] - 2 * x[1]), 3, x[0] * x[1]};
+          });
+        Array<const TinyMatrix<2>> interpolate_value =
+          InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+    }
+
+    SECTION("3D")
+    {
+      constexpr size_t Dimension = 3;
+
+      const auto& mesh_3d = MeshDataBaseForTests::get().cartesianMesh3D();
+      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 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;
+let R3_affine_3d: R^3 -> R^3, x -> (2 * x[0] + 3 * x[1] + 2, 3 * x[0] + x[1] + 2 * x[2], 2 * x[2]);
+let R3_non_linear_3d: R^3 -> R^3, x -> (2 * exp(x[0]) * sin(x[1]) + x[2] + 3, x[0] * x[2] - 2 * x[1], 3);
+let R2x2_affine_3d: R^3 -> R^2x2, x -> (2 * x[0] + 3 * x[1] + 2 * x[2] + 1, 3 * x[0] + x[1] + 2 * x[2], 2 * x[0] + x[1] + x[2], 2);
+let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3, x[0] * x[1] * x[2]);
+)";
+      TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+      auto ast = ASTBuilder::build(input);
+
+      ASTModulesImporter{*ast};
+      ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+      ASTSymbolTableBuilder{*ast};
+      ASTNodeDataTypeBuilder{*ast};
+
+      ASTNodeTypeCleaner<language::var_declaration>{*ast};
+      ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+      ASTNodeExpressionBuilder{*ast};
+
+      std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+      TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+      position.byte = data.size();   // ensure that variables are declared at this point
+
+      SECTION("scalar_affine_3d")
+      {
+        auto [i_symbol, found] = symbol_table->find("scalar_affine_3d", position);
+        REQUIRE(found);
+        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+        Array<double> cell_value{cell_id_list.size()};
+        parallel_for(
+          cell_value.size(), PUGS_LAMBDA(const size_t i) {
+            const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+            cell_value[i]                  = 2 * x[0] + 3 * x[1] + 2 * x[2] - 1;
+          });
+        Array<const double> interpolate_value =
+          InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("scalar_non_linear_3d")
+      {
+        auto [i_symbol, found] = symbol_table->find("scalar_non_linear_3d", position);
+        REQUIRE(found);
+        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+        Array<double> cell_value{cell_id_list.size()};
+        parallel_for(
+          cell_value.size(), PUGS_LAMBDA(const size_t i) {
+            const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+            cell_value[i]                  = 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
+          });
+        Array<const double> interpolate_value =
+          InterpolateItemValue<double(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R3_affine_3d")
+      {
+        auto [i_symbol, found] = symbol_table->find("R3_affine_3d", position);
+        REQUIRE(found);
+        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+        Array<TinyVector<3>> cell_value{cell_id_list.size()};
+        parallel_for(
+          cell_value.size(), PUGS_LAMBDA(const size_t i) {
+            const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+            cell_value[i] = TinyVector<3>{2 * x[0] + 3 * x[1] + 2, 3 * x[0] + x[1] + 2 * x[2], 2 * x[2]};
+          });
+        Array<const TinyVector<3>> interpolate_value =
+          InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R3_non_linear_3d")
+      {
+        auto [i_symbol, found] = symbol_table->find("R3_non_linear_3d", position);
+        REQUIRE(found);
+        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+        Array<TinyVector<3>> cell_value{cell_id_list.size()};
+        parallel_for(
+          cell_value.size(), PUGS_LAMBDA(const size_t i) {
+            const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+            cell_value[i] = TinyVector<3>{2 * exp(x[0]) * sin(x[1]) + x[2] + 3, x[0] * x[2] - 2 * x[1], 3};
+          });
+        Array<const TinyVector<3>> interpolate_value =
+          InterpolateItemValue<TinyVector<3>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R2x2_affine_3d")
+      {
+        auto [i_symbol, found] = symbol_table->find("R2x2_affine_3d", position);
+        REQUIRE(found);
+        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+        Array<TinyMatrix<2>> cell_value{cell_id_list.size()};
+        parallel_for(
+          cell_value.size(), PUGS_LAMBDA(const size_t i) {
+            const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+            cell_value[i] =
+              TinyMatrix<2>{2 * x[0] + 3 * x[1] + 2 * x[2] + 1, 3 * x[0] + x[1] + 2 * x[2], 2 * x[0] + x[1] + x[2], 2};
+          });
+        Array<const TinyMatrix<2>> interpolate_value =
+          InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+
+      SECTION("R2x2_non_linear_3d")
+      {
+        auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_3d", position);
+        REQUIRE(found);
+        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+        Array<TinyMatrix<2>> cell_value{cell_id_list.size()};
+        parallel_for(
+          cell_value.size(), PUGS_LAMBDA(const size_t i) {
+            const TinyVector<Dimension>& x = xj[cell_id_list[i]];
+            cell_value[i] = TinyMatrix<2>{2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3,
+                                          x[0] * x[1] * x[2]};
+          });
+        Array<const TinyMatrix<2>> interpolate_value =
+          InterpolateItemValue<TinyMatrix<2>(TinyVector<Dimension>)>::interpolate(function_symbol_id, xj, cell_id_list);
+
+        REQUIRE(same_cell_value(cell_value, interpolate_value));
+      }
+    }
+  }
+}
-- 
GitLab