diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp
index e9c6559e5798a24abde01d17eeb6f0828364822d..ffcda127022b7b93eae812784f89785a9691657b 100644
--- a/src/scheme/DiscreteFunctionP0.hpp
+++ b/src/scheme/DiscreteFunctionP0.hpp
@@ -26,30 +26,35 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   DiscreteFunctionDescriptorP0 m_discrete_function_descriptor;
 
  public:
+  PUGS_INLINE
   ASTNodeDataType
   dataType() const final
   {
     return ast_node_data_type_from<DataType>;
   }
 
+  PUGS_INLINE
   const CellValue<DataType>&
   cellValues() const
   {
     return m_cell_values;
   }
 
+  PUGS_INLINE
   std::shared_ptr<const IMesh>
   mesh() const
   {
     return m_mesh;
   }
 
+  PUGS_INLINE
   const IDiscreteFunctionDescriptor&
   descriptor() const final
   {
     return m_discrete_function_descriptor;
   }
 
+  PUGS_FORCEINLINE
   operator DiscreteFunctionP0<Dimension, const DataType>() const
   {
     return DiscreteFunctionP0<Dimension, const DataType>(m_mesh, m_cell_values);
@@ -63,23 +68,46 @@ class DiscreteFunctionP0 : public IDiscreteFunction
     m_cell_values.fill(data);
   }
 
-  PUGS_FORCEINLINE
-  DataType&
+  friend PUGS_INLINE DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>>
+  copy(const DiscreteFunctionP0& source)
+  {
+    return DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>>{source.m_mesh, copy(source.cellValues())};
+  }
+
+  friend PUGS_INLINE void
+  copy_to(const DiscreteFunctionP0<Dimension, DataType>& source,
+          DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>>& destination)
+  {
+    Assert(source.m_mesh == destination.m_mesh);
+    copy_to(source.m_cell_values, destination.m_cell_values);
+  }
+
+  PUGS_FORCEINLINE DataType&
   operator[](const CellId cell_id) const noexcept(NO_ASSERT)
   {
     return m_cell_values[cell_id];
   }
 
+  PUGS_INLINE DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>>
+  operator-() const
+  {
+    Assert(m_cell_values.isBuilt());
+    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> opposite = copy(*this);
+    parallel_for(
+      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { opposite[cell_id] = -opposite[cell_id]; });
+
+    return opposite;
+  }
+
   template <typename DataType2T>
   PUGS_INLINE DiscreteFunctionP0<Dimension, decltype(DataType{} + DataType2T{})>
   operator+(const DiscreteFunctionP0<Dimension, DataType2T>& g) const
   {
     const DiscreteFunctionP0& f = *this;
-    Assert(f.mesh() == g.mesh(), "functions are nor defined on the same mesh");
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(f.mesh());
-    DiscreteFunctionP0<Dimension, decltype(DataType{} + DataType2T{})> sum(mesh);
+    Assert(f.mesh() == g.mesh(), "functions are not defined on the same mesh");
+    DiscreteFunctionP0<Dimension, decltype(DataType{} + DataType2T{})> sum(f.m_mesh);
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum[cell_id] = f[cell_id] + g[cell_id]; });
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum[cell_id] = f[cell_id] + g[cell_id]; });
     return sum;
   }
 
@@ -87,11 +115,10 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   PUGS_INLINE friend DiscreteFunctionP0<Dimension, decltype(LHSDataType{} + DataType{})>
   operator+(const LHSDataType& a, const DiscreteFunctionP0& g)
   {
-    std::shared_ptr mesh  = std::dynamic_pointer_cast<const MeshType>(g.mesh());
     using ProductDataType = decltype(LHSDataType{} + DataType{});
-    DiscreteFunctionP0<Dimension, ProductDataType> sum(mesh);
+    DiscreteFunctionP0<Dimension, ProductDataType> sum(g.m_mesh);
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum[cell_id] = a + g[cell_id]; });
+      g.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum[cell_id] = a + g[cell_id]; });
     return sum;
   }
 
@@ -99,11 +126,10 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   PUGS_INLINE friend DiscreteFunctionP0<Dimension, decltype(DataType{} + RHSDataType{})>
   operator+(const DiscreteFunctionP0& f, const RHSDataType& b)
   {
-    std::shared_ptr mesh  = std::dynamic_pointer_cast<const MeshType>(f.mesh());
     using ProductDataType = decltype(DataType{} + RHSDataType{});
-    DiscreteFunctionP0<Dimension, ProductDataType> sum(mesh);
+    DiscreteFunctionP0<Dimension, ProductDataType> sum(f.m_mesh);
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum[cell_id] = f[cell_id] + b; });
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum[cell_id] = f[cell_id] + b; });
     return sum;
   }
 
@@ -112,11 +138,10 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   operator-(const DiscreteFunctionP0<Dimension, DataType2T>& g) const
   {
     const DiscreteFunctionP0& f = *this;
-    Assert(f.mesh() == g.mesh(), "functions are nor defined on the same mesh");
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(f.mesh());
-    DiscreteFunctionP0<Dimension, decltype(DataType{} - DataType2T{})> difference(mesh);
+    Assert(f.mesh() == g.mesh(), "functions are not defined on the same mesh");
+    DiscreteFunctionP0<Dimension, decltype(DataType{} - DataType2T{})> difference(f.m_mesh);
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference[cell_id] = f[cell_id] - g[cell_id]; });
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference[cell_id] = f[cell_id] - g[cell_id]; });
     return difference;
   }
 
