Issue/mesh data

src/language/modules/SchemeModule.cpp

@@ -64,7 +64,7 @@ struct GlaceScheme
   using ConnectivityType = Connectivity<Dimension>;
   using MeshType         = Mesh<ConnectivityType>;
   using MeshDataType     = MeshData<Dimension>;
-  using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
+  using UnknownsType     = FiniteVolumesEulerUnknowns<MeshType>;
 
   std::shared_ptr<const MeshType> m_mesh;
 
@@ -75,8 +75,6 @@ struct GlaceScheme
               const FunctionSymbolId& p_id)
     : m_mesh{std::dynamic_pointer_cast<const MeshType>(i_mesh)}
   {
-    MeshDataType mesh_data(*m_mesh);
-
     std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
 
     std::vector<BoundaryConditionHandler> bc_list;
@@ -115,24 +113,25 @@ struct GlaceScheme
       }
     }
 
-    UnknownsType unknowns(mesh_data);
+    UnknownsType unknowns(*m_mesh);
+
+    {
+      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
 
       unknowns.rhoj() =
-      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(rho_id, mesh_data.xj());
+        InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(rho_id,
+                                                                                                  mesh_data.xj());
 
       unknowns.pj() =
         InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(p_id, mesh_data.xj());
 
-    unknowns.uj() =
-      InterpolateItemValue<TinyVector<Dimension>(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(u_id,
-                                                                                                               mesh_data
-                                                                                                                 .xj());
+      unknowns.uj() = InterpolateItemValue<TinyVector<Dimension>(
+        TinyVector<Dimension>)>::template interpolate<ItemType::cell>(u_id, mesh_data.xj());
+    }
     unknowns.gammaj().fill(1.4);
 
     AcousticSolver acoustic_solver(m_mesh, bc_list);
 
-    const CellValue<const double>& Vj = mesh_data.Vj();
-
     const double tmax = 0.2;
     double t          = 0;
 
@@ -152,6 +151,11 @@ struct GlaceScheme
     BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
     block_eos.updateEandCFromRhoP();
 
+    {
+      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+      const CellValue<const double> Vj = mesh_data.Vj();
+
       parallel_for(
         m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { Ej[j] = ej[j] + 0.5 * (uj[j], uj[j]); });
 
@@ -160,17 +164,20 @@ struct GlaceScheme
 
       parallel_for(
         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)) {
+      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
       vtk_writer.write(m_mesh,
                        {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                        NamedItemValue{"coords", m_mesh->xr()},
+                        NamedItemValue{"coords", m_mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
                         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);
+      double dt = 0.4 * acoustic_solver.acoustic_dt(mesh_data.Vj(), cj);
       if (t + dt > tmax) {
         dt = tmax - t;
       }
@@ -190,13 +197,17 @@ struct GlaceScheme
     std::cout << rang::style::bold << "Final time=" << rang::fgB::green << t << rang::style::reset << " reached after "
               << rang::fgB::cyan << iteration << rang::style::reset << rang::style::bold << " iterations"
               << rang::style::reset << '\n';
+    {
+      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
       vtk_writer.write(m_mesh,
                        {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                      NamedItemValue{"coords", m_mesh->xr()},
+                        NamedItemValue{"coords", m_mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
                         NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()},
                         NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}},
                        t, true);   // forces last output
     }
+  }
 };
 
 SchemeModule::SchemeModule()