diff --git a/dune/tectonic/globalruinanonlinearity.hh b/dune/tectonic/globalruinanonlinearity.hh
index 4c15db46287bc18288a73765c0a9c7f26d594169..4397aa341d59caa0be6ff9366a2a51d56eb5a39c 100644
--- a/dune/tectonic/globalruinanonlinearity.hh
+++ b/dune/tectonic/globalruinanonlinearity.hh
@@ -30,7 +30,13 @@ class GlobalRuinaNonlinearity
         eta(eta),
         normalStress(normalStress),
         trivialNonlinearity(
-            new LocalNonlinearity<dim>(make_shared<TrivialFunction const>())) {}
+            new LocalNonlinearity<dim>(make_shared<TrivialFunction const>())),
+        restrictions(nodalIntegrals->size()) // TODO: can we get the size from
+                                             // another place?
+  {
+    for (auto &x : restrictions)
+      x = shared_ptr<LocalNonlinearity<dim> const>();
+  }
 
   /*
     Return a restriction of the outer function to the i'th node.
@@ -39,11 +45,14 @@ class GlobalRuinaNonlinearity
     if ((*nodalIntegrals)[i][0] == 0)
       return trivialNonlinearity;
 
-    // TODO: cache functions
+    if (restrictions[i] != NULL)
+      return restrictions[i];
+
     auto const func = make_shared<RuinaFunction const>(
         (*nodalIntegrals)[i][0], (*a)[i][0], (*mu)[i][0], (*eta)[i][0],
         (*normalStress)[i][0]);
-    return make_shared<LocalNonlinearity<dim> const>(func);
+    restrictions[i] = make_shared<LocalNonlinearity<dim> const>(func);
+    return restrictions[i];
   }
 
 private:
@@ -54,6 +63,8 @@ class GlobalRuinaNonlinearity
   shared_ptr<BlockVector<FieldVector<double, 1>> const> mu;
   shared_ptr<BlockVector<FieldVector<double, 1>> const> eta;
   shared_ptr<BlockVector<FieldVector<double, 1>> const> normalStress;
+
+  std::vector<shared_ptr<LocalNonlinearity<dim> const>> mutable restrictions;
 };
 }
 #endif