diff --git a/src/main.cpp b/src/main.cpp
index 8fab5c67e1ee868008b27d4d7ee4efbb6cd14d85..fb2da6d32b435e5257739d05668950707b3bbd20 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -226,7 +226,7 @@ int main(int argc, char *argv[])
     // class for acoustic solver test
     Kokkos::Timer timer;
     timer.reset();
-    Connectivity1D connectivity(number);
+    std::shared_ptr<Connectivity1D>connectivity( new Connectivity1D(number));
     typedef Mesh<Connectivity1D> MeshType;
     typedef MeshData<MeshType> MeshDataType;
     typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType;
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 0bb82b29eb2a5849dcedbff4f2edaf55a7050a0d..45151c589c82774f3e289c855403a9c6d8bdaf04 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -805,8 +805,8 @@ GmshReader::__proceedData()
     for (size_t j=nb_triangles; j<nb_triangles+nb_quadrangles; ++j) {
       cell_nb_nodes[j] = 4;
     }
-    m_connectivity = new Connectivity2D(cell_nb_nodes, cell_nodes);
-    Connectivity2D& connectivity = *m_connectivity;
+    Connectivity2D* p_connectivity = new Connectivity2D(cell_nb_nodes, cell_nodes);
+    Connectivity2D& connectivity = *p_connectivity;
 
     std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
     for (unsigned int e=0; e<__edges.extent(0); ++e) {
@@ -856,7 +856,8 @@ GmshReader::__proceedData()
       xr[i][0] = __vertices[i][0];
       xr[i][1] = __vertices[i][1];
     }
-    m_mesh = new MeshType(connectivity, xr);
+    m_mesh = new MeshType(std::shared_ptr<Connectivity2D>(p_connectivity),
+			  xr);
     MeshType& mesh = *m_mesh;
 
     std::ofstream gnuplot("mesh.gnu");
diff --git a/src/mesh/GmshReader.hpp b/src/mesh/GmshReader.hpp
index 5cb7e9e7b6c76d9d1a9615eb122bab673def1050..1baddb1ee729294c5a6500ab086069c8b2893b86 100644
--- a/src/mesh/GmshReader.hpp
+++ b/src/mesh/GmshReader.hpp
@@ -184,7 +184,6 @@ private:
    */
   void __readPhysicalNames2_2();
 
-  Connectivity2D* m_connectivity;
   Mesh<Connectivity2D>* m_mesh;
 public:
 
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index f7fa0cde68a6bfa0adbc52121fe22b7c99fa52ec..40936ce84b9eb656d0d74b10f2a792c5b2de2966 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -4,6 +4,8 @@
 #include <Kokkos_Core.hpp>
 #include <TinyVector.hpp>
 
+#include <memory>
+
 template <typename ConnectivityType>
 class Mesh
 {
@@ -13,29 +15,29 @@ public:
   typedef TinyVector<dimension> Rd;
 
 private:
-  const Connectivity& m_connectivity;
+  const std::shared_ptr<Connectivity> m_connectivity;
   Kokkos::View<Rd*> m_xr;
 
 public:
 
   const Connectivity& connectivity() const
   {
-    return m_connectivity;
+    return *m_connectivity;
   }
 
   const size_t& numberOfNodes() const
   {
-    return m_connectivity.numberOfNodes();
+    return m_connectivity->numberOfNodes();
   }
 
   const size_t& numberOfFaces() const
   {
-    return m_connectivity.numberOfFaces();
+    return m_connectivity->numberOfFaces();
   }
 
   const size_t& numberOfCells() const
   {
-    return m_connectivity.numberOfCells();
+    return m_connectivity->numberOfCells();
   }
 
   Kokkos::View<Rd*> xr() const
@@ -44,19 +46,19 @@ public:
   }
 
   KOKKOS_INLINE_FUNCTION
-  Mesh(const Connectivity& connectivity)
+  Mesh(const std::shared_ptr<Connectivity>& connectivity)
     : m_connectivity(connectivity),
-      m_xr("xr", connectivity.numberOfNodes())
+      m_xr("xr", connectivity->numberOfNodes())
   {
     static_assert (Connectivity::dimension ==1,"no automatic calculation of vertices in dim>1");
-    const double delta_x = 1./connectivity.numberOfCells();
-    Kokkos::parallel_for(connectivity.numberOfNodes(), KOKKOS_LAMBDA(const int& r){
+    const double delta_x = 1./connectivity->numberOfCells();
+    Kokkos::parallel_for(connectivity->numberOfNodes(), KOKKOS_LAMBDA(const int& r){
 	m_xr[r][0] = r*delta_x;
       });
   }
 
   KOKKOS_INLINE_FUNCTION
-  Mesh(const Connectivity& connectivity,
+  Mesh(const std::shared_ptr<Connectivity>& connectivity,
        Kokkos::View<Rd*>& xr)
     : m_connectivity(connectivity),
       m_xr(xr)