diff --git a/src/mesh/ItemArray.hpp b/src/mesh/ItemArray.hpp
index a38f950954084a3527717a2a858d452e3e3cf2ad..78e5f8dd64c430ec8149c3e608d444897987f108 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[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/mesh/SubItemArrayPerItem.hpp b/src/mesh/SubItemArrayPerItem.hpp
index 431c3f4bee962857ff0b361c56a9ab5e987fde35..d39bff6a3015b61fbbe853dbeb9970e237068068 100644
--- a/src/mesh/SubItemArrayPerItem.hpp
+++ b/src/mesh/SubItemArrayPerItem.hpp
@@ -8,6 +8,8 @@
 #include <mesh/ItemType.hpp>
 #include <utils/Array.hpp>
 #include <utils/PugsAssert.hpp>
+#include <utils/Table.hpp>
+
 #include <utils/SubArray.hpp>
 
 #include <memory>
@@ -34,9 +36,7 @@ class SubItemArrayPerItem
   ConnectivityPtr m_connectivity_ptr;
 
   typename ConnectivityMatrix::HostRowType m_host_row_map;
-  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 SubItemArrayPerItem<std::add_const_t<DataType>, ItemOfItem, ConnectivitySharedPtr>;
@@ -53,56 +53,6 @@ class SubItemArrayPerItem
  public:
   using ToShared = SubItemArrayPerItem<DataType, ItemOfItem, ConnectivitySharedPtr>;
 
-  class SubView
-  {
-   public:
-    using data_type = DataType;
-
-   private:
-    PUGS_RESTRICT DataType* const m_sub_arrays;
-    const size_t m_size;
-    const size_t m_size_of_arrays;
-
-   public:
-    template <typename IndexType>
-    PUGS_FORCEINLINE SubArray<DataType>
-    operator[](IndexType i) const noexcept(NO_ASSERT)
-    {
-      static_assert(std::is_integral_v<IndexType>, "SubView is indexed by integral arrays");
-      Assert(static_cast<size_t>(i) < m_size);
-      return SubArray<DataType>(m_sub_arrays, i * m_size_of_arrays, m_size_of_arrays);
-    }
-
-    template <typename IndexType>
-    PUGS_FORCEINLINE SubArray<DataType>
-    operator[](IndexType i) noexcept(NO_ASSERT)
-    {
-      static_assert(std::is_integral_v<IndexType>, "SubView is indexed by integral arrays");
-      Assert(static_cast<size_t>(i) < m_size);
-      return SubArray<DataType>(m_sub_arrays, i * m_size_of_arrays, m_size_of_arrays);
-    }
-
-    PUGS_INLINE
-    size_t
-    size() const noexcept
-    {
-      return m_size;
-    }
-
-    SubView(const SubView&) = delete;
-
-    PUGS_INLINE
-    SubView(SubView&&) noexcept = delete;
-
-    PUGS_INLINE
-    SubView(const Array<DataType>& arrays, size_t begin, size_t end, size_t size_of_arrays) noexcept(NO_ASSERT)
-      : m_sub_arrays{&(arrays[begin * size_of_arrays])}, m_size{end - begin}, m_size_of_arrays{size_of_arrays}
-    {
-      Assert(begin <= end);
-      Assert(end <= arrays.size());
-    }
-  };
-
   friend PUGS_INLINE SubItemArrayPerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityPtr>
   copy(SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& source)
   {
@@ -110,8 +60,8 @@ class SubItemArrayPerItem
 
     image.m_connectivity_ptr = source.m_connectivity_ptr;
     image.m_host_row_map     = source.m_host_row_map;
-    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;
   }
 
@@ -134,34 +84,34 @@ class SubItemArrayPerItem
   }
 
   template <typename IndexType, typename SubIndexType>
-  PUGS_FORCEINLINE SubArray<DataType>
+  PUGS_FORCEINLINE Array<DataType>
   operator()(IndexType item_id, SubIndexType i) const noexcept(NO_ASSERT)
   {
     static_assert(std::is_same_v<IndexType, ItemId>, "first index must be of the correct ItemId type");
     static_assert(std::is_integral_v<SubIndexType>, "second index must be an integral type");
     Assert(this->isBuilt());
     Assert(i + m_host_row_map(size_t{item_id}) < m_host_row_map(size_t{item_id} + 1));
-    return SubArray(m_arrays_values, (m_host_row_map(size_t{item_id}) + i) * m_size_of_arrays, m_size_of_arrays);
+    return m_values[m_host_row_map(size_t{item_id}) + i];
   }
 
   // Following Kokkos logic, these classes are view and const view does allow
   // changes in data
   template <typename ArrayIndexType>
-  DataType&
+  Array<DataType>
   operator[](const ArrayIndexType& i) const noexcept(NO_ASSERT)
   {
     static_assert(std::is_integral_v<ArrayIndexType>, "index must be an integral type");
     Assert(this->isBuilt());
-    Assert(static_cast<size_t>(i) < m_arrays_values.size());
-    return m_arrays_values[i];
+    Assert(static_cast<size_t>(i) < m_values.nbRows());
+    return m_values[i];
   }
 
   PUGS_INLINE
   size_t
-  numberOfValues() const noexcept(NO_ASSERT)
+  numberOfArrays() const noexcept(NO_ASSERT)
   {
     Assert(this->isBuilt());
-    return m_arrays_values.size();
+    return m_values.nbRows();
   }
 
   PUGS_INLINE
@@ -177,7 +127,7 @@ class SubItemArrayPerItem
   sizeOfArrays() const noexcept(NO_ASSERT)
   {
     Assert(this->isBuilt());
-    return m_size_of_arrays;
+    return m_values.nbColumns();
   }
 
   template <typename IndexType>
@@ -190,30 +140,38 @@ class SubItemArrayPerItem
     return m_host_row_map(size_t{item_id} + 1) - m_host_row_map(size_t{item_id});
   }
 
+  PUGS_INLINE
+  void
+  fill(const DataType& data) const noexcept
+  {
+    static_assert(not std::is_const_v<DataType>, "Cannot modify ItemValue of const");
+    m_values.fill(data);
+  }
+
   template <typename IndexType>
-  PUGS_INLINE SubView
-  itemArrays(IndexType item_id) noexcept(NO_ASSERT)
+  PUGS_INLINE Table<DataType>
+  itemTable(IndexType item_id) noexcept(NO_ASSERT)
   {
     static_assert(std::is_same_v<IndexType, ItemId>, "index must be an ItemId");
     Assert(this->isBuilt());
     Assert(item_id < this->numberOfItems());
     const auto& item_begin = m_host_row_map(size_t{item_id});
     const auto& item_end   = m_host_row_map(size_t{item_id} + 1);
-    return SubView(m_arrays_values, item_begin, item_end, m_size_of_arrays);
+    return subTable(m_values, item_begin, item_end - item_begin, 0, this->sizeOfArrays());
   }
 
   // Following Kokkos logic, these classes are view and const view does allow
   // changes in data
   template <typename IndexType>
-  PUGS_INLINE SubView
-  itemArrays(IndexType item_id) const noexcept(NO_ASSERT)
+  PUGS_INLINE Table<DataType>
+  itemTable(IndexType item_id) const noexcept(NO_ASSERT)
   {
     static_assert(std::is_same_v<IndexType, ItemId>, "index must be an ItemId");
     Assert(this->isBuilt());
     Assert(item_id < this->numberOfItems());
     const auto& item_begin = m_host_row_map(size_t{item_id});
     const auto& item_end   = m_host_row_map(size_t{item_id} + 1);
-    return SubView(m_arrays_values, item_begin, item_end, m_size_of_arrays);
+    return subTable(m_values, item_begin, item_end - item_begin, 0, this->sizeOfArrays());
   }
 
   template <typename DataType2, typename ConnectivityPtr2>
@@ -227,9 +185,8 @@ class SubItemArrayPerItem
     static_assert(((std::is_const_v<DataType2> and std::is_const_v<DataType>) or not std::is_const_v<DataType2>),
                   "Cannot assign SubItemArrayPerItem of const data to "
                   "SubItemArrayPerItem of non-const data");
