Use axpy

dune/tectonic/localnonlinearity.hh

@@ -88,15 +88,13 @@ template <int dimension> class LocalNonlinearity {
       A[k][k] += (h2 - h1ox) * x[k] * x[k] / normX2 + h1ox;
   }
 
-  // TODO: do not evaluate at zero
   void addGradient(VectorType const &x, VectorType &y) const {
-    if (x.two_norm() <
-        1e-8 or std::isinf(func_->regularity(x.two_norm()))) // TODO
+    double const xnorm = x.two_norm();
+    if (xnorm < 1e-8 or std::isinf(func_->regularity(xnorm))) // TODO
       return;
 
-    VectorType tmp;
-    upperGradient(x, tmp); // upper and lower gradient coincide in this case
-    y += tmp;
+    // left and right differential coincide in this case
+    y.axpy(func_->rightDifferential(xnorm) / xnorm, x);
   }
 
   void upperGradient(VectorType const &x, VectorType &ret) const {