diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index 0a02886dd673680a7c800d549b5ee3fc9aacd07a..a819af495a1dbb358a2705231e3ff30f2bcbed02 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -7,8 +7,7 @@
 #include <map>
 
 template <size_t Dimension>
-Connectivity<Dimension>::Connectivity()
-{}
+Connectivity<Dimension>::Connectivity() = default;
 
 template <size_t Dimension>
 void
@@ -22,6 +21,22 @@ Connectivity<Dimension>::_buildFrom(const ConnectivityDescriptor& descriptor)
     Assert(descriptor.face_owner_vector.size() == descriptor.face_number_vector.size());
   }
 
+  m_number_of_cells = descriptor.cell_number_vector.size();
+  m_number_of_nodes = descriptor.node_number_vector.size();
+
+  if constexpr (Dimension == 1) {
+    m_number_of_edges = m_number_of_nodes;
+    m_number_of_faces = m_number_of_nodes;
+  } else {
+    m_number_of_faces = descriptor.face_number_vector.size();
+    if constexpr (Dimension == 2) {
+      m_number_of_edges = m_number_of_faces;
+    } else {
+      static_assert(Dimension == 3, "unexpected dimension");
+      m_number_of_edges = descriptor.edge_number_vector.size();
+    }
+  }
+
   auto& cell_to_node_matrix = m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::node)];
   cell_to_node_matrix       = descriptor.cell_to_node_vector;
 
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index f1e8c9f2dcc9e883987fde4fc69baed16971dbd0..86b6606c7f24c2271569f07cbb5fe9d10a8b8fe3 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -47,6 +47,11 @@ class Connectivity final : public IConnectivity
   }
 
  private:
+  size_t m_number_of_cells;
+  size_t m_number_of_faces;
+  size_t m_number_of_edges;
+  size_t m_number_of_nodes;
+
   ConnectivityMatrix m_item_to_item_matrix[Dimension + 1][Dimension + 1];
   WeakCellValue<const CellType> m_cell_type;
 
@@ -75,21 +80,21 @@ class Connectivity final : public IConnectivity
   WeakFaceValuePerCell<const bool> m_cell_face_is_reversed;
   WeakEdgeValuePerFace<const bool> m_face_edge_is_reversed;
 
-  WeakNodeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_nodes;
-  WeakFaceValuePerCell<const unsigned short> m_cell_local_numbers_in_their_faces;
-  WeakEdgeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_edges;
+  WeakFaceValuePerCell<const uint16_t> m_cell_local_numbers_in_their_faces;
+  WeakEdgeValuePerCell<const uint16_t> m_cell_local_numbers_in_their_edges;
+  WeakNodeValuePerCell<const uint16_t> m_cell_local_numbers_in_their_nodes;
 
-  WeakCellValuePerFace<const unsigned short> m_face_local_numbers_in_their_cells;
-  WeakEdgeValuePerFace<const unsigned short> m_face_local_numbers_in_their_edges;
-  WeakNodeValuePerFace<const unsigned short> m_face_local_numbers_in_their_nodes;
+  WeakCellValuePerFace<const uint16_t> m_face_local_numbers_in_their_cells;
+  WeakEdgeValuePerFace<const uint16_t> m_face_local_numbers_in_their_edges;
+  WeakNodeValuePerFace<const uint16_t> m_face_local_numbers_in_their_nodes;
 
-  WeakCellValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_cells;
-  WeakFaceValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_faces;
-  WeakNodeValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_nodes;
+  WeakCellValuePerEdge<const uint16_t> m_edge_local_numbers_in_their_cells;
+  WeakFaceValuePerEdge<const uint16_t> m_edge_local_numbers_in_their_faces;
+  WeakNodeValuePerEdge<const uint16_t> m_edge_local_numbers_in_their_nodes;
 
-  WeakCellValuePerNode<const unsigned short> m_node_local_numbers_in_their_cells;
-  WeakFaceValuePerNode<const unsigned short> m_node_local_numbers_in_their_faces;
-  WeakEdgeValuePerNode<const unsigned short> m_node_local_numbers_in_their_edges;
+  WeakCellValuePerNode<const uint16_t> m_node_local_numbers_in_their_cells;
+  WeakFaceValuePerNode<const uint16_t> m_node_local_numbers_in_their_faces;
+  WeakEdgeValuePerNode<const uint16_t> m_node_local_numbers_in_their_edges;
 
   ConnectivityComputer m_connectivity_computer;
 
@@ -98,31 +103,18 @@ class Connectivity final : public IConnectivity
   std::vector<RefEdgeList> m_ref_edge_list_vector;
   std::vector<RefNodeList> m_ref_node_list_vector;
 
-  WeakCellValue<const double> m_inv_cell_nb_nodes;
-
   void _computeCellFaceAndFaceNodeConnectivities();
 
   void _buildIsBoundaryFace() const;
   void _buildIsBoundaryEdge() const;
   void _buildIsBoundaryNode() const;
 
-  template <typename SubItemValuePerItemType>
-  PUGS_INLINE const SubItemValuePerItemType&
-  _lazzyBuildItemNumberInTheirChild(const SubItemValuePerItemType& sub_item_value_per_item) const
-  {
-    using ReversedItemOfItem = typename SubItemValuePerItemType::ItemOfItemType::Reversed;
-    if (not sub_item_value_per_item.isBuilt()) {
-      const_cast<SubItemValuePerItemType&>(sub_item_value_per_item) =
-        m_connectivity_computer.computeLocalItemNumberInChildItem<ReversedItemOfItem>(*this);
-    }
-    return sub_item_value_per_item;
-  }
-
   friend class ConnectivityComputer;
 
+ public:
   PUGS_INLINE
   const ConnectivityMatrix&
-  _getMatrix(const ItemType& item_type_0, const ItemType& item_type_1) const
+  getMatrix(const ItemType& item_type_0, const ItemType& item_type_1) const final
   {
     const ConnectivityMatrix& connectivity_matrix = m_item_to_item_matrix[itemTId(item_type_0)][itemTId(item_type_1)];
     if (not connectivity_matrix.isBuilt()) {
@@ -132,7 +124,6 @@ class Connectivity final : public IConnectivity
     return connectivity_matrix;
   }
 
- public:
   PUGS_INLINE
   CellValue<const CellType>
   cellType() const
@@ -333,7 +324,7 @@ class Connectivity final : public IConnectivity
   PUGS_INLINE ItemToItemMatrix<source_item_type, target_item_type>
   getItemToItemMatrix() const
   {
-    return ItemToItemMatrix<source_item_type, target_item_type>(_getMatrix(source_item_type, target_item_type));
+    return ItemToItemMatrix<source_item_type, target_item_type>(getMatrix(source_item_type, target_item_type));
   }
 
   PUGS_INLINE
@@ -443,109 +434,129 @@ class Connectivity final : public IConnectivity
   }
 
   PUGS_INLINE
-  const auto&
-  cellLocalNumbersInTheirNodes() const
+  FaceValuePerCell<const uint16_t>
+  cellLocalNumbersInTheirFaces() const
   {
-    return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_nodes);
+    if (not m_cell_local_numbers_in_their_faces.isBuilt()) {
+      m_connectivity_computer.computeLocalItemNumberInChildItem<FaceOfCell>(*this);
+    }
+    return m_cell_local_numbers_in_their_faces;
   }
 
   PUGS_INLINE
-  const auto&
+  EdgeValuePerCell<const uint16_t>
   cellLocalNumbersInTheirEdges() const
   {
-    if constexpr (Dimension > 2) {
-      return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_edges);
-    } else {
-      return cellLocalNumbersInTheirFaces();
+    if (not m_cell_local_numbers_in_their_edges.isBuilt()) {
+      m_connectivity_computer.computeLocalItemNumberInChildItem<EdgeOfCell>(*this);
     }
+    return m_cell_local_numbers_in_their_edges;
   }
 
   PUGS_INLINE
-  const auto&
-  cellLocalNumbersInTheirFaces() const
+  NodeValuePerCell<const uint16_t>
+  cellLocalNumbersInTheirNodes() const
   {
-    if constexpr (Dimension > 1) {
-      return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_faces);
-    } else {
-      return cellLocalNumbersInTheirNodes();
+    if (not m_cell_local_numbers_in_their_nodes.isBuilt()) {
+      m_connectivity_computer.computeLocalItemNumberInChildItem<NodeOfCell>(*this);
     }
+    return m_cell_local_numbers_in_their_nodes;
   }
 
   PUGS_INLINE
-  const auto&
+  CellValuePerFace<const uint16_t>
   faceLocalNumbersInTheirCells() const
   {
-    if constexpr (Dimension > 1) {
-      return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_cells);
-    } else {
-      return nodeLocalNumbersInTheirCells();
+    if (not m_face_local_numbers_in_their_cells.isBuilt()) {
+      m_connectivity_computer.computeLocalItemNumberInChildItem<CellOfFace>(*this);
     }
+    return m_face_local_numbers_in_their_cells;
   }
 
   PUGS_INLINE
-  const auto&
+  EdgeValuePerFace<const uint16_t>
   faceLocalNumbersInTheirEdges() const
   {
     static_assert(Dimension > 2, "this function has no meaning in 1d or 2d");
-    return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_edges);
+    if (not m_face_local_numbers_in_their_edges.isBuilt()) {
+      m_connectivity_computer.computeLocalItemNumberInChildItem<EdgeOfFace>(*this);
+    }
+    return m_face_local_numbers_in_their_edges;
   }
 
   PUGS_INLINE
-  const auto&
+  NodeValuePerFace<const uint16_t>
   faceLocalNumbersInTheirNodes() const
   {
     static_assert(Dimension > 1, "this function has no meaning in 1d");
-    return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_nodes);
+    if (not m_face_local_numbers_in_their_nodes.isBuilt()) {
+      m_connectivity_computer.computeLocalItemNumberInChildItem<NodeOfFace>(*this);
+    }
+    return m_face_local_numbers_in_their_nodes;
   }
 
   PUGS_INLINE
-  const auto&
+  CellValuePerEdge<const uint16_t>
   edgeLocalNumbersInTheirCells() const
   {
-    if constexpr (Dimension > 2) {
-      return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_cells);
-    } else {
-      return faceLocalNumbersInTheirCells();
+    if (not m_edge_local_numbers_in_their_cells.isBuilt()) {
+      m_connectivity_computer.computeLocalItemNumberInChildItem<CellOfEdge>(*this);
     }
+    return m_edge_local_numbers_in_their_cells;
   }
 
   PUGS_INLINE
-  const auto&
+  FaceValuePerEdge<const uint16_t>
   edgeLocalNumbersInTheirFaces() const
   {
     static_assert(Dimension > 2, "this function has no meaning in 1d or 2d");
-    return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_faces);
+    if (not m_edge_local_numbers_in_their_faces.isBuilt()) {
+      m_connectivity_computer.computeLocalItemNumberInChildItem<FaceOfEdge>(*this);
+    }
+    return m_edge_local_numbers_in_their_faces;
   }
 
   PUGS_INLINE
-  const auto&
+  NodeValuePerEdge<const uint16_t>
   edgeLocalNumbersInTheirNodes() const
   {
-    static_assert(Dimension > 2, "this function has no meaning in 1d or 2d");
-    return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_nodes);
+    static_assert(Dimension > 1, "this function has no meaning in 1d");
+    if (not m_edge_local_numbers_in_their_nodes.isBuilt()) {
+      m_connectivity_computer.computeLocalItemNumberInChildItem<NodeOfEdge>(*this);
+    }
+    return m_edge_local_numbers_in_their_nodes;
   }
 
   PUGS_INLINE
-  const auto&
+  CellValuePerNode<const uint16_t>
   nodeLocalNumbersInTheirCells() const
   {
-    return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_cells);
+    if (not m_node_local_numbers_in_their_cells.isBuilt()) {
+      m_connectivity_computer.computeLocalItemNumberInChildItem<CellOfNode>(*this);
+    }
+    return m_node_local_numbers_in_their_cells;
   }
 
   PUGS_INLINE
-  const auto&
-  nodeLocalNumbersInTheirEdges() const
+  FaceValuePerNode<const uint16_t>
+  nodeLocalNumbersInTheirFaces() const
   {
-    static_assert(Dimension > 2, "this function has no meaning in 1d or 2d");
-    return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_edges);
+    static_assert(Dimension > 1, "this function has no meaning in 1d");
+    if (not m_node_local_numbers_in_their_faces.isBuilt()) {
+      m_connectivity_computer.computeLocalItemNumberInChildItem<FaceOfNode>(*this);
+    }
+    return m_node_local_numbers_in_their_faces;
   }
 
   PUGS_INLINE
-  const auto&
-  nodeLocalNumbersInTheirFaces() const
+  EdgeValuePerNode<const uint16_t>
+  nodeLocalNumbersInTheirEdges() const
   {
-    static_assert(Dimension > 1, "this function has no meaning in 1d");
-    return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_faces);
+    static_assert(Dimension > 1, "this function has no meaning in 1d or 2d");
+    if (not m_node_local_numbers_in_their_edges.isBuilt()) {
+      m_connectivity_computer.computeLocalItemNumberInChildItem<EdgeOfNode>(*this);
+    }
+    return m_node_local_numbers_in_their_edges;
   }
 
   template <ItemType item_type>
@@ -639,38 +650,28 @@ class Connectivity final : public IConnectivity
   size_t
   numberOfNodes() const final
   {
-    const auto& node_to_cell_matrix = this->_getMatrix(ItemType::node, ItemType::cell);
-    return node_to_cell_matrix.numRows();
+    return m_number_of_nodes;
   }
 
   PUGS_INLINE
   size_t
   numberOfEdges() const final
   {
-    if constexpr (Dimension == 1) {
-      return this->numberOfNodes();
-    } else if constexpr (Dimension == 2) {
-      return this->numberOfFaces();
-    } else {
-      const auto& edge_to_node_matrix = this->_getMatrix(ItemType::edge, ItemType::node);
-      return edge_to_node_matrix.numRows();
-    }
+    return m_number_of_edges;
   }
 
   PUGS_INLINE
   size_t
   numberOfFaces() const final
   {
-    const auto& face_to_node_matrix = this->_getMatrix(ItemType::face, ItemType::cell);
-    return face_to_node_matrix.numRows();
+    return m_number_of_faces;
   }
 
   PUGS_INLINE
   size_t
   numberOfCells() const final
   {
-    const auto& cell_to_node_matrix = this->_getMatrix(ItemType::cell, ItemType::node);
-    return cell_to_node_matrix.numRows();
+    return m_number_of_cells;
   }
 
   template <ItemType item_type>
@@ -690,24 +691,6 @@ class Connectivity final : public IConnectivity
     }
   }
 
-  CellValue<const double>
-  invCellNbNodes() const
-  {
-    if (not m_inv_cell_nb_nodes.isBuilt()) {
-      const auto& cell_to_node_matrix = m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::node)];
-
-      WeakCellValue<double> inv_cell_nb_nodes(*this);
-      parallel_for(
-        this->numberOfCells(), PUGS_LAMBDA(CellId j) {
-          const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
-          inv_cell_nb_nodes[j]   = 1. / cell_nodes.length;
-        });
-      const_cast<WeakCellValue<const double>&>(m_inv_cell_nb_nodes) = inv_cell_nb_nodes;
-    }
-
-    return m_inv_cell_nb_nodes;
-  }
-
   Connectivity(const Connectivity&) = delete;
 
  private:
diff --git a/src/mesh/ConnectivityComputer.cpp b/src/mesh/ConnectivityComputer.cpp
index 46deb5b95166d42bc3bb5ade5d3a02683157abfa..aa8ca779cfb965ffb3e85e4feb1b4b9feabb79ca 100644
--- a/src/mesh/ConnectivityComputer.cpp
+++ b/src/mesh/ConnectivityComputer.cpp
@@ -4,8 +4,6 @@
 #include <utils/Exceptions.hpp>
 
 #include <iostream>
-#include <map>
-#include <sstream>
 
 template <typename ConnectivityType>
 PUGS_INLINE ConnectivityMatrix
