diff --git a/doc/userdoc.org b/doc/userdoc.org
index 57da5149942fcc80cc51ee41114065ca9e822878..490db865223c800e05ee0ee31b52bd86486f8872 100644
--- a/doc/userdoc.org
+++ b/doc/userdoc.org
@@ -2831,32 +2831,35 @@ The information produced concerns
 
 This type is used to designate kinds of items (cell, face, edge or node).
 
-***** ~item_value~
+***** ~item_value~ and ~item_array~
 
-The type ~item_value~ is an abstract type use to designate arrays of
-values defined on the entities of a ~mesh~. Actually, these values are
-not associated to the mesh itself but to the *connectivity*. The values
-on the entities can be of type ~B~, ~N~, ~Z~, ~R~, ~R^1~, ~R^2~, ~R^3~, ~R^1x1~, ~R^2x2~
-and ~R^3x3~. Entities themselves can be cells, faces, edges or nodes.
+The types ~item_value~ and ~item_array~ are abstract types use to
+designate values or arrays of values defined on each entities of a
+~mesh~. Actually, these set of values (or arrays) are not associated to
+the mesh itself but to the *connectivity*. The values on the entities
+can be of type ~B~, ~N~, ~Z~, ~R~, ~R^1~, ~R^2~, ~R^3~, ~R^1x1~, ~R^2x2~ and
+~R^3x3~. Entities themselves can be cells, faces, edges or nodes.
 
 These variables are used to pass data from one function to another.
 
 #+BEGIN_warning
-By now, no mathematical operation is defined on ~item_value~ variables.
+By now, no mathematical operation is defined on ~item_value~ or
+~item_array~ variables.
 #+END_warning
 
 Moreover, ~item_value~ variables can be post processed. Observe that
 edge or face values cannot be post processed since neither ~VTK~ nor
-~Gnuplot~ can handle these data.
+~Gnuplot~ can handle these data. Also ~item_raray~ can only be post
+processed if array per item contain scalar data ( ~B~, ~N~, ~Z~, ~R~).
 
-***** ~sub_item_value~
+***** ~sub_item_value~ and ~sub_item_arrays~
 
-This abstract type handles values defined on the sub items of the
-items of a ~mesh~. Similarly these values are attached to a *connectivity*
-and not to a mesh. The values associated to the sub items can be of
-type ~B~, ~N~, ~Z~, ~R~, ~R^1~, ~R^2~, ~R^3~, ~R^1x1~, ~R^2x2~ or ~R^3x3~. An example of
-~sub_item_value~ is the $\mathbf{C}_{jr}$ vectors which are defined at each
-node of each cell.
+These abstract type handles values or arrays defined on the sub items
+of the items of a ~mesh~. Similarly these sets of values or arrays are
+attached to a *connectivity* and not to a mesh. The values associated to
+the sub items can be of type ~B~, ~N~, ~Z~, ~R~, ~R^1~, ~R^2~, ~R^3~, ~R^1x1~, ~R^2x2~
+or ~R^3x3~. An example of ~sub_item_value~ is the $\mathbf{C}_{jr}$ vectors
+which are defined at each node of each cell.
 
 These variables are used to pass data from one function to
 another. They cannot be post processed.
diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index af1fca45e7dec5898a12731430358691728c7a2a..1f2d0ce4db829694dca9875593f168ebcac47b6a 100644
--- a/src/language/modules/MeshModule.cpp
+++ b/src/language/modules/MeshModule.cpp
@@ -5,6 +5,7 @@
 #include <language/utils/BinaryOperatorProcessorBuilder.hpp>
 #include <language/utils/BuiltinFunctionEmbedder.hpp>
 #include <language/utils/FunctionTable.hpp>
+#include <language/utils/ItemArrayVariantFunctionInterpoler.hpp>
 #include <language/utils/ItemValueVariantFunctionInterpoler.hpp>
 #include <language/utils/OStream.hpp>
 #include <language/utils/OperatorRepository.hpp>
@@ -17,6 +18,7 @@
 #include <mesh/GmshReader.hpp>
 #include <mesh/IBoundaryDescriptor.hpp>
 #include <mesh/IZoneDescriptor.hpp>
+#include <mesh/ItemArrayVariant.hpp>
 #include <mesh/ItemValueVariant.hpp>
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshRelaxer.hpp>
@@ -25,6 +27,7 @@
 #include <mesh/NamedZoneDescriptor.hpp>
 #include <mesh/NumberedBoundaryDescriptor.hpp>
 #include <mesh/NumberedZoneDescriptor.hpp>
+#include <mesh/SubItemArrayPerItemVariant.hpp>
 #include <mesh/SubItemValuePerItemVariant.hpp>
 #include <utils/Exceptions.hpp>
 
@@ -39,7 +42,9 @@ MeshModule::MeshModule()
   this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const ItemType>>);
 
   this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const ItemValueVariant>>);
+  this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const ItemArrayVariant>>);
   this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const SubItemValuePerItemVariant>>);
+  this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const SubItemArrayPerItemVariant>>);
 
   this->_addBuiltinFunction("cell", std::function(
 
@@ -125,6 +130,17 @@ MeshModule::MeshModule()
 
                               ));
 
