diff --git a/dune/tectonic/localnonlinearity.hh b/dune/tectonic/localnonlinearity.hh
index ee3d81645ea043cff35ac548a535a91122911927..bb1edb7ca3fbb3fccdf0018707b60ffdddd2e260 100644
--- a/dune/tectonic/localnonlinearity.hh
+++ b/dune/tectonic/localnonlinearity.hh
@@ -21,9 +21,7 @@ template <int dimension> class LocalNonlinearity {
   LocalNonlinearity(shared_ptr<NiceFunction const> func) : func(func) {}
 
   double operator()(VectorType const &x) const {
-    double ret;
-    func->evaluate(x.two_norm(), ret);
-    return ret;
+    return func->evaluate(x.two_norm());
   }
 
   double regularity(VectorType const &x) const {
diff --git a/dune/tectonic/nicefunction.hh b/dune/tectonic/nicefunction.hh
index 2c5e793be54cef75a7a4c4d2063021b5fb3698a4..39df709bb8b639b5fc74ce788a05ff60f06bfaa2 100644
--- a/dune/tectonic/nicefunction.hh
+++ b/dune/tectonic/nicefunction.hh
@@ -23,7 +23,7 @@ class NiceFunction {
   // Whether H(|.|) is smooth at zero
   bool virtual smoothesNorm() const { return false; }
 
-  void virtual evaluate(double const &x, double &y) const {
+  double virtual evaluate(double x) const {
     DUNE_THROW(NotImplemented, "evaluation not implemented");
   }
 };
@@ -83,7 +83,7 @@ class RuinaFunction : public NiceFunction {
 
 class TrivialFunction : public NiceFunction {
 public:
-  void virtual evaluate(double const &x, double &y) const { y = 0; }
+  double virtual evaluate(double) const { return 0; }
 
   double virtual leftDifferential(double) const { return 0; }
 
diff --git a/src/test-gradient-method-nicefunction.hh b/src/test-gradient-method-nicefunction.hh
index 7ef9e64236d924fcf92b031bf83faab0c4b5e675..bb80dc7b8dc6b53f182a087d328ee8c210ea7eea 100644
--- a/src/test-gradient-method-nicefunction.hh
+++ b/src/test-gradient-method-nicefunction.hh
@@ -19,7 +19,7 @@ class MyFunction : public NiceFunction {
 };
 
 class Parabola : public MyFunction {
-  void virtual evaluate(double const &x, double &y) const { y = x * x; }
+  double virtual evaluate(double x) const { return x * x; }
 
   double virtual leftDifferential(double s) const { return 2 * s; }
 
@@ -36,9 +36,7 @@ class LinearFunction : public MyFunction {
 public:
   LinearFunction(double a) : coefficient(a) {}
 
-  void virtual evaluate(double const &x, double &y) const {
-    y = coefficient * x;
-  }
+  double virtual evaluate(double x) const { return coefficient * x; }
 
   double virtual leftDifferential(double s) const { return coefficient; }
 
@@ -54,8 +52,8 @@ class LinearFunction : public MyFunction {
 
 template <int slope> class SampleFunction : public MyFunction {
 public:
-  void virtual evaluate(double const &x, double &y) const {
-    y = (x < 1) ? x : (slope * (x - 1) + 1);
+  double virtual evaluate(double x) const {
+    return (x < 1) ? x : (slope * (x - 1) + 1);
   }
 
   double virtual leftDifferential(double s) const {
@@ -69,7 +67,7 @@ template <int slope> class SampleFunction : public MyFunction {
 
 class SteepFunction : public MyFunction {
 public:
-  void virtual evaluate(double const &x, double &y) const { y = 100 * x; }
+  double virtual evaluate(double x) const { return 100 * x; }
 
   double virtual leftDifferential(double s) const { return 100; }
 
@@ -79,10 +77,10 @@ class SteepFunction : public MyFunction {
 // slope in [n-1,n] is n
 class HorribleFunction : public MyFunction {
 public:
-  void virtual evaluate(double const &x, double &y) const {
+  double virtual evaluate(double x) const {
     double const fl = floor(x);
     double const sum = fl * (fl + 1) / 2;
-    y = sum + (fl + 1) * (x - fl);
+    return sum + (fl + 1) * (x - fl);
   }
 
   double virtual leftDifferential(double x) const {
@@ -105,14 +103,14 @@ class HorribleFunction : public MyFunction {
 // slope in [n-1,n] is log(n+1)
 class HorribleFunctionLogarithmic : public MyFunction {
 public:
-  void virtual evaluate(double const &x, double &y) const {
-    y = 0;
+  double virtual evaluate(double x) const {
+    double y = 0;
     size_t const fl = floor(x);
     for (size_t i = 1; i <= fl;)
       y += std::log(
           ++i); // factorials grow to fast so we compute this incrementally
 
-    y += std::log(fl + 2) * (x - fl);
+    return y + std::log(fl + 2) * (x - fl);
   }
 
   double virtual leftDifferential(double x) const {