@@ -15,9 +13,7 @@ ConnectivityComputer::computeConnectivityMatrix(const ConnectivityType& connecti
 {
   ConnectivityMatrix item_to_child_item_matrix;
   if (connectivity.isConnectivityMatrixBuilt(child_item_type, item_type)) {
-    const ConnectivityMatrix& child_to_item_matrix = connectivity._getMatrix(child_item_type, item_type);
-
-    std::cout << "computing connectivity " << itemName(item_type) << " -> " << itemName(child_item_type) << '\n';
+    const ConnectivityMatrix& child_to_item_matrix = connectivity.getMatrix(child_item_type, item_type);
 
     item_to_child_item_matrix = this->_computeInverse(child_to_item_matrix);
   } else {
@@ -44,127 +40,271 @@ template ConnectivityMatrix ConnectivityComputer::computeConnectivityMatrix(cons
 ConnectivityMatrix
 ConnectivityComputer::_computeInverse(const ConnectivityMatrix& item_to_child_matrix) const
 {
-  std::map<unsigned int, std::vector<unsigned int>> child_to_item_vector_map;
-  const size_t& number_of_items = item_to_child_matrix.numRows();
+  const size_t& number_of_rows = item_to_child_matrix.numberOfRows();
+
+  if ((item_to_child_matrix.values().size() > 0)) {
+    const size_t& number_of_columns = max(item_to_child_matrix.values());
+
+    Array<uint32_t> transposed_next_free_column_index(number_of_columns + 1);
+    transposed_next_free_column_index.fill(0);
 
-  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);
+    Array<uint32_t> transposed_rows_map(transposed_next_free_column_index.size() + 1);
+    transposed_rows_map.fill(0);
+    for (size_t i = 0; i < number_of_rows; ++i) {
+      for (size_t j = item_to_child_matrix.rowsMap()[i]; j < item_to_child_matrix.rowsMap()[i + 1]; ++j) {
+        transposed_rows_map[item_to_child_matrix.values()[j] + 1]++;
+      }
     }
-  }
+    for (size_t i = 1; i < transposed_rows_map.size(); ++i) {
+      transposed_rows_map[i] += transposed_rows_map[i - 1];
+    }
+    Array<uint32_t> transposed_column_indices(transposed_rows_map[transposed_rows_map.size() - 1]);
 
-  {
-    size_t i = 0;
-    for (const auto& [child_item_id, item_vector] : child_to_item_vector_map) {
-      if (child_item_id != i) {
-        throw NotImplementedError("sparse item numerotation NIY\n");
+    for (size_t i = 0; i < number_of_rows; ++i) {
+      for (size_t j = item_to_child_matrix.rowsMap()[i]; j < item_to_child_matrix.rowsMap()[i + 1]; ++j) {
+        size_t i_column_index = item_to_child_matrix.values()[j];
+        uint32_t& shift       = transposed_next_free_column_index[i_column_index];
+
+        transposed_column_indices[transposed_rows_map[i_column_index] + shift] = i;
+        Assert(shift < transposed_rows_map[i_column_index + 1]);
+        shift++;
       }
-      ++i;
     }
-  }
 
-  std::vector<std::vector<unsigned int>> child_to_items_vector(child_to_item_vector_map.size());
-  for (const auto& [child_item_id, item_vector] : child_to_item_vector_map) {
-    child_to_items_vector[child_item_id] = item_vector;
-  }
+    return ConnectivityMatrix{transposed_rows_map, transposed_column_indices};
+  } else {
+    // empty connectivity
+    Array<uint32_t> transposed_rows_map{1};
+    transposed_rows_map[0] = 0;
+    Array<uint32_t> transposed_column_indices;
 
-  return ConnectivityMatrix(child_to_items_vector);
+    return ConnectivityMatrix{transposed_rows_map, transposed_column_indices};
+  }
 }
 
 template <typename ItemOfItem, typename ConnectivityType>
-WeakSubItemValuePerItem<const unsigned short, typename ItemOfItem::Reversed>
+void
 ConnectivityComputer::computeLocalItemNumberInChildItem(const ConnectivityType& connectivity) const
 {
-  using ReversedItemOfItem = typename ItemOfItem::Reversed;
-
-  constexpr ItemType item_type       = ReversedItemOfItem::item_type;
-  constexpr ItemType child_item_type = ReversedItemOfItem::sub_item_type;
+  constexpr ItemType item_type     = ItemOfItem::item_type;
+  constexpr ItemType sub_item_type = ItemOfItem::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, sub_item_type>();
+  ItemToItemMatrix child_item_to_items_matrix = connectivity.template getItemToItemMatrix<sub_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) {
-          item_number_in_child_item(r, J) = R;
-          break;
+  Array<unsigned short> number_in_child_array{item_to_child_items_matrix.numberOfValues()};
+  {
+    WeakSubItemValuePerItem<unsigned short, ItemOfItem> item_number_in_child_item(connectivity, number_in_child_array);
+    for (ItemIdT<item_type> item_id = 0; item_id < connectivity.template numberOf<item_type>(); ++item_id) {
+      const auto& item_to_child_items = item_to_child_items_matrix[item_id];
+      for (unsigned short i_item_child = 0; i_item_child < item_to_child_items_matrix[item_id].size(); ++i_item_child) {
+        ItemIdT j                       = item_to_child_items[i_item_child];
+        const auto& child_item_to_items = child_item_to_items_matrix[j];
+
+        for (unsigned int i_item_in_child = 0; i_item_in_child < child_item_to_items.size(); ++i_item_in_child) {
+          if (child_item_to_items[i_item_in_child] == item_id) {
+            item_number_in_child_item(item_id, i_item_child) = i_item_in_child;
+            break;
+          }
         }
       }
     }
   }
 
-  return item_number_in_child_item;
+  if constexpr (ConnectivityType::Dimension == 1) {
+    if constexpr (item_type == ItemType::cell) {
+      const_cast<WeakSubItemValuePerItem<const unsigned short, FaceOfCell>&>(
+        connectivity.m_cell_local_numbers_in_their_faces) =
+        WeakSubItemValuePerItem<unsigned short, FaceOfCell>(connectivity, number_in_child_array);
+      const_cast<WeakSubItemValuePerItem<const unsigned short, EdgeOfCell>&>(
+        connectivity.m_cell_local_numbers_in_their_edges) =
+        WeakSubItemValuePerItem<unsigned short, EdgeOfCell>(connectivity, number_in_child_array);
+      const_cast<WeakSubItemValuePerItem<const unsigned short, NodeOfCell>&>(
+        connectivity.m_cell_local_numbers_in_their_nodes) =
+        WeakSubItemValuePerItem<unsigned short, NodeOfCell>(connectivity, number_in_child_array);
+    } else if constexpr (sub_item_type == ItemType::cell) {
+      const_cast<WeakSubItemValuePerItem<const unsigned short, CellOfFace>&>(
+        connectivity.m_face_local_numbers_in_their_cells) =
+        WeakSubItemValuePerItem<unsigned short, CellOfFace>(connectivity, number_in_child_array);
+      const_cast<WeakSubItemValuePerItem<const unsigned short, CellOfEdge>&>(
+        connectivity.m_edge_local_numbers_in_their_cells) =
+        WeakSubItemValuePerItem<unsigned short, CellOfEdge>(connectivity, number_in_child_array);
+      const_cast<WeakSubItemValuePerItem<const unsigned short, CellOfNode>&>(
+        connectivity.m_node_local_numbers_in_their_cells) =
+        WeakSubItemValuePerItem<unsigned short, CellOfNode>(connectivity, number_in_child_array);
+    }
+  } else if constexpr (ConnectivityType::Dimension == 2) {
+    if constexpr (item_type == ItemType::cell) {
+      if constexpr ((sub_item_type == ItemType::face) or (sub_item_type == ItemType::edge)) {
+        const_cast<WeakSubItemValuePerItem<const unsigned short, FaceOfCell>&>(
+          connectivity.m_cell_local_numbers_in_their_faces) =
+          WeakSubItemValuePerItem<unsigned short, FaceOfCell>(connectivity, number_in_child_array);
+        const_cast<WeakSubItemValuePerItem<const unsigned short, EdgeOfCell>&>(
+          connectivity.m_cell_local_numbers_in_their_edges) =
+          WeakSubItemValuePerItem<unsigned short, EdgeOfCell>(connectivity, number_in_child_array);
+      } else {
+        static_assert(sub_item_type == ItemType::node);
+        const_cast<WeakSubItemValuePerItem<const unsigned short, NodeOfCell>&>(
+          connectivity.m_cell_local_numbers_in_their_nodes) =
+          WeakSubItemValuePerItem<unsigned short, NodeOfCell>(connectivity, number_in_child_array);
+      }
+    } else if constexpr (item_type == ItemType::face or item_type == ItemType::edge) {
+      if constexpr (sub_item_type == ItemType::cell) {
+        const_cast<WeakSubItemValuePerItem<const unsigned short, CellOfFace>&>(
+          connectivity.m_face_local_numbers_in_their_cells) =
+          WeakSubItemValuePerItem<unsigned short, CellOfFace>(connectivity, number_in_child_array);
+        const_cast<WeakSubItemValuePerItem<const unsigned short, CellOfEdge>&>(
+          connectivity.m_edge_local_numbers_in_their_cells) =
+          WeakSubItemValuePerItem<unsigned short, CellOfEdge>(connectivity, number_in_child_array);
+      } else {
+        static_assert(sub_item_type == ItemType::node);
+        const_cast<WeakSubItemValuePerItem<const unsigned short, NodeOfFace>&>(
+          connectivity.m_face_local_numbers_in_their_nodes) =
+          WeakSubItemValuePerItem<unsigned short, NodeOfFace>(connectivity, number_in_child_array);
+        const_cast<WeakSubItemValuePerItem<const unsigned short, NodeOfEdge>&>(
+          connectivity.m_edge_local_numbers_in_their_nodes) =
+          WeakSubItemValuePerItem<unsigned short, NodeOfEdge>(connectivity, number_in_child_array);
+      }
+    } else {
+      static_assert(item_type == ItemType::node);
+      if constexpr (sub_item_type == ItemType::cell) {
+        const_cast<WeakSubItemValuePerItem<const unsigned short, CellOfNode>&>(
+          connectivity.m_node_local_numbers_in_their_cells) =
+          WeakSubItemValuePerItem<unsigned short, CellOfNode>(connectivity, number_in_child_array);
+      } else {
+        static_assert(sub_item_type == ItemType::face or sub_item_type == ItemType::edge);
+        const_cast<WeakSubItemValuePerItem<const unsigned short, FaceOfNode>&>(
+          connectivity.m_node_local_numbers_in_their_faces) =
+          WeakSubItemValuePerItem<unsigned short, FaceOfNode>(connectivity, number_in_child_array);
+        const_cast<WeakSubItemValuePerItem<const unsigned short, EdgeOfNode>&>(
+          connectivity.m_node_local_numbers_in_their_edges) =
+          WeakSubItemValuePerItem<unsigned short, EdgeOfNode>(connectivity, number_in_child_array);
+      }
+    }
+  } else if constexpr (ConnectivityType::Dimension == 3) {
+    if constexpr (item_type == ItemType::cell) {
+      if constexpr (sub_item_type == ItemType::face) {
+        const_cast<WeakSubItemValuePerItem<const unsigned short, FaceOfCell>&>(
+          connectivity.m_cell_local_numbers_in_their_faces) =
+          WeakSubItemValuePerItem<unsigned short, FaceOfCell>(connectivity, number_in_child_array);
+      } else if constexpr (sub_item_type == ItemType::edge) {
+        const_cast<WeakSubItemValuePerItem<const unsigned short, EdgeOfCell>&>(
+          connectivity.m_cell_local_numbers_in_their_edges) =
+          WeakSubItemValuePerItem<unsigned short, EdgeOfCell>(connectivity, number_in_child_array);
+      } else {
+        static_assert(sub_item_type == ItemType::node);
+        const_cast<WeakSubItemValuePerItem<const unsigned short, NodeOfCell>&>(
+          connectivity.m_cell_local_numbers_in_their_nodes) =
+          WeakSubItemValuePerItem<unsigned short, NodeOfCell>(connectivity, number_in_child_array);
+      }
+    } else if constexpr (item_type == ItemType::face) {
+      if constexpr (sub_item_type == ItemType::cell) {
+        const_cast<WeakSubItemValuePerItem<const unsigned short, CellOfFace>&>(
+          connectivity.m_face_local_numbers_in_their_cells) =
+          WeakSubItemValuePerItem<unsigned short, CellOfFace>(connectivity, number_in_child_array);
+      } else if constexpr (sub_item_type == ItemType::edge) {
+        const_cast<WeakSubItemValuePerItem<const unsigned short, EdgeOfFace>&>(
+          connectivity.m_face_local_numbers_in_their_edges) =
+          WeakSubItemValuePerItem<unsigned short, EdgeOfFace>(connectivity, number_in_child_array);
+      } else {
+        static_assert(sub_item_type == ItemType::node);
+        const_cast<WeakSubItemValuePerItem<const unsigned short, NodeOfFace>&>(
+          connectivity.m_face_local_numbers_in_their_nodes) =
+          WeakSubItemValuePerItem<unsigned short, NodeOfFace>(connectivity, number_in_child_array);
+      }
+    } else if constexpr (item_type == ItemType::edge) {
+      if constexpr (sub_item_type == ItemType::cell) {
+        const_cast<WeakSubItemValuePerItem<const unsigned short, CellOfEdge>&>(
+          connectivity.m_edge_local_numbers_in_their_cells) =
+          WeakSubItemValuePerItem<unsigned short, CellOfEdge>(connectivity, number_in_child_array);
+      } else if constexpr (sub_item_type == ItemType::face) {
+        const_cast<WeakSubItemValuePerItem<const unsigned short, FaceOfEdge>&>(
+          connectivity.m_edge_local_numbers_in_their_faces) =
+          WeakSubItemValuePerItem<unsigned short, FaceOfEdge>(connectivity, number_in_child_array);
+      } else {
+        static_assert(sub_item_type == ItemType::node);
+        const_cast<WeakSubItemValuePerItem<const unsigned short, NodeOfEdge>&>(
+          connectivity.m_edge_local_numbers_in_their_nodes) =
+          WeakSubItemValuePerItem<unsigned short, NodeOfEdge>(connectivity, number_in_child_array);
+      }
+    } else {
+      static_assert(item_type == ItemType::node);
+      if constexpr (sub_item_type == ItemType::cell) {
+        const_cast<WeakSubItemValuePerItem<const unsigned short, CellOfNode>&>(
+          connectivity.m_node_local_numbers_in_their_cells) =
+          WeakSubItemValuePerItem<unsigned short, CellOfNode>(connectivity, number_in_child_array);
+      } else if constexpr (sub_item_type == ItemType::face) {
+        const_cast<WeakSubItemValuePerItem<const unsigned short, FaceOfNode>&>(
+          connectivity.m_node_local_numbers_in_their_faces) =
+          WeakSubItemValuePerItem<unsigned short, FaceOfNode>(connectivity, number_in_child_array);
+      } else {
+        static_assert(sub_item_type == ItemType::edge);
+        const_cast<WeakSubItemValuePerItem<const unsigned short, EdgeOfNode>&>(
+          connectivity.m_node_local_numbers_in_their_edges) =
+          WeakSubItemValuePerItem<unsigned short, EdgeOfNode>(connectivity, number_in_child_array);
+      }
+    }
+  }
 }
 
 // 1D
 
-template WeakSubItemValuePerItem<const unsigned short, CellOfNode>
-ConnectivityComputer::computeLocalItemNumberInChildItem<NodeOfCell>(const Connectivity1D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<CellOfFace>(const Connectivity1D&) const;
 
-template WeakSubItemValuePerItem<const unsigned short, NodeOfCell>
-ConnectivityComputer::computeLocalItemNumberInChildItem<CellOfNode>(const Connectivity1D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<CellOfEdge>(const Connectivity1D&) const;
+
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<CellOfNode>(const Connectivity1D&) const;
+
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<FaceOfCell>(const Connectivity1D&) const;
+
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<EdgeOfCell>(const Connectivity1D&) const;
+
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<NodeOfCell>(const Connectivity1D&) const;
 
 // 2D
 
-template WeakSubItemValuePerItem<const unsigned short, CellOfNode>
-ConnectivityComputer::computeLocalItemNumberInChildItem<NodeOfCell>(const Connectivity2D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<CellOfFace>(const Connectivity2D&) const;
+
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<CellOfEdge>(const Connectivity2D&) const;
+
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<CellOfNode>(const Connectivity2D&) const;
+
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<FaceOfCell>(const Connectivity2D&) const;
+
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<FaceOfNode>(const Connectivity2D&) const;
 
-template WeakSubItemValuePerItem<const unsigned short, CellOfFace>
-ConnectivityComputer::computeLocalItemNumberInChildItem<FaceOfCell>(const Connectivity2D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<EdgeOfCell>(const Connectivity2D&) const;
 
-template WeakSubItemValuePerItem<const unsigned short, FaceOfNode>
-ConnectivityComputer::computeLocalItemNumberInChildItem<NodeOfFace>(const Connectivity2D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<EdgeOfNode>(const Connectivity2D&) const;
 
-template WeakSubItemValuePerItem<const unsigned short, FaceOfCell>
-ConnectivityComputer::computeLocalItemNumberInChildItem<CellOfFace>(const Connectivity2D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<NodeOfCell>(const Connectivity2D&) const;
 
-template WeakSubItemValuePerItem<const unsigned short, NodeOfFace>
-ConnectivityComputer::computeLocalItemNumberInChildItem<FaceOfNode>(const Connectivity2D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<NodeOfFace>(const Connectivity2D&) const;
 
-template WeakSubItemValuePerItem<const unsigned short, NodeOfCell>
-ConnectivityComputer::computeLocalItemNumberInChildItem<CellOfNode>(const Connectivity2D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<NodeOfEdge>(const Connectivity2D&) const;
 
 // 3D
 
-template WeakSubItemValuePerItem<const unsigned short, CellOfNode>
-ConnectivityComputer::computeLocalItemNumberInChildItem<NodeOfCell>(const Connectivity3D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<CellOfFace>(const Connectivity3D&) const;
 
-template WeakSubItemValuePerItem<const unsigned short, CellOfEdge>
-ConnectivityComputer::computeLocalItemNumberInChildItem<EdgeOfCell>(const Connectivity3D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<CellOfEdge>(const Connectivity3D&) const;
 
-template WeakSubItemValuePerItem<const unsigned short, CellOfFace>
-ConnectivityComputer::computeLocalItemNumberInChildItem<FaceOfCell>(const Connectivity3D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<CellOfNode>(const Connectivity3D&) const;
 
-template WeakSubItemValuePerItem<const unsigned short, FaceOfNode>
-ConnectivityComputer::computeLocalItemNumberInChildItem<NodeOfFace>(const Connectivity3D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<FaceOfCell>(const Connectivity3D&) const;
 
-template WeakSubItemValuePerItem<const unsigned short, FaceOfEdge>
-ConnectivityComputer::computeLocalItemNumberInChildItem<EdgeOfFace>(const Connectivity3D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<FaceOfEdge>(const Connectivity3D&) const;
 
-template WeakSubItemValuePerItem<const unsigned short, FaceOfCell>
-ConnectivityComputer::computeLocalItemNumberInChildItem<CellOfFace>(const Connectivity3D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<FaceOfNode>(const Connectivity3D&) const;
 
-template WeakSubItemValuePerItem<const unsigned short, EdgeOfNode>
-ConnectivityComputer::computeLocalItemNumberInChildItem<NodeOfEdge>(const Connectivity3D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<EdgeOfCell>(const Connectivity3D&) const;
 
-template WeakSubItemValuePerItem<const unsigned short, EdgeOfFace>
-ConnectivityComputer::computeLocalItemNumberInChildItem<FaceOfEdge>(const Connectivity3D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<EdgeOfFace>(const Connectivity3D&) const;
 
-template WeakSubItemValuePerItem<const unsigned short, EdgeOfCell>
-ConnectivityComputer::computeLocalItemNumberInChildItem<CellOfEdge>(const Connectivity3D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<EdgeOfNode>(const Connectivity3D&) const;
 
-template WeakSubItemValuePerItem<const unsigned short, NodeOfEdge>
-ConnectivityComputer::computeLocalItemNumberInChildItem<EdgeOfNode>(const Connectivity3D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<NodeOfCell>(const Connectivity3D&) const;
 
-template WeakSubItemValuePerItem<const unsigned short, NodeOfFace>
-ConnectivityComputer::computeLocalItemNumberInChildItem<FaceOfNode>(const Connectivity3D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<NodeOfFace>(const Connectivity3D&) const;
 
-template WeakSubItemValuePerItem<const unsigned short, NodeOfCell>
-ConnectivityComputer::computeLocalItemNumberInChildItem<CellOfNode>(const Connectivity3D&) const;
+template void ConnectivityComputer::computeLocalItemNumberInChildItem<NodeOfEdge>(const Connectivity3D&) const;
diff --git a/src/mesh/ConnectivityComputer.hpp b/src/mesh/ConnectivityComputer.hpp
index a4f843e218107b49a6913e3f9028bf034df7c63f..b033a25d5ccc04db84c194e79e492a081216a49e 100644
--- a/src/mesh/ConnectivityComputer.hpp
+++ b/src/mesh/ConnectivityComputer.hpp
@@ -16,8 +16,7 @@ class ConnectivityComputer
                                                ItemType child_item_type) const;
 
   template <typename ItemOfItem, typename ConnectivityType>
-  WeakSubItemValuePerItem<const unsigned short, typename ItemOfItem::Reversed> computeLocalItemNumberInChildItem(
-    const ConnectivityType& connectivity) const;
+  void computeLocalItemNumberInChildItem(const ConnectivityType& connectivity) const;
 
   ConnectivityComputer(const ConnectivityComputer&) = default;
   ConnectivityComputer()                            = default;
diff --git a/src/mesh/ConnectivityMatrix.hpp b/src/mesh/ConnectivityMatrix.hpp
index 576f7935fb11dce55ef20f9fdc5610bfbd27bbff..6b777288cd6b8e3793a9addff6629c661d45d82d 100644
--- a/src/mesh/ConnectivityMatrix.hpp
+++ b/src/mesh/ConnectivityMatrix.hpp
@@ -4,70 +4,94 @@
 #include <utils/Array.hpp>
 #include <utils/PugsUtils.hpp>
 
-#include <Kokkos_StaticCrsGraph.hpp>
-
 class ConnectivityMatrix
 {
+ public:
+  using IndexType = uint32_t;
+
  private:
-  using HostMatrix = Kokkos::StaticCrsGraph<unsigned int, Kokkos::HostSpace>;
-  HostMatrix m_host_matrix;
+  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
-  entries() const
+  const Array<const IndexType>&
+  values() 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
-  numEntries() const
+  size_t
+  numberOfValues() const
   {
-    return m_host_matrix.entries.extent(0);
+    return m_column_indices.size();
   }
 
   PUGS_INLINE
-  auto
-  numRows() const
+  size_t
+  numberOfRows() const
   {
-    return m_host_matrix.numRows();
+    return m_row_map.size() - 1;
   }
 
   PUGS_INLINE
   auto
-  rowConst(size_t j) const
+  operator[](size_t j) const
   {
-    return m_host_matrix.rowConst(j);
+    return subArrayView(m_column_indices, m_row_map[j], m_row_map[j + 1] - m_row_map[j]);
   }
 
   PUGS_INLINE
-  const auto&
-  rowMap(size_t j) const
+  ConnectivityMatrix(const Array<const uint32_t>& row_map,
+                     const Array<const uint32_t>& column_indices) noexcept(NO_ASSERT)
+    : m_row_map{row_map}, m_column_indices{column_indices}, m_is_built{true}
   {
-    return m_host_matrix.row_map[j];
+    Assert(m_column_indices.size() == m_row_map[m_row_map.size() - 1], "incompatible row map and column indices");
+    Assert(m_row_map.size() > 0, "invalid row map");
+    Assert(m_row_map[0] == 0, "row map should start with 0");
+#ifndef NDEBUG
+    for (size_t i = 1; i < m_row_map.size(); ++i) {
+      Assert(m_row_map[i] > m_row_map[i - 1], "row map values must be strictly increasing");
+    }
+#endif   // NDEBUG
   }
 
   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) noexcept : 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/IConnectivity.hpp b/src/mesh/IConnectivity.hpp
index 90c49e4d818e6caa6414992e9ebaf7bdf46f15cc..def3475c386a3e3b321ea640cdf378456f39fac6 100644
--- a/src/mesh/IConnectivity.hpp
+++ b/src/mesh/IConnectivity.hpp
@@ -9,16 +9,15 @@
 
 class IConnectivity : public std::enable_shared_from_this<IConnectivity>
 {
- protected:
+ public:
   template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
   friend class SubItemValuePerItem;
 
   template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
   friend class SubItemArrayPerItem;
 
-  virtual const ConnectivityMatrix& _getMatrix(const ItemType& item_type_0, const ItemType& item_type_1) const = 0;
+  virtual const ConnectivityMatrix& getMatrix(const ItemType& item_type_0, const ItemType& item_type_1) const = 0;
 
- public:
   virtual size_t dimension() const = 0;
 
   std::shared_ptr<const IConnectivity>
diff --git a/src/mesh/ItemArray.hpp b/src/mesh/ItemArray.hpp
index 9213bd3132d94da121cd6c78f8bc0d3bba096922..175537ef27f1d48123f83b0db9a951042b6ceb0e 100644
--- a/src/mesh/ItemArray.hpp
+++ b/src/mesh/ItemArray.hpp
@@ -43,7 +43,7 @@ class ItemArray
   friend ItemArray<std::remove_const_t<DataType>, item_type, ConnectivityWeakPtr>;
 
  public:
-  friend PUGS_INLINE ItemArray<std::remove_const_t<DataType>, item_type, ConnectivityPtr>
+  [[nodiscard]] friend PUGS_INLINE ItemArray<std::remove_const_t<DataType>, item_type, ConnectivityPtr>
   copy(const ItemArray<DataType, item_type, ConnectivityPtr>& source)
   {
     ItemArray<std::remove_const_t<DataType>, item_type, ConnectivityPtr> image;
@@ -62,19 +62,18 @@ class ItemArray
   copy_to(const ItemArray<DataType, item_type, ConnectivityPtr>& source,
           const ItemArray<std::remove_const_t<DataType>, item_type, ConnectivityPtr2>& destination)
   {
-    Assert(destination.connectivity_ptr() == source.connectivity_ptr());
-    copy_to(source.m_values, destination.m_values);
+    Assert(destination.connectivity_ptr() == source.connectivity_ptr(), "different connectivities");
+    Assert(source.sizeOfArrays() == destination.sizeOfArrays(), "incompatible size of arrays")
+      copy_to(source.m_values, destination.m_values);
   }
 
-  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>) {
@@ -92,37 +91,27 @@ class ItemArray
     m_values.fill(data);
   }
 
-  // Following Kokkos logic, these classes are view and const view does allow
-  // changes in data
-  PUGS_FORCEINLINE
-  Array<DataType>
-  operator[](const ItemId& i) const noexcept(NO_ASSERT)
+  template <ItemType item_t>
+  [[nodiscard]] PUGS_INLINE typename Table<DataType>::UnsafeRowView
+  operator[](const ItemIdT<item_t>& item_id) const noexcept(NO_ASSERT)
   {
-    Assert(this->isBuilt());
-    Assert(i < this->numberOfItems());
-    return m_values[i];
+    static_assert(item_t == item_type, "invalid ItemId type");
+    Assert(this->isBuilt(), "ItemArray is not built");
+    Assert(item_id < this->numberOfItems(), "invalid item_id");
+    return m_values[item_id];
   }
 
-  template <typename IndexType>
-  Array<DataType>
-  operator[](const IndexType&) const noexcept(NO_ASSERT)
-  {
-    static_assert(std::is_same_v<IndexType, ItemId>, "ItemArray must be indexed by ItemId");
-  }
-
-  PUGS_INLINE
-  size_t
+  [[nodiscard]] PUGS_INLINE size_t
   numberOfItems() const noexcept(NO_ASSERT)
   {
-    Assert(this->isBuilt());
+    Assert(this->isBuilt(), "ItemArray is not built");
     return m_values.numberOfRows();
   }
 
-  PUGS_INLINE
-  size_t
+  [[nodiscard]] PUGS_INLINE size_t
   sizeOfArrays() const
   {
-    Assert(this->isBuilt());
+    Assert(this->isBuilt(), "ItemArray is not built");
     return m_values.numberOfColumns();
   }
 
@@ -178,7 +167,7 @@ class ItemArray
   ItemArray(const IConnectivity& connectivity, const Table<DataType>& table) noexcept(NO_ASSERT)
     : m_connectivity_ptr{connectivity.shared_ptr()}, m_values{table}
   {
-    Assert(connectivity.numberOf<item_type>() == table.numberOfRows(), "Invalid table, wrong number of rows");
+    Assert(connectivity.numberOf<item_type>() == table.numberOfRows(), "invalid table, wrong number of rows");
   }
 
   PUGS_INLINE
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/ItemId.hpp b/src/mesh/ItemId.hpp
index e62fe7bae6b75f37285257d4211e7ad70d8b6a26..4d7d6714342c62d4c0d3fc35dbb5ae769aafdb8f 100644
--- a/src/mesh/ItemId.hpp
+++ b/src/mesh/ItemId.hpp
@@ -9,7 +9,7 @@ template <ItemType item_type>
 class ItemIdT
 {
  public:
-  using base_type = unsigned int;
+  using base_type = uint32_t;
 
  private:
   base_type m_id;
diff --git a/src/mesh/ItemToItemMatrix.hpp b/src/mesh/ItemToItemMatrix.hpp
index 1cb20d3ee595cf0963b95b4cc67f201cf2cddb5c..6cd2d798b6eaa2844d4eae0d33f8405d344593b7 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,90 +15,80 @@ class ItemToItemMatrix
   using SourceItemId = ItemIdT<SourceItemType>;
   using TargetItemId = ItemIdT<TargetItemType>;
 
-  template <typename RowType>
-  class SubItemList
+  class UnsafeSubItemArray
   {
    private:
-    const RowType m_row;
+    using IndexType = typename ConnectivityMatrix::IndexType;
+
+    const IndexType* const m_values;
+    const 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_values[j];
     }
 
     PUGS_INLINE
-    SubItemList(const RowType& row) : m_row{row}
-    {
-      ;
-    }
+    UnsafeSubItemArray(const ConnectivityMatrix& connectivity_matrix, SourceItemId source_item_id)
+      : m_values{&(connectivity_matrix.values()[connectivity_matrix.rowsMap()[source_item_id]])},
+        m_size(connectivity_matrix.rowsMap()[source_item_id + 1] - connectivity_matrix.rowsMap()[source_item_id])
+    {}
 
     PUGS_INLINE
-    SubItemList& operator=(const SubItemList&) = default;
+    UnsafeSubItemArray& operator=(const UnsafeSubItemArray&) = delete;
 
     PUGS_INLINE
-    SubItemList& operator=(SubItemList&&) = default;
+    UnsafeSubItemArray(const UnsafeSubItemArray&) = delete;
 
     PUGS_INLINE
-    SubItemList(const SubItemList&) = default;
-
-    PUGS_INLINE
-    SubItemList(SubItemList&&) = default;
-
-    PUGS_INLINE
-    ~SubItemList() = default;
+    ~UnsafeSubItemArray() = default;
   };
 
  private:
   const ConnectivityMatrix& m_connectivity_matrix;
 
  public:
-  auto
-  numberOfEntries() const
+  PUGS_INLINE
+  size_t
+  numberOfValues() const
   {
-    return m_connectivity_matrix.numEntries();
+    return m_connectivity_matrix.numberOfValues();
   }
 
+  PUGS_INLINE
   auto
-  entries() const
+  values() const
   {
-    return m_connectivity_matrix.entries();
+    return m_connectivity_matrix.values();
   }
 
   PUGS_INLINE
   auto
   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));
+    Assert(source_id < m_connectivity_matrix.numberOfRows());
+    return UnsafeSubItemArray(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/ItemValue.hpp b/src/mesh/ItemValue.hpp
index e6df4b7059a56ee7922c16f50fcf829cff8d8f4d..496c33242317e49f19842441e9c0fbec0f16ef70 100644
--- a/src/mesh/ItemValue.hpp
+++ b/src/mesh/ItemValue.hpp
@@ -43,7 +43,7 @@ class ItemValue
   friend ItemValue<std::remove_const_t<DataType>, item_type, ConnectivityWeakPtr>;
 
  public:
-  friend PUGS_INLINE ItemValue<std::remove_const_t<DataType>, item_type, ConnectivityPtr>
+  [[nodiscard]] friend PUGS_INLINE ItemValue<std::remove_const_t<DataType>, item_type, ConnectivityPtr>
   copy(const ItemValue<DataType, item_type, ConnectivityPtr>& source)
   {
     ItemValue<std::remove_const_t<DataType>, item_type, ConnectivityPtr> image;
@@ -62,19 +62,17 @@ class ItemValue
   copy_to(const ItemValue<DataType, item_type, ConnectivityPtr>& source,
           const ItemValue<std::remove_const_t<DataType>, item_type, ConnectivityPtr2>& destination)
   {
-    Assert(destination.connectivity_ptr() == source.connectivity_ptr());
+    Assert(destination.connectivity_ptr() == source.connectivity_ptr(), "different connectivities");
     copy_to(source.m_values, destination.m_values);
   }
 
-  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>) {
@@ -84,11 +82,10 @@ class ItemValue
     }
   }
 
-  PUGS_INLINE
-  size_t
+  [[nodiscard]] PUGS_INLINE size_t
   numberOfItems() const noexcept(NO_ASSERT)
   {
-    Assert(this->isBuilt());
+    Assert(this->isBuilt(), "ItemValue is not built");
     return m_values.size();
   }
 
@@ -100,21 +97,14 @@ class ItemValue
     m_values.fill(data);
   }
 
-  // Following Kokkos logic, these classes are view and const view does allow
-  // changes in data
-  PUGS_FORCEINLINE
-  DataType&
-  operator[](const ItemId& i) const noexcept(NO_ASSERT)
-  {
-    Assert(this->isBuilt());
-    return m_values[i];
-  }
-
-  template <typename IndexType>
-  DataType&
-  operator[](const IndexType&) const noexcept(NO_ASSERT)
+  template <ItemType item_t>
+  [[nodiscard]] PUGS_INLINE DataType&
+  operator[](const ItemIdT<item_t>& item_id) const noexcept(NO_ASSERT)
   {
-    static_assert(std::is_same_v<IndexType, ItemId>, "ItemValue must be indexed by ItemId");
+    static_assert(item_t == item_type, "invalid ItemId type");
+    Assert(this->isBuilt(), "ItemValue is not built");
+    Assert(item_id < this->numberOfItems(), "invalid item_id");
+    return m_values[item_id];
   }
 
   template <typename DataType2, typename ConnectivityPtr2>
@@ -168,7 +158,7 @@ class ItemValue
   ItemValue(const IConnectivity& connectivity, const Array<DataType>& values) noexcept(NO_ASSERT)
     : m_connectivity_ptr{connectivity.shared_ptr()}, m_values{values}
   {
-    Assert((m_values.size() == connectivity.numberOf<item_type>()), "Invalid values size");
+    Assert(m_values.size() == connectivity.numberOf<item_type>(), "invalid values size");
   }
 
   PUGS_INLINE
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 8842fd8ac604221d869e046a7d5e9e0975916ad8..ada0f8fbdbcf0b8d2075ebbb8ad05cafa3a63a2e 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -174,8 +174,6 @@ class MeshData : public IMeshData
     } else {
       const NodeValue<const Rd>& xr = m_mesh.xr();
 
-      const CellValue<const double>& inv_cell_nb_nodes = m_mesh.connectivity().invCellNbNodes();
-
       const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
       CellValue<Rd> cell_iso_barycenter{m_mesh.connectivity()};
       parallel_for(
@@ -185,7 +183,7 @@ class MeshData : public IMeshData
           for (size_t R = 0; R < cell_nodes.size(); ++R) {
             X += xr[cell_nodes[R]];
           }
-          cell_iso_barycenter[j] = inv_cell_nb_nodes[j] * X;
+          cell_iso_barycenter[j] = (1. / cell_nodes.size()) * X;
         });
       m_cell_iso_barycenter = cell_iso_barycenter;
     }
diff --git a/src/mesh/PrimalToDiamondDualConnectivityDataMapper.hpp b/src/mesh/PrimalToDiamondDualConnectivityDataMapper.hpp
index ce695a4d3eb4f261677c1e5f2d2ff191890efdc1..a803732dbb59be501f2a6ae995c57505a6dc4b0a 100644
--- a/src/mesh/PrimalToDiamondDualConnectivityDataMapper.hpp
+++ b/src/mesh/PrimalToDiamondDualConnectivityDataMapper.hpp
@@ -116,8 +116,8 @@ class PrimalToDiamondDualConnectivityDataMapper : public IPrimalToDualConnectivi
 
   template <typename OriginDataType, typename DestinationDataType, ItemType destination_face_type>
   void
-  fromDualCell(const CellValue<DestinationDataType>& dual_cell_value,
-               const ItemValue<OriginDataType, destination_face_type>& primal_face_value) const
+  fromDualCell(const CellValue<OriginDataType>& dual_cell_value,
+               const ItemValue<DestinationDataType, destination_face_type>& primal_face_value) const
   {
     static_assert(not std::is_const_v<DestinationDataType>, "destination data type must not be constant");
     static_assert(std::is_same_v<std::remove_const_t<OriginDataType>, DestinationDataType>, "incompatible types");
diff --git a/src/mesh/PrimalToDual1DConnectivityDataMapper.hpp b/src/mesh/PrimalToDual1DConnectivityDataMapper.hpp
index 723025ca793cac1b558a18673a8179298c0ac330..22570810227cbff09c1a72c680f5d34c0f522a44 100644
--- a/src/mesh/PrimalToDual1DConnectivityDataMapper.hpp
+++ b/src/mesh/PrimalToDual1DConnectivityDataMapper.hpp
@@ -142,8 +142,8 @@ class PrimalToDual1DConnectivityDataMapper : public IPrimalToDualConnectivityDat
 
   template <typename OriginDataType, typename DestinationDataType, ItemType destination_node_type>
   void
-  fromDualCell(const CellValue<DestinationDataType>& dual_cell_value,
-               const ItemValue<OriginDataType, destination_node_type>& primal_node_value) const
+  fromDualCell(const CellValue<OriginDataType>& dual_cell_value,
+               const ItemValue<DestinationDataType, destination_node_type>& primal_node_value) const
   {
     static_assert(is_node_in_1d_v<destination_node_type>, "invalid destination node type");
 
diff --git a/src/mesh/PrimalToMedianDualConnectivityDataMapper.hpp b/src/mesh/PrimalToMedianDualConnectivityDataMapper.hpp
index 8f11d66520d975b9c98e336751f5097cec1e4a96..366508eb3e6b09026529b9ff75b9d822e98e4b1f 100644
--- a/src/mesh/PrimalToMedianDualConnectivityDataMapper.hpp
+++ b/src/mesh/PrimalToMedianDualConnectivityDataMapper.hpp
@@ -159,8 +159,8 @@ class PrimalToMedianDualConnectivityDataMapper : public IPrimalToDualConnectivit
 
   template <typename OriginDataType, typename DestinationDataType>
   void
-  fromDualCell(const CellValue<DestinationDataType>& dual_cell_value,
-               const NodeValue<OriginDataType>& primal_node_value) const
+  fromDualCell(const CellValue<OriginDataType>& dual_cell_value,
+               const NodeValue<DestinationDataType>& primal_node_value) const
   {
     static_assert(not std::is_const_v<DestinationDataType>, "destination data type must not be constant");
     static_assert(std::is_same_v<std::remove_const_t<OriginDataType>, DestinationDataType>, "incompatible types");
diff --git a/src/mesh/SubItemArrayPerItem.hpp b/src/mesh/SubItemArrayPerItem.hpp
index 7b204d2193aa6cdb450f7e8d8b23f91d0a47e2c3..a98c039e085cc9539b077bf3a717f94c9fda7667 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,64 +78,69 @@ 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(this->isBuilt(), "SubItemArrayPerItem is not built");
+    Assert(item_id < this->numberOfItems(), "invalid item_id");
+
+    Assert(i + m_row_map[size_t{item_id}] < m_row_map[size_t{item_id} + 1], "invalid index");
+    return m_values[m_row_map[size_t{item_id}] + i];
+  }
+
+  template <ItemType item_t>
+  [[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");
+    return this->itemTable(item_id);
   }
 
-  // Following Kokkos logic, these classes are view and const view does allow
-  // changes in data
   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");
-    Assert(this->isBuilt());
-    Assert(static_cast<size_t>(i) < m_values.numberOfRows());
+    Assert(this->isBuilt(), "SubItemArrayPerItem is not built");
+    Assert(static_cast<size_t>(i) < m_values.numberOfRows(), "invalid index");
     return m_values[i];
   }
 
-  PUGS_INLINE
-  size_t
+  [[nodiscard]] PUGS_INLINE size_t
   numberOfArrays() const noexcept(NO_ASSERT)
   {
-    Assert(this->isBuilt());
+    Assert(this->isBuilt(), "SubItemArrayPerItem is not built");
     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(this->isBuilt(), "SubItemArrayPerItem is not built");
+    Assert(m_row_map.size() > 0, "invalid row map");
+    return m_row_map.size() - 1;
   }
 
-  PUGS_INLINE size_t
+  [[nodiscard]] PUGS_INLINE size_t
   sizeOfArrays() const noexcept(NO_ASSERT)
   {
-    Assert(this->isBuilt());
+    Assert(this->isBuilt(), "SubItemArrayPerItem is not built");
     return m_values.numberOfColumns();
   }
 
   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});
+    Assert(this->isBuilt(), "SubItemArrayPerItem is not built");
+    Assert(item_id < this->numberOfItems(), "invalid 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");
@@ -145,29 +148,15 @@ class SubItemArrayPerItem
   }
 
   template <typename IndexType>
-  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 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 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());
+    Assert(this->isBuilt(), "SubItemArrayPerItem is not built");
+    Assert(item_id < this->numberOfItems(), "invalid item_id");
+    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>
@@ -181,8 +170,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>) {
@@ -210,12 +199,19 @@ class SubItemArrayPerItem
     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);
+    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.numberOfValues(), size_of_arrays);
   }
 
+  SubItemArrayPerItem(const IConnectivity& connectivity, const Table<DataType>& table) noexcept(NO_ASSERT)
+    : m_connectivity_ptr{connectivity.shared_ptr()}, m_values{table}
+  {
+    Assert(m_values.numberOfRows() == connectivity.getMatrix(item_type, sub_item_type).numberOfValues(),
+           "invalid size of provided table");
+    m_row_map = connectivity.getMatrix(item_type, sub_item_type).rowsMap();
+  }
   ~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 ea0e900681279b306ee82a9615e843650cfb1ea7..f5ce39d5f7c67c448c327baf8fea9719d17021fd 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
@@ -48,26 +48,24 @@ class SubItemValuePerItem
   friend SubItemValuePerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityWeakPtr>;
 
  public:
-  friend PUGS_INLINE SubItemValuePerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityPtr>
+  [[nodiscard]] friend PUGS_INLINE SubItemValuePerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityPtr>
   copy(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& source)
   {
     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;
   }
 
-  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>) {
@@ -78,53 +76,58 @@ class SubItemValuePerItem
   }
 
   template <typename IndexType, typename SubIndexType>
-  PUGS_FORCEINLINE DataType&
+  [[nodiscard]] PUGS_INLINE 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 m_values[m_host_row_map(size_t{item_id}) + i];
+    Assert(this->isBuilt(), "SubItemValuePerItem is not built");
+    Assert(item_id < this->numberOfItems(), "invalid item_id");
+    Assert(i + m_row_map[size_t{item_id}] < m_row_map[size_t{item_id} + 1], "invalid index");
+    return m_values[m_row_map[size_t{item_id}] + i];
+  }
+
+  template <ItemType item_t>
+  [[nodiscard]] PUGS_INLINE typename Array<DataType>::UnsafeArrayView
+  operator[](const ItemIdT<item_t>& item_id) const noexcept(NO_ASSERT)
+  {
+    static_assert(item_t == item_type, "invalid ItemId type");
+    return this->itemArray(item_id);
   }
 
-  // Following Kokkos logic, these classes are view and const view does allow
-  // changes in data
   template <typename ArrayIndexType>
-  DataType&
+  [[nodiscard]] PUGS_INLINE 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_values.size());
+    Assert(this->isBuilt(), "SubItemValuePerItem is not built");
+    Assert(static_cast<size_t>(i) < m_values.size(), "invalid index");
     return m_values[i];
   }
 
-  PUGS_INLINE
-  size_t
+  [[nodiscard]] PUGS_INLINE size_t
   numberOfValues() const noexcept(NO_ASSERT)
   {
-    Assert(this->isBuilt());
+    Assert(this->isBuilt(), "SubItemValuePerItem is not built");
     return m_values.size();
   }
 
-  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(this->isBuilt(), "SubItemValuePerItem is not built");
+    Assert(m_row_map.size() > 0, "invalid row map");
+    return m_row_map.size() - 1;
   }
 
   template <typename IndexType>
-  PUGS_INLINE size_t
+  [[nodiscard]] PUGS_INLINE size_t
   numberOfSubValues(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});
+    Assert(this->isBuilt(), "SubItemValuePerItem is not built");
+    Assert(item_id < this->numberOfItems(), "invalid item_id");
+    return m_row_map[size_t{item_id} + 1] - m_row_map[size_t{item_id}];
   }
 
   PUGS_INLINE
@@ -136,29 +139,15 @@ class SubItemValuePerItem
   }
 
   template <typename IndexType>
-  PUGS_INLINE Array<DataType>
-  itemArray(IndexType item_id) noexcept(NO_ASSERT)
+  [[nodiscard]] 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);
-  }
-
-  // Following Kokkos logic, these classes are view and const view does allow
-  // changes in data
-  template <typename IndexType>
-  PUGS_INLINE Array<DataType>
-  itemArray(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);
+    Assert(this->isBuilt(), "SubItemValuePerItem is not built");
+    Assert(item_id < this->numberOfItems(), "invalid item_id");
+    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>
@@ -172,8 +161,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>) {
@@ -200,10 +189,18 @@ class SubItemValuePerItem
     static_assert(not std::is_const_v<DataType>, "Cannot allocate SubItemValuePerItem of const data: only "
                                                  "view is supported");
 
-    ConnectivityMatrix connectivity_matrix = connectivity._getMatrix(item_type, sub_item_type);
+    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<DataType>(connectivity_matrix.numberOfValues());
+  }
+
+  SubItemValuePerItem(const IConnectivity& connectivity, const Array<DataType>& values) noexcept(NO_ASSERT)
+    : m_connectivity_ptr{connectivity.shared_ptr()}, m_values{values}
+  {
+    Assert(m_values.size() == connectivity.getMatrix(item_type, sub_item_type).numberOfValues(),
+           "invalid size of provided values");
+    m_row_map = connectivity.getMatrix(item_type, sub_item_type).rowsMap();
   }
 
   ~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..98173a1402b588ceb9c9584bf5dbe84c136fc8e3 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
@@ -456,7 +456,7 @@ VTKWriter::_write(const std::shared_ptr<const MeshType>& mesh,
     {
       const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
 
-      _write_array(fout, "connectivity", cell_to_node_matrix.entries(), serialize_data_list);
+      _write_array(fout, "connectivity", cell_to_node_matrix.values(), serialize_data_list);
     }
 
     {
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index db0a20dca4e640878a1c6e842b025bd32676a8ff..1ef88c67af2968ba891874da558e206a3dd86964 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -298,21 +298,13 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
 
             bc_list.emplace_back(PressureBoundaryCondition{mesh_node_boundary, node_values});
           } else {
-            constexpr ItemType FaceType = [] {
-              if constexpr (Dimension > 1) {
-                return ItemType::face;
-              } else {
-                return ItemType::node;
-              }
-            }();
-
             MeshFaceBoundary<Dimension> mesh_face_boundary =
               getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
 
             MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
             Array<const double> face_values =
-              InterpolateItemValue<double(Rd)>::template interpolate<FaceType>(pressure_id, mesh_data.xl(),
-                                                                               mesh_face_boundary.faceList());
+              InterpolateItemValue<double(Rd)>::template interpolate<ItemType::face>(pressure_id, mesh_data.xl(),
+                                                                                     mesh_face_boundary.faceList());
             bc_list.emplace_back(PressureBoundaryCondition{mesh_face_boundary, face_values});
           }
 
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_Connectivity.cpp b/tests/test_Connectivity.cpp
index 9abb8fdb3aa0f81bc6123daf70dd467006657f54..f7d2ab7de77a04e3f5eb46efd458d5299740fc36 100644
--- a/tests/test_Connectivity.cpp
+++ b/tests/test_Connectivity.cpp
@@ -12,6 +12,120 @@
 
 TEST_CASE("Connectivity", "[mesh]")
 {
+  SECTION("numberOfItems")
+  {
+    SECTION("1D")
+    {
+      SECTION("unordered 1D mesh")
+      {
+        const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh();
+
+        if (parallel::size() == 1) {
+          REQUIRE(mesh.numberOfNodes() == 35);
+          REQUIRE(mesh.numberOfEdges() == 35);
+          REQUIRE(mesh.numberOfFaces() == 35);
+          REQUIRE(mesh.numberOfCells() == 34);
+        }
+
+        REQUIRE(mesh.numberOfNodes() == mesh.numberOf<ItemType::node>());
+        REQUIRE(mesh.numberOfEdges() == mesh.numberOf<ItemType::edge>());
+        REQUIRE(mesh.numberOfFaces() == mesh.numberOf<ItemType::face>());
+        REQUIRE(mesh.numberOfCells() == mesh.numberOf<ItemType::cell>());
+      }
+
+      SECTION("cartesian 1D mesh")
+      {
+        const auto& mesh = *MeshDataBaseForTests::get().cartesian1DMesh();
+
+        if (parallel::size() == 1) {
+          REQUIRE(mesh.numberOfNodes() == 24);
+          REQUIRE(mesh.numberOfEdges() == 24);
+          REQUIRE(mesh.numberOfFaces() == 24);
+          REQUIRE(mesh.numberOfCells() == 23);
+        }
+
+        REQUIRE(mesh.numberOfNodes() == mesh.numberOf<ItemType::node>());
+        REQUIRE(mesh.numberOfEdges() == mesh.numberOf<ItemType::edge>());
+        REQUIRE(mesh.numberOfFaces() == mesh.numberOf<ItemType::face>());
+        REQUIRE(mesh.numberOfCells() == mesh.numberOf<ItemType::cell>());
+      }
+    }
+
+    SECTION("2D")
+    {
+      SECTION("hybrid 2D mesh")
+      {
+        const auto& mesh = *MeshDataBaseForTests::get().hybrid2DMesh();
+
+        if (parallel::size() == 1) {
+          REQUIRE(mesh.numberOfNodes() == 53);
+          REQUIRE(mesh.numberOfEdges() == 110);
+          REQUIRE(mesh.numberOfFaces() == 110);
+          REQUIRE(mesh.numberOfCells() == 58);
+        }
+
+        REQUIRE(mesh.numberOfNodes() == mesh.numberOf<ItemType::node>());
+        REQUIRE(mesh.numberOfEdges() == mesh.numberOf<ItemType::edge>());
+        REQUIRE(mesh.numberOfFaces() == mesh.numberOf<ItemType::face>());
+        REQUIRE(mesh.numberOfCells() == mesh.numberOf<ItemType::cell>());
+      }
+
+      SECTION("cartesian 2D mesh")
+      {
+        const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh();
+
+        if (parallel::size() == 1) {
+          REQUIRE(mesh.numberOfNodes() == 56);
+          REQUIRE(mesh.numberOfEdges() == 97);
+          REQUIRE(mesh.numberOfFaces() == 97);
+          REQUIRE(mesh.numberOfCells() == 42);
+        }
+
+        REQUIRE(mesh.numberOfNodes() == mesh.numberOf<ItemType::node>());
+        REQUIRE(mesh.numberOfEdges() == mesh.numberOf<ItemType::edge>());
+        REQUIRE(mesh.numberOfFaces() == mesh.numberOf<ItemType::face>());
+        REQUIRE(mesh.numberOfCells() == mesh.numberOf<ItemType::cell>());
+      }
+    }
+
+    SECTION("3D")
+    {
+      SECTION("hybrid 3D mesh")
+      {
+        const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh();
+
+        if (parallel::size() == 1) {
+          REQUIRE(mesh.numberOfNodes() == 132);
+          REQUIRE(mesh.numberOfEdges() == 452);
+          REQUIRE(mesh.numberOfFaces() == 520);
+          REQUIRE(mesh.numberOfCells() == 199);
+        }
+
+        REQUIRE(mesh.numberOfNodes() == mesh.numberOf<ItemType::node>());
+        REQUIRE(mesh.numberOfEdges() == mesh.numberOf<ItemType::edge>());
+        REQUIRE(mesh.numberOfFaces() == mesh.numberOf<ItemType::face>());
+        REQUIRE(mesh.numberOfCells() == mesh.numberOf<ItemType::cell>());
+      }
+
+      SECTION("cartesian 3D mesh")
+      {
+        const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh();
+
+        if (parallel::size() == 1) {
+          REQUIRE(mesh.numberOfNodes() == 280);
+          REQUIRE(mesh.numberOfEdges() == 709);
+          REQUIRE(mesh.numberOfFaces() == 598);
+          REQUIRE(mesh.numberOfCells() == 168);
+        }
+
+        REQUIRE(mesh.numberOfNodes() == mesh.numberOf<ItemType::node>());
+        REQUIRE(mesh.numberOfEdges() == mesh.numberOf<ItemType::edge>());
+        REQUIRE(mesh.numberOfFaces() == mesh.numberOf<ItemType::face>());
+        REQUIRE(mesh.numberOfCells() == mesh.numberOf<ItemType::cell>());
+      }
+    }
+  }
+
   SECTION("isOwned")
   {
     SECTION("1D")
@@ -555,4 +669,305 @@ TEST_CASE("Connectivity", "[mesh]")
       }
     }
   }
+
+  SECTION("ItemLocalNumbersInTheirSubItems")
+  {
+    auto check_item_local_numbers_in_their_subitems = [](auto item_to_subitem_matrix, auto subitem_to_item_matrix,
+                                                         auto item_local_numbers_in_subitems) -> bool {
+      using ItemId    = typename decltype(item_to_subitem_matrix)::SourceItemId;
+      using SubItemId = typename decltype(item_to_subitem_matrix)::TargetItemId;
+
+      static_assert(std::is_same_v<typename decltype(item_to_subitem_matrix)::SourceItemId,
+                                   typename decltype(subitem_to_item_matrix)::TargetItemId>);
+      static_assert(std::is_same_v<typename decltype(item_to_subitem_matrix)::TargetItemId,
+                                   typename decltype(subitem_to_item_matrix)::SourceItemId>);
+      static_assert(std::is_same_v<typename decltype(item_to_subitem_matrix)::SourceItemId,
+                                   typename decltype(item_local_numbers_in_subitems)::ItemId>);
+
+      for (ItemId item_id = 0; item_id < item_local_numbers_in_subitems.numberOfItems(); ++item_id) {
+        auto item_subitem_list = item_to_subitem_matrix[item_id];
+        auto item_numbers      = item_local_numbers_in_subitems.itemArray(item_id);
+
+        for (size_t i_item_subitem = 0; i_item_subitem < item_subitem_list.size(); ++i_item_subitem) {
+          const SubItemId node_cell_id = item_subitem_list[i_item_subitem];
+          const size_t i_local_item    = item_numbers[i_item_subitem];
+          if (subitem_to_item_matrix[node_cell_id][i_local_item] != item_id) {
+            return false;
+          }
+        }
+      }
+      return true;
+    };
+
+    SECTION("1D")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh                           = named_mesh.mesh();
+          const Connectivity<1>& connectivity = mesh->connectivity();
+
+          SECTION("node <-> cell")
+          {
+            auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+            auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
+
+            auto node_local_numbers_in_their_cells = connectivity.nodeLocalNumbersInTheirCells();
+            REQUIRE(check_item_local_numbers_in_their_subitems(node_to_cell_matrix, cell_to_node_matrix,
+                                                               node_local_numbers_in_their_cells));
+
+            auto cell_local_numbers_in_their_nodes = connectivity.cellLocalNumbersInTheirNodes();
+            REQUIRE(check_item_local_numbers_in_their_subitems(cell_to_node_matrix, node_to_cell_matrix,
+                                                               cell_local_numbers_in_their_nodes));
+          }
+
+          SECTION("edge <-> cell")
+          {
+            auto edge_to_cell_matrix = connectivity.edgeToCellMatrix();
+            auto cell_to_edge_matrix = connectivity.cellToEdgeMatrix();
+
+            auto edge_local_numbers_in_their_cells = connectivity.edgeLocalNumbersInTheirCells();
+            REQUIRE(check_item_local_numbers_in_their_subitems(edge_to_cell_matrix, cell_to_edge_matrix,
+                                                               edge_local_numbers_in_their_cells));
+
+            auto cell_local_numbers_in_their_edges = connectivity.cellLocalNumbersInTheirEdges();
+            REQUIRE(check_item_local_numbers_in_their_subitems(cell_to_edge_matrix, edge_to_cell_matrix,
+                                                               cell_local_numbers_in_their_edges));
+          }
+
+          SECTION("face <-> cell")
+          {
+            auto face_to_cell_matrix = connectivity.faceToCellMatrix();
+            auto cell_to_face_matrix = connectivity.cellToFaceMatrix();
+
+            auto face_local_numbers_in_their_cells = connectivity.faceLocalNumbersInTheirCells();
+            REQUIRE(check_item_local_numbers_in_their_subitems(face_to_cell_matrix, cell_to_face_matrix,
+                                                               face_local_numbers_in_their_cells));
+
+            auto cell_local_numbers_in_their_faces = connectivity.cellLocalNumbersInTheirFaces();
+            REQUIRE(check_item_local_numbers_in_their_subitems(cell_to_face_matrix, face_to_cell_matrix,
+                                                               cell_local_numbers_in_their_faces));
+          }
+
+          SECTION("shared data")
+          {
+            auto face_local_numbers_in_their_cells = connectivity.faceLocalNumbersInTheirCells();
+            auto edge_local_numbers_in_their_cells = connectivity.edgeLocalNumbersInTheirCells();
+            auto node_local_numbers_in_their_cells = connectivity.nodeLocalNumbersInTheirCells();
+
+            REQUIRE(&(face_local_numbers_in_their_cells[0]) == &(edge_local_numbers_in_their_cells[0]));
+            REQUIRE(&(face_local_numbers_in_their_cells[0]) == &(node_local_numbers_in_their_cells[0]));
+
+            auto cell_local_numbers_in_their_faces = connectivity.cellLocalNumbersInTheirFaces();
+            auto cell_local_numbers_in_their_edges = connectivity.cellLocalNumbersInTheirEdges();
+            auto cell_local_numbers_in_their_nodes = connectivity.cellLocalNumbersInTheirNodes();
+
+            REQUIRE(&(cell_local_numbers_in_their_faces[0]) == &(cell_local_numbers_in_their_edges[0]));
+            REQUIRE(&(cell_local_numbers_in_their_faces[0]) == &(cell_local_numbers_in_their_nodes[0]));
+          }
+        }
+      }
+    }
+
+    SECTION("2D")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh                           = named_mesh.mesh();
+          const Connectivity<2>& connectivity = mesh->connectivity();
+
+          SECTION("node <-> cell")
+          {
+            auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+            auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
+
+            auto node_local_numbers_in_their_cells = connectivity.nodeLocalNumbersInTheirCells();
+            REQUIRE(check_item_local_numbers_in_their_subitems(node_to_cell_matrix, cell_to_node_matrix,
+                                                               node_local_numbers_in_their_cells));
+
+            auto cell_local_numbers_in_their_nodes = connectivity.cellLocalNumbersInTheirNodes();
+            REQUIRE(check_item_local_numbers_in_their_subitems(cell_to_node_matrix, node_to_cell_matrix,
+                                                               cell_local_numbers_in_their_nodes));
+          }
+
+          SECTION("edge <-> cell")
+          {
+            auto edge_to_cell_matrix = connectivity.edgeToCellMatrix();
+            auto cell_to_edge_matrix = connectivity.cellToEdgeMatrix();
+
+            auto edge_local_numbers_in_their_cells = connectivity.edgeLocalNumbersInTheirCells();
+            REQUIRE(check_item_local_numbers_in_their_subitems(edge_to_cell_matrix, cell_to_edge_matrix,
+                                                               edge_local_numbers_in_their_cells));
+
+            auto cell_local_numbers_in_their_edges = connectivity.cellLocalNumbersInTheirEdges();
+            REQUIRE(check_item_local_numbers_in_their_subitems(cell_to_edge_matrix, edge_to_cell_matrix,
+                                                               cell_local_numbers_in_their_edges));
+          }
+
+          SECTION("face <-> cell")
+          {
+            auto face_to_cell_matrix = connectivity.faceToCellMatrix();
+            auto cell_to_face_matrix = connectivity.cellToFaceMatrix();
+
+            auto face_local_numbers_in_their_cells = connectivity.faceLocalNumbersInTheirCells();
+            REQUIRE(check_item_local_numbers_in_their_subitems(face_to_cell_matrix, cell_to_face_matrix,
+                                                               face_local_numbers_in_their_cells));
+
+            auto cell_local_numbers_in_their_faces = connectivity.cellLocalNumbersInTheirFaces();
+            REQUIRE(check_item_local_numbers_in_their_subitems(cell_to_face_matrix, face_to_cell_matrix,
+                                                               cell_local_numbers_in_their_faces));
+          }
+
+          SECTION("node <-> face")
+          {
+            auto node_to_face_matrix = connectivity.nodeToFaceMatrix();
+            auto face_to_node_matrix = connectivity.faceToNodeMatrix();
+
+            auto node_local_numbers_in_their_faces = connectivity.nodeLocalNumbersInTheirFaces();
+            REQUIRE(check_item_local_numbers_in_their_subitems(node_to_face_matrix, face_to_node_matrix,
+                                                               node_local_numbers_in_their_faces));
+
+            auto face_local_numbers_in_their_nodes = connectivity.faceLocalNumbersInTheirNodes();
+            REQUIRE(check_item_local_numbers_in_their_subitems(face_to_node_matrix, node_to_face_matrix,
+                                                               face_local_numbers_in_their_nodes));
+          }
+
+          SECTION("node <-> edge")
+          {
+            auto node_to_edge_matrix = connectivity.nodeToEdgeMatrix();
+            auto edge_to_node_matrix = connectivity.edgeToNodeMatrix();
+
+            auto node_local_numbers_in_their_edges = connectivity.nodeLocalNumbersInTheirEdges();
+            REQUIRE(check_item_local_numbers_in_their_subitems(node_to_edge_matrix, edge_to_node_matrix,
+                                                               node_local_numbers_in_their_edges));
+
+            auto edge_local_numbers_in_their_nodes = connectivity.edgeLocalNumbersInTheirNodes();
+            REQUIRE(check_item_local_numbers_in_their_subitems(edge_to_node_matrix, node_to_edge_matrix,
+                                                               edge_local_numbers_in_their_nodes));
+          }
+
+          SECTION("shared data")
+          {
+            auto face_local_numbers_in_their_cells = connectivity.faceLocalNumbersInTheirCells();
+            auto edge_local_numbers_in_their_cells = connectivity.edgeLocalNumbersInTheirCells();
+            auto face_local_numbers_in_their_nodes = connectivity.faceLocalNumbersInTheirNodes();
+            auto edge_local_numbers_in_their_nodes = connectivity.edgeLocalNumbersInTheirNodes();
+
+            REQUIRE(&(face_local_numbers_in_their_cells[0]) == &(edge_local_numbers_in_their_cells[0]));
+            REQUIRE(&(face_local_numbers_in_their_nodes[0]) == &(edge_local_numbers_in_their_nodes[0]));
+
+            auto cell_local_numbers_in_their_faces = connectivity.cellLocalNumbersInTheirFaces();
+            auto cell_local_numbers_in_their_edges = connectivity.cellLocalNumbersInTheirEdges();
+            auto node_local_numbers_in_their_faces = connectivity.nodeLocalNumbersInTheirFaces();
+            auto node_local_numbers_in_their_edges = connectivity.nodeLocalNumbersInTheirEdges();
+
+            REQUIRE(&(cell_local_numbers_in_their_faces[0]) == &(cell_local_numbers_in_their_edges[0]));
+            REQUIRE(&(node_local_numbers_in_their_faces[0]) == &(node_local_numbers_in_their_edges[0]));
+          }
+        }
+      }
+    }
+
+    SECTION("3D")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh                           = named_mesh.mesh();
+          const Connectivity<3>& connectivity = mesh->connectivity();
+
+          SECTION("node <-> cell")
+          {
+            auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+            auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
+
+            auto node_local_numbers_in_their_cells = connectivity.nodeLocalNumbersInTheirCells();
+            REQUIRE(check_item_local_numbers_in_their_subitems(node_to_cell_matrix, cell_to_node_matrix,
+                                                               node_local_numbers_in_their_cells));
+
+            auto cell_local_numbers_in_their_nodes = connectivity.cellLocalNumbersInTheirNodes();
+            REQUIRE(check_item_local_numbers_in_their_subitems(cell_to_node_matrix, node_to_cell_matrix,
+                                                               cell_local_numbers_in_their_nodes));
+          }
+
+          SECTION("edge <-> cell")
+          {
+            auto edge_to_cell_matrix = connectivity.edgeToCellMatrix();
+            auto cell_to_edge_matrix = connectivity.cellToEdgeMatrix();
+
+            auto edge_local_numbers_in_their_cells = connectivity.edgeLocalNumbersInTheirCells();
+            REQUIRE(check_item_local_numbers_in_their_subitems(edge_to_cell_matrix, cell_to_edge_matrix,
+                                                               edge_local_numbers_in_their_cells));
+
+            auto cell_local_numbers_in_their_edges = connectivity.cellLocalNumbersInTheirEdges();
+            REQUIRE(check_item_local_numbers_in_their_subitems(cell_to_edge_matrix, edge_to_cell_matrix,
+                                                               cell_local_numbers_in_their_edges));
+          }
+
+          SECTION("face <-> cell")
+          {
+            auto face_to_cell_matrix = connectivity.faceToCellMatrix();
+            auto cell_to_face_matrix = connectivity.cellToFaceMatrix();
+
+            auto face_local_numbers_in_their_cells = connectivity.faceLocalNumbersInTheirCells();
+            REQUIRE(check_item_local_numbers_in_their_subitems(face_to_cell_matrix, cell_to_face_matrix,
+                                                               face_local_numbers_in_their_cells));
+
+            auto cell_local_numbers_in_their_faces = connectivity.cellLocalNumbersInTheirFaces();
+            REQUIRE(check_item_local_numbers_in_their_subitems(cell_to_face_matrix, face_to_cell_matrix,
+                                                               cell_local_numbers_in_their_faces));
+          }
+
+          SECTION("node <-> face")
+          {
+            auto node_to_face_matrix = connectivity.nodeToFaceMatrix();
+            auto face_to_node_matrix = connectivity.faceToNodeMatrix();
+
+            auto node_local_numbers_in_their_faces = connectivity.nodeLocalNumbersInTheirFaces();
+            REQUIRE(check_item_local_numbers_in_their_subitems(node_to_face_matrix, face_to_node_matrix,
+                                                               node_local_numbers_in_their_faces));
+
+            auto face_local_numbers_in_their_nodes = connectivity.faceLocalNumbersInTheirNodes();
+            REQUIRE(check_item_local_numbers_in_their_subitems(face_to_node_matrix, node_to_face_matrix,
+                                                               face_local_numbers_in_their_nodes));
+          }
+
+          SECTION("edge <-> face")
+          {
+            auto edge_to_face_matrix = connectivity.edgeToFaceMatrix();
+            auto face_to_edge_matrix = connectivity.faceToEdgeMatrix();
+
+            auto edge_local_numbers_in_their_faces = connectivity.edgeLocalNumbersInTheirFaces();
+            REQUIRE(check_item_local_numbers_in_their_subitems(edge_to_face_matrix, face_to_edge_matrix,
+                                                               edge_local_numbers_in_their_faces));
+
+            auto face_local_numbers_in_their_edges = connectivity.faceLocalNumbersInTheirEdges();
+            REQUIRE(check_item_local_numbers_in_their_subitems(face_to_edge_matrix, edge_to_face_matrix,
+                                                               face_local_numbers_in_their_edges));
+          }
+
+          SECTION("node <-> edge")
+          {
+            auto node_to_edge_matrix = connectivity.nodeToEdgeMatrix();
+            auto edge_to_node_matrix = connectivity.edgeToNodeMatrix();
+
+            auto node_local_numbers_in_their_edges = connectivity.nodeLocalNumbersInTheirEdges();
+            REQUIRE(check_item_local_numbers_in_their_subitems(node_to_edge_matrix, edge_to_node_matrix,
+                                                               node_local_numbers_in_their_edges));
+
+            auto edge_local_numbers_in_their_nodes = connectivity.edgeLocalNumbersInTheirNodes();
+            REQUIRE(check_item_local_numbers_in_their_subitems(edge_to_node_matrix, node_to_edge_matrix,
+                                                               edge_local_numbers_in_their_nodes));
+          }
+        }
+      }
+    }
+  }
 }