+  this->_addBuiltinFunction(
+    "interpolate_array",
+    std::function(
+
+      [](std::shared_ptr<const IMesh> mesh, std::shared_ptr<const ItemType> item_type,
+         const std::vector<FunctionSymbolId>& function_id_list) -> std::shared_ptr<const ItemArrayVariant> {
+        return ItemArrayVariantFunctionInterpoler{mesh, *item_type, function_id_list}.interpolate();
+      }
+
+      ));
+
   this->_addBuiltinFunction("transform", std::function(
 
                                            [](std::shared_ptr<const IMesh> p_mesh,
diff --git a/src/language/modules/MeshModule.hpp b/src/language/modules/MeshModule.hpp
index bf9f0b714b327a1216f7ec3e54fc98a3ec65d70f..580f25aa839e2c1cc691247dd7359d8a6c1aa036 100644
--- a/src/language/modules/MeshModule.hpp
+++ b/src/language/modules/MeshModule.hpp
@@ -29,11 +29,21 @@ template <>
 inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const ItemValueVariant>> =
   ASTNodeDataType::build<ASTNodeDataType::type_id_t>("item_value");
 
+class ItemArrayVariant;
+template <>
+inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const ItemArrayVariant>> =
+  ASTNodeDataType::build<ASTNodeDataType::type_id_t>("item_array");
+
 class SubItemValuePerItemVariant;
 template <>
 inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const SubItemValuePerItemVariant>> =
   ASTNodeDataType::build<ASTNodeDataType::type_id_t>("sub_item_value");
 
+class SubItemArrayPerItemVariant;
+template <>
+inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const SubItemArrayPerItemVariant>> =
+  ASTNodeDataType::build<ASTNodeDataType::type_id_t>("sub_item_array");
+
 class MeshModule : public BuiltinModule
 {
  public:
diff --git a/src/language/modules/WriterModule.cpp b/src/language/modules/WriterModule.cpp
index aa479c0eb94005d59000d98dd460d8fedb89ab21..af7d0cb4b1cf50de79462a56d07f262cf1bf95ad 100644
--- a/src/language/modules/WriterModule.cpp
+++ b/src/language/modules/WriterModule.cpp
@@ -4,12 +4,14 @@
 #include <language/utils/TypeDescriptor.hpp>
 #include <mesh/Connectivity.hpp>
 #include <mesh/GmshReader.hpp>
+#include <mesh/ItemArrayVariant.hpp>
 #include <mesh/ItemValueVariant.hpp>
 #include <mesh/Mesh.hpp>
 #include <output/GnuplotWriter.hpp>
 #include <output/GnuplotWriter1D.hpp>
 #include <output/IWriter.hpp>
 #include <output/NamedDiscreteFunction.hpp>
+#include <output/NamedItemArrayVariant.hpp>
 #include <output/NamedItemValueVariant.hpp>
 #include <output/VTKWriter.hpp>
 #include <scheme/DiscreteFunctionVariant.hpp>
@@ -80,6 +82,7 @@ WriterModule::WriterModule()
                                              }
 
                                              ));
+
   this->_addBuiltinFunction("name_output", std::function(
 
                                              [](std::shared_ptr<const ItemValueVariant> item_value_variant,
@@ -90,6 +93,16 @@ WriterModule::WriterModule()
 
                                              ));
 
+  this->_addBuiltinFunction("name_output", std::function(
+
+                                             [](std::shared_ptr<const ItemArrayVariant> item_array_variant,
+                                                const std::string& name) -> std::shared_ptr<const INamedDiscreteData> {
+                                               return std::make_shared<const NamedItemArrayVariant>(item_array_variant,
+                                                                                                    name);
+                                             }
+
+                                             ));
+
   this->_addBuiltinFunction("write_mesh",
                             std::function(
 
diff --git a/src/language/utils/CMakeLists.txt b/src/language/utils/CMakeLists.txt
index a81ffa8aecdf7af0085296b181f2b3ec4d63238f..db3e7856b758b6a4822d97ac0330958b3b367176 100644
--- a/src/language/utils/CMakeLists.txt
+++ b/src/language/utils/CMakeLists.txt
@@ -28,6 +28,7 @@ add_library(PugsLanguageUtils
   FunctionSymbolId.cpp
   IncDecOperatorRegisterForN.cpp
   IncDecOperatorRegisterForZ.cpp
+  ItemArrayVariantFunctionInterpoler.cpp
   ItemValueVariantFunctionInterpoler.cpp
   OFStream.cpp
   OperatorRepository.cpp
diff --git a/src/language/utils/ItemArrayVariantFunctionInterpoler.cpp b/src/language/utils/ItemArrayVariantFunctionInterpoler.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..02d47d399ed9ec0fdb570fc3877b0b384f7a5429
--- /dev/null
+++ b/src/language/utils/ItemArrayVariantFunctionInterpoler.cpp
@@ -0,0 +1,159 @@
+#include <language/utils/ItemArrayVariantFunctionInterpoler.hpp>
+
+#include <language/utils/InterpolateItemArray.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemArrayVariant.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <utils/Exceptions.hpp>
+
+#include <memory>
+
+template <size_t Dimension, typename DataType, typename ArrayType>
+std::shared_ptr<ItemArrayVariant>
+ItemArrayVariantFunctionInterpoler::_interpolate() const
+{
+  std::shared_ptr p_mesh = std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(m_mesh);
+  using MeshDataType     = MeshData<Dimension>;
+
+  switch (m_item_type) {
+  case ItemType::cell: {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);
+    return std::make_shared<ItemArrayVariant>(
+      InterpolateItemArray<DataType(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(m_function_id_list,
+                                                                                                  mesh_data.xj()));
+  }
+  case ItemType::face: {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);
+    return std::make_shared<ItemArrayVariant>(
+      InterpolateItemArray<DataType(TinyVector<Dimension>)>::template interpolate<ItemType::face>(m_function_id_list,
+                                                                                                  mesh_data.xl()));
+  }
+  case ItemType::edge: {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);
+    return std::make_shared<ItemArrayVariant>(
+      InterpolateItemArray<DataType(TinyVector<Dimension>)>::template interpolate<ItemType::edge>(m_function_id_list,
+                                                                                                  mesh_data.xe()));
+  }
+  case ItemType::node: {
+    return std::make_shared<ItemArrayVariant>(
+      InterpolateItemArray<DataType(TinyVector<Dimension>)>::template interpolate<ItemType::node>(m_function_id_list,
+                                                                                                  p_mesh->xr()));
+  }
+    // LCOV_EXCL_START
+  default: {
+    throw UnexpectedError("invalid item type");
+  }
+    // LCOV_EXCL_STOP
+  }
+}
+
+template <size_t Dimension>
+std::shared_ptr<ItemArrayVariant>
+ItemArrayVariantFunctionInterpoler::_interpolate() const
+{
+  const ASTNodeDataType data_type = [&] {
+    const auto& function0_descriptor = m_function_id_list[0].descriptor();
+    Assert(function0_descriptor.domainMappingNode().children[1]->m_data_type == ASTNodeDataType::typename_t);
+
+    ASTNodeDataType data_type = function0_descriptor.domainMappingNode().children[1]->m_data_type.contentType();
+
+    for (size_t i = 1; i < m_function_id_list.size(); ++i) {
+      const auto& function_descriptor = m_function_id_list[i].descriptor();
+      Assert(function_descriptor.domainMappingNode().children[1]->m_data_type == ASTNodeDataType::typename_t);
+      if (data_type != function_descriptor.domainMappingNode().children[1]->m_data_type.contentType()) {
+        throw NormalError("functions must have the same type");
+      }
+    }
+
+    return data_type;
+  }();
+
+  switch (data_type) {
+  case ASTNodeDataType::bool_t: {
+    return this->_interpolate<Dimension, bool>();
+  }
+  case ASTNodeDataType::unsigned_int_t: {
+    return this->_interpolate<Dimension, uint64_t>();
+  }
+  case ASTNodeDataType::int_t: {
+    return this->_interpolate<Dimension, int64_t>();
+  }
+  case ASTNodeDataType::double_t: {
+    return this->_interpolate<Dimension, double>();
+  }
+  case ASTNodeDataType::vector_t: {
+    switch (data_type.dimension()) {
+    case 1: {
+      return this->_interpolate<Dimension, TinyVector<1>>();
+    }
+    case 2: {
+      return this->_interpolate<Dimension, TinyVector<2>>();
+    }
+    case 3: {
+      return this->_interpolate<Dimension, TinyVector<3>>();
+    }
+      // LCOV_EXCL_START
+    default: {
+      std::ostringstream os;
+      os << "invalid vector dimension " << rang::fgB::red << data_type.dimension() << rang::style::reset;
+
+      throw UnexpectedError(os.str());
+    }
+      // LCOV_EXCL_STOP
+    }
+  }
+  case ASTNodeDataType::matrix_t: {
+    Assert(data_type.numberOfColumns() == data_type.numberOfRows(), "undefined matrix type");
+    switch (data_type.numberOfColumns()) {
+    case 1: {
+      return this->_interpolate<Dimension, TinyMatrix<1>>();
+    }
+    case 2: {
+      return this->_interpolate<Dimension, TinyMatrix<2>>();
+    }
+    case 3: {
+      return this->_interpolate<Dimension, TinyMatrix<3>>();
+    }
+      // LCOV_EXCL_START
+    default: {
+      std::ostringstream os;
+      os << "invalid vector dimension " << rang::fgB::red << data_type.dimension() << rang::style::reset;
+
+      throw UnexpectedError(os.str());
+    }
+      // LCOV_EXCL_STOP
+    }
+  }
+    // LCOV_EXCL_START
+  default: {
+    std::ostringstream os;
+    os << "invalid interpolation array type: " << rang::fgB::red << dataTypeName(data_type) << rang::style::reset;
+
+    throw UnexpectedError(os.str());
+  }
+    // LCOV_EXCL_STOP
+  }
+}
+
+std::shared_ptr<ItemArrayVariant>
+ItemArrayVariantFunctionInterpoler::interpolate() const
+{
+  switch (m_mesh->dimension()) {
+  case 1: {
+    return this->_interpolate<1>();
+  }
+  case 2: {
+    return this->_interpolate<2>();
+  }
+  case 3: {
+    return this->_interpolate<3>();
+  }
+    // LCOV_EXCL_START
+  default: {
+    throw UnexpectedError("invalid dimension");
+  }
+    // LCOV_EXCL_STOP
+  }
+}
diff --git a/src/language/utils/ItemArrayVariantFunctionInterpoler.hpp b/src/language/utils/ItemArrayVariantFunctionInterpoler.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..eac61b5fb067e83c2cf8f132060e68884316cad4
--- /dev/null
+++ b/src/language/utils/ItemArrayVariantFunctionInterpoler.hpp
@@ -0,0 +1,38 @@
+#ifndef ITEM_ARRAY_VARIANT_FUNCTION_INTERPOLER_HPP
+#define ITEM_ARRAY_VARIANT_FUNCTION_INTERPOLER_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IMesh.hpp>
+#include <mesh/IZoneDescriptor.hpp>
+#include <mesh/ItemArrayVariant.hpp>
+#include <mesh/ItemType.hpp>
+
+class ItemArrayVariantFunctionInterpoler
+{
+ private:
+  std::shared_ptr<const IMesh> m_mesh;
+  const ItemType m_item_type;
+  const std::vector<FunctionSymbolId> m_function_id_list;
+
+  template <size_t Dimension, typename DataType, typename ArrayType = DataType>
+  std::shared_ptr<ItemArrayVariant> _interpolate() const;
+
+  template <size_t Dimension>
+  std::shared_ptr<ItemArrayVariant> _interpolate() const;
+
+ public:
+  std::shared_ptr<ItemArrayVariant> interpolate() const;
+
+  ItemArrayVariantFunctionInterpoler(const std::shared_ptr<const IMesh>& mesh,
+                                     const ItemType& item_type,
+                                     const std::vector<FunctionSymbolId>& function_id_list)
+    : m_mesh{mesh}, m_item_type{item_type}, m_function_id_list{function_id_list}
+  {}
+
+  ItemArrayVariantFunctionInterpoler(const ItemArrayVariantFunctionInterpoler&) = delete;
+  ItemArrayVariantFunctionInterpoler(ItemArrayVariantFunctionInterpoler&&)      = delete;
+
+  ~ItemArrayVariantFunctionInterpoler() = default;
+};
+
+#endif   // ITEM_ARRAY_VARIANT_FUNCTION_INTERPOLER_HPP
diff --git a/src/language/utils/ItemValueVariantFunctionInterpoler.cpp b/src/language/utils/ItemValueVariantFunctionInterpoler.cpp
index e10e834e55a126823414316ce5d527716befd761..98dbe1658b45b8b8613892d0cfa72e67ccc67700 100644
--- a/src/language/utils/ItemValueVariantFunctionInterpoler.cpp
+++ b/src/language/utils/ItemValueVariantFunctionInterpoler.cpp
@@ -60,13 +60,13 @@ ItemValueVariantFunctionInterpoler::_interpolate() const
 
   switch (data_type) {
   case ASTNodeDataType::bool_t: {
-    return this->_interpolate<Dimension, bool, double>();
+    return this->_interpolate<Dimension, bool>();
   }
   case ASTNodeDataType::unsigned_int_t: {
-    return this->_interpolate<Dimension, uint64_t, double>();
+    return this->_interpolate<Dimension, uint64_t>();
   }
   case ASTNodeDataType::int_t: {
-    return this->_interpolate<Dimension, int64_t, double>();
+    return this->_interpolate<Dimension, int64_t>();
   }
   case ASTNodeDataType::double_t: {
     return this->_interpolate<Dimension, double>();
diff --git a/src/mesh/ItemArrayVariant.hpp b/src/mesh/ItemArrayVariant.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0284cce20965d6568c6272358163e866e70ed393
--- /dev/null
+++ b/src/mesh/ItemArrayVariant.hpp
@@ -0,0 +1,113 @@
+#ifndef ITEM_ARRAY_VARIANT_HPP
+#define ITEM_ARRAY_VARIANT_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <mesh/ItemArray.hpp>
+#include <utils/Exceptions.hpp>
+
+class ItemArrayVariant
+{
+ private:
+  using Variant = std::variant<NodeArray<const bool>,
+                               NodeArray<const long int>,
+                               NodeArray<const unsigned long int>,
+                               NodeArray<const double>,
+                               NodeArray<const TinyVector<1, double>>,
+                               NodeArray<const TinyVector<2, double>>,
+                               NodeArray<const TinyVector<3, double>>,
+                               NodeArray<const TinyMatrix<1, 1, double>>,
+                               NodeArray<const TinyMatrix<2, 2, double>>,
+                               NodeArray<const TinyMatrix<3, 3, double>>,
+
+                               EdgeArray<const bool>,
+                               EdgeArray<const long int>,
+                               EdgeArray<const unsigned long int>,
+                               EdgeArray<const double>,
+                               EdgeArray<const TinyVector<1, double>>,
+                               EdgeArray<const TinyVector<2, double>>,
+                               EdgeArray<const TinyVector<3, double>>,
+                               EdgeArray<const TinyMatrix<1, 1, double>>,
+                               EdgeArray<const TinyMatrix<2, 2, double>>,
+                               EdgeArray<const TinyMatrix<3, 3, double>>,
+
+                               FaceArray<const bool>,
+                               FaceArray<const long int>,
+                               FaceArray<const unsigned long int>,
+                               FaceArray<const double>,
+                               FaceArray<const TinyVector<1, double>>,
+                               FaceArray<const TinyVector<2, double>>,
+                               FaceArray<const TinyVector<3, double>>,
+                               FaceArray<const TinyMatrix<1, 1, double>>,
+                               FaceArray<const TinyMatrix<2, 2, double>>,
+                               FaceArray<const TinyMatrix<3, 3, double>>,
+
+                               CellArray<const bool>,
+                               CellArray<const long int>,
+                               CellArray<const unsigned long int>,
+                               CellArray<const double>,
+                               CellArray<const TinyVector<1, double>>,
+                               CellArray<const TinyVector<2, double>>,
+                               CellArray<const TinyVector<3, double>>,
+                               CellArray<const TinyMatrix<1, 1, double>>,
+                               CellArray<const TinyMatrix<2, 2, double>>,
+                               CellArray<const TinyMatrix<3, 3, double>>>;
+
+  Variant m_item_array;
+
+ public:
+  PUGS_INLINE
+  const Variant&
+  itemArray() const
+  {
+    return m_item_array;
+  }
+
+  template <typename ItemArrayT>
+  PUGS_INLINE auto
+  get() const
+  {
+    using DataType               = typename ItemArrayT::data_type;
+    constexpr ItemType item_type = ItemArrayT::item_t;
+
+    if constexpr (std::is_same_v<ItemArrayT, ItemArray<DataType, item_type>> or
+                  std::is_same_v<ItemArrayT, ItemArray<const DataType, item_type>> or
+                  std::is_same_v<ItemArrayT, WeakItemArray<DataType, item_type>> or
+                  std::is_same_v<ItemArrayT, WeakItemArray<const DataType, item_type>>) {
+      if (not std::holds_alternative<ItemArray<const DataType, item_type>>(this->m_item_array)) {
+        throw NormalError("invalid ItemArray type");
+      }
+      return std::get<ItemArray<const DataType, item_type>>(this->itemArray());
+    } else {
+      static_assert(std::is_same_v<ItemArrayT, ItemArrayT>, "invalid template argument");
+    }
+  }
+
+  template <typename DataType, ItemType item_type>
+  ItemArrayVariant(const ItemArray<DataType, item_type>& item_array)
+    : m_item_array{ItemArray<const DataType, item_type>{item_array}}
+  {
+    static_assert(std::is_same_v<std::remove_const_t<DataType>, bool> or                         //
+                    std::is_same_v<std::remove_const_t<DataType>, long int> or                   //
+                    std::is_same_v<std::remove_const_t<DataType>, unsigned long int> or          //
+                    std::is_same_v<std::remove_const_t<DataType>, double> or                     //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyVector<1, double>> or      //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyVector<2, double>> or      //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyVector<3, double>> or      //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<1, 1, double>> or   //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<2, 2, double>> or   //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<3, 3, double>>,
+                  "ItemArray with this DataType is not allowed in variant");
+  }
+
+  ItemArrayVariant& operator=(ItemArrayVariant&&) = default;
+  ItemArrayVariant& operator=(const ItemArrayVariant&) = default;
+
+  ItemArrayVariant(const ItemArrayVariant&) = default;
+  ItemArrayVariant(ItemArrayVariant&&)      = default;
+
+  ItemArrayVariant()  = delete;
+  ~ItemArrayVariant() = default;
+};
+
+#endif   // ITEM_ARRAY_VARIANT_HPP
diff --git a/src/mesh/SubItemArrayPerItemVariant.hpp b/src/mesh/SubItemArrayPerItemVariant.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ffb0ea673241236a7e6fd1a3ee0566439abcf765
--- /dev/null
+++ b/src/mesh/SubItemArrayPerItemVariant.hpp
@@ -0,0 +1,202 @@
+#ifndef SUB_ITEM_ARRAY_PER_ITEM_VARIANT_HPP
+#define SUB_ITEM_ARRAY_PER_ITEM_VARIANT_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <mesh/SubItemArrayPerItem.hpp>
+#include <utils/Exceptions.hpp>
+
+class SubItemArrayPerItemVariant
+{
+ private:
+  using Variant = std::variant<NodeArrayPerEdge<const bool>,
+                               NodeArrayPerEdge<const long int>,
+                               NodeArrayPerEdge<const unsigned long int>,
+                               NodeArrayPerEdge<const double>,
+                               NodeArrayPerEdge<const TinyVector<1, double>>,
+                               NodeArrayPerEdge<const TinyVector<2, double>>,
+                               NodeArrayPerEdge<const TinyVector<3, double>>,
+                               NodeArrayPerEdge<const TinyMatrix<1, 1, double>>,
+                               NodeArrayPerEdge<const TinyMatrix<2, 2, double>>,
+                               NodeArrayPerEdge<const TinyMatrix<3, 3, double>>,
+
+                               NodeArrayPerFace<const bool>,
+                               NodeArrayPerFace<const long int>,
+                               NodeArrayPerFace<const unsigned long int>,
+                               NodeArrayPerFace<const double>,
+                               NodeArrayPerFace<const TinyVector<1, double>>,
+                               NodeArrayPerFace<const TinyVector<2, double>>,
+                               NodeArrayPerFace<const TinyVector<3, double>>,
+                               NodeArrayPerFace<const TinyMatrix<1, 1, double>>,
+                               NodeArrayPerFace<const TinyMatrix<2, 2, double>>,
+                               NodeArrayPerFace<const TinyMatrix<3, 3, double>>,
+
+                               NodeArrayPerCell<const bool>,
+                               NodeArrayPerCell<const long int>,
+                               NodeArrayPerCell<const unsigned long int>,
+                               NodeArrayPerCell<const double>,
+                               NodeArrayPerCell<const TinyVector<1, double>>,
+                               NodeArrayPerCell<const TinyVector<2, double>>,
+                               NodeArrayPerCell<const TinyVector<3, double>>,
+                               NodeArrayPerCell<const TinyMatrix<1, 1, double>>,
+                               NodeArrayPerCell<const TinyMatrix<2, 2, double>>,
+                               NodeArrayPerCell<const TinyMatrix<3, 3, double>>,
+
+                               EdgeArrayPerNode<const bool>,
+                               EdgeArrayPerNode<const long int>,
+                               EdgeArrayPerNode<const unsigned long int>,
+                               EdgeArrayPerNode<const double>,
+                               EdgeArrayPerNode<const TinyVector<1, double>>,
+                               EdgeArrayPerNode<const TinyVector<2, double>>,
+                               EdgeArrayPerNode<const TinyVector<3, double>>,
+                               EdgeArrayPerNode<const TinyMatrix<1, 1, double>>,
+                               EdgeArrayPerNode<const TinyMatrix<2, 2, double>>,
+                               EdgeArrayPerNode<const TinyMatrix<3, 3, double>>,
+
+                               EdgeArrayPerFace<const bool>,
+                               EdgeArrayPerFace<const long int>,
+                               EdgeArrayPerFace<const unsigned long int>,
+                               EdgeArrayPerFace<const double>,
+                               EdgeArrayPerFace<const TinyVector<1, double>>,
+                               EdgeArrayPerFace<const TinyVector<2, double>>,
+                               EdgeArrayPerFace<const TinyVector<3, double>>,
+                               EdgeArrayPerFace<const TinyMatrix<1, 1, double>>,
+                               EdgeArrayPerFace<const TinyMatrix<2, 2, double>>,
+                               EdgeArrayPerFace<const TinyMatrix<3, 3, double>>,
+
+                               EdgeArrayPerCell<const bool>,
+                               EdgeArrayPerCell<const long int>,
+                               EdgeArrayPerCell<const unsigned long int>,
+                               EdgeArrayPerCell<const double>,
+                               EdgeArrayPerCell<const TinyVector<1, double>>,
+                               EdgeArrayPerCell<const TinyVector<2, double>>,
+                               EdgeArrayPerCell<const TinyVector<3, double>>,
+                               EdgeArrayPerCell<const TinyMatrix<1, 1, double>>,
+                               EdgeArrayPerCell<const TinyMatrix<2, 2, double>>,
+                               EdgeArrayPerCell<const TinyMatrix<3, 3, double>>,
+
+                               FaceArrayPerNode<const bool>,
+                               FaceArrayPerNode<const long int>,
+                               FaceArrayPerNode<const unsigned long int>,
+                               FaceArrayPerNode<const double>,
+                               FaceArrayPerNode<const TinyVector<1, double>>,
+                               FaceArrayPerNode<const TinyVector<2, double>>,
+                               FaceArrayPerNode<const TinyVector<3, double>>,
+                               FaceArrayPerNode<const TinyMatrix<1, 1, double>>,
+                               FaceArrayPerNode<const TinyMatrix<2, 2, double>>,
+                               FaceArrayPerNode<const TinyMatrix<3, 3, double>>,
+
+                               FaceArrayPerEdge<const bool>,
+                               FaceArrayPerEdge<const long int>,
+                               FaceArrayPerEdge<const unsigned long int>,
+                               FaceArrayPerEdge<const double>,
+                               FaceArrayPerEdge<const TinyVector<1, double>>,
+                               FaceArrayPerEdge<const TinyVector<2, double>>,
+                               FaceArrayPerEdge<const TinyVector<3, double>>,
+                               FaceArrayPerEdge<const TinyMatrix<1, 1, double>>,
+                               FaceArrayPerEdge<const TinyMatrix<2, 2, double>>,
+                               FaceArrayPerEdge<const TinyMatrix<3, 3, double>>,
+
+                               FaceArrayPerCell<const bool>,
+                               FaceArrayPerCell<const long int>,
+                               FaceArrayPerCell<const unsigned long int>,
+                               FaceArrayPerCell<const double>,
+                               FaceArrayPerCell<const TinyVector<1, double>>,
+                               FaceArrayPerCell<const TinyVector<2, double>>,
+                               FaceArrayPerCell<const TinyVector<3, double>>,
+                               FaceArrayPerCell<const TinyMatrix<1, 1, double>>,
+                               FaceArrayPerCell<const TinyMatrix<2, 2, double>>,
+                               FaceArrayPerCell<const TinyMatrix<3, 3, double>>,
+
+                               CellArrayPerNode<const bool>,
+                               CellArrayPerNode<const long int>,
+                               CellArrayPerNode<const unsigned long int>,
+                               CellArrayPerNode<const double>,
+                               CellArrayPerNode<const TinyVector<1, double>>,
+                               CellArrayPerNode<const TinyVector<2, double>>,
+                               CellArrayPerNode<const TinyVector<3, double>>,
+                               CellArrayPerNode<const TinyMatrix<1, 1, double>>,
+                               CellArrayPerNode<const TinyMatrix<2, 2, double>>,
+                               CellArrayPerNode<const TinyMatrix<3, 3, double>>,
+
+                               CellArrayPerEdge<const bool>,
+                               CellArrayPerEdge<const long int>,
+                               CellArrayPerEdge<const unsigned long int>,
+                               CellArrayPerEdge<const double>,
+                               CellArrayPerEdge<const TinyVector<1, double>>,
+                               CellArrayPerEdge<const TinyVector<2, double>>,
+                               CellArrayPerEdge<const TinyVector<3, double>>,
+                               CellArrayPerEdge<const TinyMatrix<1, 1, double>>,
+                               CellArrayPerEdge<const TinyMatrix<2, 2, double>>,
+                               CellArrayPerEdge<const TinyMatrix<3, 3, double>>,
+
+                               CellArrayPerFace<const bool>,
+                               CellArrayPerFace<const long int>,
+                               CellArrayPerFace<const unsigned long int>,
+                               CellArrayPerFace<const double>,
+                               CellArrayPerFace<const TinyVector<1, double>>,
+                               CellArrayPerFace<const TinyVector<2, double>>,
+                               CellArrayPerFace<const TinyVector<3, double>>,
+                               CellArrayPerFace<const TinyMatrix<1, 1, double>>,
+                               CellArrayPerFace<const TinyMatrix<2, 2, double>>,
+                               CellArrayPerFace<const TinyMatrix<3, 3, double>>>;
+
+  Variant m_sub_item_array_per_item;
+
+ public:
+  PUGS_INLINE
+  const Variant&
+  itemArray() const
+  {
+    return m_sub_item_array_per_item;
+  }
+
+  template <typename SubItemArrayPerItemT>
+  PUGS_INLINE auto
+  get() const
+  {
+    using DataType        = typename SubItemArrayPerItemT::data_type;
+    using ItemOfItemTypeT = typename SubItemArrayPerItemT::ItemOfItemType;
+
+    if constexpr (std::is_same_v<SubItemArrayPerItemT, SubItemArrayPerItem<DataType, ItemOfItemTypeT>> or
+                  std::is_same_v<SubItemArrayPerItemT, SubItemArrayPerItem<const DataType, ItemOfItemTypeT>> or
+                  std::is_same_v<SubItemArrayPerItemT, WeakSubItemArrayPerItem<DataType, ItemOfItemTypeT>> or
+                  std::is_same_v<SubItemArrayPerItemT, WeakSubItemArrayPerItem<const DataType, ItemOfItemTypeT>>) {
+      if (not std::holds_alternative<SubItemArrayPerItem<const DataType, ItemOfItemTypeT>>(
+            this->m_sub_item_array_per_item)) {
+        throw NormalError("invalid SubItemArrayPerItem type");
+      }
+      return std::get<SubItemArrayPerItem<const DataType, ItemOfItemTypeT>>(this->m_sub_item_array_per_item);
+    } else {
+      static_assert(std::is_same_v<SubItemArrayPerItemT, SubItemArrayPerItemT>, "invalid template argument");
+    }
+  }
+
+  template <typename DataType, typename ItemOfItemTypeT>
+  SubItemArrayPerItemVariant(const SubItemArrayPerItem<DataType, ItemOfItemTypeT>& sub_item_array_per_item)
+    : m_sub_item_array_per_item{SubItemArrayPerItem<const DataType, ItemOfItemTypeT>{sub_item_array_per_item}}
+  {
+    static_assert(std::is_same_v<std::remove_const_t<DataType>, bool> or                         //
+                    std::is_same_v<std::remove_const_t<DataType>, long int> or                   //
+                    std::is_same_v<std::remove_const_t<DataType>, unsigned long int> or          //
+                    std::is_same_v<std::remove_const_t<DataType>, double> or                     //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyVector<1, double>> or      //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyVector<2, double>> or      //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyVector<3, double>> or      //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<1, 1, double>> or   //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<2, 2, double>> or   //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<3, 3, double>>,
+                  "SubItemArrayPerItem with this DataType is not allowed in variant");
+  }
+
+  SubItemArrayPerItemVariant& operator=(SubItemArrayPerItemVariant&&) = default;
+  SubItemArrayPerItemVariant& operator=(const SubItemArrayPerItemVariant&) = default;
+
+  SubItemArrayPerItemVariant(const SubItemArrayPerItemVariant&) = default;
+  SubItemArrayPerItemVariant(SubItemArrayPerItemVariant&&)      = default;
+
+  SubItemArrayPerItemVariant()  = delete;
+  ~SubItemArrayPerItemVariant() = default;
+};
+
+#endif   // SUB_ITEM_ARRAY_PER_ITEM_VARIANT_HPP
diff --git a/src/output/INamedDiscreteData.hpp b/src/output/INamedDiscreteData.hpp
index 5d22e82c2e9d599cb695970ec9ba7049efdf08bf..b4fae2b34a5bbbd14ca2127b5cad973beb6d9747 100644
--- a/src/output/INamedDiscreteData.hpp
+++ b/src/output/INamedDiscreteData.hpp
@@ -8,6 +8,7 @@ class INamedDiscreteData
  public:
   enum class Type
   {
+    item_array,
     item_value,
     discrete_function
   };
diff --git a/src/output/NamedItemArrayVariant.hpp b/src/output/NamedItemArrayVariant.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4b92d5585d56aa74067c29c7d3e6945262f642e8
--- /dev/null
+++ b/src/output/NamedItemArrayVariant.hpp
@@ -0,0 +1,47 @@
+#ifndef NAMED_ITEM_ARRAY_VARIANT_HPP
+#define NAMED_ITEM_ARRAY_VARIANT_HPP
+
+#include <output/INamedDiscreteData.hpp>
+
+#include <memory>
+#include <string>
+
+class ItemArrayVariant;
+
+class NamedItemArrayVariant final : public INamedDiscreteData
+{
+ private:
+  std::shared_ptr<const ItemArrayVariant> m_item_array_variant;
+  std::string m_name;
+
+ public:
+  Type
+  type() const final
+  {
+    return INamedDiscreteData::Type::item_array;
+  }
+
+  const std::string&
+  name() const final
+  {
+    return m_name;
+  }
+
+  const std::shared_ptr<const ItemArrayVariant>
+  itemArrayVariant() const
+  {
+    return m_item_array_variant;
+  }
+
+  NamedItemArrayVariant(const std::shared_ptr<const ItemArrayVariant>& item_array_variant, const std::string& name)
+    : m_item_array_variant{item_array_variant}, m_name{name}
+  {}
+
+  NamedItemArrayVariant(const NamedItemArrayVariant&) = default;
+  NamedItemArrayVariant(NamedItemArrayVariant&&)      = default;
+
+  NamedItemArrayVariant()  = default;
+  ~NamedItemArrayVariant() = default;
+};
+
+#endif   // NAMED_ITEM_ARRAY_VARIANT_HPP
diff --git a/src/output/WriterBase.cpp b/src/output/WriterBase.cpp
index b5cbee4c412d485fcdd12d6497e4b1fec20706f1..a96da1e4e0a5e1eab52c4c3859dc9f0ae2304f08 100644
--- a/src/output/WriterBase.cpp
+++ b/src/output/WriterBase.cpp
@@ -1,8 +1,10 @@
 #include <output/WriterBase.hpp>
 
 #include <mesh/IMesh.hpp>
+#include <mesh/ItemArrayVariant.hpp>
 #include <mesh/ItemValueVariant.hpp>
 #include <output/NamedDiscreteFunction.hpp>
+#include <output/NamedItemArrayVariant.hpp>
 #include <output/NamedItemValueVariant.hpp>
 #include <output/OutputNamedItemValueSet.hpp>
 #include <scheme/DiscreteFunctionP0.hpp>
@@ -94,6 +96,21 @@ WriterBase::_checkSignature(
         named_item_value.itemValueVariant()->itemValue());
       break;
     }
+    case INamedDiscreteData::Type::item_array: {
+      const NamedItemArrayVariant& named_item_array = dynamic_cast<const NamedItemArrayVariant&>(*named_discrete_data);
+      std::visit(
+        [&](auto&& item_value) {
+          using ItemValueT = std::decay_t<decltype(item_value)>;
+          using DataType   = std::decay_t<typename ItemValueT::data_type>;
+
+          std::ostringstream type_name;
+          type_name << "item_array(" << dataTypeName(ast_node_data_type_from<DataType>) << ')';
+
+          name_type_map[named_discrete_data->name()] = type_name.str();
+        },
+        named_item_array.itemArrayVariant()->itemArray());
+      break;
+    }
     case INamedDiscreteData::Type::discrete_function: {
       const NamedDiscreteFunction& named_discrete_function =
         dynamic_cast<const NamedDiscreteFunction&>(*named_discrete_data);
@@ -107,9 +124,6 @@ WriterBase::_checkSignature(
         named_discrete_function.discreteFunctionVariant()->discreteFunction());
       break;
     }
-    default: {
-      throw UnexpectedError("unexpected discrete data type");
-    }
     }
   }
 
@@ -204,6 +218,17 @@ WriterBase::_getMesh(const std::vector<std::shared_ptr<const INamedDiscreteData>
                    auto&&
                      item_value) { connectivity_set[item_value.connectivity_ptr()] = named_item_value_variant.name(); },
                  named_item_value_variant.itemValueVariant()->itemValue());
+      break;
+    }
+    case INamedDiscreteData::Type::item_array: {
+      const NamedItemArrayVariant& named_item_array_variant =
+        dynamic_cast<const NamedItemArrayVariant&>(*named_discrete_data);
+
+      std::visit([&](
+                   auto&&
+                     item_array) { connectivity_set[item_array.connectivity_ptr()] = named_item_array_variant.name(); },
+                 named_item_array_variant.itemArrayVariant()->itemArray());
+      break;
     }
     }
   }
