diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 9237f5db7f60d323e2bad3127628952dc43a8919..ef2669559feb7e81180de8f0f74d8395ee41d79e 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -38,8 +38,11 @@
 #include <scheme/HypoelasticSolver.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/InflowBoundaryConditionDescriptor.hpp>
 #include <scheme/NeumannBoundaryConditionDescriptor.hpp>
+#include <scheme/OutflowBoundaryConditionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <scheme/VariableBCDescriptor.hpp>
 #include <utils/Socket.hpp>
 
 #include <memory>
@@ -51,6 +54,7 @@ SchemeModule::SchemeModule()
   this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IQuadratureDescriptor>>);
 
   this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IBoundaryConditionDescriptor>>);
+  this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const VariableBCDescriptor>>);
 
   this->_addBuiltinFunction("P0", std::function(
 
@@ -76,6 +80,16 @@ SchemeModule::SchemeModule()
 
                                        ));
 
+  this->_addBuiltinFunction("variable_bc", std::function(
+
+                                             [](const std::shared_ptr<const DiscreteFunctionVariant>& v,
+                                                const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                                  bc_descriptor_list) -> std::shared_ptr<const VariableBCDescriptor> {
+                                               return std::make_shared<VariableBCDescriptor>(v, bc_descriptor_list);
+                                             }
+
+                                             ));
+
   this->_addBuiltinFunction("GaussLobatto", std::function(
 
                                               [](uint64_t degree) -> std::shared_ptr<const IQuadratureDescriptor> {
@@ -368,6 +382,26 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this->_addBuiltinFunction("inflow", std::function(
+
+                                        [](std::shared_ptr<const IBoundaryDescriptor> boundary,
+                                           const FunctionSymbolId& function_id)
+                                          -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+                                          return std::make_shared<InflowBoundaryConditionDescriptor>(boundary,
+                                                                                                     function_id);
+                                        }
+
+                                        ));
+
+  this->_addBuiltinFunction("outflow", std::function(
+
+                                         [](std::shared_ptr<const IBoundaryDescriptor> boundary)
+                                           -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+                                           return std::make_shared<OutflowBoundaryConditionDescriptor>(boundary);
+                                         }
+
+                                         ));
+
   this->_addBuiltinFunction("glace_fluxes", std::function(
 
                                               [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
@@ -765,6 +799,7 @@ SchemeModule::SchemeModule()
                                                  }
 
                                                  ));
+
   this->_addBuiltinFunction("hypoelastic_dt", std::function(
 
                                                 [](const std::shared_ptr<const DiscreteFunctionVariant>& c) -> double {
@@ -773,16 +808,16 @@ SchemeModule::SchemeModule()
 
                                                 ));
 
-  this->_addBuiltinFunction("fluxing_advection",
-                            std::function(
+  this->_addBuiltinFunction("fluxing_advection", std::function(
 
-                              [](const std::shared_ptr<const IMesh> new_mesh,
-                                 const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& remapped_variables)
-                                -> std::vector<std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return advectByFluxing(new_mesh, remapped_variables);
-                              }
+                                                   [](const std::shared_ptr<const IMesh> new_mesh,
+                                                      const std::vector<std::shared_ptr<const VariableBCDescriptor>>&
+                                                        remapped_quantity_with_bc)
+                                                     -> std::vector<std::shared_ptr<const DiscreteFunctionVariant>> {
+                                                     return advectByFluxing(new_mesh, remapped_quantity_with_bc);
+                                                   }
 
-                              ));
+                                                   ));
 
   MathFunctionRegisterForVh{*this};
 }
diff --git a/src/language/modules/SchemeModule.hpp b/src/language/modules/SchemeModule.hpp
index d56a5ec74069bd23d8169f25ff8a829c7c4a53ff..68cffcf18dc9438cfa5607bdf70739bf223a175d 100644
--- a/src/language/modules/SchemeModule.hpp
+++ b/src/language/modules/SchemeModule.hpp
@@ -10,6 +10,11 @@ template <>
 inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const IBoundaryConditionDescriptor>> =
   ASTNodeDataType::build<ASTNodeDataType::type_id_t>("boundary_condition");
 
+class VariableBCDescriptor;
+template <>
+inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const VariableBCDescriptor>> =
+  ASTNodeDataType::build<ASTNodeDataType::type_id_t>("variable_boundary_condition");
+
 class DiscreteFunctionVariant;
 template <>
 inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const DiscreteFunctionVariant>> =
diff --git a/src/language/utils/EvaluateArrayAtPoints.hpp b/src/language/utils/EvaluateArrayAtPoints.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..558cdff30cc355287cd2c1233884f765b1cac1b7
--- /dev/null
+++ b/src/language/utils/EvaluateArrayAtPoints.hpp
@@ -0,0 +1,72 @@
+#ifndef EVALUATE_ARRAY_AT_POINTS_HPP
+#define EVALUATE_ARRAY_AT_POINTS_HPP
+
+#include <language/utils/PugsFunctionAdapter.hpp>
+#include <utils/Array.hpp>
+#include <utils/Table.hpp>
+
+class FunctionSymbolId;
+
+template <typename T>
+class EvaluateArrayAtPoints;
+template <typename OutputType, typename InputType>
+class EvaluateArrayAtPoints<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)>
+{
+  using Adapter = PugsFunctionAdapter<OutputType(InputType)>;
+
+ public:
+  template <typename InputArrayT, typename OutputTableT>
+  static PUGS_INLINE void
+  evaluateTo(const FunctionSymbolId& function_symbol_id, const InputArrayT& position, OutputTableT& table)
+  {
+    static_assert(std::is_same_v<std::remove_const_t<typename InputArrayT::data_type>, InputType>,
+                  "invalid input data type");
+    static_assert(std::is_same_v<std::remove_const_t<typename OutputTableT::data_type>, OutputType>,
+                  "invalid output data type");
+
+    auto& expression    = Adapter::getFunctionExpression(function_symbol_id);
+    auto convert_result = Adapter::getArrayResultConverter(expression.m_data_type);
+
+    auto context_list = Adapter::getContextList(expression);
+
+    using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
+    Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
+
+    if constexpr (std::is_arithmetic_v<OutputType>) {
+      table.fill(0);
+    } else if constexpr (is_tiny_vector_v<OutputType> or is_tiny_matrix_v<OutputType>) {
+      table.fill(zero);
+    } else {
+      static_assert(std::is_same_v<OutputType, double>, "unexpected output type");
+    }
+
+    parallel_for(size(position), [=, &expression, &tokens](typename InputArrayT::index_type i) {
+      const int32_t t = tokens.acquire();
+
+      auto& execution_policy = context_list[t];
+
+      Adapter::convertArgs(execution_policy.currentContext(), position[i]);
+      auto result  = expression.execute(execution_policy);
+      auto&& array = convert_result(std::move(result));
+
+      for (size_t j = 0; j < array.size(); ++j) {
+        table[i][j] = array[j];
+      }
+
+      tokens.release(t);
+    });
+  }
+
+  template <class InputArrayT>
+  static PUGS_INLINE Table<OutputType>
+  evaluate(const FunctionSymbolId& function_symbol_id, const InputArrayT& position)
+  {
+    static_assert(std::is_same_v<std::remove_const_t<typename InputArrayT::data_type>, InputType>,
+                  "invalid input data type");
+    Table<OutputType> value(size(position));
+    evaluateArrayTo(function_symbol_id, position, value);
+    return value;
+  }
+};
+
+#endif   // EVALUATE_ARRAY_AT_POINTS_HPP
diff --git a/src/language/utils/InterpolateItemArray.hpp b/src/language/utils/InterpolateItemArray.hpp
index 17b72de6d62e0ac68062343010154257c43d3006..fde135d6bcbb4c8995395660226e8d3435c599b4 100644
--- a/src/language/utils/InterpolateItemArray.hpp
+++ b/src/language/utils/InterpolateItemArray.hpp
@@ -1,6 +1,7 @@
 #ifndef INTERPOLATE_ITEM_ARRAY_HPP
 #define INTERPOLATE_ITEM_ARRAY_HPP
 
+#include <language/utils/EvaluateArrayAtPoints.hpp>
 #include <language/utils/InterpolateItemValue.hpp>
 #include <mesh/ItemArray.hpp>
 #include <mesh/ItemType.hpp>
