diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index 1717b48c78deaedc7cabd9f60feb215e29a80ebe..8e0bbd4db39f20977afa3725177f728d7094710d 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -43,6 +43,18 @@ class [[nodiscard]] TinyMatrix
     }
   }
 
+  template <typename... Args>
+  PUGS_FORCEINLINE constexpr void
+  _unpackVariadicInput(const TinyVector<M, T>& t, Args&&... args) noexcept
+  {
+    for (size_t j = 0; j < M; ++j) {
+      m_values[_index(j, N - 1 - sizeof...(args))] = t[j];
+    }
+    if constexpr (sizeof...(args) > 0) {
+      this->_unpackVariadicInput(std::forward<Args>(args)...);
+    }
+  }
+
  public:
   PUGS_INLINE
   constexpr bool
@@ -322,6 +334,13 @@ class [[nodiscard]] TinyMatrix
     this->_unpackVariadicInput(t, std::forward<Args>(args)...);
   }
 
+  template <typename... Args>
+  PUGS_INLINE explicit constexpr TinyMatrix(const TinyVector<M, T>& t, Args&&... args) noexcept
+  {
+    static_assert(sizeof...(args) == N - 1, "wrong number of parameters");
+    this->_unpackVariadicInput(t, std::forward<Args>(args)...);
+  }
+
   // One does not use the '=default' constructor to avoid (unexpected)
   // performances issues
   PUGS_INLINE
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 1c888821f8d6c9c62e8ddc58ce23cd97c4f05b32..813bde93a055a514e35a1e949ca6ea2241609ad2 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -16,6 +16,7 @@
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/MeshRandomizer.hpp>
+#include <mesh/MeshSmoother.hpp>
 #include <scheme/AcousticSolver.hpp>
 #include <scheme/AxisBoundaryConditionDescriptor.hpp>
 #include <scheme/DirichletBoundaryConditionDescriptor.hpp>
@@ -30,13 +31,17 @@
 #include <scheme/DiscreteFunctionVectorInterpoler.hpp>
 #include <scheme/ExternalBoundaryConditionDescriptor.hpp>
 #include <scheme/FixedBoundaryConditionDescriptor.hpp>
+#include <scheme/FluxingAdvectionSolver.hpp>
 #include <scheme/FourierBoundaryConditionDescriptor.hpp>
 #include <scheme/FreeBoundaryConditionDescriptor.hpp>
 #include <scheme/HyperelasticSolver.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>
