diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 06e241bf3f59e4f04f8a3e553fbc81a83cf3cf76..9a367ef65a4770ac51b4b701b5fe39fe70848f59 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -63,19 +63,19 @@ struct GlaceScheme
 {
   using ConnectivityType = Connectivity<Dimension>;
   using MeshType         = Mesh<ConnectivityType>;
-  using MeshDataType     = MeshDataLegacy<MeshType>;
+  using MeshDataType     = MeshData<Dimension>;
   using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
 
-  const MeshType& m_mesh;
+  std::shared_ptr<const MeshType> m_mesh;
 
-  GlaceScheme(const IMesh& mesh,
+  GlaceScheme(std::shared_ptr<const IMesh> i_mesh,
               const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
               const FunctionSymbolId& rho_id,
               const FunctionSymbolId& u_id,
               const FunctionSymbolId& p_id)
-    : m_mesh{dynamic_cast<const MeshType&>(mesh)}
+    : m_mesh{std::dynamic_pointer_cast<const MeshType>(i_mesh)}
   {
-    MeshDataType mesh_data(m_mesh);
+    MeshDataType mesh_data(*m_mesh);
 
     std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
 
@@ -95,8 +95,8 @@ struct GlaceScheme
           const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
             dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
           for (size_t i_ref_face_list = 0;
-               i_ref_face_list < m_mesh.connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
-            const auto& ref_face_list = m_mesh.connectivity().template refItemList<FaceType>(i_ref_face_list);
+               i_ref_face_list < m_mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+            const auto& ref_face_list = m_mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
             const RefId& ref          = ref_face_list.refId();
             if (ref == sym_bc_descriptor.boundaryDescriptor()) {
               SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
@@ -129,7 +129,7 @@ struct GlaceScheme
                                                                                                                  .xj());
     unknowns.gammaj().fill(1.4);
 
-    AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
+    AcousticSolver acoustic_solver(m_mesh, bc_list);
 
     const CellValue<const double>& Vj = mesh_data.Vj();
 
@@ -153,22 +153,22 @@ struct GlaceScheme
     block_eos.updateEandCFromRhoP();
 
     parallel_for(
-      m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { Ej[j] = ej[j] + 0.5 * (uj[j], uj[j]); });
+      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { Ej[j] = ej[j] + 0.5 * (uj[j], uj[j]); });
 
     parallel_for(
-      m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { mj[j] = rhoj[j] * Vj[j]; });
+      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { mj[j] = rhoj[j] * Vj[j]; });
 
     parallel_for(
-      m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { inv_mj[j] = 1. / mj[j]; });
+      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { inv_mj[j] = 1. / mj[j]; });
 
     VTKWriter vtk_writer("mesh_" + std::to_string(Dimension), 0.01);
 
     while ((t < tmax) and (iteration < itermax)) {
       vtk_writer.write(m_mesh,
                        {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                        NamedItemValue{"coords", m_mesh.xr()},
-                        NamedItemValue{"cell_owner", m_mesh.connectivity().cellOwner()},
-                        NamedItemValue{"node_owner", m_mesh.connectivity().nodeOwner()}},
+                        NamedItemValue{"coords", m_mesh->xr()},
+                        NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()},
+                        NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}},
                        t);
       double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
       if (t + dt > tmax) {
@@ -192,9 +192,9 @@ struct GlaceScheme
               << rang::style::reset << '\n';
     vtk_writer.write(m_mesh,
                      {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                      NamedItemValue{"coords", m_mesh.xr()},
-                      NamedItemValue{"cell_owner", m_mesh.connectivity().cellOwner()},
-                      NamedItemValue{"node_owner", m_mesh.connectivity().nodeOwner()}},
+                      NamedItemValue{"coords", m_mesh->xr()},
+                      NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()},
+                      NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}},
                      t, true);   // forces last output
   }
 };