@@ -12,24 +13,47 @@ class InterpolateItemArray<OutputType(InputType)>
 {
   static constexpr size_t Dimension = OutputType::Dimension;
 
+ private:
+  PUGS_INLINE static bool
+  _isSingleTupleFunction(const std::vector<FunctionSymbolId>& function_symbol_id_list)
+  {
+    if (function_symbol_id_list.size() > 1) {
+      return false;
+    } else {
+      Assert(function_symbol_id_list.size() == 1);
+      const FunctionSymbolId& function_symbol_id = function_symbol_id_list[0];
+      return (function_symbol_id.descriptor().domainMappingNode().children[1]->m_data_type == ASTNodeDataType::tuple_t);
+    }
+  }
+
  public:
   template <ItemType item_type>
   PUGS_INLINE static ItemArray<OutputType, item_type>
   interpolate(const std::vector<FunctionSymbolId>& function_symbol_id_list,
               const ItemValue<const InputType, item_type>& position)
   {
-    ItemArray<OutputType, item_type> item_array{*position.connectivity_ptr(), function_symbol_id_list.size()};
+    if (_isSingleTupleFunction(function_symbol_id_list)) {
+      const FunctionSymbolId& function_symbol_id = function_symbol_id_list[0];
+      const size_t table_size = function_symbol_id.descriptor().definitionNode().children[1]->children.size();
 
-    for (size_t i_function_symbol = 0; i_function_symbol < function_symbol_id_list.size(); ++i_function_symbol) {
-      const FunctionSymbolId& function_symbol_id = function_symbol_id_list[i_function_symbol];
-      ItemValue<OutputType, item_type> item_value =
-        InterpolateItemValue<OutputType(InputType)>::interpolate(function_symbol_id, position);
-      parallel_for(
-        item_value.numberOfItems(),
-        PUGS_LAMBDA(ItemIdT<item_type> item_id) { item_array[item_id][i_function_symbol] = item_value[item_id]; });
-    }
+      ItemArray<OutputType, item_type> item_array{*position.connectivity_ptr(), table_size};
+      EvaluateArrayAtPoints<OutputType(const InputType)>::evaluateTo(function_symbol_id, position, item_array);
+
+      return item_array;
+    } else {
+      ItemArray<OutputType, item_type> item_array{*position.connectivity_ptr(), function_symbol_id_list.size()};
+
+      for (size_t i_function_symbol = 0; i_function_symbol < function_symbol_id_list.size(); ++i_function_symbol) {
+        const FunctionSymbolId& function_symbol_id = function_symbol_id_list[i_function_symbol];
+        ItemValue<OutputType, item_type> item_value =
+          InterpolateItemValue<OutputType(InputType)>::interpolate(function_symbol_id, position);
+        parallel_for(
+          item_value.numberOfItems(),
+          PUGS_LAMBDA(ItemIdT<item_type> item_id) { item_array[item_id][i_function_symbol] = item_value[item_id]; });
+      }
 
-    return item_array;
+      return item_array;
+    }
   }
 
   template <ItemType item_type>
@@ -38,18 +62,36 @@ class InterpolateItemArray<OutputType(InputType)>
               const ItemValue<const InputType, item_type>& position,
               const Array<const ItemIdT<item_type>>& list_of_items)
   {
-    Table<OutputType> table{list_of_items.size(), function_symbol_id_list.size()};
+    if (_isSingleTupleFunction(function_symbol_id_list)) {
+      Array<InputType> item_position{list_of_items.size()};
+      using ItemId = ItemIdT<item_type>;
+      parallel_for(
+        list_of_items.size(), PUGS_LAMBDA(size_t i_item) {
+          ItemId item_id        = list_of_items[i_item];
+          item_position[i_item] = position[item_id];
+        });
 
-    for (size_t i_function_symbol = 0; i_function_symbol < function_symbol_id_list.size(); ++i_function_symbol) {
-      const FunctionSymbolId& function_symbol_id = function_symbol_id_list[i_function_symbol];
-      Array<OutputType> array =
-        InterpolateItemValue<OutputType(InputType)>::interpolate(function_symbol_id, position, list_of_items);
+      const FunctionSymbolId& function_symbol_id = function_symbol_id_list[0];
+      const size_t table_size = function_symbol_id.descriptor().definitionNode().children[1]->children.size();
 
-      parallel_for(
-        array.size(), PUGS_LAMBDA(size_t i) { table[i][i_function_symbol] = array[i]; });
-    }
+      Table<OutputType> table{list_of_items.size(), table_size};
+      EvaluateArrayAtPoints<OutputType(const InputType)>::evaluateTo(function_symbol_id, item_position, table);
 
-    return table;
+      return table;
+    } else {
+      Table<OutputType> table{list_of_items.size(), function_symbol_id_list.size()};
+
+      for (size_t i_function_symbol = 0; i_function_symbol < function_symbol_id_list.size(); ++i_function_symbol) {
+        const FunctionSymbolId& function_symbol_id = function_symbol_id_list[i_function_symbol];
+        Array<OutputType> array =
+          InterpolateItemValue<OutputType(InputType)>::interpolate(function_symbol_id, position, list_of_items);
+
+        parallel_for(
+          array.size(), PUGS_LAMBDA(size_t i) { table[i][i_function_symbol] = array[i]; });
+      }
+
+      return table;
+    }
   }
 
   template <ItemType item_type>
diff --git a/src/language/utils/ItemArrayVariantFunctionInterpoler.cpp b/src/language/utils/ItemArrayVariantFunctionInterpoler.cpp
index 02d47d399ed9ec0fdb570fc3877b0b384f7a5429..48b51efba5aa8a9a9905247629e9e398a70ab17a 100644
--- a/src/language/utils/ItemArrayVariantFunctionInterpoler.cpp
+++ b/src/language/utils/ItemArrayVariantFunctionInterpoler.cpp
@@ -55,7 +55,8 @@ ItemArrayVariantFunctionInterpoler::_interpolate() const
 {
   const ASTNodeDataType data_type = [&] {
     const auto& function0_descriptor = m_function_id_list[0].descriptor();
-    Assert(function0_descriptor.domainMappingNode().children[1]->m_data_type == ASTNodeDataType::typename_t);
+    Assert(function0_descriptor.domainMappingNode().children[1]->m_data_type == ASTNodeDataType::typename_t or
+           function0_descriptor.domainMappingNode().children[1]->m_data_type == ASTNodeDataType::tuple_t);
 
     ASTNodeDataType data_type = function0_descriptor.domainMappingNode().children[1]->m_data_type.contentType();
 
diff --git a/src/language/utils/PugsFunctionAdapter.hpp b/src/language/utils/PugsFunctionAdapter.hpp
index de5e884bf63bc0ac3726864c4c4571c0a0fb9f55..a4d1f76bd45285294008bfbaa24b6172a2b0b774 100644
--- a/src/language/utils/PugsFunctionAdapter.hpp
+++ b/src/language/utils/PugsFunctionAdapter.hpp
@@ -288,6 +288,92 @@ class PugsFunctionAdapter<OutputType(InputType...)>
     }
   }
 