@@ -48,6 +53,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(
 
@@ -73,6 +79,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> {
@@ -238,6 +254,58 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this->_addBuiltinFunction("smoothMesh", std::function(
+
+                                            [](std::shared_ptr<const IMesh> p_mesh,
+                                               const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                                 bc_descriptor_list) -> std::shared_ptr<const IMesh> {
+                                              MeshSmootherHandler handler;
+                                              return handler.getSmoothedMesh(p_mesh, bc_descriptor_list);
+                                            }
+
+                                            ));
+
+  this->_addBuiltinFunction("smoothMesh",
+                            std::function(
+
+                              [](std::shared_ptr<const IMesh> p_mesh,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const FunctionSymbolId& function_symbol_id) -> std::shared_ptr<const IMesh> {
+                                MeshSmootherHandler handler;
+                                return handler.getSmoothedMesh(p_mesh, bc_descriptor_list, function_symbol_id);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("smoothMesh",
+                            std::function(
+
+                              [](std::shared_ptr<const IMesh> p_mesh,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const std::vector<std::shared_ptr<const IZoneDescriptor>>& smoothing_zone_list)
+                                -> std::shared_ptr<const IMesh> {
+                                MeshSmootherHandler handler;
+                                return handler.getSmoothedMesh(p_mesh, bc_descriptor_list, smoothing_zone_list);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("smoothMesh", std::function(
+
+                                            [](std::shared_ptr<const IMesh> p_mesh,
+                                               const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                                 bc_descriptor_list,
+                                               const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>&
+                                                 discrete_function_variant_list) -> std::shared_ptr<const IMesh> {
+                                              MeshSmootherHandler handler;
+                                              return handler.getSmoothedMesh(p_mesh, bc_descriptor_list,
+                                                                             discrete_function_variant_list);
+                                            }
+
+                                            ));
+
   this->_addBuiltinFunction("fixed", std::function(
 
                                        [](std::shared_ptr<const IBoundaryDescriptor> boundary)
@@ -313,6 +381,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,
@@ -590,6 +678,17 @@ SchemeModule::SchemeModule()
 
                                                  ));
 
+  this->_addBuiltinFunction("fluxing_advection", std::function(
+
+                                                   [](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/BuiltinFunctionEmbedder.hpp b/src/language/utils/BuiltinFunctionEmbedder.hpp
index 460ec3e866be479633e5bf9a9a416c1f84fc09e6..251c39e6d5eb5dc46423e75dedea28947d1a29be 100644
--- a/src/language/utils/BuiltinFunctionEmbedder.hpp
+++ b/src/language/utils/BuiltinFunctionEmbedder.hpp
@@ -114,10 +114,21 @@ class BuiltinFunctionEmbedderBase<FX(Args...)> : public IBuiltinFunctionEmbedder
   }
 
   template <typename T>
-  PUGS_INLINE std::shared_ptr<IDataHandler>
+  PUGS_INLINE EmbeddedData
   _createHandler(std::shared_ptr<T> data) const
   {
-    return std::make_shared<DataHandler<T>>(data);
+    return EmbeddedData{std::make_shared<DataHandler<T>>(data)};
+  }
+
+  template <typename T>
+  PUGS_INLINE std::vector<EmbeddedData>
+  _createHandler(std::vector<std::shared_ptr<T>> data) const
+  {
+    std::vector<EmbeddedData> embedded(data.size());
+    for (size_t i_data = 0; i_data < data.size(); ++i_data) {
+      embedded[i_data] = EmbeddedData(std::make_shared<DataHandler<T>>(data[i_data]));
+    }
+    return embedded;
   }
 
   template <typename ResultT>
@@ -127,7 +138,7 @@ class BuiltinFunctionEmbedderBase<FX(Args...)> : public IBuiltinFunctionEmbedder
     if constexpr (is_data_variant_v<std::decay_t<ResultT>>) {
       return std::move(result);
     } else {
-      return EmbeddedData(_createHandler(std::move(result)));
+      return _createHandler(std::move(result));
     }
   }
 
@@ -327,7 +338,7 @@ class BuiltinFunctionEmbedder<FX(Args...)> : public BuiltinFunctionEmbedderBase<
       std::apply(m_f, t);
       return {};
     } else {
-      return EmbeddedData(this->template _createHandler(std::apply(m_f, t)));
+      return this->template _createHandler(std::apply(m_f, t));
     }
   }
 
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/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 62dbd50536ad2e725e67a25d9cddf481a62b40c6..39f31282296315f1dea4647bcb9d02aa43169707 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -22,6 +22,7 @@ add_library(
   MedianDualMeshBuilder.cpp
   MeshBuilderBase.cpp
   MeshCellZone.cpp
+  MeshData.cpp
   MeshDataManager.cpp
   MeshEdgeBoundary.cpp
   MeshFaceBoundary.cpp
@@ -34,5 +35,6 @@ add_library(
   MeshLineNodeBoundary.cpp
   MeshNodeBoundary.cpp
   MeshRandomizer.cpp
+  MeshSmoother.cpp
   MeshTransformer.cpp
   SynchronizerManager.cpp)
diff --git a/src/mesh/ItemArrayUtils.hpp b/src/mesh/ItemArrayUtils.hpp
index 6afa39c167c9d2fa95ab5ef53e502df27359dce7..abba5cfe2631bf6c472f3a4c60a5c0eec8703acb 100644
--- a/src/mesh/ItemArrayUtils.hpp
+++ b/src/mesh/ItemArrayUtils.hpp
@@ -308,7 +308,7 @@ sum(const ItemArray<DataType, item_type, ConnectivityPtr>& item_value)
 
 template <typename DataType, ItemType item_type, typename ConnectivityPtr>
 void
-synchronize(ItemArray<DataType, item_type, ConnectivityPtr>& item_array)
+synchronize(ItemArray<DataType, item_type, ConnectivityPtr> item_array)
 {
   static_assert(not std::is_const_v<DataType>, "cannot synchronize ItemArray of const data");
   if (parallel::size() > 1) {
diff --git a/src/mesh/ItemValueUtils.hpp b/src/mesh/ItemValueUtils.hpp
index ef99f59e44f9f71a01267ca0c424220241573f60..3bff4868c30677ef8d11db3a94cb493c5bcb07a0 100644
--- a/src/mesh/ItemValueUtils.hpp
+++ b/src/mesh/ItemValueUtils.hpp
@@ -296,7 +296,7 @@ sum(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
 
 template <typename DataType, ItemType item_type, typename ConnectivityPtr>
 void
-synchronize(ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
+synchronize(ItemValue<DataType, item_type, ConnectivityPtr> item_value)
 {
   static_assert(not std::is_const_v<DataType>, "cannot synchronize ItemValue of const data");
   if (parallel::size() > 1) {
@@ -309,7 +309,7 @@ synchronize(ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
 
 template <typename DataType, ItemType item_type, typename ConnectivityPtr>
 bool
-isSynchronized(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
+isSynchronized(ItemValue<DataType, item_type, ConnectivityPtr> item_value)
 {
   bool is_synchronized = true;
 
diff --git a/src/mesh/MeshData.cpp b/src/mesh/MeshData.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2fdddb387554258120bb36eba31e57a940e2576c
--- /dev/null
+++ b/src/mesh/MeshData.cpp
@@ -0,0 +1,25 @@
+#include <mesh/MeshData.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemValueVariant.hpp>
+#include <mesh/Mesh.hpp>
+#include <output/NamedItemValueVariant.hpp>
+#include <output/VTKWriter.hpp>
+#include <utils/Exceptions.hpp>
+
+template <size_t Dimension>
+void
+MeshData<Dimension>::_storeBadMesh()
+{
+  VTKWriter writer("bad_mesh");
+  writer.writeOnMesh(std::make_shared<MeshType>(m_mesh.shared_connectivity(), m_mesh.xr()),
+                     {std::make_shared<NamedItemValueVariant>(std::make_shared<ItemValueVariant>(m_Vj), "volume")});
+  std::ostringstream error_msg;
+  error_msg << "mesh contains cells of non-positive volume (see " << rang::fgB::yellow << "bad_mesh.pvd"
+            << rang::fg::reset << " file).";
+  throw NormalError(error_msg.str());
+}
+
+template void MeshData<1>::_storeBadMesh();
+template void MeshData<2>::_storeBadMesh();
+template void MeshData<3>::_storeBadMesh();
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index a448862a01391a7c51bd1e7e9d75aff9e3df996e..78b9b81e69daf5e3fc0010bd8c737b7303759353 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -5,14 +5,9 @@
 #include <mesh/IMeshData.hpp>
 #include <mesh/ItemValue.hpp>
 #include <mesh/SubItemValuePerItem.hpp>
-#include <utils/Exceptions.hpp>
 #include <utils/Messenger.hpp>
 #include <utils/PugsUtils.hpp>
 
-#include <output/VTKWriter.hpp>
-
-#include <map>
-
 template <size_t Dimension>
 class Connectivity;
 
@@ -49,6 +44,8 @@ class MeshData : public IMeshData
   FaceValue<const double> m_ll;
   EdgeValue<const Rd> m_xe;
 
+  void _storeBadMesh();
+
   PUGS_INLINE
   void
   _computeNl()
@@ -537,12 +534,7 @@ class MeshData : public IMeshData
     }();
 
     if (not parallel::allReduceAnd(is_valid)) {
-      VTKWriter writer("bad_mesh");
-      writer.writeMesh(m_mesh);
-      std::ostringstream error_msg;
-      error_msg << "mesh contains cells of non-positive volume (see " << rang::fgB::yellow << "bad_mesh.pvd"
-                << rang::fg::reset << " file).";
-      throw NormalError(error_msg.str());
+      this->_storeBadMesh();
     }
   }
 
diff --git a/src/mesh/MeshSmoother.cpp b/src/mesh/MeshSmoother.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a9af80c0129fdb20267b3c62e4430535f00f154f
--- /dev/null
+++ b/src/mesh/MeshSmoother.cpp
@@ -0,0 +1,541 @@
+#include <mesh/MeshSmoother.hpp>
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshCellZone.hpp>
+#include <mesh/MeshFlatNodeBoundary.hpp>
+#include <mesh/MeshLineNodeBoundary.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <scheme/AxisBoundaryConditionDescriptor.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/FixedBoundaryConditionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <utils/RandomEngine.hpp>
+
+#include <variant>
+
+template <size_t Dimension>
+class MeshSmootherHandler::MeshSmoother
+{
+ private:
+  using Rd               = TinyVector<Dimension>;
+  using Rdxd             = TinyMatrix<Dimension>;
+  using ConnectivityType = Connectivity<Dimension>;
+  using MeshType         = Mesh<ConnectivityType>;
+
+  const MeshType& m_given_mesh;
+
+  class AxisBoundaryCondition;
+  class FixedBoundaryCondition;
+  class SymmetryBoundaryCondition;
+
+  using BoundaryCondition = std::variant<AxisBoundaryCondition, FixedBoundaryCondition, SymmetryBoundaryCondition>;
+
+  using BoundaryConditionList = std::vector<BoundaryCondition>;
+  BoundaryConditionList m_boundary_condition_list;
+
+  BoundaryConditionList
+  _getBCList(const MeshType& mesh,
+             const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+  {
+    BoundaryConditionList bc_list;
+
+    for (const auto& bc_descriptor : bc_descriptor_list) {
+      switch (bc_descriptor->type()) {
+      case IBoundaryConditionDescriptor::Type::axis: {
+        if constexpr (Dimension == 1) {
+          bc_list.emplace_back(FixedBoundaryCondition{getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
+        } else {
+          bc_list.emplace_back(
+            AxisBoundaryCondition{getMeshLineNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
+        }
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::symmetry: {
+        bc_list.emplace_back(
+          SymmetryBoundaryCondition{getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::fixed: {
+        bc_list.emplace_back(FixedBoundaryCondition{getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
+        break;
+      }
+      default: {
+        std::ostringstream error_msg;
+        error_msg << *bc_descriptor << " is an invalid boundary condition for mesh smoother";
+        throw NormalError(error_msg.str());
+      }
+      }
+    }
+
+    return bc_list;
+  }
+
+  void
+  _applyBC(NodeValue<Rd>& shift) const
+  {
+    for (auto&& boundary_condition : m_boundary_condition_list) {
+      std::visit(
+        [&](auto&& bc) {
+          using BCType = std::decay_t<decltype(bc)>;
+          if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
+            const Rd& n = bc.outgoingNormal();
+
+            const Rdxd I   = identity;
+            const Rdxd nxn = tensorProduct(n, n);
+            const Rdxd P   = I - nxn;
+
+            const Array<const NodeId>& node_list = bc.nodeList();
+            parallel_for(
+              node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+                const NodeId node_id = node_list[i_node];
+
+                shift[node_id] = P * shift[node_id];
+              });
+
+          } else if constexpr (std::is_same_v<BCType, AxisBoundaryCondition>) {
+            if constexpr (Dimension > 1) {
+              const Rd& t = bc.direction();
+
+              const Rdxd txt = tensorProduct(t, t);
+
+              const Array<const NodeId>& node_list = bc.nodeList();
+              parallel_for(
+                node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+                  const NodeId node_id = node_list[i_node];
+
+                  shift[node_id] = txt * shift[node_id];
+                });
+            } else {
+              throw UnexpectedError("AxisBoundaryCondition make no sense in dimension 1");
+            }
+
+          } else if constexpr (std::is_same_v<BCType, FixedBoundaryCondition>) {
+            const Array<const NodeId>& node_list = bc.nodeList();
+            parallel_for(
+              node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+                const NodeId node_id = node_list[i_node];
+                shift[node_id]       = zero;
+              });
+
+          } else {
+            throw UnexpectedError("invalid boundary condition type");
+          }
+        },
+        boundary_condition);
+    }
+  }
+
+  NodeValue<Rd>
+  _getDisplacement() const
+  {
+    const ConnectivityType& connectivity = m_given_mesh.connectivity();
+    NodeValue<const Rd> given_xr         = m_given_mesh.xr();
+
+    auto node_to_cell_matrix        = connectivity.nodeToCellMatrix();
+    auto cell_to_node_matrix        = connectivity.cellToNodeMatrix();
+    auto node_number_in_their_cells = connectivity.nodeLocalNumbersInTheirCells();
+
+    NodeValue<double> max_delta_xr{connectivity};
+    parallel_for(
+      connectivity.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+        const Rd& x0 = given_xr[node_id];
+
+        const auto& node_cell_list = node_to_cell_matrix[node_id];
+        double min_distance_2      = std::numeric_limits<double>::max();
+
+        for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
+          const size_t i_cell_node = node_number_in_their_cells(node_id, i_cell);
+
+          const CellId cell_id       = node_cell_list[i_cell];
+          const auto& cell_node_list = cell_to_node_matrix[cell_id];
+
+          for (size_t i_node = 0; i_node < cell_node_list.size(); ++i_node) {
+            if (i_node != i_cell_node) {
+              const NodeId cell_node_id = cell_node_list[i_node];
+              const Rd delta            = x0 - given_xr[cell_node_id];
+              min_distance_2            = std::min(min_distance_2, dot(delta, delta));
+            }
+          }
+        }
+        double max_delta = std::sqrt(min_distance_2);
+
+        max_delta_xr[node_id] = max_delta;
+      });
+
+    NodeValue<Rd> shift_r{connectivity};
+
+    parallel_for(
+      m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+        const auto& node_cell_list = node_to_cell_matrix[node_id];
+        Rd mean_position(zero);
+        size_t number_of_neighbours = 0;
+
+        for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
+          const size_t i_cell_node = node_number_in_their_cells(node_id, i_cell);
+
+          const CellId cell_id       = node_cell_list[i_cell];
+          const auto& cell_node_list = cell_to_node_matrix[cell_id];
+          for (size_t i_node = 0; i_node < cell_node_list.size(); ++i_node) {
+            if (i_node != i_cell_node) {
+              const NodeId cell_node_id = cell_node_list[i_node];
+              mean_position += given_xr[cell_node_id];
+              number_of_neighbours++;
+            }
+          }
+        }
+        mean_position    = 1. / number_of_neighbours * mean_position;
+        shift_r[node_id] = mean_position - given_xr[node_id];
+      });
+
+    this->_applyBC(shift_r);
+
+    synchronize(shift_r);
+
+    return shift_r;
+  }
+
+ public:
+  std::shared_ptr<const IMesh>
+  getSmoothedMesh() const
+  {
+    NodeValue<const Rd> given_xr = m_given_mesh.xr();
+
+    NodeValue<Rd> xr = this->_getDisplacement();
+
+    parallel_for(
+      m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { xr[node_id] += given_xr[node_id]; });
+
+    return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+  }
+
+  std::shared_ptr<const IMesh>
+  getSmoothedMesh(const FunctionSymbolId& function_symbol_id) const
+  {
+    NodeValue<const Rd> given_xr = m_given_mesh.xr();
+
+    NodeValue<const bool> is_displaced =
+      InterpolateItemValue<bool(const Rd)>::interpolate(function_symbol_id, given_xr);
+
+    NodeValue<Rd> xr = this->_getDisplacement();
+
+    parallel_for(
+      m_given_mesh.numberOfNodes(),
+      PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id]; });
+
+    return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+  }
+
+  std::shared_ptr<const IMesh>
+  getSmoothedMesh(const std::vector<std::shared_ptr<const IZoneDescriptor>>& zone_descriptor_list) const
+  {
+    NodeValue<const Rd> given_xr = m_given_mesh.xr();
+
+    auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix();
+
+    NodeValue<bool> is_displaced{m_given_mesh.connectivity()};
+    is_displaced.fill(false);
+
+    for (size_t i_zone = 0; i_zone < zone_descriptor_list.size(); ++i_zone) {
+      MeshCellZone<Dimension> cell_zone = getMeshCellZone(m_given_mesh, *zone_descriptor_list[i_zone]);
+      const auto cell_list              = cell_zone.cellList();
+      CellValue<bool> is_zone_cell{m_given_mesh.connectivity()};
+      is_zone_cell.fill(false);
+      parallel_for(
+        cell_list.size(), PUGS_LAMBDA(const size_t i_cell) { is_zone_cell[cell_list[i_cell]] = true; });
+      parallel_for(
+        m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+          auto node_cell_list = node_to_cell_matrix[node_id];
+          bool displace       = true;
+          for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) {
+            const CellId cell_id = node_cell_list[i_node_cell];
+            displace &= is_zone_cell[cell_id];
+          }
+          if (displace) {
+            is_displaced[node_id] = true;
+          }
+        });
+    }
+    synchronize(is_displaced);
+    NodeValue<Rd> xr = this->_getDisplacement();
+
+    parallel_for(
+      m_given_mesh.numberOfNodes(),
+      PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id]; });
+
+    return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+  }
+
+  std::shared_ptr<const IMesh>
+  getSmoothedMesh(
+    const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
+  {
+    NodeValue<const Rd> given_xr = m_given_mesh.xr();
+
+    auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix();
+
+    NodeValue<bool> is_displaced{m_given_mesh.connectivity()};
+    is_displaced.fill(false);
+
+    for (size_t i_zone = 0; i_zone < discrete_function_variant_list.size(); ++i_zone) {
+      auto is_zone_cell = discrete_function_variant_list[i_zone]->get<DiscreteFunctionP0<Dimension, const double>>();
+
+      parallel_for(
+        m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+          auto node_cell_list = node_to_cell_matrix[node_id];
+          bool displace       = true;
+          for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) {
+            const CellId cell_id = node_cell_list[i_node_cell];
+            displace &= (is_zone_cell[cell_id] != 0);
+          }
+          if (displace) {
+            is_displaced[node_id] = true;
+          }
+        });
+    }
+    synchronize(is_displaced);
+    NodeValue<Rd> xr = this->_getDisplacement();
+
+    parallel_for(
+      m_given_mesh.numberOfNodes(),
+      PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id]; });
+
+    return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+  }
+
+  MeshSmoother(const MeshSmoother&) = delete;
+  MeshSmoother(MeshSmoother&&)      = delete;
+
+  MeshSmoother(const MeshType& given_mesh,
+               const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+    : m_given_mesh(given_mesh), m_boundary_condition_list(this->_getBCList(given_mesh, bc_descriptor_list))
+  {}
+
+  ~MeshSmoother() = default;
+};
+
+template <size_t Dimension>
+class MeshSmootherHandler::MeshSmoother<Dimension>::AxisBoundaryCondition
+{
+ public:
+  using Rd = TinyVector<Dimension, double>;
+
+ private:
+  const MeshLineNodeBoundary<Dimension> m_mesh_line_node_boundary;
+
+ public:
+  const Rd&
+  direction() const
+  {
+    return m_mesh_line_node_boundary.direction();
+  }
+
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_line_node_boundary.nodeList();
+  }
+
+  AxisBoundaryCondition(MeshLineNodeBoundary<Dimension>&& mesh_line_node_boundary)
+    : m_mesh_line_node_boundary(mesh_line_node_boundary)
+  {
+    ;
+  }
+
+  ~AxisBoundaryCondition() = default;
+};
+
+template <>
+class MeshSmootherHandler::MeshSmoother<1>::AxisBoundaryCondition
+{
+ public:
+  AxisBoundaryCondition()  = default;
+  ~AxisBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class MeshSmootherHandler::MeshSmoother<Dimension>::FixedBoundaryCondition
+{
+ private:
+  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+
+ public:
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_node_boundary.nodeList();
+  }
+
+  FixedBoundaryCondition(MeshNodeBoundary<Dimension>&& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}
+
+  ~FixedBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class MeshSmootherHandler::MeshSmoother<Dimension>::SymmetryBoundaryCondition
+{
+ public:
+  using Rd = TinyVector<Dimension, double>;
+
+ private:
+  const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
+
+ public:
+  const Rd&
+  outgoingNormal() const
+  {
+    return m_mesh_flat_node_boundary.outgoingNormal();
+  }
+
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_flat_node_boundary.nodeList();
+  }
+
+  SymmetryBoundaryCondition(MeshFlatNodeBoundary<Dimension>&& mesh_flat_node_boundary)
+    : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
+  {
+    ;
+  }
+
+  ~SymmetryBoundaryCondition() = default;
+};
+
+std::shared_ptr<const IMesh>
+MeshSmootherHandler::getSmoothedMesh(
+  const std::shared_ptr<const IMesh>& mesh,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+{
+  switch (mesh->dimension()) {
+  case 1: {
+    constexpr size_t Dimension = 1;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh();
+  }
+  case 2: {
+    constexpr size_t Dimension = 2;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh();
+  }
+  case 3: {
+    constexpr size_t Dimension = 3;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh();
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+std::shared_ptr<const IMesh>
+MeshSmootherHandler::getSmoothedMesh(
+  const std::shared_ptr<const IMesh>& mesh,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+  const FunctionSymbolId& function_symbol_id) const
+{
+  switch (mesh->dimension()) {
+  case 1: {
+    constexpr size_t Dimension = 1;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh(function_symbol_id);
+  }
+  case 2: {
+    constexpr size_t Dimension = 2;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh(function_symbol_id);
+  }
+  case 3: {
+    constexpr size_t Dimension = 3;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh(function_symbol_id);
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+std::shared_ptr<const IMesh>
+MeshSmootherHandler::getSmoothedMesh(
+  const std::shared_ptr<const IMesh>& mesh,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+  const std::vector<std::shared_ptr<const IZoneDescriptor>>& smoothing_zone_list) const
+{
+  switch (mesh->dimension()) {
+  case 1: {
+    constexpr size_t Dimension = 1;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh(smoothing_zone_list);
+  }
+  case 2: {
+    constexpr size_t Dimension = 2;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh(smoothing_zone_list);
+  }
+  case 3: {
+    constexpr size_t Dimension = 3;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh(smoothing_zone_list);
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+std::shared_ptr<const IMesh>
+MeshSmootherHandler::getSmoothedMesh(
+  const std::shared_ptr<const IMesh>& mesh,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+  const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
+{
+  if (not hasSameMesh(discrete_function_variant_list)) {
+    throw NormalError("discrete functions are not defined on the same mesh");
+  }
+
+  std::shared_ptr<const IMesh> common_mesh = getCommonMesh(discrete_function_variant_list);
+
+  if (common_mesh != mesh) {
+    throw NormalError("discrete functions are not defined on the smoothed mesh");
+  }
+
+  switch (mesh->dimension()) {
+  case 1: {
+    constexpr size_t Dimension = 1;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh(discrete_function_variant_list);
+  }
+  case 2: {
+    constexpr size_t Dimension = 2;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh(discrete_function_variant_list);
+  }
+  case 3: {
+    constexpr size_t Dimension = 3;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshSmoother smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+    return smoother.getSmoothedMesh(discrete_function_variant_list);
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
diff --git a/src/mesh/MeshSmoother.hpp b/src/mesh/MeshSmoother.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..99db9bf76bf9ea771b1ebbc41c1b948c6d451632
--- /dev/null
+++ b/src/mesh/MeshSmoother.hpp
@@ -0,0 +1,45 @@
+#ifndef MESH_SMOOTHER_HPP
+#define MESH_SMOOTHER_HPP
+
+#include <mesh/IMesh.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+
+#include <memory>
+#include <vector>
+
+class FunctionSymbolId;
+class IZoneDescriptor;
+class DiscreteFunctionVariant;
+
+class MeshSmootherHandler
+{
+ private:
+  template <size_t Dimension>
+  class MeshSmoother;
+
+ public:
+  std::shared_ptr<const IMesh> getSmoothedMesh(
+    const std::shared_ptr<const IMesh>& mesh,
+    const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const;
+
+  std::shared_ptr<const IMesh> getSmoothedMesh(
+    const std::shared_ptr<const IMesh>& mesh,
+    const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+    const FunctionSymbolId& function_symbol_id) const;
+
+  std::shared_ptr<const IMesh> getSmoothedMesh(
+    const std::shared_ptr<const IMesh>& mesh,
+    const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+    const std::vector<std::shared_ptr<const IZoneDescriptor>>& smoothing_zone_list) const;
+
+  std::shared_ptr<const IMesh> getSmoothedMesh(
+    const std::shared_ptr<const IMesh>& mesh,
+    const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+    const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& smoothing_zone_list) const;
+
+  MeshSmootherHandler()                      = default;
+  MeshSmootherHandler(MeshSmootherHandler&&) = default;
+  ~MeshSmootherHandler()                     = default;
+};
+
+#endif   // MESH_RANDOMIZER_HPP
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index 349f5cfd968c63b98fd6b010bd43dedf2caaafcf..1379aea67566c559a17f19a5315335b4248e74c2 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -8,4 +8,5 @@ add_library(
   DiscreteFunctionInterpoler.cpp
   DiscreteFunctionUtils.cpp
   DiscreteFunctionVectorIntegrator.cpp
-  DiscreteFunctionVectorInterpoler.cpp)
+  DiscreteFunctionVectorInterpoler.cpp
+  FluxingAdvectionSolver.cpp)
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
new file mode 100644
index 0000000000000000000000000000000000000000..264ecef042956e6b77230e00f49cd321b50dcc6a
--- /dev/null
+++ b/src/scheme/FluxingAdvectionSolver.cpp
@@ -0,0 +1,891 @@
+#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>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/Mesh.hpp>
+#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
+{
+ private:
+  using Rd = TinyVector<Dimension>;
+
+  using MeshType     = Mesh<Connectivity<Dimension>>;
+  using MeshDataType = MeshData<Dimension>;
+
+  const std::shared_ptr<const MeshType> m_old_mesh;
+  const std::shared_ptr<const MeshType> m_new_mesh;
+
+  using RemapVariant = std::variant<CellValue<double>,
+                                    CellValue<TinyVector<1>>,
+                                    CellValue<TinyVector<2>>,
+                                    CellValue<TinyVector<3>>,
+                                    CellValue<TinyMatrix<1>>,
+                                    CellValue<TinyMatrix<2>>,
+                                    CellValue<TinyMatrix<3>>,
+
+                                    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;
+  size_t m_number_of_cycles;
+
+  FaceValue<double> _computeAlgebraicFluxingVolume();
+  void _computeDonorCells(FaceValue<const double> algebraic_fluxing_volumes);
+  FaceValue<double> _computeFluxingVolume(FaceValue<double> algebraic_fluxing_volumes);
+  void _computeCycleNumber(FaceValue<double> fluxing_volumes);
+  void _computeGeometricalData();
+
+  template <typename DataType>
+  void
+  _storeValues(const DiscreteFunctionP0<Dimension, const DataType>& old_q)
+  {
+    m_remapped_list.emplace_back(copy(old_q.cellValues()));
+  }
+
+  void
+  _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,
+                 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 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)},
+      m_new_mesh{std::dynamic_pointer_cast<const MeshType>(i_new_mesh)}
+  {
+    if ((m_old_mesh.use_count() == 0) or (m_new_mesh.use_count() == 0)) {
+      throw NormalError("old and new meshes must be of same type");
+    }
+
+    if (m_new_mesh->shared_connectivity() != m_old_mesh->shared_connectivity()) {
+      throw NormalError("old and new meshes must share the same connectivity");
+    }
+
+    this->_computeGeometricalData();
+  }
+
+  ~FluxingAdvectionSolver() = default;
+};
+
+template <size_t Dimension>
+void
+FluxingAdvectionSolver<Dimension>::_computeDonorCells(FaceValue<const double> algebraic_fluxing_volumes)
+{
+  m_donnor_cell = [&] {
+    const FaceValuePerCell<const bool> cell_face_is_reversed = m_new_mesh->connectivity().cellFaceIsReversed();
+    const auto face_to_cell_matrix                           = m_new_mesh->connectivity().faceToCellMatrix();
+
+    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 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 {
+          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();
+          }
+        }
+      });
+
+    return donnor_cell;
+  }();
+}
+
+template <>
+void
+FluxingAdvectionSolver<1>::_computeDonorCells(FaceValue<const double> algebraic_fluxing_volumes)
+{
+  m_donnor_cell = [&] {
+    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) {
+          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 {
+          if ((algebraic_fluxing_volumes[face_id] <= 0) xor (cell_to_face_matrix[cell_id][0] == face_id)) {
+            donnor_cell[face_id] = cell_id;
+          } else {
+            donnor_cell[face_id] = face_to_cell[1];
+          }
+        }
+      });
+
+    return donnor_cell;
+  }();
+}
+
+template <>
+FaceValue<double>
+FluxingAdvectionSolver<1>::_computeAlgebraicFluxingVolume()
+{
+  Array<double> fluxing_volumes{m_new_mesh->numberOfNodes()};
+  NodeValue<double> nodal_fluxing_volume(m_new_mesh->connectivity(), fluxing_volumes);
+  auto old_xr = m_old_mesh->xr();
+  auto new_xr = m_new_mesh->xr();
+
+  parallel_for(
+    m_new_mesh->numberOfNodes(),
+    PUGS_LAMBDA(NodeId node_id) { nodal_fluxing_volume[node_id] = new_xr[node_id][0] - old_xr[node_id][0]; });
+
+  FaceValue<double> algebraic_fluxing_volumes(m_new_mesh->connectivity(), fluxing_volumes);
+
+  synchronize(algebraic_fluxing_volumes);
+  return algebraic_fluxing_volumes;
+}
+
+template <>
+FaceValue<double>
+FluxingAdvectionSolver<2>::_computeAlgebraicFluxingVolume()
+{
+  const auto face_to_node_matrix = m_old_mesh->connectivity().faceToNodeMatrix();
+  FaceValue<double> algebraic_fluxing_volume(m_new_mesh->connectivity());
+  auto old_xr = m_old_mesh->xr();
+  auto new_xr = m_new_mesh->xr();
+  parallel_for(
+    m_new_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+      const auto& face_to_node = face_to_node_matrix[face_id];
+
+      const Rd& x0 = old_xr[face_to_node[0]];
+      const Rd& x1 = old_xr[face_to_node[1]];
+      const Rd& x2 = new_xr[face_to_node[1]];
+      const Rd& x3 = new_xr[face_to_node[0]];
+
+      TinyMatrix<2> M(x3[0] - x1[0], x2[0] - x0[0],   //
+                      x3[1] - x1[1], x2[1] - x0[1]);
+
+      algebraic_fluxing_volume[face_id] = 0.5 * det(M);
+    });
+
+  synchronize(algebraic_fluxing_volume);
+  return algebraic_fluxing_volume;
+}
+
+template <>
+FaceValue<double>
+FluxingAdvectionSolver<3>::_computeAlgebraicFluxingVolume()
+{
+  const auto face_to_node_matrix = m_old_mesh->connectivity().faceToNodeMatrix();
+  FaceValue<double> algebraic_fluxing_volume(m_new_mesh->connectivity());
+  auto old_xr = m_old_mesh->xr();
+  auto new_xr = m_new_mesh->xr();
+  parallel_for(
+    m_new_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+      const auto& face_to_node = face_to_node_matrix[face_id];
+      if (face_to_node.size() == 4) {
+        const Rd& x0 = old_xr[face_to_node[0]];
+        const Rd& x1 = old_xr[face_to_node[1]];
+        const Rd& x2 = old_xr[face_to_node[2]];
+        const Rd& x3 = old_xr[face_to_node[3]];
+
+        const Rd& x4 = new_xr[face_to_node[0]];
+        const Rd& x5 = new_xr[face_to_node[1]];
+        const Rd& x6 = new_xr[face_to_node[2]];
+        const Rd& x7 = new_xr[face_to_node[3]];
+
+        const Rd& a1 = x6 - x1;
+        const Rd& a2 = x6 - x3;
+        const Rd& a3 = x6 - x4;
+
+        const Rd& b1 = x7 - x0;
+        const Rd& b2 = x5 - x0;
+        const Rd& b3 = x2 - x0;
+
+        TinyMatrix<3> M1(a1 + b1, a2, a3);
+        TinyMatrix<3> M2(b1, a2 + b2, a3);
+        TinyMatrix<3> M3(a1, b2, a3 + b3);
+
+        algebraic_fluxing_volume[face_id] = (det(M1) + det(M2) + det(M3)) / 12;
+      } else if (face_to_node.size() == 3) {
+        const Rd& x0 = old_xr[face_to_node[0]];
+        const Rd& x1 = old_xr[face_to_node[1]];
+        const Rd& x2 = old_xr[face_to_node[2]];
+
+        const Rd& x3 = new_xr[face_to_node[0]];
+        const Rd& x4 = new_xr[face_to_node[1]];
+        const Rd& x5 = new_xr[face_to_node[2]];
+
+        const Rd& a1 = x5 - x1;
+        const Rd& a2 = x5 - x2;
+        const Rd& a3 = x5 - x3;
+
+        const Rd& b1 = x5 - x0;
+        const Rd& b2 = x4 - x0;
+        const Rd& b3 = x2 - x0;
+
+        TinyMatrix<3> M1(a1 + b1, a2, a3);
+        TinyMatrix<3> M2(b1, a2 + b2, a3);
+        TinyMatrix<3> M3(a1, b2, a3 + b3);
+
+        algebraic_fluxing_volume[face_id] = (det(M1) + det(M2) + det(M3)) / 12;
+      } else {
+        throw NotImplementedError("Not implemented for non quad faces");
+      }
+    });
+
+  return algebraic_fluxing_volume;
+}
+
+template <size_t Dimension>
+FaceValue<double>
+FluxingAdvectionSolver<Dimension>::_computeFluxingVolume(FaceValue<double> algebraic_fluxing_volumes)
+{
+  Assert(m_donnor_cell.isBuilt());
+  // Now that donnor cells are clearly defined, we consider the
+  // non-algebraic volumes of fluxing
+  parallel_for(
+    algebraic_fluxing_volumes.numberOfItems(), PUGS_LAMBDA(const FaceId face_id) {
+      algebraic_fluxing_volumes[face_id] = std::abs(algebraic_fluxing_volumes[face_id]);
+    });
+
+  return algebraic_fluxing_volumes;
+}
+
+template <size_t Dimension>
+void
+FluxingAdvectionSolver<Dimension>::_computeCycleNumber(FaceValue<double> fluxing_volumes)
+{
+  const auto cell_to_face_matrix = m_old_mesh->connectivity().cellToFaceMatrix();
+
+  const CellValue<double> total_negative_flux(m_old_mesh->connectivity());
+  total_negative_flux.fill(0);
+
+  parallel_for(
+    m_old_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+      const auto& cell_to_face = cell_to_face_matrix[cell_id];
+      for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+        FaceId face_id = cell_to_face[i_face];
+        if (cell_id == m_donnor_cell[face_id]) {
+          total_negative_flux[cell_id] += fluxing_volumes[face_id];
+        }
+      }
+    });
+
+  MeshData<Dimension>& mesh_data   = MeshDataManager::instance().getMeshData(*m_old_mesh);
+  const CellValue<const double> Vj = mesh_data.Vj();
+  CellValue<size_t> ratio(m_old_mesh->connectivity());
+
+  parallel_for(
+    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);
+
+  if (number_of_cycles > 1) {
+    const double cycle_ratio = 1. / number_of_cycles;
+
+    parallel_for(
+      fluxing_volumes.numberOfItems(), PUGS_LAMBDA(const FaceId face_id) { fluxing_volumes[face_id] *= cycle_ratio; });
+  }
+
+  m_number_of_cycles     = number_of_cycles;
+  m_cycle_fluxing_volume = fluxing_volumes;
+}
+
+template <size_t Dimension>
+void
+FluxingAdvectionSolver<Dimension>::_computeGeometricalData()
+{
+  auto fluxing_volumes = this->_computeAlgebraicFluxingVolume();
+  this->_computeDonorCells(fluxing_volumes);
+  fluxing_volumes = this->_computeFluxingVolume(fluxing_volumes);
+  this->_computeCycleNumber(fluxing_volumes);
+}
+
+template <size_t Dimension>
+template <typename CellDataType>
+void
+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();
+  const auto face_to_cell_matrix = m_new_mesh->connectivity().faceToCellMatrix();
+
+  auto new_q = copy(old_q);
+
+  if constexpr (is_item_value_v<CellDataType>) {
+    parallel_for(
+      m_new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { new_q[cell_id] *= step_Vj[cell_id]; });
+  } else if constexpr (is_item_array_v<CellDataType>) {
+    parallel_for(
+      m_new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        auto new_array      = new_q[cell_id];
+        const double volume = step_Vj[cell_id];
+
+        for (size_t i = 0; i < new_array.size(); ++i) {
+          new_array[i] *= volume;
+        }
+      });
+  }
+
+  // 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];
+      for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+        const FaceId face_id        = cell_to_face[i_face];
+        const double fluxing_volume = m_cycle_fluxing_volume[face_id];
+
+        const auto& face_to_cell = face_to_cell_matrix[face_id];
+
+        if (face_to_cell.size() == 1) {
+          continue;
+        }
+
+        CellId donnor_id = m_donnor_cell[face_id];
+
+        if constexpr (is_item_value_v<CellDataType>) {
+          auto fluxed_q = old_q[donnor_id];
+          fluxed_q *= ((donnor_id == cell_id) ? -1 : 1) * fluxing_volume;
+
+          new_q[cell_id] += fluxed_q;
+        } else if constexpr (is_item_array_v<CellDataType>) {
+          const double sign   = ((donnor_id == cell_id) ? -1 : 1);
+          auto old_cell_array = old_q[donnor_id];
+          auto new_cell_array = new_q[cell_id];
+          for (size_t i = 0; i < new_cell_array.size(); ++i) {
+            new_cell_array[i] += (sign * fluxing_volume) * old_cell_array[i];
+          }
+        }
+      }
+    });
+
+  // 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;
+}
+
+template <size_t Dimension>
+void
+FluxingAdvectionSolver<Dimension>::_remapAllQuantities()
+{
+  const auto cell_to_face_matrix              = m_new_mesh->connectivity().cellToFaceMatrix();
+  const auto face_local_number_in_their_cells = m_new_mesh->connectivity().faceLocalNumbersInTheirCells();
+
+  MeshData<Dimension>& old_mesh_data = MeshDataManager::instance().getMeshData(*m_old_mesh);
+
+  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 (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(
+      m_new_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        const auto& cell_to_face = cell_to_face_matrix[cell_id];
+        for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+          const FaceId face_id = cell_to_face[i_face];
+          CellId donnor_id     = m_donnor_cell[face_id];
+
+          double flux = ((donnor_id == cell_id) ? -1 : 1) * m_cycle_fluxing_volume[face_id];
+          step_Vj[cell_id] += flux;
+        }
+      });
+
+    synchronize(step_Vj);
+
+    CellValue<double> inv_Vj(m_old_mesh->connectivity());
+    parallel_for(
+      m_new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { inv_Vj[cell_id] = 1 / step_Vj[cell_id]; });
+
+    for (auto& remapped_q : m_remapped_list) {
+      std::visit(
+        [&](auto&& new_q) {
+          using CellDataType = std::decay_t<decltype(new_q)>;
+          static_assert(is_item_value_v<CellDataType> or is_item_array_v<CellDataType>, "invalid data type");
+
+          if constexpr (is_item_value_v<CellDataType>) {
+            parallel_for(
+              m_new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { new_q[cell_id] *= inv_Vj[cell_id]; });
+          } else if constexpr (is_item_array_v<CellDataType>) {
+            parallel_for(
+              m_new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                auto array              = new_q[cell_id];
+                const double inv_volume = inv_Vj[cell_id];
+
+                for (size_t i = 0; i < array.size(); ++i) {
+                  array[i] *= inv_volume;
+                }
+              });
+          }
+        },
+        remapped_q);
+    }
+  }
+}
+
+template <size_t Dimension>
+std::vector<std::shared_ptr<const DiscreteFunctionVariant>>
+FluxingAdvectionSolver<Dimension>::remap(
+  const std::vector<std::shared_ptr<const VariableBCDescriptor>>& quantity_list_with_bc)
+{
+  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");
+        }
+      },
+      quantity_v->discreteFunction());
+  }
+
+  this->_remapAllQuantities();
+
+  std::vector<std::shared_ptr<const DiscreteFunctionVariant>> new_variables;
+
+  for (size_t i = 0; i < quantity_list_with_bc.size(); ++i) {
+    std::visit(
+      [&](auto&& variable) {
+        using DiscreteFunctionT = std::decay_t<decltype(variable)>;
+        using DataType          = std::decay_t<typename DiscreteFunctionT::data_type>;
+
+        if constexpr (std::is_same_v<MeshType, typename DiscreteFunctionT::MeshType>) {
+          if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
+            new_variables.push_back(std::make_shared<DiscreteFunctionVariant>(
+              DiscreteFunctionT(m_new_mesh, std::get<CellValue<DataType>>(m_remapped_list[i]))));
+          } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
+            new_variables.push_back(std::make_shared<DiscreteFunctionVariant>(
+              DiscreteFunctionT(m_new_mesh, std::get<CellArray<DataType>>(m_remapped_list[i]))));
+          } else {
+            throw UnexpectedError("invalid discrete function type");
+          }
+        } else {
+          throw UnexpectedError("incompatible mesh types");
+        }
+      },
+      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 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");
+  }
+
+  const std::shared_ptr<const IMesh> i_old_mesh = getCommonMesh(remapped_variables);
+
+  switch (i_old_mesh->dimension()) {
+  case 1: {
+    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_with_bc);
+  }
+  case 3: {
+    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
new file mode 100644
index 0000000000000000000000000000000000000000..2cbb653cbeed76000d8eb1e3ba01ed61c35bc8a3
--- /dev/null
+++ b/src/scheme/FluxingAdvectionSolver.hpp
@@ -0,0 +1,18 @@
+#ifndef FLUXING_ADVECION_SOLVER_HPP
+#define FLUXING_ADVECION_SOLVER_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/VariableBCDescriptor.hpp>
+
+#include <vector>
+
+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");
+        }
       }
     }
   }
