diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 16d1b550921be0d7a307924583c7b9f22f70136c..06e241bf3f59e4f04f8a3e553fbc81a83cf3cf76 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -63,7 +63,7 @@ struct GlaceScheme
 {
   using ConnectivityType = Connectivity<Dimension>;
   using MeshType         = Mesh<ConnectivityType>;
-  using MeshDataType     = MeshData<MeshType>;
+  using MeshDataType     = MeshDataLegacy<MeshType>;
   using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
 
   const MeshType& m_mesh;
diff --git a/src/main.cpp b/src/main.cpp
index 7224ff3965d7ae986c3adb7999a158e7a410d96b..2fabfe269a9a88ad13ef8965cdf144d33ad044a6 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -78,7 +78,7 @@ main(int argc, char* argv[])
 
       using ConnectivityType = Connectivity1D;
       using MeshType         = Mesh<ConnectivityType>;
-      using MeshDataType     = MeshData<MeshType>;
+      using MeshDataType     = MeshDataLegacy<MeshType>;
       using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
 
       const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
@@ -197,7 +197,7 @@ main(int argc, char* argv[])
 
       using ConnectivityType = Connectivity2D;
       using MeshType         = Mesh<ConnectivityType>;
-      using MeshDataType     = MeshData<MeshType>;
+      using MeshDataType     = MeshDataLegacy<MeshType>;
       using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
 
       const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
@@ -303,7 +303,7 @@ main(int argc, char* argv[])
 
       using ConnectivityType = Connectivity3D;
       using MeshType         = Mesh<ConnectivityType>;
-      using MeshDataType     = MeshData<MeshType>;
+      using MeshDataType     = MeshDataLegacy<MeshType>;
       using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
 
       const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 657b3b37500304446627038514ede886ad205da4..2a7853937a3821dc7370b923f6a5d05a23268f8a 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -12,6 +12,7 @@ add_library(
   GmshReader.cpp
   LogicalConnectivityBuilder.cpp
   MeshBuilderBase.cpp
+  MeshDataManager.cpp
   SynchronizerManager.cpp)
 
 # Additional dependencies
diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualMeshBuilder.cpp
index c5416939c4c81dd67114b30a720ee8a2a258a0c8..c0141d7015d600dcb1a4f06bdf2a7e6cf10b4590 100644
--- a/src/mesh/DiamondDualMeshBuilder.cpp
+++ b/src/mesh/DiamondDualMeshBuilder.cpp
@@ -26,7 +26,7 @@ DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(
   NodeValue<TinyVector<Dimension>> diamond_xr{diamond_connectivity};
 
   const NodeValue<const TinyVector<Dimension>> primal_xr = primal_mesh.xr();
-  MeshData<MeshType> primal_mesh_data{primal_mesh};
+  MeshDataLegacy<MeshType> primal_mesh_data{primal_mesh};
   const CellValue<const TinyVector<Dimension>> primal_xj = primal_mesh_data.xj();
 
   NodeId i_node = 0;
@@ -57,7 +57,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();
-  MeshData<MeshType> primal_mesh_data{primal_mesh};
+  MeshDataLegacy<MeshType> primal_mesh_data{primal_mesh};
   const CellValue<const TinyVector<1>> primal_xj = primal_mesh_data.xj();
 
   NodeId next_node_id = 0;
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 9fab2cf9aa8f782125ae21683634b1c6fff68dee..1b209e11dab47036fc1e7ee145b3b0b067e78aff 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -11,7 +11,7 @@
 #include <map>
 
 template <typename M>
-class MeshData
+class [[deprecated]] MeshDataLegacy
 {
  public:
   using MeshType = M;
@@ -32,8 +32,7 @@ class MeshData
   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();
@@ -68,8 +67,7 @@ class MeshData
   }
 
   PUGS_INLINE
-  void
-  _updateVolume()
+  void _updateVolume()
   {
     const NodeValue<const Rd>& xr   = m_mesh.xr();
     const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
@@ -89,8 +87,7 @@ class MeshData
   }
 
   PUGS_INLINE
-  void
-  _updateCjr()
+  void _updateCjr()
   {
     if constexpr (Dimension == 1) {
       // Cjr/njr/ljr are constant overtime
@@ -218,8 +215,7 @@ class MeshData
     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) {
@@ -237,49 +233,42 @@ class MeshData
 
  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;
   }
 
   PUGS_INLINE
-  const NodeValuePerCell<const double>&
-  ljr() const
+  const NodeValuePerCell<const double>& ljr() const
   {
     return m_ljr;
   }
 
   PUGS_INLINE
-  const NodeValuePerCell<const Rd>&
-  njr() const
+  const NodeValuePerCell<const Rd>& njr() const
   {
     return m_njr;
   }
 
   PUGS_INLINE
-  const CellValue<const Rd>&
-  xj() const
+  const CellValue<const Rd>& xj() const
   {
     return m_xj;
   }
 
   PUGS_INLINE
-  const CellValue<const double>&
-  Vj() const
+  const CellValue<const double>& Vj() const
   {
     return m_Vj;
   }
 
-  void
-  updateAllData()
+  void updateAllData()
   {
     this->_updateCjr();
     this->_updateCenter();
@@ -287,7 +276,7 @@ class MeshData
     this->_checkCellVolume();
   }
 
-  MeshData(const MeshType& mesh) : m_mesh(mesh)
+  MeshDataLegacy(const MeshType& mesh) : m_mesh(mesh)
   {
     if constexpr (Dimension == 1) {
       // in 1d Cjr are computed once for all
@@ -312,10 +301,10 @@ class MeshData
     this->updateAllData();
   }
 
-  MeshData(const MeshData&) = delete;
-  MeshData(MeshData&&)      = delete;
+  MeshDataLegacy(const MeshDataLegacy&) = delete;
+  MeshDataLegacy(MeshDataLegacy &&)     = delete;
 
-  ~MeshData()
+  ~MeshDataLegacy()
   {
     ;
   }
diff --git a/src/mesh/MeshDataManager.cpp b/src/mesh/MeshDataManager.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a72fda8307d2a20033c206d31998a44931cbf23c
--- /dev/null
+++ b/src/mesh/MeshDataManager.cpp
@@ -0,0 +1,55 @@
+#include <utils/PugsAssert.hpp>
+
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+
+class MeshData
+{
+ public:
+  MeshData() = default;
+};
+
+MeshDataManager* MeshDataManager::m_instance{nullptr};
+
+void
+MeshDataManager::create()
+{
+  Assert(m_instance == nullptr, "MeshDataManager is already created");
+  m_instance = new MeshDataManager;
+}
+
+void
+MeshDataManager::destroy()
+{
+  Assert(m_instance != nullptr, "MeshDataManager was not created!");
+
+  if (m_instance->m_mesh_mesh_data_map.size() > 0) {
+    std::stringstream error;
+    error << ": some mesh data is still registered\n";
+    for (const auto& i_connectivity_synchronizer : m_instance->m_mesh_mesh_data_map) {
+      error << " - mesh data " << rang::fgB::magenta << i_connectivity_synchronizer.first << rang::style::reset << '\n';
+    }
+    throw UnexpectedError(error.str());
+  }
+  delete m_instance;
+  m_instance = nullptr;
+}
+
+void
+MeshDataManager::deleteMeshData(const IMesh* mesh_data)
+{
+  m_mesh_mesh_data_map.erase(mesh_data);
+}
+
+MeshData&
+MeshDataManager::getMeshData(const IMesh* mesh)
+{
+  if (auto connectivity_synchronizer = m_mesh_mesh_data_map.find(mesh);
+      connectivity_synchronizer != m_mesh_mesh_data_map.end()) {
+    return (*connectivity_synchronizer->second);
+  } else {
+    std::shared_ptr mesh_data  = std::make_shared<MeshData>();
+    m_mesh_mesh_data_map[mesh] = mesh_data;
+    return *mesh_data;
+  }
+}
diff --git a/src/mesh/MeshDataManager.hpp b/src/mesh/MeshDataManager.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..66041bec3525d83de7adf8190f2702dd08544dea
--- /dev/null
+++ b/src/mesh/MeshDataManager.hpp
@@ -0,0 +1,42 @@
+#ifndef MESH_DATA_MANAGER_HPP
+#define MESH_DATA_MANAGER_HPP
+
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+
+#include <map>
+#include <memory>
+
+class IMesh;
+class MeshData;
+
+class MeshDataManager
+{
+ private:
+  std::map<const IMesh*, std::shared_ptr<MeshData>> m_mesh_mesh_data_map;
+
+  static MeshDataManager* m_instance;
+
+  MeshDataManager(const MeshDataManager&) = delete;
+  MeshDataManager(MeshDataManager&&)      = delete;
+
+  MeshDataManager()  = default;
+  ~MeshDataManager() = default;
+
+ public:
+  static void create();
+  static void destroy();
+
+  PUGS_INLINE
+  static MeshDataManager&
+  instance()
+  {
+    Assert(m_instance != nullptr, "MeshDataManager was not created!");
+    return *m_instance;
+  }
+
+  void deleteMeshData(const IMesh*);
+  MeshData& getMeshData(const IMesh*);
+};
+
+#endif   // MESH_DATA_MANAGER_HPP