+  [[nodiscard]] PUGS_INLINE static std::function<std::vector<OutputType>(DataVariant&& result)>
+  getArrayResultConverter(const ASTNodeDataType& data_type)
+  {
+    Assert(data_type == ASTNodeDataType::list_t);
+
+    if constexpr (std::is_arithmetic_v<OutputType>) {
+      return [&](DataVariant&& result) -> std::vector<OutputType> {
+        return std::visit(
+          [&](auto&& value) -> std::vector<OutputType> {
+            using ValueType = std::decay_t<decltype(value)>;
+            if constexpr (std::is_same_v<ValueType, AggregateDataVariant>) {
+              std::vector<OutputType> array(value.size());
+
+              for (size_t i = 0; i < value.size(); ++i) {
+                array[i] = std::visit(
+                  [&](auto&& value_i) -> OutputType {
+                    using Value_I_Type = std::decay_t<decltype(value_i)>;
+                    if constexpr (std::is_arithmetic_v<Value_I_Type>) {
+                      return value_i;
+                    } else {
+                      // LCOV_EXCL_START
+                      throw UnexpectedError("expecting arithmetic type");
+                      // LCOV_EXCL_STOP
+                    }
+                  },
+                  value[i]);
+              }
+
+              return array;
+            } else {
+              // LCOV_EXCL_START
+              throw UnexpectedError("invalid DataVariant");
+              // LCOV_EXCL_STOP
+            }
+          },
+          result);
+      };
+    } else if constexpr (is_tiny_vector_v<OutputType> or (is_tiny_matrix_v<OutputType>)) {
+      return [&](DataVariant&& result) -> std::vector<OutputType> {
+        return std::visit(
+          [&](auto&& value) -> std::vector<OutputType> {
+            using ValueType = std::decay_t<decltype(value)>;
+            if constexpr (std::is_same_v<ValueType, AggregateDataVariant>) {
+              std::vector<OutputType> array(value.size());
+
+              for (size_t i = 0; i < value.size(); ++i) {
+                array[i] = std::visit(
+                  [&](auto&& value_i) -> OutputType {
+                    using Value_I_Type = std::decay_t<decltype(value_i)>;
+                    if constexpr (std::is_same_v<Value_I_Type, OutputType>) {
+                      return value_i;
+                    } else if constexpr (OutputType::Dimension == 1) {
+                      if constexpr (std::is_arithmetic_v<Value_I_Type>) {
+                        return OutputType(value_i);
+                      } else {
+                        // LCOV_EXCL_START
+                        throw UnexpectedError("expecting arithmetic type");
+                        // LCOV_EXCL_STOP
+                      }
+                    } else if constexpr (std::is_integral_v<Value_I_Type>) {
+                      // reaching this point, it must be a null vector
+                      // or a null matrix
+                      return OutputType{zero};
+                    } else {
+                      // LCOV_EXCL_START
+                      throw UnexpectedError("expecting arithmetic type");
+                      // LCOV_EXCL_STOP
+                    }
+                  },
+                  value[i]);
+              }
+
+              return array;
+            } else {
+              // LCOV_EXCL_START
+              throw UnexpectedError("invalid DataVariant");
+              // LCOV_EXCL_STOP
+            }
+          },
+          result);
+      };
+    } else {
+      throw UnexpectedError("non-arithmetic tuple type");
+    }
+  }
+
   PugsFunctionAdapter() = delete;
 };
 
diff --git a/src/scheme/DiscreteFunctionVectorInterpoler.cpp b/src/scheme/DiscreteFunctionVectorInterpoler.cpp
index b416400c65a719614be50a53f167a44d859c9497..57724b7722df4a8410993afee713c35f8d1c73e1 100644
--- a/src/scheme/DiscreteFunctionVectorInterpoler.cpp
+++ b/src/scheme/DiscreteFunctionVectorInterpoler.cpp
@@ -96,11 +96,7 @@ template <size_t Dimension>
 DiscreteFunctionVariant
 DiscreteFunctionVectorInterpoler::_interpolate() const
 {
-  for (const auto& function_id : m_function_id_list) {
-    const auto& function_descriptor = function_id.descriptor();
-    Assert(function_descriptor.domainMappingNode().children[1]->m_data_type == ASTNodeDataType::typename_t);
-    const ASTNodeDataType& data_type = function_descriptor.domainMappingNode().children[1]->m_data_type.contentType();
-
+  auto check_data_type = [](const ASTNodeDataType& data_type) {
     switch (data_type) {
     case ASTNodeDataType::bool_t:
     case ASTNodeDataType::unsigned_int_t:
@@ -115,7 +111,32 @@ DiscreteFunctionVectorInterpoler::_interpolate() const
       throw NormalError(os.str());
     }
     }
+  };
+
+  if (m_function_id_list.size() == 1) {
+    const auto& function_descriptor = m_function_id_list[0].descriptor();
+    if ((function_descriptor.domainMappingNode().children[1]->m_data_type == ASTNodeDataType::typename_t) or
+        (function_descriptor.domainMappingNode().children[1]->m_data_type == ASTNodeDataType::tuple_t)) {
+      const ASTNodeDataType& data_type = function_descriptor.domainMappingNode().children[1]->m_data_type.contentType();
+
+      check_data_type(data_type);
+    } else {
+      throw UnexpectedError("incorrect function");
+    }
+  } else {
+    for (const auto& function_id : m_function_id_list) {
+      const auto& function_descriptor = function_id.descriptor();
+      if (function_descriptor.domainMappingNode().children[1]->m_data_type == ASTNodeDataType::typename_t) {
+        const ASTNodeDataType& data_type =
+          function_descriptor.domainMappingNode().children[1]->m_data_type.contentType();
+
+        check_data_type(data_type);
+      } else {
+        throw UnexpectedError("incorrect function");
+      }
+    }
   }
+
   return this->_interpolate<Dimension, double>();
 }
 
diff --git a/src/scheme/FluxingAdvectionSolver.cpp b/src/scheme/FluxingAdvectionSolver.cpp
index 57f0d4842c82c5bcea68cb83de687d82a5b1bca7..264ecef042956e6b77230e00f49cd321b50dcc6a 100644
--- a/src/scheme/FluxingAdvectionSolver.cpp
+++ b/src/scheme/FluxingAdvectionSolver.cpp
@@ -1,6 +1,8 @@
 #include <scheme/FluxingAdvectionSolver.hpp>
 
 #include <language/utils/EvaluateAtPoints.hpp>
+#include <language/utils/InterpolateItemArray.hpp>
+#include <language/utils/InterpolateItemValue.hpp>
 #include <mesh/Connectivity.hpp>
 #include <mesh/IMesh.hpp>
 #include <mesh/ItemArrayUtils.hpp>
@@ -9,10 +11,18 @@
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/MeshFaceBoundary.hpp>
+#include <mesh/MeshFlatFaceBoundary.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
 #include <mesh/SubItemValuePerItem.hpp>
 #include <scheme/DiscreteFunctionP0.hpp>
 #include <scheme/DiscreteFunctionUtils.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/InflowBoundaryConditionDescriptor.hpp>
+#include <scheme/OutflowBoundaryConditionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+
+#include <variant>
+#include <vector>
 
 template <size_t Dimension>
 class FluxingAdvectionSolver
@@ -36,7 +46,25 @@ class FluxingAdvectionSolver
 
                                     CellArray<double>>;
 
+  template <typename DataType>
+  class InflowValueBoundaryCondition;
+  class InflowArrayBoundaryCondition;
+  class OutflowBoundaryCondition;
+  class SymmetryBoundaryCondition;
+
+  using BoundaryConditionVariant = std::variant<InflowValueBoundaryCondition<double>,   //
+                                                InflowValueBoundaryCondition<TinyVector<1>>,
+                                                InflowValueBoundaryCondition<TinyVector<2>>,
+                                                InflowValueBoundaryCondition<TinyVector<3>>,
+                                                InflowValueBoundaryCondition<TinyMatrix<1>>,
+                                                InflowValueBoundaryCondition<TinyMatrix<2>>,
+                                                InflowValueBoundaryCondition<TinyMatrix<3>>,
+                                                InflowArrayBoundaryCondition,
+                                                OutflowBoundaryCondition,
+                                                SymmetryBoundaryCondition>;
+
   std::vector<RemapVariant> m_remapped_list;
+  std::vector<std::vector<BoundaryConditionVariant>> m_boundary_condition_list;
 
   FaceValue<const CellId> m_donnor_cell;
   FaceValue<const double> m_cycle_fluxing_volume;
@@ -55,21 +83,111 @@ class FluxingAdvectionSolver
     m_remapped_list.emplace_back(copy(old_q.cellValues()));
   }
 
-  template <typename DataType>
   void
-  _storeValues(const DiscreteFunctionP0Vector<Dimension, const DataType>& old_q)
+  _storeValues(const DiscreteFunctionP0Vector<Dimension, const double>& old_q)
   {
     m_remapped_list.emplace_back(copy(old_q.cellArrays()));
   }
 
