diff --git a/src/Makefile.am b/src/Makefile.am
index 89b6e40d91526597e9a4c465dc04fe0ef693cc66..2e35b6deb2bfd344b5fa87baddd29a4b16b5e9c3 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -7,6 +7,7 @@ bin_PROGRAMS = \
 
 SOURCES = \
 	assemblers.cc \
+	friction_writer.cc \
 	state/compute_state_dieterich_euler.cc \
 	state/compute_state_ruina.cc \
 	solverfactory.cc \
diff --git a/src/friction_writer.cc b/src/friction_writer.cc
new file mode 100644
index 0000000000000000000000000000000000000000..b9087fdc47e21ee15168d7d615126e161724c538
--- /dev/null
+++ b/src/friction_writer.cc
@@ -0,0 +1,52 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "friction_writer.hh"
+
+double computeCOF(FrictionData const &fd, double V, double log_state) {
+  double const mu = fd.mu0 + fd.a * std::log(V / fd.V0) +
+                    fd.b * (log_state + std::log(fd.V0 / fd.L));
+  return (mu > 0) ? mu : 0;
+}
+
+template <class BitVectorType>
+FrictionWriter<BitVectorType>::FrictionWriter(
+    FrictionData const &_fd, BitVectorType const &_boundaryNodes)
+    : coefficientWriter("coefficients", std::fstream::out),
+      displacementWriter("displacements", std::fstream::out),
+      stateWriter("states", std::fstream::out),
+      velocityWriter("velocities", std::fstream::out),
+      fd(_fd),
+      boundaryNodes(_boundaryNodes) {}
+
+template <class BitVectorType>
+FrictionWriter<BitVectorType>::~FrictionWriter() {
+  stateWriter.close();
+  displacementWriter.close();
+  velocityWriter.close();
+  coefficientWriter.close();
+}
+
+template <class BitVectorType>
+template <class SingletonVectorType, class VectorType>
+void FrictionWriter<BitVectorType>::writeInfo(SingletonVectorType const &alpha,
+                                              VectorType const &u,
+                                              VectorType const &v) {
+  for (size_t i = 0; i < boundaryNodes.size(); ++i) {
+    if (!boundaryNodes[i][0])
+      continue;
+
+    coefficientWriter << computeCOF(fd, v[i].two_norm(), alpha[i]) << " ";
+    displacementWriter << u[i][0] << " ";
+    stateWriter << alpha[i][0] << " ";
+    velocityWriter << v[i][0] << " ";
+  }
+
+  stateWriter << std::endl;
+  displacementWriter << std::endl;
+  velocityWriter << std::endl;
+  coefficientWriter << std::endl;
+}
+
+#include "friction_writer_tmpl.cc"
diff --git a/src/friction_writer.hh b/src/friction_writer.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a4902e228d344428d045fefbc68e7970e3c9496b
--- /dev/null
+++ b/src/friction_writer.hh
@@ -0,0 +1,23 @@
+//#include <iostream>
+#include <fstream>
+
+#include <dune/tectonic/frictiondata.hh>
+
+template <class BitVectorType> class FrictionWriter {
+public:
+  FrictionWriter(FrictionData const &_fd, BitVectorType const &_boundaryNodes);
+
+  ~FrictionWriter();
+
+  template <class SingletonVectorType, class VectorType>
+  void writeInfo(SingletonVectorType const &alpha, VectorType const &u,
+                 VectorType const &v);
+
+private:
+  std::fstream coefficientWriter;
+  std::fstream displacementWriter;
+  std::fstream stateWriter;
+  std::fstream velocityWriter;
+  FrictionData const &fd;
+  BitVectorType const &boundaryNodes;
+};
diff --git a/src/friction_writer_tmpl.cc b/src/friction_writer_tmpl.cc
new file mode 100644
index 0000000000000000000000000000000000000000..b37ab8c01e9a7c2698ce5708aab4a5e7a5995df8
--- /dev/null
+++ b/src/friction_writer_tmpl.cc
@@ -0,0 +1,16 @@
+#include <dune/common/bitsetvector.hh>
+#include <dune/istl/bvector.hh>
+
+using BitVectorType = Dune::BitSetVector<1>;
+using SingletonVectorType = Dune::BlockVector<Dune::FieldVector<double, 1>>;
+using VectorType2 = Dune::BlockVector<Dune::FieldVector<double, 2>>;
+using VectorType3 = Dune::BlockVector<Dune::FieldVector<double, 3>>;
+
+template class FrictionWriter<BitVectorType>;
+
+template void FrictionWriter<BitVectorType>::writeInfo(
+    SingletonVectorType const &alpha, VectorType2 const &u,
+    VectorType2 const &v);
+template void FrictionWriter<BitVectorType>::writeInfo(
+    SingletonVectorType const &alpha, VectorType3 const &u,
+    VectorType3 const &v);
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 22f432e9841b0742908e53f48d78bfa7788c933a..f0aa09dd5186bb8950fd92ee9962c81a465dfaa0 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -82,6 +82,7 @@
 #include <dune/tectonic/globalnonlinearity.hh>
 
 #include "assemblers.hh"
