diff --git a/src/mesh/Connectivity1D.hpp b/src/mesh/Connectivity1D.hpp
index b78a8bf330eca6213ede9e20f1efd18d52bf50c2..21d89450c7ea14d0d9ac8f8467c5e929d98a057f 100644
--- a/src/mesh/Connectivity1D.hpp
+++ b/src/mesh/Connectivity1D.hpp
@@ -14,6 +14,9 @@ class Connectivity1D
 {
 public:
   static constexpr size_t dimension = 1;
+
+  ConnectivityMatrix m_cell_to_node_matrix;
+
   ConnectivityMatrix m_node_to_cell_matrix;
   ConnectivityMatrixShort m_node_to_cell_local_node_matrix;
 
@@ -22,12 +25,12 @@ private:
 
   const size_t m_number_of_cells;
 
-  const Kokkos::View<const unsigned int**> m_cell_nodes;
-  const Kokkos::View<const unsigned int**>& m_cell_faces = m_cell_nodes;
+  Kokkos::View<const unsigned int**> m_cell_nodes;
+  Kokkos::View<const unsigned int**>& m_cell_faces = m_cell_nodes;
 
-  const Kokkos::View<const unsigned short*> m_cell_nb_nodes;
+  Kokkos::View<const unsigned short*> m_cell_nb_nodes;
   Kokkos::View<double*> m_inv_cell_nb_nodes;
-  const Kokkos::View<const unsigned short*>& m_cell_nb_faces = m_cell_nb_nodes;
+  Kokkos::View<const unsigned short*>& m_cell_nb_faces = m_cell_nb_nodes;
 
   size_t m_max_nb_node_per_cell;
 
@@ -139,13 +142,30 @@ public:
                                       m_node_to_cell_local_node_matrix);
   }
 
-  Connectivity1D(const Kokkos::View<const unsigned short*> cell_nb_nodes,
-                 const Kokkos::View<const unsigned int**> cell_nodes)
-    : m_number_of_cells (cell_nb_nodes.extent(0)),
-      m_cell_nodes (cell_nodes),
-      m_cell_nb_nodes (cell_nb_nodes),
-      m_inv_cell_nb_nodes ("inv_cell_nb_nodes", m_number_of_cells)
-  {
+  Connectivity1D(const std::vector<std::vector<unsigned int>>& cell_by_node_vector)
+      : m_number_of_cells (cell_by_node_vector.size()),
+        m_inv_cell_nb_nodes ("inv_cell_nb_nodes", m_number_of_cells)
+  {
+    m_cell_to_node_matrix
+        = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("cell_to_node_matrix",
+                                                            cell_by_node_vector);
+    {
+      Kokkos::View<unsigned short*> cell_nb_nodes("cell_nb_nodes", cell_by_node_vector.size());
+      for (size_t j=0; j<cell_by_node_vector.size(); ++j) {
+        cell_nb_nodes[j] = cell_by_node_vector[j].size();
+      }
+      m_cell_nb_nodes = cell_nb_nodes;
+    }
+    {
+      Kokkos::View<unsigned int**> cell_nodes("cell_nodes", cell_by_node_vector.size(), 2);
+      for (size_t j=0; j<cell_by_node_vector.size(); ++j) {
+        for (size_t R=0; R<cell_by_node_vector[j].size(); ++R) {
+          cell_nodes(j,R) = cell_by_node_vector[j][R];
+        }
+      }
+      m_cell_nodes = cell_nodes;
+    }
+
     Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const size_t& j) {
         m_inv_cell_nb_nodes[j] = 1./m_cell_nb_nodes[j];
       });
diff --git a/src/mesh/Connectivity2D.hpp b/src/mesh/Connectivity2D.hpp
index 4eecf748545721b785106a18c516a66032b39a01..0bab4bf58d349575eb42f2baa920bce85f71c382 100644
--- a/src/mesh/Connectivity2D.hpp
+++ b/src/mesh/Connectivity2D.hpp
@@ -19,11 +19,14 @@ class Connectivity2D
  public:
   static constexpr size_t dimension = 2;
 
-  ConnectivityMatrix m_node_to_cell_matrix;
-  ConnectivityMatrixShort m_node_to_cell_local_node_matrix;
+  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;
+
  private:
   std::vector<RefFaceList> m_ref_face_list;
   std::vector<RefNodeList> m_ref_node_list;
