diff --git a/src/Makefile.am b/src/Makefile.am
index 6e924008f481947a41f5659c85078610b68f963a..b97aa8cea74a9d4ed0d91259f051e49790b4a121 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -23,7 +23,7 @@ run-one-body-sample-gdb: one-body-sample
 	libtool --mode execute gdb ./one-body-sample
 
 one_body_sample_SOURCES = \
-	one-body-sample.cc assemblers.cc LambertW.cc compute_state.cc vtk.cc
+	one-body-sample.cc assemblers.cc LambertW.cc compute_state.cc vtk.cc print.cc
 one_body_sample_CPPFLAGS = \
 	$(AM_CPPFLAGS) -Dsrcdir=\"$(srcdir)\"
 one_body_sample_LDADD = \
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 4d08d013d4fd0ff028035fe68613f70cb2f59ad6..06ec0ceacc542da23c4be4b77b6517e1a252747c 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -64,6 +64,7 @@
 
 #include "assemblers.hh"
 #include "compute_state.hh"
+#include "print.hh"
 #include "vtk.hh"
 
 int const dim = 2;
@@ -310,23 +311,10 @@ int main(int argc, char *argv[]) {
           }
         }
 
-        if (parset.get<bool>("printEvolution")) {
-          // Find the first node that belongs to the frictional boundary
-          for (size_t i = 0; i < frictionalNodes.size(); ++i) {
-            if (!frictionalNodes[i][0])
-              continue;
-
-            boost::format const formatter("s[%03d] = %+3e, "
-                                          "%|40t|u[%03d] = %+3e");
-            std::cout << boost::format(formatter) % run % (*s4_new)[i] % run %
-                             u4[i] << std::endl;
-            double out;
-            functions.get("sampleFunction").evaluate(h * run, out);
-            octave_writer << (*s4_new)[i] << " " << u4[i][0] * 1e6 << " " << out
-                          << std::endl;
-            break;
-          }
-        }
+        if (parset.get<bool>("printEvolution"))
+          print_evolution<SingletonVectorType, VectorType>(
+              frictionalNodes, *s4_new, u4, functions.get("sampleFunction"),
+              run, h * run, octave_writer);
       }
 
       u4 += u4_diff;
diff --git a/src/print.cc b/src/print.cc
new file mode 100644
index 0000000000000000000000000000000000000000..4441bf3d18ce96afac746e537c9c2175308c2488
--- /dev/null
+++ b/src/print.cc
@@ -0,0 +1,30 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <boost/format.hpp>
+
+#include "print.hh"
+
+template <class SingletonVectorType, class VectorType>
+void print_evolution(Dune::BitSetVector<1> const &frictionalNodes,
+                     SingletonVectorType const &state,
+                     VectorType const &displacement,
+                     Dune::VirtualFunction<double, double> const &function,
+                     int run, double time, std::ostream &octave_writer) {
+  // Find the first node that belongs to the frictional boundary
+  for (size_t i = 0; i < frictionalNodes.size(); ++i) {
+    if (!frictionalNodes[i][0])
+      continue;
+
+    double out;
+    function.evaluate(time, out);
+    std::cout << boost::format("s[%03d] = %+3e, %|20t|u[%03d] = %+3e") % run %
+                     state[i] % run % displacement[i] << std::endl;
+    octave_writer << state[i][0] << " " << displacement[i][0] * 1e6 << " "
+                  << out << std::endl;
+    break;
+  }
+}
+
+#include "print_tmpl.cc"
diff --git a/src/print.hh b/src/print.hh
new file mode 100644
index 0000000000000000000000000000000000000000..79d52066569ca18f239593e10fb958b4a5cf5659
--- /dev/null
+++ b/src/print.hh
@@ -0,0 +1,14 @@
+#ifndef PRINT_HH
+#define PRINT_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/bitsetvector.hh>
+
+template <class SingletonVectorType, class VectorType>
+void print_evolution(Dune::BitSetVector<1> const &frictionalNodes,
+                     SingletonVectorType const &state,
+                     VectorType const &displacement,
+                     Dune::VirtualFunction<double, double> const &function,
+                     int run, double time, std::ostream &octave_writer);
+
+#endif
diff --git a/src/print_tmpl.cc b/src/print_tmpl.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ae475a813ce56d01705e68e30aeae58b9b749e6a
--- /dev/null
+++ b/src/print_tmpl.cc
@@ -0,0 +1,26 @@
+#include <dune/common/fvector.hh>
+#include <dune/istl/bvector.hh>
+
+typedef Dune::BlockVector<Dune::FieldVector<double, 1>> SingletonVectorType;
+
+// {{{ 2D
+typedef Dune::FieldVector<double, 2> SmallVector2;
+typedef Dune::BlockVector<SmallVector2> VectorType2;
+
+template void print_evolution<SingletonVectorType, VectorType2>(
+    Dune::BitSetVector<1> const &frictionalNodes,
+    SingletonVectorType const &state, VectorType2 const &displacement,
+    Dune::VirtualFunction<double, double> const &function, int run, double time,
+    std::ostream &octave_writer);
+// }}}
+
+// {{{ 3D
+typedef Dune::FieldVector<double, 3> SmallVector3;
+typedef Dune::BlockVector<SmallVector3> VectorType3;
+
+template void print_evolution<SingletonVectorType, VectorType3>(
+    Dune::BitSetVector<1> const &frictionalNodes,
+    SingletonVectorType const &state, VectorType3 const &displacement,
+    Dune::VirtualFunction<double, double> const &function, int run, double time,
+    std::ostream &octave_writer);
+// }}}