Fix MeshData and MeshDataManager

MeshData constructor can now only be called by MeshDataManager! This
ensures only one construction of the data storage for each mesh

Also fix following issues
- volume calculation and checking was dangerously nested
- few iterators were poorly named (due to copy/paste in MeshDataManager)

src/language/modules/SchemeModule.cpp

@@ -116,7 +116,7 @@ struct GlaceScheme
     UnknownsType unknowns(*m_mesh);
 
     {
-      MeshDataType mesh_data(*m_mesh);
+      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
 
       unknowns.rhoj() =
         InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(rho_id,
@@ -152,7 +152,7 @@ struct GlaceScheme
     block_eos.updateEandCFromRhoP();
 
     {
-      MeshDataType mesh_data(*m_mesh);
+      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
 
       const CellValue<const double> Vj = mesh_data.Vj();
 
@@ -169,7 +169,7 @@ struct GlaceScheme
     VTKWriter vtk_writer("mesh_" + std::to_string(Dimension), 0.01);
 
     while ((t < tmax) and (iteration < itermax)) {
-      MeshDataType mesh_data(*m_mesh);
+      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
 
       vtk_writer.write(m_mesh,
                        {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
@@ -198,7 +198,7 @@ struct GlaceScheme
               << rang::fgB::cyan << iteration << rang::style::reset << rang::style::bold << " iterations"
               << rang::style::reset << '\n';
     {
-      MeshDataType mesh_data(*m_mesh);
+      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
 
       vtk_writer.write(m_mesh,
                        {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},

src/mesh/MeshData.hpp

@@ -73,7 +73,6 @@ class MeshData : public IMeshData
         });
       m_xj = std::make_shared<CellValue<const Rd>>(xj);
     }
-    std::cout << "- computed xj for mesh " << &m_mesh << '\n';
   }
 
   PUGS_INLINE
@@ -97,8 +96,6 @@ class MeshData : public IMeshData
         Vj[j] = inv_Dimension * sum_cjr_xr;
       });
     m_Vj = std::make_shared<CellValue<const double>>(Vj);
-
-    std::cout << "- computed Vj for mesh " << &m_mesh << '\n';
   }
 
   PUGS_INLINE
@@ -275,7 +272,8 @@ class MeshData : public IMeshData
   void
   _checkCellVolume()
   {
-    auto Vj = this->Vj();
+    Assert(m_Vj);
+    auto& Vj = *m_Vj;
 
     bool is_valid = [&] {
       for (CellId j = 0; j < m_mesh.numberOfCells(); ++j) {
@@ -360,17 +358,18 @@ class MeshData : public IMeshData
     return *m_Vj;
   }
 
+ private:
+  // MeshData **must** be constructed through MeshDataManager
+  friend class MeshDataManager;
   MeshData(const MeshType& mesh) : m_mesh(mesh) {}
 
+ public:
   MeshData() = delete;
 
-  MeshData(const MeshData&) = default;
-  MeshData(MeshData&&)      = default;
+  MeshData(const MeshData&) = delete;
+  MeshData(MeshData&&)      = delete;
 
-  ~MeshData()
-  {
-    ;
-  }
+  ~MeshData() {}
 };
 
 #endif   // MESH_DATA_HPP

src/mesh/MeshDataManager.cpp

@@ -25,8 +25,8 @@ MeshDataManager::destroy()
   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';
+    for (const auto& i_mesh_data : m_instance->m_mesh_mesh_data_map) {
+      error << " - mesh data " << rang::fgB::magenta << i_mesh_data.first << rang::style::reset << '\n';
     }
     throw UnexpectedError(error.str());
   }
@@ -46,11 +46,12 @@ MeshDataManager::getMeshData(const Mesh<Connectivity<Dimension>>& 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 dynamic_cast<MeshData<Dimension>&>(*connectivity_synchronizer->second);
+  if (auto i_mesh_data = m_mesh_mesh_data_map.find(p_mesh); i_mesh_data != m_mesh_mesh_data_map.end()) {
+    return dynamic_cast<MeshData<Dimension>&>(*i_mesh_data->second);
   } else {
-    std::shared_ptr mesh_data    = std::make_shared<MeshData<Dimension>>(mesh);
+    // **cannot** use make_shared since MeshData constructor is **private**
+    std::shared_ptr<MeshData<Dimension>> mesh_data{new MeshData<Dimension>(mesh)};
+
     m_mesh_mesh_data_map[p_mesh] = mesh_data;
     return *mesh_data;
   }