diff --git a/src/scheme/DiscreteFunctionP0Vector.hpp b/src/scheme/DiscreteFunctionP0Vector.hpp
index 77ddf70f1d7f14a1c4dcab4cb352d4c6280d1715..f4ccc7d3c3395d5c2c070f7ab3464cd0561b9344 100644
--- a/src/scheme/DiscreteFunctionP0Vector.hpp
+++ b/src/scheme/DiscreteFunctionP0Vector.hpp
@@ -31,6 +31,7 @@ class DiscreteFunctionP0Vector : public IDiscreteFunction
   DiscreteFunctionDescriptorP0Vector m_discrete_function_descriptor;
 
  public:
+  PUGS_INLINE
   ASTNodeDataType
   dataType() const final
   {
@@ -44,18 +45,21 @@ class DiscreteFunctionP0Vector : public IDiscreteFunction
     return m_cell_arrays.sizeOfArrays();
   }
 
+  PUGS_INLINE
   const CellArray<DataType>&
   cellArrays() const
   {
     return m_cell_arrays;
   }
 
+  PUGS_INLINE
   std::shared_ptr<const IMesh>
   mesh() const
   {
     return m_mesh;
   }
 
+  PUGS_INLINE
   const IDiscreteFunctionDescriptor&
   descriptor() const final
   {
@@ -63,21 +67,75 @@ class DiscreteFunctionP0Vector : public IDiscreteFunction
   }
 
   PUGS_FORCEINLINE
-  Array<double>
+  operator DiscreteFunctionP0Vector<Dimension, const DataType>() const
+  {
+    return DiscreteFunctionP0Vector<Dimension, const DataType>(m_mesh, m_cell_arrays);
+  }
+
+  PUGS_INLINE
+  void
+  fill(const DataType& data) const noexcept
+  {
+    static_assert(not std::is_const_v<DataType>, "Cannot modify ItemValue of const");
+    m_cell_arrays.fill(data);
+  }
+
+  PUGS_INLINE DiscreteFunctionP0Vector
+  operator=(const DiscreteFunctionP0Vector& f)
+  {
+    Assert(m_mesh == f.m_mesh);
+    Assert(this->size() == f.size());
+    m_cell_arrays = f.m_cell_arrays;
+    return *this;
+  }
+
+  friend PUGS_INLINE DiscreteFunctionP0Vector<Dimension, std::remove_const_t<DataType>>
+  copy(const DiscreteFunctionP0Vector& source)
+  {
+    return DiscreteFunctionP0Vector<Dimension, std::remove_const_t<DataType>>{source.m_mesh, copy(source.cellArrays())};
+  }
+
+  friend PUGS_INLINE void
+  copy_to(const DiscreteFunctionP0Vector<Dimension, DataType>& source,
+          DiscreteFunctionP0Vector<Dimension, std::remove_const_t<DataType>>& destination)
+  {
+    Assert(source.m_mesh == destination.m_mesh);
+    copy_to(source.m_cell_arrays, destination.m_cell_arrays);
+  }
+
+  PUGS_FORCEINLINE
+  Array<DataType>
   operator[](const CellId cell_id) const noexcept(NO_ASSERT)
   {
     return m_cell_arrays[cell_id];
   }
 
-  friend DiscreteFunctionP0Vector
+  PUGS_INLINE DiscreteFunctionP0Vector<Dimension, std::remove_const_t<DataType>>
+  operator-() const
+  {
+    Assert(m_cell_arrays.isBuilt());
+    DiscreteFunctionP0Vector<Dimension, std::remove_const_t<DataType>> opposite = copy(*this);
+
+    const size_t size_of_arrays = this->size();
+    parallel_for(
+      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        auto array = opposite[cell_id];
+        for (size_t i = 0; i < size_of_arrays; ++i) {
+          array[i] = -array[i];
+        }
+      });
+
+    return opposite;
+  }
+
+  friend DiscreteFunctionP0Vector<Dimension, std::remove_const_t<DataType>>
   operator+(const DiscreteFunctionP0Vector& f, const DiscreteFunctionP0Vector& g)
   {
-    Assert(f.mesh() == g.mesh(), "functions are nor defined on the same mesh");
+    Assert(f.m_mesh == g.m_mesh, "functions are nor defined on the same mesh");
     Assert(f.size() == g.size(), "P0 vectors have different sizes");
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(f.mesh());
-    DiscreteFunctionP0Vector sum(mesh, f.size());
+    DiscreteFunctionP0Vector<Dimension, std::remove_const_t<DataType>> sum(f.m_mesh, f.size());
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
         for (size_t i = 0; i < f.size(); ++i) {
           sum[cell_id][i] = f[cell_id][i] + g[cell_id][i];
         }
@@ -85,15 +143,14 @@ class DiscreteFunctionP0Vector : public IDiscreteFunction
     return sum;
   }
 