+  template <typename DataType>
+  void
+  _storeValueBCList(const size_t i_quantity,
+                    const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_list)
+  {
+    const Mesh<Connectivity<Dimension>>& mesh = *m_old_mesh;
+
+    for (auto& i_bc : bc_list) {
+      switch (i_bc->type()) {
+      case IBoundaryConditionDescriptor::Type::symmetry: {
+        m_boundary_condition_list[i_quantity].emplace_back(
+          SymmetryBoundaryCondition(getMeshFlatFaceBoundary(mesh, i_bc->boundaryDescriptor())));
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::outflow: {
+        m_boundary_condition_list[i_quantity].emplace_back(
+          OutflowBoundaryCondition(getMeshFaceBoundary(mesh, i_bc->boundaryDescriptor())));
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::inflow: {
+        const InflowBoundaryConditionDescriptor& inflow_bc_descriptor =
+          dynamic_cast<const InflowBoundaryConditionDescriptor&>(*i_bc);
+
+        MeshFaceBoundary<Dimension> mesh_face_boundary =
+          getMeshFaceBoundary(mesh, inflow_bc_descriptor.boundaryDescriptor());
+
+        FaceValue<const Rd> xl = MeshDataManager::instance().getMeshData(*m_old_mesh).xl();
+        Array<const DataType> value_list =
+          InterpolateItemValue<DataType(Rd)>::template interpolate<ItemType::face>(inflow_bc_descriptor
+                                                                                     .functionSymbolId(),
+                                                                                   xl, mesh_face_boundary.faceList());
+
+        m_boundary_condition_list[i_quantity].emplace_back(
+          InflowValueBoundaryCondition<DataType>{mesh_face_boundary, value_list});
+        break;
+      }
+      default: {
+        std::ostringstream error_msg;
+        error_msg << "invalid boundary condition for advection: " << *i_bc;
+        throw NormalError(error_msg.str());
+      }
+      }
+    }
+  }
+
+  void
+  _storeArrayBCList(const size_t i_quantity,
+                    const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_list)
+  {
+    const Mesh<Connectivity<Dimension>>& mesh = *m_old_mesh;
+
+    for (auto& i_bc : bc_list) {
+      switch (i_bc->type()) {
+      case IBoundaryConditionDescriptor::Type::symmetry: {
+        m_boundary_condition_list[i_quantity].emplace_back(
+          SymmetryBoundaryCondition(getMeshFlatFaceBoundary(mesh, i_bc->boundaryDescriptor())));
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::outflow: {
+        m_boundary_condition_list[i_quantity].emplace_back(
+          OutflowBoundaryCondition(getMeshFaceBoundary(mesh, i_bc->boundaryDescriptor())));
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::inflow: {
+        const InflowBoundaryConditionDescriptor& inflow_bc_descriptor =
+          dynamic_cast<const InflowBoundaryConditionDescriptor&>(*i_bc);
+
+        MeshFaceBoundary<Dimension> mesh_face_boundary =
+          getMeshFaceBoundary(mesh, inflow_bc_descriptor.boundaryDescriptor());
+
+        FaceValue<const Rd> xl = MeshDataManager::instance().getMeshData(*m_old_mesh).xl();
+        Table<const double> array_list =
+          InterpolateItemArray<double(Rd)>::template interpolate<ItemType::face>({inflow_bc_descriptor
+                                                                                    .functionSymbolId()},
+                                                                                 xl, mesh_face_boundary.faceList());
+
+        m_boundary_condition_list[i_quantity].emplace_back(
+          InflowArrayBoundaryCondition{mesh_face_boundary, array_list});
+        break;
+      }
+      default: {
+        std::ostringstream error_msg;
+        error_msg << "invalid boundary condition for advection: " << *i_bc;
+        throw NormalError(error_msg.str());
+      }
+      }
+    }
+  }
+
   template <typename CellDataType>
-  void _remapOne(const CellValue<const double>& step_Vj, CellDataType& old_q);
+  void _remapOne(const CellValue<const double>& step_Vj,
+                 const std::vector<BoundaryConditionVariant>& q_bc_list,
+                 CellDataType& old_q);
 
   void _remapAllQuantities();
 
  public:
   std::vector<std::shared_ptr<const DiscreteFunctionVariant>>   //
-  remap(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& quantities);
+  remap(const std::vector<std::shared_ptr<const VariableBCDescriptor>>& quantity_list_with_bc);
 
   FluxingAdvectionSolver(const std::shared_ptr<const IMesh> i_old_mesh, const std::shared_ptr<const IMesh> i_new_mesh)
     : m_old_mesh{std::dynamic_pointer_cast<const MeshType>(i_old_mesh)},
@@ -102,16 +220,16 @@ FluxingAdvectionSolver<Dimension>::_computeDonorCells(FaceValue<const double> al
     FaceValue<CellId> donnor_cell{m_old_mesh->connectivity()};
     parallel_for(
       m_new_mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
-        const auto& face_to_cell = face_to_cell_matrix[face_id];
-        if (face_to_cell.size() == 1) {
-          donnor_cell[face_id] = face_to_cell[0];
+        const auto& face_to_cell    = face_to_cell_matrix[face_id];
+        const size_t i_face_in_cell = face_local_number_in_their_cells[face_id][0];
+        const CellId cell_id        = face_to_cell[0];
+        if (cell_face_is_reversed[cell_id][i_face_in_cell] xor (algebraic_fluxing_volumes[face_id] <= 0)) {
+          donnor_cell[face_id] = cell_id;
         } else {
-          const CellId cell_id        = face_to_cell[0];
-          const size_t i_face_in_cell = face_local_number_in_their_cells[face_id][0];
-          if (cell_face_is_reversed[cell_id][i_face_in_cell] xor (algebraic_fluxing_volumes[face_id] <= 0)) {
-            donnor_cell[face_id] = cell_id;
-          } else {
+          if (face_to_cell.size() == 2) {
             donnor_cell[face_id] = face_to_cell[1];
+          } else {
+            donnor_cell[face_id] = std::numeric_limits<CellId::base_type>::max();
           }
         }
       });
@@ -128,14 +246,20 @@ FluxingAdvectionSolver<1>::_computeDonorCells(FaceValue<const double> algebraic_
     const auto face_to_cell_matrix = m_new_mesh->connectivity().faceToCellMatrix();
     const auto cell_to_face_matrix = m_new_mesh->connectivity().cellToFaceMatrix();
 
+    const auto face_local_number_in_their_cells = m_new_mesh->connectivity().faceLocalNumbersInTheirCells();
+
     FaceValue<CellId> donnor_cell{m_old_mesh->connectivity()};
     parallel_for(
       m_new_mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
         const auto& face_to_cell = face_to_cell_matrix[face_id];
+        const CellId cell_id     = face_to_cell[0];
         if (face_to_cell.size() == 1) {
-          donnor_cell[face_id] = face_to_cell[0];
+          if ((algebraic_fluxing_volumes[face_id] <= 0) xor (cell_to_face_matrix[cell_id][1] == face_id)) {
+            donnor_cell[face_id] = std::numeric_limits<CellId::base_type>::max();
+          } else {
+            donnor_cell[face_id] = cell_id;
+          }
         } else {
-          const CellId cell_id = face_to_cell[0];
           if ((algebraic_fluxing_volumes[face_id] <= 0) xor (cell_to_face_matrix[cell_id][0] == face_id)) {
             donnor_cell[face_id] = cell_id;
           } else {
@@ -299,8 +423,9 @@ FluxingAdvectionSolver<Dimension>::_computeCycleNumber(FaceValue<double> fluxing
   CellValue<size_t> ratio(m_old_mesh->connectivity());
 
   parallel_for(
-    m_old_mesh->numberOfCells(),
-    PUGS_LAMBDA(CellId cell_id) { ratio[cell_id] = std::ceil(total_negative_flux[cell_id] / Vj[cell_id]); });
+    m_old_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+      ratio[cell_id] = (1. / 0.6) * std::ceil(total_negative_flux[cell_id] / Vj[cell_id]);
+    });
   synchronize(ratio);
 
   size_t number_of_cycles = max(ratio);
@@ -329,8 +454,12 @@ FluxingAdvectionSolver<Dimension>::_computeGeometricalData()
 template <size_t Dimension>
 template <typename CellDataType>
 void
-FluxingAdvectionSolver<Dimension>::_remapOne(const CellValue<const double>& step_Vj, CellDataType& old_q)
+FluxingAdvectionSolver<Dimension>::_remapOne(const CellValue<const double>& step_Vj,
+                                             const std::vector<BoundaryConditionVariant>& q_bc_list,
+                                             CellDataType& old_q)
 {
+  using DataType = std::decay_t<typename CellDataType::data_type>;
+
   static_assert(is_item_value_v<CellDataType> or is_item_array_v<CellDataType>, "invalid data type");
 
   const auto cell_to_face_matrix = m_new_mesh->connectivity().cellToFaceMatrix();
@@ -353,6 +482,7 @@ FluxingAdvectionSolver<Dimension>::_remapOne(const CellValue<const double>& step
       });
   }
 
+  // First we deal with inner faces
   parallel_for(
     m_new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
       const auto& cell_to_face = cell_to_face_matrix[cell_id];
@@ -384,6 +514,119 @@ FluxingAdvectionSolver<Dimension>::_remapOne(const CellValue<const double>& step
       }
     });
 
+  // Now we deal with boundary faces
+  for (const auto& bc_descriptor : q_bc_list) {
+    std::visit(
+      [&](auto&& bc) {
+        using TypeOfBC = std::decay_t<decltype(bc)>;
+
+        if constexpr (std::is_same_v<TypeOfBC, SymmetryBoundaryCondition>) {
+          auto face_list = bc.faceList();
+          for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+            const FaceId face_id        = face_list[i_face];
+            const double fluxing_volume = m_cycle_fluxing_volume[face_id];
+            const CellId face_cell_id   = face_to_cell_matrix[face_id][0];
+            if (fluxing_volume > 1E-12 * step_Vj[face_cell_id]) {
+              std::ostringstream error_msg;
+              error_msg << "invalid symmetry for face " << face_id
+                        << " (number=" << m_old_mesh->connectivity().faceNumber()[face_id] << ")\n"
+                        << "face has non zero fluxing volume " << fluxing_volume
+                        << " (cell volume = " << step_Vj[face_cell_id] << ")\n";
+              throw NormalError(error_msg.str());
+            }
+          }
+        } else if constexpr (std::is_same_v<TypeOfBC, OutflowBoundaryCondition>) {
+          auto face_list = bc.faceList();
+          for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+            const FaceId face_id        = face_list[i_face];
+            const CellId cell_id        = m_donnor_cell[face_id];
+            const double fluxing_volume = m_cycle_fluxing_volume[face_id];
+            if (cell_id != std::numeric_limits<CellId::base_type>::max()) {
+              if constexpr (is_item_array_v<CellDataType>) {
+                for (size_t i = 0; i < new_q[cell_id].size(); ++i) {
+                  new_q[cell_id][i] -= fluxing_volume * old_q[cell_id][i];
+                }
+              } else {
+                new_q[cell_id] -= fluxing_volume * old_q[cell_id];
+              }
+            } else {
+              const CellId face_cell_id = face_to_cell_matrix[face_id][0];
+              if (fluxing_volume > 1E-12 * step_Vj[face_cell_id]) {
+                std::ostringstream error_msg;
+                error_msg << "invalid outflow for face " << face_id
+                          << " (number=" << m_old_mesh->connectivity().faceNumber()[face_id] << ")"
+                          << "face has inflow fluxing volume " << fluxing_volume
+                          << " (cell volume = " << step_Vj[face_cell_id] << ")\n";
+                throw NormalError(error_msg.str());
+              }
+            }
+          }
+        } else if constexpr (std::is_same_v<TypeOfBC, InflowValueBoundaryCondition<DataType>>) {
+          if constexpr (is_item_value_v<CellDataType>) {
+            auto face_list  = bc.faceList();
+            auto value_list = bc.valueList();
+            for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+              const FaceId face_id        = face_list[i_face];
+              const CellId cell_id        = m_donnor_cell[face_id];
+              const double fluxing_volume = m_cycle_fluxing_volume[face_id];
+              if (cell_id == std::numeric_limits<CellId::base_type>::max()) {
+                Assert(face_to_cell_matrix[face_id].size() == 1, "boundary face must be connected to a single cell");
+                const CellId face_cell_id = face_to_cell_matrix[face_id][0];
+                new_q[face_cell_id] += fluxing_volume * value_list[i_face];
+              } else {
+                const CellId face_cell_id = face_to_cell_matrix[face_id][0];
+                if (fluxing_volume > 1E-12 * step_Vj[face_cell_id]) {
+                  std::ostringstream error_msg;
+                  error_msg << "invalid inflow for face " << face_id
+                            << " (number=" << m_old_mesh->connectivity().faceNumber()[face_id] << ")"
+                            << "face has outflow fluxing volume " << fluxing_volume
+                            << " (cell volume = " << step_Vj[face_cell_id] << ")\n";
+                  throw NormalError(error_msg.str());
+                }
+              }
+            }
+          } else {
+            throw UnexpectedError("invalid boundary condition for fluxing advection");
+          }
+        } else if constexpr (std::is_same_v<TypeOfBC, InflowArrayBoundaryCondition>) {
+          if constexpr (is_item_array_v<CellDataType>) {
+            auto face_list  = bc.faceList();
+            auto array_list = bc.arrayList();
+            for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+              const FaceId face_id        = face_list[i_face];
+              const CellId cell_id        = m_donnor_cell[face_id];
+              const double fluxing_volume = m_cycle_fluxing_volume[face_id];
+              if (cell_id == std::numeric_limits<CellId::base_type>::max()) {
+                Assert(face_to_cell_matrix[face_id].size() == 1, "boundary face must be connected to a single cell");
+                const CellId face_cell_id = face_to_cell_matrix[face_id][0];
+                auto new_array            = new_q[face_cell_id];
+                const auto& bc_array      = array_list[i_face];
+
+                for (size_t i = 0; i < new_array.size(); ++i) {
+                  new_array[i] += fluxing_volume * bc_array[i];
+                }
+              } else {
+                const CellId face_cell_id = face_to_cell_matrix[face_id][0];
+                if (fluxing_volume > 1E-12 * step_Vj[face_cell_id]) {
+                  std::ostringstream error_msg;
+                  error_msg << "invalid inflow for face " << face_id
+                            << " (number=" << m_old_mesh->connectivity().faceNumber()[face_id] << ")"
+                            << "face has outflow fluxing volume " << fluxing_volume
+                            << " (cell volume = " << step_Vj[face_cell_id] << ")\n";
+                  throw NormalError(error_msg.str());
+                }
+              }
+            }
+          } else {
+            throw UnexpectedError("invalid boundary condition for fluxing advection");
+          }
+        } else {
+          throw UnexpectedError("invalid boundary condition for fluxing advection");
+        }
+      },
+      bc_descriptor);
+  }
+
   synchronize(new_q);
   old_q = new_q;
 }
@@ -399,10 +642,12 @@ FluxingAdvectionSolver<Dimension>::_remapAllQuantities()
 
   const CellValue<const double> old_Vj = old_mesh_data.Vj();
   const CellValue<double> step_Vj      = copy(old_Vj);
-
   for (size_t jstep = 0; jstep < m_number_of_cycles; ++jstep) {
-    for (auto& remapped_q : m_remapped_list) {
-      std::visit([&](auto&& old_q) { this->_remapOne(step_Vj, old_q); }, remapped_q);
+    for (size_t i = 0; i < m_remapped_list.size(); ++i) {
+      Assert(m_remapped_list.size() == m_boundary_condition_list.size());
+      auto& remapped_q = m_remapped_list[i];
+      auto& q_bc_list  = m_boundary_condition_list[i];
+      std::visit([&](auto&& old_q) { this->_remapOne(step_Vj, q_bc_list, old_q); }, remapped_q);
     }
 
     parallel_for(
@@ -452,26 +697,38 @@ FluxingAdvectionSolver<Dimension>::_remapAllQuantities()
 template <size_t Dimension>
 std::vector<std::shared_ptr<const DiscreteFunctionVariant>>
 FluxingAdvectionSolver<Dimension>::remap(
-  const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& quantity_list)
+  const std::vector<std::shared_ptr<const VariableBCDescriptor>>& quantity_list_with_bc)
 {
-  for (auto&& variable_v : quantity_list) {
+  m_boundary_condition_list.resize(quantity_list_with_bc.size());
+  for (size_t i = 0; i < quantity_list_with_bc.size(); ++i) {
+    const auto& quantity_v_with_bc = quantity_list_with_bc[i];
+    const auto& quantity_v         = quantity_v_with_bc->discreteFunctionVariant();
+    const auto& bc_list            = quantity_v_with_bc->bcDescriptorList();
     std::visit(
       [&](auto&& variable) {
         using DiscreteFunctionT = std::decay_t<decltype(variable)>;
         if constexpr (std::is_same_v<MeshType, typename DiscreteFunctionT::MeshType>) {
           this->_storeValues(variable);
+          if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
+            using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
+            this->_storeValueBCList<DataType>(i, bc_list);
+          } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
+            this->_storeArrayBCList(i, bc_list);
+          } else {
+            throw UnexpectedError("invalid discrete function type");
+          }
         } else {
           throw UnexpectedError("incompatible mesh types");
         }
       },
-      variable_v->discreteFunction());
+      quantity_v->discreteFunction());
   }
 
   this->_remapAllQuantities();
 
   std::vector<std::shared_ptr<const DiscreteFunctionVariant>> new_variables;
 
-  for (size_t i = 0; i < quantity_list.size(); ++i) {
+  for (size_t i = 0; i < quantity_list_with_bc.size(); ++i) {
     std::visit(
       [&](auto&& variable) {
         using DiscreteFunctionT = std::decay_t<decltype(variable)>;
@@ -491,16 +748,126 @@ FluxingAdvectionSolver<Dimension>::remap(
           throw UnexpectedError("incompatible mesh types");
         }
       },
-      quantity_list[i]->discreteFunction());
+      quantity_list_with_bc[i]->discreteFunctionVariant()->discreteFunction());
   }
 
   return new_variables;
 }
 
+template <size_t Dimension>
+template <typename DataType>
+class FluxingAdvectionSolver<Dimension>::InflowValueBoundaryCondition
+{
+ private:
+  const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
+
+  Array<const DataType> m_value_list;
+
+ public:
+  const Array<const FaceId>&
+  faceList() const
+  {
+    return m_mesh_face_boundary.faceList();
+  }
+
+  const Array<const DataType>&
+  valueList() const
+  {
+    return m_value_list;
+  }
+
+  InflowValueBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary, Array<const DataType> value_list)
+    : m_mesh_face_boundary(mesh_face_boundary), m_value_list(value_list)
+  {
+    ;
+  }
+
+  ~InflowValueBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class FluxingAdvectionSolver<Dimension>::InflowArrayBoundaryCondition
+{
+ private:
+  const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
+
+  Table<const double> m_array_list;
+
+ public:
+  const Array<const FaceId>&
+  faceList() const
+  {
+    return m_mesh_face_boundary.faceList();
+  }
+
+  const Table<const double>&
+  arrayList() const
+  {
+    return m_array_list;
+  }
+
+  InflowArrayBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary, Table<const double> array_list)
+    : m_mesh_face_boundary(mesh_face_boundary), m_array_list(array_list)
+  {
+    ;
+  }
+
+  ~InflowArrayBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class FluxingAdvectionSolver<Dimension>::OutflowBoundaryCondition
+{
+ private:
+  const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
+
+ public:
+  const Array<const FaceId>&
+  faceList() const
+  {
+    return m_mesh_face_boundary.faceList();
+  }
+
+  OutflowBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary)
+    : m_mesh_face_boundary(mesh_face_boundary)
+  {
+    ;
+  }
+
+  ~OutflowBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class FluxingAdvectionSolver<Dimension>::SymmetryBoundaryCondition
+{
+ private:
+  const MeshFlatFaceBoundary<Dimension> m_mesh_flat_face_boundary;
+
+ public:
+  const Array<const FaceId>&
+  faceList() const
+  {
+    return m_mesh_flat_face_boundary.faceList();
+  }
+
+  SymmetryBoundaryCondition(const MeshFlatFaceBoundary<Dimension>& mesh_flat_face_boundary)
+    : m_mesh_flat_face_boundary(mesh_flat_face_boundary)
+  {
+    ;
+  }
+
+  ~SymmetryBoundaryCondition() = default;
+};
+
 std::vector<std::shared_ptr<const DiscreteFunctionVariant>>
 advectByFluxing(const std::shared_ptr<const IMesh> i_new_mesh,
-                const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& remapped_variables)
+                const std::vector<std::shared_ptr<const VariableBCDescriptor>>& remapped_variables_with_bc)
 {
+  std::vector<std::shared_ptr<const DiscreteFunctionVariant>> remapped_variables;
+  for (auto i_remapped_var_with_bc : remapped_variables_with_bc) {
+    remapped_variables.push_back(i_remapped_var_with_bc->discreteFunctionVariant());
+  }
+
   if (not hasSameMesh(remapped_variables)) {
     throw NormalError("remapped quantities are not defined on the same mesh");
   }
@@ -509,13 +876,13 @@ advectByFluxing(const std::shared_ptr<const IMesh> i_new_mesh,
 
   switch (i_old_mesh->dimension()) {
   case 1: {
-    return FluxingAdvectionSolver<1>{i_old_mesh, i_new_mesh}.remap(remapped_variables);
+    return FluxingAdvectionSolver<1>{i_old_mesh, i_new_mesh}.remap(remapped_variables_with_bc);
   }
   case 2: {
-    return FluxingAdvectionSolver<2>{i_old_mesh, i_new_mesh}.remap(remapped_variables);
+    return FluxingAdvectionSolver<2>{i_old_mesh, i_new_mesh}.remap(remapped_variables_with_bc);
   }
   case 3: {
-    return FluxingAdvectionSolver<3>{i_old_mesh, i_new_mesh}.remap(remapped_variables);
+    return FluxingAdvectionSolver<3>{i_old_mesh, i_new_mesh}.remap(remapped_variables_with_bc);
   }
   default: {
     throw UnexpectedError("Invalid mesh dimension");
diff --git a/src/scheme/FluxingAdvectionSolver.hpp b/src/scheme/FluxingAdvectionSolver.hpp
index 05dce50551c6caf20eb4efb25ee8939125c567d0..2cbb653cbeed76000d8eb1e3ba01ed61c35bc8a3 100644
--- a/src/scheme/FluxingAdvectionSolver.hpp
+++ b/src/scheme/FluxingAdvectionSolver.hpp
@@ -3,6 +3,7 @@
 
 #include <language/utils/FunctionSymbolId.hpp>
 #include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/VariableBCDescriptor.hpp>
 
 #include <vector>
 
@@ -10,4 +11,8 @@ std::vector<std::shared_ptr<const DiscreteFunctionVariant>> advectByFluxing(
   const std::shared_ptr<const IMesh> new_mesh,
   const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& remapped_variables);
 
+std::vector<std::shared_ptr<const DiscreteFunctionVariant>> advectByFluxing(
+  const std::shared_ptr<const IMesh> new_mesh,
+  const std::vector<std::shared_ptr<const VariableBCDescriptor>>& remapped_variables_with_bc);
+
 #endif   // FLUXING_ADVECION_SOLVER_HPP
diff --git a/src/scheme/IBoundaryConditionDescriptor.hpp b/src/scheme/IBoundaryConditionDescriptor.hpp
index f99eff9355f6fdf6fbd65ad1898f3cba655af284..1350aa4dba4f7f6739b7de6689c10554ceec6971 100644
--- a/src/scheme/IBoundaryConditionDescriptor.hpp
+++ b/src/scheme/IBoundaryConditionDescriptor.hpp
@@ -16,7 +16,9 @@ class IBoundaryConditionDescriptor
     fourier,
     fixed,
     free,
+    inflow,
     neumann,
+    outflow,
     symmetry
   };
 
diff --git a/src/scheme/InflowBoundaryConditionDescriptor.hpp b/src/scheme/InflowBoundaryConditionDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e648c1593cb4316f0cc77a4f18946b0457d65b5e
--- /dev/null
+++ b/src/scheme/InflowBoundaryConditionDescriptor.hpp
@@ -0,0 +1,55 @@
+#ifndef INFLOW_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#define INFLOW_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IBoundaryDescriptor.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+
+#include <memory>
+
+class InflowBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+{
+ private:
+  std::ostream&
+  _write(std::ostream& os) const final
+  {
+    os << "inflow(" << ',' << *m_boundary_descriptor << ")";
+    return os;
+  }
+
+  std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
+  const FunctionSymbolId m_function_symbol_id;
+
+ public:
+  FunctionSymbolId
+  functionSymbolId() const
+  {
+    return m_function_symbol_id;
+  }
+
+  const IBoundaryDescriptor&
+  boundaryDescriptor() const final
+  {
+    return *m_boundary_descriptor;
+  }
+
+  Type
+  type() const final
+  {
+    return Type::inflow;
+  }
+
+  InflowBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
+                                    const FunctionSymbolId& function_symbol_id)
+    : m_boundary_descriptor(boundary_descriptor), m_function_symbol_id{function_symbol_id}
+  {
+    ;
+  }
+
+  InflowBoundaryConditionDescriptor(const InflowBoundaryConditionDescriptor&) = delete;
+  InflowBoundaryConditionDescriptor(InflowBoundaryConditionDescriptor&&)      = delete;
+
+  ~InflowBoundaryConditionDescriptor() = default;
+};
+
+#endif   // INFLOW_BOUNDARY_CONDITION_DESCRIPTOR_HPP
diff --git a/src/scheme/OutflowBoundaryConditionDescriptor.hpp b/src/scheme/OutflowBoundaryConditionDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..83ef7a775a81269a5d3930ed2e07d5634d99959a
--- /dev/null
+++ b/src/scheme/OutflowBoundaryConditionDescriptor.hpp
@@ -0,0 +1,46 @@
+#ifndef OUTFLOW_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#define OUTFLOW_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+
+#include <mesh/IBoundaryDescriptor.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+
+#include <memory>
+
+class OutflowBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+{
+ private:
+  std::ostream&
+  _write(std::ostream& os) const final
+  {
+    os << "outflow(" << *m_boundary_descriptor << ")";
+    return os;
+  }
+
+  std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
+
+ public:
+  const IBoundaryDescriptor&
+  boundaryDescriptor() const final
+  {
+    return *m_boundary_descriptor;
+  }
+
+  Type
+  type() const final
+  {
+    return Type::outflow;
+  }
+
+  OutflowBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor)
+    : m_boundary_descriptor(boundary_descriptor)
+  {
+    ;
+  }
+
+  OutflowBoundaryConditionDescriptor(const OutflowBoundaryConditionDescriptor&) = delete;
+  OutflowBoundaryConditionDescriptor(OutflowBoundaryConditionDescriptor&&)      = delete;
+
+  ~OutflowBoundaryConditionDescriptor() = default;
+};
+
+#endif   // OUTFLOW_BOUNDARY_CONDITION_DESCRIPTOR_HPP
diff --git a/src/scheme/VariableBCDescriptor.hpp b/src/scheme/VariableBCDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..18b35925433bc98148a87485e97e8edddfa20fbc
--- /dev/null
+++ b/src/scheme/VariableBCDescriptor.hpp
@@ -0,0 +1,41 @@
+#ifndef VARIABLE_BC_DESCRIPTOR_HPP
+#define VARIABLE_BC_DESCRIPTOR_HPP
+
+class DiscreteFunctionVariant;
+class IBoundaryConditionDescriptor;
+
+#include <memory>
+#include <vector>
+
+class VariableBCDescriptor
+{
+ private:
+  std::shared_ptr<const DiscreteFunctionVariant> m_discrete_function_variant;
+  std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>> m_bc_descriptor_list;
+
+ public:
+  const std::shared_ptr<const DiscreteFunctionVariant>&
+  discreteFunctionVariant() const
+  {
+    return m_discrete_function_variant;
+  }
+
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+  bcDescriptorList() const
+  {
+    return m_bc_descriptor_list;
+  }
+
+  VariableBCDescriptor(const std::shared_ptr<const DiscreteFunctionVariant>& discrete_function_variant,
+                       const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+    : m_discrete_function_variant{discrete_function_variant}, m_bc_descriptor_list{bc_descriptor_list}
+  {}
+
+  VariableBCDescriptor(const VariableBCDescriptor&) = default;
+  VariableBCDescriptor(VariableBCDescriptor&&)      = default;
+
+  VariableBCDescriptor()  = default;
+  ~VariableBCDescriptor() = default;
+};
+
+#endif   // VARIABLE_BC_DESCRIPTOR_HPP
diff --git a/tests/test_BuiltinFunctionEmbedder.cpp b/tests/test_BuiltinFunctionEmbedder.cpp
index 2e56fe228bc4b5d0ae863ae177b3633efbfa6e90..31cde714627687bb0aee8f62ff122cf2e5ae0a1b 100644
--- a/tests/test_BuiltinFunctionEmbedder.cpp
+++ b/tests/test_BuiltinFunctionEmbedder.cpp
@@ -310,6 +310,64 @@ TEST_CASE("BuiltinFunctionEmbedder", "[language]")
     REQUIRE(*data_type.contentTypeList()[2] == ast_node_data_type_from<std::shared_ptr<const double>>);
   }
 
+  SECTION("(R^2) -> (R) BuiltinFunctionEmbedder")
+  {
+    std::function c = [](const std::vector<TinyVector<2>>& x) -> std::vector<double> {
+      std::vector<double> sum(x.size());
+      for (size_t i = 0; i < x.size(); ++i) {
+        sum[i] = x[i][0] + x[i][1];
+      }
+      return sum;
+    };
+
+    std::unique_ptr<IBuiltinFunctionEmbedder> i_embedded_c =
+      std::make_unique<BuiltinFunctionEmbedder<std::vector<double>(const std::vector<TinyVector<2>>&)>>(c);
+
+    REQUIRE(i_embedded_c->numberOfParameters() == 1);
+    REQUIRE(i_embedded_c->getParameterDataTypes().size() == 1);
+    REQUIRE(i_embedded_c->getParameterDataTypes()[0] == ASTNodeDataType::tuple_t);
+    REQUIRE(i_embedded_c->getReturnDataType() == ASTNodeDataType::tuple_t);
+
+    std::vector<TinyVector<2>> x = {TinyVector<2>{1, 2}, TinyVector<2>{3, 4}, TinyVector<2>{-2, 4}};
+    std::vector<double> result   = c(x);
+
+    const std::vector value = std::get<std::vector<double>>(i_embedded_c->apply(std::vector<DataVariant>{x}));
+
+    for (size_t i = 0; i < result.size(); ++i) {
+      REQUIRE(value[i] == result[i]);
+    }
+  }
+
+  SECTION("std::vector<EmbeddedData>(EmbeddedData, EmbeddedData) BuiltinFunctionEmbedder")
+  {
+    std::function sum = [&](const uint64_t& i, const uint64_t& j) -> std::vector<std::shared_ptr<const uint64_t>> {
+      std::vector<std::shared_ptr<const uint64_t>> x = {std::make_shared<const uint64_t>(i),
+                                                        std::make_shared<const uint64_t>(j)};
+      return x;
+    };
+
+    std::unique_ptr<IBuiltinFunctionEmbedder> i_embedded_c = std::make_unique<
+      BuiltinFunctionEmbedder<std::vector<std::shared_ptr<const uint64_t>>(const uint64_t&, const uint64_t&)>>(sum);
+
+    REQUIRE(i_embedded_c->numberOfParameters() == 2);
+    REQUIRE(i_embedded_c->getParameterDataTypes().size() == 2);
+    REQUIRE(i_embedded_c->getParameterDataTypes()[0] == ASTNodeDataType::unsigned_int_t);
+    REQUIRE(i_embedded_c->getParameterDataTypes()[1] == ASTNodeDataType::unsigned_int_t);
+    REQUIRE(i_embedded_c->getReturnDataType() == ASTNodeDataType::tuple_t);
+
+    std::vector<uint64_t> values{3, 4};
+
+    std::vector<EmbeddedData> embedded_data_list =
+      std::get<std::vector<EmbeddedData>>(i_embedded_c->apply(std::vector<DataVariant>{values[0], values[1]}));
+
+    for (size_t i = 0; i < embedded_data_list.size(); ++i) {
+      const EmbeddedData& embedded_data       = embedded_data_list[i];
+      const IDataHandler& handled_data        = embedded_data.get();
+      const DataHandler<const uint64_t>& data = dynamic_cast<const DataHandler<const uint64_t>&>(handled_data);
+      REQUIRE(*data.data_ptr() == values[i]);
+    }
+  }
+
   SECTION("void -> N*R*shared_const_double BuiltinFunctionEmbedder")
   {
     std::function c = [](void) -> std::tuple<uint64_t, double, std::shared_ptr<const double>> {
diff --git a/tests/test_InterpolateItemArray.cpp b/tests/test_InterpolateItemArray.cpp
index ccee6e77877eadc359cc2c11313a7114081cc8c0..6d020282ad0d0a953ca7e815357d5461f5dd3a76 100644
--- a/tests/test_InterpolateItemArray.cpp
+++ b/tests/test_InterpolateItemArray.cpp
@@ -48,69 +48,192 @@ 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));
+            }
+
+            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;
+              });
 
