diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index fde78fbbff3ee3eeddd9c25952d5bbc2d908c227..c52dab2062ad14637fbeb435f90b3277d9173c66 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -158,11 +158,9 @@ void Connectivity<3>::_computeFaceCellConnectivities()
   m_connectivity_computer.computeInverseConnectivityMatrix(cell_to_face_matrix,
                                                            face_to_cell_matrix);
 
-  CellValuePerFace<unsigned short> face_to_cell_local_face(*this);
-  m_connectivity_computer.computeLocalChildItemNumberInItem(cell_to_face_matrix,
-                                                            face_to_cell_matrix,
-                                                            face_to_cell_local_face);
-  m_face_to_cell_local_face = face_to_cell_local_face;
+  m_face_to_cell_local_face
+      = m_connectivity_computer.computeLocalItemNumberInChildItem<TypeOfItem::face,
+                                                                  TypeOfItem::cell>(*this);
 
   const auto& face_to_node_matrix
       = m_item_to_item_matrix[itemId(TypeOfItem::face)][itemId(TypeOfItem::node)];
@@ -253,3 +251,50 @@ void Connectivity<2>::_computeFaceCellConnectivities()
     face_to_cell_matrix = face_to_cell_vector;
   }
 }
+
+
+template<size_t Dimension>
+Connectivity<Dimension>::
+Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector)
+{
+  auto& cell_to_node_matrix
+      = m_item_to_item_matrix[itemId(TypeOfItem::cell)][itemId(TypeOfItem::node)];
+  cell_to_node_matrix = cell_by_node_vector;
+
+  Assert(this->numberOfCells()>0);
+
+  {
+    Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", this->numberOfCells());
+    Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const int& j){
+        const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+        inv_cell_nb_nodes[j] = 1./cell_nodes.length;
+      });
+    m_inv_cell_nb_nodes = inv_cell_nb_nodes;
+  }
+
+  auto& node_to_cell_matrix
+      = m_item_to_item_matrix[itemId(TypeOfItem::node)][itemId(TypeOfItem::cell)];
+
+  m_connectivity_computer.computeInverseConnectivityMatrix(cell_to_node_matrix,
+                                                           node_to_cell_matrix);
+
+  m_node_to_cell_local_node
+      = m_connectivity_computer.computeLocalItemNumberInChildItem<TypeOfItem::node, TypeOfItem::cell>(*this);
+
+  m_cell_to_node_local_cell
+      = m_connectivity_computer.computeLocalItemNumberInChildItem<TypeOfItem::cell, TypeOfItem::node>(*this);
+
+  if constexpr (Dimension>1) {
+    this->_computeFaceCellConnectivities();
+  }
+}
+
+
+template Connectivity1D::
+Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector);
+
+template Connectivity2D::
+Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector);
+
+template Connectivity3D::
+Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector);
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index 566ce48505c23b895466acb98b173098230fde3b..dc30bf95076f8d9fe685e9bc82b7c215860f2356 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -6,9 +6,11 @@
 
 #include <Kokkos_Core.hpp>
 
+#include <IConnectivity.hpp>
+
 #include <ConnectivityMatrix.hpp>
-#include <ConnectivityComputer.hpp>
 #include <SubItemValuePerItem.hpp>
+#include <ConnectivityComputer.hpp>
 
 #include <vector>
 #include <unordered_map>
@@ -219,7 +221,8 @@ class ConnectivityFace<3>
 };
 
 template <size_t Dimension>
-class Connectivity
+class Connectivity final
+    : public IConnectivity
 {
  private:
   constexpr static auto& itemId = ItemId<Dimension>::itemId;
@@ -233,11 +236,19 @@ class Connectivity
     return m_item_to_item_matrix[itemId(item_type_0)][itemId(item_type_1)];
   };
 
+  const ConnectivityMatrix& getMatrix(const TypeOfItem& item_type_0,
+                                      const TypeOfItem& item_type_1) const
+  {
+    return m_item_to_item_matrix[itemId(item_type_0)][itemId(item_type_1)];
+  }
+
   const auto& cellFaceIsReversed() const
   {
+    static_assert(dimension>1, "reversed faces makes no sense in dimension 1");
     return m_cell_face_is_reversed;
   }
 
