The descent direction is just "-gradient"

dune/tectonic/ellipticenergy.hh

@@ -29,46 +29,6 @@ template <int dim> class EllipticEnergy {
     return y * v + (*phi)(v); // <1/2 Av - b,v> + H(|v|)
   }
 
-  void descentAtZero(SmallVector &ret) const {
-    SmallVector const zero(0);
-    // If there is a direction of descent, this is it
-    smoothGradient(zero, ret);
-    ret *= -1;
-  }
-
-  bool descentDirection(SmallVector const &x, SmallVector &ret) const {
-    // Check the squared norm rather than each component because
-    // complementaryProjection() divides by it
-    if (x.two_norm2() == 0.0) {
-      descentAtZero(ret);
-      return true;
-    }
-
-    SmallVector pg;
-    gradient(x, pg);
-    SmallVector mg;
-    gradient(x, mg);
-    double const pgx = pg * x;
-    double const mgx = mg * x;
-    if (pgx >= 0 && mgx >= 0) {
-      ret = pg;
-      ret *= -1;
-      return true;
-    } else if (pgx <= 0 && mgx <= 0) {
-      ret = mg;
-      ret *= -1;
-      return true;
-    } else {
-      // Includes the case that pg points in direction x and mg
-      // points in direction -x. The projection will then be zero.
-      SmallVector d;
-      smoothGradient(x, d);
-      complementaryProjection(d, x, ret);
-      ret *= -1;
-      return false;
-    }
-  }
-
   SmallMatrix const &A;
   SmallVector const &b;
   shared_ptr<NonlinearityType const> const phi;
@@ -76,32 +36,14 @@ template <int dim> class EllipticEnergy {
                     // to dim-1; the special value dim means that no
                     // dimension should be ignored
 
-private:
-  // Gradient of the smooth part
-  void smoothGradient(SmallVector const &x, SmallVector &y) const {
-    A.mv(x, y); // y = Av
-    y -= b;     // y = Av - b
-    if (ignore != dim)
-      y[ignore] = 0;
-  }
-
   void gradient(SmallVector const &x, SmallVector &y) const {
-    phi->gradient(x, y);
-    SmallVector z;
-    smoothGradient(x, z);
-    y += z;
+    A.mv(x, y);
+    y -= b;
+    phi->addGradient(x, y);
+
     if (ignore != dim)
       y[ignore] = 0;
   }
-
-  // y = (id - P)(d) where P is the projection onto the line t*x
-  void complementaryProjection(SmallVector const &d, SmallVector const &x,
-                               SmallVector &y) const {
-    double const dx = d * x;
-    double const xx = x.two_norm2();
-    y = d;
-    y.axpy(-dx / xx, x);
-  }
 };
 }
 #endif

dune/tectonic/localfriction.hh

@@ -114,12 +114,8 @@ template <int dimension> class LocalFriction {
     if (std::isinf(func->regularity(xnorm)))
       return;
 
-    y.axpy(func->differential(xnorm) / xnorm, x);
-  }
-
-  void gradient(VectorType const &x, VectorType &ret) const {
-    ret = 0;
-    addGradient(x, ret);
+    if (xnorm >= smp)
+      y.axpy(func->differential(xnorm) / xnorm, x);
   }
 
   void directionalDomain(VectorType const &, VectorType const &,

dune/tectonic/minimisation.hh

@@ -76,7 +76,8 @@ void minimise(Functional const &J, typename Functional::SmallVector &x,
 
   for (size_t step = 0; step < steps; ++step) {
     SmallVector descDir;
-    J.descentDirection(x, descDir);
+    J.gradient(x, descDir);
+    descDir *= -1;
 
     if (descDir.two_norm() < 1e-14) // TODO: Make controllable
       return;