@@ -31,8 +34,8 @@ class Connectivity2D
   const size_t  m_number_of_cells;
   size_t m_number_of_faces;
 
-  const Kokkos::View<const unsigned short*> m_cell_nb_nodes;
-  const Kokkos::View<const unsigned int**> m_cell_nodes;
+  Kokkos::View<const unsigned short*> m_cell_nb_nodes;
+  Kokkos::View<const unsigned int**> m_cell_nodes;
   Kokkos::View<double*> m_inv_cell_nb_nodes;
 
   Kokkos::View<const unsigned short*> m_cell_nb_faces;
@@ -260,12 +263,30 @@ class Connectivity2D
 
   Connectivity2D(const Connectivity2D&) = delete;
 
-  Connectivity2D(const Kokkos::View<const unsigned short*> cell_nb_nodes,
-                 const Kokkos::View<const unsigned int**> cell_nodes)
-      : m_number_of_cells (cell_nodes.extent(0)),
-        m_cell_nb_nodes(cell_nb_nodes),
-        m_cell_nodes (cell_nodes)
+  Connectivity2D(std::vector<std::vector<unsigned int>> cell_by_node_vector)
+      : m_number_of_cells (cell_by_node_vector.size())
   {
+    m_cell_to_node_matrix
+        = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("cell_to_node_matrix",
+                                                            cell_by_node_vector);
+
+    {
+      Kokkos::View<unsigned short*> cell_nb_nodes("cell_nb_nodes", cell_by_node_vector.size());
+      for (size_t j=0; j<cell_by_node_vector.size(); ++j) {
+        cell_nb_nodes[j] = cell_by_node_vector[j].size();
+      }
+      m_cell_nb_nodes = cell_nb_nodes;
+    }
+    {
+      Kokkos::View<unsigned int**> cell_nodes("cell_nodes", cell_by_node_vector.size(), 8);
+      for (size_t j=0; j<cell_by_node_vector.size(); ++j) {
+        for (size_t R=0; R<cell_by_node_vector[j].size(); ++R) {
+          cell_nodes(j,R) = cell_by_node_vector[j][R];
+        }
+      }
+      m_cell_nodes = cell_nodes;
+    }
+
     {
       Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", m_number_of_cells);
       Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const int&j){
diff --git a/src/mesh/Connectivity3D.hpp b/src/mesh/Connectivity3D.hpp
index c509e8e5bcdf029bfd00bd9932a9c91d794a462f..3a92bcfc9f3b88b17cdc6d550cf9677020380c1f 100644
--- a/src/mesh/Connectivity3D.hpp
+++ b/src/mesh/Connectivity3D.hpp
@@ -22,6 +22,8 @@ class Connectivity3D
 public:
   static constexpr size_t dimension = 3;
 
+  ConnectivityMatrix m_cell_to_node_matrix;
+
   ConnectivityMatrix m_cell_to_face_matrix;
   ConnectivityMatrixShort m_cell_to_face_is_reversed_matrix;
 
@@ -39,8 +41,8 @@ private:
   const size_t m_number_of_cells;
   size_t m_number_of_faces;
 
-  const Kokkos::View<const unsigned short*> m_cell_nb_nodes;
-  const Kokkos::View<const unsigned int**> m_cell_nodes;
+  Kokkos::View<const unsigned short*> m_cell_nb_nodes;
+  Kokkos::View<const unsigned int**> m_cell_nodes;
   Kokkos::View<double*> m_inv_cell_nb_nodes;
 
   Kokkos::View<const unsigned short*> m_node_nb_faces;
@@ -477,12 +479,28 @@ private:
 
   Connectivity3D(const Connectivity3D&) = delete;
 
-  Connectivity3D(const Kokkos::View<const unsigned short*> cell_nb_nodes,
-                 const Kokkos::View<const unsigned int**> cell_nodes)
-    : m_number_of_cells (cell_nodes.extent(0)),
-      m_cell_nb_nodes(cell_nb_nodes),
-      m_cell_nodes (cell_nodes)
+  Connectivity3D(const std::vector<std::vector<unsigned int>>& cell_by_node_vector)
+      : m_number_of_cells (cell_by_node_vector.size())
   {
+    m_cell_to_node_matrix
+        = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("cell_to_node_matrix",
+                                                            cell_by_node_vector);
+    {
+      Kokkos::View<unsigned short*> cell_nb_nodes("cell_nb_nodes", cell_by_node_vector.size());
+      for (size_t j=0; j<cell_by_node_vector.size(); ++j) {
+        cell_nb_nodes[j] = cell_by_node_vector[j].size();
+      }
+      m_cell_nb_nodes = cell_nb_nodes;
+    }
+    {
+      Kokkos::View<unsigned int**> cell_nodes("cell_nodes", cell_by_node_vector.size(), 8);
+      for (size_t j=0; j<cell_by_node_vector.size(); ++j) {
+        for (size_t R=0; R<cell_by_node_vector[j].size(); ++R) {
+          cell_nodes(j,R) = cell_by_node_vector[j][R];
+        }
+      }
+      m_cell_nodes = cell_nodes;
+    }
     {
       Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", m_number_of_cells);
       Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const int& j){
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index c8c2e09e7a28bf484cac0840b41d964f87fc994a..11516ed001bf8bf182ca55c4f00218f183770a31 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -799,29 +799,25 @@ GmshReader::__proceedData()
     if (elementNumber[4] > 0) {
       max_nb_node_per_cell = 8;
     }
-    const Kokkos::View<unsigned int**> cell_nodes("cell_nodes", nb_cells, max_nb_node_per_cell);
+
+    std::vector<std::vector<unsigned int>> cell_by_node_vector(nb_cells);
     const size_t nb_tetrahedra = __tetrahedra.extent(0);
     for (size_t j=0; j<nb_tetrahedra; ++j) {
+      cell_by_node_vector[j].resize(4);
       for (int r=0; r<4; ++r) {
-        cell_nodes(j,r) = __tetrahedra[j][r];
+        cell_by_node_vector[j][r] = __tetrahedra[j][r];
       }
     }
     const size_t nb_hexahedra = __hexahedra.extent(0);
     for (size_t j=0; j<nb_hexahedra; ++j) {
-      const size_t jh = j+nb_tetrahedra;
+      const size_t jh = nb_tetrahedra+j;
+      cell_by_node_vector[jh].resize(8);
       for (int r=0; r<8; ++r) {
-        cell_nodes(jh,r) = __hexahedra[j][r];
+        cell_by_node_vector[jh][r] = __hexahedra[j][r];
       }
     }
-    const Kokkos::View<unsigned short*> cell_nb_nodes("cell_nb_nodes", nb_cells);
-    for (size_t j=0; j<nb_tetrahedra; ++j) {
-      cell_nb_nodes[j] = 4;
-    }
-    for (size_t j=nb_tetrahedra; j<nb_tetrahedra+nb_hexahedra; ++j) {
-      cell_nb_nodes[j] = 8;
-    }
 
-    std::shared_ptr<Connectivity3D> p_connectivity(new Connectivity3D(cell_nb_nodes, cell_nodes));
+    std::shared_ptr<Connectivity3D> p_connectivity(new Connectivity3D(cell_by_node_vector));
     Connectivity3D& connectivity = *p_connectivity;
 
     std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
@@ -847,16 +843,6 @@ GmshReader::__proceedData()
       connectivity.addRefFaceList(RefFaceList(physical_ref_id.refId(), face_list));
     }
 
-    // for (size_t j=0; j<nb_cells; ++j) {
-    //   std::cout << std::setw(5) << j << ": ";
-    //   for (size_t r=0; r<cell_nb_nodes[j]; ++r) {
-    //     std::cout << ' ' << cell_nodes(j,r);
-    //   }
-    //   std::cout << '\n';
-    // }
-
-    // std::cerr << "*** using a 3d mesh (NIY)\n";
-    // std::exit(0);
     typedef Mesh<Connectivity3D> MeshType;
     typedef TinyVector<3, double> Rd;
 
@@ -874,30 +860,49 @@ GmshReader::__proceedData()
     if (elementNumber[2] > 0) {
       max_nb_node_per_cell = 4;
     }
-    const Kokkos::View<unsigned int**> cell_nodes("cell_nodes", nb_cells, max_nb_node_per_cell);
+
+    std::vector<std::vector<unsigned int>> cell_by_node_vector(nb_cells);
     const size_t nb_triangles = __triangles.extent(0);
     for (size_t j=0; j<nb_triangles; ++j) {
-      cell_nodes(j,0) = __triangles[j][0];
-      cell_nodes(j,1) = __triangles[j][1];
-      cell_nodes(j,2) = __triangles[j][2];
+      cell_by_node_vector[j].resize(3);
+      for (int r=0; r<3; ++r) {
+        cell_by_node_vector[j][r] = __triangles[j][r];
+      }
     }
 
     const size_t nb_quadrangles = __quadrangles.extent(0);
     for (size_t j=0; j<nb_quadrangles; ++j) {
       const size_t jq = j+nb_triangles;
-      cell_nodes(jq,0) = __quadrangles[j][0];
-      cell_nodes(jq,1) = __quadrangles[j][1];
-      cell_nodes(jq,2) = __quadrangles[j][2];
-      cell_nodes(jq,3) = __quadrangles[j][3];
-    }
-    const Kokkos::View<unsigned short*> cell_nb_nodes("cell_nb_nodes", nb_cells);
-    for (size_t j=0; j<nb_triangles; ++j) {
-      cell_nb_nodes[j] = 3;
-    }
-    for (size_t j=nb_triangles; j<nb_triangles+nb_quadrangles; ++j) {
-      cell_nb_nodes[j] = 4;
+      cell_by_node_vector[jq].resize(4);
+      for (int r=0; r<4; ++r) {
+        cell_by_node_vector[jq][r] = __quadrangles[j][r];
+      }
     }
-    std::shared_ptr<Connectivity2D> p_connectivity(new Connectivity2D(cell_nb_nodes, cell_nodes));
+
+    // const Kokkos::View<unsigned int**> cell_nodes("cell_nodes", nb_cells, max_nb_node_per_cell);
+
+    // for (size_t j=0; j<nb_triangles; ++j) {
+    //   cell_nodes(j,0) = __triangles[j][0];
+    //   cell_nodes(j,1) = __triangles[j][1];
+    //   cell_nodes(j,2) = __triangles[j][2];
+    // }
+
+    // const size_t nb_quadrangles = __quadrangles.extent(0);
+    // for (size_t j=0; j<nb_quadrangles; ++j) {
+    //   const size_t jq = j+nb_triangles;
+    //   cell_nodes(jq,0) = __quadrangles[j][0];
+    //   cell_nodes(jq,1) = __quadrangles[j][1];
+    //   cell_nodes(jq,2) = __quadrangles[j][2];
+    //   cell_nodes(jq,3) = __quadrangles[j][3];
+    // }
+    // const Kokkos::View<unsigned short*> cell_nb_nodes("cell_nb_nodes", nb_cells);
+    // for (size_t j=0; j<nb_triangles; ++j) {
+    //   cell_nb_nodes[j] = 3;
+    // }
+    // for (size_t j=nb_triangles; j<nb_triangles+nb_quadrangles; ++j) {
+    //   cell_nb_nodes[j] = 4;
+    // }
+    std::shared_ptr<Connectivity2D> p_connectivity(new Connectivity2D(cell_by_node_vector));
     Connectivity2D& connectivity = *p_connectivity;
 
     std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
@@ -945,20 +950,16 @@ GmshReader::__proceedData()
   } else if ((dimension1_mask, elementNumber)>0) {
 
     const size_t nb_cells = (dimension1_mask, elementNumber);
-    const size_t max_nb_node_per_cell=2;
-
-    const Kokkos::View<unsigned int**> cell_nodes("cell_nodes", nb_cells, max_nb_node_per_cell);
-    for (size_t j=0; j<nb_cells; ++j) {
-      cell_nodes(j,0) = __edges[j][0];
-      cell_nodes(j,1) = __edges[j][1];
-    }
 
-    const Kokkos::View<unsigned short*> cell_nb_nodes("cell_nb_nodes", nb_cells);
+    std::vector<std::vector<unsigned int>> cell_by_node_vector(nb_cells);
     for (size_t j=0; j<nb_cells; ++j) {
-      cell_nb_nodes[j] = 2;
+      cell_by_node_vector[j].resize(2);
+      for (int r=0; r<2; ++r) {
+        cell_by_node_vector[j][r] = __edges[j][r];
+      }
     }
 
-    std::shared_ptr<Connectivity1D> p_connectivity(new Connectivity1D(cell_nb_nodes, cell_nodes));
+    std::shared_ptr<Connectivity1D> p_connectivity(new Connectivity1D(cell_by_node_vector));
     Connectivity1D& connectivity = *p_connectivity;
 
     std::map<unsigned int, std::vector<unsigned int>> ref_points_map;