@@ -249,15 +249,15 @@ SchemeModule::SchemeModule()
                                  const FunctionSymbolId& p_id) -> void {
                                 switch (p_mesh->dimension()) {
                                 case 1: {
-                                  GlaceScheme<1>{*p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
+                                  GlaceScheme<1>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
                                   break;
                                 }
                                 case 2: {
-                                  GlaceScheme<2>{*p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
+                                  GlaceScheme<2>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
                                   break;
                                 }
                                 case 3: {
-                                  GlaceScheme<3>{*p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
+                                  GlaceScheme<3>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
                                   break;
                                 }
                                 default: {
diff --git a/src/language/modules/VTKModule.cpp b/src/language/modules/VTKModule.cpp
index 30d2c11e48b580622ae24e30c9d71bb2d46b007d..f815a78e763b387ee5707d9ac65ec1f038a36b14 100644
--- a/src/language/modules/VTKModule.cpp
+++ b/src/language/modules/VTKModule.cpp
@@ -20,22 +20,25 @@ VTKModule::VTKModule()
 
                                 switch (p_mesh->dimension()) {
                                 case 1: {
-                                  using MeshType       = Mesh<Connectivity<1>>;
-                                  const MeshType& mesh = dynamic_cast<const MeshType&>(*p_mesh);
+                                  using MeshType = Mesh<Connectivity<1>>;
+                                  const std::shared_ptr<const MeshType> mesh =
+                                    std::dynamic_pointer_cast<const MeshType>(p_mesh);
 
                                   writer.write(mesh, OutputNamedItemValueSet{}, time, true);
                                   break;
                                 }
                                 case 2: {
-                                  using MeshType       = Mesh<Connectivity<2>>;
-                                  const MeshType& mesh = dynamic_cast<const MeshType&>(*p_mesh);
+                                  using MeshType = Mesh<Connectivity<2>>;
+                                  const std::shared_ptr<const MeshType> mesh =
+                                    std::dynamic_pointer_cast<const MeshType>(p_mesh);
 
                                   writer.write(mesh, OutputNamedItemValueSet{}, time, true);
                                   break;
                                 }
                                 case 3: {
-                                  using MeshType       = Mesh<Connectivity<3>>;
-                                  const MeshType& mesh = dynamic_cast<const MeshType&>(*p_mesh);
+                                  using MeshType = Mesh<Connectivity<3>>;
+                                  const std::shared_ptr<const MeshType> mesh =
+                                    std::dynamic_pointer_cast<const MeshType>(p_mesh);
 
                                   writer.write(mesh, OutputNamedItemValueSet{}, time, true);
                                   break;
diff --git a/src/main.cpp b/src/main.cpp
index 2fabfe269a9a88ad13ef8965cdf144d33ad044a6..aa9e2e3ff6378c1f1e20fe27388ea6019e765b08 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -44,7 +44,9 @@ main(int argc, char* argv[])
   std::regex gmsh_regex("(.*).msh");
   if (not std::regex_match(filename, gmsh_regex)) {
     SynchronizerManager::create();
+    MeshDataManager::create();
     parser(filename);
+    MeshDataManager::destroy();
     SynchronizerManager::destroy();
     finalize();
     return 0;
@@ -52,6 +54,7 @@ main(int argc, char* argv[])
 
   std::map<std::string, double> method_cost_map;
 
+  MeshDataManager::create();
   SynchronizerManager::create();
 
   if (filename != "") {
@@ -78,14 +81,14 @@ main(int argc, char* argv[])
 
       using ConnectivityType = Connectivity1D;
       using MeshType         = Mesh<ConnectivityType>;
-      using MeshDataType     = MeshDataLegacy<MeshType>;
+      using MeshDataType     = MeshData<1>;
       using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
 
-      const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
+      std::shared_ptr<const MeshType> mesh = std::dynamic_pointer_cast<const MeshType>(gmsh_reader.mesh());
 
       Timer timer;
       timer.reset();
-      MeshDataType mesh_data(mesh);
+      MeshDataType mesh_data(*mesh);
 
       std::vector<BoundaryConditionHandler> bc_list;
       {
@@ -95,8 +98,8 @@ main(int argc, char* argv[])
             const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
               dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
             for (size_t i_ref_node_list = 0;
-                 i_ref_node_list < mesh.connectivity().numberOfRefItemList<ItemType::node>(); ++i_ref_node_list) {
-              const RefNodeList& ref_node_list = mesh.connectivity().refItemList<ItemType::node>(i_ref_node_list);
+                 i_ref_node_list < mesh->connectivity().numberOfRefItemList<ItemType::node>(); ++i_ref_node_list) {
+              const RefNodeList& ref_node_list = mesh->connectivity().refItemList<ItemType::node>(i_ref_node_list);
               const RefId& ref                 = ref_node_list.refId();
               if (ref == sym_bc_descriptor.boundaryDescriptor()) {
                 SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
@@ -119,7 +122,7 @@ main(int argc, char* argv[])
 
       unknowns.initializeSod();
 
-      AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
+      AcousticSolver acoustic_solver(mesh, bc_list);
 
       using Rd = TinyVector<MeshType::Dimension>;
 
@@ -144,9 +147,9 @@ main(int argc, char* argv[])
       while ((t < tmax) and (iteration < itermax)) {
         vtk_writer.write(mesh,
                          {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                          NamedItemValue{"coords", mesh.xr()},
-                          NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
-                          NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
+                          NamedItemValue{"coords", mesh->xr()},
+                          NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
+                          NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
                          t);
         double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
         if (t + dt > tmax) {
@@ -161,9 +164,9 @@ main(int argc, char* argv[])
       }
       vtk_writer.write(mesh,
                        {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                        NamedItemValue{"coords", mesh.xr()},
-                        NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
-                        NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
+                        NamedItemValue{"coords", mesh->xr()},
+                        NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
+                        NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
                        t, true);   // forces last output
 
       std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
@@ -175,7 +178,7 @@ main(int argc, char* argv[])
         const CellValue<const Rd>& xj       = mesh_data.xj();
         const CellValue<const double>& rhoj = unknowns.rhoj();
         std::ofstream fout("rho");
-        for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+        for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
           fout << xj[j][0] << ' ' << rhoj[j] << '\n';
         }
       }
@@ -197,14 +200,14 @@ main(int argc, char* argv[])
 
       using ConnectivityType = Connectivity2D;
       using MeshType         = Mesh<ConnectivityType>;
-      using MeshDataType     = MeshDataLegacy<MeshType>;
+      using MeshDataType     = MeshData<2>;
       using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
 
-      const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
+      std::shared_ptr<const MeshType> mesh = std::dynamic_pointer_cast<const MeshType>(gmsh_reader.mesh());
 
       Timer timer;
       timer.reset();
-      MeshDataType mesh_data(mesh);
+      MeshDataType mesh_data(*mesh);
 
       std::vector<BoundaryConditionHandler> bc_list;
       {
@@ -214,8 +217,8 @@ main(int argc, char* argv[])
             const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
               dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
             for (size_t i_ref_face_list = 0;
-                 i_ref_face_list < mesh.connectivity().numberOfRefItemList<ItemType::face>(); ++i_ref_face_list) {
-              const RefFaceList& ref_face_list = mesh.connectivity().refItemList<ItemType::face>(i_ref_face_list);
+                 i_ref_face_list < mesh->connectivity().numberOfRefItemList<ItemType::face>(); ++i_ref_face_list) {
+              const RefFaceList& ref_face_list = mesh->connectivity().refItemList<ItemType::face>(i_ref_face_list);
               const RefId& ref                 = ref_face_list.refId();
               if (ref == sym_bc_descriptor.boundaryDescriptor()) {
                 SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
@@ -238,7 +241,7 @@ main(int argc, char* argv[])
 
       unknowns.initializeSod();
 
-      AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
+      AcousticSolver acoustic_solver(mesh, bc_list);
 
       const CellValue<const double>& Vj = mesh_data.Vj();
 
@@ -261,9 +264,9 @@ main(int argc, char* argv[])
       while ((t < tmax) and (iteration < itermax)) {
         vtk_writer.write(mesh,
                          {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                          NamedItemValue{"coords", mesh.xr()},
-                          NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
-                          NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
+                          NamedItemValue{"coords", mesh->xr()},
+                          NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
+                          NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
                          t);
         double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
         if (t + dt > tmax) {
@@ -278,9 +281,9 @@ main(int argc, char* argv[])
       }
       vtk_writer.write(mesh,
                        {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                        NamedItemValue{"coords", mesh.xr()},
-                        NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
-                        NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
+                        NamedItemValue{"coords", mesh->xr()},
+                        NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
+                        NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
                        t, true);   // forces last output
 
       std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
@@ -303,14 +306,14 @@ main(int argc, char* argv[])
 
       using ConnectivityType = Connectivity3D;
       using MeshType         = Mesh<ConnectivityType>;
-      using MeshDataType     = MeshDataLegacy<MeshType>;
+      using MeshDataType     = MeshData<3>;
       using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
 
-      const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
+      std::shared_ptr<const MeshType> mesh = std::dynamic_pointer_cast<const MeshType>(gmsh_reader.mesh());
 
       Timer timer;
       timer.reset();
-      MeshDataType mesh_data(mesh);
+      MeshDataType mesh_data(*mesh);
 
       std::vector<BoundaryConditionHandler> bc_list;
       {
@@ -320,8 +323,8 @@ main(int argc, char* argv[])
             const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
               dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
             for (size_t i_ref_face_list = 0;
-                 i_ref_face_list < mesh.connectivity().numberOfRefItemList<ItemType::face>(); ++i_ref_face_list) {
-              const RefFaceList& ref_face_list = mesh.connectivity().refItemList<ItemType::face>(i_ref_face_list);
+                 i_ref_face_list < mesh->connectivity().numberOfRefItemList<ItemType::face>(); ++i_ref_face_list) {
+              const RefFaceList& ref_face_list = mesh->connectivity().refItemList<ItemType::face>(i_ref_face_list);
               const RefId& ref                 = ref_face_list.refId();
               if (ref == sym_bc_descriptor.boundaryDescriptor()) {
                 SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
@@ -344,7 +347,7 @@ main(int argc, char* argv[])
 
       unknowns.initializeSod();
 
-      AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
+      AcousticSolver acoustic_solver(mesh, bc_list);
 
       const CellValue<const double>& Vj = mesh_data.Vj();
 
@@ -367,9 +370,9 @@ main(int argc, char* argv[])
       while ((t < tmax) and (iteration < itermax)) {
         vtk_writer.write(mesh,
                          {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                          NamedItemValue{"coords", mesh.xr()},
-                          NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
-                          NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
+                          NamedItemValue{"coords", mesh->xr()},
+                          NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
+                          NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
                          t);
         double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
         if (t + dt > tmax) {
@@ -383,9 +386,9 @@ main(int argc, char* argv[])
       }
       vtk_writer.write(mesh,
                        {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                        NamedItemValue{"coords", mesh.xr()},
-                        NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
-                        NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
+                        NamedItemValue{"coords", mesh->xr()},
+                        NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
+                        NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
                        t, true);   // forces last output
 
       std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
@@ -403,6 +406,7 @@ main(int argc, char* argv[])
     throw NormalError("Connectivity1D defined by number of nodes no more implemented\n");
   }
 
+  MeshDataManager::destroy();
   SynchronizerManager::destroy();
 
   finalize();
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 2a7853937a3821dc7370b923f6a5d05a23268f8a..92c865c599767db0d5e5fed84e7d16d0cbec8579 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -10,6 +10,7 @@ add_library(
   DiamondDualConnectivityBuilder.cpp
   DiamondDualMeshBuilder.cpp
   GmshReader.cpp
+  IMesh.cpp
   LogicalConnectivityBuilder.cpp
   MeshBuilderBase.cpp
   MeshDataManager.cpp
diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualMeshBuilder.cpp
index c0141d7015d600dcb1a4f06bdf2a7e6cf10b4590..2218541b338e196a913b3c0859b882ad668d252e 100644
--- a/src/mesh/DiamondDualMeshBuilder.cpp
+++ b/src/mesh/DiamondDualMeshBuilder.cpp
@@ -8,6 +8,7 @@
 #include <mesh/ItemValueUtils.hpp>
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
 #include <mesh/RefId.hpp>
 #include <utils/Array.hpp>
 #include <utils/Messenger.hpp>
@@ -26,7 +27,7 @@ DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(
   NodeValue<TinyVector<Dimension>> diamond_xr{diamond_connectivity};
 
   const NodeValue<const TinyVector<Dimension>> primal_xr = primal_mesh.xr();
-  MeshDataLegacy<MeshType> primal_mesh_data{primal_mesh};
+  MeshData<Dimension>& primal_mesh_data                  = MeshDataManager::instance().getMeshData(primal_mesh);
   const CellValue<const TinyVector<Dimension>> primal_xj = primal_mesh_data.xj();
 
   NodeId i_node = 0;
@@ -57,7 +58,7 @@ DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(const Mesh<Connectivity<1>>& p
   NodeValue<TinyVector<1>> diamond_xr{diamond_connectivity};
 
   const NodeValue<const TinyVector<1>> primal_xr = primal_mesh.xr();
-  MeshDataLegacy<MeshType> primal_mesh_data{primal_mesh};
+  MeshData<1>& primal_mesh_data                  = MeshDataManager::instance().getMeshData(primal_mesh);
   const CellValue<const TinyVector<1>> primal_xj = primal_mesh_data.xj();
 
   NodeId next_node_id = 0;
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 5e65fcbb74b655acae2cbd428df35f35fc675973..71284c13ed850c7ecd0da8c3bcab938ef93c5bf9 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -8,7 +8,6 @@
 #include <mesh/ConnectivityDispatcher.hpp>
 #include <mesh/ItemValueUtils.hpp>
 #include <mesh/Mesh.hpp>
-#include <mesh/MeshData.hpp>
 #include <mesh/RefItemList.hpp>
 #include <utils/ArrayUtils.hpp>
 #include <utils/Exceptions.hpp>
diff --git a/src/mesh/IMesh.cpp b/src/mesh/IMesh.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cb2accc141ed09c778237587373494797fee726d
--- /dev/null
+++ b/src/mesh/IMesh.cpp
@@ -0,0 +1,8 @@
+#include <mesh/IMesh.hpp>
+
+#include <mesh/MeshDataManager.hpp>
+
+IMesh::~IMesh()
+{
+  MeshDataManager::instance().deleteMeshData(*this);
+}
diff --git a/src/mesh/IMesh.hpp b/src/mesh/IMesh.hpp
index 3e22eb85bb750a5083f658cad83a32ea7626eabb..cb524849aaf0a7f994d8fb34fe2e97508c5e4f8a 100644
--- a/src/mesh/IMesh.hpp
+++ b/src/mesh/IMesh.hpp
@@ -3,7 +3,7 @@
 
 #include <cstddef>
 
-struct IMesh
+class IMesh
 {
  public:
   virtual size_t dimension() const = 0;
@@ -11,8 +11,8 @@ struct IMesh
   IMesh(const IMesh&) = delete;
   IMesh(IMesh&&)      = delete;
 
-  IMesh()  = default;
-  ~IMesh() = default;
+  IMesh() = default;
+  virtual ~IMesh();
 };
 
 #endif   // I_MESH_HPP
diff --git a/src/mesh/IMeshData.hpp b/src/mesh/IMeshData.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8cad9586d201de6b383bbbfaf94f09e5797078e5
--- /dev/null
+++ b/src/mesh/IMeshData.hpp
@@ -0,0 +1,10 @@
+#ifndef I_MESH_DATA_HPP
+#define I_MESH_DATA_HPP
+
+class IMeshData
+{
+ public:
+  virtual ~IMeshData() = default;
+};
+
+#endif   // I_MESH_DATA_HPP
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 1b209e11dab47036fc1e7ee145b3b0b067e78aff..f8cd50f30bed85e3f905aa6efc877b2612e669e0 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -2,6 +2,7 @@
 #define MESH_DATA_HPP
 
 #include <algebra/TinyVector.hpp>
+#include <mesh/IMeshData.hpp>
 #include <mesh/ItemValue.hpp>
 #include <mesh/SubItemValuePerItem.hpp>
 #include <utils/Exceptions.hpp>
@@ -10,29 +11,35 @@
 
 #include <map>
 
-template <typename M>
-class [[deprecated]] MeshDataLegacy
+template <size_t Dimension>
+class Connectivity;
+
+template <typename ConnectivityType>
+class Mesh;
+
+template <size_t Dimension>
+class MeshData : public IMeshData
 {
  public:
-  using MeshType = M;
-
-  static constexpr size_t Dimension = MeshType::Dimension;
   static_assert(Dimension > 0, "dimension must be strictly positive");
 
+  using MeshType = Mesh<Connectivity<Dimension>>;
+
   using Rd = TinyVector<Dimension>;
 
   static constexpr double inv_Dimension = 1. / Dimension;
 
  private:
   const MeshType& m_mesh;
-  NodeValuePerCell<const Rd> m_Cjr;
-  NodeValuePerCell<const double> m_ljr;
-  NodeValuePerCell<const Rd> m_njr;
-  CellValue<const Rd> m_xj;
-  CellValue<const double> m_Vj;
+  std::shared_ptr<NodeValuePerCell<const Rd>> m_Cjr;
+  std::shared_ptr<NodeValuePerCell<const double>> m_ljr;
+  std::shared_ptr<NodeValuePerCell<const Rd>> m_njr;
+  std::shared_ptr<CellValue<const Rd>> m_xj;
+  std::shared_ptr<CellValue<const double>> m_Vj;
 
   PUGS_INLINE
-  void _updateCenter()
+  void
+  _updateCenter()
   {   // Computes vertices isobarycenter
     if constexpr (Dimension == 1) {
       const NodeValue<const Rd>& xr = m_mesh.xr();
@@ -45,7 +52,7 @@ class [[deprecated]] MeshDataLegacy
           const auto& cell_nodes = cell_to_node_matrix[j];
           xj[j]                  = 0.5 * (xr[cell_nodes[0]] + xr[cell_nodes[1]]);
         });
-      m_xj = xj;
+      m_xj = std::make_shared<CellValue<const Rd>>(xj);
     } else {
       const NodeValue<const Rd>& xr = m_mesh.xr();
 
@@ -62,12 +69,13 @@ class [[deprecated]] MeshDataLegacy
           }
           xj[j] = inv_cell_nb_nodes[j] * X;
         });
-      m_xj = xj;
+      m_xj = std::make_shared<CellValue<const Rd>>(xj);
     }
   }
 
   PUGS_INLINE
-  void _updateVolume()
+  void
+  _updateVolume()
   {
     const NodeValue<const Rd>& xr   = m_mesh.xr();
     const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
@@ -79,15 +87,16 @@ class [[deprecated]] MeshDataLegacy
         const auto& cell_nodes = cell_to_node_matrix[j];
 
         for (size_t R = 0; R < cell_nodes.size(); ++R) {
-          sum_cjr_xr += (xr[cell_nodes[R]], m_Cjr(j, R));
+          sum_cjr_xr += (xr[cell_nodes[R]], (*m_Cjr)(j, R));
         }
         Vj[j] = inv_Dimension * sum_cjr_xr;
       });
-    m_Vj = Vj;
+    m_Vj = std::make_shared<CellValue<const double>>(Vj);
   }
 
   PUGS_INLINE
-  void _updateCjr()
+  void
+  _updateCjr()
   {
     if constexpr (Dimension == 1) {
       // Cjr/njr/ljr are constant overtime
@@ -107,21 +116,21 @@ class [[deprecated]] MeshDataLegacy
               Cjr(j, R)       = Rd{half_xrp_xrm[1], -half_xrp_xrm[0]};
             }
           });
-        m_Cjr = Cjr;
+        m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr);
       }
 
       {
         NodeValuePerCell<double> ljr(m_mesh.connectivity());
         parallel_for(
-          m_Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = l2Norm(m_Cjr[jr]); });
-        m_ljr = ljr;
+          m_Cjr->numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = l2Norm((*m_Cjr)[jr]); });
+        m_ljr = std::make_shared<NodeValuePerCell<const double>>(ljr);
       }
 
       {
         NodeValuePerCell<Rd> njr(m_mesh.connectivity());
         parallel_for(
-          m_Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { njr[jr] = (1. / m_ljr[jr]) * m_Cjr[jr]; });
-        m_njr = njr;
+          m_Cjr->numberOfValues(), PUGS_LAMBDA(size_t jr) { njr[jr] = (1. / (*m_ljr)[jr]) * (*m_Cjr)[jr]; });
+        m_njr = std::make_shared<NodeValuePerCell<const Rd>>(njr);
       }
     } else if (Dimension == 3) {
       const NodeValue<const Rd>& xr = m_mesh.xr();
@@ -195,31 +204,32 @@ class [[deprecated]] MeshDataLegacy
             }
           });
 
-        m_Cjr = Cjr;
+        m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr);
       }
 
       {
         NodeValuePerCell<double> ljr(m_mesh.connectivity());
         parallel_for(
-          m_Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = l2Norm(m_Cjr[jr]); });
-        m_ljr = ljr;
+          m_Cjr->numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = l2Norm((*m_Cjr)[jr]); });
+        m_ljr = std::make_shared<NodeValuePerCell<const double>>(ljr);
       }
 
       {
         NodeValuePerCell<Rd> njr(m_mesh.connectivity());
         parallel_for(
-          m_Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { njr[jr] = (1. / m_ljr[jr]) * m_Cjr[jr]; });
-        m_njr = njr;
+          m_Cjr->numberOfValues(), PUGS_LAMBDA(size_t jr) { njr[jr] = (1. / (*m_ljr)[jr]) * (*m_Cjr)[jr]; });
+        m_njr = std::make_shared<NodeValuePerCell<const Rd>>(njr);
       }
     }
     static_assert((Dimension <= 3), "only 1d, 2d and 3d are implemented");
   }
 
-  void _checkCellVolume() const
+  void
+  _checkCellVolume() const
   {
     bool is_valid = [&] {
       for (CellId j = 0; j < m_mesh.numberOfCells(); ++j) {
-        if (m_Vj[j] <= 0) {
+        if ((*m_Vj)[j] <= 0) {
           return false;
         }
       }
@@ -233,50 +243,60 @@ class [[deprecated]] MeshDataLegacy
 
  public:
   PUGS_INLINE
-  const MeshType& mesh() const
+  const MeshType&
+  mesh() const
   {
     return m_mesh;
   }
 
   PUGS_INLINE
-  const NodeValuePerCell<const Rd>& Cjr() const
+  const NodeValuePerCell<const Rd>
+  Cjr() const
   {
-    return m_Cjr;
+    return *m_Cjr;
   }
 
   PUGS_INLINE
-  const NodeValuePerCell<const double>& ljr() const
+  const NodeValuePerCell<const double>
+  ljr() const
   {
-    return m_ljr;
+    return *m_ljr;
   }
 
   PUGS_INLINE
-  const NodeValuePerCell<const Rd>& njr() const
+  const NodeValuePerCell<const Rd>
+  njr() const
   {
-    return m_njr;
+    return *m_njr;
   }
 
   PUGS_INLINE
-  const CellValue<const Rd>& xj() const
+  const CellValue<const Rd>
+  xj()
   {
-    return m_xj;
+    if (not m_xj) {
+      this->_updateCenter();
+    }
+    return *m_xj;
   }
 
   PUGS_INLINE
-  const CellValue<const double>& Vj() const
+  const CellValue<const double>
+  Vj() const
   {
-    return m_Vj;
+    return *m_Vj;
   }
 
-  void updateAllData()
+  void
+  updateAllData()
   {
     this->_updateCjr();
-    this->_updateCenter();
+    //    this->_updateCenter();
     this->_updateVolume();
     this->_checkCellVolume();
   }
 
-  MeshDataLegacy(const MeshType& mesh) : m_mesh(mesh)
+  MeshData(const MeshType& mesh) : m_mesh(mesh)
   {
     if constexpr (Dimension == 1) {
       // in 1d Cjr are computed once for all
@@ -287,7 +307,7 @@ class [[deprecated]] MeshDataLegacy
             Cjr(j, 0) = -1;
             Cjr(j, 1) = 1;
           });
-        m_Cjr = Cjr;
+        m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr);
       }
       // in 1d njr=Cjr (here performs a shallow copy)
       m_njr = m_Cjr;
@@ -295,16 +315,18 @@ class [[deprecated]] MeshDataLegacy
         NodeValuePerCell<double> ljr(m_mesh.connectivity());
         parallel_for(
           ljr.numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = 1; });
-        m_ljr = ljr;
+        m_ljr = std::make_shared<NodeValuePerCell<const double>>(ljr);
       }
     }
     this->updateAllData();
   }
 
-  MeshDataLegacy(const MeshDataLegacy&) = delete;
-  MeshDataLegacy(MeshDataLegacy &&)     = delete;
+  MeshData() = delete;
+
+  MeshData(const MeshData&) = default;
+  MeshData(MeshData&&)      = default;
 
-  ~MeshDataLegacy()
+  ~MeshData()
   {
     ;
   }
diff --git a/src/mesh/MeshDataManager.cpp b/src/mesh/MeshDataManager.cpp
index a72fda8307d2a20033c206d31998a44931cbf23c..ae77ee8255f62447f8d49cd8f9a85a7aec6ddd10 100644
--- a/src/mesh/MeshDataManager.cpp
+++ b/src/mesh/MeshDataManager.cpp
@@ -1,13 +1,12 @@
 #include <utils/PugsAssert.hpp>
 
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
+#include <utils/Exceptions.hpp>
 
-class MeshData
-{
- public:
-  MeshData() = default;
-};
+#include <sstream>
 
 MeshDataManager* MeshDataManager::m_instance{nullptr};
 
@@ -36,20 +35,27 @@ MeshDataManager::destroy()
 }
 
 void