diff --git a/tests/test_ItemArray.cpp b/tests/test_ItemArray.cpp
index b98c12ba2528f2933a6ec9f7475fa559fa7d055b..997e327a1237752463e5b8ace7ef9bb1a167fcf7 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);
         }
@@ -297,16 +297,16 @@ TEST_CASE("ItemArray", "[mesh]")
     SECTION("checking for build ItemArray")
     {
       CellArray<int> cell_array;
-      REQUIRE_THROWS_AS(cell_array[CellId{0}], AssertError);
+      REQUIRE_THROWS_WITH(cell_array[CellId{0}], "ItemArray is not built");
 
       FaceArray<int> face_array;
-      REQUIRE_THROWS_AS(face_array[FaceId{0}], AssertError);
+      REQUIRE_THROWS_WITH(face_array[FaceId{0}], "ItemArray is not built");
 
       EdgeArray<int> edge_array;
-      REQUIRE_THROWS_AS(edge_array[EdgeId{0}], AssertError);
+      REQUIRE_THROWS_WITH(edge_array[EdgeId{0}], "ItemArray is not built");
 
       NodeArray<int> node_array;
-      REQUIRE_THROWS_AS(node_array[NodeId{0}], AssertError);
+      REQUIRE_THROWS_WITH(node_array[NodeId{0}], "ItemArray is not built");
     }
 
     SECTION("checking for bounds violation")