@@ -271,8 +296,26 @@ WriterBase::_getOutputNamedItemDataSet(
                  item_value_variant.itemValue());
       break;
     }
-    default: {
-      throw UnexpectedError("invalid discrete data type");
+    case INamedDiscreteData::Type::item_array: {
+      const NamedItemArrayVariant& named_item_array_variant =
+        dynamic_cast<const NamedItemArrayVariant&>(*named_discrete_data);
+
+      const std::string& name = named_item_array_variant.name();
+
+      const ItemArrayVariant& item_value_variant = *named_item_array_variant.itemArrayVariant();
+
+      std::visit(
+        [&](auto&& item_array) {
+          using ItemArrayType = std::decay_t<decltype(item_array)>;
+          using DataType      = std::decay_t<typename ItemArrayType::data_type>;
+          if constexpr (std::is_arithmetic_v<DataType>) {
+            named_item_data_set.add(NamedItemData{name, item_array});
+          } else {
+            throw NormalError("can only write item_array containing scalar values");
+          }
+        },
+        item_value_variant.itemArray());
+      break;
     }
     }
   }
@@ -340,8 +383,8 @@ WriterBase::writeOnMesh(const std::shared_ptr<const IMesh>& mesh,
   if (m_period_manager.has_value()) {
     throw NormalError("this writer requires time value");
   } else {
-    _checkMesh(mesh, named_discrete_data_list);
-    _checkConnectivity(mesh, named_discrete_data_list);
+    this->_checkMesh(mesh, named_discrete_data_list);
+    this->_checkConnectivity(mesh, named_discrete_data_list);
     this->_write(*mesh, named_discrete_data_list);
   }
 }