-MeshDataManager::deleteMeshData(const IMesh* mesh_data)
+MeshDataManager::deleteMeshData(const IMesh& mesh)
 {
-  m_mesh_mesh_data_map.erase(mesh_data);
+  m_mesh_mesh_data_map.erase(&mesh);
 }
 
-MeshData&
-MeshDataManager::getMeshData(const IMesh* mesh)
+template <size_t Dimension>
+MeshData<Dimension>&
+MeshDataManager::getMeshData(const Mesh<Connectivity<Dimension>>& mesh)
 {
-  if (auto connectivity_synchronizer = m_mesh_mesh_data_map.find(mesh);
+  const IMesh* p_mesh = &mesh;
+
+  if (auto connectivity_synchronizer = m_mesh_mesh_data_map.find(p_mesh);
       connectivity_synchronizer != m_mesh_mesh_data_map.end()) {
-    return (*connectivity_synchronizer->second);
+    return dynamic_cast<MeshData<Dimension>&>(*connectivity_synchronizer->second);
   } else {
-    std::shared_ptr mesh_data  = std::make_shared<MeshData>();
-    m_mesh_mesh_data_map[mesh] = mesh_data;
+    std::shared_ptr mesh_data    = std::make_shared<MeshData<Dimension>>(mesh);
+    m_mesh_mesh_data_map[p_mesh] = mesh_data;
     return *mesh_data;
   }
 }
+
+template MeshData<1>& MeshDataManager::getMeshData(const Mesh<Connectivity<1>>&);
+template MeshData<2>& MeshDataManager::getMeshData(const Mesh<Connectivity<2>>&);
+template MeshData<3>& MeshDataManager::getMeshData(const Mesh<Connectivity<3>>&);
diff --git a/src/mesh/MeshDataManager.hpp b/src/mesh/MeshDataManager.hpp
index 66041bec3525d83de7adf8190f2702dd08544dea..7b127148b889bc585082d626b79c435293fa4cd7 100644
--- a/src/mesh/MeshDataManager.hpp
+++ b/src/mesh/MeshDataManager.hpp
@@ -1,6 +1,7 @@
 #ifndef MESH_DATA_MANAGER_HPP
 #define MESH_DATA_MANAGER_HPP
 
+#include <mesh/IMeshData.hpp>
 #include <utils/PugsAssert.hpp>
 #include <utils/PugsMacros.hpp>
 
@@ -8,12 +9,20 @@
 #include <memory>
 
 class IMesh;
+
+template <size_t>
+class Connectivity;
+
+template <typename ConnectivityType>
+class Mesh;
+
+template <size_t Dimension>
 class MeshData;
 
 class MeshDataManager
 {
  private:
-  std::map<const IMesh*, std::shared_ptr<MeshData>> m_mesh_mesh_data_map;
+  std::map<const IMesh*, std::shared_ptr<IMeshData>> m_mesh_mesh_data_map;
 
   static MeshDataManager* m_instance;
 
@@ -35,8 +44,10 @@ class MeshDataManager
     return *m_instance;
   }
 
-  void deleteMeshData(const IMesh*);
-  MeshData& getMeshData(const IMesh*);
+  void deleteMeshData(const IMesh&);
+
+  template <size_t Dimension>
+  MeshData<Dimension>& getMeshData(const Mesh<Connectivity<Dimension>>&);
 };
 
 #endif   // MESH_DATA_MANAGER_HPP
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index e02c9b7520db65081977e4bb1774a348c9cd6e25..a49db8b4b358f52e7dd199a64a55a43b1123de2e 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -35,10 +35,10 @@ class MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
   }
 
   template <typename MeshType>
