diff --git a/src/mesh/ItemArray.hpp b/src/mesh/ItemArray.hpp
index a38f950954084a3527717a2a858d452e3e3cf2ad..ccfff9143c6d7e0f2ff9ff6eeea36533e016c7f3 100644
--- a/src/mesh/ItemArray.hpp
+++ b/src/mesh/ItemArray.hpp
@@ -4,9 +4,8 @@
 #include <mesh/IConnectivity.hpp>
 #include <mesh/ItemId.hpp>
 #include <mesh/ItemType.hpp>
-#include <utils/Array.hpp>
 #include <utils/PugsAssert.hpp>
-#include <utils/SubArray.hpp>
+#include <utils/Table.hpp>
 
 #include <memory>
 
@@ -29,9 +28,7 @@ class ItemArray
 
   ConnectivityPtr m_connectivity_ptr;
 
-  Array<DataType> m_arrays_values;
-
-  size_t m_size_of_arrays;
+  Table<DataType> m_values;
 
   // Allow const std:shared_ptr version to access our data
   friend ItemArray<std::add_const_t<DataType>, item_type, ConnectivitySharedPtr>;
@@ -46,8 +43,7 @@ class ItemArray
     ItemArray<std::remove_const_t<DataType>, item_type, ConnectivityPtr> image;
 
     image.m_connectivity_ptr = source.m_connectivity_ptr;
-    image.m_arrays_values    = copy(source.m_arrays_values);
-    image.m_size_of_arrays   = source.m_size_of_arrays;
+    image.m_values           = copy(source.m_values);
     return image;
   }
 
@@ -69,34 +65,27 @@ class ItemArray
     }
   }
 
-  PUGS_INLINE
-  size_t
-  numberOfValues() const noexcept(NO_ASSERT)
-  {
-    Assert(this->isBuilt());
-    return m_arrays_values.size();
-  }
-
   PUGS_INLINE
   void
   fill(const DataType& data) const noexcept
   {
     static_assert(not std::is_const_v<DataType>, "Cannot modify ItemArray of const");
-    m_arrays_values.fill(data);
+    m_values.fill(data);
   }
 
   // Following Kokkos logic, these classes are view and const view does allow
   // changes in data
   PUGS_FORCEINLINE
-  SubArray<DataType>
+  Array<DataType>
   operator[](const ItemId& i) const noexcept(NO_ASSERT)
   {
     Assert(this->isBuilt());
-    return SubArray{m_arrays_values, i * m_size_of_arrays, m_size_of_arrays};
+    Assert(i < this->numberOfItems());
+    return m_values.row(i);
   }
 
   template <typename IndexType>
