diff --git a/src/sand-wedge-data/special_writer.hh b/src/sand-wedge-data/special_writer.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a9f3adedd8a671f24716dfefe45d5d3daee743da
--- /dev/null
+++ b/src/sand-wedge-data/special_writer.hh
@@ -0,0 +1,84 @@
+#ifndef SPECIAL_WRITER_HH
+#define SPECIAL_WRITER_HH
+
+#include <fstream>
+#include <utility>
+
+#include <dune/common/fvector.hh>
+#include <dune/grid/utility/hierarchicsearch.hh>
+
+#include "mygeometry.hh"
+
+template <class GridView, int dimension> class SpecialWriter {
+  using LocalVector = Dune::FieldVector<double, dimension>;
+  using Element = typename GridView::Grid::template Codim<0>::Entity;
+  using ElementPointer =
+      typename GridView::Grid::template Codim<0>::EntityPointer;
+
+  void writeHorizontal(LocalVector const &v) {
+    writer_ << MyGeometry::horizontalProjection(v) << " ";
+  }
+  void writeVertical(LocalVector const &v) {
+    writer_ << MyGeometry::verticalProjection(v) << " ";
+  }
+
+  std::pair<ElementPointer, LocalVector> globalToLocal(LocalVector const &x)
+      const {
+    auto const element = hsearch_.findEntity(x);
+    return std::pair<ElementPointer, LocalVector>(element,
+                                                  element->geometry().local(x));
+  }
+
+  std::fstream writer_;
+
+  Dune::HierarchicSearch<typename GridView::Grid, GridView> const hsearch_;
+
+  std::pair<ElementPointer, LocalVector> const G;
+  std::pair<ElementPointer, LocalVector> const H;
+  std::pair<ElementPointer, LocalVector> const J;
+  std::pair<ElementPointer, LocalVector> const I;
+  std::pair<ElementPointer, LocalVector> const U;
+  std::pair<ElementPointer, LocalVector> const Z;
+
+public:
+  SpecialWriter(std::string filename, GridView const &gridView)
+      : writer_(filename, std::fstream::out),
+        hsearch_(gridView.grid(), gridView),
+        G(globalToLocal(MyGeometry::G)),
+        H(globalToLocal(MyGeometry::H)),
+        J(globalToLocal(MyGeometry::J)),
+        I(globalToLocal(MyGeometry::I)),
+        U(globalToLocal(MyGeometry::U)),
+        Z(globalToLocal(MyGeometry::Z)) {
+    writer_ << "Gh Hh Jh Ih Uv Uh Zv Zh" << std::endl;
+  }
+
+  void write(VirtualGridFunction<typename GridView::Grid, LocalVector> const &
+                 specialField) {
+    LocalVector value;
+
+    specialField.evaluateLocal(*G.first, G.second, value);
+    writeHorizontal(value);
+
+    specialField.evaluateLocal(*H.first, H.second, value);
+    writeHorizontal(value);
+
+    specialField.evaluateLocal(*J.first, J.second, value);
+    writeHorizontal(value);
+
+    specialField.evaluateLocal(*I.first, I.second, value);
+    writeHorizontal(value);
+
+    specialField.evaluateLocal(*U.first, U.second, value);
+    writeVertical(value);
+    writeHorizontal(value);
+
+    specialField.evaluateLocal(*Z.first, Z.second, value);
+    writeVertical(value);
+    writeHorizontal(value);
+
+    writer_ << std::endl;
+  }
+};
+
+#endif
diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
index b1c6178d04c75e3038beb047f6bf392edcab35f3..e2e7ecb162707868cbe85c0e5316b46d92128834 100644
--- a/src/sand-wedge.cc
+++ b/src/sand-wedge.cc
@@ -80,6 +80,7 @@
 #include "sand-wedge-data/myglobalfrictiondata.hh"
 #include "sand-wedge-data/mygrid.cc" // FIXME
 #include "sand-wedge-data/mygrid.hh"