-  MeshNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
+  MeshNodeBoundary(const std::shared_ptr<const MeshType>& mesh, const RefFaceList& ref_face_list)
   {
     static_assert(Dimension == MeshType::Dimension);
-    const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
+    const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix();
 
     const Array<const FaceId>& face_list = ref_face_list.list();
     parallel_for(
@@ -52,7 +52,7 @@ class MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
     Kokkos::vector<unsigned int> node_ids;
     // not enough but should reduce significantly the number of resizing
     node_ids.reserve(Dimension * face_list.size());
-    const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
+    const auto& face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
 
     for (size_t l = 0; l < face_list.size(); ++l) {
       const FaceId face_number = face_list[l];
@@ -73,7 +73,8 @@ class MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
   }
 
   template <typename MeshType>
-  MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list) : m_node_list(ref_node_list.list())
+  MeshNodeBoundary(const std::shared_ptr<const MeshType>&, const RefNodeList& ref_node_list)
+    : m_node_list(ref_node_list.list())
   {
     static_assert(Dimension == MeshType::Dimension);
   }
@@ -96,16 +97,16 @@ class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension>   // clazy:exclu
   const Rd m_outgoing_normal;
 
   template <typename MeshType>
-  PUGS_INLINE Rd _getNormal(const MeshType& mesh);
+  PUGS_INLINE Rd _getNormal(const std::shared_ptr<const MeshType>& mesh);
 
   template <typename MeshType>
   PUGS_INLINE void _checkBoundaryIsFlat(TinyVector<2, double> normal,
                                         TinyVector<2, double> xmin,
                                         TinyVector<2, double> xmax,
-                                        const MeshType& mesh) const;
+                                        const std::shared_ptr<const MeshType>& mesh) const;
 
   template <typename MeshType>