-    m_host_row_map   = sub_item_array_per_item.m_host_row_map;
-    m_arrays_values  = sub_item_array_per_item.m_arrays_values;
-    m_size_of_arrays = sub_item_array_per_item.m_size_of_arrays;
+    m_host_row_map = sub_item_array_per_item.m_host_row_map;
+    m_values       = sub_item_array_per_item.m_values;
 
     if constexpr (std::is_same_v<ConnectivityPtr, ConnectivitySharedPtr> and
                   std::is_same_v<ConnectivityPtr2, ConnectivityWeakPtr>) {
@@ -251,16 +208,16 @@ class SubItemArrayPerItem
 
   SubItemArrayPerItem() = default;
 
-  SubItemArrayPerItem(const IConnectivity& connectivity, size_t size_of_array) noexcept
-    : m_connectivity_ptr{connectivity.shared_ptr()}, m_size_of_arrays(size_of_array)
+  SubItemArrayPerItem(const IConnectivity& connectivity, size_t size_of_arrays) noexcept
+    : m_connectivity_ptr{connectivity.shared_ptr()}
   {
     static_assert(not std::is_const_v<DataType>, "Cannot allocate SubItemArrayPerItem of const data: only "
                                                  "view is supported");
 
     ConnectivityMatrix connectivity_matrix = connectivity._getMatrix(item_type, sub_item_type);
 
-    m_host_row_map  = connectivity_matrix.rowsMap();
-    m_arrays_values = Array<std::remove_const_t<DataType>>(connectivity_matrix.numEntries() * m_size_of_arrays);
+    m_host_row_map = connectivity_matrix.rowsMap();
+    m_values       = Table<std::remove_const_t<DataType>>(connectivity_matrix.numEntries(), size_of_arrays);
   }
 
   ~SubItemArrayPerItem() = default;
diff --git a/src/mesh/SubItemValuePerItem.hpp b/src/mesh/SubItemValuePerItem.hpp
index 63bff8a35848746535df9d05c7df0a12c899c3a6..b86b075ecbb44ff0f95a1a28b1f8f2468d98b2e3 100644
--- a/src/mesh/SubItemValuePerItem.hpp
+++ b/src/mesh/SubItemValuePerItem.hpp
@@ -50,55 +50,6 @@ class SubItemValuePerItem
  public:
   using ToShared = SubItemValuePerItem<DataType, ItemOfItem, ConnectivitySharedPtr>;
 
-  class SubView
-  {
-   public:
-    using data_type = DataType;
-
-   private:
-    PUGS_RESTRICT DataType* const m_sub_values;
-    const size_t m_size;
-
-   public:
-    template <typename IndexType>
-    PUGS_FORCEINLINE DataType&
-    operator[](IndexType i) const noexcept(NO_ASSERT)
-    {
-      static_assert(std::is_integral_v<IndexType>, "SubView is indexed by integral values");
-      Assert(i < m_size);
-      return m_sub_values[i];
-    }
-
-    template <typename IndexType>
-    PUGS_FORCEINLINE DataType&
-    operator[](IndexType i) noexcept(NO_ASSERT)
-    {
-      static_assert(std::is_integral_v<IndexType>, "SubView is indexed by integral values");
-      Assert(i < m_size);
-      return m_sub_values[i];
-    }
-
-    PUGS_INLINE
-    size_t
-    size() const noexcept
-    {
-      return m_size;
-    }
-
-    SubView(const SubView&)     = delete;
-    SubView(SubView&&) noexcept = delete;
-
-    PUGS_INLINE
-    SubView(const Array<DataType>& values, size_t begin, size_t end) noexcept(NO_ASSERT)
-      : m_sub_values(&(values[begin])), m_size(end - begin)
-    {
-      Assert(begin <= end);
-      Assert(end <= values.size());
-    }
-
-    ~SubView() = default;
-  };
-
   friend PUGS_INLINE SubItemValuePerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityPtr>
   copy(SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& source)
   {
@@ -178,8 +129,16 @@ class SubItemValuePerItem
     return m_host_row_map(size_t{item_id} + 1) - m_host_row_map(size_t{item_id});
   }
 
+  PUGS_INLINE
+  void
+  fill(const DataType& data) const noexcept
+  {
+    static_assert(not std::is_const_v<DataType>, "Cannot modify ItemValue of const");
+    m_values.fill(data);
+  }
+
   template <typename IndexType>
-  PUGS_INLINE SubView
+  PUGS_INLINE Array<DataType>
   itemValues(IndexType item_id) noexcept(NO_ASSERT)
   {
     static_assert(std::is_same_v<IndexType, ItemId>, "index must be an ItemId");
@@ -187,13 +146,13 @@ class SubItemValuePerItem
     Assert(item_id < this->numberOfItems());
     const auto& item_begin = m_host_row_map(size_t{item_id});
     const auto& item_end   = m_host_row_map(size_t{item_id} + 1);
-    return SubView(m_values, item_begin, item_end);
+    return subArray(m_values, item_begin, item_end - item_begin);
   }
 
   // Following Kokkos logic, these classes are view and const view does allow
   // changes in data
   template <typename IndexType>
-  PUGS_INLINE SubView
+  PUGS_INLINE Array<DataType>
   itemValues(IndexType item_id) const noexcept(NO_ASSERT)
   {
     static_assert(std::is_same_v<IndexType, ItemId>, "index must be an ItemId");
@@ -201,7 +160,7 @@ class SubItemValuePerItem
     Assert(item_id < this->numberOfItems());
     const auto& item_begin = m_host_row_map(size_t{item_id});
     const auto& item_end   = m_host_row_map(size_t{item_id} + 1);
-    return SubView(m_values, item_begin, item_end);
+    return subArray(m_values, item_begin, item_end - item_begin);
   }
 
   template <typename DataType2, typename ConnectivityPtr2>
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..67f544b997f32086131a549237de091e75ba60c1
--- /dev/null
+++ b/src/utils/Table.hpp
@@ -0,0 +1,144 @@
+#ifndef TABLE_HPP
+#define TABLE_HPP
+
+#include <Kokkos_CopyViews.hpp>
+#include <utils/Array.hpp>
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+#include <utils/PugsUtils.hpp>
+
+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);
+  }
+
+  PUGS_INLINE
+  Array<DataType> operator[](index_type i) const
+  {
+    Assert(i < this->nbRows());
+    return encapsulate(Kokkos::View<DataType*>(m_values, i, Kokkos::ALL));
+  }
+
+  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);
+
+  template <typename DataType2>
+  friend PUGS_INLINE Table<DataType2> subTable(const Table<DataType2>& table,
+                                               typename Table<DataType2>::index_type row_begin,
+                                               typename Table<DataType2>::index_type row_size,
+                                               typename Table<DataType2>::index_type column_begin,
+                                               typename Table<DataType2>::index_type column_size);
+
+  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;
+}
+
+template <typename DataType>
+PUGS_INLINE Table<DataType>
+subTable(const Table<DataType>& table,
+         typename Table<DataType>::index_type row_begin,
+         typename Table<DataType>::index_type row_size,
+         typename Table<DataType>::index_type column_begin,
+         typename Table<DataType>::index_type column_size)
+{
+  Assert(row_begin < table.nbRows());
+  Assert(row_begin + row_size <= table.nbRows());
+  Assert(column_begin < table.nbColumns());
+  Assert(column_begin + column_size <= table.nbColumns());
+  return encapsulate(Kokkos::View<DataType**>(table.m_values, std::make_pair(row_begin, row_begin + row_size),
+                                              std::make_pair(column_begin, column_begin + column_size)));
+}
+
+#endif   // TABLE_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index c99154ff0dee78cf7aadf95c70cccfbeb36ffa75..c36830918d3d769870250d212ea6d04e5d8316fc 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -85,8 +85,8 @@ add_executable (unit_tests
   test_PugsUtils.cpp
   test_RevisionInfo.cpp
   test_SparseMatrixDescriptor.cpp
-  test_SubArray.cpp
   test_SymbolTable.cpp