+#include "friction_writer.hh"
 #include "solverfactory.hh"
 #include "vtk.hh"
 
@@ -152,12 +153,6 @@ void initPython() {
   Python::run("sys.path.append('" srcdir "')");
 }
 
-double computeCOF(FrictionData const &fd, double V, double log_state) {
-  double const mu = fd.mu0 + fd.a * std::log(V / fd.V0) +
-                    fd.b * (log_state + std::log(fd.V0 / fd.L));
-  return (mu > 0) ? mu : 0;
-}
-
 int main(int argc, char *argv[]) {
   try {
     Dune::Timer timer;
@@ -442,6 +437,8 @@ int main(int argc, char *argv[]) {
       initialAccelerationProblemSolver.solve();
     }
     // }}}
+    FrictionWriter<Dune::BitSetVector<1>> writer(frictionData, frictionalNodes);
+    writer.writeInfo(alpha_initial, u_initial, v_initial);
 
     // Set up TNNMG solver
     using NonlinearFactoryType = SolverFactory<
@@ -460,11 +457,7 @@ int main(int argc, char *argv[]) {
           coordinateWriter << coordinates[i] << std::endl;
       coordinateWriter.close();
     }
-    std::fstream stateWriter("states", std::fstream::out),
-        displacementWriter("displacements", std::fstream::out),
-        velocityWriter("velocities", std::fstream::out),
-        coefficientWriter("coefficients", std::fstream::out),
-        iterationWriter("iterations", std::fstream::out),
+    std::fstream iterationWriter("iterations", std::fstream::out),
         relaxationWriter("relaxation", std::fstream::out);
 
     auto timeSteppingScheme =
@@ -579,19 +572,7 @@ int main(int argc, char *argv[]) {
       if (printProgress)
         std::cerr << std::endl;
 
-      for (size_t i = 0; i < frictionalNodes.size(); ++i)
-        if (frictionalNodes[i][0]) {
-          stateWriter << alpha[i][0] << " ";
-          displacementWriter << u[i][0] << " ";
-          velocityWriter << v[i][0] << " ";
-          coefficientWriter << computeCOF(frictionData, v[i].two_norm(),
-                                          alpha[i]) << " ";
-        }
-
-      stateWriter << std::endl;
-      displacementWriter << std::endl;
-      velocityWriter << std::endl;
-      coefficientWriter << std::endl;
+      writer.writeInfo(alpha, u, v);
       iterationWriter << std::endl;
       relaxationWriter << std::endl;
 
@@ -612,10 +593,6 @@ int main(int argc, char *argv[]) {
       std::cerr << std::endl << "Making " << timesteps << " time steps took "
                 << timer.elapsed() << "s" << std::endl;
 
-    stateWriter.close();
-    displacementWriter.close();
-    velocityWriter.close();
-    coefficientWriter.close();
     iterationWriter.close();
     relaxationWriter.close();