-            function_symbol_id_list.push_back(
-              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            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();
+
+            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"};
 
-            function_symbol_id_list.push_back(
-              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            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));
+            }
+
+            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));
           }
+        }
+      }
 
-          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;
-            });
+      SECTION("from -> (R^1)")
+      {
+        for (const auto& named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_1d = named_mesh.mesh();
 
-          CellArray<const double> interpolate_array =
-            InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj);
+            auto xj = MeshDataManager::instance().getMeshData(*mesh_1d).xj();
 
-          REQUIRE(same_cell_array(cell_array, interpolate_array));
+            std::string_view data = R"(
+import math;
+let f_1d: R^1 -> (R^1), 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));
+            }
+
+            using R1 = TinyVector<1>;
+            CellArray<R1> 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] = R1{2 * x[0] + 2};
+                cell_array[cell_id][1] = R1{2 * exp(x[0]) + 3};
+              });
+
+            CellArray<const R1> interpolate_array =
+              InterpolateItemArray<R1(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj);
+
+            REQUIRE(same_cell_array(cell_array, interpolate_array));
+          }
         }
       }
     }
@@ -121,69 +244,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));
+            }
+
+            {
+              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));
+            }
+
+            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();
 
-          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;
-            });
+            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);
+
+            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<const double> interpolate_array =
-            InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj);
+            std::vector<FunctionSymbolId> function_symbol_id_list;
 