-  friend DiscreteFunctionP0Vector
+  friend DiscreteFunctionP0Vector<Dimension, std::remove_const_t<DataType>>
   operator-(const DiscreteFunctionP0Vector& f, const DiscreteFunctionP0Vector& g)
   {
     Assert(f.mesh() == g.mesh(), "functions are nor defined on the same mesh");
     Assert(f.size() == g.size(), "P0 vectors have different sizes");
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(f.mesh());
-    DiscreteFunctionP0Vector difference(mesh, f.size());
+    DiscreteFunctionP0Vector<Dimension, std::remove_const_t<DataType>> difference(f.m_mesh, f.size());
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
         for (size_t i = 0; i < f.size(); ++i) {
           difference[cell_id][i] = f[cell_id][i] - g[cell_id][i];
         }
@@ -102,32 +159,33 @@ class DiscreteFunctionP0Vector : public IDiscreteFunction
   }
 
   template <typename DataType2>
-  friend DiscreteFunctionP0Vector
+  friend DiscreteFunctionP0Vector<Dimension, std::remove_const_t<DataType>>
   operator*(const DiscreteFunctionP0<Dimension, DataType2>& f, const DiscreteFunctionP0Vector& g)
   {
     static_assert(std::is_arithmetic_v<DataType2>, "left operand must be arithmetic");
 
     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());
-    DiscreteFunctionP0Vector product(mesh, g.size());
+    DiscreteFunctionP0Vector<Dimension, std::remove_const_t<DataType>> product(g.m_mesh, g.size());
+    const size_t size_of_arrays = g.size();
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
-        const double f_value = f[cell_id];
-        for (size_t i = 0; i < g.size(); ++i) {
+      g.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        const DataType2& f_value = f[cell_id];
+        for (size_t i = 0; i < size_of_arrays; ++i) {
           product[cell_id][i] = f_value * g[cell_id][i];
         }
       });
     return product;
   }
 
-  friend DiscreteFunctionP0Vector
-  operator*(double a, const DiscreteFunctionP0Vector& f)
+  template <typename DataType2>
+  friend DiscreteFunctionP0Vector<Dimension, std::remove_const_t<DataType>>
+  operator*(const DataType2& a, const DiscreteFunctionP0Vector& f)
   {
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(f.mesh());
-    DiscreteFunctionP0Vector product(mesh, f.size());
+    DiscreteFunctionP0Vector<Dimension, std::remove_const_t<DataType>> product(f.m_mesh, f.size());
+    const size_t size_of_arrays = f.size();
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
-        for (size_t i = 0; i < f.size(); ++i) {
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        for (size_t i = 0; i < size_of_arrays; ++i) {
           product[cell_id][i] = a * f[cell_id][i];
         }
       });
@@ -149,16 +207,10 @@ class DiscreteFunctionP0Vector : public IDiscreteFunction
     : m_mesh{mesh}, m_cell_arrays{mesh->connectivity(), size}
   {}
 
-  template <typename DataType2>
-  DiscreteFunctionP0Vector(const std::shared_ptr<const MeshType>& mesh, const CellArray<DataType2>& cell_arrays)
-    : m_mesh{mesh}
+  DiscreteFunctionP0Vector(const std::shared_ptr<const MeshType>& mesh, const CellArray<DataType>& cell_arrays)
+    : m_mesh{mesh}, m_cell_arrays{cell_arrays}
   {
-    static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>);
-    static_assert(std::is_const_v<DataType> or not std::is_const_v<DataType2>);
-
-    Assert(mesh->connectivity_ptr() == cell_arrays.connectivity_ptr());
-
-    m_cell_arrays = cell_arrays;
+    Assert(mesh->shared_connectivity() == cell_arrays.connectivity_ptr());
   }
 
   DiscreteFunctionP0Vector(const DiscreteFunctionP0Vector&) noexcept = default;