@@ -322,19 +322,19 @@ TEST_CASE("ItemArray", "[mesh]")
 
           CellArray<int> cell_array{connectivity, 1};
           CellId invalid_cell_id = connectivity.numberOfCells();
-          REQUIRE_THROWS_AS(cell_array[invalid_cell_id], AssertError);
+          REQUIRE_THROWS_WITH(cell_array[invalid_cell_id], "invalid item_id");
 
           FaceArray<int> face_array{connectivity, 2};
           FaceId invalid_face_id = connectivity.numberOfFaces();
-          REQUIRE_THROWS_AS(face_array[invalid_face_id], AssertError);
+          REQUIRE_THROWS_WITH(face_array[invalid_face_id], "invalid item_id");
 
           EdgeArray<int> edge_array{connectivity, 1};
           EdgeId invalid_edge_id = connectivity.numberOfEdges();
-          REQUIRE_THROWS_AS(edge_array[invalid_edge_id], AssertError);
+          REQUIRE_THROWS_WITH(edge_array[invalid_edge_id], "invalid item_id");
 
           NodeArray<int> node_array{connectivity, 0};
           NodeId invalid_node_id = connectivity.numberOfNodes();
-          REQUIRE_THROWS_AS(node_array[invalid_node_id], AssertError);
+          REQUIRE_THROWS_WITH(node_array[invalid_node_id], "invalid item_id");
         }
       }
     }