+  test_Table.cpp
   test_Timer.cpp
   test_TinyMatrix.cpp
   test_TinyVector.cpp
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) {
diff --git a/tests/test_SubItemArrayPerItem.cpp b/tests/test_SubItemArrayPerItem.cpp
index 5438398b0c1d5ba0559214473f41d13dba3cb8b4..7636ad420dc87ddf029d21b154ea0067a4eb5684 100644
--- a/tests/test_SubItemArrayPerItem.cpp
+++ b/tests/test_SubItemArrayPerItem.cpp
@@ -51,15 +51,13 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
   SECTION("dimensions")
   {
-    auto number_of_values = [](const auto& sub_item_array_per_item) -> size_t {
+    auto number_of_arrays = [](const auto& sub_item_array_per_item) -> size_t {
       using SubItemArrayPerItemType = std::decay_t<decltype(sub_item_array_per_item)>;
       using ItemId                  = typename SubItemArrayPerItemType::ItemId;
 
       size_t number = 0;
       for (ItemId item_id = 0; item_id < sub_item_array_per_item.numberOfItems(); ++item_id) {
-        for (size_t i = 0; i < sub_item_array_per_item.numberOfSubArrays(item_id); ++i) {
-          number += sub_item_array_per_item(item_id, i).size();
-        }
+        number += sub_item_array_per_item.numberOfSubArrays(item_id);
       }
       return number;
     };
@@ -73,7 +71,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
       {
         NodeArrayPerCell<int> node_array_per_cell{connectivity, 3};
         REQUIRE(node_array_per_cell.numberOfItems() == connectivity.numberOfCells());
-        REQUIRE(node_array_per_cell.numberOfValues() == number_of_values(node_array_per_cell));
+        REQUIRE(node_array_per_cell.numberOfArrays() == number_of_arrays(node_array_per_cell));
         REQUIRE(node_array_per_cell.sizeOfArrays() == 3);
 
         auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
@@ -82,7 +80,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
           for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
             is_correct &=
               (cell_to_node_matrix[cell_id].size() == node_array_per_cell.numberOfSubArrays(cell_id)) and
-              (node_array_per_cell.itemArrays(cell_id).size() == node_array_per_cell.numberOfSubArrays(cell_id));
+              (node_array_per_cell.itemTable(cell_id).nbRows() == node_array_per_cell.numberOfSubArrays(cell_id));
           }
           REQUIRE(is_correct);
         }
@@ -92,14 +90,14 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
           bool is_correct = true;
           for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
             is_correct &=
-              (const_node_array_per_cell.itemArrays(cell_id).size() == node_array_per_cell.numberOfSubArrays(cell_id));
+              (const_node_array_per_cell.itemTable(cell_id).nbRows() == node_array_per_cell.numberOfSubArrays(cell_id));
           }
           REQUIRE(is_correct);
         }
 
         EdgeArrayPerCell<int> edge_array_per_cell{connectivity, 4};
         REQUIRE(edge_array_per_cell.numberOfItems() == connectivity.numberOfCells());
-        REQUIRE(edge_array_per_cell.numberOfValues() == number_of_values(edge_array_per_cell));
+        REQUIRE(edge_array_per_cell.numberOfArrays() == number_of_arrays(edge_array_per_cell));
         REQUIRE(edge_array_per_cell.sizeOfArrays() == 4);
 
         auto cell_to_edge_matrix = connectivity.cellToEdgeMatrix();
@@ -113,7 +111,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
         FaceArrayPerCell<int> face_array_per_cell{connectivity, 2};
         REQUIRE(face_array_per_cell.numberOfItems() == connectivity.numberOfCells());
-        REQUIRE(face_array_per_cell.numberOfValues() == number_of_values(face_array_per_cell));
+        REQUIRE(face_array_per_cell.numberOfArrays() == number_of_arrays(face_array_per_cell));
         REQUIRE(face_array_per_cell.sizeOfArrays() == 2);
 
         auto cell_to_face_matrix = connectivity.cellToFaceMatrix();
