diff --git a/src/main.cpp b/src/main.cpp
index a2e5222d87c245ead518eb9709cbfa37e9344be9..34ad22fea1716fa49262894add80484df6662583 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -12,7 +12,6 @@
 // #include <AcousticSolverTest.hpp>
 
 #include <Connectivity1D.hpp>
-#include <Connectivity2D.hpp>
 #include <Connectivity3D.hpp>
 
 #include <Mesh.hpp>
diff --git a/src/mesh/Connectivity2D.hpp b/src/mesh/Connectivity2D.hpp
deleted file mode 100644
index b7f2412ffd87dba5f78ca2338f4ec6cf8a14963c..0000000000000000000000000000000000000000
--- a/src/mesh/Connectivity2D.hpp
+++ /dev/null
@@ -1,274 +0,0 @@
-#ifndef CONNECTIVITY_2D_HPP
-#define CONNECTIVITY_2D_HPP
-
-#include <Kokkos_Core.hpp>
-#include <PastisAssert.hpp>
-#include <TinyVector.hpp>
-
-#include <ConnectivityUtils.hpp>
-#include <vector>
-#include <map>
-#include <algorithm>
-
-#include <TypeOfItem.hpp>
-#include <RefId.hpp>
-#include <RefNodeList.hpp>
-#include <RefFaceList.hpp>
-
-class Connectivity2D
-{
- public:
-  static constexpr size_t dimension = 2;
-
-  ConnectivityMatrix m_cell_to_node_matrix;
-
-  ConnectivityMatrix m_face_to_cell_matrix;
-  ConnectivityMatrix m_face_to_node_matrix;
-
-  ConnectivityMatrix m_node_to_cell_matrix;
-  ConnectivityMatrixShort m_node_to_cell_local_node_matrix;
-
-  // Stores numbering of nodes of each cell.
-  //
-  // This is different from m_cell_to_node_matrix which return the global id of
-  // a local node in a cell
-  ConnectivityMatrix m_node_id_per_cell_matrix;
-
-  inline ConnectivityMatrix subItemIdPerItemMatrix() const
-  {
-    return m_node_id_per_cell_matrix;
-  }
-
- private:
-  std::vector<RefFaceList> m_ref_face_list;
-  std::vector<RefNodeList> m_ref_node_list;
-
-  Kokkos::View<double*> m_inv_cell_nb_nodes;
-
-  Kokkos::View<unsigned short**> m_face_cell_local_face;
-  Kokkos::View<unsigned short**> m_face_node_local_face;
-
-  struct Face
-  {
-    const unsigned int m_node0_id;
-    const unsigned int m_node1_id;
-
-    KOKKOS_INLINE_FUNCTION
-    bool operator<(const Face& f) const
-    {
-      return ((m_node0_id<f.m_node0_id) or
-              ((m_node0_id == f.m_node0_id) and
-               (m_node1_id<f.m_node1_id)));
-    }
-
-    KOKKOS_INLINE_FUNCTION
-    Face& operator=(const Face&) = default;
-
-    KOKKOS_INLINE_FUNCTION
-    Face& operator=(Face&&) = default;
-
-    KOKKOS_INLINE_FUNCTION
-    Face(const Face&) = default;
-
-    KOKKOS_INLINE_FUNCTION
-    Face(Face&&) = default;
-
-    KOKKOS_INLINE_FUNCTION
-    Face(unsigned int node0_id,
-         unsigned int node1_id)
-        : m_node0_id(node0_id),
-          m_node1_id(node1_id)
-    {
-      ;
-    }
-
-    KOKKOS_INLINE_FUNCTION
-    ~Face() = default;
-  };
-
-  void _computeFaceCellConnectivities()
-  {
-    // In 2D faces are simply define
-    typedef std::pair<unsigned int, unsigned short> CellFaceId;
-    std::map<Face, std::vector<CellFaceId>> face_cells_map;
-    for (unsigned int j=0; j<this->numberOfCells(); ++j) {
-      const auto& cell_nodes = m_cell_to_node_matrix.rowConst(j);
-      for (unsigned short r=0; r<cell_nodes.length; ++r) {
-        unsigned int node0_id = cell_nodes(r);
-        unsigned int node1_id = cell_nodes((r+1)%cell_nodes.length);
-        if (node1_id<node0_id) {
-          std::swap(node0_id, node1_id);
-        }
-        face_cells_map[Face(node0_id, node1_id)].push_back(std::make_pair(j, r));
-      }
-    }
-
-    {
-      std::vector<std::vector<unsigned int>> face_to_node_vector(face_cells_map.size());
-      int l=0;
-      for (const auto& face_info : face_cells_map) {
-        const Face& face = face_info.first;
-        face_to_node_vector[l] = {face.m_node0_id, face.m_node1_id};
-        ++l;
-      }
-      m_face_to_node_matrix
-          = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("face_to_node_matrix", face_to_node_vector);
-    }
-
-    {
-      std::vector<std::vector<unsigned int>> face_to_cell_vector(face_cells_map.size());
-      int l=0;
-      for (const auto& face_cells_vector : face_cells_map) {
-        const auto& [face, cell_info_vector] = face_cells_vector;
-        for (const auto& cell_info : cell_info_vector) {
-          face_to_cell_vector[l].push_back(cell_info.second);
-        }
-        ++l;
-      }
-      m_face_to_cell_matrix
-          = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("face_to_cell_matrix", face_to_cell_vector);
-    }
-
-    // Kokkos::View<unsigned short**> face_cell_local_face("face_cell_local_face",
-    //                                                     face_cells_map.size(), m_max_nb_node_per_cell);
-    // {
-    //   int l=0;
-    //   for (const auto& face_cells_vector : face_cells_map) {
-    //     const auto& cells_vector = face_cells_vector.second;
-    //     for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
-    //       unsigned short cell_local_face = cells_vector[lj].second;
-    //       face_cell_local_face(l,lj) = cell_local_face;
-    //     }
-    //     ++l;
-    //   }
-    // }
-    // m_face_cell_local_face = face_cell_local_face;
-  }
-
- public:
-  void addRefFaceList(const RefFaceList& ref_face_list)
-  {
-    m_ref_face_list.push_back(ref_face_list);
-  }
-
-  size_t numberOfRefFaceList() const
-  {
-    return m_ref_face_list.size();
-  }
-
-  const RefFaceList& refFaceList(const size_t& i) const
-  {
-    return m_ref_face_list[i];
-  }
-
-  void addRefNodeList(const RefNodeList& ref_node_list)
-  {
-    m_ref_node_list.push_back(ref_node_list);
-  }
-
-  size_t numberOfRefNodeList() const
-  {
-    return m_ref_node_list.size();
-  }
-
-  const RefNodeList& refNodeList(const size_t& i) const
-  {
-    return m_ref_node_list[i];
-  }
-
-  KOKKOS_INLINE_FUNCTION
-  size_t numberOfNodes() const
-  {
-    return m_node_to_cell_matrix.numRows();
-  }
-
-  KOKKOS_INLINE_FUNCTION
-  size_t numberOfFaces() const
-  {
-    return m_face_to_cell_matrix.numRows();
-  }
-
-  KOKKOS_INLINE_FUNCTION
-  size_t numberOfCells() const
-  {
-    return m_cell_to_node_matrix.numRows();
-  }
-
-  const Kokkos::View<const double*> invCellNbNodes() const
-  {
-    return m_inv_cell_nb_nodes;
-  }
-
-  // const Kokkos::View<const unsigned short**> faceCellLocalFace() const
-  // {
-  //   return m_face_cell_local_face;
-  // }
-
-  unsigned int getFaceNumber(const unsigned int node0_id,
-                             const unsigned int node1_id) const
-  {
-#warning rewrite
-    const unsigned int n0_id = std::min(node0_id, node1_id);
-    const unsigned int n1_id = std::max(node0_id, node1_id);
-    // worst code ever
-    for (unsigned int l=0; l<this->numberOfFaces(); ++l) {
-      const auto& face_nodes = m_face_to_node_matrix.rowConst(l);
-      if ((face_nodes(0) == n0_id) and (face_nodes(1) == n1_id)) {
-        return l;
-      }
-    }
-    std::cerr << "Face not found!\n";
-    std::exit(0);
-    return 0;
-  }
-
-  Connectivity2D(const Connectivity2D&) = delete;
-
-  Connectivity2D(std::vector<std::vector<unsigned int>> cell_by_node_vector)
-  {
-    Assert(cell_by_node_vector.size()>0);
-
-    m_cell_to_node_matrix
-        = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("cell_to_node_matrix",
-                                                            cell_by_node_vector);
-
-    {
-      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 = m_cell_to_node_matrix.rowConst(j);
-          inv_cell_nb_nodes[j] = 1./cell_nodes.length;
-        });
-      m_inv_cell_nb_nodes = inv_cell_nb_nodes;
-    }
-
-    {
-      std::vector<std::vector<unsigned int>> node_id_per_cell_vector(this->numberOfCells());
-      unsigned int id=0;
-      for (unsigned int j=0; j<this->numberOfCells(); ++j) {
-        const auto& cell_to_node = m_cell_to_node_matrix.rowConst(j);
-        auto& node_id_per_cell = node_id_per_cell_vector[j];
-        node_id_per_cell.resize(cell_to_node.length);
-        for (size_t r=0; r<cell_to_node.length; ++r) {
-          node_id_per_cell[r] = id++;
-        }
-      }
-      m_node_id_per_cell_matrix
-          = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("node_id_per_cell_matrix",
-                                                              cell_by_node_vector);
-    }
-
-    ConnectivityUtils utils;
-    utils.computeNodeCellConnectivity(m_cell_to_node_matrix,
-                                      m_node_to_cell_matrix,
-                                      m_node_to_cell_local_node_matrix);
-
-    this->_computeFaceCellConnectivities();
-  }
-
-  ~Connectivity2D()
-  {
-    ;
-  }
-};
-
-#endif // CONNECTIVITY_2D_HPP
diff --git a/src/mesh/Connectivity3D.hpp b/src/mesh/Connectivity3D.hpp
index 0eb1c631898a37582d44e4634f7df6739fba0f7c..10f8f2d1b6b74fb284573bc280dd9176a4c462ec 100644
--- a/src/mesh/Connectivity3D.hpp
+++ b/src/mesh/Connectivity3D.hpp
@@ -18,369 +18,245 @@
 
 #include <tuple>
 
