[Extend] Handle descent at zero properly

dune/tectonic/ellipticenergy.hh

@@ -6,6 +6,7 @@
 #include <dune/common/stdstreams.hh>
 
 #include <dune/fufem/arithmetic.hh>
+#include <dune/solvers/common/interval.hh>
 
 #include "localfriction.hh"
 
@@ -30,14 +31,33 @@ template <size_t dim> class EllipticEnergy {
   std::shared_ptr<Nonlinearity const> const phi;
   typename Dune::BitSetVector<dim>::const_reference const ignore;
 
-  void gradient(LocalVector const &x, LocalVector &y) const {
-    A.mv(x, y);
-    y -= b;
-    phi->addGradient(x, y);
-
-    for (size_t i = 0; i < dim; ++i)
-      if (ignore[i])
-        y[i] = 0;
+  void descentDirection(LocalVector const &x, LocalVector &y) const {
+    if (x.two_norm() >= 0.0) {
+      A.mv(x, y);
+      y -= b;
+      phi->addGradient(x, y);
+
+      for (size_t i = 0; i < dim; ++i)
+        if (ignore[i])
+          y[i] = 0;
+      y *= -1;
+    } else {
+      A.mv(x, y);
+      y -= b;
+
+      for (size_t i = 0; i < dim; ++i)
+        if (ignore[i])
+          y[i] = 0;
+      y *= -1;
+
+      // The contribution of phi to the directional derivative is the
+      // same for all directions here! Choose the one that is best of
+      // the smooth part and check if we have overall decline
+      Dune::Solvers::Interval<double> D;
+      phi->directionalSubDiff(x, y, D);
+      if (D[1] - y.two_norm2() >= 0.0)
+        y = 0.0;
+    }
   }
 };
 #endif

dune/tectonic/minimisation.hh

@@ -39,12 +39,12 @@ void minimise(Functional const &J, typename Functional::LocalVector &x,
 
   for (size_t step = 0; step < steps; ++step) {
     LocalVector v;
-    J.gradient(x, v);
+    J.descentDirection(x, v);
 
     double const vnorm = v.two_norm();
     if (vnorm <= 0.0)
       return;
-    v /= -vnorm; // normalise for numerical stability; note the minus
+    v /= vnorm;
 
     double const alpha = lineSearch(J, x, v, bisection);
     Arithmetic::addProduct(x, alpha, v);