diff --git a/tests/test_DiscreteFunctionP0.cpp b/tests/test_DiscreteFunctionP0.cpp
index e445d86d9b30683b86a886be305cd1d3d64e5806..1f60554b2b64f60e3ce6a1f440de587662e95eb7 100644
--- a/tests/test_DiscreteFunctionP0.cpp
+++ b/tests/test_DiscreteFunctionP0.cpp
@@ -135,6 +135,82 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
     }
   }
 
+  SECTION("fill")
+  {
+    auto all_values_equal = [](const auto& f, const auto& g) {
+      size_t number_of_cells = f.cellValues().numberOfItems();
+      for (CellId cell_id = 0; cell_id < number_of_cells; ++cell_id) {
+        if (f[cell_id] != g) {
+          return false;
+        }
+      }
+      return true;
+    };
+
+    SECTION("1D")
+    {
+      std::shared_ptr mesh       = MeshDataBaseForTests::get().cartesianMesh1D();
+      constexpr size_t Dimension = 1;
+
+      DiscreteFunctionP0<Dimension, double> f{mesh};
+      f.fill(3);
+
+      REQUIRE(all_values_equal(f, 3));
+
+      DiscreteFunctionP0<Dimension, TinyVector<3>> v{mesh};
+      v.fill(TinyVector<3>{1, 2, 3});
+
+      REQUIRE(all_values_equal(v, TinyVector<3>{1, 2, 3}));
+
+      DiscreteFunctionP0<Dimension, TinyMatrix<3>> A{mesh};
+      A.fill(TinyMatrix<3>{1, 2, 3, 4, 5, 6, 7, 8, 9});
+
+      REQUIRE(all_values_equal(A, TinyMatrix<3>{1, 2, 3, 4, 5, 6, 7, 8, 9}));
+    }
+
+    SECTION("2D")
+    {
+      std::shared_ptr mesh       = MeshDataBaseForTests::get().cartesianMesh2D();
+      constexpr size_t Dimension = 2;
+
+      DiscreteFunctionP0<Dimension, double> f{mesh};
+      f.fill(3);
+
+      REQUIRE(all_values_equal(f, 3));
+
+      DiscreteFunctionP0<Dimension, TinyVector<3>> v{mesh};
+      v.fill(TinyVector<3>{1, 2, 3});
+
+      REQUIRE(all_values_equal(v, TinyVector<3>{1, 2, 3}));
+
+      DiscreteFunctionP0<Dimension, TinyMatrix<3>> A{mesh};
+      A.fill(TinyMatrix<3>{1, 2, 3, 4, 5, 6, 7, 8, 9});
+
+      REQUIRE(all_values_equal(A, TinyMatrix<3>{1, 2, 3, 4, 5, 6, 7, 8, 9}));
+    }
+
+    SECTION("3D")
+    {
+      std::shared_ptr mesh       = MeshDataBaseForTests::get().cartesianMesh3D();
+      constexpr size_t Dimension = 3;
+
+      DiscreteFunctionP0<Dimension, double> f{mesh};
+      f.fill(3);
+
+      REQUIRE(all_values_equal(f, 3));
+
+      DiscreteFunctionP0<Dimension, TinyVector<3>> v{mesh};
+      v.fill(TinyVector<3>{1, 2, 3});
+
+      REQUIRE(all_values_equal(v, TinyVector<3>{1, 2, 3}));
+
+      DiscreteFunctionP0<Dimension, TinyMatrix<3>> A{mesh};
+      A.fill(TinyMatrix<3>{1, 2, 3, 4, 5, 6, 7, 8, 9});
+
+      REQUIRE(all_values_equal(A, TinyMatrix<3>{1, 2, 3, 4, 5, 6, 7, 8, 9}));
+    }
+  }
+
   SECTION("binary operators")
   {
     SECTION("1D")