-          REQUIRE(same_cell_array(cell_array, interpolate_array));
+            {
+              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 +379,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);
 
-          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;
-            });
+            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_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);
+
+            REQUIRE(same_cell_array(cell_array, interpolate_array));
+          }
         }
       }
     }
@@ -264,7 +511,7 @@ let scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
 
   SECTION("interpolate on items list")
   {
-    auto same_cell_value = [](auto interpolated, auto reference) -> bool {
+    auto same_cell_array = [](auto interpolated, auto reference) -> bool {
       for (size_t i = 0; i < interpolated.numberOfRows(); ++i) {
         for (size_t j = 0; j < interpolated.numberOfColumns(); ++j) {
           if (interpolated[i][j] != reference[i][j]) {
@@ -281,77 +528,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_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();
+
+            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};
 
-          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;
-            });
+            ASTSymbolTableBuilder{*ast};
+            ASTNodeDataTypeBuilder{*ast};
 
-          Table<const double> interpolate_array =
-            InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj, cell_id_list);
+            ASTNodeTypeCleaner<language::var_declaration>{*ast};
+            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+            ASTNodeExpressionBuilder{*ast};
 
-          REQUIRE(same_cell_value(cell_array, interpolate_array));
+            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<const double> interpolate_array =
+              InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj,
+                                                                               cell_id_list);
+
+            REQUIRE(same_cell_array(cell_array, interpolate_array));
+          }
         }
       }
     }