+#include "sand-wedge-data/special_writer.hh"
 #include "solverfactory.hh"
 #include "state.hh"
 #include "timestepping.hh"
@@ -94,78 +95,6 @@ void initPython() {
   Python::run("sys.path.append('" datadir "')");
 }
 
-template <class GridView> class SpecialWriter {
-  using LocalVector = Dune::FieldVector<double, dims>;
-  using Element = typename GridView::Grid::template Codim<0>::Entity;
-  using ElementPointer =
-      typename GridView::Grid::template Codim<0>::EntityPointer;
-
-  void writeHorizontal(LocalVector const &v) {
-    writer_ << MyGeometry::horizontalProjection(v) << " ";
-  }
-  void writeVertical(LocalVector const &v) {
-    writer_ << MyGeometry::verticalProjection(v) << " ";
-  }
-
-  std::pair<ElementPointer, LocalVector> globalToLocal(LocalVector const &x)
-      const {
-    auto const element = hsearch_.findEntity(x);
-    return std::pair<ElementPointer, LocalVector>(element,
-                                                  element->geometry().local(x));
-  }
-
-  std::fstream writer_;
-
-  Dune::HierarchicSearch<typename GridView::Grid, GridView> const hsearch_;
-
-  std::pair<ElementPointer, LocalVector> const G;
-  std::pair<ElementPointer, LocalVector> const H;
-  std::pair<ElementPointer, LocalVector> const J;
-  std::pair<ElementPointer, LocalVector> const I;
-  std::pair<ElementPointer, LocalVector> const U;
-  std::pair<ElementPointer, LocalVector> const Z;
-
-public:
-  SpecialWriter(std::string filename, GridView const &gridView)
-      : writer_(filename, std::fstream::out),
-        hsearch_(gridView.grid(), gridView),
-        G(globalToLocal(MyGeometry::G)),
-        H(globalToLocal(MyGeometry::H)),
-        J(globalToLocal(MyGeometry::J)),
-        I(globalToLocal(MyGeometry::I)),
-        U(globalToLocal(MyGeometry::U)),
-        Z(globalToLocal(MyGeometry::Z)) {
-    writer_ << "Gh Hh Jh Ih Uv Uh Zv Zh" << std::endl;
-  }
-
-  void write(VirtualGridFunction<typename GridView::Grid, LocalVector> const &
-                 specialField) {
-    LocalVector value;
-
-    specialField.evaluateLocal(*G.first, G.second, value);
-    writeHorizontal(value);
-
-    specialField.evaluateLocal(*H.first, H.second, value);
-    writeHorizontal(value);
-
-    specialField.evaluateLocal(*J.first, J.second, value);
-    writeHorizontal(value);
-
-    specialField.evaluateLocal(*I.first, I.second, value);
-    writeHorizontal(value);
-
-    specialField.evaluateLocal(*U.first, U.second, value);
-    writeVertical(value);
-    writeHorizontal(value);
-
-    specialField.evaluateLocal(*Z.first, Z.second, value);
-    writeVertical(value);
-    writeHorizontal(value);
-
-    writer_ << std::endl;
-  }
-};
-
 int main(int argc, char *argv[]) {
   try {
     Dune::ParameterTree parset;
@@ -387,10 +316,10 @@ int main(int argc, char *argv[]) {
       vtkWriter.write(0, u_initial, v_initial, alpha_initial, stress);
     }
 
-    SpecialWriter<GridView> specialVelocityWriter("specialVelocities",
-                                                  leafView);
-    SpecialWriter<GridView> specialDisplacementWriter("specialDisplacements",
-                                                      leafView);
+    SpecialWriter<GridView, dims> specialVelocityWriter("specialVelocities",
+                                                        leafView);
+    SpecialWriter<GridView, dims> specialDisplacementWriter(
+        "specialDisplacements", leafView);
 
     // Set up TNNMG solver
     using NonlinearFactory =