@@ -351,7 +351,39 @@ TEST_CASE("ItemArray", "[mesh]")
           const Connectivity<3>& connectivity = mesh_3d->connectivity();
 
           Table<size_t> values{connectivity.numberOfCells() + 3, 3};
-          REQUIRE_THROWS_WITH(CellArray<size_t>(connectivity, values), "Invalid table, wrong number of rows");
+          REQUIRE_THROWS_WITH(CellArray<size_t>(connectivity, values), "invalid table, wrong number of rows");
+        }
+      }
+    }
+
+    SECTION("invalid copy_to")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_2d = named_mesh.mesh();
+
+          const Connectivity<2>& connectivity_2d = mesh_2d->connectivity();
+
+          std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+          for (auto named_mesh : mesh_list) {
+            SECTION(named_mesh.name())
+            {
+              auto mesh_3d = named_mesh.mesh();
+
+              const Connectivity<3>& connectivity_3d = mesh_3d->connectivity();
+
+              CellArray<int> cell_2d_array{connectivity_2d, 3};
+              CellArray<int> cell_3d_array{connectivity_3d, 3};
+              REQUIRE_THROWS_WITH(copy_to(cell_2d_array, cell_3d_array), "different connectivities");
+
+              CellArray<int> cell_2d_array2{connectivity_2d, 2};
+              REQUIRE_THROWS_WITH(copy_to(cell_2d_array, cell_2d_array2), "incompatible size of arrays");
+            }
+          }
         }
       }
     }
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_ItemValue.cpp b/tests/test_ItemValue.cpp
index be1307593d0ed8bff2c916d3162cccc34fb88536..702110f7a906654bb5e60676e1e0e2d4c7a7fc03 100644
--- a/tests/test_ItemValue.cpp
+++ b/tests/test_ItemValue.cpp
@@ -241,16 +241,16 @@ TEST_CASE("ItemValue", "[mesh]")
     SECTION("checking for build ItemValue")
     {
       CellValue<int> cell_value;
-      REQUIRE_THROWS_AS(cell_value[CellId{0}], AssertError);
+      REQUIRE_THROWS_WITH(cell_value[CellId{0}], "ItemValue is not built");
 
       FaceValue<int> face_value;
-      REQUIRE_THROWS_AS(face_value[FaceId{0}], AssertError);
+      REQUIRE_THROWS_WITH(face_value[FaceId{0}], "ItemValue is not built");
 
       EdgeValue<int> edge_value;
-      REQUIRE_THROWS_AS(edge_value[EdgeId{0}], AssertError);
+      REQUIRE_THROWS_WITH(edge_value[EdgeId{0}], "ItemValue is not built");
 
       NodeValue<int> node_value;
-      REQUIRE_THROWS_AS(node_value[NodeId{0}], AssertError);
+      REQUIRE_THROWS_WITH(node_value[NodeId{0}], "ItemValue is not built");
     }
 
     SECTION("checking for bounds violation")
