diff --git a/src/language/ast/ASTNodeBuiltinFunctionExpressionBuilder.cpp b/src/language/ast/ASTNodeBuiltinFunctionExpressionBuilder.cpp
index 487259530b7c3e4872371e9623517c0546d31e4f..412906b9d67291ed7901d3ac00943cc8a9107449 100644
--- a/src/language/ast/ASTNodeBuiltinFunctionExpressionBuilder.cpp
+++ b/src/language/ast/ASTNodeBuiltinFunctionExpressionBuilder.cpp
@@ -231,6 +231,9 @@ ASTNodeBuiltinFunctionExpressionBuilder::_getArgumentConverter(const ASTNodeData
       }
       case ASTNodeDataType::double_t: {
         return std::make_unique<FunctionTupleArgumentConverter<ParameterContentT, double>>(argument_number);
+      }
+      case ASTNodeDataType::function_t: {
+        return std::make_unique<FunctionTupleArgumentConverter<ParameterContentT, FunctionSymbolId>>(argument_number);
       }
         // LCOV_EXCL_START
       default: {
@@ -240,7 +243,15 @@ ASTNodeBuiltinFunctionExpressionBuilder::_getArgumentConverter(const ASTNodeData
       }
     }
     case ASTNodeDataType::list_t: {
-      return std::make_unique<FunctionListArgumentConverter<ParameterContentT, ParameterContentT>>(argument_number);
+      if constexpr (std::is_same_v<ParameterContentT, FunctionSymbolId>) {
+        const ASTNode& parent_node = argument_node_sub_data_type.m_parent_node;
+        auto symbol_table          = parent_node.m_symbol_table;
+
+        return std::make_unique<FunctionListArgumentConverter<FunctionSymbolId, FunctionSymbolId>>(argument_number,
+                                                                                                   symbol_table);
+      } else {
+        return std::make_unique<FunctionListArgumentConverter<ParameterContentT, ParameterContentT>>(argument_number);
+      }
     }
     case ASTNodeDataType::type_id_t: {
       return std::make_unique<FunctionTupleArgumentConverter<ParameterContentT, EmbeddedData>>(argument_number);
@@ -257,6 +268,12 @@ ASTNodeBuiltinFunctionExpressionBuilder::_getArgumentConverter(const ASTNodeData
     case ASTNodeDataType::double_t: {
       return std::make_unique<FunctionTupleArgumentConverter<ParameterContentT, double>>(argument_number);
     }
+    case ASTNodeDataType::function_t: {
+      const ASTNode& parent_node = argument_node_sub_data_type.m_parent_node;
+      auto symbol_table          = parent_node.m_symbol_table;
+
+      return std::make_unique<FunctionArgumentToTupleFunctionSymbolIdConverter>(argument_number, symbol_table);
+    }
     case ASTNodeDataType::vector_t: {
       switch (arg_data_type.dimension()) {
       case 1: {
@@ -398,6 +415,9 @@ ASTNodeBuiltinFunctionExpressionBuilder::_getArgumentConverter(const ASTNodeData
       case ASTNodeDataType::double_t: {
         return get_function_argument_to_tuple_converter(double{});
       }
+      case ASTNodeDataType::function_t: {
+        return get_function_argument_to_tuple_converter(FunctionSymbolId{});
+      }
       case ASTNodeDataType::vector_t: {
         switch (parameter_type.contentType().dimension()) {
         case 1: {
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 4d1a866326a46c63f71b9b45b5188d7aa11746ca..7c0d55cf1f23e9f5b186f220fbe2608e2ada7b7f 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -11,6 +11,7 @@
 #include <scheme/DiscreteFunctionDescriptorP0.hpp>
 #include <scheme/DiscreteFunctionInterpoler.hpp>
 #include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/DiscreteFunctionVectorInterpoler.hpp>
 #include <scheme/FourierBoundaryConditionDescriptor.hpp>
 #include <scheme/FreeBoundaryConditionDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
@@ -42,6 +43,19 @@ SchemeModule::SchemeModule()
 
                                     ));
 
+  this->_addBuiltinFunction(
+    "interpolate",
+    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+      const IDiscreteFunction>(std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunctionDescriptor>,
+                               const std::vector<FunctionSymbolId>&)>>(
+      [](std::shared_ptr<const IMesh> mesh,
+         std::shared_ptr<const IDiscreteFunctionDescriptor> discrete_function_descriptor,
+         const std::vector<FunctionSymbolId>& function_id_list) -> std::shared_ptr<const IDiscreteFunction> {
+        return DiscreteFunctionVectorInterpoler{mesh, discrete_function_descriptor, function_id_list}.interpolate();
+      }
+
+      ));
+
   this->_addBuiltinFunction(
     "interpolate",
     std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
diff --git a/src/language/node_processor/FunctionArgumentConverter.hpp b/src/language/node_processor/FunctionArgumentConverter.hpp
index af08e7806a4fd5daa01a7aecc98a84e52291040f..b13eb222b08cd2c0f5dde5e72d4de5626ef7758b 100644
--- a/src/language/node_processor/FunctionArgumentConverter.hpp
+++ b/src/language/node_processor/FunctionArgumentConverter.hpp
@@ -242,6 +242,7 @@ class FunctionListArgumentConverter final : public IFunctionArgumentConverter
   DataVariant
   convert(ExecutionPolicy& exec_policy, DataVariant&& value)
   {
+    static_assert(not std::is_same_v<ContentType, FunctionSymbolId>);
     using TupleType = std::vector<ContentType>;
     if constexpr (std::is_same_v<ContentType, ProvidedValueType>) {
       std::visit(
@@ -304,6 +305,57 @@ class FunctionListArgumentConverter final : public IFunctionArgumentConverter
   FunctionListArgumentConverter(size_t argument_id) : m_argument_id{argument_id} {}
 };
 
+template <>
+class FunctionListArgumentConverter<FunctionSymbolId, FunctionSymbolId> final : public IFunctionArgumentConverter
+{
+ private:
+  size_t m_argument_id;
+  std::shared_ptr<SymbolTable> m_symbol_table;
+
+ public:
+  DataVariant
+  convert(ExecutionPolicy& exec_policy, DataVariant&& value)
+  {
+    using TupleType = std::vector<FunctionSymbolId>;
+    std::visit(
+      [&](auto&& v) {
+        using ValueT = std::decay_t<decltype(v)>;
+        if constexpr (std::is_same_v<ValueT, AggregateDataVariant>) {
+          TupleType list_value;
+          list_value.reserve(v.size());
+          for (size_t i = 0; i < v.size(); ++i) {
+            std::visit(
+              [&](auto&& vi) {
+                using Vi_T = std::decay_t<decltype(vi)>;
+                if constexpr (std::is_same_v<Vi_T, uint64_t>) {
+                  list_value.emplace_back(FunctionSymbolId{vi, m_symbol_table});
+                } else {
+                  // LCOV_EXCL_START
+                  throw UnexpectedError(std::string{"invalid conversion of '"} + demangle<Vi_T>() + "' to '" +
+                                        demangle<FunctionSymbolId>() + "'");
+                  // LCOV_EXCL_STOP
+                }
+              },
+              (v[i]));
+          }
+          exec_policy.currentContext()[m_argument_id] = std::move(list_value);
+        } else {
+          // LCOV_EXCL_START
+          throw UnexpectedError(std::string{"invalid conversion of '"} + demangle<ValueT>() + "' to '" +
+                                demangle<FunctionSymbolId>() + "' list");
+          // LCOV_EXCL_STOP
+        }
+      },
+      value);
+
+    return {};
+  }
+
+  FunctionListArgumentConverter(size_t argument_id, const std::shared_ptr<SymbolTable>& symbol_table)
+    : m_argument_id{argument_id}, m_symbol_table{symbol_table}
+  {}
+};
+
 class FunctionArgumentToFunctionSymbolIdConverter final : public IFunctionArgumentConverter
 {
  private:
@@ -324,4 +376,25 @@ class FunctionArgumentToFunctionSymbolIdConverter final : public IFunctionArgume
   {}
 };
 
+class FunctionArgumentToTupleFunctionSymbolIdConverter final : public IFunctionArgumentConverter
+{
+ private:
+  size_t m_argument_id;
+  std::shared_ptr<SymbolTable> m_symbol_table;
+
+ public:
+  DataVariant
+  convert(ExecutionPolicy& exec_policy, DataVariant&& value)
+  {
+    exec_policy.currentContext()[m_argument_id] =
+      std::vector{FunctionSymbolId{std::get<uint64_t>(value), m_symbol_table}};
+
+    return {};
+  }
+
+  FunctionArgumentToTupleFunctionSymbolIdConverter(size_t argument_id, const std::shared_ptr<SymbolTable>& symbol_table)
+    : m_argument_id{argument_id}, m_symbol_table{symbol_table}
+  {}
+};
+
 #endif   // FUNCTION_ARGUMENT_CONVERTER_HPP
diff --git a/src/language/utils/BuiltinFunctionEmbedderUtils.cpp b/src/language/utils/BuiltinFunctionEmbedderUtils.cpp
index 9899fee3e1855a2ded4e306dd5dcf041d84e051b..ffe1112d0543fa644aa952a9599d9fb12eb5e383 100644
--- a/src/language/utils/BuiltinFunctionEmbedderUtils.cpp
+++ b/src/language/utils/BuiltinFunctionEmbedderUtils.cpp
@@ -233,12 +233,20 @@ getBuiltinFunctionEmbedder(ASTNode& n)
       return callable_id_list[0];
     }
     default: {
+      auto& builtin_function_embedder_table = n.m_symbol_table->builtinFunctionEmbedderTable();
+      for (auto callable_id : callable_id_list) {
+        std::shared_ptr builtin_function_embedder = builtin_function_embedder_table[callable_id];
+        // If has exact match return it
+        if (dataTypeName(args_node.m_data_type) == dataTypeName(builtin_function_embedder->getParameterDataTypes())) {
+          return callable_id;
+        }
+      }
+      // else this is an ambiguous call
       std::ostringstream error_msg;
       error_msg << "ambiguous function call " << rang::fgB::red << builtin_function_name << rang::style::reset
                 << rang::style::bold << ": " << rang::fgB::yellow << dataTypeName(args_node.m_data_type)
                 << rang::style::reset << rang::style::bold << "\nnote: candidates are";
 
-      auto& builtin_function_embedder_table = n.m_symbol_table->builtinFunctionEmbedderTable();
       for (auto callable_id : callable_id_list) {
         std::shared_ptr builtin_function_embedder = builtin_function_embedder_table[callable_id];
 
diff --git a/src/language/utils/DataVariant.hpp b/src/language/utils/DataVariant.hpp
index 6add9e52853bcbf58ee42e4f2869649dab43d37d..964044c6121478b4247c089a45eaf3ebe025195e 100644
--- a/src/language/utils/DataVariant.hpp
+++ b/src/language/utils/DataVariant.hpp
@@ -40,7 +40,8 @@ using DataVariant = std::variant<std::monostate,
                                  std::vector<TinyMatrix<3>>,
                                  std::vector<EmbeddedData>,
                                  AggregateDataVariant,
-                                 FunctionSymbolId>;
+                                 FunctionSymbolId,
+                                 std::vector<FunctionSymbolId>>;
 
 template <typename T, typename...>
 inline constexpr bool is_data_variant_v = is_variant<T, DataVariant>::value;
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionOperators.cpp b/src/language/utils/EmbeddedIDiscreteFunctionOperators.cpp
index db7a609b9449d27380c2030c75e938554e6a095b..af0ac957a430cecc6bcdd8b643064816a2402b94 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionOperators.cpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionOperators.cpp
@@ -2,6 +2,7 @@
 
 #include <language/node_processor/BinaryExpressionProcessor.hpp>
 #include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
 #include <scheme/IDiscreteFunction.hpp>
 #include <utils/Exceptions.hpp>
 
@@ -16,14 +17,29 @@ template <>
 PUGS_INLINE std::string
 name(const IDiscreteFunction& f)
 {
-  return "Vh(" + dataTypeName(f.dataType()) + ")";
+  return "Vh(" + name(f.descriptor().type()) + ":" + dataTypeName(f.dataType()) + ")";
 }
 
 template <>
 PUGS_INLINE std::string
 name(const std::shared_ptr<const IDiscreteFunction>& f)
 {
-  return "Vh(" + dataTypeName(f->dataType()) + ")";
+  return "Vh(" + name(f->descriptor().type()) + ":" + dataTypeName(f->dataType()) + ")";
+}
+
+PUGS_INLINE
+bool
+isSameDiscretization(const IDiscreteFunction& f, const IDiscreteFunction& g)
+{
+  return (f.dataType() == g.dataType()) and (f.descriptor().type() == g.descriptor().type());
+}
+
+PUGS_INLINE
+bool
+isSameDiscretization(const std::shared_ptr<const IDiscreteFunction>& f,
+                     const std::shared_ptr<const IDiscreteFunction>& g)
+{
+  return (f->dataType() == g->dataType()) and (f->descriptor().type() == g->descriptor().type());
 }
 
 template <typename LHS_T, typename RHS_T>
@@ -56,16 +72,36 @@ innerCompositionLaw(const std::shared_ptr<const IDiscreteFunction>& f,
                     const std::shared_ptr<const IDiscreteFunction>& g)
 {
   Assert(f->mesh() == g->mesh());
-  Assert(f->dataType() == g->dataType());
+  Assert(isSameDiscretization(f, g));
 
   switch (f->dataType()) {
   case ASTNodeDataType::double_t: {
-    auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*f);
-    auto gh = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*g);
+    if (f->descriptor().type() == DiscreteFunctionType::P0) {
+      auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*f);
+      auto gh = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*g);
 
-    return innerCompositionLaw<BinOperatorT>(fh, gh);
+      return innerCompositionLaw<BinOperatorT>(fh, gh);
+
+    } else if (f->descriptor().type() == DiscreteFunctionType::P0Vector) {
+      if constexpr (std::is_same_v<BinOperatorT, language::plus_op> or
+                    std::is_same_v<BinOperatorT, language::minus_op>) {
+        auto fh = dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*f);
+        auto gh = dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*g);
+
+        if (fh.size() != gh.size()) {
+          throw NormalError(name(f) + " spaces have different sizes");
+        }
+
+        return innerCompositionLaw<BinOperatorT>(fh, gh);
+      } else {
+        throw NormalError(invalid_operands(f, g));
+      }
+    } else {
+      throw UnexpectedError(invalid_operands(f, g));
+    }
   }
   case ASTNodeDataType::vector_t: {
+    Assert(f->descriptor().type() == DiscreteFunctionType::P0);
     switch (f->dataType().dimension()) {
     case 1: {
       auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<1>>&>(*f);
@@ -91,6 +127,7 @@ innerCompositionLaw(const std::shared_ptr<const IDiscreteFunction>& f,
     }
   }
   case ASTNodeDataType::matrix_t: {
+    Assert(f->descriptor().type() == DiscreteFunctionType::P0);
     Assert(f->dataType().nbRows() == f->dataType().nbColumns());
     switch (f->dataType().nbRows()) {
     case 1: {
@@ -130,7 +167,7 @@ innerCompositionLaw(const std::shared_ptr<const IDiscreteFunction>& f,
   if (f->mesh() != g->mesh()) {
     throw NormalError("discrete functions defined on different meshes");
   }
-  if (f->dataType() != g->dataType()) {
+  if (not isSameDiscretization(f, g)) {
     throw NormalError(invalid_operands(f, g));
   }
 
@@ -155,10 +192,8 @@ std::shared_ptr<const IDiscreteFunction>
 applyBinaryOperation(const LeftDiscreteFunctionT& lhs, const RightDiscreteFunctionT& rhs)
 {
   Assert(lhs.mesh() == rhs.mesh());
-  using lhs_data_type = typename LeftDiscreteFunctionT::data_type;
-  using rhs_data_type = typename RightDiscreteFunctionT::data_type;
 
-  static_assert(not std::is_same_v<rhs_data_type, lhs_data_type>,
+  static_assert(not std::is_same_v<LeftDiscreteFunctionT, RightDiscreteFunctionT>,
                 "use innerCompositionLaw when data types are the same");
 
   return std::make_shared<decltype(BinOp<BinOperatorT>{}.eval(lhs, rhs))>(BinOp<BinOperatorT>{}.eval(lhs, rhs));
@@ -169,7 +204,7 @@ std::shared_ptr<const IDiscreteFunction>
 applyBinaryOperation(const DiscreteFunctionT& fh, const std::shared_ptr<const IDiscreteFunction>& g)
 {
   Assert(fh.mesh() == g->mesh());
-  Assert(fh.dataType() != g->dataType());
+  Assert(not isSameDiscretization(fh, *g));
   using lhs_data_type = std::decay_t<typename DiscreteFunctionT::data_type>;
 
   switch (g->dataType()) {
@@ -178,6 +213,15 @@ applyBinaryOperation(const DiscreteFunctionT& fh, const std::shared_ptr<const ID
       if constexpr (not is_tiny_matrix_v<lhs_data_type>) {
         auto gh = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*g);
 
+        return applyBinaryOperation<BinOperatorT>(fh, gh);
+      } else {
+        throw NormalError(invalid_operands(fh, g));
+      }
+    } else if constexpr (std::is_same_v<BinOperatorT, language::multiply_op> and
+                         std::is_same_v<DiscreteFunctionT, DiscreteFunctionP0<Dimension, double>>) {
+      if (g->descriptor().type() == DiscreteFunctionType::P0Vector) {
+        auto gh = dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*g);
+
         return applyBinaryOperation<BinOperatorT>(fh, gh);
       } else {
         throw NormalError(invalid_operands(fh, g));
@@ -266,7 +310,7 @@ applyBinaryOperation(const std::shared_ptr<const IDiscreteFunction>& f,
                      const std::shared_ptr<const IDiscreteFunction>& g)
 {
   Assert(f->mesh() == g->mesh());
-  Assert(f->dataType() != g->dataType());
+  Assert(not isSameDiscretization(f, g));
 
   switch (f->dataType()) {
   case ASTNodeDataType::double_t: {
@@ -312,7 +356,7 @@ applyBinaryOperation(const std::shared_ptr<const IDiscreteFunction>& f,
     throw NormalError("functions defined on different meshes");
   }
 
-  Assert(f->dataType() != g->dataType(), "should call inner composition instead");
+  Assert(not isSameDiscretization(f, g), "should call inner composition instead");
 
   switch (f->mesh()->dimension()) {
   case 1: {
@@ -345,7 +389,7 @@ operator-(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_p
 std::shared_ptr<const IDiscreteFunction>
 operator*(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<const IDiscreteFunction>& g)
 {
-  if (f->dataType() == g->dataType()) {
+  if (isSameDiscretization(f, g)) {
     return innerCompositionLaw<language::multiply_op>(f, g);
   } else {
     return applyBinaryOperation<language::multiply_op>(f, g);
@@ -355,7 +399,7 @@ operator*(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_p
 std::shared_ptr<const IDiscreteFunction>
 operator/(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<const IDiscreteFunction>& g)
 {
-  if (f->dataType() == g->dataType()) {
+  if (isSameDiscretization(f, g)) {
     return innerCompositionLaw<language::divide_op>(f, g);
   } else {
     return applyBinaryOperation<language::divide_op>(f, g);
@@ -392,8 +436,15 @@ applyBinaryOperationWithLeftConstant(const DataType& a, const std::shared_ptr<co
   case ASTNodeDataType::unsigned_int_t:
   case ASTNodeDataType::int_t:
   case ASTNodeDataType::double_t: {
-    auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*f);
-    return applyBinaryOperationWithLeftConstant<BinOperatorT>(a, fh);
+    if (f->descriptor().type() == DiscreteFunctionType::P0) {
+      auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*f);
+      return applyBinaryOperationWithLeftConstant<BinOperatorT>(a, fh);
+    } else if (f->descriptor().type() == DiscreteFunctionType::P0Vector) {
+      auto fh = dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*f);
+      return applyBinaryOperationWithLeftConstant<BinOperatorT>(a, fh);
+    } else {
+      throw NormalError(invalid_operands(a, f));
+    }
   }
   case ASTNodeDataType::vector_t: {
     if constexpr (is_tiny_matrix_v<DataType>) {
@@ -528,6 +579,8 @@ template <typename BinOperatorT, typename DataType, typename DiscreteFunctionT>
 std::shared_ptr<const IDiscreteFunction>
 applyBinaryOperationWithRightConstant(const DiscreteFunctionT& f, const DataType& a)
 {
+  Assert(f.descriptor().type() != DiscreteFunctionType::P0);
+
   using lhs_data_type = std::decay_t<typename DiscreteFunctionT::data_type>;
   using rhs_data_type = std::decay_t<DataType>;
 
@@ -549,6 +602,10 @@ template <typename BinOperatorT, size_t Dimension, typename DataType>
 std::shared_ptr<const IDiscreteFunction>
 applyBinaryOperationWithRightConstant(const std::shared_ptr<const IDiscreteFunction>& f, const DataType& a)
 {
+  if (f->descriptor().type() != DiscreteFunctionType::P0) {
+    throw NormalError(invalid_operands(f, a));
+  }
+
   switch (f->dataType()) {
   case ASTNodeDataType::bool_t:
   case ASTNodeDataType::unsigned_int_t:
diff --git a/src/language/utils/InterpolateItemArray.hpp b/src/language/utils/InterpolateItemArray.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..fcef1db90a2e17b07e54fe55966a96973d443b9e
--- /dev/null
+++ b/src/language/utils/InterpolateItemArray.hpp
@@ -0,0 +1,58 @@
+#ifndef INTERPOLATE_ITEM_ARRAY_HPP
+#define INTERPOLATE_ITEM_ARRAY_HPP
+
+#include <language/utils/InterpolateItemValue.hpp>
+#include <language/utils/PugsFunctionAdapter.hpp>
+#include <mesh/ItemArray.hpp>
+#include <mesh/ItemType.hpp>
+
+template <typename T>
+class InterpolateItemArray;
+template <typename OutputType, typename InputType>
+class InterpolateItemArray<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)>
+{
+  static constexpr size_t Dimension = OutputType::Dimension;
+  using Adapter                     = PugsFunctionAdapter<OutputType(InputType)>;
+
+ public:
+  template <ItemType item_type>
+  static inline 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()};
+
+    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;
+  }
+
+  template <ItemType item_type>
+  static inline Table<OutputType>
+  interpolate(const std::vector<FunctionSymbolId>& function_symbol_id_list,
+              const ItemArray<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()};
+
+    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;
+  }
+};
+
+#endif   // INTERPOLATE_ITEM_ARRAY_HPP
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index 5ca484b2258adfbe64b11cb6913c126442846447..20caa238cb25019df43eccd6c05c3a65192a0a94 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -658,11 +658,12 @@ class Connectivity final : public IConnectivity
   ~Connectivity() = default;
 };
 
-template <size_t Dim>
-PUGS_INLINE std::shared_ptr<Connectivity<Dim>>
-Connectivity<Dim>::build(const ConnectivityDescriptor& descriptor)
+template <size_t Dimension>
+PUGS_INLINE std::shared_ptr<Connectivity<Dimension>>
+Connectivity<Dimension>::build(const ConnectivityDescriptor& descriptor)
 {
-  std::shared_ptr<Connectivity<Dim>> connectivity_ptr(new Connectivity<Dim>);
+  // Cannot use std::make_shared in this context since the default constructor is private
+  std::shared_ptr<Connectivity<Dimension>> connectivity_ptr(new Connectivity<Dimension>);
   connectivity_ptr->_buildFrom(descriptor);
 
   return connectivity_ptr;
diff --git a/src/output/GnuplotWriter.cpp b/src/output/GnuplotWriter.cpp
index 78b1e25ce8277888819395cafe92118a46d992fb..5cb59553a54df6f1b7c173066079110058c84330 100644
--- a/src/output/GnuplotWriter.cpp
+++ b/src/output/GnuplotWriter.cpp
@@ -41,7 +41,7 @@ GnuplotWriter::_getFilename() const
 
 template <size_t Dimension>
 void
-GnuplotWriter::_writePreamble(const OutputNamedItemValueSet& output_named_item_value_set, std::ostream& fout) const
+GnuplotWriter::_writePreamble(const OutputNamedItemDataSet& output_named_item_data_set, std::ostream& fout) const
 {
   fout << "# list of data\n";
   fout << "# 1:x";
@@ -49,37 +49,49 @@ GnuplotWriter::_writePreamble(const OutputNamedItemValueSet& output_named_item_v
     fout << " 2:y";
   }
   uint64_t i = Dimension + 1;
-  for (const auto& i_named_item_value : output_named_item_value_set) {
-    const std::string name         = i_named_item_value.first;
-    const auto& item_value_variant = i_named_item_value.second;
+  for (const auto& i_named_item_data : output_named_item_data_set) {
+    const std::string name        = i_named_item_data.first;
+    const auto& item_data_variant = i_named_item_data.second;
     std::visit(
-      [&](auto&& item_value) {
-        using ItemValueType = std::decay_t<decltype(item_value)>;
-        using DataType      = std::decay_t<typename ItemValueType::data_type>;
-        if constexpr (std::is_arithmetic_v<DataType>) {
-          fout << ' ' << i++ << ':' << name;
-        } else if constexpr (is_tiny_vector_v<DataType>) {
-          for (size_t j = 0; j < DataType{}.dimension(); ++j) {
-            fout << ' ' << i++ << ':' << name << '[' << j << ']';
+      [&](auto&& item_data) {
+        using ItemDataType = std::decay_t<decltype(item_data)>;
+        using DataType     = std::decay_t<typename ItemDataType::data_type>;
+        if constexpr (is_item_value_v<ItemDataType>) {
+          if constexpr (std::is_arithmetic_v<DataType>) {
+            fout << ' ' << i++ << ':' << name;
+          } else if constexpr (is_tiny_vector_v<DataType>) {
+            for (size_t j = 0; j < DataType{}.dimension(); ++j) {
+              fout << ' ' << i++ << ':' << name << '[' << j << ']';
+            }
+          } else if constexpr (is_tiny_matrix_v<DataType>) {
+            for (size_t j = 0; j < DataType{}.nbRows(); ++j) {
+              for (size_t k = 0; k < DataType{}.nbColumns(); ++k) {
+                fout << ' ' << i++ << ':' << name << '(' << j << ',' << k << ')';
+              }
+            }
+          } else {
+            throw UnexpectedError("invalid data type");
           }
-        } else if constexpr (is_tiny_matrix_v<DataType>) {
-          for (size_t j = 0; j < DataType{}.nbRows(); ++j) {
-            for (size_t k = 0; k < DataType{}.nbColumns(); ++k) {
-              fout << ' ' << i++ << ':' << name << '(' << j << ',' << k << ')';
+        } else if constexpr (is_item_array_v<ItemDataType>) {
+          if constexpr (std::is_arithmetic_v<DataType>) {
+            for (size_t j = 0; j < item_data.sizeOfArrays(); ++j) {
+              fout << ' ' << i++ << ':' << name << '[' << j << ']';
             }
+          } else {
+            throw UnexpectedError("invalid data type");
           }
         } else {
-          throw UnexpectedError("invalid data type");
+          throw UnexpectedError("invalid ItemData type");
         }
       },
-      item_value_variant);
+      item_data_variant);
   }
   fout << "\n\n";
 }
 
 template <typename DataType>
 void
-GnuplotWriter::_writeCellValue(const CellValue<DataType>& cell_value, CellId cell_id, std::ostream& fout) const
+GnuplotWriter::_writeCellData(const CellValue<DataType>& cell_value, CellId cell_id, std::ostream& fout) const
 {
   const auto& value = cell_value[cell_id];
   if constexpr (std::is_arithmetic_v<DataType>) {
@@ -99,17 +111,31 @@ GnuplotWriter::_writeCellValue(const CellValue<DataType>& cell_value, CellId cel
   }
 }
 
-template <typename DataType, ItemType item_type>
+template <typename DataType>
+void
+GnuplotWriter::_writeCellData(const CellArray<DataType>& cell_array, CellId cell_id, std::ostream& fout) const
+{
+  const auto& array = cell_array[cell_id];
+  if constexpr (std::is_arithmetic_v<DataType>) {
+    for (size_t i = 0; i < array.size(); ++i) {
+      fout << ' ' << array[i];
+    }
+  } else {
+    throw UnexpectedError("invalid data type for cell value output: " + demangle<DataType>());
+  }
+}
+
+template <typename ItemDataT>
 void
-GnuplotWriter::_writeValue(const ItemValue<DataType, item_type>& item_value,
-                           [[maybe_unused]] CellId cell_id,
-                           [[maybe_unused]] NodeId node_id,
-                           std::ostream& fout) const
+GnuplotWriter::_writeData(const ItemDataT& item_data,
+                          [[maybe_unused]] CellId cell_id,
+                          [[maybe_unused]] NodeId node_id,
+                          std::ostream& fout) const
 {
-  if constexpr (item_type == ItemType::cell) {
-    this->_writeCellValue(item_value, cell_id, fout);
-  } else if constexpr (item_type == ItemType::node) {
-    this->_writeNodeValue(item_value, node_id, fout);
+  if constexpr (ItemDataT::item_t == ItemType::cell) {
+    this->_writeCellData(item_data, cell_id, fout);
+  } else if constexpr (ItemDataT::item_t == ItemType::node) {
+    this->_writeNodeData(item_data, node_id, fout);
   } else {
     throw UnexpectedError{"invalid item type"};
   }
@@ -117,9 +143,9 @@ GnuplotWriter::_writeValue(const ItemValue<DataType, item_type>& item_value,
 
 template <typename MeshType>
 void
-GnuplotWriter::_writeValues(const std::shared_ptr<const MeshType>& mesh,
-                            const OutputNamedItemValueSet& output_named_item_value_set,
-                            std::ostream& fout) const
+GnuplotWriter::_writeDataAtNodes(const std::shared_ptr<const MeshType>& mesh,
+                                 const OutputNamedItemDataSet& output_named_item_data_set,
+                                 std::ostream& fout) const
 {
   if constexpr (MeshType::Dimension == 1) {
     auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
@@ -132,8 +158,8 @@ GnuplotWriter::_writeValues(const std::shared_ptr<const MeshType>& mesh,
           const NodeId& node_id                     = cell_nodes[i_node];
           const TinyVector<MeshType::Dimension>& xr = mesh->xr()[node_id];
           fout << xr[0];
-          for (auto [name, item_value] : output_named_item_value_set) {
-            std::visit([&](auto&& cell_value) { _writeValue(cell_value, cell_id, node_id, fout); }, item_value);
+          for (auto [name, item_data_variant] : output_named_item_data_set) {
+            std::visit([&](auto&& item_data) { _writeData(item_data, cell_id, node_id, fout); }, item_data_variant);
           }
           fout << '\n';
         }
@@ -152,16 +178,16 @@ GnuplotWriter::_writeValues(const std::shared_ptr<const MeshType>& mesh,
           const NodeId& node_id                     = cell_nodes[i_node];
           const TinyVector<MeshType::Dimension>& xr = mesh->xr()[node_id];
           fout << xr[0] << ' ' << xr[1];
-          for (auto [name, item_value] : output_named_item_value_set) {
-            std::visit([&](auto&& cell_value) { _writeValue(cell_value, cell_id, node_id, fout); }, item_value);
+          for (auto [name, item_data_variant] : output_named_item_data_set) {
+            std::visit([&](auto&& item_data) { _writeData(item_data, cell_id, node_id, fout); }, item_data_variant);
           }
           fout << '\n';
         }
         const NodeId& node_id                     = cell_nodes[0];
         const TinyVector<MeshType::Dimension>& xr = mesh->xr()[node_id];
         fout << xr[0] << ' ' << xr[1];
-        for (auto [name, item_value] : output_named_item_value_set) {
-          std::visit([&](auto&& cell_value) { _writeValue(cell_value, cell_id, node_id, fout); }, item_value);
+        for (auto [name, item_data_variant] : output_named_item_data_set) {
+          std::visit([&](auto&& item_data) { _writeData(item_data, cell_id, node_id, fout); }, item_data_variant);
         }
         fout << "\n\n\n";
       }
@@ -171,32 +197,46 @@ GnuplotWriter::_writeValues(const std::shared_ptr<const MeshType>& mesh,
   }
 }
 
-template <typename ValueType>
+template <typename DataType>
 void
-GnuplotWriter::_writeNodeValue(const NodeValue<ValueType>& node_value, NodeId node_id, std::ostream& fout) const
+GnuplotWriter::_writeNodeData(const NodeValue<DataType>& node_value, NodeId node_id, std::ostream& fout) const
 {
   const auto& value = node_value[node_id];
-  if constexpr (std::is_arithmetic_v<ValueType>) {
+  if constexpr (std::is_arithmetic_v<DataType>) {
     fout << ' ' << value;
-  } else if constexpr (is_tiny_vector_v<std::decay_t<ValueType>>) {
+  } else if constexpr (is_tiny_vector_v<std::decay_t<DataType>>) {
     for (size_t i = 0; i < value.dimension(); ++i) {
       fout << ' ' << value[i];
     }
-  } else if constexpr (is_tiny_matrix_v<std::decay_t<ValueType>>) {
+  } else if constexpr (is_tiny_matrix_v<std::decay_t<DataType>>) {
     for (size_t i = 0; i < value.nbRows(); ++i) {
       for (size_t j = 0; j < value.nbColumns(); ++j) {
         fout << ' ' << value(i, j);
       }
     }
   } else {
-    throw UnexpectedError("invalid data type for cell value output: " + demangle<ValueType>());
+    throw UnexpectedError("invalid data type for cell value output: " + demangle<DataType>());
+  }
+}
+
+template <typename DataType>
+void
+GnuplotWriter::_writeNodeData(const NodeArray<DataType>& node_array, NodeId node_id, std::ostream& fout) const
+{
+  const auto& array = node_array[node_id];
+  if constexpr (std::is_arithmetic_v<DataType>) {
+    for (size_t i = 0; i < array.size(); ++i) {
+      fout << ' ' << array[i];
+    }
+  } else {
+    throw UnexpectedError("invalid data type for cell value output: " + demangle<DataType>());
   }
 }
 
 template <typename MeshType>
 void
 GnuplotWriter::_write(const std::shared_ptr<const MeshType>& mesh,
-                      const OutputNamedItemValueSet& output_named_item_value_set,
+                      const OutputNamedItemDataSet& output_named_item_data_set,
                       double time) const
 {
   if (parallel::rank() == 0) {
@@ -207,7 +247,7 @@ GnuplotWriter::_write(const std::shared_ptr<const MeshType>& mesh,
     fout << this->_getDateAndVersionComment();
     fout << "# time = " << time << "\n\n";
 
-    this->_writePreamble<MeshType::Dimension>(output_named_item_value_set, fout);
+    this->_writePreamble<MeshType::Dimension>(output_named_item_data_set, fout);
   }
 
   for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
@@ -216,7 +256,7 @@ GnuplotWriter::_write(const std::shared_ptr<const MeshType>& mesh,
       fout.precision(15);
       fout.setf(std::ios_base::scientific);
 
-      this->_writeValues(mesh, output_named_item_value_set, fout);
+      this->_writeDataAtNodes(mesh, output_named_item_data_set, fout);
     }
     parallel::barrier();
   }
@@ -225,15 +265,15 @@ GnuplotWriter::_write(const std::shared_ptr<const MeshType>& mesh,
 void
 GnuplotWriter::writeMesh(const std::shared_ptr<const IMesh>& p_mesh) const
 {
-  OutputNamedItemValueSet output_named_item_value_set{};
+  OutputNamedItemDataSet output_named_item_data_set{};
 
   switch (p_mesh->dimension()) {
   case 1: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(p_mesh), output_named_item_value_set, 0);
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(p_mesh), output_named_item_data_set, 0);
     break;
   }
   case 2: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(p_mesh), output_named_item_value_set, 0);
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(p_mesh), output_named_item_data_set, 0);
     break;
   }
   default: {
@@ -248,15 +288,15 @@ GnuplotWriter::write(const std::vector<std::shared_ptr<const NamedDiscreteFuncti
 {
   std::shared_ptr mesh = this->_getMesh(named_discrete_function_list);
 
-  OutputNamedItemValueSet output_named_item_value_set = this->_getOutputNamedItemValueSet(named_discrete_function_list);
+  OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_function_list);
 
   switch (mesh->dimension()) {
   case 1: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_value_set, time);
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_data_set, time);
     break;
   }
   case 2: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_value_set, time);
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_data_set, time);
     break;
   }
   default: {
diff --git a/src/output/GnuplotWriter.hpp b/src/output/GnuplotWriter.hpp
index bf6487c998c29d8a5c12177a574eea587a1fb6ea..608d830152cb9c90d6f8b4014bbf2d8054099fe7 100644
--- a/src/output/GnuplotWriter.hpp
+++ b/src/output/GnuplotWriter.hpp
@@ -19,28 +19,34 @@ class GnuplotWriter : public WriterBase
   std::string _getFilename() const;
 
   template <typename DataType>
-  void _writeCellValue(const CellValue<DataType>& cell_value, CellId cell_id, std::ostream& fout) const;
+  void _writeCellData(const CellValue<DataType>& cell_value, CellId cell_id, std::ostream& fout) const;
 
   template <typename DataType>
-  void _writeNodeValue(const NodeValue<DataType>& node_value, NodeId cell_id, std::ostream& fout) const;
+  void _writeCellData(const CellArray<DataType>& cell_value, CellId cell_id, std::ostream& fout) const;
+
+  template <typename DataType>
+  void _writeNodeData(const NodeValue<DataType>& node_value, NodeId cell_id, std::ostream& fout) const;
+
+  template <typename DataType>
+  void _writeNodeData(const NodeArray<DataType>& node_value, NodeId cell_id, std::ostream& fout) const;
 
   template <size_t Dimension>
-  void _writePreamble(const OutputNamedItemValueSet& output_named_item_value_set, std::ostream& fout) const;
+  void _writePreamble(const OutputNamedItemDataSet& output_named_item_data_set, std::ostream& fout) const;
 
-  template <typename DataType, ItemType item_type>
-  void _writeValue(const ItemValue<DataType, item_type>& item_value,
-                   CellId cell_id,
-                   NodeId node_id,
-                   std::ostream& fout) const;
+  template <typename ItemDataT>
+  void _writeData(const ItemDataT& item_data,
+                  [[maybe_unused]] CellId cell_id,
+                  [[maybe_unused]] NodeId node_id,
+                  std::ostream& fout) const;
 
   template <typename MeshType>
-  void _writeValues(const std::shared_ptr<const MeshType>& mesh,
-                    const OutputNamedItemValueSet& output_named_item_value_set,
-                    std::ostream& fout) const;
+  void _writeDataAtNodes(const std::shared_ptr<const MeshType>& mesh,
+                         const OutputNamedItemDataSet& output_named_item_data_set,
+                         std::ostream& fout) const;
 
   template <typename MeshType>
   void _write(const std::shared_ptr<const MeshType>& mesh,
-              const OutputNamedItemValueSet& output_named_item_value_set,
+              const OutputNamedItemDataSet& output_named_item_data_set,
               double time) const;
 
  public:
diff --git a/src/output/GnuplotWriter1D.cpp b/src/output/GnuplotWriter1D.cpp
index 79423d5da05b699b995d472313e134a8c9dd4a52..31320054e6f244e9e0771190c539eecd733d66ee 100644
--- a/src/output/GnuplotWriter1D.cpp
+++ b/src/output/GnuplotWriter1D.cpp
@@ -39,71 +39,69 @@ GnuplotWriter1D::_getFilename() const
   return sout.str();
 }
 
-template <typename DataType>
+template <typename ItemDataType>
 bool
-GnuplotWriter1D::_is_cell_value(const CellValue<const DataType>&) const
+GnuplotWriter1D::_is_cell_data(const ItemDataType&) const
 {
-  return true;
+  return ItemDataType::item_t == ItemType::cell;
 }
 
-template <typename DataType>
+template <typename ItemDataType>
 bool
-GnuplotWriter1D::_is_cell_value(const NodeValue<const DataType>&) const
+GnuplotWriter1D::_is_node_data(const ItemDataType&) const
 {
-  return false;
-}
-
-template <typename DataType>
-bool
-GnuplotWriter1D::_is_node_value(const CellValue<const DataType>&) const
-{
-  return false;
-}
-
-template <typename DataType>
-bool
-GnuplotWriter1D::_is_node_value(const NodeValue<const DataType>&) const
-{
-  return true;
+  return ItemDataType::item_t == ItemType::node;
 }
 
 void
-GnuplotWriter1D::_writePreamble(const OutputNamedItemValueSet& output_named_item_value_set, std::ostream& fout) const
+GnuplotWriter1D::_writePreamble(const OutputNamedItemDataSet& output_named_item_data_set, std::ostream& fout) const
 {
   fout << "# list of data\n";
   fout << "# 1:x";
   uint64_t i = 2;
-  for (const auto& i_named_item_value : output_named_item_value_set) {
-    const std::string name         = i_named_item_value.first;
-    const auto& item_value_variant = i_named_item_value.second;
+  for (const auto& i_named_item_data : output_named_item_data_set) {
+    const std::string name        = i_named_item_data.first;
+    const auto& item_data_variant = i_named_item_data.second;
     std::visit(
-      [&](auto&& item_value) {
-        using ItemValueType = std::decay_t<decltype(item_value)>;
-        using DataType      = std::decay_t<typename ItemValueType::data_type>;
-        if constexpr (std::is_arithmetic_v<DataType>) {
-          fout << ' ' << i++ << ':' << name;
-        } else if constexpr (is_tiny_vector_v<DataType>) {
-          for (size_t j = 0; j < DataType{}.dimension(); ++j) {
-            fout << ' ' << i++ << ':' << name << '[' << j << ']';
+      [&](auto&& item_data) {
+        using ItemDataType = std::decay_t<decltype(item_data)>;
+        using DataType     = std::decay_t<typename ItemDataType::data_type>;
+        if constexpr (is_item_value_v<ItemDataType>) {
+          if constexpr (std::is_arithmetic_v<DataType>) {
+            fout << ' ' << i++ << ':' << name;
+          } else if constexpr (is_tiny_vector_v<DataType>) {
+            for (size_t j = 0; j < DataType{}.dimension(); ++j) {
+              fout << ' ' << i++ << ':' << name << '[' << j << ']';
+            }
+          } else if constexpr (is_tiny_matrix_v<DataType>) {
+            for (size_t j = 0; j < DataType{}.nbRows(); ++j) {
+              for (size_t k = 0; k < DataType{}.nbColumns(); ++k) {
+                fout << ' ' << i++ << ':' << name << '(' << j << ',' << k << ')';
+              }
+            }
+          } else {
+            throw UnexpectedError("invalid data type");
           }
-        } else if constexpr (is_tiny_matrix_v<DataType>) {
-          for (size_t j = 0; j < DataType{}.nbRows(); ++j) {
-            for (size_t k = 0; k < DataType{}.nbColumns(); ++k) {
-              fout << ' ' << i++ << ':' << name << '(' << j << ',' << k << ')';
+        } else if constexpr (is_item_array_v<ItemDataType>) {
+          if constexpr (std::is_arithmetic_v<DataType>) {
+            for (size_t j = 0; j < item_data.sizeOfArrays(); ++j) {
+              fout << ' ' << i++ << ':' << name << '[' << j << ']';
             }
+          } else {
+            throw UnexpectedError("invalid data type");
           }
         } else {
-          throw UnexpectedError("invalid data type");
+          throw UnexpectedError("invalid ItemData type");
         }
       },
-      item_value_variant);
+      item_data_variant);
   }
   fout << "\n\n";
 }
 
 template <typename DataType, ItemType item_type>
 size_t
-GnuplotWriter1D::_itemValueNbRow(const ItemValue<DataType, item_type>&) const
+GnuplotWriter1D::_itemDataNbRow(const ItemValue<DataType, item_type>&) const
 {
   if constexpr (std::is_arithmetic_v<DataType>) {
     return 1;
@@ -116,18 +114,25 @@ GnuplotWriter1D::_itemValueNbRow(const ItemValue<DataType, item_type>&) const
   }
 }
 
+template <typename DataType, ItemType item_type>
+size_t
+GnuplotWriter1D::_itemDataNbRow(const ItemArray<DataType, item_type>& item_array) const
+{
+  return item_array.sizeOfArrays();
+}
+
 template <typename MeshType, ItemType item_type>
 void
-GnuplotWriter1D::_writeItemValues(const std::shared_ptr<const MeshType>& mesh,
-                                  const OutputNamedItemValueSet& output_named_item_value_set,
-                                  std::ostream& fout) const
+GnuplotWriter1D::_writeItemDatas(const std::shared_ptr<const MeshType>& mesh,
+                                 const OutputNamedItemDataSet& output_named_item_data_set,
+                                 std::ostream& fout) const
 {
   using ItemId = ItemIdT<item_type>;
 
   const size_t& number_of_columns = [&] {
     size_t number_of_columns = 1;
-    for (auto [name, item_value] : output_named_item_value_set) {
-      std::visit([&](auto&& value) { number_of_columns += _itemValueNbRow(value); }, item_value);
+    for (auto [name, item_data] : output_named_item_data_set) {
+      std::visit([&](auto&& value) { number_of_columns += _itemDataNbRow(value); }, item_data);
     }
     return number_of_columns;
   }();
@@ -175,37 +180,53 @@ GnuplotWriter1D::_writeItemValues(const std::shared_ptr<const MeshType>& mesh,
   }
 
   size_t column_number = 1;
-  for (auto [name, output_item_value] : output_named_item_value_set) {
+  for (auto [name, output_item_data] : output_named_item_data_set) {
     std::visit(
-      [&](auto&& item_value) {
-        using ItemValueT = std::decay_t<decltype(item_value)>;
-        if constexpr (ItemValueT::item_t == item_type) {
-          using DataT  = std::decay_t<typename ItemValueT::data_type>;
-          size_t index = 0;
-          for (ItemId item_id = 0; item_id < item_value.numberOfItems(); ++item_id) {
-            if (is_owned[item_id]) {
-              if constexpr (std::is_arithmetic_v<DataT>) {
-                values[number_of_columns * index + column_number] = item_value[item_id];
-              } else if constexpr (is_tiny_vector_v<DataT>) {
-                const size_t k = number_of_columns * index + column_number;
-                for (size_t j = 0; j < DataT::Dimension; ++j) {
-                  values[k + j] = item_value[item_id][j];
+      [&](auto&& item_data) {
+        using ItemDataT = std::decay_t<decltype(item_data)>;
+        if constexpr (ItemDataT::item_t == item_type) {
+          if constexpr (is_item_value_v<ItemDataT>) {
+            using DataT  = std::decay_t<typename ItemDataT::data_type>;
+            size_t index = 0;
+            for (ItemId item_id = 0; item_id < item_data.numberOfItems(); ++item_id) {
+              if (is_owned[item_id]) {
+                if constexpr (std::is_arithmetic_v<DataT>) {
+                  values[number_of_columns * index + column_number] = item_data[item_id];
+                } else if constexpr (is_tiny_vector_v<DataT>) {
+                  const size_t k = number_of_columns * index + column_number;
+                  for (size_t j = 0; j < DataT::Dimension; ++j) {
+                    values[k + j] = item_data[item_id][j];
+                  }
+                } else if constexpr (is_tiny_matrix_v<DataT>) {
+                  size_t k = number_of_columns * index + column_number;
+                  for (size_t i = 0; i < DataT{}.nbRows(); ++i) {
+                    for (size_t j = 0; j < DataT{}.nbColumns(); ++j) {
+                      values[k++] = item_data[item_id](i, j);
+                    }
+                  }
                 }
-              } else if constexpr (is_tiny_matrix_v<DataT>) {
-                size_t k = number_of_columns * index + column_number;
-                for (size_t i = 0; i < DataT{}.nbRows(); ++i) {
-                  for (size_t j = 0; j < DataT{}.nbColumns(); ++j) {
-                    values[k++] = item_value[item_id](i, j);
+                ++index;
+              }
+            }
+          } else {
+            using DataT  = std::decay_t<typename ItemDataT::data_type>;
+            size_t index = 0;
+            for (ItemId item_id = 0; item_id < item_data.numberOfItems(); ++item_id) {
+              if (is_owned[item_id]) {
+                if constexpr (std::is_arithmetic_v<DataT>) {
+                  const size_t k = number_of_columns * index + column_number;
+                  for (size_t j = 0; j < item_data.sizeOfArrays(); ++j) {
+                    values[k + j] = item_data[item_id][j];
                   }
                 }
+                ++index;
               }
-              ++index;
             }
           }
         }
-        column_number += _itemValueNbRow(item_value);
+        column_number += _itemDataNbRow(item_data);
       },
-      output_item_value);
+      output_item_data);
   }
 
   if (parallel::size() > 1) {
@@ -236,23 +257,21 @@ GnuplotWriter1D::_writeItemValues(const std::shared_ptr<const MeshType>& mesh,
 template <typename MeshType>
 void
 GnuplotWriter1D::_write(const std::shared_ptr<const MeshType>& mesh,
-                        const OutputNamedItemValueSet& output_named_item_value_set,
+                        const OutputNamedItemDataSet& output_named_item_data_set,
                         double time) const
 {
-  bool has_cell_value = false;
-  for (const auto& [name, item_value_variant] : output_named_item_value_set) {
-    has_cell_value |=
-      std::visit([&](auto&& item_value) { return this->_is_cell_value(item_value); }, item_value_variant);
+  bool has_cell_data = false;
+  for (const auto& [name, item_data_variant] : output_named_item_data_set) {
+    has_cell_data |= std::visit([&](auto&& item_data) { return this->_is_cell_data(item_data); }, item_data_variant);
   }
 
-  bool has_node_value = false;
-  for (const auto& [name, item_value_variant] : output_named_item_value_set) {
-    has_node_value |=
-      std::visit([&](auto&& item_value) { return this->_is_node_value(item_value); }, item_value_variant);
+  bool has_node_data = false;
+  for (const auto& [name, item_data_variant] : output_named_item_data_set) {
+    has_node_data |= std::visit([&](auto&& item_data) { return this->_is_node_data(item_data); }, item_data_variant);
   }
 
-  if (has_cell_value and has_node_value) {
-    throw NormalError("cannot store both node and cell values in a gnuplot file");
+  if (has_cell_data and has_node_data) {
+    throw NormalError("cannot store both node and cell data in the same gnuplot file");
   }
 
   std::ofstream fout;
@@ -265,13 +284,13 @@ GnuplotWriter1D::_write(const std::shared_ptr<const MeshType>& mesh,
 
     fout << "# time = " << time << "\n\n";
 
-    _writePreamble(output_named_item_value_set, fout);
+    _writePreamble(output_named_item_data_set, fout);
   }
 
-  if (has_cell_value) {
-    this->_writeItemValues<MeshType, ItemType::cell>(mesh, output_named_item_value_set, fout);
+  if (has_cell_data) {
+    this->_writeItemDatas<MeshType, ItemType::cell>(mesh, output_named_item_data_set, fout);
   } else {   // has_node_value
-    this->_writeItemValues<MeshType, ItemType::node>(mesh, output_named_item_value_set, fout);
+    this->_writeItemDatas<MeshType, ItemType::node>(mesh, output_named_item_data_set, fout);
   }
 }
 
@@ -287,11 +306,11 @@ GnuplotWriter1D::write(const std::vector<std::shared_ptr<const NamedDiscreteFunc
 {
   std::shared_ptr mesh = this->_getMesh(named_discrete_function_list);
 
-  OutputNamedItemValueSet output_named_item_value_set = this->_getOutputNamedItemValueSet(named_discrete_function_list);
+  OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_function_list);
 
   switch (mesh->dimension()) {
   case 1: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_value_set, time);
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_data_set, time);
     break;
   }
   case 2: {
diff --git a/src/output/GnuplotWriter1D.hpp b/src/output/GnuplotWriter1D.hpp
index 6402a6002a105e1ef8b85b60783dce7a1cba82bd..b003d1756f8917dadf7350adada701052494b581 100644
--- a/src/output/GnuplotWriter1D.hpp
+++ b/src/output/GnuplotWriter1D.hpp
@@ -18,31 +18,28 @@ class GnuplotWriter1D : public WriterBase
 
   std::string _getFilename() const;
 
-  template <typename DataType>
-  bool _is_cell_value(const CellValue<const DataType>&) const;
+  template <typename ItemDataType>
+  bool _is_cell_data(const ItemDataType&) const;
 
-  template <typename DataType>
-  bool _is_cell_value(const NodeValue<const DataType>&) const;
+  template <typename ItemDataType>
+  bool _is_node_data(const ItemDataType&) const;
 
-  template <typename DataType>
-  bool _is_node_value(const CellValue<const DataType>&) const;
-
-  template <typename DataType>
-  bool _is_node_value(const NodeValue<const DataType>&) const;
+  template <typename DataType, ItemType item_type>
+  size_t _itemDataNbRow(const ItemValue<DataType, item_type>&) const;
 
   template <typename DataType, ItemType item_type>
-  size_t _itemValueNbRow(const ItemValue<DataType, item_type>&) const;
+  size_t _itemDataNbRow(const ItemArray<DataType, item_type>&) const;
 
   template <typename MeshType, ItemType item_type>
-  void _writeItemValues(const std::shared_ptr<const MeshType>& mesh,
-                        const OutputNamedItemValueSet& output_named_item_value_set,
-                        std::ostream& fout) const;
+  void _writeItemDatas(const std::shared_ptr<const MeshType>& mesh,
+                       const OutputNamedItemDataSet& output_named_item_data_set,
+                       std::ostream& fout) const;
 
-  void _writePreamble(const OutputNamedItemValueSet& output_named_item_value_set, std::ostream& fout) const;
+  void _writePreamble(const OutputNamedItemDataSet& output_named_item_value_set, std::ostream& fout) const;
 
   template <typename MeshType>
   void _write(const std::shared_ptr<const MeshType>& mesh,
-              const OutputNamedItemValueSet& output_named_item_value_set,
+              const OutputNamedItemDataSet& output_named_item_value_set,
               double time) const;
 
  public:
diff --git a/src/output/OutputNamedItemValueSet.hpp b/src/output/OutputNamedItemValueSet.hpp
index c7464135e3531c36ba3b1d3a801088cba002fbd7..bf90ddc3bad9636f2339695e87a446b11d81ae1a 100644
--- a/src/output/OutputNamedItemValueSet.hpp
+++ b/src/output/OutputNamedItemValueSet.hpp
@@ -1,6 +1,7 @@
-#ifndef OUTPUT_NAMED_ITEM_VALUE_SET_HPP
-#define OUTPUT_NAMED_ITEM_VALUE_SET_HPP
+#ifndef OUTPUT_NAMED_ITEM_DATA_SET_HPP
+#define OUTPUT_NAMED_ITEM_DATA_SET_HPP
 
+#include <mesh/ItemArray.hpp>
 #include <mesh/ItemValue.hpp>
 
 #include <algebra/TinyMatrix.hpp>
@@ -10,12 +11,17 @@
 #include <string>
 #include <variant>
 
-template <typename DataType, ItemType item_type>
-class NamedItemValue
+template <typename DataType,
+          ItemType item_type,
+          typename ConnectivityPtr,
+          template <typename DataTypeT, ItemType item_type_t, typename ConnectivityPtrT>
+          typename ItemDataT>
+class NamedItemData
 {
  private:
   std::string m_name;
-  ItemValue<const DataType, item_type> m_item_value;
+  using ItemDataType = ItemDataT<const DataType, item_type, std::shared_ptr<const IConnectivity>>;
+  ItemDataType m_item_data;
 
  public:
   constexpr const std::string&
@@ -24,77 +30,99 @@ class NamedItemValue
     return m_name;
   }
 
-  constexpr const ItemValue<const DataType, item_type>&
-  itemValue() const
+  constexpr const ItemDataType&
+  itemData() const
   {
-    return m_item_value;
+    return m_item_data;
   }
 
-  NamedItemValue& operator=(const NamedItemValue&) = default;
-  NamedItemValue& operator=(NamedItemValue&&) = default;
+  NamedItemData& operator=(const NamedItemData&) = default;
+  NamedItemData& operator=(NamedItemData&&) = default;
 
-  template <typename ConnectivityPtr>
-  NamedItemValue(const std::string& name, const ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
-    : m_name(name), m_item_value(item_value)
+  NamedItemData(const std::string& name, const ItemDataT<DataType, item_type, ConnectivityPtr>& item_data)
+    : m_name(name), m_item_data(item_data)
   {
     ;
   }
 
-  NamedItemValue(const std::string& name, const ItemValue<const DataType, item_type>& item_value)
-    : m_name(name), m_item_value(item_value)
+  NamedItemData(const std::string& name, const ItemDataT<const DataType, item_type, ConnectivityPtr>& item_data)
+    : m_name(name), m_item_data(item_data)
   {
     ;
   }
 
-  NamedItemValue(const NamedItemValue&) = default;
-  NamedItemValue(NamedItemValue&&)      = default;
-  ~NamedItemValue()                     = default;
+  NamedItemData(const NamedItemData&) = default;
+  NamedItemData(NamedItemData&&)      = default;
+  ~NamedItemData()                    = default;
 };
 
-class OutputNamedItemValueSet
+class OutputNamedItemDataSet
 {
  public:
-  using ItemValueVariant = std::variant<NodeValue<const bool>,
-                                        NodeValue<const int>,
-                                        NodeValue<const long int>,
-                                        NodeValue<const unsigned long int>,
-                                        NodeValue<const double>,
-                                        NodeValue<const TinyVector<1, double>>,
-                                        NodeValue<const TinyVector<2, double>>,
-                                        NodeValue<const TinyVector<3, double>>,
-                                        NodeValue<const TinyMatrix<1, double>>,
-                                        NodeValue<const TinyMatrix<2, double>>,
-                                        NodeValue<const TinyMatrix<3, double>>,
-
-                                        CellValue<const bool>,
-                                        CellValue<const int>,
-                                        CellValue<const long int>,
-                                        CellValue<const unsigned long int>,
-                                        CellValue<const double>,
-                                        CellValue<const TinyVector<1, double>>,
-                                        CellValue<const TinyVector<2, double>>,
-                                        CellValue<const TinyVector<3, double>>,
-                                        CellValue<const TinyMatrix<1, double>>,
-                                        CellValue<const TinyMatrix<2, double>>,
-                                        CellValue<const TinyMatrix<3, double>>>;
+  using ItemDataVariant = std::variant<NodeValue<const bool>,
+                                       NodeValue<const int>,
+                                       NodeValue<const long int>,
+                                       NodeValue<const unsigned long int>,
+                                       NodeValue<const double>,
+                                       NodeValue<const TinyVector<1, double>>,
+                                       NodeValue<const TinyVector<2, double>>,
+                                       NodeValue<const TinyVector<3, double>>,
+                                       NodeValue<const TinyMatrix<1, double>>,
+                                       NodeValue<const TinyMatrix<2, double>>,
+                                       NodeValue<const TinyMatrix<3, double>>,
+
+                                       CellValue<const bool>,
+                                       CellValue<const int>,
+                                       CellValue<const long int>,
+                                       CellValue<const unsigned long int>,
+                                       CellValue<const double>,
+                                       CellValue<const TinyVector<1, double>>,
+                                       CellValue<const TinyVector<2, double>>,
+                                       CellValue<const TinyVector<3, double>>,
+                                       CellValue<const TinyMatrix<1, double>>,
+                                       CellValue<const TinyMatrix<2, double>>,
+                                       CellValue<const TinyMatrix<3, double>>,
+
+                                       NodeArray<const bool>,
+                                       NodeArray<const int>,
+                                       NodeArray<const long int>,
+                                       NodeArray<const unsigned long int>,
+                                       NodeArray<const double>,
+
+                                       CellArray<const bool>,
+                                       CellArray<const int>,
+                                       CellArray<const long int>,
+                                       CellArray<const unsigned long int>,
+                                       CellArray<const double>>;
 
  private:
-  std::map<std::string, ItemValueVariant> m_name_itemvariant_map;
-
-  template <typename DataType, ItemType item_type>
+  std::map<std::string, ItemDataVariant> m_name_itemvariant_map;
+
+  template <typename DataType,
+            ItemType item_type,
+            typename ConnectivityPtr,
+            template <typename DataTypeT, ItemType item_type_t, typename ConnectivityPtrT>
+            typename ItemDataT,
+            typename... Args>
   PUGS_FORCEINLINE constexpr void
-  _doInsert(const NamedItemValue<DataType, item_type>& named_itemvalue)
+  _doInsert(const NamedItemData<DataType, item_type, ConnectivityPtr, ItemDataT>& named_item_data)
   {
-    if (m_name_itemvariant_map.find(named_itemvalue.name()) == m_name_itemvariant_map.end()) {
-      m_name_itemvariant_map[named_itemvalue.name()] = named_itemvalue.itemValue();
+    if (m_name_itemvariant_map.find(named_item_data.name()) == m_name_itemvariant_map.end()) {
+      m_name_itemvariant_map[named_item_data.name()] = named_item_data.itemData();
     }
   }
 
-  template <typename DataType, ItemType item_type, typename... Args>
+  template <typename DataType,
+            ItemType item_type,
+            typename ConnectivityPtr,
+            template <typename DataTypeT, ItemType item_type_t, typename ConnectivityPtrT>
+            typename ItemDataT,
+            typename... Args>
   PUGS_FORCEINLINE constexpr void
-  _unpackVariadicInput(const NamedItemValue<DataType, item_type>& named_itemvalue, Args&&... args)
+  _unpackVariadicInput(const NamedItemData<DataType, item_type, ConnectivityPtr, ItemDataT>& named_item_data,
+                       Args&&... args)
   {
-    _doInsert(named_itemvalue);
+    _doInsert(named_item_data);
     if constexpr (sizeof...(args) > 0) {
       this->_unpackVariadicInput(std::forward<Args>(args)...);
     }
@@ -113,22 +141,30 @@ class OutputNamedItemValueSet
     return m_name_itemvariant_map.end();
   }
 
-  template <typename DataType, ItemType item_type>
+  template <typename DataType,
+            ItemType item_type,
+            typename ConnectivityPtr,
+            template <typename DataTypeT, ItemType item_type_t, typename ConnectivityPtrT>
+            typename ItemDataT>
   void
-  add(const NamedItemValue<DataType, item_type>& named_itemvalue)
+  add(const NamedItemData<DataType, item_type, ConnectivityPtr, ItemDataT>& named_item_data)
   {
-    _doInsert(named_itemvalue);
+    _doInsert(named_item_data);
   }
 
-  template <typename... DataType, ItemType... item_type>
-  OutputNamedItemValueSet(NamedItemValue<DataType, item_type>... named_itemvalue)
+  template <typename... DataType,
+            ItemType... item_type,
+            typename... ConnectivityPtr,
+            template <typename DataTypeT, ItemType item_type_t, typename ConnectivityPtrT>
+            typename... ItemDataT>
+  OutputNamedItemDataSet(NamedItemData<DataType, item_type, ConnectivityPtr, ItemDataT>... named_item_data)
   {
-    _unpackVariadicInput(named_itemvalue...);
+    _unpackVariadicInput(named_item_data...);
   }
 
-  OutputNamedItemValueSet(const OutputNamedItemValueSet&) = default;
-  OutputNamedItemValueSet()                               = default;
-  ~OutputNamedItemValueSet()                              = default;
+  OutputNamedItemDataSet(const OutputNamedItemDataSet&) = default;
+  OutputNamedItemDataSet()                              = default;
+  ~OutputNamedItemDataSet()                             = default;
 };
 
-#endif   // OUTPUT_NAMED_ITEM_VALUE_SET_HPP
+#endif   // OUTPUT_NAMED_ITEM_DATA_SET_HPP
diff --git a/src/output/VTKWriter.cpp b/src/output/VTKWriter.cpp
index f65a4acfa89830c54e131539252e80a08de67ff3..0bbdde090daab2cb49a086d359dbbe3b4e53739a 100644
--- a/src/output/VTKWriter.cpp
+++ b/src/output/VTKWriter.cpp
@@ -95,6 +95,23 @@ class VTKWriter::SerializedDataList
     this->add(array);
   }
 
+  template <typename DataT, ItemType item_type, typename ConnectivityT>
+  void
+  add(const ItemArray<DataT, item_type, ConnectivityT>& item_array)
+  {
+    Array<std::remove_const_t<DataT>> array(item_array.numberOfItems() * item_array.sizeOfArrays());
+
+    parallel_for(
+      item_array.numberOfItems(), PUGS_LAMBDA(ItemIdT<item_type> item_id) {
+        const size_t begin = item_id * item_array.sizeOfArrays();
+        for (size_t j = 0; j < item_array.sizeOfArrays(); ++j) {
+          array[begin + j] = item_array[item_id][j];
+        }
+      });
+
+    this->add(array);
+  }
+
   void
   write(std::ostream& os) const
   {
@@ -144,70 +161,48 @@ VTKWriter::_getFilenameVTU(int rank_number) const
   return sout.str();
 }
 
-template <typename DataType>
-void
-VTKWriter::_write_node_pvtu(std::ofstream& os, const std::string& name, const NodeValue<const DataType>&) const
-{
-  os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\"/>\n";
-}
-
-template <size_t N, typename DataType>
+template <typename ItemDataT>
 void
-VTKWriter::_write_node_pvtu(std::ofstream& os,
-                            const std::string& name,
-                            const NodeValue<const TinyVector<N, DataType>>&) const
+VTKWriter::_write_item_pvtu(std::ofstream& os, const std::string& name, const ItemDataT& item_data) const
 {
-  os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N
-     << "\"/>\n";
-}
-
-template <size_t N, typename DataType>
-void
-VTKWriter::_write_node_pvtu(std::ofstream& os,
-                            const std::string& name,
-                            const NodeValue<const TinyMatrix<N, DataType>>&) const
-{
-  os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N * N
-     << "\"/>\n";
-}
-
-template <typename DataType>
-void
-VTKWriter::_write_node_pvtu(std::ofstream&, const std::string&, const CellValue<const DataType>&) const
-{}
+  if constexpr (is_item_value_v<ItemDataT>) {
+    using DataType = std::decay_t<typename ItemDataT::data_type>;
+    if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) {
+      using SubDataType = typename DataType::data_type;
+      os << "<PDataArray type=\"" << VTKType<SubDataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\""
+         << DataType{}.dimension() << "\"/>\n";
+    } else {
+      static_assert(std::is_arithmetic_v<DataType>);
+      os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\"/>\n";
+    }
+  } else {
+    static_assert(is_item_array_v<ItemDataT>);
+    using DataType = typename ItemDataT::data_type;
+    static_assert(std::is_arithmetic_v<DataType>);
 
-template <typename DataType>
-void
-VTKWriter::_write_cell_pvtu(std::ofstream& os, const std::string& name, const CellValue<const DataType>&) const
-{
-  os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\"/>\n";
+    os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\""
+       << item_data.sizeOfArrays() << "\"/>\n";
+  }
 }
 
-template <size_t N, typename DataType>
+template <typename ItemDataT>
 void
-VTKWriter::_write_cell_pvtu(std::ofstream& os,
-                            const std::string& name,
-                            const CellValue<const TinyVector<N, DataType>>&) const
+VTKWriter::_write_node_pvtu(std::ofstream& os, const std::string& name, const ItemDataT& item_data) const
 {
-  os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N
-     << "\"/>\n";
+  if constexpr (ItemDataT::item_t == ItemType::node) {
+    _write_item_pvtu(os, name, item_data);
+  }
 }
 
-template <size_t N, typename DataType>
+template <typename ItemDataT>
 void
-VTKWriter::_write_cell_pvtu(std::ofstream& os,
-                            const std::string& name,
-                            const CellValue<const TinyMatrix<N, DataType>>&) const
+VTKWriter::_write_cell_pvtu(std::ofstream& os, const std::string& name, const ItemDataT& item_data) const
 {
-  os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N * N
-     << "\"/>\n";
+  if constexpr (ItemDataT::item_t == ItemType::cell) {
+    _write_item_pvtu(os, name, item_data);
+  }
 }
 
-template <typename DataType>
-void
-VTKWriter::_write_cell_pvtu(std::ofstream&, const std::string&, const NodeValue<const DataType>&) const
-{}
-
 template <typename DataType>
 struct VTKWriter::VTKType
 {
@@ -230,24 +225,24 @@ template <typename DataType>
 void
 VTKWriter::_write_array(std::ofstream& os,
                         const std::string& name,
-                        const Array<DataType>& item_value,
+                        const Array<DataType>& array,
                         SerializedDataList& serialized_data_list) const
 {
   os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" offset=\""
      << serialized_data_list.offset() << "\" format=\"appended\"/>\n";
-  serialized_data_list.add(item_value);
+  serialized_data_list.add(array);
 }
 
 template <size_t N, typename DataType>
 void
 VTKWriter::_write_array(std::ofstream& os,
                         const std::string& name,
-                        const Array<TinyVector<N, DataType>>& item_value,
+                        const Array<TinyVector<N, DataType>>& array,
                         SerializedDataList& serialized_data_list) const
 {
   os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N
      << "\" offset=\"" << serialized_data_list.offset() << "\" format=\"appended\"/>\n";
-  serialized_data_list.add(item_value);
+  serialized_data_list.add(array);
 }
 
 template <typename DataType>
@@ -294,61 +289,74 @@ VTKWriter::_write_node_value(std::ofstream&,
                              SerializedDataList&) const
 {}
 
-template <typename DataType>
+template <typename ItemDataT>
 void
-VTKWriter::_write_cell_value(std::ofstream& os,
-                             const std::string& name,
-                             const CellValue<const DataType>& item_value,
-                             SerializedDataList& serialized_data_list) const
+VTKWriter::_write_item_data(std::ofstream& os,
+                            const std::string& name,
+                            const ItemDataT& item_data,
+                            SerializedDataList& serialized_data_list) const
 {
-  os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" offset=\""
-     << serialized_data_list.offset() << "\" format=\"appended\"/>\n";
-  serialized_data_list.add(item_value);
+  using DataType = std::decay_t<typename ItemDataT::data_type>;
+  if constexpr (is_item_value_v<ItemDataT>) {
+    if constexpr (std::is_arithmetic_v<DataType>) {
+      os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" offset=\""
+         << serialized_data_list.offset() << "\" format=\"appended\"/>\n";
+    } else if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) {
+      using SubDataType = typename DataType::data_type;
+      os << "<DataArray type=\"" << VTKType<SubDataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\""
+         << DataType{}.dimension() << "\" offset=\"" << serialized_data_list.offset() << "\" format=\"appended\"/>\n";
+    } else {
+      static_assert(std::is_arithmetic_v<DataType>, "unexpected data type");
+    }
+    serialized_data_list.add(item_data);
+
+  } else {
+    static_assert(is_item_array_v<ItemDataT>);
+    if constexpr (std::is_arithmetic_v<DataType>) {
+      os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\""
+         << item_data.sizeOfArrays() << "\" offset=\"" << serialized_data_list.offset() << "\" format=\"appended\"/>\n";
+    } else {
+      static_assert(std::is_arithmetic_v<DataType>, "unexpected data type");
+    }
+    serialized_data_list.add(item_data);
+  }
 }
 
-template <size_t N, typename DataType>
+template <typename ItemDataT>
 void
-VTKWriter::_write_cell_value(std::ofstream& os,
-                             const std::string& name,
-                             const CellValue<const TinyVector<N, DataType>>& item_value,
-                             SerializedDataList& serialized_data_list) const
+VTKWriter::_write_cell_data(std::ofstream& os,
+                            const std::string& name,
+                            const ItemDataT& item_data,
+                            SerializedDataList& serialized_data_list) const
 {
-  os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N
-     << "\" offset=\"" << serialized_data_list.offset() << "\" format=\"appended\"/>\n";
-  serialized_data_list.add(item_value);
+  if constexpr (ItemDataT::item_t == ItemType::cell) {
+    _write_item_data(os, name, item_data, serialized_data_list);
+  }
 }
 
-template <size_t N, typename DataType>
+template <typename ItemDataT>
 void
-VTKWriter::_write_cell_value(std::ofstream& os,
-                             const std::string& name,
-                             const CellValue<const TinyMatrix<N, DataType>>& item_value,
-                             SerializedDataList& serialized_data_list) const
+VTKWriter::_write_node_data(std::ofstream& os,
+                            const std::string& name,
+                            const ItemDataT& item_data,
+                            SerializedDataList& serialized_data_list) const
 {
-  os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N * N
-     << "\" offset=\"" << serialized_data_list.offset() << "\" format=\"appended\"/>\n";
-  serialized_data_list.add(item_value);
+  if constexpr (ItemDataT::item_t == ItemType::node) {
+    _write_item_data(os, name, item_data, serialized_data_list);
+  }
 }
 
-template <typename DataType>
-void
-VTKWriter::_write_cell_value(std::ofstream&,
-                             const std::string&,
-                             const NodeValue<const DataType>&,
-                             SerializedDataList&) const
-{}
-
 template <typename MeshType>
 void
 VTKWriter::_write(const std::shared_ptr<const MeshType>& mesh,
-                  const OutputNamedItemValueSet& given_output_named_item_value_set,
+                  const OutputNamedItemDataSet& given_output_named_item_data_set,
                   double time) const
 {
-  OutputNamedItemValueSet output_named_item_value_set{given_output_named_item_value_set};
+  OutputNamedItemDataSet output_named_item_data_set{given_output_named_item_data_set};
   // Adding basic mesh information
-  output_named_item_value_set.add(NamedItemValue{"cell_number", mesh->connectivity().cellNumber()});
-  output_named_item_value_set.add(NamedItemValue{"node_number", mesh->connectivity().nodeNumber()});
-  output_named_item_value_set.add(NamedItemValue{"cell_center", MeshDataManager::instance().getMeshData(*mesh).xj()});
+  output_named_item_data_set.add(NamedItemData{"cell_number", mesh->connectivity().cellNumber()});
+  output_named_item_data_set.add(NamedItemData{"node_number", mesh->connectivity().nodeNumber()});
+  output_named_item_data_set.add(NamedItemData{"cell_center", MeshDataManager::instance().getMeshData(*mesh).xj()});
 
   if (parallel::rank() == 0) {   // write PVTK file
     std::ofstream fout(_getFilenamePVTU());
@@ -376,14 +384,14 @@ VTKWriter::_write(const std::shared_ptr<const MeshType>& mesh,
     fout << "</PCells>\n";
 
     fout << "<PPointData>\n";
-    for (const auto& [name, item_value_variant] : output_named_item_value_set) {
+    for (const auto& [name, item_value_variant] : output_named_item_data_set) {
       std::visit([&, name = name](auto&& item_value) { return this->_write_node_pvtu(fout, name, item_value); },
                  item_value_variant);
     }
     fout << "</PPointData>\n";
 
     fout << "<PCellData>\n";
-    for (const auto& [name, item_value_variant] : output_named_item_value_set) {
+    for (const auto& [name, item_value_variant] : output_named_item_data_set) {
       std::visit([&, name = name](auto&& item_value) { return this->_write_cell_pvtu(fout, name, item_value); },
                  item_value_variant);
     }
@@ -414,16 +422,16 @@ VTKWriter::_write(const std::shared_ptr<const MeshType>& mesh,
     fout << "<Piece NumberOfPoints=\"" << mesh->numberOfNodes() << "\" NumberOfCells=\"" << mesh->numberOfCells()
          << "\">\n";
     fout << "<CellData>\n";
-    for (const auto& [name, item_value_variant] : output_named_item_value_set) {
+    for (const auto& [name, item_value_variant] : output_named_item_data_set) {
       std::visit([&, name = name](
-                   auto&& item_value) { return this->_write_cell_value(fout, name, item_value, serialize_data_list); },
+                   auto&& item_value) { return this->_write_cell_data(fout, name, item_value, serialize_data_list); },
                  item_value_variant);
     }
     fout << "</CellData>\n";
     fout << "<PointData>\n";
-    for (const auto& [name, item_value_variant] : output_named_item_value_set) {
+    for (const auto& [name, item_value_variant] : output_named_item_data_set) {
       std::visit([&, name = name](
-                   auto&& item_value) { return this->_write_node_value(fout, name, item_value, serialize_data_list); },
+                   auto&& item_value) { return this->_write_node_data(fout, name, item_value, serialize_data_list); },
                  item_value_variant);
     }
     fout << "</PointData>\n";
@@ -617,19 +625,19 @@ VTKWriter::_write(const std::shared_ptr<const MeshType>& mesh,
 void
 VTKWriter::writeMesh(const std::shared_ptr<const IMesh>& mesh) const
 {
-  OutputNamedItemValueSet output_named_item_value_set;
+  OutputNamedItemDataSet output_named_item_data_set;
 
   switch (mesh->dimension()) {
   case 1: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_value_set, 0);
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_data_set, 0);
     break;
   }
   case 2: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_value_set, 0);
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_data_set, 0);
     break;
   }
   case 3: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(mesh), output_named_item_value_set, 0);
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(mesh), output_named_item_data_set, 0);
     break;
   }
   default: {
@@ -644,19 +652,19 @@ VTKWriter::write(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>
 {
   std::shared_ptr mesh = this->_getMesh(named_discrete_function_list);
 
-  OutputNamedItemValueSet output_named_item_value_set = this->_getOutputNamedItemValueSet(named_discrete_function_list);
+  OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_function_list);
 
   switch (mesh->dimension()) {
   case 1: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_value_set, time);
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_data_set, time);
     break;
   }
   case 2: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_value_set, time);
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_data_set, time);
     break;
   }
   case 3: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(mesh), output_named_item_value_set, time);
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(mesh), output_named_item_data_set, time);
     break;
   }
   default: {
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index 7f38119222474af2df46508e96fb27d7a40053f9..7baede4556fb37e2b3b29f449dfbdaf5dfa8b143 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -22,37 +22,14 @@ class VTKWriter : public WriterBase
 
   std::string _getFilenameVTU(int rank_number) const;
 
-  template <typename DataType>
-  void _write_node_pvtu(std::ofstream& os, const std::string& name, const NodeValue<const DataType>&) const;
-
-  template <size_t N, typename DataType>
-  void _write_node_pvtu(std::ofstream& os,
-                        const std::string& name,
-                        const NodeValue<const TinyVector<N, DataType>>&) const;
-
-  template <size_t N, typename DataType>
-  void _write_node_pvtu(std::ofstream& os,
-                        const std::string& name,
-                        const NodeValue<const TinyMatrix<N, DataType>>&) const;
-
-  template <typename DataType>
-  void _write_node_pvtu(std::ofstream&, const std::string&, const CellValue<const DataType>&) const;
-
-  template <typename DataType>
-  void _write_cell_pvtu(std::ofstream& os, const std::string& name, const CellValue<const DataType>&) const;
+  template <typename ItemDataT>
+  void _write_item_pvtu(std::ofstream& os, const std::string& name, const ItemDataT&) const;
 
-  template <size_t N, typename DataType>
-  void _write_cell_pvtu(std::ofstream& os,
-                        const std::string& name,
-                        const CellValue<const TinyVector<N, DataType>>&) const;
+  template <typename ItemDataT>
+  void _write_node_pvtu(std::ofstream& os, const std::string& name, const ItemDataT&) const;
 
-  template <size_t N, typename DataType>
-  void _write_cell_pvtu(std::ofstream& os,
-                        const std::string& name,
-                        const CellValue<const TinyMatrix<N, DataType>>&) const;
-
-  template <typename DataType>
-  void _write_cell_pvtu(std::ofstream&, const std::string&, const NodeValue<const DataType>&) const;
+  template <typename ItemDataT>
+  void _write_cell_pvtu(std::ofstream& os, const std::string& name, const ItemDataT&) const;
 
   template <typename DataType>
   struct VTKType;
@@ -60,13 +37,13 @@ class VTKWriter : public WriterBase
   template <typename DataType>
   void _write_array(std::ofstream& os,
                     const std::string& name,
-                    const Array<DataType>& item_value,
+                    const Array<DataType>& array,
                     SerializedDataList& serialized_data_list) const;
 
   template <size_t N, typename DataType>
   void _write_array(std::ofstream& os,
                     const std::string& name,
-                    const Array<TinyVector<N, DataType>>& item_value,
+                    const Array<TinyVector<N, DataType>>& array,
                     SerializedDataList& serialized_data_list) const;
 
   template <typename DataType>
@@ -99,27 +76,27 @@ class VTKWriter : public WriterBase
                          const CellValue<const DataType>& item_value,
                          SerializedDataList& serialized_data_list) const;
 
-  template <size_t N, typename DataType>
-  void _write_cell_value(std::ofstream& os,
-                         const std::string& name,
-                         const CellValue<const TinyVector<N, DataType>>& item_value,
-                         SerializedDataList& serialized_data_list) const;
+  template <typename ItemDataT>
+  void _write_item_data(std::ofstream& os,
+                        const std::string& name,
+                        const ItemDataT& item_data,
+                        SerializedDataList& serialized_data_list) const;
 
-  template <size_t N, typename DataType>
-  void _write_cell_value(std::ofstream& os,
-                         const std::string& name,
-                         const CellValue<const TinyMatrix<N, DataType>>& item_value,
-                         SerializedDataList& serialized_data_list) const;
+  template <typename ItemDataT>
+  void _write_cell_data(std::ofstream& os,
+                        const std::string& name,
+                        const ItemDataT& item_data,
+                        SerializedDataList& serialized_data_list) const;
 
-  template <typename DataType>
-  void _write_cell_value(std::ofstream&,
-                         const std::string&,
-                         const NodeValue<const DataType>&,
-                         SerializedDataList&) const;
+  template <typename ItemDataT>
+  void _write_node_data(std::ofstream& os,
+                        const std::string& name,
+                        const ItemDataT& item_data,
+                        SerializedDataList& serialized_data_list) const;
 
   template <typename MeshType>
   void _write(const std::shared_ptr<const MeshType>& mesh,
-              const OutputNamedItemValueSet& output_named_item_value_set,
+              const OutputNamedItemDataSet& output_named_item_data_set,
               double time) const;
 
  public:
diff --git a/src/output/WriterBase.cpp b/src/output/WriterBase.cpp
index 4a26605cd330301ddf2fc3346c544228c34f4f9b..40d08b305094c6c33576181a9d730e579656497e 100644
--- a/src/output/WriterBase.cpp
+++ b/src/output/WriterBase.cpp
@@ -3,83 +3,112 @@
 #include <mesh/IMesh.hpp>
 #include <output/OutputNamedItemValueSet.hpp>
 #include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
 #include <scheme/IDiscreteFunction.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 #include <utils/Exceptions.hpp>
 
-template <size_t Dimension, typename DataType>
+template <typename DiscreteFunctionType>
 void
-WriterBase::registerDiscreteFunctionP0(const std::string& name,
-                                       const IDiscreteFunction& i_discrete_function,
-                                       OutputNamedItemValueSet& named_item_value_set)
+WriterBase::_register(const std::string& name,
+                      const DiscreteFunctionType& discrete_function,
+                      OutputNamedItemDataSet& named_item_data_set)
 {
-  const DiscreteFunctionP0<Dimension, DataType>& discrete_function =
-    dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(i_discrete_function);
-  named_item_value_set.add(NamedItemValue{name, discrete_function.cellValues()});
+  if constexpr (DiscreteFunctionType::handled_data_type == IDiscreteFunction::HandledItemDataType::value) {
+    named_item_data_set.add(NamedItemData{name, discrete_function.cellValues()});
+  } else {
+    named_item_data_set.add(NamedItemData{name, discrete_function.cellArrays()});
+  }
 }
 
-template <size_t Dimension>
+template <size_t Dimension, template <size_t DimensionT, typename DataTypeT> typename DiscreteFunctionType>
 void
-WriterBase::registerDiscreteFunctionP0(const std::string& name,
-                                       const IDiscreteFunction& i_discrete_function,
-                                       OutputNamedItemValueSet& named_item_value_set)
+WriterBase::_register(const std::string& name,
+                      const IDiscreteFunction& i_discrete_function,
+                      OutputNamedItemDataSet& named_item_data_set)
 {
   const ASTNodeDataType& data_type = i_discrete_function.dataType();
   switch (data_type) {
   case ASTNodeDataType::bool_t: {
-    registerDiscreteFunctionP0<Dimension, bool>(name, i_discrete_function, named_item_value_set);
+    _register(name, dynamic_cast<const DiscreteFunctionType<Dimension, bool>&>(i_discrete_function),
+              named_item_data_set);
     break;
   }
   case ASTNodeDataType::unsigned_int_t: {
-    registerDiscreteFunctionP0<Dimension, uint64_t>(name, i_discrete_function, named_item_value_set);
+    _register(name, dynamic_cast<const DiscreteFunctionType<Dimension, uint64_t>&>(i_discrete_function),
+              named_item_data_set);
     break;
   }
   case ASTNodeDataType::int_t: {
-    registerDiscreteFunctionP0<Dimension, int64_t>(name, i_discrete_function, named_item_value_set);
+    _register(name, dynamic_cast<const DiscreteFunctionType<Dimension, int64_t>&>(i_discrete_function),
+              named_item_data_set);
     break;
   }
   case ASTNodeDataType::double_t: {
-    registerDiscreteFunctionP0<Dimension, double>(name, i_discrete_function, named_item_value_set);
+    _register(name, dynamic_cast<const DiscreteFunctionType<Dimension, double>&>(i_discrete_function),
+              named_item_data_set);
     break;
   }
   case ASTNodeDataType::vector_t: {
-    switch (data_type.dimension()) {
-    case 1: {
-      registerDiscreteFunctionP0<Dimension, TinyVector<1, double>>(name, i_discrete_function, named_item_value_set);
-      break;
-    }
-    case 2: {
-      registerDiscreteFunctionP0<Dimension, TinyVector<2, double>>(name, i_discrete_function, named_item_value_set);
-      break;
-    }
-    case 3: {
-      registerDiscreteFunctionP0<Dimension, TinyVector<3, double>>(name, i_discrete_function, named_item_value_set);
-      break;
-    }
-    default: {
-      throw UnexpectedError("invalid vector dimension");
-    }
+    if constexpr (DiscreteFunctionType<Dimension, double>::handled_data_type ==
+                  IDiscreteFunction::HandledItemDataType::vector) {
+      throw UnexpectedError("invalid data type for vector data");
+    } else {
+      switch (data_type.dimension()) {
+      case 1: {
+        _register(name,
+                  dynamic_cast<const DiscreteFunctionType<Dimension, TinyVector<1, double>>&>(i_discrete_function),
+                  named_item_data_set);
+        break;
+      }
+      case 2: {
+        _register(name,
+                  dynamic_cast<const DiscreteFunctionType<Dimension, TinyVector<2, double>>&>(i_discrete_function),
+                  named_item_data_set);
+        break;
+      }
+      case 3: {
+        _register(name,
+                  dynamic_cast<const DiscreteFunctionType<Dimension, TinyVector<3, double>>&>(i_discrete_function),
+                  named_item_data_set);
+        break;
+      }
+      default: {
+        throw UnexpectedError("invalid vector dimension");
+      }
+      }
     }
     break;
   }
   case ASTNodeDataType::matrix_t: {
-    Assert(data_type.nbRows() == data_type.nbColumns(), "invalid matrix dimensions");
-    switch (data_type.nbRows()) {
-    case 1: {
-      registerDiscreteFunctionP0<Dimension, TinyMatrix<1, double>>(name, i_discrete_function, named_item_value_set);
-      break;
-    }
-    case 2: {
-      registerDiscreteFunctionP0<Dimension, TinyMatrix<2, double>>(name, i_discrete_function, named_item_value_set);
-      break;
-    }
-    case 3: {
-      registerDiscreteFunctionP0<Dimension, TinyMatrix<3, double>>(name, i_discrete_function, named_item_value_set);
-      break;
-    }
-    default: {
-      throw UnexpectedError("invalid matrix dimension");
-    }
+    if constexpr (DiscreteFunctionType<Dimension, double>::handled_data_type ==
+                  IDiscreteFunction::HandledItemDataType::vector) {
+      throw UnexpectedError("invalid data type for vector data");
+    } else {
+      Assert(data_type.nbRows() == data_type.nbColumns(), "invalid matrix dimensions");
+      switch (data_type.nbRows()) {
+      case 1: {
+        _register(name,
+                  dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<1, double>>&>(i_discrete_function),
+                  named_item_data_set);
+        break;
+      }
+      case 2: {
+        _register(name,
+                  dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<2, double>>&>(i_discrete_function),
+                  named_item_data_set);
+        break;
+      }
+      case 3: {
+        _register(name,
+                  dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<3, double>>&>(i_discrete_function),
+                  named_item_data_set);
+        break;
+      }
+      default: {
+        throw UnexpectedError("invalid matrix dimension");
+      }
+      }
     }
     break;
   }
@@ -89,23 +118,23 @@ WriterBase::registerDiscreteFunctionP0(const std::string& name,
   }
 }
 
+template <template <size_t Dimension, typename DataType> typename DiscreteFunctionType>
 void
-WriterBase::registerDiscreteFunctionP0(const NamedDiscreteFunction& named_discrete_function,
-                                       OutputNamedItemValueSet& named_item_value_set)
+WriterBase::_register(const NamedDiscreteFunction& named_discrete_function, OutputNamedItemDataSet& named_item_data_set)
 {
   const IDiscreteFunction& i_discrete_function = *named_discrete_function.discreteFunction();
   const std::string& name                      = named_discrete_function.name();
   switch (i_discrete_function.mesh()->dimension()) {
   case 1: {
-    registerDiscreteFunctionP0<1>(name, i_discrete_function, named_item_value_set);
+    _register<1, DiscreteFunctionType>(name, i_discrete_function, named_item_data_set);
     break;
   }
   case 2: {
-    registerDiscreteFunctionP0<2>(name, i_discrete_function, named_item_value_set);
+    _register<2, DiscreteFunctionType>(name, i_discrete_function, named_item_data_set);
     break;
   }
   case 3: {
-    registerDiscreteFunctionP0<3>(name, i_discrete_function, named_item_value_set);
+    _register<3, DiscreteFunctionType>(name, i_discrete_function, named_item_data_set);
     break;
   }
   default: {
@@ -135,18 +164,22 @@ WriterBase::_getMesh(
   return mesh;
 }
 
-OutputNamedItemValueSet
-WriterBase::_getOutputNamedItemValueSet(
+OutputNamedItemDataSet
+WriterBase::_getOutputNamedItemDataSet(
   const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const
 {
-  OutputNamedItemValueSet named_item_value_set;
+  OutputNamedItemDataSet named_item_data_set;
 
   for (auto& named_discrete_function : named_discrete_function_list) {
     const IDiscreteFunction& i_discrete_function = *named_discrete_function->discreteFunction();
 
     switch (i_discrete_function.descriptor().type()) {
     case DiscreteFunctionType::P0: {
-      WriterBase::registerDiscreteFunctionP0(*named_discrete_function, named_item_value_set);
+      WriterBase::_register<DiscreteFunctionP0>(*named_discrete_function, named_item_data_set);
+      break;
+    }
+    case DiscreteFunctionType::P0Vector: {
+      WriterBase::_register<DiscreteFunctionP0Vector>(*named_discrete_function, named_item_data_set);
       break;
     }
     default: {
@@ -158,7 +191,7 @@ WriterBase::_getOutputNamedItemValueSet(
     }
   }
 
-  return named_item_value_set;
+  return named_item_data_set;
 }
 
 double
diff --git a/src/output/WriterBase.hpp b/src/output/WriterBase.hpp
index babdec90733416cf9007e55c9652f2c380d8899d..3ba91af2867b50078037a0d699f914e42ed63f3e 100644
--- a/src/output/WriterBase.hpp
+++ b/src/output/WriterBase.hpp
@@ -6,7 +6,7 @@
 #include <string>
 
 class IMesh;
-class OutputNamedItemValueSet;
+class OutputNamedItemDataSet;
 
 class WriterBase : public IWriter
 {
@@ -18,19 +18,20 @@ class WriterBase : public IWriter
   mutable std::vector<double> m_saved_times;
 
  private:
-  template <size_t Dimension, typename DataType>
-  static void registerDiscreteFunctionP0(const std::string& name, const IDiscreteFunction&, OutputNamedItemValueSet&);
+  template <typename DiscreteFunctionType>
+  static void _register(const std::string& name, const DiscreteFunctionType&, OutputNamedItemDataSet&);
 
-  template <size_t Dimension>
-  static void registerDiscreteFunctionP0(const std::string& name, const IDiscreteFunction&, OutputNamedItemValueSet&);
+  template <size_t Dimension, template <size_t DimensionT, typename DataTypeT> typename DiscreteFunctionType>
+  static void _register(const std::string& name, const IDiscreteFunction&, OutputNamedItemDataSet&);
 
-  static void registerDiscreteFunctionP0(const NamedDiscreteFunction&, OutputNamedItemValueSet&);
+  template <template <size_t DimensionT, typename DataTypeT> typename DiscreteFunctionType>
+  static void _register(const NamedDiscreteFunction&, OutputNamedItemDataSet&);
 
  protected:
   std::shared_ptr<const IMesh> _getMesh(
     const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const;
 
-  OutputNamedItemValueSet _getOutputNamedItemValueSet(
+  OutputNamedItemDataSet _getOutputNamedItemDataSet(
     const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const;
 
   virtual void write(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index b7522b9bbd08a90ee1bd95a56747753d902e10ed..b9a24259f7fcf800f4215f7b291115c0e110f41b 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -5,4 +5,5 @@ add_library(
   AcousticSolver.cpp
   DiscreteFunctionInterpoler.cpp
   DiscreteFunctionUtils.cpp
+  DiscreteFunctionVectorInterpoler.cpp
   PerfectGas.cpp)
diff --git a/src/scheme/DiscreteFunctionDescriptorP0.hpp b/src/scheme/DiscreteFunctionDescriptorP0.hpp
index 97bfbbf729d4224dbb7af41bf61d489c5a3080cb..142ecefdd0f6212e532e1a16180679efda848ba5 100644
--- a/src/scheme/DiscreteFunctionDescriptorP0.hpp
+++ b/src/scheme/DiscreteFunctionDescriptorP0.hpp
@@ -1,5 +1,5 @@
-#ifndef DISCRETE_FUNCTION_DESCRIPTOR_P_0HPP
-#define DISCRETE_FUNCTION_DESCRIPTOR_P_0HPP
+#ifndef DISCRETE_FUNCTION_DESCRIPTOR_P0_HPP
+#define DISCRETE_FUNCTION_DESCRIPTOR_P0_HPP
 
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 
@@ -21,4 +21,4 @@ class DiscreteFunctionDescriptorP0 : public IDiscreteFunctionDescriptor
   ~DiscreteFunctionDescriptorP0() noexcept = default;
 };
 
-#endif   // DISCRETE_FUNCTION_DESCRIPTOR_P_0HPP
+#endif   // DISCRETE_FUNCTION_DESCRIPTOR_P0_HPP
diff --git a/src/scheme/DiscreteFunctionDescriptorP0Vector.hpp b/src/scheme/DiscreteFunctionDescriptorP0Vector.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..14cb8ad017137e7018a5317f938b0e6a74689b9e
--- /dev/null
+++ b/src/scheme/DiscreteFunctionDescriptorP0Vector.hpp
@@ -0,0 +1,24 @@
+#ifndef DISCRETE_FUNCTION_DESCRIPTOR_P0_VECTOR_HPP
+#define DISCRETE_FUNCTION_DESCRIPTOR_P0_VECTOR_HPP
+
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+
+class DiscreteFunctionDescriptorP0Vector : public IDiscreteFunctionDescriptor
+{
+ public:
+  DiscreteFunctionType
+  type() const final
+  {
+    return DiscreteFunctionType::P0Vector;
+  }
+
+  DiscreteFunctionDescriptorP0Vector() noexcept = default;
+
+  DiscreteFunctionDescriptorP0Vector(const DiscreteFunctionDescriptorP0Vector&) = default;
+
+  DiscreteFunctionDescriptorP0Vector(DiscreteFunctionDescriptorP0Vector&&) noexcept = default;
+
+  ~DiscreteFunctionDescriptorP0Vector() noexcept = default;
+};
+
+#endif   // DISCRETE_FUNCTION_DESCRIPTOR_P0_VECTOR_HPP
diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp
index 53037ea72b28813cd2b2e619574f8faefa7484f0..d41d98db99a0fa74bbe8e05496d5324f33db0634 100644
--- a/src/scheme/DiscreteFunctionP0.hpp
+++ b/src/scheme/DiscreteFunctionP0.hpp
@@ -17,6 +17,8 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   using data_type = DataType;
   using MeshType  = Mesh<Connectivity<Dimension>>;
 
+  static constexpr HandledItemDataType handled_data_type = HandledItemDataType::value;
+
  private:
   std::shared_ptr<const MeshType> m_mesh;
   CellValue<DataType> m_cell_values;
@@ -50,7 +52,7 @@ class DiscreteFunctionP0 : public IDiscreteFunction
 
   PUGS_FORCEINLINE
   DataType&
-  operator[](const CellId& cell_id) const noexcept(NO_ASSERT)
+  operator[](const CellId cell_id) const noexcept(NO_ASSERT)
   {
     return m_cell_values[cell_id];
   }
@@ -131,15 +133,12 @@ class DiscreteFunctionP0 : public IDiscreteFunction
                                                                                                   mesh_data.xj());
   }
 
-  DiscreteFunctionP0(const std::shared_ptr<const MeshType>& mesh) : m_mesh{mesh}, m_cell_values{mesh->connectivity()}
-  {
-    ;
-  }
+  DiscreteFunctionP0(const std::shared_ptr<const MeshType>& mesh) : m_mesh{mesh}, m_cell_values{mesh->connectivity()} {}
 
   DiscreteFunctionP0(const std::shared_ptr<const MeshType>& mesh, const CellValue<DataType>& cell_value)
     : m_mesh{mesh}, m_cell_values{cell_value}
   {
-    ;
+    Assert(mesh->connectivity().shared_ptr() == cell_value.connectivity_ptr());
   }
 
   DiscreteFunctionP0(const DiscreteFunctionP0&) noexcept = default;
diff --git a/src/scheme/DiscreteFunctionP0Vector.hpp b/src/scheme/DiscreteFunctionP0Vector.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..77ddf70f1d7f14a1c4dcab4cb352d4c6280d1715
--- /dev/null
+++ b/src/scheme/DiscreteFunctionP0Vector.hpp
@@ -0,0 +1,170 @@
+#ifndef DISCRETE_FUNCTION_P0_VECTOR_HPP
+#define DISCRETE_FUNCTION_P0_VECTOR_HPP
+
+#include <scheme/IDiscreteFunction.hpp>
+
+#include <language/utils/InterpolateItemArray.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemArray.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <utils/Exceptions.hpp>
+
+template <size_t Dimension, typename DataType>
+class DiscreteFunctionP0Vector : public IDiscreteFunction
+{
+ public:
+  using data_type = DataType;
+  using MeshType  = Mesh<Connectivity<Dimension>>;
+
+  static constexpr HandledItemDataType handled_data_type = HandledItemDataType::vector;
+
+  static_assert(std::is_arithmetic_v<DataType>, "DiscreteFunctionP0Vector are only defined for arithmetic data type");
+
+ private:
+  std::shared_ptr<const MeshType> m_mesh;
+  CellArray<DataType> m_cell_arrays;
+
+  DiscreteFunctionDescriptorP0Vector m_discrete_function_descriptor;
+
+ public:
+  ASTNodeDataType
+  dataType() const final
+  {
+    return ast_node_data_type_from<data_type>;
+  }
+
+  PUGS_INLINE
+  size_t
+  size() const
+  {
+    return m_cell_arrays.sizeOfArrays();
+  }
+
+  const CellArray<DataType>&
+  cellArrays() const
+  {
+    return m_cell_arrays;
+  }
+
+  std::shared_ptr<const IMesh>
+  mesh() const
+  {
+    return m_mesh;
+  }
+
+  const IDiscreteFunctionDescriptor&
+  descriptor() const final
+  {
+    return m_discrete_function_descriptor;
+  }
+
+  PUGS_FORCEINLINE
+  Array<double>
+  operator[](const CellId cell_id) const noexcept(NO_ASSERT)
+  {
+    return m_cell_arrays[cell_id];
+  }
+
+  friend DiscreteFunctionP0Vector
+  operator+(const DiscreteFunctionP0Vector& f, const DiscreteFunctionP0Vector& g)
+  {
+    Assert(f.mesh() == g.mesh(), "functions are nor defined on the same mesh");
+    Assert(f.size() == g.size(), "P0 vectors have different sizes");
+    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(f.mesh());
+    DiscreteFunctionP0Vector sum(mesh, f.size());
+    parallel_for(
+      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        for (size_t i = 0; i < f.size(); ++i) {
+          sum[cell_id][i] = f[cell_id][i] + g[cell_id][i];
+        }
+      });
+    return sum;
+  }
+
+  friend DiscreteFunctionP0Vector
+  operator-(const DiscreteFunctionP0Vector& f, const DiscreteFunctionP0Vector& g)
+  {
+    Assert(f.mesh() == g.mesh(), "functions are nor defined on the same mesh");
+    Assert(f.size() == g.size(), "P0 vectors have different sizes");
+    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(f.mesh());
+    DiscreteFunctionP0Vector difference(mesh, f.size());
+    parallel_for(
+      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        for (size_t i = 0; i < f.size(); ++i) {
+          difference[cell_id][i] = f[cell_id][i] - g[cell_id][i];
+        }
+      });
+    return difference;
+  }
+
+  template <typename DataType2>
+  friend DiscreteFunctionP0Vector
+  operator*(const DiscreteFunctionP0<Dimension, DataType2>& f, const DiscreteFunctionP0Vector& g)
+  {
+    static_assert(std::is_arithmetic_v<DataType2>, "left operand must be arithmetic");
+
+    Assert(f.mesh() == g.mesh(), "functions are nor defined on the same mesh");
+    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(f.mesh());
+    DiscreteFunctionP0Vector product(mesh, g.size());
+    parallel_for(
+      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        const double f_value = f[cell_id];
+        for (size_t i = 0; i < g.size(); ++i) {
+          product[cell_id][i] = f_value * g[cell_id][i];
+        }
+      });
+    return product;
+  }
+
+  friend DiscreteFunctionP0Vector
+  operator*(double a, const DiscreteFunctionP0Vector& f)
+  {
+    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(f.mesh());
+    DiscreteFunctionP0Vector product(mesh, f.size());
+    parallel_for(
+      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        for (size_t i = 0; i < f.size(); ++i) {
+          product[cell_id][i] = a * f[cell_id][i];
+        }
+      });
+    return product;
+  }
+
+  DiscreteFunctionP0Vector(const std::shared_ptr<const MeshType>& mesh,
+                           const std::vector<FunctionSymbolId>& function_symbol_id_list)
+    : m_mesh(mesh)
+  {
+    using MeshDataType      = MeshData<Dimension>;
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+    m_cell_arrays = InterpolateItemArray<DataType(
+      TinyVector<Dimension>)>::template interpolate<ItemType::cell>(function_symbol_id_list, mesh_data.xj());
+  }
+
+  DiscreteFunctionP0Vector(const std::shared_ptr<const MeshType>& mesh, size_t size)
+    : m_mesh{mesh}, m_cell_arrays{mesh->connectivity(), size}
+  {}
+
+  template <typename DataType2>
+  DiscreteFunctionP0Vector(const std::shared_ptr<const MeshType>& mesh, const CellArray<DataType2>& cell_arrays)
+    : m_mesh{mesh}
+  {
+    static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>);
+    static_assert(std::is_const_v<DataType> or not std::is_const_v<DataType2>);
+
+    Assert(mesh->connectivity_ptr() == cell_arrays.connectivity_ptr());
+
+    m_cell_arrays = cell_arrays;
+  }
+
+  DiscreteFunctionP0Vector(const DiscreteFunctionP0Vector&) noexcept = default;
+  DiscreteFunctionP0Vector(DiscreteFunctionP0Vector&&) noexcept      = default;
+
+  ~DiscreteFunctionP0Vector() = default;
+};
+
+#endif   // DISCRETE_FUNCTION_P0_VECTOR_HPP
diff --git a/src/scheme/DiscreteFunctionType.hpp b/src/scheme/DiscreteFunctionType.hpp
index debf863499df64cd161ac0b7c20a1a5ce4657107..8532a25e9111d40a8fcfbd4d0fe8cb6805d525c9 100644
--- a/src/scheme/DiscreteFunctionType.hpp
+++ b/src/scheme/DiscreteFunctionType.hpp
@@ -1,9 +1,27 @@
 #ifndef DISCRETE_FUNCTION_TYPE_HPP
 #define DISCRETE_FUNCTION_TYPE_HPP
 
+#include <utils/PugsMacros.hpp>
+
+#include <string>
+
 enum class DiscreteFunctionType
 {
-  P0
+  P0,
+  P0Vector
 };
 
+PUGS_INLINE
+std::string
+name(DiscreteFunctionType type)
+{
+  switch (type) {
+  case DiscreteFunctionType::P0:
+    return "P0";
+  case DiscreteFunctionType::P0Vector:
+    return "P0Vector";
+  }
+  return {};
+}
+
 #endif   // DISCRETE_FUNCTION_TYPE_HPP
diff --git a/src/scheme/DiscreteFunctionVectorInterpoler.cpp b/src/scheme/DiscreteFunctionVectorInterpoler.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1a5077377dd1858f5f3c2fa495621598784e8360
--- /dev/null
+++ b/src/scheme/DiscreteFunctionVectorInterpoler.cpp
@@ -0,0 +1,60 @@
+#include <scheme/DiscreteFunctionVectorInterpoler.hpp>
+
+#include <scheme/DiscreteFunctionP0Vector.hpp>
+#include <utils/Exceptions.hpp>
+
+template <size_t Dimension, typename DataType>
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionVectorInterpoler::_interpolate() const
+{
+  std::shared_ptr mesh = std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(m_mesh);
+  return std::make_shared<DiscreteFunctionP0Vector<Dimension, DataType>>(mesh, m_function_id_list);
+}
+
+template <size_t Dimension>
+std::shared_ptr<IDiscreteFunction>
+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();
+
+    switch (data_type) {
+    case ASTNodeDataType::bool_t:
+    case ASTNodeDataType::unsigned_int_t:
+    case ASTNodeDataType::int_t:
+    case ASTNodeDataType::double_t: {
+      break;
+    }
+    default: {
+      std::ostringstream os;
+      os << "vector functions require scalar value type.\n"
+         << "Invalid interpolation value type:" << rang::fgB::red << dataTypeName(data_type) << rang::style::reset;
+      throw NormalError(os.str());
+    }
+    }
+  }
+  return this->_interpolate<Dimension, double>();
+}
+
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionVectorInterpoler::interpolate() const
+{
+  std::shared_ptr<IDiscreteFunction> discrete_function;
+  switch (m_mesh->dimension()) {
+  case 1: {
+    return this->_interpolate<1>();
+  }
+  case 2: {
+    return this->_interpolate<2>();
+  }
+  case 3: {
+    return this->_interpolate<3>();
+  }
+  default: {
+    throw UnexpectedError("invalid dimension");
+  }
+  }
+  return nullptr;
+}
diff --git a/src/scheme/DiscreteFunctionVectorInterpoler.hpp b/src/scheme/DiscreteFunctionVectorInterpoler.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..093f290c40b99dc7c3c8f6db1ab32cbae10dfcb2
--- /dev/null
+++ b/src/scheme/DiscreteFunctionVectorInterpoler.hpp
@@ -0,0 +1,41 @@
+#ifndef DISCRETE_FUNCTION_VECTOR_INTERPOLER_HPP
+#define DISCRETE_FUNCTION_VECTOR_INTERPOLER_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IMesh.hpp>
+#include <scheme/IDiscreteFunction.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+
+#include <memory>
+#include <vector>
+
+class DiscreteFunctionVectorInterpoler
+{
+ private:
+  std::shared_ptr<const IMesh> m_mesh;
+  std::shared_ptr<const IDiscreteFunctionDescriptor> m_discrete_function_descriptor;
+  const std::vector<FunctionSymbolId> m_function_id_list;
+
+  template <size_t Dimension, typename DataType>
+  std::shared_ptr<IDiscreteFunction> _interpolate() const;
+
+  template <size_t Dimension>
+  std::shared_ptr<IDiscreteFunction> _interpolate() const;
+
+ public:
+  std::shared_ptr<IDiscreteFunction> interpolate() const;
+
+  DiscreteFunctionVectorInterpoler(
+    const std::shared_ptr<const IMesh>& mesh,
+    const std::shared_ptr<const IDiscreteFunctionDescriptor>& discrete_function_descriptor,
+    const std::vector<FunctionSymbolId>& function_id_list)
+    : m_mesh{mesh}, m_discrete_function_descriptor{discrete_function_descriptor}, m_function_id_list{function_id_list}
+  {}
+
+  DiscreteFunctionVectorInterpoler(const DiscreteFunctionVectorInterpoler&) = delete;
+  DiscreteFunctionVectorInterpoler(DiscreteFunctionVectorInterpoler&&)      = delete;
+
+  ~DiscreteFunctionVectorInterpoler() = default;
+};
+
+#endif   // DISCRETE_FUNCTION_VECTOR_INTERPOLER_HPP
diff --git a/src/scheme/IDiscreteFunction.hpp b/src/scheme/IDiscreteFunction.hpp
index a8bf77f43efb48b03a8967eb16d81966fde86570..be6ca66278de95265e9a301ebc8f2fe0fa7ee5db 100644
--- a/src/scheme/IDiscreteFunction.hpp
+++ b/src/scheme/IDiscreteFunction.hpp
@@ -10,6 +10,12 @@ class IDiscreteFunctionDescriptor;
 class IDiscreteFunction
 {
  public:
+  enum class HandledItemDataType
+  {
+    value,
+    vector,
+  };
+
   virtual std::shared_ptr<const IMesh> mesh() const             = 0;
   virtual const IDiscreteFunctionDescriptor& descriptor() const = 0;
   virtual ASTNodeDataType dataType() const                      = 0;
diff --git a/src/utils/PugsTraits.hpp b/src/utils/PugsTraits.hpp
index e1dc8b48d6b15fa05ac9b9b9378385e351fadb2b..504d694ad86a08811a4a87e2c8098dab77cd607d 100644
--- a/src/utils/PugsTraits.hpp
+++ b/src/utils/PugsTraits.hpp
@@ -8,11 +8,18 @@
 #include <variant>
 #include <vector>
 
+#include <mesh/ItemType.hpp>
+
 template <size_t N, typename T>
 class TinyVector;
 template <size_t N, typename T>
 class TinyMatrix;
 
+template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+class ItemValue;
+template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+class ItemArray;
+
 // Traits is_trivially_castable
 
 template <typename T>
@@ -101,6 +108,22 @@ inline constexpr bool is_tiny_matrix_v = false;
 template <size_t N, typename T>
 inline constexpr bool is_tiny_matrix_v<TinyMatrix<N, T>> = true;
 
+// Trais is ItemValue
+
+template <typename T>
+constexpr inline bool is_item_value_v = false;
+
+template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+constexpr inline bool is_item_value_v<ItemValue<DataType, item_type, ConnectivityPtr>> = true;
+
+// Trais is ItemValue
+
+template <typename T>
+constexpr inline bool is_item_array_v = false;
+
+template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+constexpr inline bool is_item_array_v<ItemArray<DataType, item_type, ConnectivityPtr>> = true;
+
 // helper to check if a type is part of a variant
 
 template <typename T, typename V>