diff --git a/ProcessLib/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlowProcess.h
index c4db0b63676e4e9bc16c6f0f74797c24e45d938f..fb06b0dd0124399d49f1a6947ed5e674a5099f36 100644
--- a/ProcessLib/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlowProcess.h
@@ -17,8 +17,6 @@
 #include "AssemblerLib/LocalAssemblerBuilder.h"
 #include "AssemblerLib/LocalDataInitializer.h"
 
-#include "FileIO/VtkIO/VtuInterface.h"
-
 #include "GroundwaterFlowFEM.h"
 #include "Parameter.h"
 #include "Process.h"
@@ -169,57 +167,6 @@ public:
         return true;
     }
 
-    void post(std::string const& file_name) override
-    {
-        DBUG("Postprocessing GroundwaterFlowProcess.");
-
-        std::string const property_name = "Result";
-
-        // Get or create a property vector for results.
-        boost::optional<MeshLib::PropertyVector<double>&> result;
-        if (this->_mesh.getProperties().hasPropertyVector(property_name))
-        {
-            result = this->_mesh.getProperties().template
-                getPropertyVector<double>(property_name);
-        }
-        else
-        {
-            result = this->_mesh.getProperties().template
-                createNewPropertyVector<double>(property_name,
-                    MeshLib::MeshItemType::Node);
-            result->resize(this->_x->size());
-        }
-        assert(result && result->size() == this->_x->size());
-
-#ifdef USE_PETSC
-        std::unique_ptr<double[]> u(new double[this->_x->size()]);
-        this->_x->getGlobalVector(u.get());  // get the global solution
-
-        std::size_t const n = this->_mesh.getNNodes();
-        for (std::size_t i = 0; i < n; ++i)
-        {
-            MeshLib::Location const l(this->_mesh.getID(),
-                                      MeshLib::MeshItemType::Node, i);
-            auto const global_index = std::abs(  // 0 is the component id.
-                this->_local_to_global_index_map->getGlobalIndex(l, 0));
-            (*result)[i] = u[global_index];
-        }
-#else
-        // Copy result
-        for (std::size_t i = 0; i < this->_x->size(); ++i)
-            (*result)[i] = (*this->_x)[i];
-#endif
-
-        // Write output file
-        FileIO::VtuInterface vtu_interface(&this->_mesh, vtkXMLWriter::Binary, true);
-        vtu_interface.writeToFile(file_name);
-    }
-
-    void postTimestep(std::string const& file_name, const unsigned /*timestep*/) override
-    {
-        post(file_name);
-    }
-
     ~GroundwaterFlowProcess()
     {
         for (auto p : _local_assemblers)
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index fc86eb09239410556da60805634f067c2e64dc10..f5d20b033db8227103ce20a862e69a89e464cc62 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -20,6 +20,7 @@
 #include "AssemblerLib/LocalToGlobalIndexMap.h"
 #include "AssemblerLib/VectorMatrixAssembler.h"
 #include "BaseLib/ConfigTreeNew.h"
+#include "FileIO/VtkIO/VtuInterface.h"
 #include "MathLib/LinAlg/ApplyKnownSolution.h"
 #include "MathLib/LinAlg/SetMatrixSparsity.h"
 #include "MeshGeoToolsLib/MeshNodeSearcher.h"
@@ -64,9 +65,11 @@ public:
 
 	/// Postprocessing after solve().
 	/// The file_name is indicating the name of possible output file.
-	virtual void post(std::string const& file_name) = 0;
-	virtual void postTimestep(std::string const& file_name,
-	                          const unsigned timestep) = 0;
+	void postTimestep(std::string const& file_name, const unsigned /*timestep*/)
+	{
+		post();
+		output(file_name);
+	}
 
 	void initialize()
 	{
@@ -128,6 +131,8 @@ public:
 	}
 
 protected:
+	virtual void post() { };
+
 	/// Set linear solver options; called by the derived process which is
 	/// parsing the configuration.
 	void setLinearSolverOptions(BaseLib::ConfigTreeNew&& config)
@@ -231,6 +236,52 @@ private:
 		    *_local_to_global_index_map, _mesh));
 	}
 
+	void output(std::string const& file_name)
+	{
+		DBUG("Process output.");
+
+		std::string const property_name = "Result";
+
+		// Get or create a property vector for results.
+		boost::optional<MeshLib::PropertyVector<double>&> result;
+		if (_mesh.getProperties().hasPropertyVector(property_name))
+		{
+			result = _mesh.getProperties().template getPropertyVector<double>(
+			    property_name);
+		}
+		else
+		{
+			result =
+			    _mesh.getProperties().template createNewPropertyVector<double>(
+			        property_name, MeshLib::MeshItemType::Node);
+			result->resize(_x->size());
+		}
+		assert(result && result->size() == _x->size());
+
+#ifdef USE_PETSC
+		std::unique_ptr<double[]> u(new double[_x->size()]);
+		_x->getGlobalVector(u.get());  // get the global solution
+
+		std::size_t const n = _mesh.getNNodes();
+		for (std::size_t i = 0; i < n; ++i)
+		{
+			MeshLib::Location const l(_mesh.getID(),
+			                          MeshLib::MeshItemType::Node, i);
+			auto const global_index = std::abs(  // 0 is the component id.
+			    _local_to_global_index_map->getGlobalIndex(l, 0));
+			(*result)[i] = u[global_index];
+		}
+#else
+		// Copy result
+		for (std::size_t i = 0; i < _x->size(); ++i)
+			(*result)[i] = (*_x)[i];
+#endif
+
+		// Write output file
+		FileIO::VtuInterface vtu_interface(&_mesh, vtkXMLWriter::Binary, true);
+		vtu_interface.writeToFile(file_name);
+	}
+
 protected:
 	unsigned const _integration_order = 2;