diff --git a/tests/test_TinyMatrix.cpp b/tests/test_TinyMatrix.cpp
index f6f09cf55dd2ad1ffef55c9275b989c203598f7e..03d3bbe5dc0ad69655cda5d3ae764d8e770b0d60 100644
--- a/tests/test_TinyMatrix.cpp
+++ b/tests/test_TinyMatrix.cpp
@@ -32,10 +32,14 @@ TEST_CASE("TinyMatrix", "[algebra]")
   REQUIRE(TinyMatrix<5, 4, int>::NumberOfColumns == 4);
 
   TinyMatrix<3, 4, int> A(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12);
+  TinyVector<3, int> x1(1, 5, 9);
+  TinyVector<3, int> x2(2, 6, 10);
+  TinyVector<3, int> x3(3, 7, 11);
+  TinyVector<3, int> x4(4, 8, 12);
   REQUIRE(((A(0, 0) == 1) and (A(0, 1) == 2) and (A(0, 2) == 3) and (A(0, 3) == 4) and   //
            (A(1, 0) == 5) and (A(1, 1) == 6) and (A(1, 2) == 7) and (A(1, 3) == 8) and   //
            (A(2, 0) == 9) and (A(2, 1) == 10) and (A(2, 2) == 11) and (A(2, 3) == 12)));
-
+  REQUIRE(A == TinyMatrix<3, 4, int>(x1, x2, x3, x4));
   TinyMatrix<3, 4, int> B(6, 5, 3, 8, 34, 6, 35, 6, 7, 1, 3, 6);
 
   SECTION("checking for opposed matrix")