diff --git a/src/language/modules/VTKModule.cpp b/src/language/modules/VTKModule.cpp
index ce84141ee866f7ac6fdbb8448db78d0d1b64de88..3144460fb5729ba29db2f20f864d6a5552e74a9b 100644
--- a/src/language/modules/VTKModule.cpp
+++ b/src/language/modules/VTKModule.cpp
@@ -6,10 +6,11 @@
 #include <mesh/GmshReader.hpp>
 #include <mesh/Mesh.hpp>
 #include <output/VTKWriter.hpp>
+#include <scheme/IDiscreteFunction.hpp>
 
 VTKModule::VTKModule()
 {
-  this->_addBuiltinFunction("writeVTK",
+  this->_addBuiltinFunction("writeMeshVTK",
                             std::make_shared<
                               BuiltinFunctionEmbedder<void(std::shared_ptr<const IMesh>, const std::string&)>>(
 
@@ -57,5 +58,65 @@ VTKModule::VTKModule()
                                 time++;
                               }
 
+                              ));
+
+  this->_addBuiltinFunction("writeVTK",
+                            std::make_shared<BuiltinFunctionEmbedder<
+                              void(const std::vector<std::shared_ptr<const IDiscreteFunction>>&, const std::string&)>>(
+
+                              [](const std::vector<std::shared_ptr<const IDiscreteFunction>>& discrete_function_list,
+                                 const std::string& filename) -> void {
+                                Assert(discrete_function_list.size() > 0);
+
+                                VTKWriter writer(filename, 0.1);
+
+                                static double time = 0;
+
+                                std::shared_ptr p_mesh = discrete_function_list[0]->mesh();
+                                for (size_t i = 1; i < discrete_function_list.size(); ++i) {
+                                  if (p_mesh != discrete_function_list[i]->mesh()) {
+                                    throw NormalError("discrete functions must be defined on the same mesh!");
+                                  }
+                                }
+
+                                switch (p_mesh->dimension()) {
+                                case 1: {
+                                  using MeshType = Mesh<Connectivity<1>>;
+                                  const std::shared_ptr<const MeshType> mesh =
+                                    std::dynamic_pointer_cast<const MeshType>(p_mesh);
+
+                                  writer.write(mesh,
+                                               OutputNamedItemValueSet{
+                                                 NamedItemValue{"cell_number", mesh->connectivity().cellNumber()}},
+                                               time, true);
+                                  break;
+                                }
+                                case 2: {
+                                  using MeshType = Mesh<Connectivity<2>>;
+                                  const std::shared_ptr<const MeshType> mesh =
+                                    std::dynamic_pointer_cast<const MeshType>(p_mesh);
+
+                                  writer.write(mesh,
+                                               OutputNamedItemValueSet{
+                                                 NamedItemValue{"cell_number", mesh->connectivity().cellNumber()}},
+                                               time, true);
+                                  break;
+                                }
+                                case 3: {
+                                  using MeshType = Mesh<Connectivity<3>>;
+                                  const std::shared_ptr<const MeshType> mesh =
+                                    std::dynamic_pointer_cast<const MeshType>(p_mesh);
+
+                                  writer.write(mesh,
+                                               OutputNamedItemValueSet{
+                                                 NamedItemValue{"cell_number", mesh->connectivity().cellNumber()}},
+                                               time, true);
+                                  break;
+                                }
+                                }
+
+                                time++;
+                              }
+
                               ));
 }
diff --git a/src/scheme/DiscreteFunctionInterpoler.cpp b/src/scheme/DiscreteFunctionInterpoler.cpp
index 07e113b986b4dc8bda4cb1a902d27e7b0a938e5c..60b481cdeffcbae4571c5a297a5521099d01686d 100644
--- a/src/scheme/DiscreteFunctionInterpoler.cpp
+++ b/src/scheme/DiscreteFunctionInterpoler.cpp
@@ -17,6 +17,12 @@ class DiscreteFunctionP0 : public IDiscreteFunction
   CellValue<DataType> m_cell_values;
 
  public:
+  std::shared_ptr<const IMesh>
+  mesh() const
+  {
+    return m_mesh;
+  }
+
   PUGS_FORCEINLINE
   DataType&
   operator[](const CellId& cell_id) const noexcept(NO_ASSERT)
diff --git a/src/scheme/DiscreteFunctionType.hpp b/src/scheme/DiscreteFunctionType.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..debf863499df64cd161ac0b7c20a1a5ce4657107
--- /dev/null
+++ b/src/scheme/DiscreteFunctionType.hpp
@@ -0,0 +1,9 @@
+#ifndef DISCRETE_FUNCTION_TYPE_HPP
+#define DISCRETE_FUNCTION_TYPE_HPP
+
+enum class DiscreteFunctionType
+{
+  P0
+};
+
+#endif   // DISCRETE_FUNCTION_TYPE_HPP
diff --git a/src/scheme/IDiscreteFunction.hpp b/src/scheme/IDiscreteFunction.hpp
index 3eea07af9530bc1f52213b4610fff028a7d9a306..eafea0e961c03b4acae7032bf96c3334e8048810 100644
--- a/src/scheme/IDiscreteFunction.hpp
+++ b/src/scheme/IDiscreteFunction.hpp
@@ -1,9 +1,17 @@
 #ifndef I_DISCRETE_FUNCTION_HPP
 #define I_DISCRETE_FUNCTION_HPP
 
+#include <scheme/DiscreteFunctionType.hpp>
+
+#include <memory>
+
+class IMesh;
+
 class IDiscreteFunction
 {
  public:
+  virtual std::shared_ptr<const IMesh> mesh() const = 0;
+
   IDiscreteFunction() = default;
 
   IDiscreteFunction(const IDiscreteFunction&) = default;
diff --git a/src/scheme/IDiscreteFunctionDescriptor.hpp b/src/scheme/IDiscreteFunctionDescriptor.hpp
index c2b795d6dc924a595d95f19de17c1e3da025683e..fcd9fdd1aa8e07973aa9d3503424035d6ce9c646 100644
--- a/src/scheme/IDiscreteFunctionDescriptor.hpp
+++ b/src/scheme/IDiscreteFunctionDescriptor.hpp
@@ -1,15 +1,10 @@
 #ifndef I_DISCRETE_FUNCTION_DESCRIPTOR_HPP
 #define I_DISCRETE_FUNCTION_DESCRIPTOR_HPP
 
-enum class DiscreteFunctionDescriptorType
-{
-  P0
-};
-
 class IDiscreteFunctionDescriptor
 {
  public:
-  virtual DiscreteFunctionDescriptorType type() const = 0;
+  virtual DiscreteFunctionType type() const = 0;
 
   IDiscreteFunctionDescriptor() noexcept = default;
 
@@ -23,10 +18,10 @@ class IDiscreteFunctionDescriptor
 class DiscreteFunctionDescriptorP0 : public IDiscreteFunctionDescriptor
 {
  public:
-  DiscreteFunctionDescriptorType
+  DiscreteFunctionType
   type() const final
   {
-    return DiscreteFunctionDescriptorType::P0;
+    return DiscreteFunctionType::P0;
   }
 
   DiscreteFunctionDescriptorP0() noexcept = default;