diff --git a/src/mesh/ConnectivityComputer.cpp b/src/mesh/ConnectivityComputer.cpp
index 46deb5b95166d42bc3bb5ade5d3a02683157abfa..e08a638d750bb70a3e41709a4c70c3539199e048 100644
--- a/src/mesh/ConnectivityComputer.cpp
+++ b/src/mesh/ConnectivityComputer.cpp
@@ -48,9 +48,9 @@ ConnectivityComputer::_computeInverse(const ConnectivityMatrix& item_to_child_ma
   const size_t& number_of_items = item_to_child_matrix.numRows();
 
   for (unsigned int j = 0; j < number_of_items; ++j) {
-    const auto& item_to_child = item_to_child_matrix.rowConst(j);
-    for (unsigned int r = 0; r < item_to_child.length; ++r) {
-      child_to_item_vector_map[item_to_child(r)].push_back(j);
+    const auto& item_to_child = item_to_child_matrix[j];
+    for (unsigned int r = 0; r < item_to_child.size(); ++r) {
+      child_to_item_vector_map[item_to_child[r]].push_back(j);
     }
   }
 
@@ -81,18 +81,18 @@ ConnectivityComputer::computeLocalItemNumberInChildItem(const ConnectivityType&
   constexpr ItemType item_type       = ReversedItemOfItem::item_type;
   constexpr ItemType child_item_type = ReversedItemOfItem::sub_item_type;
 
-  const ConnectivityMatrix& child_item_to_items_matrix = connectivity._getMatrix(child_item_type, item_type);
-  const ConnectivityMatrix& item_to_child_items_matrix = connectivity._getMatrix(item_type, child_item_type);
+  ItemToItemMatrix item_to_child_items_matrix = connectivity.template getItemToItemMatrix<item_type, child_item_type>();
+  ItemToItemMatrix child_item_to_items_matrix = connectivity.template getItemToItemMatrix<child_item_type, item_type>();
 
   WeakSubItemValuePerItem<unsigned short, ReversedItemOfItem> item_number_in_child_item(connectivity);
-  for (ItemIdT<item_type> r = 0; r < item_to_child_items_matrix.numRows(); ++r) {
-    const auto& item_to_child_items = item_to_child_items_matrix.rowConst(r);
-    for (unsigned short J = 0; J < item_to_child_items.length; ++J) {
-      const unsigned int j            = item_to_child_items(J);
-      const auto& child_item_to_items = child_item_to_items_matrix.rowConst(j);
-
-      for (unsigned int R = 0; R < child_item_to_items.length; ++R) {
-        if (child_item_to_items(R) == r) {
+  for (ItemIdT<item_type> r = 0; r < connectivity.template numberOf<item_type>(); ++r) {
+    const auto& item_to_child_items = item_to_child_items_matrix[r];
+    for (unsigned short J = 0; J < item_to_child_items_matrix[r].size(); ++J) {
+      ItemIdT j                       = item_to_child_items[J];
+      const auto& child_item_to_items = child_item_to_items_matrix[j];
+
+      for (unsigned int R = 0; R < child_item_to_items.size(); ++R) {
+        if (child_item_to_items[R] == r) {
           item_number_in_child_item(r, J) = R;
           break;
         }
diff --git a/src/mesh/ConnectivityMatrix.hpp b/src/mesh/ConnectivityMatrix.hpp
index 576f7935fb11dce55ef20f9fdc5610bfbd27bbff..11df5978bbeb3eced43fa57fed878bdf60327916 100644
--- a/src/mesh/ConnectivityMatrix.hpp
+++ b/src/mesh/ConnectivityMatrix.hpp
@@ -4,70 +4,78 @@
 #include <utils/Array.hpp>
 #include <utils/PugsUtils.hpp>
 
-#include <Kokkos_StaticCrsGraph.hpp>
-
 class ConnectivityMatrix
 {
  private:
-  using HostMatrix = Kokkos::StaticCrsGraph<unsigned int, Kokkos::HostSpace>;
-  HostMatrix m_host_matrix;
+  using IndexType = uint32_t;
+
+  Array<const IndexType> m_row_map;
+  Array<const IndexType> m_column_indices;
 
   bool m_is_built{false};
 
  public:
-  using HostRowType = HostMatrix::row_map_type;
-
   const bool&
   isBuilt() const
   {
     return m_is_built;
   }
 
-  auto
+  const Array<const IndexType>&
   entries() const
   {
-    return encapsulate(m_host_matrix.entries);
+    return m_column_indices;
   }
 
   const auto&
   rowsMap() const
   {
-    return m_host_matrix.row_map;
+    return m_row_map;
   }
 
   PUGS_INLINE
-  auto
+  size_t
   numEntries() const
   {
-    return m_host_matrix.entries.extent(0);
+    return m_column_indices.size();
   }
 
   PUGS_INLINE
-  auto
+  size_t
   numRows() const
   {
-    return m_host_matrix.numRows();
+    return m_row_map.size() - 1;
   }
 
   PUGS_INLINE
   auto
-  rowConst(size_t j) const
-  {
-    return m_host_matrix.rowConst(j);
-  }
-
-  PUGS_INLINE
-  const auto&
-  rowMap(size_t j) const
+  operator[](size_t j) const
   {
-    return m_host_matrix.row_map[j];
+    return subArrayView(m_column_indices, m_row_map[j], m_row_map[j + 1] - m_row_map[j]);
   }
 
   PUGS_INLINE
-  ConnectivityMatrix(const std::vector<std::vector<unsigned int>>& initializer)
-    : m_host_matrix{Kokkos::create_staticcrsgraph<HostMatrix>("connectivity_matrix", initializer)}, m_is_built{true}
+  ConnectivityMatrix(const std::vector<std::vector<unsigned int>>& initializer) : m_is_built{true}
   {
-    ;
+    m_row_map = [&] {
+      Array<uint32_t> row_map(initializer.size() + 1);
+      row_map[0] = 0;
+      for (size_t i = 0; i < initializer.size(); ++i) {
+        row_map[i + 1] = row_map[i] + initializer[i].size();
+      }
+      return row_map;
+    }();
+
+    m_column_indices = [&] {
+      Array<uint32_t> column_indices(m_row_map[m_row_map.size() - 1]);
+      size_t index = 0;
+      for (const auto& row : initializer) {
+        for (const auto& col_index : row) {
+          column_indices[index++] = col_index;
+        }
+      }
+      return column_indices;
+    }();
   }
 
   ConnectivityMatrix& operator=(const ConnectivityMatrix&) = default;
diff --git a/src/mesh/ItemArray.hpp b/src/mesh/ItemArray.hpp
index 90ba27629874251874801129dd5d250335c393f2..b1d4fcb1ae384fc78c07165542bb54170b3b47b4 100644
--- a/src/mesh/ItemArray.hpp
+++ b/src/mesh/ItemArray.hpp
@@ -93,7 +93,7 @@ class ItemArray
   }
 
   template <ItemType item_t>
-  PUGS_INLINE Array<DataType>
+  PUGS_INLINE typename Table<DataType>::UnsafeRowView
   operator[](const ItemIdT<item_t>& item_id) const noexcept(NO_ASSERT)
   {
     static_assert(item_t == item_type, "invalid ItemId type");
diff --git a/src/mesh/ItemArrayUtils.hpp b/src/mesh/ItemArrayUtils.hpp
index ca339de9a10809d939c8924844c507442c83b8ef..6afa39c167c9d2fa95ab5ef53e502df27359dce7 100644
--- a/src/mesh/ItemArrayUtils.hpp
+++ b/src/mesh/ItemArrayUtils.hpp
@@ -42,7 +42,7 @@ min(const ItemArray<DataType, item_type, ConnectivityPtr>& item_value)
     operator()(const index_type& i_item, data_type& data) const
     {
       if (m_is_owned[i_item]) {
-        Array array = m_item_value[i_item];
+        auto array = m_item_value[i_item];
         for (size_t i = 0; i < m_item_value.sizeOfArrays(); ++i) {
           if (array[i] < data) {
             data = array[i];
@@ -141,7 +141,7 @@ max(const ItemArray<DataType, item_type, ConnectivityPtr>& item_value)
     operator()(const index_type& i_item, data_type& data) const
     {
       if (m_is_owned[i_item]) {
-        Array array = m_item_value[i_item];
+        auto array = m_item_value[i_item];
         for (size_t i = 0; i < m_item_value.sizeOfArrays(); ++i) {
           if (array[i] > data) {
             data = array[i];
@@ -239,7 +239,7 @@ sum(const ItemArray<DataType, item_type, ConnectivityPtr>& item_value)
     operator()(const index_type& i_item, data_type& data) const
     {
       if (m_is_owned[i_item]) {
-        Array array = m_item_value[i_item];
+        auto array = m_item_value[i_item];
         for (size_t i = 0; i < m_item_value.sizeOfArrays(); ++i) {
           data += array[i];
         }
diff --git a/src/mesh/ItemToItemMatrix.hpp b/src/mesh/ItemToItemMatrix.hpp
index 1cb20d3ee595cf0963b95b4cc67f201cf2cddb5c..acee8ecd8e85f081ca0e640e2f408a541ac9157a 100644
--- a/src/mesh/ItemToItemMatrix.hpp
+++ b/src/mesh/ItemToItemMatrix.hpp
@@ -5,6 +5,7 @@
 #include <mesh/ItemId.hpp>
 #include <mesh/ItemType.hpp>
 
+#include <utils/Exceptions.hpp>
 #include <utils/PugsUtils.hpp>
 
 template <ItemType SourceItemType, ItemType TargetItemType>
@@ -14,33 +15,35 @@ class ItemToItemMatrix
   using SourceItemId = ItemIdT<SourceItemType>;
   using TargetItemId = ItemIdT<TargetItemType>;
 
-  template <typename RowType>
   class SubItemList
   {
    private:
-    const RowType m_row;
+    const ConnectivityMatrix& m_connectivity_matrix;
+    SourceItemId m_source_item_id;
+    size_t m_size;
 
    public:
     PUGS_INLINE
     size_t
     size() const
     {
-      return m_row.length;
+      return m_size;
     }
 
-    PUGS_INLINE
-    TargetItemId
+    PUGS_INLINE TargetItemId
     operator[](size_t j) const
     {
-      Assert(j < m_row.length);
-      return m_row(j);
+      Assert(j < m_size);
+      return m_connectivity_matrix.entries()[m_connectivity_matrix.rowsMap()[m_source_item_id] + j];
     }
 
     PUGS_INLINE
-    SubItemList(const RowType& row) : m_row{row}
-    {
-      ;
-    }
+    SubItemList(const ConnectivityMatrix& connectivity_matrix, SourceItemId source_item_id)
+      : m_connectivity_matrix{connectivity_matrix},
+        m_source_item_id(source_item_id),
+        m_size(m_connectivity_matrix.rowsMap()[m_source_item_id + 1] -
+               m_connectivity_matrix.rowsMap()[m_source_item_id])
+    {}
 
     PUGS_INLINE
     SubItemList& operator=(const SubItemList&) = default;
@@ -62,12 +65,6 @@ class ItemToItemMatrix
   const ConnectivityMatrix& m_connectivity_matrix;
 
  public:
-  auto
-  numberOfEntries() const
-  {
-    return m_connectivity_matrix.numEntries();
-  }
-
   auto
   entries() const
   {
@@ -79,25 +76,18 @@ class ItemToItemMatrix
   operator[](const SourceItemId& source_id) const
   {
     Assert(source_id < m_connectivity_matrix.numRows());
-    using RowType = decltype(m_connectivity_matrix.rowConst(source_id));
-    return SubItemList<RowType>(m_connectivity_matrix.rowConst(source_id));
+    return SubItemList(m_connectivity_matrix, source_id);
   }
 
   template <typename IndexType>
-  PUGS_INLINE const auto&
-  operator[](const IndexType& source_id) const
+  PUGS_INLINE void
+  operator[](const IndexType&) const
   {
-    Assert(source_id < m_connectivity_matrix.numRows());
     static_assert(std::is_same_v<IndexType, SourceItemId>, "ItemToItemMatrix must be indexed using correct ItemId");
-    using RowType = decltype(m_connectivity_matrix.rowConst(source_id));
-    return SubItemList<RowType>(m_connectivity_matrix.rowConst(source_id));
   }
 
   PUGS_INLINE
-  ItemToItemMatrix(const ConnectivityMatrix& connectivity_matrix) : m_connectivity_matrix{connectivity_matrix}
-  {
-    ;
-  }
+  ItemToItemMatrix(const ConnectivityMatrix& connectivity_matrix) : m_connectivity_matrix{connectivity_matrix} {}
 
   PUGS_INLINE
   ItemToItemMatrix& operator=(const ItemToItemMatrix&) = default;
diff --git a/src/mesh/SubItemArrayPerItem.hpp b/src/mesh/SubItemArrayPerItem.hpp
index e0bc1e611c10ccecb088c272c75d1d5e2c7d9d38..cb88e85d76f758262316b547140e85e000b669ca 100644
--- a/src/mesh/SubItemArrayPerItem.hpp
+++ b/src/mesh/SubItemArrayPerItem.hpp
@@ -33,7 +33,7 @@ class SubItemArrayPerItem
 
   ConnectivityPtr m_connectivity_ptr;
 
-  typename ConnectivityMatrix::HostRowType m_host_row_map;
+  Array<const uint32_t> m_row_map;
   Table<DataType> m_values;
 
   // Allow const std:shared_ptr version to access our data
@@ -49,27 +49,25 @@ class SubItemArrayPerItem
   friend SubItemArrayPerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityWeakPtr>;
 
  public:
-  friend PUGS_INLINE SubItemArrayPerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityPtr>
+  [[nodiscard]] friend PUGS_INLINE SubItemArrayPerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityPtr>
   copy(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& source)
   {
     SubItemArrayPerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityPtr> image;
 
     image.m_connectivity_ptr = source.m_connectivity_ptr;
-    image.m_host_row_map     = source.m_host_row_map;
+    image.m_row_map          = source.m_row_map;
     image.m_values           = copy(source.m_values);
 
     return image;
   }
 
-  PUGS_INLINE
-  bool
+  [[nodiscard]] PUGS_INLINE bool
   isBuilt() const noexcept
   {
     return m_connectivity_ptr.use_count() != 0;
   }
 
-  PUGS_INLINE
-  std::shared_ptr<const IConnectivity>
+  [[nodiscard]] PUGS_INLINE std::shared_ptr<const IConnectivity>
   connectivity_ptr() const noexcept
   {
     if constexpr (std::is_same_v<ConnectivityPtr, ConnectivitySharedPtr>) {
@@ -80,18 +78,18 @@ class SubItemArrayPerItem
   }
 
   template <typename IndexType, typename SubIndexType>
-  PUGS_FORCEINLINE Array<DataType>
+  [[nodiscard]] PUGS_INLINE typename Table<DataType>::UnsafeRowView
   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 m_values[m_host_row_map(size_t{item_id}) + i];
+    Assert(i + m_row_map[size_t{item_id}] < m_row_map[size_t{item_id} + 1]);
+    return m_values[m_row_map[size_t{item_id}] + i];
   }
 
   template <ItemType item_t>
-  Table<DataType>
+  [[nodiscard]] PUGS_INLINE typename Table<DataType>::UnsafeTableView
   operator[](const ItemIdT<item_t>& item_id) const noexcept(NO_ASSERT)
   {
     static_assert(item_t == item_type, "invalid ItemId type");
@@ -99,7 +97,7 @@ class SubItemArrayPerItem
   }
 
   template <typename ArrayIndexType>
-  Array<DataType>
+  [[nodiscard]] PUGS_INLINE typename Table<DataType>::UnsafeRowView
   operator[](const ArrayIndexType& i) const noexcept(NO_ASSERT)
   {
     static_assert(std::is_integral_v<ArrayIndexType>, "index must be an integral type");
@@ -108,24 +106,22 @@ class SubItemArrayPerItem
     return m_values[i];
   }
 
-  PUGS_INLINE
-  size_t
+  [[nodiscard]] PUGS_INLINE size_t
   numberOfArrays() const noexcept(NO_ASSERT)
   {
     Assert(this->isBuilt());
     return m_values.numberOfRows();
   }
 
-  PUGS_INLINE
-  size_t
+  [[nodiscard]] PUGS_INLINE size_t
   numberOfItems() const noexcept(NO_ASSERT)
   {
     Assert(this->isBuilt());
-    Assert(m_host_row_map.extent(0) > 0);
-    return m_host_row_map.extent(0) - 1;
+    Assert(m_row_map.size() > 0);
+    return m_row_map.size() - 1;
   }
 
-  PUGS_INLINE size_t
+  [[nodiscard]] PUGS_INLINE size_t
   sizeOfArrays() const noexcept(NO_ASSERT)
   {
     Assert(this->isBuilt());
@@ -133,17 +129,16 @@ class SubItemArrayPerItem
   }
 
   template <typename IndexType>
-  PUGS_INLINE size_t
+  [[nodiscard]] PUGS_INLINE size_t
   numberOfSubArrays(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());
-    return m_host_row_map(size_t{item_id} + 1) - m_host_row_map(size_t{item_id});
+    return m_row_map[size_t{item_id} + 1] - m_row_map[size_t{item_id}];
   }
 
-  PUGS_INLINE
-  void
+  PUGS_INLINE void
   fill(const DataType& data) const noexcept
   {
     static_assert(not std::is_const_v<DataType>, "Cannot modify SubItemArrayPerItem of const");
@@ -151,15 +146,15 @@ class SubItemArrayPerItem
   }
 
   template <typename IndexType>
-  PUGS_INLINE Table<DataType>
+  [[nodiscard]] PUGS_INLINE typename Table<DataType>::UnsafeTableView
   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 subTable(m_values, item_begin, item_end - item_begin, 0, this->sizeOfArrays());
+    const auto& item_begin = m_row_map[size_t{item_id}];
+    const auto& item_end   = m_row_map[size_t{item_id} + 1];
+    return subTableView(m_values, item_begin, item_end - item_begin, 0, this->sizeOfArrays());
   }
 
   template <typename DataType2, typename ConnectivityPtr2>
@@ -173,8 +168,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_values       = sub_item_array_per_item.m_values;
+    m_row_map = sub_item_array_per_item.m_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>) {
@@ -204,8 +199,8 @@ class SubItemArrayPerItem
 
     ConnectivityMatrix connectivity_matrix = connectivity._getMatrix(item_type, sub_item_type);
 
-    m_host_row_map = connectivity_matrix.rowsMap();
-    m_values       = Table<std::remove_const_t<DataType>>(connectivity_matrix.numEntries(), size_of_arrays);
+    m_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/SubItemArrayPerItemUtils.hpp b/src/mesh/SubItemArrayPerItemUtils.hpp
index 5e4f122e3d8dc7a06858091a2e8f8baa402c95e4..1c5ad24683c5527e8028896348ba73ff8160af74 100644
--- a/src/mesh/SubItemArrayPerItemUtils.hpp
+++ b/src/mesh/SubItemArrayPerItemUtils.hpp
@@ -44,7 +44,7 @@ min(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& sub_item_a
     operator()(const index_type& item_id, data_type& data) const
     {
       if (m_is_owned[item_id]) {
-        Table sub_item_table = m_sub_item_array_per_item.itemTable(item_id);
+        auto sub_item_table = m_sub_item_array_per_item.itemTable(item_id);
         for (size_t i = 0; i < sub_item_table.numberOfRows(); ++i) {
           for (size_t j = 0; j < sub_item_table.numberOfColumns(); ++j) {
             if (sub_item_table(i, j) < data) {
@@ -146,7 +146,7 @@ max(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& sub_item_a
     operator()(const index_type& item_id, data_type& data) const
     {
       if (m_is_owned[item_id]) {
-        Table sub_item_table = m_sub_item_array_per_item.itemTable(item_id);
+        auto sub_item_table = m_sub_item_array_per_item.itemTable(item_id);
         for (size_t i = 0; i < sub_item_table.numberOfRows(); ++i) {
           for (size_t j = 0; j < sub_item_table.numberOfColumns(); ++j) {
             if (sub_item_table(i, j) > data) {
@@ -247,7 +247,7 @@ sum(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& item_array
     operator()(const index_type& item_id, data_type& data) const
     {
       if (m_is_owned[item_id]) {
-        Table sub_item_table = m_sub_item_array_per_item.itemTable(item_id);
+        auto sub_item_table = m_sub_item_array_per_item.itemTable(item_id);
         for (size_t i = 0; i < sub_item_table.numberOfRows(); ++i) {
           for (size_t j = 0; j < sub_item_table.numberOfColumns(); ++j) {
             data += sub_item_table(i, j);
@@ -341,8 +341,8 @@ isSynchronized(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>&
     synchronize(sub_item_array_per_item_copy);
 
     for (size_t i = 0; i < sub_item_array_per_item.numberOfArrays(); ++i) {
-      Array sub_item_array      = sub_item_array_per_item[i];
-      Array sub_item_array_copy = sub_item_array_per_item_copy[i];
+      auto sub_item_array      = sub_item_array_per_item[i];
+      auto sub_item_array_copy = sub_item_array_per_item_copy[i];
       for (size_t j = 0; j < sub_item_array.size(); ++j) {
         if (sub_item_array_copy[j] != sub_item_array[j]) {
           is_synchronized = false;
diff --git a/src/mesh/SubItemValuePerItem.hpp b/src/mesh/SubItemValuePerItem.hpp
index 90ef57b2ca0221752a5890f0ede7f481a6fe4f70..2b6ad94604551cabf495e70e758b7eadf00833c6 100644
--- a/src/mesh/SubItemValuePerItem.hpp
+++ b/src/mesh/SubItemValuePerItem.hpp
@@ -32,7 +32,7 @@ class SubItemValuePerItem
 
   ConnectivityPtr m_connectivity_ptr;
 
-  typename ConnectivityMatrix::HostRowType m_host_row_map;
+  Array<const uint32_t> m_row_map;
   Array<DataType> m_values;
 
   // Allow const std:shared_ptr version to access our data
@@ -54,7 +54,7 @@ class SubItemValuePerItem
     SubItemValuePerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityPtr> image;
 
     image.m_connectivity_ptr = source.m_connectivity_ptr;
-    image.m_host_row_map     = source.m_host_row_map;
+    image.m_row_map          = source.m_row_map;
     image.m_values           = copy(source.m_values);
     return image;
   }
@@ -84,12 +84,12 @@ class SubItemValuePerItem
     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 m_values[m_host_row_map(size_t{item_id}) + i];
+    Assert(i + m_row_map[size_t{item_id}] < m_row_map[size_t{item_id} + 1]);
+    return m_values[m_row_map[size_t{item_id}] + i];
   }
 
   template <ItemType item_t>
-  Array<DataType>
+  typename Array<DataType>::UnsafeArrayView
   operator[](const ItemIdT<item_t>& item_id) const noexcept(NO_ASSERT)
   {
     static_assert(item_t == item_type, "invalid ItemId type");
@@ -119,8 +119,8 @@ class SubItemValuePerItem
   numberOfItems() const noexcept(NO_ASSERT)
   {
     Assert(this->isBuilt());
-    Assert(m_host_row_map.extent(0) > 0);
-    return m_host_row_map.extent(0) - 1;
+    Assert(m_row_map.size() > 0);
+    return m_row_map.size() - 1;
   }
 
   template <typename IndexType>
@@ -130,7 +130,7 @@ class SubItemValuePerItem
     static_assert(std::is_same_v<IndexType, ItemId>, "index must be an ItemId");
     Assert(this->isBuilt());
     Assert(item_id < this->numberOfItems());
-    return m_host_row_map(size_t{item_id} + 1) - m_host_row_map(size_t{item_id});
+    return m_row_map[size_t{item_id} + 1] - m_row_map[size_t{item_id}];
   }
 
   PUGS_INLINE
@@ -142,15 +142,15 @@ class SubItemValuePerItem
   }
 
   template <typename IndexType>
-  PUGS_INLINE Array<DataType>
+  PUGS_INLINE typename Array<DataType>::UnsafeArrayView
   itemArray(const 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 subArray(m_values, item_begin, item_end - item_begin);
+    const auto& item_begin = m_row_map[size_t{item_id}];
+    const auto& item_end   = m_row_map[size_t{item_id} + 1];
+    return subArrayView(m_values, item_begin, item_end - item_begin);
   }
 
   template <typename DataType2, typename ConnectivityPtr2>
@@ -164,8 +164,8 @@ class SubItemValuePerItem
     static_assert(((std::is_const_v<DataType2> and std::is_const_v<DataType>) or not std::is_const_v<DataType2>),
                   "Cannot assign SubItemValuePerItem of const data to "
                   "SubItemValuePerItem of non-const data");
-    m_host_row_map = sub_item_value_per_item.m_host_row_map;
-    m_values       = sub_item_value_per_item.m_values;
+    m_row_map = sub_item_value_per_item.m_row_map;
+    m_values  = sub_item_value_per_item.m_values;
 
     if constexpr (std::is_same_v<ConnectivityPtr, ConnectivitySharedPtr> and
                   std::is_same_v<ConnectivityPtr2, ConnectivityWeakPtr>) {
@@ -194,8 +194,8 @@ class SubItemValuePerItem
 
     ConnectivityMatrix connectivity_matrix = connectivity._getMatrix(item_type, sub_item_type);
 
-    m_host_row_map = connectivity_matrix.rowsMap();
-    m_values       = Array<std::remove_const_t<DataType>>(connectivity_matrix.numEntries());
+    m_row_map = connectivity_matrix.rowsMap();
+    m_values  = Array<std::remove_const_t<DataType>>(connectivity_matrix.numEntries());
   }
 
   ~SubItemValuePerItem() = default;
diff --git a/src/mesh/SubItemValuePerItemUtils.hpp b/src/mesh/SubItemValuePerItemUtils.hpp
index 88f36d1cc15704a947d0990f0a2190c9c375c53c..8443047745680161747998a8cd114ffac8fe2647 100644
--- a/src/mesh/SubItemValuePerItemUtils.hpp
+++ b/src/mesh/SubItemValuePerItemUtils.hpp
@@ -44,7 +44,7 @@ min(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& sub_item_v
     operator()(const index_type& item_id, data_type& data) const
     {
       if (m_is_owned[item_id]) {
-        Array sub_item_array = m_sub_item_value_per_item.itemArray(item_id);
+        auto sub_item_array = m_sub_item_value_per_item.itemArray(item_id);
         for (size_t i = 0; i < sub_item_array.size(); ++i) {
           if (sub_item_array[i] < data) {
             data = sub_item_array[i];
@@ -144,7 +144,7 @@ max(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& sub_item_v
     operator()(const index_type& item_id, data_type& data) const
     {
       if (m_is_owned[item_id]) {
-        Array sub_item_array = m_sub_item_value_per_item.itemArray(item_id);
+        auto sub_item_array = m_sub_item_value_per_item.itemArray(item_id);
         for (size_t i = 0; i < sub_item_array.size(); ++i) {
           if (sub_item_array[i] > data) {
             data = sub_item_array[i];
@@ -243,7 +243,7 @@ sum(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& item_value
     operator()(const index_type& item_id, data_type& data) const
     {
       if (m_is_owned[item_id]) {
-        Array sub_item_array = m_sub_item_value_per_item.itemArray(item_id);
+        auto sub_item_array = m_sub_item_value_per_item.itemArray(item_id);
         for (size_t i = 0; i < sub_item_array.size(); ++i) {
           data += sub_item_array[i];
         }
diff --git a/src/output/VTKWriter.cpp b/src/output/VTKWriter.cpp
index 7b407256ac7cfac5e354330a99a081b73e822540..5e8b84f956d31e76635f51fa78dcd9df7ecef233 100644
--- a/src/output/VTKWriter.cpp
+++ b/src/output/VTKWriter.cpp
@@ -30,7 +30,7 @@ class ICharArrayEmbedder
 template <typename InputDataT>
 class CharArrayEmbedder : public ICharArrayEmbedder
 {
-  CastArray<InputDataT, char> m_char_cast_array;
+  CastArray<InputDataT, const char> m_char_cast_array;
 
  public:
   size_t
diff --git a/src/scheme/DiscreteFunctionP0Vector.hpp b/src/scheme/DiscreteFunctionP0Vector.hpp
index 5f00d5964d3b0e060549a6337b62e40fe2b97f21..8eaa5d0a7580ddc65cb17e82d6f5551f86114e1c 100644
--- a/src/scheme/DiscreteFunctionP0Vector.hpp
+++ b/src/scheme/DiscreteFunctionP0Vector.hpp
@@ -107,7 +107,7 @@ class DiscreteFunctionP0Vector : public IDiscreteFunction
   }
 
   PUGS_FORCEINLINE
-  Array<DataType>
+  typename Table<DataType>::UnsafeRowView
   operator[](const CellId cell_id) const noexcept(NO_ASSERT)
   {
     return m_cell_arrays[cell_id];
@@ -202,8 +202,17 @@ class DiscreteFunctionP0Vector : public IDiscreteFunction
     Assert(f.size() == g.size());
     DiscreteFunctionP0<Dimension, double> result{f.m_mesh};
     parallel_for(
-      f.m_mesh->numberOfCells(),
-      PUGS_LAMBDA(CellId cell_id) { result[cell_id] = dot(Vector{f[cell_id]}, Vector{g[cell_id]}); });
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        const auto& f_cell_id = f[cell_id];
+        const auto& g_cell_id = g[cell_id];
+
+        double sum = 0;
+        for (size_t i = 0; i < f.size(); ++i) {
+          sum += f_cell_id[i] * g_cell_id[i];
+        }
+
+        result[cell_id] = sum;
+      });
 
     return result;
   }
diff --git a/src/utils/Array.hpp b/src/utils/Array.hpp
index ca8a6ea6b349b2918c54ccd8de354bff8e999a6f..d641dbdeffae9fa0b6434ba9c8ef2b30853c0ff0 100644
--- a/src/utils/Array.hpp
+++ b/src/utils/Array.hpp
@@ -18,6 +18,49 @@ class [[nodiscard]] Array
   using data_type  = DataType;
   using index_type = size_t;
 
+  class [[nodiscard]] UnsafeArrayView
+  {
+   private:
+    DataType* const m_values;
+    const size_t m_size;
+
+   public:
+    [[nodiscard]] PUGS_INLINE size_t size() const
+    {
+      return m_size;
+    }
+
+    [[nodiscard]] PUGS_INLINE DataType& operator[](size_t i) const
+    {
+      Assert(i < m_size, "invalid index");
+      return m_values[i];
+    }
+
+    PUGS_INLINE void fill(const DataType& data) const
+    {
+      for (size_t i = 0; i < m_size; ++i) {
+        m_values[i] = data;
+      }
+    }
+
+    UnsafeArrayView& operator=(const UnsafeArrayView&) = delete;
+    UnsafeArrayView& operator=(UnsafeArrayView&&) = delete;
+
+    UnsafeArrayView(const Array<DataType>& array, index_type begin, index_type size)
+      : m_values{&array[begin]}, m_size{size}
+    {
+      Assert((begin < array.size()) and (begin + size <= array.size()), "required view is not contained in the Array");
+    }
+
+    // To try to keep these views close to the initial array one
+    // forbids copy constructor and take benefits of C++-17 copy
+    // elisions.
+    UnsafeArrayView(const UnsafeArrayView&) = delete;
+
+    UnsafeArrayView()  = delete;
+    ~UnsafeArrayView() = default;
+  };
+
  private:
   Kokkos::View<DataType*> m_values;
 
@@ -25,12 +68,12 @@ class [[nodiscard]] Array
   friend Array<std::add_const_t<DataType>>;
 
  public:
-  PUGS_INLINE size_t size() const noexcept
+  [[nodiscard]] PUGS_INLINE size_t size() const noexcept
   {
     return m_values.extent(0);
   }
 
-  friend PUGS_INLINE Array<std::remove_const_t<DataType>> copy(const Array<DataType>& source)
+  [[nodiscard]] friend PUGS_INLINE Array<std::remove_const_t<DataType>> copy(const Array<DataType>& source)
   {
     Array<std::remove_const_t<DataType>> image(source.size());
     Kokkos::deep_copy(image.m_values, source.m_values);
@@ -41,7 +84,7 @@ class [[nodiscard]] Array
   friend PUGS_INLINE void copy_to(const Array<DataType>& source,
                                   const Array<std::remove_const_t<DataType>>& destination)
   {
-    Assert(source.size() == destination.size());
+    Assert(source.size() == destination.size(), "incompatible Array sizes");
     Kokkos::deep_copy(destination.m_values, source.m_values);
   }
 
@@ -49,13 +92,13 @@ class [[nodiscard]] Array
   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 typename Array<DataType2>::UnsafeArrayView subArrayView(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)
+  [[nodiscard]] PUGS_INLINE DataType& operator[](index_type i) const noexcept(NO_ASSERT)
   {
-    Assert(i < m_values.extent(0));
+    Assert(i < m_values.extent(0), "invalid index");
     return m_values[i];
   }
 
@@ -158,14 +201,12 @@ encapsulate(const Kokkos::View<DataType*, RT...>& values)
 }
 
 template <typename DataType>
-[[nodiscard]] PUGS_INLINE Array<DataType>
-subArray(const Array<DataType>& array,
-         typename Array<DataType>::index_type begin,
-         typename Array<DataType>::index_type size)
+[[nodiscard]] PUGS_INLINE typename Array<DataType>::UnsafeArrayView
+subArrayView(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)));
+  return typename Array<DataType>::UnsafeArrayView(array, begin, size);
 }
 
 // map, multimap, unordered_map and stack cannot be copied this way
diff --git a/src/utils/Table.hpp b/src/utils/Table.hpp
index cd380661b5b59d33f20812f969f9df7b93cd5a73..6171f1a242d0cce43137c4a24458f4c00059269c 100644
--- a/src/utils/Table.hpp
+++ b/src/utils/Table.hpp
@@ -1,8 +1,8 @@
 #ifndef TABLE_HPP
 #define TABLE_HPP
 
-#include <utils/Array.hpp>
 #include <utils/InvalidData.hpp>
+#include <utils/NaNHelper.hpp>
 #include <utils/PugsAssert.hpp>
 #include <utils/PugsMacros.hpp>
 #include <utils/PugsUtils.hpp>
@@ -25,24 +25,178 @@ class [[nodiscard]] Table
   friend Table<std::add_const_t<DataType>>;
 
  public:
-  PUGS_INLINE size_t numberOfRows() const noexcept
+  class [[nodiscard]] UnsafeRowView
+  {
+   private:
+    // @note for performances reasons it is very important to keep
+    // m_table as a reference. Kokkos' reference counting mechanism
+    // seems quite expensive
+    const Table<DataType>& m_table;
+    const size_t m_row;
+
+   public:
+    PUGS_INLINE size_t size() const noexcept
+    {
+      return m_table.numberOfColumns();
+    }
+
+    PUGS_INLINE
+    DataType& operator[](size_t i) const
+    {
+      Assert(i < m_table.numberOfColumns(), "invalid index");
+      return m_table(m_row, i);
+    }
+
+    PUGS_INLINE void fill(const DataType& data) const
+    {
+      for (size_t i = 0; i < this->size(); ++i) {
+        m_table(m_row, i) = data;
+      }
+    }
+
+    UnsafeRowView& operator=(const UnsafeRowView&) = delete;
+    UnsafeRowView& operator=(UnsafeRowView&&) = delete;
+
+    UnsafeRowView(const Table<DataType>& table, index_type row) : m_table{table}, m_row{row}
+    {
+      Assert(row < table.numberOfRows(), "required row view is not contained in the Table");
+    }
+
+    // To try to keep these views close to the initial array one
+    // forbids copy constructor and take benefits of C++-17 copy
+    // elisions.
+    UnsafeRowView(const UnsafeRowView&) = delete;
+
+    UnsafeRowView()  = delete;
+    ~UnsafeRowView() = default;
+  };
+
+  class [[nodiscard]] UnsafeTableView
+  {
+   private:
+    // @note for performances reasons it is very important to keep
+    // m_table as a reference. Kokkos' reference counting mechanism
+    // seems quite expensive
+    const Table<DataType>& m_table;
+    const size_t m_row_begin;
+    const size_t m_column_begin;
+    const size_t m_row_size;
+    const size_t m_column_size;
+
+   public:
+    class [[nodiscard]] RowView
+    {
+     private:
+      const UnsafeTableView& m_table_view;
+      const size_t m_row;
+
+     public:
+      [[nodiscard]] PUGS_INLINE size_t size() const
+      {
+        return m_table_view.numberOfColumns();
+      }
+
+      [[nodiscard]] PUGS_INLINE DataType& operator[](size_t i) const
+      {
+        Assert(i < m_table_view.numberOfColumns(), "invalid index");
+        return m_table_view(m_row, i);
+      }
+
+      PUGS_INLINE void fill(const DataType& data) const
+      {
+        for (size_t i = 0; i < this->size(); ++i) {
+          m_table_view(m_row, i) = data;
+        }
+      }
+
+      RowView(const UnsafeTableView& table_view, index_type row) : m_table_view{table_view}, m_row{row}
+      {
+        Assert(row < m_table_view.numberOfRows(), "required row view is not contained in the Table view");
+      }
+
+      // To try to keep these views close to the initial array one
+      // forbids copy constructor and take benefits of C++-17 copy
+      // elisions.
+      RowView(const RowView&) = delete;
+
+      RowView()  = delete;
+      ~RowView() = default;
+    };
+
+    [[nodiscard]] PUGS_INLINE size_t numberOfRows() const noexcept
+    {
+      return m_row_size;
+    }
+
+    [[nodiscard]] PUGS_INLINE size_t numberOfColumns() const noexcept
+    {
+      return m_column_size;
+    }
+
+    [[nodiscard]] PUGS_INLINE RowView operator[](size_t i) const
+    {
+      Assert(i < this->numberOfRows(), "invalid index");
+      return RowView(*this, i);
+    }
+
+    [[nodiscard]] PUGS_INLINE DataType& operator()(size_t i, size_t j) const
+    {
+      Assert(i < m_row_size, "invalid row index");
+      Assert(j < m_column_size, "invalid column index");
+
+      return m_table(m_row_begin + i, m_column_begin + j);
+    }
+
+    PUGS_INLINE void fill(const DataType& data) const
+    {
+      for (size_t i = 0; i < m_row_size; ++i) {
+        for (size_t j = 0; j < m_column_size; ++j) {
+          m_table(m_row_begin + i, m_column_begin + j) = data;
+        }
+      }
+    }
+    UnsafeTableView& operator=(const UnsafeTableView&) = delete;
+    UnsafeTableView& operator=(UnsafeTableView&&) = delete;
+
+    UnsafeTableView(const Table<DataType>& table, index_type row_begin, index_type row_size, index_type column_begin,
+                    index_type column_size)
+      : m_table{table},
+        m_row_begin{row_begin},
+        m_column_begin{column_begin},
+        m_row_size{row_size},
+        m_column_size{column_size}
+    {
+      Assert((row_begin < table.numberOfRows()) and (row_begin + row_size <= table.numberOfRows()),
+             "required view has rows not contained in the Table");
+      Assert((column_begin < table.numberOfColumns()) and (column_begin + column_size <= table.numberOfColumns()),
+             "required view has columns not contained in the Table");
+    }
+
+    // To try to keep these views close to the initial array one
+    // forbids copy constructor and take benefits of C++-17 copy
+    // elisions.
+    UnsafeTableView(const UnsafeTableView&) = delete;
+
+    UnsafeTableView()  = delete;
+    ~UnsafeTableView() = default;
+  };
+
+  [[nodiscard]] PUGS_INLINE size_t numberOfRows() const noexcept
   {
     return m_values.extent(0);
   }
 
-  PUGS_INLINE size_t numberOfColumns() const noexcept
+  [[nodiscard]] PUGS_INLINE size_t numberOfColumns() const noexcept
   {
     return m_values.extent(1);
   }
 
-  PUGS_INLINE
-  Array<DataType> operator[](index_type i) const
+  [[nodiscard]] PUGS_INLINE Table<DataType>::UnsafeRowView operator[](index_type i) const
   {
-    Assert(i < this->numberOfRows());
-    return encapsulate(Kokkos::View<DataType*>(m_values, i, Kokkos::ALL));
+    return UnsafeRowView(*this, i);
   }
 
-  friend PUGS_INLINE Table<std::remove_const_t<DataType>> copy(const Table<DataType>& source)
+  [[nodiscard]] friend PUGS_INLINE Table<std::remove_const_t<DataType>> copy(const Table<DataType>& source)
   {
     Table<std::remove_const_t<DataType>> image(source.numberOfRows(), source.numberOfColumns());
     Kokkos::deep_copy(image.m_values, source.m_values);
@@ -53,8 +207,8 @@ class [[nodiscard]] Table
   friend PUGS_INLINE void copy_to(const Table<DataType>& source,
                                   const Table<std::remove_const_t<DataType>>& destination)
   {
-    Assert(source.numberOfRows() == destination.numberOfRows());
-    Assert(source.numberOfColumns() == destination.numberOfColumns());
+    Assert(source.numberOfRows() == destination.numberOfRows(), "incompatible number of rows");
+    Assert(source.numberOfColumns() == destination.numberOfColumns(), "incompatible number of columns");
     Kokkos::deep_copy(destination.m_values, source.m_values);
   }
 
@@ -62,16 +216,15 @@ class [[nodiscard]] Table
   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);
+  friend PUGS_INLINE typename Table<DataType2>::UnsafeTableView
+  subTableView(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)
+  [[nodiscard]] PUGS_INLINE DataType& operator()(index_type i, index_type j) const noexcept(NO_ASSERT)
   {
-    Assert(i < this->numberOfRows());
-    Assert(j < this->numberOfColumns());
+    Assert(i < this->numberOfRows(), "invalid row index");
+    Assert(j < this->numberOfColumns(), "invalid column index");
     return m_values(i, j);
   }
 
@@ -160,7 +313,7 @@ class [[nodiscard]] Table
 };
 
 template <typename DataType, typename... RT>
-PUGS_INLINE Table<DataType>
+[[nodiscard]] PUGS_INLINE Table<DataType>
 encapsulate(const Kokkos::View<DataType**, RT...>& values)
 {
   Table<DataType> table;
@@ -169,19 +322,14 @@ encapsulate(const Kokkos::View<DataType**, RT...>& values)
 }
 
 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)
+[[nodiscard]] PUGS_INLINE typename Table<DataType>::UnsafeTableView
+subTableView(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.numberOfRows());
-  Assert(row_begin + row_size <= table.numberOfRows());
-  Assert(column_begin < table.numberOfColumns());
-  Assert(column_begin + column_size <= table.numberOfColumns());
-  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)));
+  return typename Table<DataType>::UnsafeTableView(table, row_begin, row_size, column_begin, column_size);
 }
 
 #endif   // TABLE_HPP
diff --git a/tests/test_Array.cpp b/tests/test_Array.cpp
index fc5a76735fcb78c81ca7635751f39e85f8d24040..f7d121df0ab11e229bae06f505e8f12361fe00da 100644
--- a/tests/test_Array.cpp
+++ b/tests/test_Array.cpp
@@ -319,6 +319,51 @@ TEST_CASE("Array", "[utils]")
     }
   }
 
+  SECTION("checking for subArrayView")
+  {
+    Array<int> array{10};
+    for (size_t i = 0; i < array.size(); ++i) {
+      array[i] = 2 * i + 1;
+    }
+
+    auto view = subArrayView(array, 3, 6);
+
+    REQUIRE(view.size() == 6);
+
+    bool use_same_memory = true;
+    for (size_t i = 0; i < view.size(); ++i) {
+      use_same_memory &= (&(view[i]) == &(array[i + 3]));
+    }
+    REQUIRE(use_same_memory);
+
+    view.fill(-3);
+    for (size_t i = 0; i < view.size(); ++i) {
+      REQUIRE(view[i] == -3);
+    }
+
+    for (size_t i = 0; i < array.size(); ++i) {
+      int int_i = i;
+      if ((i >= 3) and (i < 9)) {
+        REQUIRE(array[i] == -3);
+      } else {
+        REQUIRE(array[i] == 2 * int_i + 1);
+      }
+    }
+
+    for (size_t i = 0; i < view.size(); ++i) {
+      view[i] = 3 - 2 * i;
+    }
+
+    for (size_t i = 0; i < array.size(); ++i) {
+      int int_i = i;
+      if ((i >= 3) and (i < 9)) {
+        REQUIRE(array[i] == 3 - 2 * (int_i - 3));
+      } else {
+        REQUIRE(array[i] == 2 * int_i + 1);
+      }
+    }
+  }
+
 #ifndef NDEBUG
 
   SECTION("output with signaling NaN")
@@ -340,13 +385,13 @@ TEST_CASE("Array", "[utils]")
 
   SECTION("checking for bounds violation")
   {
-    REQUIRE_THROWS_AS(a[10], AssertError);
+    REQUIRE_THROWS_WITH(a[10], "invalid index");
   }
 
   SECTION("invalid copy_to")
   {
     Array<int> b{2 * a.size()};
-    REQUIRE_THROWS_AS(copy_to(a, b), AssertError);
+    REQUIRE_THROWS_WITH(copy_to(a, b), "incompatible Array sizes");
   }
 
   SECTION("checking for nan initialization")
@@ -366,5 +411,21 @@ TEST_CASE("Array", "[utils]")
       REQUIRE(array[i] == std::numeric_limits<int>::max() / 2);
     }
   }
+
+  SECTION("invalid UnsafeArrayView")
+  {
+    Array<int> array(10);
+
+    REQUIRE_THROWS_WITH(subArrayView(array, 3, 8), "required view is not contained in the Array");
+    REQUIRE_THROWS_WITH(subArrayView(array, 9, 2), "required view is not contained in the Array");
+  }
+
+  SECTION("check for bounds in UnsafeArrayView")
+  {
+    Array<int> array(10);
+    auto view = subArrayView(array, 2, 8);
+
+    REQUIRE_THROWS_WITH(view[8], "invalid index");
+  }
 #endif   // NDEBUG
 }
diff --git a/tests/test_ItemArray.cpp b/tests/test_ItemArray.cpp
index b98c12ba2528f2933a6ec9f7475fa559fa7d055b..8e60b44c1fddc10b6c80d6cd3e11f0bae0c45a0d 100644
--- a/tests/test_ItemArray.cpp
+++ b/tests/test_ItemArray.cpp
@@ -160,7 +160,7 @@ TEST_CASE("ItemArray", "[mesh]")
         auto is_same = [](const CellArray<size_t>& cell_array, const Table<size_t>& table) {
           bool is_same = true;
           for (CellId cell_id = 0; cell_id < cell_array.numberOfItems(); ++cell_id) {
-            Array sub_array = cell_array[cell_id];
+            auto 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));
             }
@@ -179,7 +179,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) {
-        Array sub_array = cell_array[cell_id];
+        auto sub_array = cell_array[cell_id];
         for (size_t i = 0; i < sub_array.size(); ++i) {
           is_same &= (sub_array[i] == value);
         }
diff --git a/tests/test_ItemArrayUtils.cpp b/tests/test_ItemArrayUtils.cpp
index 99958839ce3ce945945d0156ff3243ae1ebd5b85..c08167d7ecef9bbb0592838c075624464c5c4a4a 100644
--- a/tests/test_ItemArrayUtils.cpp
+++ b/tests/test_ItemArrayUtils.cpp
@@ -33,7 +33,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
             auto node_number = connectivity.nodeNumber();
 
             for (NodeId i_node = 0; i_node < mesh_1d->numberOfNodes(); ++i_node) {
-              Array array         = weak_node_array[i_node];
+              auto 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;
@@ -51,7 +51,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) {
-                Array array         = node_array[i_node];
+                auto 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]) {
@@ -79,7 +79,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
               bool is_synchronized = true;
               for (NodeId i_node = 0; i_node < mesh_1d->numberOfNodes(); ++i_node) {
-                Array array         = node_array[i_node];
+                auto 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) {
@@ -98,7 +98,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
             auto edge_number = connectivity.edgeNumber();
 
             for (EdgeId i_edge = 0; i_edge < mesh_1d->numberOfEdges(); ++i_edge) {
-              Array array         = weak_edge_array[i_edge];
+              auto 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;
@@ -116,7 +116,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) {
-                Array array         = edge_array[i_edge];
+                auto 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]) {
@@ -144,7 +144,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
               bool is_synchronized = true;
               for (EdgeId i_edge = 0; i_edge < mesh_1d->numberOfEdges(); ++i_edge) {
-                Array array         = edge_array[i_edge];
+                auto 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) {
@@ -163,7 +163,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
             auto face_number = connectivity.faceNumber();
 
             for (FaceId i_face = 0; i_face < mesh_1d->numberOfFaces(); ++i_face) {
-              Array array         = weak_face_array[i_face];
+              auto 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;
@@ -181,7 +181,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) {
-                Array array         = face_array[i_face];
+                auto 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]) {
@@ -209,7 +209,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
               bool is_synchronized = true;
               for (FaceId i_face = 0; i_face < mesh_1d->numberOfFaces(); ++i_face) {
-                Array array         = face_array[i_face];
+                auto 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) {
@@ -228,7 +228,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
             auto cell_number = connectivity.cellNumber();
 
             for (CellId i_cell = 0; i_cell < mesh_1d->numberOfCells(); ++i_cell) {
-              Array array         = weak_cell_array[i_cell];
+              auto 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;
@@ -246,7 +246,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) {
-                Array array         = cell_array[i_cell];
+                auto 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]) {
@@ -274,7 +274,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
               bool is_synchronized = true;
               for (CellId i_cell = 0; i_cell < mesh_1d->numberOfCells(); ++i_cell) {
-                Array array         = cell_array[i_cell];
+                auto 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) {
@@ -307,7 +307,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
             auto node_number = connectivity.nodeNumber();
 
             for (NodeId i_node = 0; i_node < mesh_2d->numberOfNodes(); ++i_node) {
-              Array array         = weak_node_array[i_node];
+              auto 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;
@@ -325,7 +325,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) {
-                Array array         = node_array[i_node];
+                auto 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]) {
@@ -353,7 +353,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
               bool is_synchronized = true;
               for (NodeId i_node = 0; i_node < mesh_2d->numberOfNodes(); ++i_node) {
-                Array array         = node_array[i_node];
+                auto 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) {
@@ -372,7 +372,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
             auto edge_number = connectivity.edgeNumber();
 
             for (EdgeId i_edge = 0; i_edge < mesh_2d->numberOfEdges(); ++i_edge) {
-              Array array         = weak_edge_array[i_edge];
+              auto 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;
@@ -390,7 +390,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) {
-                Array array         = edge_array[i_edge];
+                auto 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]) {
@@ -418,7 +418,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
               bool is_synchronized = true;
               for (EdgeId i_edge = 0; i_edge < mesh_2d->numberOfEdges(); ++i_edge) {
-                Array array         = edge_array[i_edge];
+                auto 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) {
@@ -437,7 +437,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
             auto face_number = connectivity.faceNumber();
 
             for (FaceId i_face = 0; i_face < mesh_2d->numberOfFaces(); ++i_face) {
-              Array array         = weak_face_array[i_face];
+              auto 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;
@@ -455,7 +455,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) {
-                Array array         = face_array[i_face];
+                auto 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]) {
@@ -483,7 +483,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
               bool is_synchronized = true;
               for (FaceId i_face = 0; i_face < mesh_2d->numberOfFaces(); ++i_face) {
-                Array array         = face_array[i_face];
+                auto 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) {
@@ -502,7 +502,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
             auto cell_number = connectivity.cellNumber();
 
             for (CellId i_cell = 0; i_cell < mesh_2d->numberOfCells(); ++i_cell) {
-              Array array         = weak_cell_array[i_cell];
+              auto 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;
@@ -520,7 +520,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) {
-                Array array         = cell_array[i_cell];
+                auto 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]) {
@@ -548,7 +548,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
               bool is_synchronized = true;
               for (CellId i_cell = 0; i_cell < mesh_2d->numberOfCells(); ++i_cell) {
-                Array array         = cell_array[i_cell];
+                auto 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) {
@@ -581,7 +581,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
             auto node_number = connectivity.nodeNumber();
 
             for (NodeId i_node = 0; i_node < mesh_3d->numberOfNodes(); ++i_node) {
-              Array array         = weak_node_array[i_node];
+              auto 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;
@@ -599,7 +599,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) {
-                Array array         = node_array[i_node];
+                auto 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]) {
@@ -627,7 +627,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
               bool is_synchronized = true;
               for (NodeId i_node = 0; i_node < mesh_3d->numberOfNodes(); ++i_node) {
-                Array array         = node_array[i_node];
+                auto 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) {
@@ -646,7 +646,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
             auto edge_number = connectivity.edgeNumber();
 
             for (EdgeId i_edge = 0; i_edge < mesh_3d->numberOfEdges(); ++i_edge) {
-              Array array         = weak_edge_array[i_edge];
+              auto 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;
@@ -664,7 +664,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) {
-                Array array         = edge_array[i_edge];
+                auto 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]) {
@@ -692,7 +692,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
               bool is_synchronized = true;
               for (EdgeId i_edge = 0; i_edge < mesh_3d->numberOfEdges(); ++i_edge) {
-                Array array         = edge_array[i_edge];
+                auto 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) {
@@ -711,7 +711,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
             auto face_number = connectivity.faceNumber();
 
             for (FaceId i_face = 0; i_face < mesh_3d->numberOfFaces(); ++i_face) {
-              Array array         = weak_face_array[i_face];
+              auto 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;
@@ -729,7 +729,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) {
-                Array array         = face_array[i_face];
+                auto 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]) {
@@ -757,7 +757,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
               bool is_synchronized = true;
               for (FaceId i_face = 0; i_face < mesh_3d->numberOfFaces(); ++i_face) {
-                Array array         = face_array[i_face];
+                auto 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) {
@@ -776,7 +776,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
             auto cell_number = connectivity.cellNumber();
 
             for (CellId i_cell = 0; i_cell < mesh_3d->numberOfCells(); ++i_cell) {
-              Array array         = weak_cell_array[i_cell];
+              auto 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;
@@ -794,7 +794,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) {
-                Array array         = cell_array[i_cell];
+                auto 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]) {
@@ -822,7 +822,7 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
 
               bool is_synchronized = true;
               for (CellId i_cell = 0; i_cell < mesh_3d->numberOfCells(); ++i_cell) {
-                Array array         = cell_array[i_cell];
+                auto 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_Table.cpp b/tests/test_Table.cpp
index a89344b75262a2eb3322589b613d6add8ae5764e..5ce90fe9ce460041d26ebcfbc88ff891da0ffdde 100644
--- a/tests/test_Table.cpp
+++ b/tests/test_Table.cpp
@@ -196,6 +196,173 @@ TEST_CASE("Table", "[utils]")
     REQUIRE(table_ost.str() == ref_ost.str());
   }
 
+  SECTION("UnsafeTableView")
+  {
+    Table<int> table(3, 4);
+    for (size_t i = 0; i < table.numberOfRows(); ++i) {
+      for (size_t j = 0; j < table.numberOfColumns(); ++j) {
+        table(i, j) = 2 * i + 3 * j;
+      }
+    }
+
+    auto sub_table_view = subTableView(table, 1, 2, 1, 3);
+    REQUIRE(sub_table_view.numberOfRows() == 2);
+    REQUIRE(sub_table_view.numberOfColumns() == 3);
+
+    {
+      bool is_same_memory = true;
+      for (size_t i = 0; i < sub_table_view.numberOfRows(); ++i) {
+        for (size_t j = 0; j < sub_table_view.numberOfColumns(); ++j) {
+          is_same_memory &= (&(sub_table_view(i, j)) == &(table(i + 1, j + 1)));
+        }
+      }
+      REQUIRE(is_same_memory);
+    }
+
+    for (size_t i = 0; i < sub_table_view.numberOfRows(); ++i) {
+      int ip1 = i + 1;
+      for (size_t j = 0; j < sub_table_view.numberOfColumns(); ++j) {
+        int jp1 = j + 1;
+        REQUIRE(sub_table_view(i, j) == 2 * ip1 + 3 * jp1);
+      }
+    }
+
+    sub_table_view.fill(-2);
+    for (size_t i = 0; i < table.numberOfRows(); ++i) {
+      int int_i = i;
+      for (size_t j = 0; j < table.numberOfColumns(); ++j) {
+        int int_j = j;
+        if ((i > 0) and (j > 0)) {
+          REQUIRE(table(i, j) == -2);
+        } else {
+          REQUIRE(table(i, j) == 2 * int_i + 3 * int_j);
+        }
+      }
+    }
+
+    for (size_t i = 0; i < sub_table_view.numberOfRows(); ++i) {
+      for (size_t j = 0; j < sub_table_view.numberOfColumns(); ++j) {
+        sub_table_view(i, j) = -3 * i + 2 * j;
+      }
+    }
+
+    for (size_t i = 0; i < table.numberOfRows(); ++i) {
+      int int_i = i;
+      for (size_t j = 0; j < table.numberOfColumns(); ++j) {
+        int int_j = j;
+        if ((i > 0) and (j > 0)) {
+          REQUIRE(table(i, j) == -3 * (int_i - 1) + 2 * (int_j - 1));
+        } else {
+          REQUIRE(table(i, j) == 2 * int_i + 3 * int_j);
+        }
+      }
+    }
+  }
+
+  SECTION("UnsafeTableView::RowView")
+  {
+    Table<int> table(3, 4);
+    for (size_t i = 0; i < table.numberOfRows(); ++i) {
+      for (size_t j = 0; j < table.numberOfColumns(); ++j) {
+        table(i, j) = 2 * i + 3 * j;
+      }
+    }
+
+    auto sub_table_view = subTableView(table, 1, 2, 1, 3);
+
+    auto sub_table_row_view = sub_table_view[1];
+
+    REQUIRE(sub_table_row_view.size() == 3);
+
+    {
+      bool is_same_memory = true;
+      for (size_t i = 0; i < sub_table_row_view.size(); ++i) {
+        is_same_memory &= (&(sub_table_row_view[i]) == &(table(2, i + 1)));
+      }
+      REQUIRE(is_same_memory);
+    }
+
+    sub_table_row_view.fill(-5);
+    for (size_t i = 0; i < table.numberOfRows(); ++i) {
+      int int_i = i;
+      for (size_t j = 0; j < table.numberOfColumns(); ++j) {
+        int int_j = j;
+        if ((i == 2) and (j > 0)) {
+          REQUIRE(table(i, j) == -5);
+        } else {
+          REQUIRE(table(i, j) == 2 * int_i + 3 * int_j);
+        }
+      }
+    }
+
+    for (size_t i = 0; i < sub_table_row_view.size(); ++i) {
+      sub_table_row_view[i] = -3 * i + 2;
+    }
+
+    for (size_t i = 0; i < table.numberOfRows(); ++i) {
+      int int_i = i;
+      for (size_t j = 0; j < table.numberOfColumns(); ++j) {
+        int int_j = j;
+        if ((i == 2) and (j > 0)) {
+          REQUIRE(table(i, j) == -3 * (int_j - 1) + 2);
+        } else {
+          REQUIRE(table(i, j) == 2 * int_i + 3 * int_j);
+        }
+      }
+    }
+  }
+
+  SECTION("UnsafeRowView")
+  {
+    Table<int> table(3, 4);
+    for (size_t i = 0; i < table.numberOfRows(); ++i) {
+      for (size_t j = 0; j < table.numberOfColumns(); ++j) {
+        table(i, j) = 2 * i + 3 * j;
+      }
+    }
+
+    auto table_row_view = table[1];
+
+    REQUIRE(table_row_view.size() == table.numberOfColumns());
+
+    {
+      bool is_same_memory = true;
+      for (size_t i = 0; i < table_row_view.size(); ++i) {
+        is_same_memory &= (&(table_row_view[i]) == &(table(1, i)));
+      }
+      REQUIRE(is_same_memory);
+    }
+
+    table_row_view.fill(-5);
+    for (size_t i = 0; i < table.numberOfRows(); ++i) {
+      int int_i = i;
+      for (size_t j = 0; j < table.numberOfColumns(); ++j) {
+        int int_j = j;
+        if (i == 1) {
+          REQUIRE(table(i, j) == -5);
+        } else {
+          REQUIRE(table(i, j) == 2 * int_i + 3 * int_j);
+        }
+      }
+    }
+
+    for (size_t i = 0; i < table_row_view.size(); ++i) {
+      table_row_view[i] = -3 * i + 2;
+    }
+
+    for (size_t i = 0; i < table.numberOfRows(); ++i) {
+      int int_i = i;
+      for (size_t j = 0; j < table.numberOfColumns(); ++j) {
+        int int_j = j;
+        if (i == 1) {
+          REQUIRE(table(i, j) == -3 * int_j + 2);
+        } else {
+          REQUIRE(table(i, j) == 2 * int_i + 3 * int_j);
+        }
+      }
+    }
+  }
+
 #ifndef NDEBUG
   SECTION("output with signaling NaN")
   {
@@ -216,8 +383,8 @@ TEST_CASE("Table", "[utils]")
 
   SECTION("checking for bounds violation")
   {
-    REQUIRE_THROWS_AS(a(4, 0), AssertError);
-    REQUIRE_THROWS_AS(a(0, 3), AssertError);
+    REQUIRE_THROWS_WITH(a(4, 0), "invalid row index");
+    REQUIRE_THROWS_WITH(a(0, 3), "invalid column index");
   }
 
   SECTION("invalid copy_to")
@@ -225,13 +392,13 @@ TEST_CASE("Table", "[utils]")
     SECTION("wrong row number")
     {
       Table<int> b{2 * a.numberOfRows(), a.numberOfColumns()};
-      REQUIRE_THROWS_AS(copy_to(a, b), AssertError);
+      REQUIRE_THROWS_WITH(copy_to(a, b), "incompatible number of rows");
     }
 
     SECTION("wrong column number")
     {
       Table<int> c{a.numberOfRows(), 2 * a.numberOfColumns()};
-      REQUIRE_THROWS_AS(copy_to(a, c), AssertError);
+      REQUIRE_THROWS_WITH(copy_to(a, c), "incompatible number of columns");
     }
   }
 
@@ -256,5 +423,59 @@ TEST_CASE("Table", "[utils]")
       }
     }
   }
+
+  SECTION("invalid UnsafeTableView")
+  {
+    Table<int> table(3, 4);
+
+    REQUIRE_THROWS_WITH(subTableView(table, 2, 2, 2, 2), "required view has rows not contained in the Table");
+    REQUIRE_THROWS_WITH(subTableView(table, 4, 2, 2, 2), "required view has rows not contained in the Table");
+
+    REQUIRE_THROWS_WITH(subTableView(table, 1, 2, 2, 4), "required view has columns not contained in the Table");
+    REQUIRE_THROWS_WITH(subTableView(table, 1, 2, 4, 2), "required view has columns not contained in the Table");
+  }
+
+  SECTION("check for bounds in UnsafeTableView")
+  {
+    Table<int> table(3, 4);
+    auto view = subTableView(table, 1, 2, 1, 3);
+
+    REQUIRE_THROWS_WITH(view(2, 1), "invalid row index");
+    REQUIRE_THROWS_WITH(view(1, 3), "invalid column index");
+  }
+
+  SECTION("invalid UnsafeTableView")
+  {
+    Table<int> table(3, 4);
+    auto sub_table_view = subTableView(table, 1, 2, 1, 3);
+
+    REQUIRE_THROWS_WITH(sub_table_view[2], "invalid index");
+  }
+
+  SECTION("check for bounds in UnsafeTableView::RowView")
+  {
+    Table<int> table(3, 4);
+    auto sub_table_view = subTableView(table, 1, 2, 1, 3);
+
+    auto sub_table_row_view = sub_table_view[1];
+
+    REQUIRE_THROWS_WITH(sub_table_row_view[3], "invalid index");
+  }
+
+  SECTION("invalid UnsafeRowView")
+  {
+    Table<int> table(3, 4);
+
+    REQUIRE_THROWS_WITH(table[3], "required row view is not contained in the Table");
+  }
+
+  SECTION("check for bounds in UnsafeRowView")
+  {
+    Table<int> table(3, 4);
+    auto table_row_view = table[1];
+
+    REQUIRE_THROWS_WITH(table_row_view[4], "invalid index");
+  }
+
 #endif   // NDEBUG
 }