@@ -352,10 +395,12 @@ WriterBase::writeOnMeshIfNeeded(const std::shared_ptr<const IMesh>& mesh,
                                 double time) const
 {
   if (m_period_manager.has_value()) {
-    if (time == m_period_manager->getLastTime())
+    if (time == m_period_manager->getLastTime()) {
       return;   // output already performed
-    _checkMesh(mesh, named_discrete_data_list);
-    _checkConnectivity(mesh, named_discrete_data_list);
+    }
+
+    this->_checkMesh(mesh, named_discrete_data_list);
+    this->_checkConnectivity(mesh, named_discrete_data_list);
     this->_writeAtTime(*mesh, named_discrete_data_list, time);
     m_period_manager->setSaveTime(time);
   } else {
@@ -369,10 +414,11 @@ WriterBase::writeOnMeshForced(const std::shared_ptr<const IMesh>& mesh,
                               double time) const
 {
   if (m_period_manager.has_value()) {
-    if (time == m_period_manager->getLastTime())
+    if (time == m_period_manager->getLastTime()) {
       return;   // output already performed
-    _checkMesh(mesh, named_discrete_data_list);
-    _checkConnectivity(mesh, named_discrete_data_list);
+    }
+    this->_checkMesh(mesh, named_discrete_data_list);
+    this->_checkConnectivity(mesh, named_discrete_data_list);
     this->_writeAtTime(*mesh, named_discrete_data_list, time);
     m_period_manager->setSaveTime(time);
   } else {
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index b5301f36ae05dd0ad9e034a4953f7d7fa95d17ec..9f1f7d9005d6aa8a583283bf2bfda8dc6e23eedd 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -178,6 +178,8 @@ add_executable (mpi_unit_tests
   test_InterpolateItemValue.cpp
   test_ItemArray.cpp
   test_ItemArrayUtils.cpp
+  test_ItemArrayVariant.cpp
+  test_ItemArrayVariantFunctionInterpoler.cpp
   test_ItemValue.cpp
   test_ItemValueUtils.cpp
   test_ItemValueVariant.cpp
@@ -195,6 +197,7 @@ add_executable (mpi_unit_tests
   test_OFStream.cpp
   test_Partitioner.cpp
   test_RandomEngine.cpp
+  test_SubItemArrayPerItemVariant.cpp
   test_SubItemValuePerItem.cpp
   test_SubItemValuePerItemVariant.cpp
   test_SubItemValuePerItemUtils.cpp
diff --git a/tests/test_ItemArrayVariant.cpp b/tests/test_ItemArrayVariant.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f3b529126f6552cff212790f024728b1d3a2cda9
--- /dev/null
+++ b/tests/test_ItemArrayVariant.cpp
@@ -0,0 +1,101 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemArrayVariant.hpp>
+#include <mesh/Mesh.hpp>
+#include <utils/Messenger.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("ItemArrayVariant", "[mesh]")
+{
+  std::shared_ptr mesh = MeshDataBaseForTests::get().hybrid2DMesh();
+
+  const Connectivity<2>& connectivity = *mesh->shared_connectivity();
+
+  using R1 = TinyVector<1>;
+  using R2 = TinyVector<2>;
+  using R3 = TinyVector<3>;
+
+  using R1x1 = TinyMatrix<1>;
+  using R2x2 = TinyMatrix<2>;
+  using R3x3 = TinyMatrix<3>;
+
+  SECTION("NodeArray<double>")
+  {
+    NodeArray<double> node_array{connectivity, 2};
+    ItemArrayVariant v(node_array);
+    REQUIRE_NOTHROW(v.get<NodeArray<const double>>());
+    REQUIRE_THROWS_WITH(v.get<NodeArray<int64_t>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<NodeArray<uint64_t>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<NodeArray<const R1>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<NodeArray<const R2>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<NodeArray<const R3>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<NodeArray<const R1x1>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<NodeArray<const R2x2>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<NodeArray<const R3x3>>(), "error: invalid ItemArray type");
+
+    REQUIRE_THROWS_WITH(v.get<EdgeArray<const double>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<FaceArray<const double>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<CellArray<const double>>(), "error: invalid ItemArray type");
+  }
+
+  SECTION("EdgeArray<R3>")
+  {
+    EdgeArray<R3> node_array{connectivity, 3};
+    ItemArrayVariant v(node_array);
+    REQUIRE_THROWS_WITH(v.get<EdgeArray<int64_t>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArray<uint64_t>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArray<double>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArray<const R1>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArray<const R2>>(), "error: invalid ItemArray type");
+    REQUIRE_NOTHROW(v.get<EdgeArray<const R3>>());
+    REQUIRE_THROWS_WITH(v.get<EdgeArray<const R1x1>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArray<const R2x2>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArray<const R3x3>>(), "error: invalid ItemArray type");
+
+    REQUIRE_THROWS_WITH(v.get<NodeArray<const R3>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<FaceArray<const R3>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<CellArray<const R3>>(), "error: invalid ItemArray type");
+  }
+
+  SECTION("FaceArray<R3x3>")
+  {
+    FaceArray<R3x3> node_array{connectivity, 2};
+    ItemArrayVariant v(node_array);
+    REQUIRE_THROWS_WITH(v.get<FaceArray<int64_t>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<FaceArray<uint64_t>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<FaceArray<double>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<FaceArray<const R1>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<FaceArray<const R2>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<FaceArray<const R3>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<FaceArray<const R1x1>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<FaceArray<const R2x2>>(), "error: invalid ItemArray type");
+    REQUIRE_NOTHROW(v.get<FaceArray<const R3x3>>());
+
+    REQUIRE_THROWS_WITH(v.get<NodeArray<const R3x3>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArray<const R3x3>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<CellArray<const R3x3>>(), "error: invalid ItemArray type");
+  }
+
+  SECTION("CellArray<int64_t>")
+  {
+    CellArray<int64_t> node_array{connectivity, 2};
+    ItemArrayVariant v(node_array);
+    REQUIRE_NOTHROW(v.get<CellArray<const int64_t>>());
+    REQUIRE_THROWS_WITH(v.get<CellArray<const uint64_t>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<CellArray<const double>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<CellArray<const R1>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<CellArray<const R2>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<CellArray<const R3>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<CellArray<const R1x1>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<CellArray<const R2x2>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<CellArray<const R3x3>>(), "error: invalid ItemArray type");
+
+    REQUIRE_THROWS_WITH(v.get<NodeArray<const int64_t>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArray<const int64_t>>(), "error: invalid ItemArray type");
+    REQUIRE_THROWS_WITH(v.get<FaceArray<const int64_t>>(), "error: invalid ItemArray type");
+  }
+}
diff --git a/tests/test_ItemArrayVariantFunctionInterpoler.cpp b/tests/test_ItemArrayVariantFunctionInterpoler.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..640e3e77aa5e6ede05fe0f217d7587f7839e4d56
--- /dev/null
+++ b/tests/test_ItemArrayVariantFunctionInterpoler.cpp
@@ -0,0 +1,648 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <language/ast/ASTBuilder.hpp>
+#include <language/ast/ASTModulesImporter.hpp>
+#include <language/ast/ASTNodeDataTypeBuilder.hpp>
+#include <language/ast/ASTNodeExpressionBuilder.hpp>
+#include <language/ast/ASTNodeFunctionEvaluationExpressionBuilder.hpp>
+#include <language/ast/ASTNodeFunctionExpressionBuilder.hpp>
+#include <language/ast/ASTNodeTypeCleaner.hpp>
+#include <language/ast/ASTSymbolTableBuilder.hpp>
+#include <language/utils/PugsFunctionAdapter.hpp>
+#include <language/utils/SymbolTable.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+
+#include <language/utils/ItemArrayVariantFunctionInterpoler.hpp>
+
+#include <pegtl/string_input.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("ItemArrayVariantFunctionInterpoler", "[scheme]")
+{
+  auto same_item_array = [](auto f, auto g) -> bool {
+    using ItemIdType = typename decltype(f)::index_type;
+    if (f.sizeOfArrays() != g.sizeOfArrays()) {
+      return false;
+    }
+
+    for (ItemIdType item_id = 0; item_id < f.numberOfItems(); ++item_id) {
+      for (size_t i = 0; i < f.sizeOfArrays(); ++i) {
+        if (f[item_id][i] != g[item_id][i]) {
+          return false;
+        }
+      }
+    }
+
+    return true;
+  };
+
+  SECTION("1D")
+  {
+    constexpr size_t Dimension = 1;
+
+    std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+    for (const auto& named_mesh : mesh_list) {
+      SECTION(named_mesh.name())
+      {
+        auto mesh_1d = named_mesh.mesh();
+
+        auto xj = MeshDataManager::instance().getMeshData(*mesh_1d).xj();
+        auto xr = mesh_1d->xr();
+
+        std::string_view data = R"(
+import math;
+let B_scalar_non_linear1_1d: R^1 -> B, x -> (exp(2 * x[0]) + 3 > 4);
+let B_scalar_non_linear2_1d: R^1 -> B, x -> (exp(2 * x[0]) + 3 < 4);
+
+let N_scalar_non_linear1_1d: R^1 -> N, x -> floor(3 * x[0] * x[0] + 2);
+let N_scalar_non_linear2_1d: R^1 -> N, x -> floor(2 * x[0] * x[0]);
+
+let Z_scalar_non_linear1_1d: R^1 -> Z, x -> floor(exp(2 * x[0]) - 1);
+let Z_scalar_non_linear2_1d: R^1 -> Z, x -> floor(cos(2 * x[0]) + 0.5);
+
+let R_scalar_non_linear1_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+let R_scalar_non_linear2_1d: R^1 -> R, x -> 2 * sin(x[0]) + 1;
+let R_scalar_non_linear3_1d: R^1 -> R, x -> x[0] * sin(x[0]);
+
+let R1_non_linear1_1d: R^1 -> R^1, x -> 2 * exp(x[0]);
+let R1_non_linear2_1d: R^1 -> R^1, x -> 2 * exp(x[0])*x[0];
+
+let R2_non_linear_1d: R^1 -> R^2, x -> [2 * exp(x[0]), -3*x[0]];
+
+let R3_non_linear_1d: R^1 -> R^3, x -> [2 * exp(x[0]) + 3, x[0] - 2, 3];
+
+let R1x1_non_linear_1d: R^1 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[0]) + 3);
+
+let R2x2_non_linear_1d: R^1 -> R^2x2, x -> [[2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0])], [3, x[0] * x[0]]];
+
+let R3x3_non_linear_1d: R^1 -> R^3x3, x -> [[2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3], [x[0] * x[0], -4*x[0], 2*x[0]+1], [3, -6*x[0], exp(x[0])]];
+)";
+        TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+        auto ast = ASTBuilder::build(input);
+
+        ASTModulesImporter{*ast};
+        ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+        ASTSymbolTableBuilder{*ast};
+        ASTNodeDataTypeBuilder{*ast};
+
+        ASTNodeTypeCleaner<language::var_declaration>{*ast};
+        ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+        ASTNodeExpressionBuilder{*ast};
+
+        std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+        TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+        position.byte = data.size();   // ensure that variables are declared at this point
+
+        SECTION("B_scalar_non_linear_1d")
+        {
+          auto [i_symbol_f1, found_f1] = symbol_table->find("B_scalar_non_linear1_1d", position);
+          REQUIRE(found_f1);
+          REQUIRE(i_symbol_f1->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function1_symbol_id(std::get<uint64_t>(i_symbol_f1->attributes().value()), symbol_table);
+
+          auto [i_symbol_f2, found_f2] = symbol_table->find("B_scalar_non_linear2_1d", position);
+          REQUIRE(found_f2);
+          REQUIRE(i_symbol_f2->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function2_symbol_id(std::get<uint64_t>(i_symbol_f2->attributes().value()), symbol_table);
+
+          CellArray<bool> cell_array{mesh_1d->connectivity(), 2};
+          parallel_for(
+            cell_array.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_array[cell_id][0]         = std::exp(2 * x[0]) + 3 > 4;
+              cell_array[cell_id][1]         = std::exp(2 * x[0]) + 3 < 4;
+            });
+
+          ItemArrayVariantFunctionInterpoler interpoler(mesh_1d, ItemType::cell,
+                                                        {function1_symbol_id, function2_symbol_id});
+          std::shared_ptr item_data_variant = interpoler.interpolate();
+
+          REQUIRE(same_item_array(cell_array, item_data_variant->get<CellArray<bool>>()));
+        }
+
+        SECTION("N_scalar_non_linear_1d")
+        {
+          auto [i_symbol_f1, found_f1] = symbol_table->find("N_scalar_non_linear1_1d", position);
+          REQUIRE(found_f1);
+          REQUIRE(i_symbol_f1->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function1_symbol_id(std::get<uint64_t>(i_symbol_f1->attributes().value()), symbol_table);
+
+          auto [i_symbol_f2, found_f2] = symbol_table->find("N_scalar_non_linear2_1d", position);
+          REQUIRE(found_f2);
+          REQUIRE(i_symbol_f2->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function2_symbol_id(std::get<uint64_t>(i_symbol_f2->attributes().value()), symbol_table);
+
+          NodeArray<uint64_t> node_array{mesh_1d->connectivity(), 2};
+          parallel_for(
+            node_array.numberOfItems(), PUGS_LAMBDA(const NodeId node_id) {
+              const TinyVector<Dimension>& x = xr[node_id];
+              node_array[node_id][0]         = std::floor(3 * x[0] * x[0] + 2);
+              node_array[node_id][1]         = std::floor(2 * x[0] * x[0]);
+            });
+
+          ItemArrayVariantFunctionInterpoler interpoler(mesh_1d, ItemType::node,
+                                                        {function1_symbol_id, function2_symbol_id});
+          std::shared_ptr item_data_variant = interpoler.interpolate();
+
+          REQUIRE(same_item_array(node_array, item_data_variant->get<NodeArray<uint64_t>>()));
+        }
+
+        SECTION("Z_scalar_non_linear_1d")
+        {
+          auto [i_symbol_f1, found_f1] = symbol_table->find("Z_scalar_non_linear1_1d", position);
+          REQUIRE(found_f1);
+          REQUIRE(i_symbol_f1->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function1_symbol_id(std::get<uint64_t>(i_symbol_f1->attributes().value()), symbol_table);
+
+          auto [i_symbol_f2, found_f2] = symbol_table->find("Z_scalar_non_linear2_1d", position);
+          REQUIRE(found_f2);
+          REQUIRE(i_symbol_f2->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function2_symbol_id(std::get<uint64_t>(i_symbol_f2->attributes().value()), symbol_table);
+
+          CellArray<int64_t> cell_array{mesh_1d->connectivity(), 2};
+          parallel_for(
+            cell_array.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_array[cell_id][0]         = std::floor(std::exp(2 * x[0]) - 1);
+              cell_array[cell_id][1]         = std::floor(cos(2 * x[0]) + 0.5);
+            });
+
+          ItemArrayVariantFunctionInterpoler interpoler(mesh_1d, ItemType::cell,
+                                                        {function1_symbol_id, function2_symbol_id});
+          std::shared_ptr item_array_variant = interpoler.interpolate();
+
+          REQUIRE(same_item_array(cell_array, item_array_variant->get<CellArray<int64_t>>()));
+        }
+
+        SECTION("R_scalar_non_linear_1d")
+        {
+          auto [i_symbol_f1, found_f1] = symbol_table->find("R_scalar_non_linear1_1d", position);
+          REQUIRE(found_f1);
+          REQUIRE(i_symbol_f1->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function1_symbol_id(std::get<uint64_t>(i_symbol_f1->attributes().value()), symbol_table);
+
+          auto [i_symbol_f2, found_f2] = symbol_table->find("R_scalar_non_linear2_1d", position);
+          REQUIRE(found_f2);
+          REQUIRE(i_symbol_f2->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function2_symbol_id(std::get<uint64_t>(i_symbol_f2->attributes().value()), symbol_table);
+
+          auto [i_symbol_f3, found_f3] = symbol_table->find("R_scalar_non_linear3_1d", position);
+          REQUIRE(found_f3);
+          REQUIRE(i_symbol_f3->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function3_symbol_id(std::get<uint64_t>(i_symbol_f3->attributes().value()), symbol_table);
+
+          CellArray<double> cell_array{mesh_1d->connectivity(), 3};
+          parallel_for(
+            cell_array.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+
+              cell_array[cell_id][0] = 2 * std::exp(x[0]) + 3;
+              cell_array[cell_id][1] = 2 * std::sin(x[0]) + 1;
+              cell_array[cell_id][2] = x[0] * std::sin(x[0]);
+            });
+
+          ItemArrayVariantFunctionInterpoler interpoler(mesh_1d, ItemType::cell,
+                                                        {function1_symbol_id, function2_symbol_id,
+                                                         function3_symbol_id});
+          std::shared_ptr item_array_variant = interpoler.interpolate();
+
+          REQUIRE(same_item_array(cell_array, item_array_variant->get<CellArray<double>>()));
+        }
+
+        SECTION("R1_non_linear_1d")
+        {
+          using DataType = TinyVector<1>;
+
+          auto [i_symbol_f1, found_f1] = symbol_table->find("R1_non_linear1_1d", position);
+          REQUIRE(found_f1);
+          REQUIRE(i_symbol_f1->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function1_symbol_id(std::get<uint64_t>(i_symbol_f1->attributes().value()), symbol_table);
+
+          auto [i_symbol_f2, found_f2] = symbol_table->find("R1_non_linear2_1d", position);
+          REQUIRE(found_f2);
+          REQUIRE(i_symbol_f2->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function2_symbol_id(std::get<uint64_t>(i_symbol_f2->attributes().value()), symbol_table);
+
+          CellArray<DataType> cell_array{mesh_1d->connectivity(), 2};
+          parallel_for(
+            cell_array.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+
+              cell_array[cell_id][0] = DataType{2 * std::exp(x[0])};
+              cell_array[cell_id][1] = DataType{2 * std::exp(x[0]) * x[0]};
+            });
+
+          ItemArrayVariantFunctionInterpoler interpoler(mesh_1d, ItemType::cell,
+                                                        {function1_symbol_id, function2_symbol_id});
+          std::shared_ptr item_array_variant = interpoler.interpolate();
+
+          REQUIRE(same_item_array(cell_array, item_array_variant->get<CellArray<DataType>>()));
+        }
+
+        SECTION("R2_non_linear_1d")
+        {
+          using DataType = TinyVector<2>;
+
+          auto [i_symbol, found] = symbol_table->find("R2_non_linear_1d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellArray<DataType> cell_array{mesh_1d->connectivity(), 1};
+          parallel_for(
+            cell_array.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+
+              cell_array[cell_id][0] = DataType{2 * std::exp(x[0]), -3 * x[0]};
+            });
+
+          ItemArrayVariantFunctionInterpoler interpoler(mesh_1d, ItemType::cell, {function_symbol_id});
+          std::shared_ptr item_array_variant = interpoler.interpolate();
+
+          REQUIRE(same_item_array(cell_array, item_array_variant->get<CellArray<DataType>>()));
+        }
+
+        SECTION("R3_non_linear_1d")
+        {
+          using DataType = TinyVector<3>;
+
+          auto [i_symbol, found] = symbol_table->find("R3_non_linear_1d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellArray<DataType> cell_array{mesh_1d->connectivity(), 1};
+          parallel_for(
+            cell_array.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_array[cell_id][0]         = DataType{2 * std::exp(x[0]) + 3, x[0] - 2, 3};
+            });
+
+          ItemArrayVariantFunctionInterpoler interpoler(mesh_1d, ItemType::cell, {function_symbol_id});
+          std::shared_ptr item_array_variant = interpoler.interpolate();
+
+          REQUIRE(same_item_array(cell_array, item_array_variant->get<CellArray<DataType>>()));
+        }
+
+        SECTION("R1x1_non_linear_1d")
+        {
+          using DataType = TinyMatrix<1>;
+
+          auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_1d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellArray<DataType> cell_array{mesh_1d->connectivity(), 1};
+          parallel_for(
+            cell_array.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_array[cell_id][0]         = DataType{2 * std::exp(x[0]) * std::sin(x[0]) + 3};
+            });
+
+          ItemArrayVariantFunctionInterpoler interpoler(mesh_1d, ItemType::cell, {function_symbol_id});
+          std::shared_ptr item_array_variant = interpoler.interpolate();
+
+          REQUIRE(same_item_array(cell_array, item_array_variant->get<CellArray<DataType>>()));
+        }
+
+        SECTION("R2x2_non_linear_1d")
+        {
+          using DataType = TinyMatrix<2>;
+
+          auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_1d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellArray<DataType> cell_array{mesh_1d->connectivity(), 1};
+          parallel_for(
+            cell_array.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_array[cell_id][0] =
+                DataType{2 * std::exp(x[0]) * std::sin(x[0]) + 3, std::sin(x[0] - 2 * x[0]), 3, x[0] * x[0]};
+            });
+
+          ItemArrayVariantFunctionInterpoler interpoler(mesh_1d, ItemType::cell, {function_symbol_id});
+          std::shared_ptr item_array_variant = interpoler.interpolate();
+
+          REQUIRE(same_item_array(cell_array, item_array_variant->get<CellArray<DataType>>()));
+        }
+
+        SECTION("R3x3_non_linear_1d")
+        {
+          using DataType = TinyMatrix<3>;
+
+          auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_1d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellArray<DataType> cell_array{mesh_1d->connectivity(), 1};
+          parallel_for(
+            cell_array.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+
+              cell_array[cell_id][0] = DataType{2 * exp(x[0]) * std::sin(x[0]) + 3,
+                                                std::sin(x[0] - 2 * x[0]),
+                                                3,
+                                                x[0] * x[0],
+                                                -4 * x[0],
+                                                2 * x[0] + 1,
+                                                3,
+                                                -6 * x[0],
+                                                std::exp(x[0])};
+            });
+
+          ItemArrayVariantFunctionInterpoler interpoler(mesh_1d, ItemType::cell, {function_symbol_id});
+          std::shared_ptr item_array_variant = interpoler.interpolate();
+
+          REQUIRE(same_item_array(cell_array, item_array_variant->get<CellArray<DataType>>()));
+        }
+      }
+    }
+  }
+
+  SECTION("2D")
+  {
+    constexpr size_t Dimension = 2;
+
+    std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+    for (const auto& named_mesh : mesh_list) {
+      SECTION(named_mesh.name())
+      {
+        auto mesh_2d = named_mesh.mesh();
+
+        auto xl = MeshDataManager::instance().getMeshData(*mesh_2d).xl();
+
+        std::string_view data = R"(
+import math;
+let B_scalar_non_linear1_2d: R^2 -> B, x -> (exp(2 * x[0])< 2*x[1]);
+let B_scalar_non_linear2_2d: R^2 -> B, x -> (sin(2 * x[0])< x[1]);
+
+let R2_non_linear_2d: R^2 -> R^2, x -> [2 * exp(x[0]), -3*x[1]];
+
+let R3x3_non_linear_2d: R^2 -> R^3x3, x -> [[2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[0]), 3],
+                                           [x[1] * x[0], -4*x[1], 2*x[0]+1],
+                                           [3, -6*x[0], exp(x[1])]];
+  )";
+        TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+        auto ast = ASTBuilder::build(input);
+
+        ASTModulesImporter{*ast};
+        ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+        ASTSymbolTableBuilder{*ast};
+        ASTNodeDataTypeBuilder{*ast};
+
+        ASTNodeTypeCleaner<language::var_declaration>{*ast};
+        ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+        ASTNodeExpressionBuilder{*ast};
+
+        std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+        TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+        position.byte = data.size();   // ensure that variables are declared at this point
+
+        SECTION("B_scalar_non_linear_2d")
+        {
+          using DataType = bool;
+
+          auto [i_symbol_f1, found_f1] = symbol_table->find("B_scalar_non_linear1_2d", position);
+          REQUIRE(found_f1);
+          REQUIRE(i_symbol_f1->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function1_symbol_id(std::get<uint64_t>(i_symbol_f1->attributes().value()), symbol_table);
+
+          auto [i_symbol_f2, found_f2] = symbol_table->find("B_scalar_non_linear2_2d", position);
+          REQUIRE(found_f2);
+          REQUIRE(i_symbol_f2->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function2_symbol_id(std::get<uint64_t>(i_symbol_f2->attributes().value()), symbol_table);
+
+          FaceArray<DataType> face_array{mesh_2d->connectivity(), 2};
+          parallel_for(
+            face_array.numberOfItems(), PUGS_LAMBDA(const FaceId face_id) {
+              const TinyVector<Dimension>& x = xl[face_id];
+              face_array[face_id][0]         = std::exp(2 * x[0]) < 2 * x[1];
+              face_array[face_id][1]         = std::sin(2 * x[0]) < x[1];
+            });
+
+          ItemArrayVariantFunctionInterpoler interpoler(mesh_2d, ItemType::face,
+                                                        {function1_symbol_id, function2_symbol_id});
+          std::shared_ptr item_array_variant = interpoler.interpolate();
+
+          REQUIRE(same_item_array(face_array, item_array_variant->get<FaceArray<DataType>>()));
+        }
+
+        SECTION("R2_non_linear_2d")
+        {
+          using DataType = TinyVector<2>;
+
+          auto [i_symbol, found] = symbol_table->find("R2_non_linear_2d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          FaceArray<DataType> face_array{mesh_2d->connectivity(), 1};
+          parallel_for(
+            face_array.numberOfItems(), PUGS_LAMBDA(const FaceId face_id) {
+              const TinyVector<Dimension>& x = xl[face_id];
+
+              face_array[face_id][0] = DataType{2 * std::exp(x[0]), -3 * x[1]};
+            });
+
+          ItemArrayVariantFunctionInterpoler interpoler(mesh_2d, ItemType::face, {function_symbol_id});
+          std::shared_ptr item_array_variant = interpoler.interpolate();
+
+          REQUIRE(same_item_array(face_array, item_array_variant->get<FaceArray<DataType>>()));
+        }
+
+        SECTION("R3x3_non_linear_2d")
+        {
+          using DataType = TinyMatrix<3>;
+
+          auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_2d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          FaceArray<DataType> face_array{mesh_2d->connectivity(), 1};
+          parallel_for(
+            face_array.numberOfItems(), PUGS_LAMBDA(const FaceId face_id) {
+              const TinyVector<Dimension>& x = xl[face_id];
+
+              face_array[face_id][0] = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3,
+                                                std::sin(x[1] - 2 * x[0]),
+                                                3,
+                                                x[1] * x[0],
+                                                -4 * x[1],
+                                                2 * x[0] + 1,
+                                                3,
+                                                -6 * x[0],
+                                                std::exp(x[1])};
+            });
+
+          ItemArrayVariantFunctionInterpoler interpoler(mesh_2d, ItemType::face, {function_symbol_id});
+          std::shared_ptr item_array_variant = interpoler.interpolate();
+
+          REQUIRE(same_item_array(face_array, item_array_variant->get<FaceArray<DataType>>()));
+        }
+      }
+    }
+  }
+
+  SECTION("3D")
+  {
+    constexpr size_t Dimension = 3;
+
+    std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+    for (const auto& named_mesh : mesh_list) {
+      SECTION(named_mesh.name())
+      {
+        auto mesh_3d = named_mesh.mesh();
+
+        auto xe = MeshDataManager::instance().getMeshData(*mesh_3d).xe();
+
+        std::string_view data = R"(
+import math;
+let R_scalar_non_linear1_3d: R^3 -> R, x -> 2 * exp(x[0]+x[2]) + 3 * x[1];
+let R_scalar_non_linear2_3d: R^3 -> R, x -> 3 * sin(x[0]+x[2]) + 2 * x[1];
+
+let R3_non_linear_3d: R^3 -> R^3, x -> [2 * exp(x[0]) + 3, x[1] - 2, 3 * x[2]];
+let R2x2_non_linear_3d: R^3 -> R^2x2,
+                          x -> [[2 * exp(x[0]) * sin(x[1]) + 3, sin(x[2] - 2 * x[0])],
+                                [3, x[1] * x[0] - x[2]]];
+  )";
+        TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+        auto ast = ASTBuilder::build(input);
+
+        ASTModulesImporter{*ast};
+        ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+        ASTSymbolTableBuilder{*ast};
+        ASTNodeDataTypeBuilder{*ast};
+
+        ASTNodeTypeCleaner<language::var_declaration>{*ast};
+        ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+        ASTNodeExpressionBuilder{*ast};
+
+        std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+        TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+        position.byte = data.size();   // ensure that variables are declared at this point
+
+        SECTION("R_scalar_non_linear_3d")
+        {
+          auto [i_symbol_f1, found_f1] = symbol_table->find("R_scalar_non_linear1_3d", position);
+          REQUIRE(found_f1);
+          REQUIRE(i_symbol_f1->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function1_symbol_id(std::get<uint64_t>(i_symbol_f1->attributes().value()), symbol_table);
+
+          auto [i_symbol_f2, found_f2] = symbol_table->find("R_scalar_non_linear2_3d", position);
+          REQUIRE(found_f2);
+          REQUIRE(i_symbol_f2->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function2_symbol_id(std::get<uint64_t>(i_symbol_f2->attributes().value()), symbol_table);
+
+          EdgeArray<double> edge_array{mesh_3d->connectivity(), 2};
+          parallel_for(
+            edge_array.numberOfItems(), PUGS_LAMBDA(const EdgeId edge_id) {
+              const TinyVector<Dimension>& x = xe[edge_id];
+
+              edge_array[edge_id][0] = 2 * std::exp(x[0] + x[2]) + 3 * x[1];
+              edge_array[edge_id][1] = 3 * std::sin(x[0] + x[2]) + 2 * x[1];
+            });
+
+          ItemArrayVariantFunctionInterpoler interpoler(mesh_3d, ItemType::edge,
+                                                        {function1_symbol_id, function2_symbol_id});
+          std::shared_ptr item_array_variant = interpoler.interpolate();
+
+          REQUIRE(same_item_array(edge_array, item_array_variant->get<EdgeArray<double>>()));
+        }
+
+        SECTION("R3_non_linear_3d")
+        {
+          using DataType = TinyVector<3>;
+
+          auto [i_symbol, found] = symbol_table->find("R3_non_linear_3d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          EdgeArray<DataType> edge_array{mesh_3d->connectivity(), 1};
+          parallel_for(
+            edge_array.numberOfItems(), PUGS_LAMBDA(const EdgeId edge_id) {
+              const TinyVector<Dimension>& x = xe[edge_id];
+
+              edge_array[edge_id][0] = DataType{2 * std::exp(x[0]) + 3, x[1] - 2, 3 * x[2]};
+            });
+
+          ItemArrayVariantFunctionInterpoler interpoler(mesh_3d, ItemType::edge, {function_symbol_id});
+          std::shared_ptr item_array_variant = interpoler.interpolate();
+
+          REQUIRE(same_item_array(edge_array, item_array_variant->get<EdgeArray<DataType>>()));
+        }
+
+        SECTION("R2x2_non_linear_3d")
+        {
+          using DataType = TinyMatrix<2>;
+
+          auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_3d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          EdgeArray<DataType> edge_array{mesh_3d->connectivity(), 1};
+          parallel_for(
+            edge_array.numberOfItems(), PUGS_LAMBDA(const EdgeId edge_id) {
+              const TinyVector<Dimension>& x = xe[edge_id];
+              edge_array[edge_id][0] =
+                DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3, std::sin(x[2] - 2 * x[0]), 3, x[1] * x[0] - x[2]};
+            });
+
+          ItemArrayVariantFunctionInterpoler interpoler(mesh_3d, ItemType::edge, {function_symbol_id});
+          std::shared_ptr item_array_variant = interpoler.interpolate();
+
+          REQUIRE(same_item_array(edge_array, item_array_variant->get<EdgeArray<DataType>>()));
+        }
+      }
+    }
+  }
+}
diff --git a/tests/test_SubItemArrayPerItemVariant.cpp b/tests/test_SubItemArrayPerItemVariant.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6e3413aec062f9c785592884650a706706f21088
--- /dev/null
+++ b/tests/test_SubItemArrayPerItemVariant.cpp
@@ -0,0 +1,204 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/SubItemArrayPerItemVariant.hpp>
+#include <utils/Messenger.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("SubItemArrayPerItemVariant", "[mesh]")
+{
+  std::shared_ptr mesh = MeshDataBaseForTests::get().hybrid3DMesh();
+
+  const Connectivity<3>& connectivity = *mesh->shared_connectivity();
+
+  using R1   = TinyVector<1>;
+  using R2   = TinyVector<2>;
+  using R3   = TinyVector<3>;
+  using R1x1 = TinyMatrix<1>;
+  using R2x2 = TinyMatrix<2>;
+  using R3x3 = TinyMatrix<3>;
+
+  SECTION("NodeArrayPerCell<double>")
+  {
+    NodeArrayPerCell<double> node_array{connectivity, 2};
+    SubItemArrayPerItemVariant v(node_array);
+    REQUIRE_NOTHROW(v.get<NodeArrayPerCell<double>>());
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerCell<int64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerCell<uint64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerCell<R1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerCell<R2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerCell<R3>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerCell<R1x1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerCell<R2x2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerCell<R3x3>>(), "error: invalid SubItemArrayPerItem type");
+  }
+
+  SECTION("NodeArrayPerFace<R1>")
+  {
+    NodeArrayPerFace<R1> node_array{connectivity, 2};
+    SubItemArrayPerItemVariant v(node_array);
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerFace<const double>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerFace<const int64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerFace<const uint64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_NOTHROW(v.get<NodeArrayPerFace<const R1>>());
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerFace<const R2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerFace<const R3>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerFace<const R1x1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerFace<const R2x2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerFace<const R3x3>>(), "error: invalid SubItemArrayPerItem type");
+  }
+
+  SECTION("NodeArrayPerEdge<int64_t>")
+  {
+    NodeArrayPerEdge<int64_t> node_array{connectivity, 3};
+    SubItemArrayPerItemVariant v(node_array);
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerEdge<double>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_NOTHROW(v.get<NodeArrayPerEdge<int64_t>>());
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerEdge<uint64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerEdge<R1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerEdge<R2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerEdge<R3>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerEdge<R1x1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerEdge<R2x2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerEdge<R3x3>>(), "error: invalid SubItemArrayPerItem type");
+  }
+
+  SECTION("EdgeArrayPerCell<R2>")
+  {
+    EdgeArrayPerCell<R2> node_array{connectivity, 3};
+    SubItemArrayPerItemVariant v(node_array);
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerCell<double>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerCell<int64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerCell<uint64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerCell<R1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_NOTHROW(v.get<EdgeArrayPerCell<R2>>());
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerCell<R3>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerCell<R1x1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerCell<R2x2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerCell<R3x3>>(), "error: invalid SubItemArrayPerItem type");
+  }
+
+  SECTION("EdgeArrayPerFace<R1>")
+  {
+    EdgeArrayPerFace<R1> node_array{connectivity, 1};
+    SubItemArrayPerItemVariant v(node_array);
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerFace<const double>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerFace<int64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerFace<uint64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_NOTHROW(v.get<EdgeArrayPerFace<R1>>());
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerFace<R2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerFace<R3>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerFace<R1x1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerFace<R2x2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerFace<R3x3>>(), "error: invalid SubItemArrayPerItem type");
+  }
+
+  SECTION("EdgeArrayPerNode<double>")
+  {
+    EdgeArrayPerNode<double> node_array{connectivity, 2};
+    SubItemArrayPerItemVariant v(node_array);
+    REQUIRE_NOTHROW(v.get<EdgeArrayPerNode<double>>());
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerNode<int64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerNode<uint64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerNode<R1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerNode<R2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerNode<R3>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerNode<R1x1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerNode<R2x2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<EdgeArrayPerNode<R3x3>>(), "error: invalid SubItemArrayPerItem type");
+  }
+
+  SECTION("FaceArrayPerCell<R3x3>")
+  {
+    FaceArrayPerCell<R3x3> node_array{connectivity, 1};
+    SubItemArrayPerItemVariant v(node_array);
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerCell<double>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerCell<int64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerCell<uint64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerCell<R1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerCell<R2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerCell<R3>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerCell<R1x1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerCell<R2x2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_NOTHROW(v.get<FaceArrayPerCell<R3x3>>());
+  }
+
+  SECTION("FaceArrayPerEdge<R2x2>")
+  {
+    FaceArrayPerEdge<R2x2> node_array{connectivity, 2};
+    SubItemArrayPerItemVariant v(node_array);
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerEdge<double>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerEdge<int64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerEdge<uint64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerEdge<R1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerEdge<R2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerEdge<R3>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerEdge<R1x1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_NOTHROW(v.get<FaceArrayPerEdge<R2x2>>());
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerEdge<R3x3>>(), "error: invalid SubItemArrayPerItem type");
+  }
+
+  SECTION("FaceArrayPerNode<uint64_t>")
+  {
+    FaceArrayPerNode<uint64_t> node_array{connectivity, 2};
+    SubItemArrayPerItemVariant v(node_array);
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerNode<double>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerNode<int64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_NOTHROW(v.get<FaceArrayPerNode<uint64_t>>());
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerNode<R1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerNode<R2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerNode<R3>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerNode<R1x1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerNode<R2x2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<FaceArrayPerNode<R3x3>>(), "error: invalid SubItemArrayPerItem type");
+  }
+
+  SECTION("NodeArrayPerCell<R1x1>")
+  {
+    NodeArrayPerCell<R1x1> node_array{connectivity, 1};
+    SubItemArrayPerItemVariant v(node_array);
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerCell<double>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerCell<int64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerCell<uint64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerCell<R1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerCell<R2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerCell<R3>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_NOTHROW(v.get<NodeArrayPerCell<R1x1>>());
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerCell<R2x2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerCell<R3x3>>(), "error: invalid SubItemArrayPerItem type");
+  }
+
+  SECTION("NodeArrayPerFace<R3>")
+  {
+    NodeArrayPerFace<R3> node_array{connectivity, 3};
+    SubItemArrayPerItemVariant v(node_array);
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerFace<double>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerFace<int64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerFace<uint64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerFace<R1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerFace<R2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_NOTHROW(v.get<NodeArrayPerFace<R3>>());
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerFace<R1x1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerFace<R2x2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerFace<R3x3>>(), "error: invalid SubItemArrayPerItem type");
+  }
+
+  SECTION("NodeArrayPerEdge<double>")
+  {
+    NodeArrayPerEdge<double> node_array{connectivity, 2};
+    SubItemArrayPerItemVariant v(node_array);
+    REQUIRE_NOTHROW(v.get<NodeArrayPerEdge<double>>());
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerEdge<int64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerEdge<uint64_t>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerEdge<R1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerEdge<R2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerEdge<R3>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerEdge<R1x1>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerEdge<R2x2>>(), "error: invalid SubItemArrayPerItem type");
+    REQUIRE_THROWS_WITH(v.get<NodeArrayPerEdge<R3x3>>(), "error: invalid SubItemArrayPerItem type");
+  }
+}