@@ -362,74 +681,213 @@ 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));
+            }
 
-            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));
+            }
+
+            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_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();
+
+            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"};
+
+            auto ast = ASTBuilder::build(input);
+
+            ASTModulesImporter{*ast};
+            ASTNodeTypeCleaner<language::import_instruction>{*ast};
 
-            function_symbol_id_list.push_back(
-              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+            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_array(cell_array, interpolate_array));
           }
+        }
+      }
+
+      SECTION("from -> (R^2x2)")
+      {
+        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();
 
-          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;
-            });
+            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^2x2), x -> ([[x[0],0],[2-x[1], x[0]*x[1]]], [[2*x[0], x[1]],[2, -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
+
+            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<const double> interpolate_array =
-            InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj, cell_id_list);
+            using R2x2 = TinyMatrix<2>;
 
-          REQUIRE(same_cell_value(cell_array, interpolate_array));
+            Table<R2x2> 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] = R2x2{x[0], 0,   //
+                                        2 - x[1], x[0] * x[1]};
+
+                cell_array[i][1] = R2x2{2 * x[0], x[1],   //
+                                        2, -x[1]};
+              });
+
+            Table<const R2x2> interpolate_array =
+              InterpolateItemArray<R2x2(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj, cell_id_list);
+
+            REQUIRE(same_cell_array(cell_array, interpolate_array));
+          }
         }
       }
     }