@@ -266,19 +266,19 @@ TEST_CASE("ItemValue", "[mesh]")
 
           CellValue<int> cell_value{connectivity};
           CellId invalid_cell_id = connectivity.numberOfCells();
-          REQUIRE_THROWS_AS(cell_value[invalid_cell_id], AssertError);
+          REQUIRE_THROWS_WITH(cell_value[invalid_cell_id], "invalid item_id");
 
           FaceValue<int> face_value{connectivity};
           FaceId invalid_face_id = connectivity.numberOfFaces();
-          REQUIRE_THROWS_AS(face_value[invalid_face_id], AssertError);
+          REQUIRE_THROWS_WITH(face_value[invalid_face_id], "invalid item_id");
 
           EdgeValue<int> edge_value{connectivity};
           EdgeId invalid_edge_id = connectivity.numberOfEdges();
-          REQUIRE_THROWS_AS(edge_value[invalid_edge_id], AssertError);
+          REQUIRE_THROWS_WITH(edge_value[invalid_edge_id], "invalid item_id");
 
           NodeValue<int> node_value{connectivity};
           NodeId invalid_node_id = connectivity.numberOfNodes();
-          REQUIRE_THROWS_AS(node_value[invalid_node_id], AssertError);
+          REQUIRE_THROWS_WITH(node_value[invalid_node_id], "invalid item_id");
         }
       }
     }
@@ -295,7 +295,7 @@ TEST_CASE("ItemValue", "[mesh]")
           const Connectivity<3>& connectivity = mesh_3d->connectivity();
 
           Array<size_t> values{3 + connectivity.numberOfCells()};
-          REQUIRE_THROWS_WITH(CellValue<size_t>(connectivity, values), "Invalid values size");
+          REQUIRE_THROWS_WITH(CellValue<size_t>(connectivity, values), "invalid values size");
         }
       }
     }
@@ -322,7 +322,7 @@ TEST_CASE("ItemValue", "[mesh]")
 
               CellValue<int> cell_2d_value{connectivity_2d};
               CellValue<int> cell_3d_value{connectivity_3d};
-              REQUIRE_THROWS_AS(copy_to(cell_2d_value, cell_3d_value), AssertError);
+              REQUIRE_THROWS_WITH(copy_to(cell_2d_value, cell_3d_value), "different connectivities");
             }
           }
         }
diff --git a/tests/test_SubItemArrayPerItem.cpp b/tests/test_SubItemArrayPerItem.cpp
index 71a226f01379cccd86ca48c27683ca0b467f80fa..aa8930269fcd03756cf69495ef4f77bc44869e62 100644
--- a/tests/test_SubItemArrayPerItem.cpp
+++ b/tests/test_SubItemArrayPerItem.cpp
@@ -1030,40 +1030,57 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
     SECTION("checking for build SubItemArrayPerItem")
     {
       CellArrayPerNode<int> cell_array_per_node;
-      REQUIRE_THROWS_AS(cell_array_per_node[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.numberOfArrays(), AssertError);
-      REQUIRE_THROWS_AS(cell_array_per_node.numberOfItems(), AssertError);
-      REQUIRE_THROWS_AS(cell_array_per_node.numberOfSubArrays(NodeId{0}), AssertError);
+      REQUIRE_THROWS_WITH(cell_array_per_node[0], "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(cell_array_per_node.itemTable(NodeId{0}), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(cell_array_per_node(NodeId{0}, 0), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(cell_array_per_node.sizeOfArrays(), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(cell_array_per_node.numberOfArrays(), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(cell_array_per_node.numberOfItems(), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(cell_array_per_node.numberOfSubArrays(NodeId{0}), "SubItemArrayPerItem is not built");
 
       FaceArrayPerCell<int> face_array_per_cell;
-      REQUIRE_THROWS_AS(face_array_per_cell[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.numberOfArrays(), AssertError);
-      REQUIRE_THROWS_AS(face_array_per_cell.numberOfItems(), AssertError);
-      REQUIRE_THROWS_AS(face_array_per_cell.numberOfSubArrays(CellId{0}), AssertError);
+      REQUIRE_THROWS_WITH(face_array_per_cell[0], "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(face_array_per_cell.itemTable(CellId{0}), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(face_array_per_cell(CellId{0}, 0), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(face_array_per_cell.sizeOfArrays(), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(face_array_per_cell.numberOfArrays(), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(face_array_per_cell.numberOfItems(), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(face_array_per_cell.numberOfSubArrays(CellId{0}), "SubItemArrayPerItem is not built");
 
       CellArrayPerEdge<int> cell_array_per_edge;
-      REQUIRE_THROWS_AS(cell_array_per_edge[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.numberOfArrays(), AssertError);
-      REQUIRE_THROWS_AS(cell_array_per_edge.numberOfItems(), AssertError);
-      REQUIRE_THROWS_AS(cell_array_per_edge.numberOfSubArrays(EdgeId{0}), AssertError);
+      REQUIRE_THROWS_WITH(cell_array_per_edge[0], "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(cell_array_per_edge.itemTable(EdgeId{0}), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(cell_array_per_edge(EdgeId{0}, 0), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(cell_array_per_edge.sizeOfArrays(), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(cell_array_per_edge.numberOfArrays(), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(cell_array_per_edge.numberOfItems(), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(cell_array_per_edge.numberOfSubArrays(EdgeId{0}), "SubItemArrayPerItem is not built");
 
       NodeArrayPerFace<int> node_array_per_face;
-      REQUIRE_THROWS_AS(node_array_per_face[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.numberOfArrays(), AssertError);
-      REQUIRE_THROWS_AS(node_array_per_face.numberOfItems(), AssertError);
-      REQUIRE_THROWS_AS(node_array_per_face.numberOfSubArrays(FaceId{0}), AssertError);
+      REQUIRE_THROWS_WITH(node_array_per_face[0], "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(node_array_per_face.itemTable(FaceId{0}), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(node_array_per_face(FaceId{0}, 0), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(node_array_per_face.sizeOfArrays(), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(node_array_per_face.numberOfArrays(), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(node_array_per_face.numberOfItems(), "SubItemArrayPerItem is not built");
+      REQUIRE_THROWS_WITH(node_array_per_face.numberOfSubArrays(FaceId{0}), "SubItemArrayPerItem is not built");
+    }
+
+    SECTION("checking invalid table size in constructor")
+    {
+      auto mesh_3d = MeshDataBaseForTests::get().hybrid3DMesh();
+
+      const Connectivity<3>& connectivity = mesh_3d->connectivity();
+
+      {
+        Table<int> table{connectivity.getMatrix(ItemType::node, ItemType::cell).numberOfValues() + 3, 2};
+        REQUIRE_THROWS_WITH(CellArrayPerNode<int>(connectivity, table), "invalid size of provided table");
+      }
+
+      {
+        Table<int> table{connectivity.getMatrix(ItemType::face, ItemType::edge).numberOfValues() + 3, 2};
+        REQUIRE_THROWS_WITH(EdgeArrayPerFace<int>(connectivity, table), "invalid size of provided table");
+      }
     }
 
     SECTION("checking for bounds violation")
@@ -1080,69 +1097,73 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
           CellArrayPerFace<int> cell_array_per_face{connectivity, 3};
           {
             FaceId invalid_face_id = connectivity.numberOfFaces();
-            REQUIRE_THROWS_AS(cell_array_per_face(invalid_face_id, 0), AssertError);
+            REQUIRE_THROWS_WITH(cell_array_per_face(invalid_face_id, 0), "invalid item_id");
+            REQUIRE_THROWS_WITH(cell_array_per_face[invalid_face_id], "invalid item_id");
           }
           if (connectivity.numberOfFaces() > 0) {
             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.numberOfRows()), AssertError);
-            REQUIRE_THROWS_AS(face_table[face_table.numberOfRows()], AssertError);
-            REQUIRE_THROWS_AS(cell_array_per_face.itemTable(face_id)[face_table.numberOfRows()], 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);
+            REQUIRE_THROWS_WITH(cell_array_per_face(face_id, face_table.numberOfRows()), "invalid index");
+            REQUIRE_THROWS_WITH(face_table[face_table.numberOfRows()], "invalid index");
+            REQUIRE_THROWS_WITH(cell_array_per_face.itemTable(face_id)[face_table.numberOfRows()], "invalid index");
+            REQUIRE_THROWS_WITH(cell_array_per_face.itemTable(face_id)[0][cell_array_per_face.sizeOfArrays()],
+                                "invalid index");
+            REQUIRE_THROWS_WITH(cell_array_per_face.itemTable(face_id)[0][cell_array_per_face.sizeOfArrays()] = 2,
+                                "invalid index");
           }
 
           FaceArrayPerNode<int> face_array_per_node{connectivity, 5};
           {
             NodeId invalid_node_id = connectivity.numberOfNodes();
-            REQUIRE_THROWS_AS(face_array_per_node(invalid_node_id, 0), AssertError);
+            REQUIRE_THROWS_WITH(face_array_per_node(invalid_node_id, 0), "invalid item_id");
+            REQUIRE_THROWS_WITH(face_array_per_node[invalid_node_id], "invalid item_id");
           }
           if (connectivity.numberOfNodes() > 0) {
             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.numberOfRows()), AssertError);
-            REQUIRE_THROWS_AS(node_table[node_table.numberOfRows()], AssertError);
-            REQUIRE_THROWS_AS(face_array_per_node.itemTable(node_id)[node_table.numberOfRows()], 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);
+            REQUIRE_THROWS_WITH(face_array_per_node(node_id, node_table.numberOfRows()), "invalid index");
+            REQUIRE_THROWS_WITH(node_table[node_table.numberOfRows()], "invalid index");
+            REQUIRE_THROWS_WITH(face_array_per_node.itemTable(node_id)[node_table.numberOfRows()], "invalid index");
+            REQUIRE_THROWS_WITH(face_array_per_node.itemTable(node_id)[0][face_array_per_node.sizeOfArrays()],
+                                "invalid index");
+            REQUIRE_THROWS_WITH(face_array_per_node.itemTable(node_id)[0][face_array_per_node.sizeOfArrays()] = 2,
+                                "invalid index");
           }
 
           EdgeArrayPerCell<int> edge_array_per_cell{connectivity, 3};
           {
             CellId invalid_cell_id = connectivity.numberOfCells();
-            REQUIRE_THROWS_AS(edge_array_per_cell(invalid_cell_id, 0), AssertError);
+            REQUIRE_THROWS_WITH(edge_array_per_cell(invalid_cell_id, 0), "invalid item_id");
+            REQUIRE_THROWS_WITH(edge_array_per_cell[invalid_cell_id], "invalid item_id");
           }
           if (connectivity.numberOfCells() > 0) {
             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.numberOfRows()), AssertError);
-            REQUIRE_THROWS_AS(cell_table[cell_table.numberOfRows()], AssertError);
-            REQUIRE_THROWS_AS(edge_array_per_cell.itemTable(cell_id)[cell_table.numberOfRows()], 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);
+            REQUIRE_THROWS_WITH(edge_array_per_cell(cell_id, cell_table.numberOfRows()), "invalid index");
+            REQUIRE_THROWS_WITH(cell_table[cell_table.numberOfRows()], "invalid index");
+            REQUIRE_THROWS_WITH(edge_array_per_cell.itemTable(cell_id)[cell_table.numberOfRows()], "invalid index");
+            REQUIRE_THROWS_WITH(edge_array_per_cell.itemTable(cell_id)[0][edge_array_per_cell.sizeOfArrays()],
+                                "invalid index");
+            REQUIRE_THROWS_WITH(edge_array_per_cell.itemTable(cell_id)[0][edge_array_per_cell.sizeOfArrays()] == 2,
+                                "invalid index");
           }
 
           NodeArrayPerEdge<int> node_array_per_edge{connectivity, 3};
           {
             EdgeId invalid_edge_id = connectivity.numberOfEdges();
-            REQUIRE_THROWS_AS(node_array_per_edge(invalid_edge_id, 0), AssertError);
+            REQUIRE_THROWS_WITH(node_array_per_edge(invalid_edge_id, 0), "invalid item_id");
+            REQUIRE_THROWS_WITH(node_array_per_edge[invalid_edge_id], "invalid item_id");
           }
           if (connectivity.numberOfEdges() > 0) {
             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.numberOfRows()), AssertError);
-            REQUIRE_THROWS_AS(edge_table[edge_table.numberOfRows()], AssertError);
-            REQUIRE_THROWS_AS(node_array_per_edge.itemTable(edge_id)[edge_table.numberOfRows()], 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);
+            REQUIRE_THROWS_WITH(node_array_per_edge(edge_id, edge_table.numberOfRows()), "invalid index");
+            REQUIRE_THROWS_WITH(edge_table[edge_table.numberOfRows()], "invalid index");
+            REQUIRE_THROWS_WITH(node_array_per_edge.itemTable(edge_id)[edge_table.numberOfRows()], "invalid index");
+            REQUIRE_THROWS_WITH(node_array_per_edge.itemTable(edge_id)[0][node_array_per_edge.sizeOfArrays()],
+                                "invalid index");
+            REQUIRE_THROWS_WITH(node_array_per_edge.itemTable(edge_id)[0][node_array_per_edge.sizeOfArrays()] = 2,
+                                "invalid index");
           }
         }
       }
diff --git a/tests/test_SubItemValuePerItem.cpp b/tests/test_SubItemValuePerItem.cpp
index dda38918acdb156f8a7c5de552d1f82319463c97..08d421fc1333ae8344402f9806446ab22c8de3b4 100644
--- a/tests/test_SubItemValuePerItem.cpp
+++ b/tests/test_SubItemValuePerItem.cpp
@@ -62,6 +62,15 @@ TEST_CASE("SubItemValuePerItem", "[mesh]")
       return number;
     };
 