@@ -124,11 +149,10 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   PUGS_INLINE friend DiscreteFunctionP0<Dimension, decltype(LHSDataType{} - DataType{})>
   operator-(const LHSDataType& a, const DiscreteFunctionP0& g)
   {
-    std::shared_ptr mesh  = std::dynamic_pointer_cast<const MeshType>(g.mesh());
     using ProductDataType = decltype(LHSDataType{} - DataType{});
-    DiscreteFunctionP0<Dimension, ProductDataType> difference(mesh);
+    DiscreteFunctionP0<Dimension, ProductDataType> difference(g.m_mesh);
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference[cell_id] = a - g[cell_id]; });
+      g.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference[cell_id] = a - g[cell_id]; });
     return difference;
   }
 
@@ -136,11 +160,10 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   PUGS_INLINE friend DiscreteFunctionP0<Dimension, decltype(DataType{} - RHSDataType{})>
   operator-(const DiscreteFunctionP0& f, const RHSDataType& b)
   {
-    std::shared_ptr mesh  = std::dynamic_pointer_cast<const MeshType>(f.mesh());
     using ProductDataType = decltype(DataType{} - RHSDataType{});
-    DiscreteFunctionP0<Dimension, ProductDataType> difference(mesh);
+    DiscreteFunctionP0<Dimension, ProductDataType> difference(f.m_mesh);
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference[cell_id] = f[cell_id] - b; });
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference[cell_id] = f[cell_id] - b; });
     return difference;
   }
 
@@ -149,11 +172,10 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   operator*(const DiscreteFunctionP0<Dimension, DataType2T>& g) const
   {
     const DiscreteFunctionP0& f = *this;
-    Assert(f.mesh() == g.mesh(), "functions are nor defined on the same mesh");
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(f.mesh());
-    DiscreteFunctionP0<Dimension, decltype(DataType{} * DataType2T{})> product(mesh);
+    Assert(f.mesh() == g.mesh(), "functions are not defined on the same mesh");
+    DiscreteFunctionP0<Dimension, decltype(DataType{} * DataType2T{})> product(f.m_mesh);
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product[cell_id] = f[cell_id] * g[cell_id]; });
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product[cell_id] = f[cell_id] * g[cell_id]; });
     return product;
   }
 
@@ -161,11 +183,10 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   PUGS_INLINE friend DiscreteFunctionP0<Dimension, decltype(LHSDataType{} * DataType{})>
   operator*(const LHSDataType& a, const DiscreteFunctionP0& f)
   {
-    std::shared_ptr mesh  = std::dynamic_pointer_cast<const MeshType>(f.mesh());
     using ProductDataType = decltype(LHSDataType{} * DataType{});
-    DiscreteFunctionP0<Dimension, ProductDataType> product(mesh);
+    DiscreteFunctionP0<Dimension, ProductDataType> product(f.m_mesh);
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product[cell_id] = a * f[cell_id]; });
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product[cell_id] = a * f[cell_id]; });
     return product;
   }
 
@@ -173,11 +194,10 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   PUGS_INLINE friend DiscreteFunctionP0<Dimension, decltype(DataType{} * RHSDataType{})>
   operator*(const DiscreteFunctionP0& f, const RHSDataType& b)
   {
-    std::shared_ptr mesh  = std::dynamic_pointer_cast<const MeshType>(f.mesh());
     using ProductDataType = decltype(DataType{} * RHSDataType{});
-    DiscreteFunctionP0<Dimension, ProductDataType> product(mesh);
+    DiscreteFunctionP0<Dimension, ProductDataType> product(f.m_mesh);
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product[cell_id] = f[cell_id] * b; });
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product[cell_id] = f[cell_id] * b; });
     return product;
   }
 
@@ -186,11 +206,10 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   operator/(const DiscreteFunctionP0<Dimension, DataType2T>& g) const
   {
     const DiscreteFunctionP0& f = *this;
-    Assert(f.mesh() == g.mesh(), "functions are nor defined on the same mesh");
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(f.mesh());
-    DiscreteFunctionP0<Dimension, decltype(DataType{} / DataType2T{})> ratio(mesh);
+    Assert(f.mesh() == g.mesh(), "functions are not defined on the same mesh");
+    DiscreteFunctionP0<Dimension, decltype(DataType{} / DataType2T{})> ratio(f.m_mesh);
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio[cell_id] = f[cell_id] / g[cell_id]; });
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio[cell_id] = f[cell_id] / g[cell_id]; });
     return ratio;
   }
 
@@ -198,11 +217,10 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   PUGS_INLINE friend DiscreteFunctionP0<Dimension, decltype(LHSDataType{} / DataType{})>
   operator/(const LHSDataType& a, const DiscreteFunctionP0& f)
   {
-    std::shared_ptr mesh  = std::dynamic_pointer_cast<const MeshType>(f.mesh());
     using ProductDataType = decltype(LHSDataType{} / DataType{});
-    DiscreteFunctionP0<Dimension, ProductDataType> ratio(mesh);
+    DiscreteFunctionP0<Dimension, ProductDataType> ratio(f.m_mesh);
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio[cell_id] = a / f[cell_id]; });
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio[cell_id] = a / f[cell_id]; });
     return ratio;
   }
 
@@ -210,11 +228,10 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   PUGS_INLINE friend DiscreteFunctionP0<Dimension, decltype(DataType{} / RHSDataType{})>
   operator/(const DiscreteFunctionP0& f, const RHSDataType& b)
   {
-    std::shared_ptr mesh  = std::dynamic_pointer_cast<const MeshType>(f.mesh());
     using ProductDataType = decltype(DataType{} / RHSDataType{});
-    DiscreteFunctionP0<Dimension, ProductDataType> ratio(mesh);
+    DiscreteFunctionP0<Dimension, ProductDataType> ratio(f.m_mesh);
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio[cell_id] = f[cell_id] / b; });
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio[cell_id] = f[cell_id] / b; });
     return ratio;
   }