-  SubArray<DataType>
+  Array<DataType>
   operator[](const IndexType&) const noexcept(NO_ASSERT)
   {
     static_assert(std::is_same_v<IndexType, ItemId>, "ItemArray must be indexed by ItemId");
@@ -107,7 +96,7 @@ class ItemArray
   numberOfItems() const noexcept(NO_ASSERT)
   {
     Assert(this->isBuilt());
-    return m_connectivity_ptr->template numberOf<item_type>();
+    return m_values.nbRows();
   }
 
   PUGS_INLINE
@@ -115,12 +104,12 @@ class ItemArray
   sizeOfArrays() const
   {
     Assert(this->isBuilt());
-    return m_size_of_arrays;
+    return m_values.nbColumns();
   }
 
   template <typename DataType2>
   PUGS_INLINE ItemArray&
-  operator=(const Array<DataType2>& arrays) noexcept(NO_ASSERT)
+  operator=(const Table<DataType2>& table) noexcept(NO_ASSERT)
   {
     // ensures that DataType is the same as source DataType2
     static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
@@ -129,11 +118,13 @@ class ItemArray
     static_assert(((std::is_const_v<DataType2> and std::is_const_v<DataType>) or not std::is_const_v<DataType2>),
                   "Cannot assign ItemArray of const to ItemArray of non-const");
 
-    Assert((arrays.size() == 0) or this->isBuilt(), "Cannot assign array of arrays to a non-built ItemArray\n");
+    Assert((table.nbRows() * table.nbColumns() == 0) or this->isBuilt(),
+           "Cannot assign array of arrays to a non-built ItemArray\n");
 
-    Assert(m_arrays_values.size() == arrays.size(), "Cannot assign an array of arrays of a different size\n");
+    Assert(this->numberOfItems() == table.nbRows(), "Cannot assign a table of a different dimensions\n");
+    Assert(this->sizeOfArrays() == table.nbColumns(), "Cannot assign a table of a different dimensions\n");
 
-    m_arrays_values = arrays;
+    m_values = table;
 
     return *this;
   }
@@ -149,8 +140,7 @@ class ItemArray
     static_assert(((std::is_const_v<DataType2> and std::is_const_v<DataType>) or not std::is_const_v<DataType2>),
                   "Cannot assign ItemArray of const to ItemArray of non-const");
 
-    m_arrays_values  = array_per_item.m_arrays_values;
-    m_size_of_arrays = array_per_item.m_size_of_arrays;
+    m_values = array_per_item.m_values;
 
     if constexpr (std::is_same_v<ConnectivityPtr, ConnectivitySharedPtr> and
                   std::is_same_v<ConnectivityPtr2, ConnectivityWeakPtr>) {
@@ -174,13 +164,10 @@ class ItemArray
 
   PUGS_INLINE
   ItemArray(const IConnectivity& connectivity, size_t size_of_array) noexcept
-    : m_connectivity_ptr{connectivity.shared_ptr()},
-      m_arrays_values{connectivity.numberOf<item_type>() * size_of_array},
-      m_size_of_arrays{size_of_array}
+    : m_connectivity_ptr{connectivity.shared_ptr()}, m_values{connectivity.numberOf<item_type>(), size_of_array}
   {
     static_assert(not std::is_const_v<DataType>, "Cannot allocate ItemArray of const data: only view is "
                                                  "supported");
-    ;
   }
 
   PUGS_INLINE
diff --git a/src/utils/Array.hpp b/src/utils/Array.hpp
index 560f60ce7015c06fa68c80121b43f86179a05e83..1fe0aacaac25b7541acab4dc532ab93e7f65070e 100644
--- a/src/utils/Array.hpp
+++ b/src/utils/Array.hpp
@@ -39,6 +39,11 @@ class [[nodiscard]] Array
   template <typename DataType2, typename... RT>
   friend PUGS_INLINE Array<DataType2> encapsulate(const Kokkos::View<DataType2*, RT...>& values);
 
+  template <typename DataType2>
+  friend PUGS_INLINE Array<DataType2> subArray(const Array<DataType2>& array,
+                                               typename Array<DataType2>::index_type begin,
+                                               typename Array<DataType2>::index_type size);
+
   PUGS_INLINE DataType& operator[](index_type i) const noexcept(NO_ASSERT)
   {
     Assert(i < m_values.extent(0));
@@ -50,9 +55,7 @@ class [[nodiscard]] Array
   {
     static_assert(not std::is_const<DataType>(), "Cannot modify Array of const");
 
-    // could consider to use std::fill
-    parallel_for(
-      this->size(), PUGS_LAMBDA(index_type i) { m_values[i] = data; });
+    Kokkos::deep_copy(m_values, data);
   }
 
   template <typename DataType2>
@@ -109,6 +112,17 @@ encapsulate(const Kokkos::View<DataType*, RT...>& values)
   return array;
 }
 
+template <typename DataType>
+PUGS_INLINE Array<DataType>
+subArray(const Array<DataType>& array,
+         typename Array<DataType>::index_type begin,
+         typename Array<DataType>::index_type size)
+{
+  Assert(begin < array.size());
+  Assert(begin + size <= array.size());
+  return encapsulate(Kokkos::View<DataType*>(array.m_values, std::make_pair(begin, begin + size)));
+}
+
 // map, multimap, unordered_map and stack cannot be copied this way
 template <typename Container>
 PUGS_INLINE Array<std::remove_const_t<typename Container::value_type>>
diff --git a/src/utils/Table.hpp b/src/utils/Table.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..edb0712afe00ec962415db1a4ec4235f33672711
--- /dev/null
+++ b/src/utils/Table.hpp
@@ -0,0 +1,126 @@
+#ifndef TABLE_HPP
+#define TABLE_HPP
+
+#include <utils/PugsMacros.hpp>
+#include <utils/PugsUtils.hpp>
+
+#include <utils/PugsAssert.hpp>
+
+#include <Kokkos_CopyViews.hpp>
+#include <algorithm>
+
+template <typename DataType>
+class [[nodiscard]] Table
+{
+ public:
+  using data_type  = DataType;
+  using index_type = size_t;
+
+ private:
+  Kokkos::View<DataType**> m_values;
+
+  // Allows const version to access our data
+  friend Table<std::add_const_t<DataType>>;
+
+ public:
+  PUGS_INLINE size_t nbRows() const noexcept
+  {
+    return m_values.extent(0);
+  }
+
+  PUGS_INLINE size_t nbColumns() const noexcept
+  {
+    return m_values.extent(1);
+  }
+
+  Array<DataType> row(index_type i) const
+  {
+    return encapsulate(Kokkos::View<DataType*>(m_values, i, Kokkos::ALL()));
+  }
+
+  Array<DataType> column(index_type j) const
+  {
+    return encapsulate(Kokkos::View<DataType*>(m_values, Kokkos::ALL(), j));
+  }
+
+  friend PUGS_INLINE Table<std::remove_const_t<DataType>> copy(const Table<DataType>& source)
+  {
+    Table<std::remove_const_t<DataType>> image(source.nbRows(), source.nbColumns());
+    Kokkos::deep_copy(image.m_values, source.m_values);
+
+    return image;
+  }
+
+  template <typename DataType2, typename... RT>
+  friend PUGS_INLINE Table<DataType2> encapsulate(const Kokkos::View<DataType2**, RT...>& values);
+
+  PUGS_INLINE DataType& operator()(index_type i, index_type j) const noexcept(NO_ASSERT)
+  {
+    Assert(i < this->nbRows());
+    Assert(j < this->nbColumns());
+    return m_values(i, j);
+  }
+
+  PUGS_INLINE
+  void fill(const DataType& data) const
+  {
+    static_assert(not std::is_const<DataType>(), "Cannot modify Table of const");
+
+    Kokkos::deep_copy(m_values, data);
+  }
+
+  template <typename DataType2>
+  PUGS_INLINE Table& operator=(const Table<DataType2>& table) noexcept
+  {
+    // ensures that DataType is the same as source DataType2
+    static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
+                  "Cannot assign Table of different type");
+    // ensures that const is not lost through copy
+    static_assert(((std::is_const<DataType2>() and std::is_const<DataType>()) or not std::is_const<DataType2>()),
+                  "Cannot assign Table of const to Table of non-const");
+    m_values = table.m_values;
+    return *this;
+  }
+
+  PUGS_INLINE
+  Table& operator=(const Table&) = default;
+
+  PUGS_INLINE
+  Table& operator=(Table&&) = default;
+
+  PUGS_INLINE
+  explicit Table(size_t nb_lines, size_t nb_columns) : m_values("anonymous", nb_lines, nb_columns)
+  {
+    static_assert(not std::is_const<DataType>(), "Cannot allocate Table of const data: only view is "
+                                                 "supported");
+  }
+
+  PUGS_INLINE
+  Table() = default;
+
+  PUGS_INLINE
+  Table(const Table&) = default;
+
+  template <typename DataType2>
+  PUGS_INLINE Table(const Table<DataType2>& table) noexcept
+  {
+    this->operator=(table);
+  }
+
+  PUGS_INLINE
+  Table(Table &&) = default;
+
+  PUGS_INLINE
+  ~Table() = default;
+};
+
+template <typename DataType, typename... RT>
+PUGS_INLINE Table<DataType>
+encapsulate(const Kokkos::View<DataType**, RT...>& values)
+{
+  Table<DataType> table;
+  table.m_values = values;
+  return table;
+}
+
+#endif   // TABLE_HPP
diff --git a/tests/test_ItemArray.cpp b/tests/test_ItemArray.cpp
index 2177d83c2df0dd70674803e7098d0f6b694b9685..8cf05df5df90f3886771905f9a4fe2b727698150 100644
--- a/tests/test_ItemArray.cpp
+++ b/tests/test_ItemArray.cpp
@@ -50,11 +50,6 @@ TEST_CASE("ItemArray", "[mesh]")
     REQUIRE(face_value.numberOfItems() == node_value.numberOfItems());
     REQUIRE(cell_value.numberOfItems() + 1 == node_value.numberOfItems());
 
-    REQUIRE(node_value.numberOfValues() == 3 * node_value.numberOfItems());
-    REQUIRE(edge_value.numberOfValues() == 3 * edge_value.numberOfItems());
-    REQUIRE(face_value.numberOfValues() == 3 * face_value.numberOfItems());
-    REQUIRE(cell_value.numberOfValues() == 3 * cell_value.numberOfItems());
-
     REQUIRE(node_value.sizeOfArrays() == 3);
     REQUIRE(edge_value.sizeOfArrays() == 3);
     REQUIRE(face_value.sizeOfArrays() == 3);
@@ -83,11 +78,6 @@ TEST_CASE("ItemArray", "[mesh]")
 
     REQUIRE(edge_value.numberOfItems() == face_value.numberOfItems());
 
-    REQUIRE(node_value.numberOfValues() == 2 * node_value.numberOfItems());
-    REQUIRE(edge_value.numberOfValues() == 2 * edge_value.numberOfItems());
-    REQUIRE(face_value.numberOfValues() == 2 * face_value.numberOfItems());
-    REQUIRE(cell_value.numberOfValues() == 2 * cell_value.numberOfItems());
-
     REQUIRE(node_value.sizeOfArrays() == 2);
     REQUIRE(edge_value.sizeOfArrays() == 2);
     REQUIRE(face_value.sizeOfArrays() == 2);
@@ -114,11 +104,6 @@ TEST_CASE("ItemArray", "[mesh]")
     FaceArray<int> face_value{connectivity, 3};
     CellArray<int> cell_value{connectivity, 3};
 
-    REQUIRE(node_value.numberOfValues() == 3 * node_value.numberOfItems());
-    REQUIRE(edge_value.numberOfValues() == 3 * edge_value.numberOfItems());
-    REQUIRE(face_value.numberOfValues() == 3 * face_value.numberOfItems());
-    REQUIRE(cell_value.numberOfValues() == 3 * cell_value.numberOfItems());
-
     REQUIRE(node_value.sizeOfArrays() == 3);
     REQUIRE(edge_value.sizeOfArrays() == 3);
     REQUIRE(face_value.sizeOfArrays() == 3);
@@ -132,26 +117,29 @@ TEST_CASE("ItemArray", "[mesh]")
 
     CellArray<size_t> cell_array{connectivity, 3};
 
-    Array<size_t> array{cell_array.numberOfValues()};
-    for (size_t i = 0; i < array.size(); ++i) {
-      array[i] = i;
+    Table<size_t> table{cell_array.numberOfItems(), cell_array.sizeOfArrays()};
+    {
+      size_t k = 0;
+      for (size_t i = 0; i < table.nbRows(); ++i) {
+        for (size_t j = 0; j < table.nbColumns(); ++j) {
+          table(i, j) = k++;
+        }
+      }
     }
+    cell_array = table;
 
-    cell_array = array;
-
-    auto is_same = [](const CellArray<size_t>& cell_array, const Array<size_t>& array) {
+    auto is_same = [](const CellArray<size_t>& cell_array, const Table<size_t>& table) {
       bool is_same = true;
-      size_t k     = 0;
       for (CellId cell_id = 0; cell_id < cell_array.numberOfItems(); ++cell_id) {
-        SubArray sub_array = cell_array[cell_id];
-        for (size_t i = 0; i < sub_array.size(); ++i, ++k) {
-          is_same &= (sub_array[i] == array[k]);
+        Array sub_array = cell_array[cell_id];
+        for (size_t i = 0; i < sub_array.size(); ++i) {
+          is_same &= (sub_array[i] == table(cell_id, i));
         }
       }
       return is_same;
     };
 
-    REQUIRE(is_same(cell_array, array));
+    REQUIRE(is_same(cell_array, table));
   }
 
   SECTION("copy")
@@ -159,7 +147,7 @@ TEST_CASE("ItemArray", "[mesh]")
     auto is_same = [](const auto& cell_array, int value) {
       bool is_same = true;
       for (CellId cell_id = 0; cell_id < cell_array.numberOfItems(); ++cell_id) {
-        SubArray sub_array = cell_array[cell_id];
+        Array sub_array = cell_array[cell_id];
         for (size_t i = 0; i < sub_array.size(); ++i) {
           is_same &= (sub_array[i] == value);
         }
@@ -174,7 +162,8 @@ TEST_CASE("ItemArray", "[mesh]")
     cell_array.fill(parallel::rank());
 
     CellArray<const int> cell_array_const_view{cell_array};
-    REQUIRE(cell_array.numberOfValues() == cell_array_const_view.numberOfValues());
+    REQUIRE(cell_array.numberOfItems() == cell_array_const_view.numberOfItems());
+    REQUIRE(cell_array.sizeOfArrays() == cell_array_const_view.sizeOfArrays());
     REQUIRE(is_same(cell_array_const_view, static_cast<std::int64_t>(parallel::rank())));
 
     CellArray<const int> const_cell_array;
@@ -248,7 +237,7 @@ TEST_CASE("ItemArray", "[mesh]")
 
       CellArray<size_t> cell_array{connectivity, 2};
 
-      Array<size_t> values{3 + cell_array.numberOfValues()};
+      Table<size_t> values{3, connectivity.numberOfCells() + 3};
       REQUIRE_THROWS_AS(cell_array = values, AssertError);
     }
   }
diff --git a/tests/test_ItemArrayUtils.cpp b/tests/test_ItemArrayUtils.cpp
index dc5611278a166cccdb9225e1eee41de9da6c078a..a8c4b7fbf50eb3c4bb4f9af89926eedcd66c5992 100644
--- a/tests/test_ItemArrayUtils.cpp
+++ b/tests/test_ItemArrayUtils.cpp
@@ -28,7 +28,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
         auto node_number = connectivity.nodeNumber();
 
         for (NodeId i_node = 0; i_node < mesh_1d.numberOfNodes(); ++i_node) {
-          SubArray array      = weak_node_array[i_node];
+          Array array         = weak_node_array[i_node];
           const size_t number = node_number[i_node];
           for (size_t i = 0; i < array.size(); ++i) {
             array[i] = number + (parallel::rank() + 1) * i;
@@ -46,7 +46,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
           bool is_synchronized = (parallel::size() > 1);
           bool is_valid        = true;
           for (NodeId i_node = 0; i_node < mesh_1d.numberOfNodes(); ++i_node) {
-            SubArray array      = node_array[i_node];
+            Array array         = node_array[i_node];
             const size_t number = node_number[i_node];
             const size_t owner  = node_owner[i_node];
             if (node_is_owned[i_node]) {
@@ -71,7 +71,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
           bool is_synchronized = true;
           for (NodeId i_node = 0; i_node < mesh_1d.numberOfNodes(); ++i_node) {
-            SubArray array      = node_array[i_node];
+            Array array         = node_array[i_node];
             const size_t number = node_number[i_node];
             const size_t owner  = node_owner[i_node];
             for (size_t i = 0; i < array.size(); ++i) {
@@ -89,7 +89,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
         auto edge_number = connectivity.edgeNumber();
 
         for (EdgeId i_edge = 0; i_edge < mesh_1d.numberOfEdges(); ++i_edge) {
-          SubArray array      = weak_edge_array[i_edge];
+          Array array         = weak_edge_array[i_edge];
           const size_t number = edge_number[i_edge];
           for (size_t i = 0; i < array.size(); ++i) {
             array[i] = number + (parallel::rank() + 1) * i;
@@ -107,7 +107,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
           bool is_synchronized = (parallel::size() > 1);
           bool is_valid        = true;
           for (EdgeId i_edge = 0; i_edge < mesh_1d.numberOfEdges(); ++i_edge) {
-            SubArray array      = edge_array[i_edge];
+            Array array         = edge_array[i_edge];
             const size_t number = edge_number[i_edge];
             const size_t owner  = edge_owner[i_edge];
             if (edge_is_owned[i_edge]) {
@@ -132,7 +132,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
           bool is_synchronized = true;
           for (EdgeId i_edge = 0; i_edge < mesh_1d.numberOfEdges(); ++i_edge) {
-            SubArray array      = edge_array[i_edge];
+            Array array         = edge_array[i_edge];
             const size_t number = edge_number[i_edge];
             const size_t owner  = edge_owner[i_edge];
             for (size_t i = 0; i < array.size(); ++i) {
@@ -150,7 +150,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
         auto face_number = connectivity.faceNumber();
 
         for (FaceId i_face = 0; i_face < mesh_1d.numberOfFaces(); ++i_face) {
-          SubArray array      = weak_face_array[i_face];
+          Array array         = weak_face_array[i_face];
           const size_t number = face_number[i_face];
           for (size_t i = 0; i < array.size(); ++i) {
             array[i] = number + (parallel::rank() + 1) * i;
@@ -168,7 +168,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
           bool is_synchronized = (parallel::size() > 1);
           bool is_valid        = true;
           for (FaceId i_face = 0; i_face < mesh_1d.numberOfFaces(); ++i_face) {
-            SubArray array      = face_array[i_face];
+            Array array         = face_array[i_face];
             const size_t number = face_number[i_face];
             const size_t owner  = face_owner[i_face];
             if (face_is_owned[i_face]) {
@@ -193,7 +193,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
           bool is_synchronized = true;
           for (FaceId i_face = 0; i_face < mesh_1d.numberOfFaces(); ++i_face) {
-            SubArray array      = face_array[i_face];
+            Array array         = face_array[i_face];
             const size_t number = face_number[i_face];
             const size_t owner  = face_owner[i_face];
             for (size_t i = 0; i < array.size(); ++i) {
@@ -211,7 +211,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
         auto cell_number = connectivity.cellNumber();
 
         for (CellId i_cell = 0; i_cell < mesh_1d.numberOfCells(); ++i_cell) {
-          SubArray array      = weak_cell_array[i_cell];
+          Array array         = weak_cell_array[i_cell];
           const size_t number = cell_number[i_cell];
           for (size_t i = 0; i < array.size(); ++i) {
             array[i] = number + (parallel::rank() + 1) * i;
@@ -229,7 +229,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
           bool is_synchronized = (parallel::size() > 1);
           bool is_valid        = true;
           for (CellId i_cell = 0; i_cell < mesh_1d.numberOfCells(); ++i_cell) {
-            SubArray array      = cell_array[i_cell];
+            Array array         = cell_array[i_cell];
             const size_t number = cell_number[i_cell];
             const size_t owner  = cell_owner[i_cell];
             if (cell_is_owned[i_cell]) {
@@ -254,7 +254,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
           bool is_synchronized = true;
           for (CellId i_cell = 0; i_cell < mesh_1d.numberOfCells(); ++i_cell) {
-            SubArray array      = cell_array[i_cell];
+            Array array         = cell_array[i_cell];
             const size_t number = cell_number[i_cell];
             const size_t owner  = cell_owner[i_cell];
             for (size_t i = 0; i < array.size(); ++i) {
@@ -278,7 +278,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
         auto node_number = connectivity.nodeNumber();
 
         for (NodeId i_node = 0; i_node < mesh_2d.numberOfNodes(); ++i_node) {
-          SubArray array      = weak_node_array[i_node];
+          Array array         = weak_node_array[i_node];
           const size_t number = node_number[i_node];
           for (size_t i = 0; i < array.size(); ++i) {
             array[i] = number + (parallel::rank() + 1) * i;
@@ -296,7 +296,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
           bool is_synchronized = (parallel::size() > 1);
           bool is_valid        = true;
           for (NodeId i_node = 0; i_node < mesh_2d.numberOfNodes(); ++i_node) {
-            SubArray array      = node_array[i_node];
+            Array array         = node_array[i_node];
             const size_t number = node_number[i_node];
             const size_t owner  = node_owner[i_node];
             if (node_is_owned[i_node]) {
@@ -321,7 +321,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
           bool is_synchronized = true;
           for (NodeId i_node = 0; i_node < mesh_2d.numberOfNodes(); ++i_node) {
-            SubArray array      = node_array[i_node];
+            Array array         = node_array[i_node];
             const size_t number = node_number[i_node];
             const size_t owner  = node_owner[i_node];
             for (size_t i = 0; i < array.size(); ++i) {
@@ -339,7 +339,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
         auto edge_number = connectivity.edgeNumber();
 
         for (EdgeId i_edge = 0; i_edge < mesh_2d.numberOfEdges(); ++i_edge) {
-          SubArray array      = weak_edge_array[i_edge];
+          Array array         = weak_edge_array[i_edge];
           const size_t number = edge_number[i_edge];
           for (size_t i = 0; i < array.size(); ++i) {
             array[i] = number + (parallel::rank() + 1) * i;
@@ -357,7 +357,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
           bool is_synchronized = (parallel::size() > 1);
           bool is_valid        = true;
           for (EdgeId i_edge = 0; i_edge < mesh_2d.numberOfEdges(); ++i_edge) {
-            SubArray array      = edge_array[i_edge];
+            Array array         = edge_array[i_edge];
             const size_t number = edge_number[i_edge];
             const size_t owner  = edge_owner[i_edge];
             if (edge_is_owned[i_edge]) {
@@ -382,7 +382,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
           bool is_synchronized = true;
           for (EdgeId i_edge = 0; i_edge < mesh_2d.numberOfEdges(); ++i_edge) {
-            SubArray array      = edge_array[i_edge];
+            Array array         = edge_array[i_edge];
             const size_t number = edge_number[i_edge];
             const size_t owner  = edge_owner[i_edge];
             for (size_t i = 0; i < array.size(); ++i) {
@@ -400,7 +400,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
         auto face_number = connectivity.faceNumber();
 
         for (FaceId i_face = 0; i_face < mesh_2d.numberOfFaces(); ++i_face) {
-          SubArray array      = weak_face_array[i_face];
+          Array array         = weak_face_array[i_face];
           const size_t number = face_number[i_face];
           for (size_t i = 0; i < array.size(); ++i) {
             array[i] = number + (parallel::rank() + 1) * i;
@@ -418,7 +418,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
           bool is_synchronized = (parallel::size() > 1);
           bool is_valid        = true;
           for (FaceId i_face = 0; i_face < mesh_2d.numberOfFaces(); ++i_face) {
-            SubArray array      = face_array[i_face];
+            Array array         = face_array[i_face];
             const size_t number = face_number[i_face];
             const size_t owner  = face_owner[i_face];
             if (face_is_owned[i_face]) {
@@ -443,7 +443,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
           bool is_synchronized = true;
           for (FaceId i_face = 0; i_face < mesh_2d.numberOfFaces(); ++i_face) {
-            SubArray array      = face_array[i_face];
+            Array array         = face_array[i_face];
             const size_t number = face_number[i_face];
             const size_t owner  = face_owner[i_face];
             for (size_t i = 0; i < array.size(); ++i) {
@@ -461,7 +461,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
         auto cell_number = connectivity.cellNumber();
 
         for (CellId i_cell = 0; i_cell < mesh_2d.numberOfCells(); ++i_cell) {
-          SubArray array      = weak_cell_array[i_cell];
+          Array array         = weak_cell_array[i_cell];
           const size_t number = cell_number[i_cell];
           for (size_t i = 0; i < array.size(); ++i) {
             array[i] = number + (parallel::rank() + 1) * i;
@@ -479,7 +479,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
           bool is_synchronized = (parallel::size() > 1);
           bool is_valid        = true;
           for (CellId i_cell = 0; i_cell < mesh_2d.numberOfCells(); ++i_cell) {
-            SubArray array      = cell_array[i_cell];
+            Array array         = cell_array[i_cell];
             const size_t number = cell_number[i_cell];
             const size_t owner  = cell_owner[i_cell];
             if (cell_is_owned[i_cell]) {
@@ -504,7 +504,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
           bool is_synchronized = true;
           for (CellId i_cell = 0; i_cell < mesh_2d.numberOfCells(); ++i_cell) {
-            SubArray array      = cell_array[i_cell];
+            Array array         = cell_array[i_cell];
             const size_t number = cell_number[i_cell];
             const size_t owner  = cell_owner[i_cell];
             for (size_t i = 0; i < array.size(); ++i) {
@@ -528,7 +528,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
         auto node_number = connectivity.nodeNumber();
 
         for (NodeId i_node = 0; i_node < mesh_3d.numberOfNodes(); ++i_node) {
-          SubArray array      = weak_node_array[i_node];
+          Array array         = weak_node_array[i_node];
           const size_t number = node_number[i_node];
           for (size_t i = 0; i < array.size(); ++i) {
             array[i] = number + (parallel::rank() + 1) * i;
@@ -546,7 +546,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
           bool is_synchronized = (parallel::size() > 1);
           bool is_valid        = true;
           for (NodeId i_node = 0; i_node < mesh_3d.numberOfNodes(); ++i_node) {
-            SubArray array      = node_array[i_node];
+            Array array         = node_array[i_node];
             const size_t number = node_number[i_node];
             const size_t owner  = node_owner[i_node];
             if (node_is_owned[i_node]) {
@@ -571,7 +571,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
           bool is_synchronized = true;
           for (NodeId i_node = 0; i_node < mesh_3d.numberOfNodes(); ++i_node) {
-            SubArray array      = node_array[i_node];
+            Array array         = node_array[i_node];
             const size_t number = node_number[i_node];
             const size_t owner  = node_owner[i_node];
             for (size_t i = 0; i < array.size(); ++i) {
@@ -589,7 +589,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
         auto edge_number = connectivity.edgeNumber();
 
         for (EdgeId i_edge = 0; i_edge < mesh_3d.numberOfEdges(); ++i_edge) {
-          SubArray array      = weak_edge_array[i_edge];
+          Array array         = weak_edge_array[i_edge];
           const size_t number = edge_number[i_edge];
           for (size_t i = 0; i < array.size(); ++i) {
             array[i] = number + (parallel::rank() + 1) * i;
@@ -607,7 +607,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
           bool is_synchronized = (parallel::size() > 1);
           bool is_valid        = true;
           for (EdgeId i_edge = 0; i_edge < mesh_3d.numberOfEdges(); ++i_edge) {
-            SubArray array      = edge_array[i_edge];
+            Array array         = edge_array[i_edge];
             const size_t number = edge_number[i_edge];
             const size_t owner  = edge_owner[i_edge];
             if (edge_is_owned[i_edge]) {
@@ -632,7 +632,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
           bool is_synchronized = true;
           for (EdgeId i_edge = 0; i_edge < mesh_3d.numberOfEdges(); ++i_edge) {
-            SubArray array      = edge_array[i_edge];
+            Array array         = edge_array[i_edge];
             const size_t number = edge_number[i_edge];
             const size_t owner  = edge_owner[i_edge];
             for (size_t i = 0; i < array.size(); ++i) {
@@ -650,7 +650,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
         auto face_number = connectivity.faceNumber();
 
         for (FaceId i_face = 0; i_face < mesh_3d.numberOfFaces(); ++i_face) {
-          SubArray array      = weak_face_array[i_face];
+          Array array         = weak_face_array[i_face];
           const size_t number = face_number[i_face];
           for (size_t i = 0; i < array.size(); ++i) {
             array[i] = number + (parallel::rank() + 1) * i;
@@ -668,7 +668,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
           bool is_synchronized = (parallel::size() > 1);
           bool is_valid        = true;
           for (FaceId i_face = 0; i_face < mesh_3d.numberOfFaces(); ++i_face) {
-            SubArray array      = face_array[i_face];
+            Array array         = face_array[i_face];
             const size_t number = face_number[i_face];
             const size_t owner  = face_owner[i_face];
             if (face_is_owned[i_face]) {
@@ -693,7 +693,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
           bool is_synchronized = true;
           for (FaceId i_face = 0; i_face < mesh_3d.numberOfFaces(); ++i_face) {
-            SubArray array      = face_array[i_face];
+            Array array         = face_array[i_face];
             const size_t number = face_number[i_face];
             const size_t owner  = face_owner[i_face];
             for (size_t i = 0; i < array.size(); ++i) {
@@ -711,7 +711,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
         auto cell_number = connectivity.cellNumber();
 
         for (CellId i_cell = 0; i_cell < mesh_3d.numberOfCells(); ++i_cell) {
-          SubArray array      = weak_cell_array[i_cell];
+          Array array         = weak_cell_array[i_cell];
           const size_t number = cell_number[i_cell];
           for (size_t i = 0; i < array.size(); ++i) {
             array[i] = number + (parallel::rank() + 1) * i;
@@ -729,7 +729,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
           bool is_synchronized = (parallel::size() > 1);
           bool is_valid        = true;
           for (CellId i_cell = 0; i_cell < mesh_3d.numberOfCells(); ++i_cell) {
-            SubArray array      = cell_array[i_cell];
+            Array array         = cell_array[i_cell];
             const size_t number = cell_number[i_cell];
             const size_t owner  = cell_owner[i_cell];
             if (cell_is_owned[i_cell]) {
@@ -754,7 +754,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
           bool is_synchronized = true;
           for (CellId i_cell = 0; i_cell < mesh_3d.numberOfCells(); ++i_cell) {
-            SubArray array      = cell_array[i_cell];
+            Array array         = cell_array[i_cell];
             const size_t number = cell_number[i_cell];
             const size_t owner  = cell_owner[i_cell];
             for (size_t i = 0; i < array.size(); ++i) {