+    auto check_same_memory = [](const auto& sub_item_value_per_item, auto array) {
+      REQUIRE(sub_item_value_per_item.numberOfValues() == array.size());
+      bool is_same = true;
+      for (size_t i = 0; i < sub_item_value_per_item.numberOfValues(); ++i) {
+        is_same &= (&(sub_item_value_per_item[i]) == &(array[i]));
+      }
+      return is_same;
+    };
+
     SECTION("1D")
     {
       std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
@@ -77,7 +86,12 @@ TEST_CASE("SubItemValuePerItem", "[mesh]")
           {
             NodeValuePerCell<int> node_value_per_cell{connectivity};
             REQUIRE(node_value_per_cell.numberOfItems() == connectivity.numberOfCells());
-            REQUIRE(node_value_per_cell.numberOfValues() == number_of_values(node_value_per_cell));
+            const size_t nb_values = number_of_values(node_value_per_cell);
+            REQUIRE(node_value_per_cell.numberOfValues() == nb_values);
+
+            Array<int> array{nb_values};
+            NodeValuePerCell<int> node_value_per_cell_from_array{connectivity, array};
+            REQUIRE(check_same_memory(node_value_per_cell_from_array, array));
 
             auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
             {
@@ -131,7 +145,12 @@ TEST_CASE("SubItemValuePerItem", "[mesh]")
           {
             CellValuePerFace<int> cell_value_per_face{connectivity};
             REQUIRE(cell_value_per_face.numberOfItems() == connectivity.numberOfFaces());
-            REQUIRE(cell_value_per_face.numberOfValues() == number_of_values(cell_value_per_face));
+            const size_t nb_values = number_of_values(cell_value_per_face);
+            REQUIRE(cell_value_per_face.numberOfValues() == nb_values);
+
+            Array<int> array{nb_values};
+            CellValuePerFace<int> cell_value_per_face_from_array{connectivity, array};
+            REQUIRE(check_same_memory(cell_value_per_face_from_array, array));
 
             auto face_to_cell_matrix = connectivity.faceToCellMatrix();
             {
@@ -147,7 +166,12 @@ TEST_CASE("SubItemValuePerItem", "[mesh]")
           {
             CellValuePerEdge<int> cell_value_per_edge{connectivity};
             REQUIRE(cell_value_per_edge.numberOfItems() == connectivity.numberOfEdges());
-            REQUIRE(cell_value_per_edge.numberOfValues() == number_of_values(cell_value_per_edge));
+            const size_t nb_values = number_of_values(cell_value_per_edge);
+            REQUIRE(cell_value_per_edge.numberOfValues() == nb_values);
+
+            Array<int> array{nb_values};
+            CellValuePerEdge<int> cell_value_per_edge_from_array{connectivity, array};
+            REQUIRE(check_same_memory(cell_value_per_edge_from_array, array));
 
             auto edge_to_cell_matrix = connectivity.edgeToCellMatrix();
             {
@@ -163,7 +187,12 @@ TEST_CASE("SubItemValuePerItem", "[mesh]")
           {
             CellValuePerNode<int> cell_value_per_node{connectivity};
             REQUIRE(cell_value_per_node.numberOfItems() == connectivity.numberOfNodes());
-            REQUIRE(cell_value_per_node.numberOfValues() == number_of_values(cell_value_per_node));
+            const size_t nb_values = number_of_values(cell_value_per_node);
+            REQUIRE(cell_value_per_node.numberOfValues() == nb_values);
+
+            Array<int> array{nb_values};
+            CellValuePerNode<int> cell_value_per_node_from_array{connectivity, array};
+            REQUIRE(check_same_memory(cell_value_per_node_from_array, array));
 
             auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
             {
@@ -193,7 +222,12 @@ TEST_CASE("SubItemValuePerItem", "[mesh]")
           {
             NodeValuePerCell<int> node_value_per_cell{connectivity};
             REQUIRE(node_value_per_cell.numberOfItems() == connectivity.numberOfCells());
-            REQUIRE(node_value_per_cell.numberOfValues() == number_of_values(node_value_per_cell));
+            const size_t nb_values = number_of_values(node_value_per_cell);
+            REQUIRE(node_value_per_cell.numberOfValues() == nb_values);
+
+            Array<int> array{nb_values};
+            NodeValuePerCell<int> node_value_per_cell_from_array{connectivity, array};
+            REQUIRE(check_same_memory(node_value_per_cell_from_array, array));
 
             auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
             {
@@ -235,7 +269,12 @@ TEST_CASE("SubItemValuePerItem", "[mesh]")
           {
             CellValuePerFace<int> cell_value_per_face{connectivity};
             REQUIRE(cell_value_per_face.numberOfItems() == connectivity.numberOfFaces());
-            REQUIRE(cell_value_per_face.numberOfValues() == number_of_values(cell_value_per_face));
+            const size_t nb_values = number_of_values(cell_value_per_face);
+            REQUIRE(cell_value_per_face.numberOfValues() == nb_values);
+
+            Array<int> array{nb_values};
+            CellValuePerFace<int> cell_value_per_face_from_array{connectivity, array};
+            REQUIRE(check_same_memory(cell_value_per_face_from_array, array));
 
             auto face_to_cell_matrix = connectivity.faceToCellMatrix();
             {
@@ -264,7 +303,12 @@ TEST_CASE("SubItemValuePerItem", "[mesh]")
           {
             CellValuePerEdge<int> cell_value_per_edge{connectivity};
             REQUIRE(cell_value_per_edge.numberOfItems() == connectivity.numberOfEdges());
-            REQUIRE(cell_value_per_edge.numberOfValues() == number_of_values(cell_value_per_edge));
+            const size_t nb_values = number_of_values(cell_value_per_edge);
+            REQUIRE(cell_value_per_edge.numberOfValues() == nb_values);
+
+            Array<int> array{nb_values};
+            CellValuePerEdge<int> cell_value_per_edge_from_array{connectivity, array};
+            REQUIRE(check_same_memory(cell_value_per_edge_from_array, array));
 
             auto edge_to_cell_matrix = connectivity.edgeToCellMatrix();
             {
@@ -293,7 +337,12 @@ TEST_CASE("SubItemValuePerItem", "[mesh]")
           {
             EdgeValuePerNode<int> edge_value_per_node{connectivity};
             REQUIRE(edge_value_per_node.numberOfItems() == connectivity.numberOfNodes());
-            REQUIRE(edge_value_per_node.numberOfValues() == number_of_values(edge_value_per_node));
+            const size_t nb_values = number_of_values(edge_value_per_node);
+            REQUIRE(edge_value_per_node.numberOfValues() == nb_values);
+
+            Array<int> array{nb_values};
+            EdgeValuePerNode<int> cell_value_per_node_from_array{connectivity, array};
+            REQUIRE(check_same_memory(cell_value_per_node_from_array, array));
 
             auto node_to_edge_matrix = connectivity.nodeToEdgeMatrix();
             {
@@ -349,7 +398,12 @@ TEST_CASE("SubItemValuePerItem", "[mesh]")
           {
             NodeValuePerCell<int> node_value_per_cell{connectivity};
             REQUIRE(node_value_per_cell.numberOfItems() == connectivity.numberOfCells());
-            REQUIRE(node_value_per_cell.numberOfValues() == number_of_values(node_value_per_cell));
+            const size_t nb_values = number_of_values(node_value_per_cell);
+            REQUIRE(node_value_per_cell.numberOfValues() == nb_values);
+
+            Array<int> array{nb_values};
+            NodeValuePerCell<int> node_value_per_cell_from_array{connectivity, array};
+            REQUIRE(check_same_memory(node_value_per_cell_from_array, array));
 
             auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
             {
@@ -391,7 +445,12 @@ TEST_CASE("SubItemValuePerItem", "[mesh]")
           {
             CellValuePerFace<int> cell_value_per_face{connectivity};
             REQUIRE(cell_value_per_face.numberOfItems() == connectivity.numberOfFaces());
-            REQUIRE(cell_value_per_face.numberOfValues() == number_of_values(cell_value_per_face));
+            const size_t nb_values = number_of_values(cell_value_per_face);
+            REQUIRE(cell_value_per_face.numberOfValues() == nb_values);
+
+            Array<int> array{nb_values};
+            CellValuePerFace<int> cell_value_per_face_from_array{connectivity, array};
+            REQUIRE(check_same_memory(cell_value_per_face_from_array, array));
 
             auto face_to_cell_matrix = connectivity.faceToCellMatrix();
             {
@@ -433,7 +492,12 @@ TEST_CASE("SubItemValuePerItem", "[mesh]")
           {
             CellValuePerEdge<int> cell_value_per_edge{connectivity};
             REQUIRE(cell_value_per_edge.numberOfItems() == connectivity.numberOfEdges());
-            REQUIRE(cell_value_per_edge.numberOfValues() == number_of_values(cell_value_per_edge));
+            const size_t nb_values = number_of_values(cell_value_per_edge);
+            REQUIRE(cell_value_per_edge.numberOfValues() == nb_values);
+
+            Array<int> array{nb_values};
+            CellValuePerEdge<int> cell_value_per_edge_from_array{connectivity, array};
+            REQUIRE(check_same_memory(cell_value_per_edge_from_array, array));
 
             auto edge_to_cell_matrix = connectivity.edgeToCellMatrix();
             {
@@ -475,7 +539,12 @@ TEST_CASE("SubItemValuePerItem", "[mesh]")
           {
             EdgeValuePerNode<int> edge_value_per_node{connectivity};
             REQUIRE(edge_value_per_node.numberOfItems() == connectivity.numberOfNodes());
-            REQUIRE(edge_value_per_node.numberOfValues() == number_of_values(edge_value_per_node));
+            const size_t nb_values = number_of_values(edge_value_per_node);
+            REQUIRE(edge_value_per_node.numberOfValues() == nb_values);
+
+            Array<int> array{nb_values};
+            EdgeValuePerNode<int> cell_value_per_node_from_array{connectivity, array};
+            REQUIRE(check_same_memory(cell_value_per_node_from_array, array));
 
             auto node_to_edge_matrix = connectivity.nodeToEdgeMatrix();
             {
@@ -517,6 +586,529 @@ TEST_CASE("SubItemValuePerItem", "[mesh]")
     }
   }
 
+  SECTION("access functions")
+  {
+    SECTION("1D")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_1d = named_mesh.mesh();
+
+          const Connectivity<1>& connectivity = mesh_1d->connectivity();
+
+          SECTION("per cell")
+          {
+            {
+              NodeValuePerCell<int> node_value_per_cell{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
+                for (size_t i = 0; i < node_value_per_cell.numberOfSubValues(cell_id); ++i) {
+                  is_valid &= (&(node_value_per_cell[cell_id][i]) == &(node_value_per_cell.itemArray(cell_id)[i]));
+                  is_valid &= (&(node_value_per_cell[cell_id][i]) == &(node_value_per_cell(cell_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == node_value_per_cell.numberOfValues());
+            }
+
+            {
+              EdgeValuePerCell<int> edge_value_per_cell{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
+                for (size_t i = 0; i < edge_value_per_cell.numberOfSubValues(cell_id); ++i) {
+                  is_valid &= (&(edge_value_per_cell[cell_id][i]) == &(edge_value_per_cell.itemArray(cell_id)[i]));
+                  is_valid &= (&(edge_value_per_cell[cell_id][i]) == &(edge_value_per_cell(cell_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == edge_value_per_cell.numberOfValues());
+            }
+
+            {
+              FaceValuePerCell<int> face_value_per_cell{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
+                for (size_t i = 0; i < face_value_per_cell.numberOfSubValues(cell_id); ++i) {
+                  is_valid &= (&(face_value_per_cell[cell_id][i]) == &(face_value_per_cell.itemArray(cell_id)[i]));
+                  is_valid &= (&(face_value_per_cell[cell_id][i]) == &(face_value_per_cell(cell_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == face_value_per_cell.numberOfValues());
+            }
+          }
+
+          SECTION("per face")
+          {
+            CellValuePerFace<int> cell_value_per_face{connectivity};
+
+            bool is_valid  = true;
+            size_t counter = 0;
+            for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+              for (size_t i = 0; i < cell_value_per_face.numberOfSubValues(face_id); ++i) {
+                is_valid &= (&(cell_value_per_face[face_id][i]) == &(cell_value_per_face.itemArray(face_id)[i]));
+                is_valid &= (&(cell_value_per_face[face_id][i]) == &(cell_value_per_face(face_id, i)));
+                counter++;
+              }
+            }
+            REQUIRE(is_valid);
+            REQUIRE(counter == cell_value_per_face.numberOfValues());
+          }
+
+          SECTION("per edge")
+          {
+            CellValuePerEdge<int> cell_value_per_edge{connectivity};
+
+            bool is_valid  = true;
+            size_t counter = 0;
+            for (EdgeId edge_id = 0; edge_id < connectivity.numberOfEdges(); ++edge_id) {
+              for (size_t i = 0; i < cell_value_per_edge.numberOfSubValues(edge_id); ++i) {
+                is_valid &= (&(cell_value_per_edge[edge_id][i]) == &(cell_value_per_edge.itemArray(edge_id)[i]));
+                is_valid &= (&(cell_value_per_edge[edge_id][i]) == &(cell_value_per_edge(edge_id, i)));
+                counter++;
+              }
+            }
+            REQUIRE(is_valid);
+            REQUIRE(counter == cell_value_per_edge.numberOfValues());
+          }
+
+          SECTION("per node")
+          {
+            CellValuePerNode<int> cell_value_per_node{connectivity};
+
+            bool is_valid  = true;
+            size_t counter = 0;
+            for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+              for (size_t i = 0; i < cell_value_per_node.numberOfSubValues(node_id); ++i) {
+                is_valid &= (&(cell_value_per_node[node_id][i]) == &(cell_value_per_node.itemArray(node_id)[i]));
+                is_valid &= (&(cell_value_per_node[node_id][i]) == &(cell_value_per_node(node_id, i)));
+                counter++;
+              }
+            }
+            REQUIRE(is_valid);
+            REQUIRE(counter == cell_value_per_node.numberOfValues());
+          }
+        }
+      }
+    }
+
+    SECTION("2D")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_2d = named_mesh.mesh();
+
+          const Connectivity<2>& connectivity = mesh_2d->connectivity();
+
+          SECTION("per cell")
+          {
+            {
+              NodeValuePerCell<int> node_value_per_cell{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
+                for (size_t i = 0; i < node_value_per_cell.numberOfSubValues(cell_id); ++i) {
+                  is_valid &= (&(node_value_per_cell[cell_id][i]) == &(node_value_per_cell.itemArray(cell_id)[i]));
+                  is_valid &= (&(node_value_per_cell[cell_id][i]) == &(node_value_per_cell(cell_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == node_value_per_cell.numberOfValues());
+            }
+
+            {
+              EdgeValuePerCell<int> edge_value_per_cell{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
+                for (size_t i = 0; i < edge_value_per_cell.numberOfSubValues(cell_id); ++i) {
+                  is_valid &= (&(edge_value_per_cell[cell_id][i]) == &(edge_value_per_cell.itemArray(cell_id)[i]));
+                  is_valid &= (&(edge_value_per_cell[cell_id][i]) == &(edge_value_per_cell(cell_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == edge_value_per_cell.numberOfValues());
+            }
+
+            {
+              FaceValuePerCell<int> face_value_per_cell{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
+                for (size_t i = 0; i < face_value_per_cell.numberOfSubValues(cell_id); ++i) {
+                  is_valid &= (&(face_value_per_cell[cell_id][i]) == &(face_value_per_cell.itemArray(cell_id)[i]));
+                  is_valid &= (&(face_value_per_cell[cell_id][i]) == &(face_value_per_cell(cell_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == face_value_per_cell.numberOfValues());
+            }
+          }
+
+          SECTION("per face")
+          {
+            {
+              NodeValuePerFace<int> node_value_per_face{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+                for (size_t i = 0; i < node_value_per_face.numberOfSubValues(face_id); ++i) {
+                  is_valid &= (&(node_value_per_face[face_id][i]) == &(node_value_per_face.itemArray(face_id)[i]));
+                  is_valid &= (&(node_value_per_face[face_id][i]) == &(node_value_per_face(face_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == node_value_per_face.numberOfValues());
+            }
+
+            {
+              CellValuePerFace<int> cell_value_per_face{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+                for (size_t i = 0; i < cell_value_per_face.numberOfSubValues(face_id); ++i) {
+                  is_valid &= (&(cell_value_per_face[face_id][i]) == &(cell_value_per_face.itemArray(face_id)[i]));
+                  is_valid &= (&(cell_value_per_face[face_id][i]) == &(cell_value_per_face(face_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == cell_value_per_face.numberOfValues());
+            }
+          }
+
+          SECTION("per edge")
+          {
+            {
+              NodeValuePerEdge<int> node_value_per_edge{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (EdgeId edge_id = 0; edge_id < connectivity.numberOfEdges(); ++edge_id) {
+                for (size_t i = 0; i < node_value_per_edge.numberOfSubValues(edge_id); ++i) {
+                  is_valid &= (&(node_value_per_edge[edge_id][i]) == &(node_value_per_edge.itemArray(edge_id)[i]));
+                  is_valid &= (&(node_value_per_edge[edge_id][i]) == &(node_value_per_edge(edge_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == node_value_per_edge.numberOfValues());
+            }
+
+            {
+              CellValuePerEdge<int> cell_value_per_edge{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (EdgeId edge_id = 0; edge_id < connectivity.numberOfEdges(); ++edge_id) {
+                for (size_t i = 0; i < cell_value_per_edge.numberOfSubValues(edge_id); ++i) {
+                  is_valid &= (&(cell_value_per_edge[edge_id][i]) == &(cell_value_per_edge.itemArray(edge_id)[i]));
+                  is_valid &= (&(cell_value_per_edge[edge_id][i]) == &(cell_value_per_edge(edge_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == cell_value_per_edge.numberOfValues());
+            }
+          }
+
+          SECTION("per node")
+          {
+            {
+              EdgeValuePerNode<int> edge_value_per_node{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+                for (size_t i = 0; i < edge_value_per_node.numberOfSubValues(node_id); ++i) {
+                  is_valid &= (&(edge_value_per_node[node_id][i]) == &(edge_value_per_node.itemArray(node_id)[i]));
+                  is_valid &= (&(edge_value_per_node[node_id][i]) == &(edge_value_per_node(node_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == edge_value_per_node.numberOfValues());
+            }
+
+            {
+              FaceValuePerNode<int> face_value_per_node{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+                for (size_t i = 0; i < face_value_per_node.numberOfSubValues(node_id); ++i) {
+                  is_valid &= (&(face_value_per_node[node_id][i]) == &(face_value_per_node.itemArray(node_id)[i]));
+                  is_valid &= (&(face_value_per_node[node_id][i]) == &(face_value_per_node(node_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == face_value_per_node.numberOfValues());
+            }
+
+            {
+              CellValuePerNode<int> cell_value_per_node{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+                for (size_t i = 0; i < cell_value_per_node.numberOfSubValues(node_id); ++i) {
+                  is_valid &= (&(cell_value_per_node[node_id][i]) == &(cell_value_per_node.itemArray(node_id)[i]));
+                  is_valid &= (&(cell_value_per_node[node_id][i]) == &(cell_value_per_node(node_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == cell_value_per_node.numberOfValues());
+            }
+          }
+        }
+      }
+    }
+
+    SECTION("3D")
+    {
+      std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_3d = named_mesh.mesh();
+
+          const Connectivity<3>& connectivity = mesh_3d->connectivity();
+
+          SECTION("per cell")
+          {
+            {
+              NodeValuePerCell<int> node_value_per_cell{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
+                for (size_t i = 0; i < node_value_per_cell.numberOfSubValues(cell_id); ++i) {
+                  is_valid &= (&(node_value_per_cell[cell_id][i]) == &(node_value_per_cell.itemArray(cell_id)[i]));
+                  is_valid &= (&(node_value_per_cell[cell_id][i]) == &(node_value_per_cell(cell_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == node_value_per_cell.numberOfValues());
+            }
+
+            {
+              EdgeValuePerCell<int> edge_value_per_cell{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
+                for (size_t i = 0; i < edge_value_per_cell.numberOfSubValues(cell_id); ++i) {
+                  is_valid &= (&(edge_value_per_cell[cell_id][i]) == &(edge_value_per_cell.itemArray(cell_id)[i]));
+                  is_valid &= (&(edge_value_per_cell[cell_id][i]) == &(edge_value_per_cell(cell_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == edge_value_per_cell.numberOfValues());
+            }
+
+            {
+              FaceValuePerCell<int> face_value_per_cell{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
+                for (size_t i = 0; i < face_value_per_cell.numberOfSubValues(cell_id); ++i) {
+                  is_valid &= (&(face_value_per_cell[cell_id][i]) == &(face_value_per_cell.itemArray(cell_id)[i]));
+                  is_valid &= (&(face_value_per_cell[cell_id][i]) == &(face_value_per_cell(cell_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == face_value_per_cell.numberOfValues());
+            }
+          }
+
+          SECTION("per face")
+          {
+            {
+              NodeValuePerFace<int> node_value_per_face{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+                for (size_t i = 0; i < node_value_per_face.numberOfSubValues(face_id); ++i) {
+                  is_valid &= (&(node_value_per_face[face_id][i]) == &(node_value_per_face.itemArray(face_id)[i]));
+                  is_valid &= (&(node_value_per_face[face_id][i]) == &(node_value_per_face(face_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == node_value_per_face.numberOfValues());
+            }
+
+            {
+              EdgeValuePerFace<int> edge_value_per_face{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+                for (size_t i = 0; i < edge_value_per_face.numberOfSubValues(face_id); ++i) {
+                  is_valid &= (&(edge_value_per_face[face_id][i]) == &(edge_value_per_face.itemArray(face_id)[i]));
+                  is_valid &= (&(edge_value_per_face[face_id][i]) == &(edge_value_per_face(face_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == edge_value_per_face.numberOfValues());
+            }
+
+            {
+              CellValuePerFace<int> cell_value_per_face{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+                for (size_t i = 0; i < cell_value_per_face.numberOfSubValues(face_id); ++i) {
+                  is_valid &= (&(cell_value_per_face[face_id][i]) == &(cell_value_per_face.itemArray(face_id)[i]));
+                  is_valid &= (&(cell_value_per_face[face_id][i]) == &(cell_value_per_face(face_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == cell_value_per_face.numberOfValues());
+            }
+          }
+
+          SECTION("per edge")
+          {
+            {
+              NodeValuePerEdge<int> node_value_per_edge{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (EdgeId edge_id = 0; edge_id < connectivity.numberOfEdges(); ++edge_id) {
+                for (size_t i = 0; i < node_value_per_edge.numberOfSubValues(edge_id); ++i) {
+                  is_valid &= (&(node_value_per_edge[edge_id][i]) == &(node_value_per_edge.itemArray(edge_id)[i]));
+                  is_valid &= (&(node_value_per_edge[edge_id][i]) == &(node_value_per_edge(edge_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == node_value_per_edge.numberOfValues());
+            }
+
+            {
+              FaceValuePerEdge<int> face_value_per_edge{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (EdgeId edge_id = 0; edge_id < connectivity.numberOfEdges(); ++edge_id) {
+                for (size_t i = 0; i < face_value_per_edge.numberOfSubValues(edge_id); ++i) {
+                  is_valid &= (&(face_value_per_edge[edge_id][i]) == &(face_value_per_edge.itemArray(edge_id)[i]));
+                  is_valid &= (&(face_value_per_edge[edge_id][i]) == &(face_value_per_edge(edge_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == face_value_per_edge.numberOfValues());
+            }
+
+            {
+              CellValuePerEdge<int> cell_value_per_edge{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (EdgeId edge_id = 0; edge_id < connectivity.numberOfEdges(); ++edge_id) {
+                for (size_t i = 0; i < cell_value_per_edge.numberOfSubValues(edge_id); ++i) {
+                  is_valid &= (&(cell_value_per_edge[edge_id][i]) == &(cell_value_per_edge.itemArray(edge_id)[i]));
+                  is_valid &= (&(cell_value_per_edge[edge_id][i]) == &(cell_value_per_edge(edge_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == cell_value_per_edge.numberOfValues());
+            }
+          }
+
+          SECTION("per node")
+          {
+            {
+              EdgeValuePerNode<int> edge_value_per_node{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+                for (size_t i = 0; i < edge_value_per_node.numberOfSubValues(node_id); ++i) {
+                  is_valid &= (&(edge_value_per_node[node_id][i]) == &(edge_value_per_node.itemArray(node_id)[i]));
+                  is_valid &= (&(edge_value_per_node[node_id][i]) == &(edge_value_per_node(node_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == edge_value_per_node.numberOfValues());
+            }
+
+            {
+              FaceValuePerNode<int> face_value_per_node{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+                for (size_t i = 0; i < face_value_per_node.numberOfSubValues(node_id); ++i) {
+                  is_valid &= (&(face_value_per_node[node_id][i]) == &(face_value_per_node.itemArray(node_id)[i]));
+                  is_valid &= (&(face_value_per_node[node_id][i]) == &(face_value_per_node(node_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == face_value_per_node.numberOfValues());
+            }
+
+            {
+              CellValuePerNode<int> cell_value_per_node{connectivity};
+
+              bool is_valid  = true;
+              size_t counter = 0;
+              for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+                for (size_t i = 0; i < cell_value_per_node.numberOfSubValues(node_id); ++i) {
+                  is_valid &= (&(cell_value_per_node[node_id][i]) == &(cell_value_per_node.itemArray(node_id)[i]));
+                  is_valid &= (&(cell_value_per_node[node_id][i]) == &(cell_value_per_node(node_id, i)));
+                  counter++;
+                }
+              }
+              REQUIRE(is_valid);
+              REQUIRE(counter == cell_value_per_node.numberOfValues());
+            }
+          }
+        }
+      }
+    }
+  }
+
   SECTION("array view")
   {
     SECTION("1D")
@@ -789,36 +1381,57 @@ TEST_CASE("SubItemValuePerItem", "[mesh]")
     SECTION("checking for build SubItemValuePerItem")
     {
       CellValuePerNode<int> cell_value_per_node;
-      REQUIRE_THROWS_AS(cell_value_per_node[0], AssertError);
-      REQUIRE_THROWS_AS(cell_value_per_node.itemArray(NodeId{0}), AssertError);
-      REQUIRE_THROWS_AS(cell_value_per_node(NodeId{0}, 0), AssertError);
-      REQUIRE_THROWS_AS(cell_value_per_node.numberOfValues(), AssertError);
-      REQUIRE_THROWS_AS(cell_value_per_node.numberOfItems(), AssertError);
-      REQUIRE_THROWS_AS(cell_value_per_node.numberOfSubValues(NodeId{0}), AssertError);
+      REQUIRE_THROWS_WITH(cell_value_per_node[0], "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(cell_value_per_node.itemArray(NodeId{0}), "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(cell_value_per_node[NodeId{0}], "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(cell_value_per_node(NodeId{0}, 0), "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(cell_value_per_node.numberOfValues(), "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(cell_value_per_node.numberOfItems(), "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(cell_value_per_node.numberOfSubValues(NodeId{0}), "SubItemValuePerItem is not built");
 
       FaceValuePerCell<int> face_value_per_cell;
-      REQUIRE_THROWS_AS(face_value_per_cell[0], AssertError);
-      REQUIRE_THROWS_AS(face_value_per_cell.itemArray(CellId{0}), AssertError);
-      REQUIRE_THROWS_AS(face_value_per_cell(CellId{0}, 0), AssertError);
-      REQUIRE_THROWS_AS(face_value_per_cell.numberOfValues(), AssertError);
-      REQUIRE_THROWS_AS(face_value_per_cell.numberOfItems(), AssertError);
-      REQUIRE_THROWS_AS(face_value_per_cell.numberOfSubValues(CellId{0}), AssertError);
+      REQUIRE_THROWS_WITH(face_value_per_cell[0], "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(face_value_per_cell.itemArray(CellId{0}), "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(face_value_per_cell[CellId{0}], "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(face_value_per_cell(CellId{0}, 0), "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(face_value_per_cell.numberOfValues(), "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(face_value_per_cell.numberOfItems(), "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(face_value_per_cell.numberOfSubValues(CellId{0}), "SubItemValuePerItem is not built");
 
       CellValuePerEdge<int> cell_value_per_edge;
-      REQUIRE_THROWS_AS(cell_value_per_edge[0], AssertError);
-      REQUIRE_THROWS_AS(cell_value_per_edge.itemArray(EdgeId{0}), AssertError);
-      REQUIRE_THROWS_AS(cell_value_per_edge(EdgeId{0}, 0), AssertError);
-      REQUIRE_THROWS_AS(cell_value_per_edge.numberOfValues(), AssertError);
-      REQUIRE_THROWS_AS(cell_value_per_edge.numberOfItems(), AssertError);
-      REQUIRE_THROWS_AS(cell_value_per_edge.numberOfSubValues(EdgeId{0}), AssertError);
+      REQUIRE_THROWS_WITH(cell_value_per_edge[0], "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(cell_value_per_edge.itemArray(EdgeId{0}), "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(cell_value_per_edge[EdgeId{0}], "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(cell_value_per_edge(EdgeId{0}, 0), "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(cell_value_per_edge.numberOfValues(), "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(cell_value_per_edge.numberOfItems(), "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(cell_value_per_edge.numberOfSubValues(EdgeId{0}), "SubItemValuePerItem is not built");
 
       NodeValuePerFace<int> node_value_per_face;
-      REQUIRE_THROWS_AS(node_value_per_face[0], AssertError);
-      REQUIRE_THROWS_AS(node_value_per_face.itemArray(FaceId{0}), AssertError);
-      REQUIRE_THROWS_AS(node_value_per_face(FaceId{0}, 0), AssertError);
-      REQUIRE_THROWS_AS(node_value_per_face.numberOfValues(), AssertError);
-      REQUIRE_THROWS_AS(node_value_per_face.numberOfItems(), AssertError);
-      REQUIRE_THROWS_AS(node_value_per_face.numberOfSubValues(FaceId{0}), AssertError);
+      REQUIRE_THROWS_WITH(node_value_per_face[0], "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(node_value_per_face.itemArray(FaceId{0}), "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(node_value_per_face[FaceId{0}], "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(node_value_per_face(FaceId{0}, 0), "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(node_value_per_face.numberOfValues(), "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(node_value_per_face.numberOfItems(), "SubItemValuePerItem is not built");
+      REQUIRE_THROWS_WITH(node_value_per_face.numberOfSubValues(FaceId{0}), "SubItemValuePerItem is not built");
+    }
+
+    SECTION("checking invalid array size in constructor")
+    {
+      auto mesh_3d = MeshDataBaseForTests::get().hybrid3DMesh();
+
+      const Connectivity<3>& connectivity = mesh_3d->connectivity();
+
+      {
+        Array<int> array{connectivity.getMatrix(ItemType::node, ItemType::cell).numberOfValues() + 3};
+        REQUIRE_THROWS_WITH(CellValuePerNode<int>(connectivity, array), "invalid size of provided values");
+      }
+
+      {
+        Array<int> array{connectivity.getMatrix(ItemType::face, ItemType::edge).numberOfValues() + 3};
+        REQUIRE_THROWS_WITH(EdgeValuePerFace<int>(connectivity, array), "invalid size of provided values");
+      }
     }
 
     SECTION("checking for bounds violation")
@@ -835,53 +1448,61 @@ TEST_CASE("SubItemValuePerItem", "[mesh]")
           CellValuePerFace<int> cell_value_per_face{connectivity};
           {
             FaceId invalid_face_id = connectivity.numberOfFaces();
-            REQUIRE_THROWS_AS(cell_value_per_face(invalid_face_id, 0), AssertError);
+            REQUIRE_THROWS_WITH(cell_value_per_face(invalid_face_id, 0), "invalid item_id");
+            REQUIRE_THROWS_WITH(cell_value_per_face[invalid_face_id], "invalid item_id");
           }
           if (connectivity.numberOfFaces() > 0) {
             FaceId face_id          = 0;
             const auto& cell_values = cell_value_per_face.itemArray(face_id);
-            REQUIRE_THROWS_AS(cell_value_per_face(face_id, cell_values.size()), AssertError);
-            REQUIRE_THROWS_AS(cell_values[cell_values.size()], AssertError);
-            REQUIRE_THROWS_AS(cell_value_per_face.itemArray(face_id)[cell_values.size()] = 2, AssertError);
+            REQUIRE_THROWS_WITH(cell_value_per_face(face_id, cell_values.size()), "invalid index");
+            REQUIRE_THROWS_WITH(cell_values[cell_values.size()], "invalid index");
+            REQUIRE_THROWS_WITH(cell_value_per_face.itemArray(face_id)[cell_values.size()] = 2, "invalid index");
+            REQUIRE_THROWS_WITH(cell_value_per_face[face_id][cell_values.size()] = 2, "invalid index");
           }
 
           FaceValuePerNode<int> face_value_per_node{connectivity};
           {
             NodeId invalid_node_id = connectivity.numberOfNodes();
-            REQUIRE_THROWS_AS(face_value_per_node(invalid_node_id, 0), AssertError);
+            REQUIRE_THROWS_WITH(face_value_per_node(invalid_node_id, 0), "invalid item_id");
+            REQUIRE_THROWS_WITH(face_value_per_node[invalid_node_id], "invalid item_id");
           }
           if (connectivity.numberOfNodes() > 0) {
             NodeId node_id          = 0;
             const auto& face_values = face_value_per_node.itemArray(node_id);
-            REQUIRE_THROWS_AS(face_value_per_node(node_id, face_values.size()), AssertError);
-            REQUIRE_THROWS_AS(face_values[face_values.size()], AssertError);
-            REQUIRE_THROWS_AS(face_value_per_node.itemArray(node_id)[face_values.size()] = 2, AssertError);
+            REQUIRE_THROWS_WITH(face_value_per_node(node_id, face_values.size()), "invalid index");
+            REQUIRE_THROWS_WITH(face_values[face_values.size()], "invalid index");
+            REQUIRE_THROWS_WITH(face_value_per_node.itemArray(node_id)[face_values.size()] = 2, "invalid index");
+            REQUIRE_THROWS_WITH(face_value_per_node[node_id][face_values.size()] = 2, "invalid index");
           }
 
           EdgeValuePerCell<int> edge_value_per_cell{connectivity};
           {
             CellId invalid_cell_id = connectivity.numberOfCells();
-            REQUIRE_THROWS_AS(edge_value_per_cell(invalid_cell_id, 0), AssertError);
+            REQUIRE_THROWS_WITH(edge_value_per_cell(invalid_cell_id, 0), "invalid item_id");
+            REQUIRE_THROWS_WITH(edge_value_per_cell[invalid_cell_id], "invalid item_id");
           }
           if (connectivity.numberOfCells() > 0) {
             CellId cell_id          = 0;
             const auto& edge_values = edge_value_per_cell.itemArray(cell_id);
-            REQUIRE_THROWS_AS(edge_value_per_cell(cell_id, edge_values.size()), AssertError);
-            REQUIRE_THROWS_AS(edge_values[edge_values.size()], AssertError);
-            REQUIRE_THROWS_AS(edge_value_per_cell.itemArray(cell_id)[edge_values.size()] = 2, AssertError);
+            REQUIRE_THROWS_WITH(edge_value_per_cell(cell_id, edge_values.size()), "invalid index");
+            REQUIRE_THROWS_WITH(edge_values[edge_values.size()], "invalid index");
+            REQUIRE_THROWS_WITH(edge_value_per_cell.itemArray(cell_id)[edge_values.size()] = 2, "invalid index");
+            REQUIRE_THROWS_WITH(edge_value_per_cell[cell_id][edge_values.size()] = 2, "invalid index");
           }
 
           NodeValuePerEdge<int> node_value_per_edge{connectivity};
           {
             EdgeId invalid_edge_id = connectivity.numberOfEdges();
-            REQUIRE_THROWS_AS(node_value_per_edge(invalid_edge_id, 0), AssertError);
+            REQUIRE_THROWS_WITH(node_value_per_edge(invalid_edge_id, 0), "invalid item_id");
+            REQUIRE_THROWS_WITH(node_value_per_edge[invalid_edge_id], "invalid item_id");
           }
           if (connectivity.numberOfEdges() > 0) {
             EdgeId edge_id          = 0;
             const auto& node_values = node_value_per_edge.itemArray(edge_id);
-            REQUIRE_THROWS_AS(node_value_per_edge(edge_id, node_values.size()), AssertError);
-            REQUIRE_THROWS_AS(node_values[node_values.size()], AssertError);
-            REQUIRE_THROWS_AS(node_value_per_edge.itemArray(edge_id)[node_values.size()] = 2, AssertError);
+            REQUIRE_THROWS_WITH(node_value_per_edge(edge_id, node_values.size()), "invalid index");
+            REQUIRE_THROWS_WITH(node_values[node_values.size()], "invalid index");
+            REQUIRE_THROWS_WITH(node_value_per_edge.itemArray(edge_id)[node_values.size()] = 2, "invalid index");
+            REQUIRE_THROWS_WITH(node_value_per_edge[edge_id][node_values.size()] = 2, "invalid index");
           }
         }
       }
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
 }