@@ -130,7 +128,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
       {
         CellArrayPerFace<int> cell_array_per_face{connectivity, 2};
         REQUIRE(cell_array_per_face.numberOfItems() == connectivity.numberOfFaces());
-        REQUIRE(cell_array_per_face.numberOfValues() == number_of_values(cell_array_per_face));
+        REQUIRE(cell_array_per_face.numberOfArrays() == number_of_arrays(cell_array_per_face));
         REQUIRE(cell_array_per_face.sizeOfArrays() == 2);
 
         auto face_to_cell_matrix = connectivity.faceToCellMatrix();
@@ -147,7 +145,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
       {
         CellArrayPerEdge<int> cell_array_per_edge{connectivity, 3};
         REQUIRE(cell_array_per_edge.numberOfItems() == connectivity.numberOfEdges());
-        REQUIRE(cell_array_per_edge.numberOfValues() == number_of_values(cell_array_per_edge));
+        REQUIRE(cell_array_per_edge.numberOfArrays() == number_of_arrays(cell_array_per_edge));
         REQUIRE(cell_array_per_edge.sizeOfArrays() == 3);
 
         auto edge_to_cell_matrix = connectivity.edgeToCellMatrix();
@@ -164,7 +162,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
       {
         CellArrayPerNode<int> cell_array_per_node{connectivity, 4};
         REQUIRE(cell_array_per_node.numberOfItems() == connectivity.numberOfNodes());
-        REQUIRE(cell_array_per_node.numberOfValues() == number_of_values(cell_array_per_node));
+        REQUIRE(cell_array_per_node.numberOfArrays() == number_of_arrays(cell_array_per_node));
         REQUIRE(cell_array_per_node.sizeOfArrays() == 4);
 
         auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
@@ -187,7 +185,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
       {
         NodeArrayPerCell<int> node_array_per_cell{connectivity, 5};
         REQUIRE(node_array_per_cell.numberOfItems() == connectivity.numberOfCells());
-        REQUIRE(node_array_per_cell.numberOfValues() == number_of_values(node_array_per_cell));
+        REQUIRE(node_array_per_cell.numberOfArrays() == number_of_arrays(node_array_per_cell));
         REQUIRE(node_array_per_cell.sizeOfArrays() == 5);
 
         auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
@@ -201,7 +199,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
         EdgeArrayPerCell<int> edge_array_per_cell{connectivity, 4};
         REQUIRE(edge_array_per_cell.numberOfItems() == connectivity.numberOfCells());
-        REQUIRE(edge_array_per_cell.numberOfValues() == number_of_values(edge_array_per_cell));
+        REQUIRE(edge_array_per_cell.numberOfArrays() == number_of_arrays(edge_array_per_cell));
         REQUIRE(edge_array_per_cell.sizeOfArrays() == 4);
 
         auto cell_to_edge_matrix = connectivity.cellToEdgeMatrix();
@@ -215,7 +213,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
         FaceArrayPerCell<int> face_array_per_cell{connectivity, 3};
         REQUIRE(face_array_per_cell.numberOfItems() == connectivity.numberOfCells());
-        REQUIRE(face_array_per_cell.numberOfValues() == number_of_values(face_array_per_cell));
+        REQUIRE(face_array_per_cell.numberOfArrays() == number_of_arrays(face_array_per_cell));
         REQUIRE(face_array_per_cell.sizeOfArrays() == 3);
 
         auto cell_to_face_matrix = connectivity.cellToFaceMatrix();
@@ -232,7 +230,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
       {
         CellArrayPerFace<int> cell_array_per_face{connectivity, 3};
         REQUIRE(cell_array_per_face.numberOfItems() == connectivity.numberOfFaces());
-        REQUIRE(cell_array_per_face.numberOfValues() == number_of_values(cell_array_per_face));
+        REQUIRE(cell_array_per_face.numberOfArrays() == number_of_arrays(cell_array_per_face));
         REQUIRE(cell_array_per_face.sizeOfArrays() == 3);
 
         auto face_to_cell_matrix = connectivity.faceToCellMatrix();
@@ -246,7 +244,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
         NodeArrayPerFace<int> node_array_per_face{connectivity, 2};
         REQUIRE(node_array_per_face.numberOfItems() == connectivity.numberOfFaces());
-        REQUIRE(node_array_per_face.numberOfValues() == number_of_values(node_array_per_face));
+        REQUIRE(node_array_per_face.numberOfArrays() == number_of_arrays(node_array_per_face));
         REQUIRE(node_array_per_face.sizeOfArrays() == 2);
 
         auto face_to_node_matrix = connectivity.faceToNodeMatrix();
@@ -263,7 +261,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
       {
         CellArrayPerEdge<int> cell_array_per_edge{connectivity, 3};
         REQUIRE(cell_array_per_edge.numberOfItems() == connectivity.numberOfEdges());
-        REQUIRE(cell_array_per_edge.numberOfValues() == number_of_values(cell_array_per_edge));
+        REQUIRE(cell_array_per_edge.numberOfArrays() == number_of_arrays(cell_array_per_edge));
         REQUIRE(cell_array_per_edge.sizeOfArrays() == 3);
 
         auto edge_to_cell_matrix = connectivity.edgeToCellMatrix();
@@ -277,7 +275,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
         NodeArrayPerEdge<int> node_array_per_edge{connectivity, 5};
         REQUIRE(node_array_per_edge.numberOfItems() == connectivity.numberOfEdges());
-        REQUIRE(node_array_per_edge.numberOfValues() == number_of_values(node_array_per_edge));
+        REQUIRE(node_array_per_edge.numberOfArrays() == number_of_arrays(node_array_per_edge));
         REQUIRE(node_array_per_edge.sizeOfArrays() == 5);
 
         auto edge_to_node_matrix = connectivity.edgeToNodeMatrix();
@@ -294,7 +292,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
       {
         EdgeArrayPerNode<int> edge_array_per_node{connectivity, 4};
         REQUIRE(edge_array_per_node.numberOfItems() == connectivity.numberOfNodes());
-        REQUIRE(edge_array_per_node.numberOfValues() == number_of_values(edge_array_per_node));
+        REQUIRE(edge_array_per_node.numberOfArrays() == number_of_arrays(edge_array_per_node));
         REQUIRE(edge_array_per_node.sizeOfArrays() == 4);
 
         auto node_to_edge_matrix = connectivity.nodeToEdgeMatrix();
@@ -308,7 +306,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
         FaceArrayPerNode<int> face_array_per_node{connectivity, 3};
         REQUIRE(face_array_per_node.numberOfItems() == connectivity.numberOfNodes());
-        REQUIRE(face_array_per_node.numberOfValues() == number_of_values(face_array_per_node));
+        REQUIRE(face_array_per_node.numberOfArrays() == number_of_arrays(face_array_per_node));
         REQUIRE(face_array_per_node.sizeOfArrays() == 3);
 
         auto node_to_face_matrix = connectivity.nodeToFaceMatrix();
@@ -322,7 +320,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
         CellArrayPerNode<int> cell_array_per_node{connectivity, 2};
         REQUIRE(cell_array_per_node.numberOfItems() == connectivity.numberOfNodes());
-        REQUIRE(cell_array_per_node.numberOfValues() == number_of_values(cell_array_per_node));
+        REQUIRE(cell_array_per_node.numberOfArrays() == number_of_arrays(cell_array_per_node));
         REQUIRE(cell_array_per_node.sizeOfArrays() == 2);
 
         auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
@@ -344,7 +342,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
       {
         NodeArrayPerCell<int> node_array_per_cell{connectivity, 3};
         REQUIRE(node_array_per_cell.numberOfItems() == connectivity.numberOfCells());
-        REQUIRE(node_array_per_cell.numberOfValues() == number_of_values(node_array_per_cell));
+        REQUIRE(node_array_per_cell.numberOfArrays() == number_of_arrays(node_array_per_cell));
         REQUIRE(node_array_per_cell.sizeOfArrays() == 3);
 
         auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
@@ -358,7 +356,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
         EdgeArrayPerCell<int> edge_array_per_cell{connectivity, 4};
         REQUIRE(edge_array_per_cell.numberOfItems() == connectivity.numberOfCells());
-        REQUIRE(edge_array_per_cell.numberOfValues() == number_of_values(edge_array_per_cell));
+        REQUIRE(edge_array_per_cell.numberOfArrays() == number_of_arrays(edge_array_per_cell));
         REQUIRE(edge_array_per_cell.sizeOfArrays() == 4);
 
         auto cell_to_edge_matrix = connectivity.cellToEdgeMatrix();
@@ -372,7 +370,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
         FaceArrayPerCell<int> face_array_per_cell{connectivity, 3};
         REQUIRE(face_array_per_cell.numberOfItems() == connectivity.numberOfCells());
-        REQUIRE(face_array_per_cell.numberOfValues() == number_of_values(face_array_per_cell));
+        REQUIRE(face_array_per_cell.numberOfArrays() == number_of_arrays(face_array_per_cell));
         REQUIRE(face_array_per_cell.sizeOfArrays() == 3);
 
         auto cell_to_face_matrix = connectivity.cellToFaceMatrix();
@@ -389,7 +387,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
       {
         CellArrayPerFace<int> cell_array_per_face{connectivity, 5};
         REQUIRE(cell_array_per_face.numberOfItems() == connectivity.numberOfFaces());
-        REQUIRE(cell_array_per_face.numberOfValues() == number_of_values(cell_array_per_face));
+        REQUIRE(cell_array_per_face.numberOfArrays() == number_of_arrays(cell_array_per_face));
         REQUIRE(cell_array_per_face.sizeOfArrays() == 5);
 
         auto face_to_cell_matrix = connectivity.faceToCellMatrix();
@@ -403,7 +401,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
         EdgeArrayPerFace<int> edge_array_per_face{connectivity, 3};
         REQUIRE(edge_array_per_face.numberOfItems() == connectivity.numberOfFaces());
-        REQUIRE(edge_array_per_face.numberOfValues() == number_of_values(edge_array_per_face));
+        REQUIRE(edge_array_per_face.numberOfArrays() == number_of_arrays(edge_array_per_face));
         REQUIRE(edge_array_per_face.sizeOfArrays() == 3);
 
         auto face_to_edge_matrix = connectivity.faceToEdgeMatrix();
@@ -417,7 +415,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
         NodeArrayPerFace<int> node_array_per_face{connectivity, 2};
         REQUIRE(node_array_per_face.numberOfItems() == connectivity.numberOfFaces());
-        REQUIRE(node_array_per_face.numberOfValues() == number_of_values(node_array_per_face));
+        REQUIRE(node_array_per_face.numberOfArrays() == number_of_arrays(node_array_per_face));
         REQUIRE(node_array_per_face.sizeOfArrays() == 2);
 
         auto face_to_node_matrix = connectivity.faceToNodeMatrix();
@@ -434,7 +432,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
       {
         CellArrayPerEdge<int> cell_array_per_edge{connectivity, 3};
         REQUIRE(cell_array_per_edge.numberOfItems() == connectivity.numberOfEdges());
-        REQUIRE(cell_array_per_edge.numberOfValues() == number_of_values(cell_array_per_edge));
+        REQUIRE(cell_array_per_edge.numberOfArrays() == number_of_arrays(cell_array_per_edge));
         REQUIRE(cell_array_per_edge.sizeOfArrays() == 3);
 
         auto edge_to_cell_matrix = connectivity.edgeToCellMatrix();
@@ -448,7 +446,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
         FaceArrayPerEdge<int> face_array_per_edge{connectivity, 5};
         REQUIRE(face_array_per_edge.numberOfItems() == connectivity.numberOfEdges());
-        REQUIRE(face_array_per_edge.numberOfValues() == number_of_values(face_array_per_edge));
+        REQUIRE(face_array_per_edge.numberOfArrays() == number_of_arrays(face_array_per_edge));
         REQUIRE(face_array_per_edge.sizeOfArrays() == 5);
 
         auto edge_to_face_matrix = connectivity.edgeToFaceMatrix();
@@ -462,7 +460,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
         NodeArrayPerEdge<int> node_array_per_edge{connectivity, 3};
         REQUIRE(node_array_per_edge.numberOfItems() == connectivity.numberOfEdges());
-        REQUIRE(node_array_per_edge.numberOfValues() == number_of_values(node_array_per_edge));
+        REQUIRE(node_array_per_edge.numberOfArrays() == number_of_arrays(node_array_per_edge));
         REQUIRE(node_array_per_edge.sizeOfArrays() == 3);
 
         auto edge_to_node_matrix = connectivity.edgeToNodeMatrix();
@@ -479,7 +477,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
       {
         EdgeArrayPerNode<int> edge_array_per_node{connectivity, 3};
         REQUIRE(edge_array_per_node.numberOfItems() == connectivity.numberOfNodes());
-        REQUIRE(edge_array_per_node.numberOfValues() == number_of_values(edge_array_per_node));
+        REQUIRE(edge_array_per_node.numberOfArrays() == number_of_arrays(edge_array_per_node));
         REQUIRE(edge_array_per_node.sizeOfArrays() == 3);
 
         auto node_to_edge_matrix = connectivity.nodeToEdgeMatrix();
@@ -493,7 +491,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
         FaceArrayPerNode<int> face_array_per_node{connectivity, 4};
         REQUIRE(face_array_per_node.numberOfItems() == connectivity.numberOfNodes());
-        REQUIRE(face_array_per_node.numberOfValues() == number_of_values(face_array_per_node));
+        REQUIRE(face_array_per_node.numberOfArrays() == number_of_arrays(face_array_per_node));
         REQUIRE(face_array_per_node.sizeOfArrays() == 4);
 
         auto node_to_face_matrix = connectivity.nodeToFaceMatrix();
@@ -507,7 +505,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
         CellArrayPerNode<int> cell_array_per_node{connectivity, 5};
         REQUIRE(cell_array_per_node.numberOfItems() == connectivity.numberOfNodes());
-        REQUIRE(cell_array_per_node.numberOfValues() == number_of_values(cell_array_per_node));
+        REQUIRE(cell_array_per_node.numberOfArrays() == number_of_arrays(cell_array_per_node));
         REQUIRE(cell_array_per_node.sizeOfArrays() == 5);
 
         auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
@@ -542,14 +540,22 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
       }
       {
         bool is_same = true;
-        for (size_t i = 0; i < edge_arrays_per_cell.numberOfValues(); ++i) {
-          is_same &= (edge_arrays_per_cell[i] == i);
+        size_t k     = 0;
+        for (size_t i = 0; i < edge_arrays_per_cell.numberOfArrays(); ++i) {
+          for (size_t j = 0; j < edge_arrays_per_cell.sizeOfArrays(); ++j, ++k) {
+            is_same &= (edge_arrays_per_cell[i][j] == k);
+          }
         }
         REQUIRE(is_same);
       }
 
-      for (size_t i = 0; i < edge_arrays_per_cell.numberOfValues(); ++i) {
-        edge_arrays_per_cell[i] = i * i + 1;
+      {
+        size_t k = 0;
+        for (size_t i = 0; i < edge_arrays_per_cell.numberOfArrays(); ++i) {
+          for (size_t j = 0; j < edge_arrays_per_cell.sizeOfArrays(); ++j, ++k) {
+            edge_arrays_per_cell[i][j] = k * k + 1;
+          }
+        }
       }
       {
         bool is_same = true;
@@ -563,6 +569,22 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
         }
         REQUIRE(is_same);
       }
+
+      const EdgeArrayPerCell<const size_t> const_edge_arrays_per_cell = edge_arrays_per_cell;
+      {
+        bool is_same = true;
+        size_t i     = 0;
+        for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
+          const auto& cell_table = const_edge_arrays_per_cell.itemTable(cell_id);
+          for (size_t i_edge = 0; i_edge < cell_table.nbRows(); ++i_edge) {
+            const auto& array = cell_table[i_edge];
+            for (size_t l = 0; l < array.size(); ++l, ++i) {
+              is_same &= (array[l] == i * i + 1);
+            }
+          }
+        }
+        REQUIRE(is_same);
+      }
     }
 
     SECTION("2D")
@@ -583,13 +605,21 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
       }
       {
         bool is_same = true;
-        for (size_t i = 0; i < cell_arrays_per_face.numberOfValues(); ++i) {
-          is_same &= (cell_arrays_per_face[i] == i);
+        size_t k     = 0;
+        for (size_t i = 0; i < cell_arrays_per_face.numberOfArrays(); ++i) {
+          for (size_t j = 0; j < cell_arrays_per_face.sizeOfArrays(); ++j, ++k) {
+            is_same &= (cell_arrays_per_face[i][j] == k);
+          }
         }
         REQUIRE(is_same);
       }
-      for (size_t i = 0; i < cell_arrays_per_face.numberOfValues(); ++i) {
-        cell_arrays_per_face[i] = 3 * i + 1;
+      {
+        size_t k = 0;
+        for (size_t i = 0; i < cell_arrays_per_face.numberOfArrays(); ++i) {
+          for (size_t j = 0; j < cell_arrays_per_face.sizeOfArrays(); ++j, ++k) {
+            cell_arrays_per_face[i][j] = 3 * k + 1;
+          }
+        }
       }
       {
         bool is_same = true;
@@ -603,6 +633,22 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
         }
         REQUIRE(is_same);
       }
+
+      const CellArrayPerFace<const size_t> const_cell_arrays_per_face = cell_arrays_per_face;
+      {
+        bool is_same = true;
+        size_t i     = 0;
+        for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+          const auto& face_table = const_cell_arrays_per_face.itemTable(face_id);
+          for (size_t i_cell = 0; i_cell < face_table.nbRows(); ++i_cell) {
+            const auto& array = face_table[i_cell];
+            for (size_t l = 0; l < array.size(); ++l, ++i) {
+              is_same &= (array[l] == 3 * i + 1);
+            }
+          }
+        }
+        REQUIRE(is_same);
+      }
     }
 
     SECTION("3D")
@@ -622,14 +668,21 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
       }
       {
         bool is_same = true;
-        for (size_t i = 0; i < face_arrays_per_node.numberOfValues(); ++i) {
-          is_same &= (face_arrays_per_node[i] == i);
+        size_t k     = 0;
+        for (size_t i = 0; i < face_arrays_per_node.numberOfArrays(); ++i) {
+          for (size_t j = 0; j < face_arrays_per_node.sizeOfArrays(); ++j, ++k) {
+            is_same &= (face_arrays_per_node[i][j] == k);
+          }
+          REQUIRE(is_same);
         }
-        REQUIRE(is_same);
       }
-
-      for (size_t i = 0; i < face_arrays_per_node.numberOfValues(); ++i) {
-        face_arrays_per_node[i] = 3 + i * i;
+      {
+        size_t k = 0;
+        for (size_t i = 0; i < face_arrays_per_node.numberOfArrays(); ++i) {
+          for (size_t j = 0; j < face_arrays_per_node.sizeOfArrays(); ++j, ++k) {
+            face_arrays_per_node[i][j] = 3 + k * k;
+          }
+        }
       }
       {
         bool is_same = true;
@@ -643,6 +696,22 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
         }
         REQUIRE(is_same);
       }
+
+      const FaceArrayPerNode<const size_t> const_face_arrays_per_node = face_arrays_per_node;
+      {
+        bool is_same = true;
+        size_t i     = 0;
+        for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+          const auto& node_table = const_face_arrays_per_node.itemTable(node_id);
+          for (size_t i_face = 0; i_face < node_table.nbRows(); ++i_face) {
+            const auto& array = node_table[i_face];
+            for (size_t l = 0; l < array.size(); ++l, ++i) {
+              is_same &= (array[l] == 3 + i * i);
+            }
+          }
+        }
+        REQUIRE(is_same);
+      }
     }
   }
 
@@ -654,15 +723,21 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
     SECTION("classic")
     {
       NodeArrayPerCell<size_t> node_array_per_cell{connectivity, 3};
-      for (size_t i = 0; i < node_array_per_cell.numberOfValues(); ++i) {
-        node_array_per_cell[i] = i;
+      {
+        size_t k = 0;
+        for (size_t i = 0; i < node_array_per_cell.numberOfArrays(); ++i) {
+          for (size_t j = 0; j < node_array_per_cell.sizeOfArrays(); ++j, ++k) {
+            node_array_per_cell[i][j] = k;
+          }
+        }
       }
-
       NodeArrayPerCell<size_t> copy_node_array_per_cell = copy(node_array_per_cell);
       {
         bool is_same = true;
-        for (size_t i = 0; i < copy_node_array_per_cell.numberOfValues(); ++i) {
-          is_same &= (copy_node_array_per_cell[i] == node_array_per_cell[i]);
+        for (size_t i = 0; i < copy_node_array_per_cell.numberOfArrays(); ++i) {
+          for (size_t j = 0; j < node_array_per_cell.sizeOfArrays(); ++j) {
+            is_same &= (copy_node_array_per_cell[i][j] == node_array_per_cell[i][j]);
+          }
         }
 
         REQUIRE(is_same);
@@ -670,9 +745,9 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
       {
         for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
-          auto node_array_list = node_array_per_cell.itemArrays(cell_id);
+          auto cell_table = node_array_per_cell.itemTable(cell_id);
           for (size_t i_node = 0; i_node < node_array_per_cell.numberOfSubArrays(cell_id); ++i_node) {
-            auto node_array = node_array_list[i_node];
+            auto node_array = cell_table[i_node];
             for (size_t i = 0; i < node_array.size(); ++i) {
               node_array[i] = cell_id + i_node + i;
             }
@@ -682,8 +757,10 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
       {
         bool is_same = true;
-        for (size_t i = 0; i < copy_node_array_per_cell.numberOfValues(); ++i) {
-          is_same &= (copy_node_array_per_cell[i] == node_array_per_cell[i]);
+        for (size_t i = 0; i < copy_node_array_per_cell.numberOfArrays(); ++i) {
+          for (size_t j = 0; j < copy_node_array_per_cell.sizeOfArrays(); ++j) {
+            is_same &= (copy_node_array_per_cell[i][j] == node_array_per_cell[i][j]);
+          }
         }
 
         REQUIRE(not is_same);
@@ -693,15 +770,22 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
     SECTION("from weak")
     {
       WeakNodeArrayPerCell<size_t> node_array_per_cell{connectivity, 3};
-      for (size_t i = 0; i < node_array_per_cell.numberOfValues(); ++i) {
-        node_array_per_cell[i] = i;
+      {
+        size_t k = 0;
+        for (size_t i = 0; i < node_array_per_cell.numberOfArrays(); ++i) {
+          for (size_t j = 0; j < node_array_per_cell.sizeOfArrays(); ++j, ++k) {
+            node_array_per_cell[i][j] = k;
+          }
+        }
       }
 
       NodeArrayPerCell<size_t> copy_node_array_per_cell = copy(node_array_per_cell);
       {
         bool is_same = true;
-        for (size_t i = 0; i < copy_node_array_per_cell.numberOfValues(); ++i) {
-          is_same &= (copy_node_array_per_cell[i] == node_array_per_cell[i]);
+        for (size_t i = 0; i < copy_node_array_per_cell.numberOfArrays(); ++i) {
+          for (size_t j = 0; j < node_array_per_cell.sizeOfArrays(); ++j) {
+            is_same &= (copy_node_array_per_cell[i][j] == node_array_per_cell[i][j]);
+          }
         }
 
         REQUIRE(is_same);
@@ -709,9 +793,9 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
       {
         for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
-          auto node_array_list = node_array_per_cell.itemArrays(cell_id);
+          auto cell_table = node_array_per_cell.itemTable(cell_id);
           for (size_t i_node = 0; i_node < node_array_per_cell.numberOfSubArrays(cell_id); ++i_node) {
-            auto node_array = node_array_list[i_node];
+            auto node_array = cell_table[i_node];
             for (size_t i = 0; i < node_array.size(); ++i) {
               node_array[i] = cell_id + i_node + i;
             }
@@ -721,8 +805,10 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
 
       {
         bool is_same = true;
-        for (size_t i = 0; i < copy_node_array_per_cell.numberOfValues(); ++i) {
-          is_same &= (copy_node_array_per_cell[i] == node_array_per_cell[i]);
+        for (size_t i = 0; i < copy_node_array_per_cell.numberOfArrays(); ++i) {
+          for (size_t j = 0; j < node_array_per_cell.sizeOfArrays(); ++j) {
+            is_same &= (copy_node_array_per_cell[i][j] == node_array_per_cell[i][j]);
+          }
         }
 
         REQUIRE(not is_same);
@@ -736,9 +822,13 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
     const Connectivity<2>& connectivity  = mesh_2d.connectivity();
 
     WeakFaceArrayPerCell<int> weak_face_array_per_cell{connectivity, 3};
-
-    for (size_t i = 0; i < weak_face_array_per_cell.numberOfValues(); ++i) {
-      weak_face_array_per_cell[i] = i;
+    {
+      size_t k = 0;
+      for (size_t i = 0; i < weak_face_array_per_cell.numberOfArrays(); ++i) {
+        for (size_t j = 0; j < weak_face_array_per_cell.sizeOfArrays(); ++j, ++k) {
+          weak_face_array_per_cell[i][j] = k;
+        }
+      }
     }
 
     FaceArrayPerCell<const int> face_array_per_cell{weak_face_array_per_cell};
@@ -747,8 +837,10 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
     REQUIRE(face_array_per_cell.sizeOfArrays() == weak_face_array_per_cell.sizeOfArrays());
 
     bool is_same = true;
-    for (size_t i = 0; i < weak_face_array_per_cell.numberOfValues(); ++i) {
-      is_same &= (face_array_per_cell[i] == weak_face_array_per_cell[i]);
+    for (size_t i = 0; i < weak_face_array_per_cell.numberOfArrays(); ++i) {
+      for (size_t j = 0; j < weak_face_array_per_cell.sizeOfArrays(); ++j) {
+        is_same &= (face_array_per_cell[i][j] == weak_face_array_per_cell[i][j]);
+      }
     }
     REQUIRE(is_same);
   }
@@ -760,37 +852,37 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
     {
       CellArrayPerNode<int> cell_array_per_node;
       REQUIRE_THROWS_AS(cell_array_per_node[0], AssertError);
-      REQUIRE_THROWS_AS(cell_array_per_node.itemArrays(NodeId{0}), AssertError);
+      REQUIRE_THROWS_AS(cell_array_per_node.itemTable(NodeId{0}), AssertError);
       REQUIRE_THROWS_AS(cell_array_per_node(NodeId{0}, 0), AssertError);
       REQUIRE_THROWS_AS(cell_array_per_node.sizeOfArrays(), AssertError);
-      REQUIRE_THROWS_AS(cell_array_per_node.numberOfValues(), AssertError);
+      REQUIRE_THROWS_AS(cell_array_per_node.numberOfArrays(), AssertError);
       REQUIRE_THROWS_AS(cell_array_per_node.numberOfItems(), AssertError);
       REQUIRE_THROWS_AS(cell_array_per_node.numberOfSubArrays(NodeId{0}), AssertError);
 
       FaceArrayPerCell<int> face_array_per_cell;
       REQUIRE_THROWS_AS(face_array_per_cell[0], AssertError);
-      REQUIRE_THROWS_AS(face_array_per_cell.itemArrays(CellId{0}), AssertError);
+      REQUIRE_THROWS_AS(face_array_per_cell.itemTable(CellId{0}), AssertError);
       REQUIRE_THROWS_AS(face_array_per_cell(CellId{0}, 0), AssertError);
       REQUIRE_THROWS_AS(face_array_per_cell.sizeOfArrays(), AssertError);
-      REQUIRE_THROWS_AS(face_array_per_cell.numberOfValues(), AssertError);
+      REQUIRE_THROWS_AS(face_array_per_cell.numberOfArrays(), AssertError);
       REQUIRE_THROWS_AS(face_array_per_cell.numberOfItems(), AssertError);
       REQUIRE_THROWS_AS(face_array_per_cell.numberOfSubArrays(CellId{0}), AssertError);
 
       CellArrayPerEdge<int> cell_array_per_edge;
       REQUIRE_THROWS_AS(cell_array_per_edge[0], AssertError);
-      REQUIRE_THROWS_AS(cell_array_per_edge.itemArrays(EdgeId{0}), AssertError);
+      REQUIRE_THROWS_AS(cell_array_per_edge.itemTable(EdgeId{0}), AssertError);
       REQUIRE_THROWS_AS(cell_array_per_edge(EdgeId{0}, 0), AssertError);
       REQUIRE_THROWS_AS(cell_array_per_edge.sizeOfArrays(), AssertError);
-      REQUIRE_THROWS_AS(cell_array_per_edge.numberOfValues(), AssertError);
+      REQUIRE_THROWS_AS(cell_array_per_edge.numberOfArrays(), AssertError);
       REQUIRE_THROWS_AS(cell_array_per_edge.numberOfItems(), AssertError);
       REQUIRE_THROWS_AS(cell_array_per_edge.numberOfSubArrays(EdgeId{0}), AssertError);
 
       NodeArrayPerFace<int> node_array_per_face;
       REQUIRE_THROWS_AS(node_array_per_face[0], AssertError);
-      REQUIRE_THROWS_AS(node_array_per_face.itemArrays(FaceId{0}), AssertError);
+      REQUIRE_THROWS_AS(node_array_per_face.itemTable(FaceId{0}), AssertError);
       REQUIRE_THROWS_AS(node_array_per_face(FaceId{0}, 0), AssertError);
       REQUIRE_THROWS_AS(node_array_per_face.sizeOfArrays(), AssertError);
-      REQUIRE_THROWS_AS(node_array_per_face.numberOfValues(), AssertError);
+      REQUIRE_THROWS_AS(node_array_per_face.numberOfArrays(), AssertError);
       REQUIRE_THROWS_AS(node_array_per_face.numberOfItems(), AssertError);
       REQUIRE_THROWS_AS(node_array_per_face.numberOfSubArrays(FaceId{0}), AssertError);
     }
@@ -806,13 +898,13 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
         REQUIRE_THROWS_AS(cell_array_per_face(invalid_face_id, 0), AssertError);
       }
       if (connectivity.numberOfFaces() > 0) {
-        FaceId face_id          = 0;
-        const auto& cell_arrays = cell_array_per_face.itemArrays(face_id);
-        REQUIRE_THROWS_AS(cell_array_per_face(face_id, cell_arrays.size()), AssertError);
-        REQUIRE_THROWS_AS(cell_arrays[cell_arrays.size()], AssertError);
-        REQUIRE_THROWS_AS(cell_array_per_face.itemArrays(face_id)[cell_arrays.size()], AssertError);
-        REQUIRE_THROWS_AS(cell_array_per_face.itemArrays(face_id)[0][cell_array_per_face.sizeOfArrays()], AssertError);
-        REQUIRE_THROWS_AS(cell_array_per_face.itemArrays(face_id)[0][cell_array_per_face.sizeOfArrays()] = 2,
+        FaceId face_id         = 0;
+        const auto& face_table = cell_array_per_face.itemTable(face_id);
+        REQUIRE_THROWS_AS(cell_array_per_face(face_id, face_table.nbRows()), AssertError);
+        REQUIRE_THROWS_AS(face_table[face_table.nbRows()], AssertError);
+        REQUIRE_THROWS_AS(cell_array_per_face.itemTable(face_id)[face_table.nbRows()], AssertError);
+        REQUIRE_THROWS_AS(cell_array_per_face.itemTable(face_id)[0][cell_array_per_face.sizeOfArrays()], AssertError);
+        REQUIRE_THROWS_AS(cell_array_per_face.itemTable(face_id)[0][cell_array_per_face.sizeOfArrays()] = 2,
                           AssertError);
       }
 
@@ -822,13 +914,13 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
         REQUIRE_THROWS_AS(face_array_per_node(invalid_node_id, 0), AssertError);
       }
       if (connectivity.numberOfNodes() > 0) {
-        NodeId node_id          = 0;
-        const auto& face_arrays = face_array_per_node.itemArrays(node_id);
-        REQUIRE_THROWS_AS(face_array_per_node(node_id, face_arrays.size()), AssertError);
-        REQUIRE_THROWS_AS(face_arrays[face_arrays.size()], AssertError);
-        REQUIRE_THROWS_AS(face_array_per_node.itemArrays(node_id)[face_arrays.size()], AssertError);
-        REQUIRE_THROWS_AS(face_array_per_node.itemArrays(node_id)[0][face_array_per_node.sizeOfArrays()], AssertError);
-        REQUIRE_THROWS_AS(face_array_per_node.itemArrays(node_id)[0][face_array_per_node.sizeOfArrays()] = 2,
+        NodeId node_id         = 0;
+        const auto& node_table = face_array_per_node.itemTable(node_id);
+        REQUIRE_THROWS_AS(face_array_per_node(node_id, node_table.nbRows()), AssertError);
+        REQUIRE_THROWS_AS(node_table[node_table.nbRows()], AssertError);
+        REQUIRE_THROWS_AS(face_array_per_node.itemTable(node_id)[node_table.nbRows()], AssertError);
+        REQUIRE_THROWS_AS(face_array_per_node.itemTable(node_id)[0][face_array_per_node.sizeOfArrays()], AssertError);
+        REQUIRE_THROWS_AS(face_array_per_node.itemTable(node_id)[0][face_array_per_node.sizeOfArrays()] = 2,
                           AssertError);
       }
 
@@ -838,13 +930,13 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
         REQUIRE_THROWS_AS(edge_array_per_cell(invalid_cell_id, 0), AssertError);
       }
       if (connectivity.numberOfCells() > 0) {
-        CellId cell_id          = 0;
-        const auto& edge_arrays = edge_array_per_cell.itemArrays(cell_id);
-        REQUIRE_THROWS_AS(edge_array_per_cell(cell_id, edge_arrays.size()), AssertError);
-        REQUIRE_THROWS_AS(edge_arrays[edge_arrays.size()], AssertError);
-        REQUIRE_THROWS_AS(edge_array_per_cell.itemArrays(cell_id)[edge_arrays.size()], AssertError);
-        REQUIRE_THROWS_AS(edge_array_per_cell.itemArrays(cell_id)[0][edge_array_per_cell.sizeOfArrays()], AssertError);
-        REQUIRE_THROWS_AS(edge_array_per_cell.itemArrays(cell_id)[0][edge_array_per_cell.sizeOfArrays()] == 2,
+        CellId cell_id         = 0;
+        const auto& cell_table = edge_array_per_cell.itemTable(cell_id);
+        REQUIRE_THROWS_AS(edge_array_per_cell(cell_id, cell_table.nbRows()), AssertError);
+        REQUIRE_THROWS_AS(cell_table[cell_table.nbRows()], AssertError);
+        REQUIRE_THROWS_AS(edge_array_per_cell.itemTable(cell_id)[cell_table.nbRows()], AssertError);
+        REQUIRE_THROWS_AS(edge_array_per_cell.itemTable(cell_id)[0][edge_array_per_cell.sizeOfArrays()], AssertError);
+        REQUIRE_THROWS_AS(edge_array_per_cell.itemTable(cell_id)[0][edge_array_per_cell.sizeOfArrays()] == 2,
                           AssertError);
       }
 
@@ -854,13 +946,13 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
         REQUIRE_THROWS_AS(node_array_per_edge(invalid_edge_id, 0), AssertError);
       }
       if (connectivity.numberOfEdges() > 0) {
-        EdgeId edge_id          = 0;
-        const auto& node_arrays = node_array_per_edge.itemArrays(edge_id);
-        REQUIRE_THROWS_AS(node_array_per_edge(edge_id, node_arrays.size()), AssertError);
-        REQUIRE_THROWS_AS(node_arrays[node_arrays.size()], AssertError);
-        REQUIRE_THROWS_AS(node_array_per_edge.itemArrays(edge_id)[node_arrays.size()], AssertError);
-        REQUIRE_THROWS_AS(node_array_per_edge.itemArrays(edge_id)[0][node_array_per_edge.sizeOfArrays()], AssertError);
-        REQUIRE_THROWS_AS(node_array_per_edge.itemArrays(edge_id)[0][node_array_per_edge.sizeOfArrays()] = 2,
+        EdgeId edge_id         = 0;
+        const auto& edge_table = node_array_per_edge.itemTable(edge_id);
+        REQUIRE_THROWS_AS(node_array_per_edge(edge_id, edge_table.nbRows()), AssertError);
+        REQUIRE_THROWS_AS(edge_table[edge_table.nbRows()], AssertError);
+        REQUIRE_THROWS_AS(node_array_per_edge.itemTable(edge_id)[edge_table.nbRows()], AssertError);
+        REQUIRE_THROWS_AS(node_array_per_edge.itemTable(edge_id)[0][node_array_per_edge.sizeOfArrays()], AssertError);
+        REQUIRE_THROWS_AS(node_array_per_edge.itemTable(edge_id)[0][node_array_per_edge.sizeOfArrays()] = 2,
                           AssertError);
       }
     }
diff --git a/tests/test_Table.cpp b/tests/test_Table.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7b736f3d76f87129b45885fed737a361bf73562e
--- /dev/null
+++ b/tests/test_Table.cpp
@@ -0,0 +1,171 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <utils/PugsAssert.hpp>
+#include <utils/Table.hpp>
+#include <utils/Types.hpp>
+
+// Instantiate to ensure full coverage is performed
+template class Table<int>;
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("Table", "[utils]")
+{
+  Table<int> a(4, 3);
+  REQUIRE(a.nbRows() == 4);
+  REQUIRE(a.nbColumns() == 3);
+  REQUIRE(a[0].size() == 3);
+  REQUIRE(a[1].size() == 3);
+  REQUIRE(a[2].size() == 3);
+  REQUIRE(a[3].size() == 3);
+
+  for (size_t i = 0; i < a.nbRows(); ++i) {
+    for (size_t j = 0; j < a.nbColumns(); ++j) {
+      a(i, j) = 2 * i + j;
+    }
+  }
+
+  REQUIRE(((a(0, 0) == 0) and (a(1, 0) == 2) and (a(2, 0) == 4) and (a(3, 0) == 6) and   //
+           (a(0, 1) == 1) and (a(1, 1) == 3) and (a(2, 1) == 5) and (a(3, 1) == 7) and   //
+           (a(0, 2) == 2) and (a(1, 2) == 4) and (a(2, 2) == 6) and (a(3, 2) == 8)));
+
+  SECTION("checking for rows")
+  {
+    REQUIRE(((a[0][0] == 0) and (a[1][0] == 2) and (a[2][0] == 4) and (a[3][0] == 6) and   //
+             (a[0][1] == 1) and (a[1][1] == 3) and (a[2][1] == 5) and (a[3][1] == 7) and   //
+             (a[0][2] == 2) and (a[1][2] == 4) and (a[2][2] == 6) and (a[3][2] == 8)));
+
+    a[2][1] = 17;
+    REQUIRE(a[2][1] == a(2, 1));
+    a[2][1] = 5;
+  }
+
+  SECTION("checking for copies")
+  {
+    Table<const int> b{a};
+
+    REQUIRE(((b(0, 0) == 0) and (b(1, 0) == 2) and (b(2, 0) == 4) and (b(3, 0) == 6) and   //
+             (b(0, 1) == 1) and (b(1, 1) == 3) and (b(2, 1) == 5) and (b(3, 1) == 7) and   //
+             (b(0, 2) == 2) and (b(1, 2) == 4) and (b(2, 2) == 6) and (b(3, 2) == 8)));
+
+    Table<int> c{a};
+
+    REQUIRE(((c(0, 0) == 0) and (c(1, 0) == 2) and (c(2, 0) == 4) and (c(3, 0) == 6) and   //
+             (c(0, 1) == 1) and (c(1, 1) == 3) and (c(2, 1) == 5) and (c(3, 1) == 7) and   //
+             (c(0, 2) == 2) and (c(1, 2) == 4) and (c(2, 2) == 6) and (c(3, 2) == 8)));
+
+    Table<int> d = std::move(c);
+
+    REQUIRE(((d(0, 0) == 0) and (d(1, 0) == 2) and (d(2, 0) == 4) and (d(3, 0) == 6) and   //
+             (d(0, 1) == 1) and (d(1, 1) == 3) and (d(2, 1) == 5) and (d(3, 1) == 7) and   //
+             (d(0, 2) == 2) and (d(1, 2) == 4) and (d(2, 2) == 6) and (d(3, 2) == 8)));
+
+    a(2, 2) = 17;
+
+    REQUIRE(a(2, 2) == 17);
+    REQUIRE(b(2, 2) == 17);
+    REQUIRE(d(2, 2) == 17);
+
+    a(2, 2) = 6;
+
+    REQUIRE(a(2, 2) == 6);
+    REQUIRE(b(2, 2) == 6);
+    REQUIRE(d(2, 2) == 6);
+  }
+
+  SECTION("checking for fill")
+  {
+    Table<int> b(3, 2);
+    b.fill(3);
+
+    REQUIRE(((b(0, 0) == 3) and (b(1, 0) == 3) and (b(2, 0) == 3) and   //
+             (b(0, 1) == 3) and (b(1, 1) == 3) and (b(2, 1) == 3)));
+  }
+
+  SECTION("checking for affectations (shallow copy)")
+  {
+    Table<const int> b;
+    b = a;
+
+    REQUIRE(((b(0, 0) == 0) and (b(1, 0) == 2) and (b(2, 0) == 4) and (b(3, 0) == 6) and   //
+             (b(0, 1) == 1) and (b(1, 1) == 3) and (b(2, 1) == 5) and (b(3, 1) == 7) and   //
+             (b(0, 2) == 2) and (b(1, 2) == 4) and (b(2, 2) == 6) and (b(3, 2) == 8)));
+
+    Table<int> c;
+    c = a;
+
+    REQUIRE(((c(0, 0) == 0) and (c(1, 0) == 2) and (c(2, 0) == 4) and (c(3, 0) == 6) and   //
+             (c(0, 1) == 1) and (c(1, 1) == 3) and (c(2, 1) == 5) and (c(3, 1) == 7) and   //
+             (c(0, 2) == 2) and (c(1, 2) == 4) and (c(2, 2) == 6) and (c(3, 2) == 8)));
+
+    Table<int> d;
+    d = std::move(c);
+
+    REQUIRE(((d(0, 0) == 0) and (d(1, 0) == 2) and (d(2, 0) == 4) and (d(3, 0) == 6) and   //
+             (d(0, 1) == 1) and (d(1, 1) == 3) and (d(2, 1) == 5) and (d(3, 1) == 7) and   //
+             (d(0, 2) == 2) and (d(1, 2) == 4) and (d(2, 2) == 6) and (d(3, 2) == 8)));
+  }
+
+  SECTION("checking for affectations (deep copy)")
+  {
+    Table<int> b(copy(a));
+
+    REQUIRE(((b(0, 0) == 0) and (b(1, 0) == 2) and (b(2, 0) == 4) and (b(3, 0) == 6) and   //
+             (b(0, 1) == 1) and (b(1, 1) == 3) and (b(2, 1) == 5) and (b(3, 1) == 7) and   //
+             (b(0, 2) == 2) and (b(1, 2) == 4) and (b(2, 2) == 6) and (b(3, 2) == 8)));
+
+    b.fill(2);
+
+    REQUIRE(((b(0, 0) == 2) and (b(1, 0) == 2) and (b(2, 0) == 2) and (b(3, 0) == 2) and   //
+             (b(0, 1) == 2) and (b(1, 1) == 2) and (b(2, 1) == 2) and (b(3, 1) == 2) and   //
+             (b(0, 2) == 2) and (b(1, 2) == 2) and (b(2, 2) == 2) and (b(3, 2) == 2)));
+
+    REQUIRE(((a(0, 0) == 0) and (a(1, 0) == 2) and (a(2, 0) == 4) and (a(3, 0) == 6) and   //
+             (a(0, 1) == 1) and (a(1, 1) == 3) and (a(2, 1) == 5) and (a(3, 1) == 7) and   //
+             (a(0, 2) == 2) and (a(1, 2) == 4) and (a(2, 2) == 6) and (a(3, 2) == 8)));
+
+    Table<int> c;
+    c = a;
+
+    REQUIRE(((c(0, 0) == 0) and (c(1, 0) == 2) and (c(2, 0) == 4) and (c(3, 0) == 6) and   //
+             (c(0, 1) == 1) and (c(1, 1) == 3) and (c(2, 1) == 5) and (c(3, 1) == 7) and   //
+             (c(0, 2) == 2) and (c(1, 2) == 4) and (c(2, 2) == 6) and (c(3, 2) == 8)));
+
+    c = copy(b);
+
+    REQUIRE(((c(0, 0) == 2) and (c(1, 0) == 2) and (c(2, 0) == 2) and (c(3, 0) == 2) and   //
+             (c(0, 1) == 2) and (c(1, 1) == 2) and (c(2, 1) == 2) and (c(3, 1) == 2) and   //
+             (c(0, 2) == 2) and (c(1, 2) == 2) and (c(2, 2) == 2) and (c(3, 2) == 2)));
+  }
+
+  SECTION("checking for Kokkos::View encaspulation")
+  {
+    {
+      Kokkos::View<double**> kokkos_view("anonymous", 10, 3);
+      for (size_t i = 0; i < 10; ++i) {
+        for (size_t j = 0; j < 3; ++j) {
+          kokkos_view(i, j) = 3 * i + j;
+        }
+      }
+
+      Table table = encapsulate(kokkos_view);
+
+      REQUIRE(table.nbRows() == kokkos_view.extent(0));
+      REQUIRE(table.nbColumns() == kokkos_view.extent(1));
+      for (size_t i = 0; i < table.nbRows(); ++i) {
+        for (size_t j = 0; j < table.nbColumns(); ++j) {
+          REQUIRE(&table(i, j) == &kokkos_view(i, j));
+        }
+      }
+    }
+  }
+
+#ifndef NDEBUG
+  SECTION("checking for bounds violation")
+  {
+    REQUIRE_THROWS_AS(a(4, 0), AssertError);
+    REQUIRE_THROWS_AS(a(0, 3), AssertError);
+  }
+#endif   // NDEBUG
+}