Rework the way DiscreteFunctionP0 are constructed

This enforces coherency between data:
- the mesh
- the cell values defined on its connectivity

src/scheme/DiscreteFunctionInterpoler.cpp

@@ -7,49 +7,61 @@
 #include <mesh/MeshDataManager.hpp>
 #include <utils/Exceptions.hpp>
 
-template <typename DataType>
+template <size_t Dimension, typename DataType>
 class DiscreteFunctionP0 : public IDiscreteFunction
 {
  private:
-  CellValue<DataType> m_cell_value;
+  using MeshType = Mesh<Connectivity<Dimension>>;
 
- public:
-  DiscreteFunctionP0(const CellValue<DataType>& cell_value) : m_cell_value{cell_value} {}
-  DiscreteFunctionP0(CellValue<DataType>&& cell_value) : m_cell_value{std::move(cell_value)} {}
-  ~DiscreteFunctionP0() = default;
-};
+  std::shared_ptr<const MeshType> m_mesh;
+  CellValue<DataType> m_cell_values;
 
-template <size_t Dimension>
-std::shared_ptr<IDiscreteFunction>
-DiscreteFunctionInterpoler::_interpolate() const
+ public:
+  PUGS_FORCEINLINE
+  DataType&
+  operator[](const CellId& cell_id) const noexcept(NO_ASSERT)
   {
-  using MeshType       = Mesh<Connectivity<Dimension>>;
-  const MeshType& mesh = dynamic_cast<const MeshType&>(*m_mesh);
+    return m_cell_values[cell_id];
+  }
 
+  DiscreteFunctionP0(const std::shared_ptr<const MeshType>& mesh, const FunctionSymbolId& function_id) : m_mesh(mesh)
+  {
     using MeshDataType      = MeshData<Dimension>;
-  MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
-
-  return std::make_shared<DiscreteFunctionP0<double>>(
-    InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(m_function_id,
-                                                                                              mesh_data.xj()));
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+    m_cell_values =
+      InterpolateItemValue<DataType(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(function_id,
+                                                                                                  mesh_data.xj());
   }
 
+  DiscreteFunctionP0(const DiscreteFunctionP0&) noexcept = default;
+  DiscreteFunctionP0(DiscreteFunctionP0&&) noexcept      = default;
+
+  ~DiscreteFunctionP0() = default;
+};
+
 std::shared_ptr<IDiscreteFunction>
 DiscreteFunctionInterpoler::interpolate() const
 {
+  std::shared_ptr<IDiscreteFunction> discrete_function;
   switch (m_mesh->dimension()) {
   case 1: {
-    return this->_interpolate<1>();
+    std::shared_ptr mesh = std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(m_mesh);
+    discrete_function    = std::make_shared<DiscreteFunctionP0<1, double>>(mesh, m_function_id);
+    break;
   }
   case 2: {
-    return this->_interpolate<2>();
+    std::shared_ptr mesh = std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(m_mesh);
+    discrete_function    = std::make_shared<DiscreteFunctionP0<2, double>>(mesh, m_function_id);
+    break;
   }
   case 3: {
-    return this->_interpolate<3>();
+    std::shared_ptr mesh = std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(m_mesh);
+    discrete_function    = std::make_shared<DiscreteFunctionP0<3, double>>(mesh, m_function_id);
+    break;
   }
   default: {
     throw UnexpectedError("invalid dimension");
   }
   }
-  return nullptr;
+  return discrete_function;
 }