diff --git a/dune/tectonic/globalnonlinearity.hh b/dune/tectonic/globalnonlinearity.hh
index 73042d741c95c8ac8f30c4bca6e31d977f28ea50..a580e81f87ad51271a3c26e8662ce875e14832ae 100644
--- a/dune/tectonic/globalnonlinearity.hh
+++ b/dune/tectonic/globalnonlinearity.hh
@@ -26,6 +26,13 @@ class GlobalNonlinearity {
       res->addHessian(v[i], hessian[i][i]);
     }
   }
+
+  virtual void addGradient(const VectorType& v, VectorType& gradient) const {
+    for (size_t i = 0; i < v.size(); ++i) {
+      auto res = restriction(i);
+      res->addGradient(v[i], gradient[i]);
+    }
+  }
 };
 }
 #endif
diff --git a/dune/tectonic/localnonlinearity.hh b/dune/tectonic/localnonlinearity.hh
index 3dc5e53b149d2483c169e1fc9cea27ea04f5de90..285f309e95d9d4ebd05efa9435d7cf8937c756fc 100644
--- a/dune/tectonic/localnonlinearity.hh
+++ b/dune/tectonic/localnonlinearity.hh
@@ -78,6 +78,12 @@ template <int dimension> class LocalNonlinearity {
       A[k][k] += (h2 - h1ox) * x[k] * x[k] / normX2 + h1ox;
   }
 
+  void addGradient(VectorType const &x, VectorType &y) const {
+    VectorType tmp;
+    upperGradient(x, tmp); // TODO
+    y += tmp;
+  }
+
   void upperGradient(VectorType const &x, VectorType &ret) const {
     ret = x;
     ret *= func_->rightDifferential(x.two_norm()) / x.two_norm();