+  [[deprecated("These member functions are poorly named")]]
   const auto& cellToNodeLocalCell() const
   {
     return m_cell_to_node_local_cell;
@@ -268,6 +279,7 @@ private:
 
   FaceValuePerCell<const bool> m_cell_face_is_reversed;
 
+  [[deprecated("These members are poorly named")]]
   NodeValuePerCell<const unsigned short> m_cell_to_node_local_cell;
 
   CellValuePerFace<const unsigned short> m_face_to_cell_local_face;
@@ -370,45 +382,7 @@ private:
 
   Connectivity(const Connectivity&) = delete;
 
-  Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector)
-  {
-    auto& cell_to_node_matrix
-        = m_item_to_item_matrix[itemId(TypeOfItem::cell)][itemId(TypeOfItem::node)];
-    cell_to_node_matrix = cell_by_node_vector;
-
-    Assert(this->numberOfCells()>0);
-
-    {
-      Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", this->numberOfCells());
-      Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const int& j){
-          const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
-          inv_cell_nb_nodes[j] = 1./cell_nodes.length;
-        });
-      m_inv_cell_nb_nodes = inv_cell_nb_nodes;
-    }
-
-    auto& node_to_cell_matrix
-        = m_item_to_item_matrix[itemId(TypeOfItem::node)][itemId(TypeOfItem::cell)];
-
-    m_connectivity_computer.computeInverseConnectivityMatrix(cell_to_node_matrix,
-                                                             node_to_cell_matrix);
-
-    CellValuePerNode<unsigned short> node_to_cell_local_node(*this);
-    m_connectivity_computer.computeLocalChildItemNumberInItem(cell_to_node_matrix,
-                                                              node_to_cell_matrix,
-                                                              node_to_cell_local_node);
-    m_node_to_cell_local_node = node_to_cell_local_node;
-
-    NodeValuePerCell<unsigned short> cell_to_node_local_cell(*this);
-    m_connectivity_computer.computeLocalChildItemNumberInItem(node_to_cell_matrix,
-                                                              cell_to_node_matrix,
-                                                              cell_to_node_local_cell);
-    m_cell_to_node_local_cell = cell_to_node_local_cell;
-
-    if constexpr (Dimension>1) {
-      this->_computeFaceCellConnectivities();
-    }
-  }
+  Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector);
 
   ~Connectivity()
   {
diff --git a/src/mesh/ConnectivityComputer.cpp b/src/mesh/ConnectivityComputer.cpp
index 691df7776db83e44a1bf7f18b91bd84f6627cf18..fea44f9e4ee6b8067b922d2d581c3dfc2ae48773 100644
--- a/src/mesh/ConnectivityComputer.cpp
+++ b/src/mesh/ConnectivityComputer.cpp
@@ -1,3 +1,6 @@
+#include <Connectivity.hpp>
+#include <ConnectivityUtils.hpp>
+
 #include <ConnectivityComputer.hpp>
 #include <map>
 
@@ -33,11 +36,19 @@ computeInverseConnectivityMatrix(const ConnectivityMatrix& item_to_child_item_ma
   child_item_to_item_matrix = child_item_to_items_vector;
 }
 
-void ConnectivityComputer::
-computeLocalChildItemNumberInItem(const ConnectivityMatrix& item_to_child_items_matrix,
-                                  const ConnectivityMatrix& child_item_to_items_matrix,
-                                  CellValuePerNode<unsigned short>& child_item_number_in_item_matrix) const
+template <TypeOfItem child_item_type,
+          TypeOfItem item_type,
+          typename ConnectivityType>
+SubItemValuePerItem<const unsigned short, item_type, child_item_type>
+ConnectivityComputer::computeLocalItemNumberInChildItem(const ConnectivityType& connectivity) const
 {
+#warning check that these connectivities are named correctly
+  const ConnectivityMatrix& child_item_to_items_matrix
+      = getConnectivityMatrix<child_item_type, item_type>(connectivity);
+  const ConnectivityMatrix& item_to_child_items_matrix
+      = getConnectivityMatrix<item_type, child_item_type>(connectivity);
+
+  SubItemValuePerItem<unsigned short, item_type, child_item_type> item_number_in_child_item(connectivity);
   for (unsigned int r=0; r<child_item_to_items_matrix.numRows(); ++r) {
     const auto& child_item_to_items = child_item_to_items_matrix.rowConst(r);
     for (unsigned short J=0; J<child_item_to_items.length; ++J) {
@@ -46,54 +57,64 @@ computeLocalChildItemNumberInItem(const ConnectivityMatrix& item_to_child_items_
 
       for (unsigned int R=0; R<item_to_child_items.length; ++R) {
         if (item_to_child_items(R) == r) {
-          child_item_number_in_item_matrix(r,J) = R;
+          item_number_in_child_item(r,J) = R;
           break;
         }
       }
     }
   }
+
+  return item_number_in_child_item;
 }
 
+// 1D
 
-void ConnectivityComputer::
-computeLocalChildItemNumberInItem(const ConnectivityMatrix& item_to_child_items_matrix,
-                                  const ConnectivityMatrix& child_item_to_items_matrix,
-                                  NodeValuePerCell<unsigned short>& child_item_number_in_item_matrix) const
-{
-  for (unsigned int r=0; r<child_item_to_items_matrix.numRows(); ++r) {
-    const auto& child_item_to_items = child_item_to_items_matrix.rowConst(r);
-    for (unsigned short J=0; J<child_item_to_items.length; ++J) {
-      const unsigned int j = child_item_to_items(J);
-      const auto& item_to_child_items = item_to_child_items_matrix.rowConst(j);
+template SubItemValuePerItem<const unsigned short, TypeOfItem::cell, TypeOfItem::node>
+ConnectivityComputer::
+computeLocalItemNumberInChildItem<TypeOfItem::node,
+                                  TypeOfItem::cell>(const Connectivity1D& connectivity) const;
 
-      for (unsigned int R=0; R<item_to_child_items.length; ++R) {
-        if (item_to_child_items(R) == r) {
-          child_item_number_in_item_matrix(r,J) = R;
-          break;
-        }
-      }
-    }
-  }
-}
+template SubItemValuePerItem<const unsigned short, TypeOfItem::cell, TypeOfItem::face>
+ConnectivityComputer::
+computeLocalItemNumberInChildItem<TypeOfItem::face,
+                                  TypeOfItem::cell>(const Connectivity1D& connectivity) const;
 
+template SubItemValuePerItem<const unsigned short, TypeOfItem::node, TypeOfItem::cell>
+ConnectivityComputer::
+computeLocalItemNumberInChildItem<TypeOfItem::cell,
+                                  TypeOfItem::node>(const Connectivity1D& connectivity) const;
 
-void ConnectivityComputer::
-computeLocalChildItemNumberInItem(const ConnectivityMatrix& item_to_child_items_matrix,
-                                  const ConnectivityMatrix& child_item_to_items_matrix,
-                                  CellValuePerFace<unsigned short>& child_item_number_in_item_matrix) const
-{
-  for (unsigned int r=0; r<child_item_to_items_matrix.numRows(); ++r) {
-    const auto& child_item_to_items = child_item_to_items_matrix.rowConst(r);
-    for (unsigned short J=0; J<child_item_to_items.length; ++J) {
-      const unsigned int j = child_item_to_items(J);
-      const auto& item_to_child_items = item_to_child_items_matrix.rowConst(j);
+// 2D
 
-      for (unsigned int R=0; R<item_to_child_items.length; ++R) {
-        if (item_to_child_items(R) == r) {
-          child_item_number_in_item_matrix(r,J) = R;
-          break;
-        }
-      }
-    }
-  }
-}
+template SubItemValuePerItem<const unsigned short, TypeOfItem::cell, TypeOfItem::node>
+ConnectivityComputer::
+computeLocalItemNumberInChildItem<TypeOfItem::node,
+                                  TypeOfItem::cell>(const Connectivity2D& connectivity) const;
+
+template SubItemValuePerItem<const unsigned short, TypeOfItem::cell, TypeOfItem::face>
+ConnectivityComputer::
+computeLocalItemNumberInChildItem<TypeOfItem::face,
+                                  TypeOfItem::cell>(const Connectivity2D& connectivity) const;
+
+template SubItemValuePerItem<const unsigned short, TypeOfItem::node, TypeOfItem::cell>
+ConnectivityComputer::
+computeLocalItemNumberInChildItem<TypeOfItem::cell,
+                                  TypeOfItem::node>(const Connectivity2D& connectivity) const;
+
+// 3D
+
+template SubItemValuePerItem<const unsigned short, TypeOfItem::cell, TypeOfItem::node>
+ConnectivityComputer::
+computeLocalItemNumberInChildItem<TypeOfItem::node,
+                                  TypeOfItem::cell>(const Connectivity3D& connectivity) const;
+
+
+template SubItemValuePerItem<const unsigned short, TypeOfItem::cell, TypeOfItem::face>
+ConnectivityComputer::
+computeLocalItemNumberInChildItem<TypeOfItem::face,
+                                  TypeOfItem::cell>(const Connectivity3D& connectivity) const;
+
+template SubItemValuePerItem<const unsigned short, TypeOfItem::node, TypeOfItem::cell>
+ConnectivityComputer::
+computeLocalItemNumberInChildItem<TypeOfItem::cell,
+                                  TypeOfItem::node>(const Connectivity3D& connectivity) const;
diff --git a/src/mesh/ConnectivityComputer.hpp b/src/mesh/ConnectivityComputer.hpp
index f8e2e31bbf5d66664aff95d21499a9d33fef40b9..87dcc06727d034114879e15d8de2b3bcf0ec6560 100644
--- a/src/mesh/ConnectivityComputer.hpp
+++ b/src/mesh/ConnectivityComputer.hpp
@@ -9,17 +9,11 @@ struct ConnectivityComputer
   void computeInverseConnectivityMatrix(const ConnectivityMatrix& item_to_child_item_matrix,
                                         ConnectivityMatrix& child_item_to_item_matrix) const;
 
-  void computeLocalChildItemNumberInItem(const ConnectivityMatrix& item_to_child_items_matrix,
-                                         const ConnectivityMatrix& child_item_to_items_matrix,
-                                         CellValuePerNode<unsigned short>&  child_item_number_in_item_matrix) const;
-
-  void computeLocalChildItemNumberInItem(const ConnectivityMatrix& item_to_child_items_matrix,
-                                         const ConnectivityMatrix& child_item_to_items_matrix,
-                                         NodeValuePerCell<unsigned short>&  child_item_number_in_item_matrix) const;
-
-  void computeLocalChildItemNumberInItem(const ConnectivityMatrix& item_to_child_items_matrix,
-                                         const ConnectivityMatrix& child_item_to_items_matrix,
-                                         CellValuePerFace<unsigned short>&  child_item_number_in_item_matrix) const;
+  template <TypeOfItem child_item_type,
+            TypeOfItem item_type,
+            typename ConnectivityType>
+  SubItemValuePerItem<const unsigned short, item_type, child_item_type>
+  computeLocalItemNumberInChildItem(const ConnectivityType& c) const;
 };
 
 #endif // CONNECTIVITY_COMPUTER_HPP
diff --git a/src/mesh/IConnectivity.hpp b/src/mesh/IConnectivity.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..40843d70a75b8264c652a48cabba644351c962b9
--- /dev/null
+++ b/src/mesh/IConnectivity.hpp
@@ -0,0 +1,18 @@
+#ifndef ICONNECTIVITY_HPP
+#define ICONNECTIVITY_HPP
+
+#include <TypeOfItem.hpp>
+#include <ConnectivityMatrix.hpp>
+
+struct IConnectivity
+{
+  virtual const ConnectivityMatrix&
+  getMatrix(const TypeOfItem& item_type_0,
+            const TypeOfItem& item_type_1) const = 0;
+
+  IConnectivity() = default;
+  IConnectivity(const IConnectivity&) = delete;
+  ~IConnectivity() = default;
+};
+
+#endif // ICONNECTIVITY_HPP
diff --git a/src/mesh/SubItemValuePerItem.hpp b/src/mesh/SubItemValuePerItem.hpp
index 839131b9ab28bd0654f7a1da10a1c1b637592356..c578306284185c7c481f2f825f19e74838a9887e 100644
--- a/src/mesh/SubItemValuePerItem.hpp
+++ b/src/mesh/SubItemValuePerItem.hpp
@@ -6,6 +6,7 @@
 #include <PastisAssert.hpp>
 
 #include <ConnectivityMatrix.hpp>
+#include <IConnectivity.hpp>
 #include <ConnectivityUtils.hpp>
 
 template <typename DataType,
@@ -129,26 +130,35 @@ class SubItemValuePerItem
   }
 
   template <typename DataType2>