-  PUGS_INLINE Rd _getOutgoingNormal(const MeshType& mesh);
+  PUGS_INLINE Rd _getOutgoingNormal(const std::shared_ptr<const MeshType>& mesh);
 
  public:
   const Rd&
@@ -118,14 +119,14 @@ class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension>   // clazy:exclu
   MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default;
 
   template <typename MeshType>
-  MeshFlatNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
+  MeshFlatNodeBoundary(const std::shared_ptr<const MeshType>& mesh, const RefFaceList& ref_face_list)
     : MeshNodeBoundary<Dimension>(mesh, ref_face_list), m_outgoing_normal(_getOutgoingNormal(mesh))
   {
     ;
   }
 
   template <typename MeshType>
-  MeshFlatNodeBoundary(const MeshType& mesh, const RefNodeList& ref_node_list)
+  MeshFlatNodeBoundary(const std::shared_ptr<const MeshType>& mesh, const RefNodeList& ref_node_list)
     : MeshNodeBoundary<Dimension>(mesh, ref_node_list), m_outgoing_normal(_getOutgoingNormal(mesh))
   {
     ;
@@ -143,7 +144,7 @@ void
 MeshFlatNodeBoundary<2>::_checkBoundaryIsFlat(TinyVector<2, double> normal,
                                               TinyVector<2, double> xmin,
                                               TinyVector<2, double> xmax,
-                                              const MeshType& mesh) const
+                                              const std::shared_ptr<const MeshType>& mesh) const
 {
   static_assert(MeshType::Dimension == 2);
   using R2 = TinyVector<2, double>;
@@ -151,7 +152,7 @@ MeshFlatNodeBoundary<2>::_checkBoundaryIsFlat(TinyVector<2, double> normal,
   const R2 origin     = 0.5 * (xmin + xmax);
   const double length = l2Norm(xmax - xmin);
 
-  const NodeValue<const R2>& xr = mesh.xr();
+  const NodeValue<const R2>& xr = mesh->xr();
 
   parallel_for(
     m_node_list.size(), PUGS_LAMBDA(size_t r) {
@@ -165,14 +166,14 @@ MeshFlatNodeBoundary<2>::_checkBoundaryIsFlat(TinyVector<2, double> normal,
 template <>
 template <typename MeshType>
 PUGS_INLINE TinyVector<1, double>
-MeshFlatNodeBoundary<1>::_getNormal(const MeshType& mesh)
+MeshFlatNodeBoundary<1>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
 {
   static_assert(MeshType::Dimension == 1);
   using R = TinyVector<1, double>;
 
   const size_t number_of_bc_nodes = [&]() {
     size_t number_of_bc_nodes = 0;
-    auto node_is_owned        = mesh.connectivity().nodeIsOwned();
+    auto node_is_owned        = mesh->connectivity().nodeIsOwned();
     for (size_t i_node = 0; i_node < m_node_list.size(); ++i_node) {
       number_of_bc_nodes += (node_is_owned[m_node_list[i_node]]);
     }
@@ -189,12 +190,12 @@ MeshFlatNodeBoundary<1>::_getNormal(const MeshType& mesh)
 template <>
 template <typename MeshType>
 PUGS_INLINE TinyVector<2, double>
-MeshFlatNodeBoundary<2>::_getNormal(const MeshType& mesh)
+MeshFlatNodeBoundary<2>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
 {
   static_assert(MeshType::Dimension == 2);
   using R2 = TinyVector<2, double>;
 
-  const NodeValue<const R2>& xr = mesh.xr();
+  const NodeValue<const R2>& xr = mesh->xr();
 
   R2 xmin(std::numeric_limits<double>::max(), std::numeric_limits<double>::max());
 
@@ -244,7 +245,7 @@ MeshFlatNodeBoundary<2>::_getNormal(const MeshType& mesh)
 template <>
 template <typename MeshType>
 PUGS_INLINE TinyVector<3, double>
-MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh)
+MeshFlatNodeBoundary<3>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
 {
   static_assert(MeshType::Dimension == 3);
   using R3 = TinyVector<3, double>;
@@ -252,13 +253,12 @@ MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh)
   R3 xmin(std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max());
   R3 ymin = xmin;
   R3 zmin = xmin;
-  ;
 
   R3 xmax = -xmin;
   R3 ymax = xmax;
   R3 zmax = xmax;
 
-  const NodeValue<const R3>& xr = mesh.xr();
+  const NodeValue<const R3>& xr = mesh->xr();
 
   for (size_t r = 0; r < m_node_list.size(); ++r) {
     const R3& x = xr[m_node_list[r]];
@@ -365,7 +365,7 @@ MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh)
 template <>
 template <typename MeshType>
 PUGS_INLINE TinyVector<1, double>
-MeshFlatNodeBoundary<1>::_getOutgoingNormal(const MeshType& mesh)
+MeshFlatNodeBoundary<1>::_getOutgoingNormal(const std::shared_ptr<const MeshType>& mesh)
 {
   static_assert(MeshType::Dimension == 1);
   using R = TinyVector<1, double>;
@@ -375,10 +375,10 @@ MeshFlatNodeBoundary<1>::_getOutgoingNormal(const MeshType& mesh)
   double max_height = 0;
 
   if (m_node_list.size() > 0) {
-    const NodeValue<const R>& xr    = mesh.xr();
-    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+    const NodeValue<const R>& xr    = mesh->xr();
+    const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
 
-    const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+    const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
 
     const NodeId r0      = m_node_list[0];
     const CellId j0      = node_to_cell_matrix[r0][0];
@@ -410,7 +410,7 @@ MeshFlatNodeBoundary<1>::_getOutgoingNormal(const MeshType& mesh)
 template <>
 template <typename MeshType>
 PUGS_INLINE TinyVector<2, double>
-MeshFlatNodeBoundary<2>::_getOutgoingNormal(const MeshType& mesh)
+MeshFlatNodeBoundary<2>::_getOutgoingNormal(const std::shared_ptr<const MeshType>& mesh)
 {
   static_assert(MeshType::Dimension == 2);
   using R2 = TinyVector<2, double>;
@@ -420,10 +420,10 @@ MeshFlatNodeBoundary<2>::_getOutgoingNormal(const MeshType& mesh)
   double max_height = 0;
 
   if (m_node_list.size() > 0) {
-    const NodeValue<const R2>& xr   = mesh.xr();
-    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+    const NodeValue<const R2>& xr   = mesh->xr();
+    const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
 
-    const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+    const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
 
     const NodeId r0      = m_node_list[0];
     const CellId j0      = node_to_cell_matrix[r0][0];
@@ -454,7 +454,7 @@ MeshFlatNodeBoundary<2>::_getOutgoingNormal(const MeshType& mesh)
 template <>
 template <typename MeshType>
 PUGS_INLINE TinyVector<3, double>
-MeshFlatNodeBoundary<3>::_getOutgoingNormal(const MeshType& mesh)
+MeshFlatNodeBoundary<3>::_getOutgoingNormal(const std::shared_ptr<const MeshType>& mesh)
 {
   static_assert(MeshType::Dimension == 3);
   using R3 = TinyVector<3, double>;
@@ -464,10 +464,10 @@ MeshFlatNodeBoundary<3>::_getOutgoingNormal(const MeshType& mesh)
   double max_height = 0;
 
   if (m_node_list.size() > 0) {
-    const NodeValue<const R3>& xr   = mesh.xr();
-    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+    const NodeValue<const R3>& xr   = mesh->xr();
+    const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
 
-    const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+    const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
 
     const NodeId r0      = m_node_list[0];
     const CellId j0      = node_to_cell_matrix[r0][0];
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index 7e818bec46b799ff0acf113230db894890141d78..5d85b3630e2afc77ef7f2a3422b2f7d207abe5b7 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -199,7 +199,7 @@ class VTKWriter
  public:
   template <typename MeshType>
   void
-  write(const MeshType& mesh,
+  write(const std::shared_ptr<const MeshType>& mesh,
         const OutputNamedItemValueSet& output_named_item_value_set,
         double time,
         bool forced_output = false)
@@ -262,7 +262,7 @@ class VTKWriter
       fout << "<?xml version=\"1.0\"?>\n";
       fout << "<VTKFile type=\"UnstructuredGrid\">\n";
       fout << "<UnstructuredGrid>\n";
-      fout << "<Piece NumberOfPoints=\"" << mesh.numberOfNodes() << "\" NumberOfCells=\"" << mesh.numberOfCells()
+      fout << "<Piece NumberOfPoints=\"" << mesh->numberOfNodes() << "\" NumberOfCells=\"" << mesh->numberOfCells()
            << "\">\n";
       fout << "<CellData>\n";
       for (const auto& [name, item_value_variant] : output_named_item_value_set) {
@@ -279,10 +279,10 @@ class VTKWriter
       fout << "<Points>\n";
       {
         using Rd                      = TinyVector<MeshType::Dimension>;
-        const NodeValue<const Rd>& xr = mesh.xr();
-        Array<TinyVector<3>> positions(mesh.numberOfNodes());
+        const NodeValue<const Rd>& xr = mesh->xr();
+        Array<TinyVector<3>> positions(mesh->numberOfNodes());
         parallel_for(
-          mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+          mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
             for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
               positions[r][i] = xr[r][i];
             }
@@ -296,16 +296,16 @@ class VTKWriter
 
       fout << "<Cells>\n";
       {
-        const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+        const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
 
         _write_array(fout, "connectivity", cell_to_node_matrix.entries());
       }
 
       {
-        const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
-        Array<unsigned int> offsets(mesh.numberOfCells());
+        const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+        Array<unsigned int> offsets(mesh->numberOfCells());
         unsigned int offset = 0;
-        for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+        for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
           const auto& cell_nodes = cell_to_node_matrix[j];
           offset += cell_nodes.size();
           offsets[j] = offset;
@@ -314,10 +314,10 @@ class VTKWriter
       }
 
       {
-        Array<int8_t> types(mesh.numberOfCells());
-        const auto& cell_type = mesh.connectivity().cellType();
+        Array<int8_t> types(mesh->numberOfCells());
+        const auto& cell_type = mesh->connectivity().cellType();
         parallel_for(
-          mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+          mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
             switch (cell_type[j]) {
             case CellType::Line: {
               types[j] = 3;
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index fd10ea3203b02cf624cc0ff0b5a556517ba84512..0a5ef970d4a65242def2fe0fbad937e5fbab02ad 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -13,6 +13,7 @@
 
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
 
 #include <algebra/TinyMatrix.hpp>
 #include <algebra/TinyVector.hpp>
@@ -25,19 +26,18 @@
 
 #include <iostream>
 
-template <typename MeshData>
+template <typename MeshType>
 class AcousticSolver
 {
-  using MeshType     = typename MeshData::MeshType;
-  using UnknownsType = FiniteVolumesEulerUnknowns<MeshData>;
+  constexpr static size_t Dimension = MeshType::Dimension;
+
+  using MeshDataType = MeshData<Dimension>;
+  using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
 
-  MeshData& m_mesh_data;
-  const MeshType& m_mesh;
+  const std::shared_ptr<MeshType> m_mesh;
   const typename MeshType::Connectivity& m_connectivity;
   const std::vector<BoundaryConditionHandler>& m_boundary_condition_list;
 
-  constexpr static size_t Dimension = MeshType::Dimension;
-
   using Rd  = TinyVector<Dimension>;
   using Rdd = TinyMatrix<Dimension>;
 
@@ -47,7 +47,7 @@ class AcousticSolver
   computeRhoCj(const CellValue<const double>& rhoj, const CellValue<const double>& cj)
   {
     parallel_for(
-      m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { m_rhocj[j] = rhoj[j] * cj[j]; });
+      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { m_rhocj[j] = rhoj[j] * cj[j]; });
     return m_rhocj;
   }
 
@@ -59,7 +59,7 @@ class AcousticSolver
              const NodeValuePerCell<const Rd>& njr)
   {
     parallel_for(
-      m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
         const size_t& nb_nodes = m_Ajr.numberOfSubValues(j);
         const double& rho_c    = rhocj[j];
         for (size_t r = 0; r < nb_nodes; ++r) {
@@ -76,7 +76,7 @@ class AcousticSolver
     const auto& node_local_numbers_in_their_cells = m_connectivity.nodeLocalNumbersInTheirCells();
 
     parallel_for(
-      m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+      m_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
         Rdd sum                                    = zero;
         const auto& node_to_cell                   = node_to_cell_matrix[r];
         const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemValues(r);
@@ -103,7 +103,7 @@ class AcousticSolver
     const auto& node_local_numbers_in_their_cells = m_connectivity.nodeLocalNumbersInTheirCells();
 
     parallel_for(
-      m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+      m_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
         Rd& br                                     = m_br[r];
         br                                         = zero;
         const auto& node_to_cell                   = node_to_cell_matrix[r];
@@ -161,7 +161,7 @@ class AcousticSolver
     inverse(Ar, m_inv_Ar);
     const NodeValue<const Rdd> invAr = m_inv_Ar;
     parallel_for(
-      m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { m_ur[r] = invAr[r] * br[r]; });
+      m_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { m_ur[r] = invAr[r] * br[r]; });
 
     return m_ur;
   }
@@ -173,10 +173,10 @@ class AcousticSolver
              const CellValue<const Rd>& uj,
              const CellValue<const double>& pj)
   {
-    const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+    const auto& cell_to_node_matrix = m_mesh->connectivity().cellToNodeMatrix();
 
     parallel_for(
-      m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
         const auto& cell_nodes = cell_to_node_matrix[j];
 
         for (size_t r = 0; r < cell_nodes.size(); ++r) {
@@ -189,7 +189,7 @@ class AcousticSolver
   inverse(const NodeValue<const Rdd>& A, NodeValue<Rdd>& inv_A) const
   {
     parallel_for(
-      m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { inv_A[r] = ::inverse(A[r]); });
+      m_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { inv_A[r] = ::inverse(A[r]); });
   }
 
   PUGS_INLINE
@@ -229,10 +229,9 @@ class AcousticSolver
   CellValue<double> m_Vj_over_cj;
 
  public:
-  AcousticSolver(MeshData& mesh_data, const std::vector<BoundaryConditionHandler>& bc_list)
-    : m_mesh_data(mesh_data),
-      m_mesh(mesh_data.mesh()),
-      m_connectivity(m_mesh.connectivity()),
+  AcousticSolver(std::shared_ptr<MeshType>& p_mesh, const std::vector<BoundaryConditionHandler>& bc_list)
+    : m_mesh(p_mesh),
+      m_connectivity(m_mesh->connectivity()),
       m_boundary_condition_list(bc_list),
       m_br(m_connectivity),
       m_Ajr(m_connectivity),
@@ -250,11 +249,13 @@ class AcousticSolver
   double
   acoustic_dt(const CellValue<const double>& Vj, const CellValue<const double>& cj) const
   {
-    const NodeValuePerCell<const double>& ljr = m_mesh_data.ljr();
-    const auto& cell_to_node_matrix           = m_mesh.connectivity().cellToNodeMatrix();
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+    const NodeValuePerCell<const double>& ljr = mesh_data.ljr();
+    const auto& cell_to_node_matrix           = m_mesh->connectivity().cellToNodeMatrix();
 
     parallel_for(
-      m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
         const auto& cell_nodes = cell_to_node_matrix[j];
 
         double S = 0;
@@ -278,20 +279,22 @@ class AcousticSolver
     CellValue<double>& pj = unknowns.pj();
     CellValue<double>& cj = unknowns.cj();
 
-    const CellValue<const double>& Vj         = m_mesh_data.Vj();
-    const NodeValuePerCell<const Rd>& Cjr     = m_mesh_data.Cjr();
-    const NodeValuePerCell<const double>& ljr = m_mesh_data.ljr();
-    const NodeValuePerCell<const Rd>& njr     = m_mesh_data.njr();
+    const MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+    const CellValue<const double> Vj         = mesh_data.Vj();
+    const NodeValuePerCell<const Rd> Cjr     = mesh_data.Cjr();
+    const NodeValuePerCell<const double> ljr = mesh_data.ljr();
+    const NodeValuePerCell<const Rd> njr     = mesh_data.njr();
 
     computeExplicitFluxes(rhoj, uj, pj, cj, Cjr, ljr, njr);
 
     const NodeValuePerCell<Rd>& Fjr = m_Fjr;
     const NodeValue<const Rd> ur    = m_ur;
-    const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+    const auto& cell_to_node_matrix = m_mesh->connectivity().cellToNodeMatrix();
 
     const CellValue<const double> inv_mj = unknowns.invMj();
     parallel_for(
-      m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
         const auto& cell_nodes = cell_to_node_matrix[j];
 
         Rd momentum_fluxes   = zero;
@@ -306,16 +309,17 @@ class AcousticSolver
       });
 
     parallel_for(
-      m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { ej[j] = Ej[j] - 0.5 * (uj[j], uj[j]); });
+      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { ej[j] = Ej[j] - 0.5 * (uj[j], uj[j]); });
 
-    NodeValue<Rd> mutable_xr = m_mesh.mutableXr();
+    NodeValue<Rd> new_xr = m_mesh->mutableXr();
     parallel_for(
-      m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { mutable_xr[r] += dt * ur[r]; });
-    m_mesh_data.updateAllData();
+      m_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
+
+    MeshDataManager::instance().getMeshData(*m_mesh).updateAllData();
 
     const CellValue<const double> mj = unknowns.mj();
     parallel_for(
-      m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { rhoj[j] = mj[j] / Vj[j]; });
+      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { rhoj[j] = mj[j] / Vj[j]; });
   }
 };
 
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index a3fd012de3d0ea44fd8b0cfe724238bb115fa05c..406bc571c20729c49cfc26b90d436d438ebf3a15 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -16,7 +16,7 @@ class FiniteVolumesEulerUnknowns
   using Rd                          = TinyVector<Dimension>;
 
  private:
-  const MeshDataType& m_mesh_data;
+  MeshDataType& m_mesh_data;
   const MeshType& m_mesh;
 
   CellValue<double> m_rhoj;
@@ -182,7 +182,7 @@ class FiniteVolumesEulerUnknowns
       m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { m_inv_mj[j] = 1. / m_mj[j]; });
   }
 
-  FiniteVolumesEulerUnknowns(const MeshDataType& mesh_data)
+  FiniteVolumesEulerUnknowns(MeshDataType& mesh_data)
     : m_mesh_data(mesh_data),
       m_mesh(m_mesh_data.mesh()),
       m_rhoj(m_mesh.connectivity()),