@@ -440,74 +898,210 @@ 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_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();
+
+            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);
+
+            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;
 
-            function_symbol_id_list.push_back(
-              FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), 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_array(cell_array, interpolate_array));
           }
+        }
+      }
+
+      SECTION("from -> (R^3)")
+      {
+        for (const auto& named_mesh : mesh_list) {
+          SECTION(named_mesh.name())
+          {
+            auto mesh_3d = named_mesh.mesh();
 
-          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;
-            });
+            auto xj = MeshDataManager::instance().getMeshData(*mesh_3d).xj();
 
-          Table<const double> interpolate_array =
-            InterpolateItemArray<double(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj, cell_id_list);
+            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^3), x -> (2*x, [2*x[0]-x[1], 3*x[2]-x[0], x[1]+x[2]], 0);
+)";
+            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
 
-          REQUIRE(same_cell_value(cell_array, interpolate_array));
+            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_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));
+            }
+
+            using R3 = TinyVector<3>;
+            Table<R3> cell_array{cell_id_list.size(), 3};
+            parallel_for(
+              cell_id_list.size(), PUGS_LAMBDA(const size_t i) {
+                const R3& x = xj[cell_id_list[i]];
+
+                cell_array[i][0] = 2 * x;
+                cell_array[i][1] = R3{2 * x[0] - x[1], 3 * x[2] - x[0], x[1] + x[2]};
+                cell_array[i][2] = zero;
+              });
+
+            Table<const R3> interpolate_array =
+              InterpolateItemArray<R3(TinyVector<Dimension>)>::interpolate(function_symbol_id_list, xj, cell_id_list);
+
+            REQUIRE(same_cell_array(cell_array, interpolate_array));
+          }
         }
       }
     }
diff --git a/tests/test_ItemArrayVariantFunctionInterpoler.cpp b/tests/test_ItemArrayVariantFunctionInterpoler.cpp
index 640e3e77aa5e6ede05fe0f217d7587f7839e4d56..7ceaf2207ef53526f19511aa31415ce8efcd577a 100644
--- a/tests/test_ItemArrayVariantFunctionInterpoler.cpp
+++ b/tests/test_ItemArrayVariantFunctionInterpoler.cpp
@@ -386,6 +386,33 @@ let R3x3_non_linear_1d: R^1 -> R^3x3, x -> [[2 * exp(x[0]) * sin(x[0]) + 3, sin(
 
           REQUIRE(same_item_array(cell_array, item_array_variant->get<CellArray<DataType>>()));
         }
+
+        SECTION("different types")
+        {
+          auto [i_symbol_f1, found_f1] = symbol_table->find("R_scalar_non_linear1_1d", position);
+          REQUIRE(found_f1);
+          REQUIRE(i_symbol_f1->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function1_symbol_id(std::get<uint64_t>(i_symbol_f1->attributes().value()), symbol_table);
+
+          auto [i_symbol_f2, found_f2] = symbol_table->find("Z_scalar_non_linear2_1d", position);
+          REQUIRE(found_f2);
+          REQUIRE(i_symbol_f2->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function2_symbol_id(std::get<uint64_t>(i_symbol_f2->attributes().value()), symbol_table);
+
+          auto [i_symbol_f3, found_f3] = symbol_table->find("R_scalar_non_linear3_1d", position);
+          REQUIRE(found_f3);
+          REQUIRE(i_symbol_f3->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function3_symbol_id(std::get<uint64_t>(i_symbol_f3->attributes().value()), symbol_table);
+
+          ItemArrayVariantFunctionInterpoler interpoler(mesh_1d, ItemType::cell,
+                                                        {function1_symbol_id, function2_symbol_id,
+                                                         function3_symbol_id});
+
+          REQUIRE_THROWS_WITH(interpoler.interpolate(), "error: functions must have the same type");
+        }
       }
     }
   }