-class Connectivity3D
-{
-public:
-  static constexpr size_t dimension = 3;
+template <size_t Dimension>
+class Connectivity;
 
-  ConnectivityMatrix m_cell_to_node_matrix;
+template <size_t Dimension>
+class ConnectivityFace;
 
-  ConnectivityMatrix m_cell_to_face_matrix;
-  ConnectivityMatrixShort m_cell_to_face_is_reversed_matrix;
+template<>
+class ConnectivityFace<2>
+{
+ public:
+  friend struct Hash;
+  struct Hash
+  {
+    size_t operator()(const ConnectivityFace& f) const {
+      size_t hash = 0;
+      hash ^= std::hash<unsigned int>()(f.m_node0_id);
+      hash ^= std::hash<unsigned int>()(f.m_node1_id) >> 1;
+      return hash;
+    }
+  };
 
-  ConnectivityMatrix m_face_to_cell_matrix;
-  ConnectivityMatrixShort m_face_to_cell_local_face_matrix;
-  ConnectivityMatrix m_face_to_node_matrix;
+  unsigned int m_node0_id;
+  unsigned int m_node1_id;
 
-  ConnectivityMatrix m_node_to_cell_matrix;
-  ConnectivityMatrixShort m_node_to_cell_local_node_matrix;
 
-  // Stores numbering of nodes of each cell.
-  // gives an id to each node of each cell. (j,r) -> id
-  //
-  // This is different from m_cell_to_node_matrix which return the global id of
-  // a local node in a cell
-  ConnectivityMatrix m_node_id_per_cell_matrix;
+  friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f)
+  {
+    os << f.m_node0_id << ' ' << f.m_node1_id << ' ';
+    return os;
+  }
 
-  inline ConnectivityMatrix subItemIdPerItemMatrix() const
+  KOKKOS_INLINE_FUNCTION
+  bool operator==(const ConnectivityFace& f) const
   {
-    return m_node_id_per_cell_matrix;
+    return ((m_node0_id == f.m_node0_id) and
+            (m_node1_id == f.m_node1_id));
   }
 
-private:
-  std::vector<RefFaceList> m_ref_face_list;
-  std::vector<RefNodeList> m_ref_node_list;
+  KOKKOS_INLINE_FUNCTION
+  bool operator<(const ConnectivityFace& f) const
+  {
+    return ((m_node0_id<f.m_node0_id) or
+            ((m_node0_id == f.m_node0_id) and
+             (m_node1_id<f.m_node1_id)));
+  }
 
-  Kokkos::View<double*> m_inv_cell_nb_nodes;
+  KOKKOS_INLINE_FUNCTION
+  ConnectivityFace& operator=(const ConnectivityFace&) = default;
 
-  Kokkos::View<const unsigned short*> m_node_nb_faces;
-  Kokkos::View<const unsigned int**> m_node_faces;
+  KOKKOS_INLINE_FUNCTION
+  ConnectivityFace& operator=(ConnectivityFace&&) = default;
 
-  class Face
+  KOKKOS_INLINE_FUNCTION
+  ConnectivityFace(const std::vector<unsigned int>& given_node_id_list)
   {
-   public:
-    friend struct Hash;
-    struct Hash {
-      size_t operator()(const Face& f) const
-      {
-        size_t hash = 0;
-        for (size_t i=0; i<f.m_node_id_list.size(); ++i) {
-          hash ^= std::hash<unsigned int>()(f.m_node_id_list[i]) >> i;
-        }
-        return hash;
-      }
-    };
-   private:
-    bool m_reversed;
-    std::vector<unsigned int> m_node_id_list;
+    Assert(given_node_id_list.size()==2);
+#warning rework this dirty constructor
+    m_node0_id = std::min(given_node_id_list[0], given_node_id_list[1]);
+    m_node1_id = std::max(given_node_id_list[0], given_node_id_list[1]);
+  }
 
-   public:
-    friend std::ostream& operator<<(std::ostream& os, const Face& f)
-    {
-      for (auto id : f.m_node_id_list) {
-        std::cout << id << ' ';
+  KOKKOS_INLINE_FUNCTION
+  ConnectivityFace(const ConnectivityFace&) = default;
+
+  KOKKOS_INLINE_FUNCTION
+  ConnectivityFace(ConnectivityFace&&) = default;
+
+  // KOKKOS_INLINE_FUNCTION
+  // ConnectivityFace(unsigned int node0_id,
+  //                  unsigned int node1_id)
+  //     : m_node0_id(node0_id),
+  //       m_node1_id(node1_id)
+  // {
+  //   ;
+  // }
+
+  KOKKOS_INLINE_FUNCTION
+  ~ConnectivityFace() = default;
+};
+
+template <>
+class ConnectivityFace<3>
+{
+ private:
+  friend class Connectivity<3>;
+  friend struct Hash;
+  struct Hash
+  {
+    size_t operator()(const ConnectivityFace& f) const {
+      size_t hash = 0;
+      for (size_t i=0; i<f.m_node_id_list.size(); ++i) {
+        hash ^= std::hash<unsigned int>()(f.m_node_id_list[i]) >> i;
       }
-      return os;
+      return hash;
     }
+  };
 
-    KOKKOS_INLINE_FUNCTION
-    const bool& reversed() const
-    {
-      return m_reversed;
-    }
+  bool m_reversed;
+  std::vector<unsigned int> m_node_id_list;
 
-    KOKKOS_INLINE_FUNCTION
-    const std::vector<unsigned int>& nodeIdList() const
-    {
-      return m_node_id_list;
+  friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f)
+  {
+    for (auto id : f.m_node_id_list) {
+      std::cout << id << ' ';
     }
+    return os;
+  }
 
-    KOKKOS_INLINE_FUNCTION
-    std::vector<unsigned int> _sort(const std::vector<unsigned int>& node_list)
-    {
-      const auto min_id = std::min_element(node_list.begin(), node_list.end());
-      const int shift = std::distance(node_list.begin(), min_id);
-
-      std::vector<unsigned int> rotated_node_list(node_list.size());
-      if (node_list[(shift+1)%node_list.size()] > node_list[(shift+node_list.size()-1)%node_list.size()]) {
-        for (size_t i=0; i<node_list.size(); ++i) {
-          rotated_node_list[i] = node_list[(shift+node_list.size()-i)%node_list.size()];
-          m_reversed = true;
-        }
-      } else {
-        for (size_t i=0; i<node_list.size(); ++i) {
-          rotated_node_list[i] = node_list[(shift+i)%node_list.size()];
-        }
-      }
+  KOKKOS_INLINE_FUNCTION
+  const bool& reversed() const
+  {
+    return m_reversed;
+  }
 
-      return rotated_node_list;
+  KOKKOS_INLINE_FUNCTION
+  const std::vector<unsigned int>& nodeIdList() const
+  {
+    return m_node_id_list;
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  std::vector<unsigned int> _sort(const std::vector<unsigned int>& node_list)
+  {
+    const auto min_id = std::min_element(node_list.begin(), node_list.end());
+    const int shift = std::distance(node_list.begin(), min_id);
+
+    std::vector<unsigned int> rotated_node_list(node_list.size());
+    if (node_list[(shift+1)%node_list.size()] > node_list[(shift+node_list.size()-1)%node_list.size()]) {
+      for (size_t i=0; i<node_list.size(); ++i) {
+        rotated_node_list[i] = node_list[(shift+node_list.size()-i)%node_list.size()];
+        m_reversed = true;
+      }
+    } else {
+      for (size_t i=0; i<node_list.size(); ++i) {
+        rotated_node_list[i] = node_list[(shift+i)%node_list.size()];
+      }
     }
 
-    bool operator==(const Face& f) const
-    {
-      if (m_node_id_list.size() == f.nodeIdList().size()) {
-        for (size_t j=0; j<m_node_id_list.size(); ++j) {
-          if (m_node_id_list[j] != f.nodeIdList()[j]) {
-            return false;
-          }
+    return rotated_node_list;
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  ConnectivityFace(const std::vector<unsigned int>& given_node_id_list)
+      : m_reversed(false),
+        m_node_id_list(_sort(given_node_id_list))
+  {
+    ;
+  }
+
+ public:
+  bool operator==(const ConnectivityFace& f) const
+  {
+    if (m_node_id_list.size() == f.nodeIdList().size()) {
+      for (size_t j=0; j<m_node_id_list.size(); ++j) {
+        if (m_node_id_list[j] != f.nodeIdList()[j]) {
+          return false;
         }
-        return true;
       }
-      return false;
+      return true;
     }
+    return false;
+  }
 
-    KOKKOS_INLINE_FUNCTION
-    bool operator<(const Face& f) const
-    {
-      const size_t min_nb_nodes = std::min(f.m_node_id_list.size(), m_node_id_list.size());
-      for (size_t i=0; i<min_nb_nodes; ++i) {
-        if (m_node_id_list[i] <  f.m_node_id_list[i]) return true;
-        if (m_node_id_list[i] != f.m_node_id_list[i]) return false;
-      }
-      return m_node_id_list.size() < f.m_node_id_list.size();
+  KOKKOS_INLINE_FUNCTION
+  bool operator<(const ConnectivityFace& f) const
+  {
+    const size_t min_nb_nodes = std::min(f.m_node_id_list.size(), m_node_id_list.size());
+    for (size_t i=0; i<min_nb_nodes; ++i) {
+      if (m_node_id_list[i] <  f.m_node_id_list[i]) return true;
+      if (m_node_id_list[i] != f.m_node_id_list[i]) return false;
     }
+    return m_node_id_list.size() < f.m_node_id_list.size();
+  }
 
-    KOKKOS_INLINE_FUNCTION
-    Face& operator=(const Face&) = default;
+  KOKKOS_INLINE_FUNCTION
+  ConnectivityFace& operator=(const ConnectivityFace&) = default;
 
-    KOKKOS_INLINE_FUNCTION
-    Face& operator=(Face&&) = default;
+  KOKKOS_INLINE_FUNCTION
+  ConnectivityFace& operator=(ConnectivityFace&&) = default;
 
-    KOKKOS_INLINE_FUNCTION
-    Face(const Face&) = default;
+  KOKKOS_INLINE_FUNCTION
+  ConnectivityFace(const ConnectivityFace&) = default;
 
-    KOKKOS_INLINE_FUNCTION
-    Face(Face&&) = default;
+  KOKKOS_INLINE_FUNCTION
+  ConnectivityFace(ConnectivityFace&&) = default;
 
-    KOKKOS_INLINE_FUNCTION
-    Face(const std::vector<unsigned int>& given_node_id_list)
-        : m_reversed(false),
-          m_node_id_list(_sort(given_node_id_list))
-    {
-      ;
-    }
 
-    KOKKOS_INLINE_FUNCTION
-    Face() = delete;
+  KOKKOS_INLINE_FUNCTION
+  ConnectivityFace() = delete;
 
-    KOKKOS_INLINE_FUNCTION
-    ~Face() = default;
-  };
+  KOKKOS_INLINE_FUNCTION
+  ~ConnectivityFace() = default;
+};
 
-  std::unordered_map<Face, unsigned int, Face::Hash> m_face_number_map;
-
-  void _computeFaceCellConnectivities()
-  {
-    Kokkos::View<unsigned short*> cell_nb_faces("cell_nb_faces", this->numberOfCells());
-
-    typedef std::tuple<unsigned int, unsigned short, bool> CellFaceInfo;
-    std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
-    for (unsigned int j=0; j<this->numberOfCells(); ++j) {
-      const auto& cell_nodes = m_cell_to_node_matrix.rowConst(j);
-
-      switch (cell_nodes.length) {
-        case 4: { // tetrahedron
-          cell_nb_faces[j] = 4;
-          // face 0
-          Face f0({cell_nodes(1),
-                   cell_nodes(2),
-                   cell_nodes(3)});
-          face_cells_map[f0].push_back(std::make_tuple(j, 0, f0.reversed()));
-
-          // face 1
-          Face f1({cell_nodes(0),
-                   cell_nodes(3),
-                   cell_nodes(2)});
-          face_cells_map[f1].push_back(std::make_tuple(j, 1, f1.reversed()));
-
-          // face 2
-          Face f2({cell_nodes(0),
-                   cell_nodes(1),
-                   cell_nodes(3)});
-          face_cells_map[f2].push_back(std::make_tuple(j, 2, f2.reversed()));
-
-          // face 3
-          Face f3({cell_nodes(0),
-                   cell_nodes(2),
-                   cell_nodes(1)});
-          face_cells_map[f3].push_back(std::make_tuple(j, 3, f3.reversed()));
-          break;
-        }
-        case 8: { // hexahedron
-          // face 0
-          Face f0({cell_nodes(3),
-                   cell_nodes(2),
-                   cell_nodes(1),
-                   cell_nodes(0)});
-          face_cells_map[f0].push_back(std::make_tuple(j, 0, f0.reversed()));
-
-          // face 1
-          Face f1({cell_nodes(4),
-                   cell_nodes(5),
-                   cell_nodes(6),
-                   cell_nodes(7)});
-          face_cells_map[f1].push_back(std::make_tuple(j, 1, f1.reversed()));
-
-          // face 2
-          Face f2({cell_nodes(0),
-                   cell_nodes(4),
-                   cell_nodes(7),
-                   cell_nodes(3)});
-          face_cells_map[f2].push_back(std::make_tuple(j, 2, f2.reversed()));
-
-          // face 3
-          Face f3({cell_nodes(1),
-                   cell_nodes(2),
-                   cell_nodes(6),
-                   cell_nodes(5)});
-          face_cells_map[f3].push_back(std::make_tuple(j, 3, f3.reversed()));
-
-          // face 4
-          Face f4({cell_nodes(0),
-                   cell_nodes(1),
-                   cell_nodes(5),
-                   cell_nodes(4)});
-          face_cells_map[f4].push_back(std::make_tuple(j, 4, f4.reversed()));
-
-          // face 5
-          Face f5({cell_nodes(3),
-                   cell_nodes(7),
-                   cell_nodes(6),
-                   cell_nodes(2)});
-          face_cells_map[f5].push_back(std::make_tuple(j, 5, f5.reversed()));
-
-          cell_nb_faces[j] = 6;
-          break;
-        }
-        default: {
-          std::cerr << "unexpected cell type!\n";
-          std::exit(0);
-        }
-      }
-    }
+template <size_t Dimension>
+class Connectivity
+{
+public:
+  static constexpr size_t dimension = Dimension;
+  static constexpr bool has_edges = (dimension>2);
+  static constexpr bool has_face = (dimension>1);
 
-    {
-      std::vector<std::vector<unsigned int>> face_to_node_vector(face_cells_map.size());
-      int l=0;
-      for (const auto& face_info : face_cells_map) {
-        const Face& face = face_info.first;
-        face_to_node_vector[l] = face.nodeIdList();
-        ++l;
-      }
-      m_face_to_node_matrix
-          = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("face_to_node_matrix", face_to_node_vector);
-    }
+  ConnectivityMatrix m_cell_to_node_matrix;
 
-    {
-      int l=0;
-      for (const auto& face_cells_vector : face_cells_map) {
-        const Face& face = face_cells_vector.first;
-        m_face_number_map[face] = l;
-        ++l;
-      }
-    }
+  ConnectivityMatrix m_cell_to_face_matrix;
+  ConnectivityMatrixShort m_cell_to_face_is_reversed_matrix;
 
-    {
-      std::vector<std::vector<unsigned int>> face_to_cell_vector(face_cells_map.size());
-      size_t l=0;
-      for (const auto& face_cells_vector : face_cells_map) {
-        const auto& cells_info_vector = face_cells_vector.second;
-        std::vector<unsigned int>& cells_vector = face_to_cell_vector[l];
-        cells_vector.resize(cells_info_vector.size());
-        for (size_t j=0; j<cells_info_vector.size(); ++j) {
-          const auto& [cell_number, local_face_in_cell, reversed] = cells_info_vector[j];
-          cells_vector[j] = cell_number;
-        }
-        ++l;
-      }
-      m_face_to_cell_matrix
-          = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("face_to_cell_matrix", face_to_cell_vector);
-    }
+  ConnectivityMatrix m_face_to_cell_matrix;
+  ConnectivityMatrixShort m_face_to_cell_local_face_matrix;
+  ConnectivityMatrix m_face_to_node_matrix;
 
-    {
-      std::vector<std::vector<unsigned short>> face_to_cell_local_face_vector(face_cells_map.size());
-      size_t l=0;
-      for (const auto& face_cells_vector : face_cells_map) {
-        const auto& cells_info_vector = face_cells_vector.second;
-        std::vector<unsigned short>& cells_vector = face_to_cell_local_face_vector[l];
-        cells_vector.resize(cells_info_vector.size());
-        for (size_t j=0; j<cells_info_vector.size(); ++j) {
-          const auto& [cell_number, local_face_in_cell, reversed] = cells_info_vector[j];
-          cells_vector[j] = local_face_in_cell;
-        }
-        ++l;
-      }
-      m_face_to_cell_local_face_matrix
-          = Kokkos::create_staticcrsgraph<ConnectivityMatrixShort>("face_to_cell_local_face_matrix",
-                                                                   face_to_cell_local_face_vector);
-    }
+  ConnectivityMatrix m_node_to_cell_matrix;
+  ConnectivityMatrixShort m_node_to_cell_local_node_matrix;
 
-#warning check that the number of cell per faces is <=2
+  template <TypeOfItem SubItemType,
+            TypeOfItem ItemType>
+  const ConnectivityMatrix& subItemIdPerItemMatrix() const = delete;
 
-    {
-      std::vector<std::vector<unsigned int>> cell_to_face_vector(this->numberOfCells());
-      for (size_t j=0; j<cell_to_face_vector.size(); ++j) {
-        cell_to_face_vector[j].resize(cell_nb_faces[j]);
-      }
-      int l=0;
-      for (const auto& face_cells_vector : face_cells_map) {
-        const auto& cells_vector = face_cells_vector.second;
-        for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
-          const auto& [cell_number, cell_local_face, reversed] = cells_vector[lj];
-          cell_to_face_vector[cell_number][cell_local_face] = l;
-        }
-        ++l;
-      }
-      m_cell_to_face_matrix
-          = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("cell_to_face_matrix", cell_to_face_vector);
-    }
+private:
+  // Stores numbering of nodes of each cell.
+  // gives an id to each node of each cell. (j,r) -> id
+  //
+  // This is different from m_cell_to_node_matrix which return the global id of
+  // a local node in a cell
+  ConnectivityMatrix m_node_id_per_cell_matrix;
 
-    {
-      std::vector<std::vector<unsigned short>> cell_to_face_is_reversed_vector(this->numberOfCells());
-      for (size_t j=0; j<cell_to_face_is_reversed_vector.size(); ++j) {
-        cell_to_face_is_reversed_vector[j].resize(cell_nb_faces[j]);
-      }
-      int l=0;
-      for (const auto& face_cells_vector : face_cells_map) {
-        const auto& cells_vector = face_cells_vector.second;
-        for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
-          const auto& [cell_number, cell_local_face, reversed] = cells_vector[lj];
-          cell_to_face_is_reversed_vector[cell_number][cell_local_face] = reversed;
-        }
-        ++l;
-      }
+  std::vector<RefFaceList> m_ref_face_list;
+  std::vector<RefNodeList> m_ref_node_list;
 
-      m_cell_to_face_is_reversed_matrix
-          = Kokkos::create_staticcrsgraph<ConnectivityMatrixShort>("cell_to_face_is_reversed_matrix",
-                                                                   cell_to_face_is_reversed_vector);
-    }
+  Kokkos::View<double*> m_inv_cell_nb_nodes;
 
-    std::unordered_map<unsigned int, std::vector<unsigned int>> node_faces_map;
-    for (size_t l=0; l<m_face_to_node_matrix.numRows(); ++l) {
-      const auto& face_nodes = m_face_to_node_matrix.rowConst(l);
-      for (size_t lr=0; lr<face_nodes.length; ++lr) {
-        const unsigned int r = face_nodes(lr);
-        node_faces_map[r].push_back(l);
-      }
-    }
-    Kokkos::View<unsigned short*> node_nb_faces("node_nb_faces", this->numberOfNodes());
-    size_t max_nb_face_per_node = 0;
-    for (auto node_faces : node_faces_map) {
-      max_nb_face_per_node = std::max(node_faces.second.size(), max_nb_face_per_node);
-      node_nb_faces[node_faces.first] = node_faces.second.size();
-    }
-    m_node_nb_faces = node_nb_faces;
-
-    Kokkos::View<unsigned int**> node_faces("node_faces", this->numberOfNodes(), max_nb_face_per_node);
-    for (auto node_faces_vector : node_faces_map) {
-      const unsigned int r = node_faces_vector.first;
-      const std::vector<unsigned int>&  faces_vector = node_faces_vector.second;
-      for (size_t l=0; l < faces_vector.size(); ++l) {
-        node_faces(r, l) = faces_vector[l];
-      }
-    }
-    m_node_faces = node_faces;
-  }
+  Kokkos::View<const unsigned short*> m_node_nb_faces;
+  Kokkos::View<const unsigned int**> m_node_faces;
+
+  using Face = ConnectivityFace<Dimension>;
+
+  std::unordered_map<Face, unsigned int, typename Face::Hash> m_face_number_map;
+
+  void _computeFaceCellConnectivities();
 
  public:
   void addRefFaceList(const RefFaceList& ref_face_list)
@@ -442,14 +318,15 @@ private:
     auto i_face = m_face_number_map.find(face);
     if (i_face == m_face_number_map.end()) {
       std::cerr << "Face " << face << " not found!\n";
+      throw std::exception();
       std::exit(0);
     }
     return i_face->second;
   }
 
-  Connectivity3D(const Connectivity3D&) = delete;
+  Connectivity(const Connectivity&) = delete;
 
-  Connectivity3D(const std::vector<std::vector<unsigned int>>& cell_by_node_vector)
+  Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector)
   {
     m_cell_to_node_matrix
         = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("cell_to_node_matrix",
@@ -490,10 +367,300 @@ private:
     this->_computeFaceCellConnectivities();
   }
 
-  ~Connectivity3D()
+  ~Connectivity()
   {
     ;
   }
 };
 
+template<>
+inline void Connectivity<3>::_computeFaceCellConnectivities()
+{
+  Kokkos::View<unsigned short*> cell_nb_faces("cell_nb_faces", this->numberOfCells());
+
+  typedef std::tuple<unsigned int, unsigned short, bool> CellFaceInfo;
+  std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
+  for (unsigned int j=0; j<this->numberOfCells(); ++j) {
+    const auto& cell_nodes = m_cell_to_node_matrix.rowConst(j);
+
+    switch (cell_nodes.length) {
+      case 4: { // tetrahedron
+        cell_nb_faces[j] = 4;
+        // face 0
+        Face f0({cell_nodes(1),
+                 cell_nodes(2),
+                 cell_nodes(3)});
+        face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
+
+        // face 1
+        Face f1({cell_nodes(0),
+                 cell_nodes(3),
+                 cell_nodes(2)});
+        face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
+
+        // face 2
+        Face f2({cell_nodes(0),
+                 cell_nodes(1),
+                 cell_nodes(3)});
+        face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
+
+        // face 3
+        Face f3({cell_nodes(0),
+                 cell_nodes(2),
+                 cell_nodes(1)});
+        face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
+        break;
+      }
+      case 8: { // hexahedron
+        // face 0
+        Face f0({cell_nodes(3),
+                 cell_nodes(2),
+                 cell_nodes(1),
+                 cell_nodes(0)});
+        face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
+
+        // face 1
+        Face f1({cell_nodes(4),
+                 cell_nodes(5),
+                 cell_nodes(6),
+                 cell_nodes(7)});
+        face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
+
+        // face 2
+        Face f2({cell_nodes(0),
+                 cell_nodes(4),
+                 cell_nodes(7),
+                 cell_nodes(3)});
+        face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
+
+        // face 3
+        Face f3({cell_nodes(1),
+                 cell_nodes(2),
+                 cell_nodes(6),
+                 cell_nodes(5)});
+        face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
+
+        // face 4
+        Face f4({cell_nodes(0),
+                 cell_nodes(1),
+                 cell_nodes(5),
+                 cell_nodes(4)});
+        face_cells_map[f4].emplace_back(std::make_tuple(j, 4, f4.reversed()));
+
+        // face 5
+        Face f5({cell_nodes(3),
+                 cell_nodes(7),
+                 cell_nodes(6),
+                 cell_nodes(2)});
+        face_cells_map[f5].emplace_back(std::make_tuple(j, 5, f5.reversed()));
+
+        cell_nb_faces[j] = 6;
+        break;
+      }
+      default: {
+        std::cerr << "unexpected cell type!\n";
+        std::exit(0);
+      }
+    }
+  }
+
+  {
+    std::vector<std::vector<unsigned int>> face_to_node_vector(face_cells_map.size());
+    int l=0;
+    for (const auto& face_info : face_cells_map) {
+      const Face& face = face_info.first;
+      face_to_node_vector[l] = face.nodeIdList();
+      ++l;
+    }
+    m_face_to_node_matrix
+        = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("face_to_node_matrix", face_to_node_vector);
+  }
+
+  {
+    int l=0;
+    for (const auto& face_cells_vector : face_cells_map) {
+      const Face& face = face_cells_vector.first;
+      m_face_number_map[face] = l;
+      ++l;
+    }
+  }
+
+  {
+    std::vector<std::vector<unsigned int>> face_to_cell_vector(face_cells_map.size());
+    size_t l=0;
+    for (const auto& face_cells_vector : face_cells_map) {
+      const auto& cells_info_vector = face_cells_vector.second;
+      std::vector<unsigned int>& cells_vector = face_to_cell_vector[l];
+      cells_vector.resize(cells_info_vector.size());
+      for (size_t j=0; j<cells_info_vector.size(); ++j) {
+        const auto& [cell_number, local_face_in_cell, reversed] = cells_info_vector[j];
+        cells_vector[j] = cell_number;
+      }
+      ++l;
+    }
+    m_face_to_cell_matrix
+        = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("face_to_cell_matrix", face_to_cell_vector);
+  }
+
+  {
+    std::vector<std::vector<unsigned short>> face_to_cell_local_face_vector(face_cells_map.size());
+    size_t l=0;
+    for (const auto& face_cells_vector : face_cells_map) {
+      const auto& cells_info_vector = face_cells_vector.second;
+      std::vector<unsigned short>& cells_vector = face_to_cell_local_face_vector[l];
+      cells_vector.resize(cells_info_vector.size());
+      for (size_t j=0; j<cells_info_vector.size(); ++j) {
+        const auto& [cell_number, local_face_in_cell, reversed] = cells_info_vector[j];
+        cells_vector[j] = local_face_in_cell;
+      }
+      ++l;
+    }
+    m_face_to_cell_local_face_matrix
+        = Kokkos::create_staticcrsgraph<ConnectivityMatrixShort>("face_to_cell_local_face_matrix",
+                                                                 face_to_cell_local_face_vector);
+  }
+
+#warning check that the number of cell per faces is <=2
+
+  {
+    std::vector<std::vector<unsigned int>> cell_to_face_vector(this->numberOfCells());
+    for (size_t j=0; j<cell_to_face_vector.size(); ++j) {
+      cell_to_face_vector[j].resize(cell_nb_faces[j]);
+    }
+    int l=0;
+    for (const auto& face_cells_vector : face_cells_map) {
+      const auto& cells_vector = face_cells_vector.second;
+      for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
+        const auto& [cell_number, cell_local_face, reversed] = cells_vector[lj];
+        cell_to_face_vector[cell_number][cell_local_face] = l;
+      }
+      ++l;
+    }
+    m_cell_to_face_matrix
+        = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("cell_to_face_matrix", cell_to_face_vector);
+  }
+
+  {
+    std::vector<std::vector<unsigned short>> cell_to_face_is_reversed_vector(this->numberOfCells());
+    for (size_t j=0; j<cell_to_face_is_reversed_vector.size(); ++j) {
+      cell_to_face_is_reversed_vector[j].resize(cell_nb_faces[j]);
+    }
+    int l=0;
+    for (const auto& face_cells_vector : face_cells_map) {
+      const auto& cells_vector = face_cells_vector.second;
+      for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
+        const auto& [cell_number, cell_local_face, reversed] = cells_vector[lj];
+        cell_to_face_is_reversed_vector[cell_number][cell_local_face] = reversed;
+      }
+      ++l;
+    }
+
+    m_cell_to_face_is_reversed_matrix
+        = Kokkos::create_staticcrsgraph<ConnectivityMatrixShort>("cell_to_face_is_reversed_matrix",
+                                                                 cell_to_face_is_reversed_vector);
+  }
+
+  std::unordered_map<unsigned int, std::vector<unsigned int>> node_faces_map;
+  for (size_t l=0; l<m_face_to_node_matrix.numRows(); ++l) {
+    const auto& face_nodes = m_face_to_node_matrix.rowConst(l);
+    for (size_t lr=0; lr<face_nodes.length; ++lr) {
+      const unsigned int r = face_nodes(lr);
+      node_faces_map[r].emplace_back(l);
+    }
+  }
+  Kokkos::View<unsigned short*> node_nb_faces("node_nb_faces", this->numberOfNodes());
+  size_t max_nb_face_per_node = 0;
+  for (auto node_faces : node_faces_map) {
+    max_nb_face_per_node = std::max(node_faces.second.size(), max_nb_face_per_node);
+    node_nb_faces[node_faces.first] = node_faces.second.size();
+  }
+  m_node_nb_faces = node_nb_faces;
+
+  Kokkos::View<unsigned int**> node_faces("node_faces", this->numberOfNodes(), max_nb_face_per_node);
+  for (auto node_faces_vector : node_faces_map) {
+    const unsigned int r = node_faces_vector.first;
+    const std::vector<unsigned int>&  faces_vector = node_faces_vector.second;
+    for (size_t l=0; l < faces_vector.size(); ++l) {
+      node_faces(r, l) = faces_vector[l];
+    }
+  }
+  m_node_faces = node_faces;
+}
+
+template<>
+inline void Connectivity<2>::_computeFaceCellConnectivities()
+{
+  // In 2D faces are simply define
+  typedef std::pair<unsigned int, unsigned short> CellFaceId;
+  std::map<Face, std::vector<CellFaceId>> face_cells_map;
+  for (unsigned int j=0; j<this->numberOfCells(); ++j) {
+    const auto& cell_nodes = m_cell_to_node_matrix.rowConst(j);
+    for (unsigned short r=0; r<cell_nodes.length; ++r) {
+      unsigned int node0_id = cell_nodes(r);
+      unsigned int node1_id = cell_nodes((r+1)%cell_nodes.length);
+      if (node1_id<node0_id) {
+        std::swap(node0_id, node1_id);
+      }
+      face_cells_map[Face({node0_id, node1_id})].push_back(std::make_pair(j, r));
+    }
+  }
+
+  {
+    int l=0;
+    for (const auto& face_cells_vector : face_cells_map) {
+      const Face& face = face_cells_vector.first;
+      m_face_number_map[face] = l;
+      ++l;
+    }
+  }
+
+  {
+    std::vector<std::vector<unsigned int>> face_to_node_vector(face_cells_map.size());
+    int l=0;
+    for (const auto& face_info : face_cells_map) {
+      const Face& face = face_info.first;
+      face_to_node_vector[l] = {face.m_node0_id, face.m_node1_id};
+      ++l;
+    }
+    m_face_to_node_matrix
+        = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("face_to_node_matrix", face_to_node_vector);
+  }
+
+  {
+    std::vector<std::vector<unsigned int>> face_to_cell_vector(face_cells_map.size());
+    int l=0;
+    for (const auto& face_cells_vector : face_cells_map) {
+      const auto& [face, cell_info_vector] = face_cells_vector;
+      for (const auto& cell_info : cell_info_vector) {
+        face_to_cell_vector[l].push_back(cell_info.second);
+      }
+      ++l;
+    }
+    m_face_to_cell_matrix
+        = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("face_to_cell_matrix", face_to_cell_vector);
+  }
+}
+
+using Connectivity3D = Connectivity<3>;
+
+template <>
+template <>
+inline const ConnectivityMatrix&
+Connectivity<3>::subItemIdPerItemMatrix<TypeOfItem::node,
+                                        TypeOfItem::cell>() const
+{
+  return m_node_id_per_cell_matrix;
+}
+
+using Connectivity2D = Connectivity<2>;
+
+template <>
+template <>
+inline const ConnectivityMatrix&
+Connectivity<2>::subItemIdPerItemMatrix<TypeOfItem::node,
+                                        TypeOfItem::cell>() const
+{
+  return m_node_id_per_cell_matrix;
+}
+
 #endif // CONNECTIVITY_3D_HPP
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 139657d1f325ef1fb39b49f33c8835539ee8e0b1..6e422a5e8381336623a8fdc3a891c33dcb0d173e 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -6,7 +6,6 @@
 #include <rang.hpp>
 
 #include <Connectivity1D.hpp>
-#include <Connectivity2D.hpp>
 #include <Connectivity3D.hpp>
 
 #include <Mesh.hpp>
@@ -884,7 +883,7 @@ GmshReader::__proceedData()
 
     std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
     for (unsigned int e=0; e<__edges.extent(0); ++e) {
-      const unsigned int edge_number = connectivity.getFaceNumber(__edges[e][0], __edges[e][1]);
+      const unsigned int edge_number = connectivity.getFaceNumber({__edges[e][0], __edges[e][1]});
       const unsigned int& ref = __edges_ref[e];
       ref_faces_map[ref].push_back(edge_number);
     }
diff --git a/src/mesh/GmshReader.hpp b/src/mesh/GmshReader.hpp
index f02cfa16f116d2df744cc73bacc073b8964526de..8825a67df3d6cf1593ff54bb6af8b7a96b079b5b 100644
--- a/src/mesh/GmshReader.hpp
+++ b/src/mesh/GmshReader.hpp
@@ -4,7 +4,6 @@
 #include <Kokkos_Core.hpp>
 #include <TinyVector.hpp>
 
-#include <Connectivity2D.hpp>
 #include <Mesh.hpp>
 
 #include <RefId.hpp>
diff --git a/src/scheme/SubItemValuePerItem.hpp b/src/scheme/SubItemValuePerItem.hpp
index 225001939c8e3a40507d34e6d12efa5a35a2d13c..80acb50a8e17119d21f6618414984c0f826bd1bf 100644
--- a/src/scheme/SubItemValuePerItem.hpp
+++ b/src/scheme/SubItemValuePerItem.hpp
@@ -163,10 +163,13 @@ class SubItemValuePerItem
 
   template <typename ConnectivityType>
   SubItemValuePerItem(const ConnectivityType& connectivity)
-      : m_subitem_id_per_item_matrix(connectivity.subItemIdPerItemMatrix()),
-        m_values("values", m_subitem_id_per_item_matrix.entries.extent(0))
   {
-    ;
+    if constexpr (ConnectivityType::dimension > 1) {
+      m_subitem_id_per_item_matrix = connectivity.template subItemIdPerItemMatrix<SubItemType,ItemType>();
+    } else {
+      m_subitem_id_per_item_matrix = connectivity.subItemIdPerItemMatrix();
+    }
+    m_values = Kokkos::View<DataType*>("values", m_subitem_id_per_item_matrix.entries.extent(0));
   }
 
   ~SubItemValuePerItem() = default;