From ef039d44adcc8b5add814f9c4a99273b97818a27 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Fri, 31 Jan 2014 18:16:01 +0100
Subject: [PATCH] [Extend] Handle descent at zero properly
---
dune/tectonic/ellipticenergy.hh | 36 +++++++++++++++++++++++++--------
dune/tectonic/minimisation.hh | 4 ++--
2 files changed, 30 insertions(+), 10 deletions(-)
diff --git a/dune/tectonic/ellipticenergy.hh b/dune/tectonic/ellipticenergy.hh
index c00776ce..75b9985f 100644
--- a/dune/tectonic/ellipticenergy.hh
+++ b/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
diff --git a/dune/tectonic/minimisation.hh b/dune/tectonic/minimisation.hh
index 0e27a654..a7ad57bc 100644
--- a/dune/tectonic/minimisation.hh
+++ b/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);
--
GitLab