diff --git a/dune/tectonic/globallaursennonlinearity.hh b/dune/tectonic/globallaursennonlinearity.hh
new file mode 100644
index 0000000000000000000000000000000000000000..539abedb2a50fb50c1563204ab665dde5fb6ee92
--- /dev/null
+++ b/dune/tectonic/globallaursennonlinearity.hh
@@ -0,0 +1,52 @@
+/* -*- mode:c++; mode:semantic -*- */
+
+#ifndef DUNE_TECTONIC_GLOBAL_LAURSEN_NONLINEARITY_HH
+#define DUNE_TECTONIC_GLOBAL_LAURSEN_NONLINEARITY_HH
+
+#include <dune/common/fvector.hh>
+
+#include "globalnonlinearity.hh"
+
+namespace Dune {
+template <int dim, class OuterFunctionType>
+class GlobalLaursenNonlinearity
+    : public Dune::GlobalNonlinearity<dim, OuterFunctionType> {
+public:
+  GlobalLaursenNonlinearity(
+      std::vector<double> const &coefficientOfFriction,
+      std::vector<double> const &normalStress,
+      std::vector<Dune::FieldVector<double, 1>> const &nodalIntegrals)
+      : coefficientOfFriction(coefficientOfFriction),
+        normalStress(normalStress),
+        nodalIntegrals(nodalIntegrals) {}
+
+  /*
+    Return a restriction of the outer function to the i'th node. If
+    mu and sigma_n denote the coefficient of friction and the normal
+    stress, respectively, at the i'th node, this function is given
+    by
+
+    sigma_n [(1/eta \bar Gamma)* + mu id]
+
+    TODO: We chose Gamma = id, so that (\bar Gamma)* = \Gamma^{-1}
+    = id^{-1} = id. The factor 1/eta cancels in this special case, leaving us
+    with
+
+    sigma_n [id + mu id] = sigma_n (1 + mu) id
+  */
+  virtual void restriction(int i, OuterFunctionType &f) const {
+    double coefficient = nodalIntegrals[i][0];
+    coefficient *= normalStress[i];
+    coefficient *= 1 + coefficientOfFriction[i];
+    f = OuterFunctionType(coefficient);
+  }
+
+private:
+  // TODO: If we're clever, we only store one vector with the precomputed
+  // results
+  std::vector<double> coefficientOfFriction;
+  std::vector<double> normalStress;
+  std::vector<Dune::FieldVector<double, 1>> nodalIntegrals;
+};
+}
+#endif
diff --git a/dune/tectonic/globalnonlinearity.hh b/dune/tectonic/globalnonlinearity.hh
index a33e4a4eaf1da9488f371d969bfbb3d49617d5e2..a72c2bb29aa3d85ddd8a5d63dc2ae34ea907a78d 100644
--- a/dune/tectonic/globalnonlinearity.hh
+++ b/dune/tectonic/globalnonlinearity.hh
@@ -8,41 +8,10 @@
 namespace Dune {
 template <int dim, class OuterFunctionType> class GlobalNonlinearity {
 public:
-  GlobalNonlinearity(
-      std::vector<double> const &coefficientOfFriction,
-      std::vector<double> const &normalStress,
-      std::vector<Dune::FieldVector<double, 1>> const &nodalIntegrals)
-      : coefficientOfFriction(coefficientOfFriction),
-        normalStress(normalStress),
-        nodalIntegrals(nodalIntegrals) {}
-
   /*
     Return a restriction of the outer function to the i'th node. If
-    mu and sigma_n denote the coefficient of friction and the normal
-    stress, respectively, at the i'th node, this function is given
-    by
-
-    sigma_n [(1/eta \bar Gamma)* + mu id]
-
-    TODO: We chose Gamma = id, so that (\bar Gamma)* = \Gamma^{-1}
-    = id^{-1} = id. The factor 1/eta cancels in this special case, leaving us
-    with
-
-    sigma_n [id + mu id] = sigma_n (1 + mu) id
   */
-  void restriction(int i, OuterFunctionType &f) const {
-    double coefficient = nodalIntegrals[i][0];
-    coefficient *= normalStress[i];
-    coefficient *= 1 + coefficientOfFriction[i];
-    f = OuterFunctionType(coefficient);
-  }
-
-private:
-  // TODO: If we're clever, we only store one vector with the precomputed
-  // results
-  std::vector<double> coefficientOfFriction;
-  std::vector<double> normalStress;
-  std::vector<Dune::FieldVector<double, 1>> nodalIntegrals;
+  virtual void restriction(int i, OuterFunctionType &f) const = 0;
 };
 }
 #endif
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 826e7dee5960d9c13023a2bab05c6e081b5595c4..3b0b6146ff8f763aed1cdfe906e28d28977d92e7 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -41,7 +41,7 @@
 #include <dune/solvers/solvers/loopsolver.hh>
 #include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
 
-#include <dune/tectonic/globalnonlinearity.hh>
+#include <dune/tectonic/globallaursennonlinearity.hh>
 #include <dune/tectonic/myconvexproblem.hh>
 #include <dune/tectonic/myblockproblem.hh>
 
@@ -226,8 +226,8 @@ int main(int argc, char *argv[]) {
       std::fill(coefficientOfFriction.begin(), coefficientOfFriction.end(),
                 parset.get<double>("boundary.mu"));
 
-      Dune::GlobalNonlinearity<dim, Dune::LinearFunction> myGlobalNonlinearity(
-          coefficientOfFriction, normalStress, nodalIntegrals);
+      Dune::GlobalLaursenNonlinearity<dim, Dune::LinearFunction>
+      myGlobalNonlinearity(coefficientOfFriction, normalStress, nodalIntegrals);
       // }}}
 
       {