+  KOKKOS_INLINE_FUNCTION
   SubItemValuePerItem&
   operator=(const SubItemValuePerItem<DataType2, SubItemType, ItemType>& sub_item_value_per_item)
   {
-    // ensures that DataType is the same as source DataType2 or that it is
-    // decorated by a const
-    static_assert(std::is_same<DataType, DataType2>() or std::is_same<DataType, std::add_const_t<DataType2>>(),
-                  "Cannot assign const view to a non const view");
+    // ensures that DataType is the same as source DataType2
+    static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
+                  "Cannot assign SubItemValuePerItem of different type");
+    // ensures that const is not lost through copy
+    static_assert(((std::is_const<DataType2>() and std::is_const<DataType>())
+                   or not std::is_const<DataType2>()),
+                  "Cannot assign const  SubItemValuePerItem to a non const SubItemValuePerItem");
     m_host_row_map = sub_item_value_per_item.m_host_row_map;
     m_values = sub_item_value_per_item.m_values;
     return *this;
   }
 
+  template <typename DataType2>
+  KOKKOS_INLINE_FUNCTION
+  SubItemValuePerItem(const SubItemValuePerItem<DataType2, SubItemType, ItemType>& sub_item_value_per_item)
+  {
+    this->operator=(sub_item_value_per_item);
+  }
 
   SubItemValuePerItem() = default;
 
-  template <typename ConnectivityType>
-  SubItemValuePerItem(const ConnectivityType& connectivity)
+  SubItemValuePerItem(const IConnectivity& connectivity)
   {
     ConnectivityMatrix connectivity_matrix
-        = getConnectivityMatrix<ItemType,SubItemType>(connectivity);
+        = connectivity.getMatrix(ItemType,SubItemType);
 
     m_host_row_map = connectivity_matrix.rowsMap();
     m_values = Kokkos::View<DataType*>("values", connectivity_matrix.numEntries());