From aa2b5ac892ceeea2d65af540c151e80c3aabe5e1 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Fri, 31 May 2013 17:58:46 +0200
Subject: [PATCH] The descent direction is just "-gradient"
---
dune/tectonic/ellipticenergy.hh | 66 ++-------------------------------
dune/tectonic/localfriction.hh | 8 +---
dune/tectonic/minimisation.hh | 3 +-
3 files changed, 8 insertions(+), 69 deletions(-)
diff --git a/dune/tectonic/ellipticenergy.hh b/dune/tectonic/ellipticenergy.hh
index 36d5ec3c..d078f7c7 100644
--- a/dune/tectonic/ellipticenergy.hh
+++ b/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
diff --git a/dune/tectonic/localfriction.hh b/dune/tectonic/localfriction.hh
index 02b5b2e2..5edabfca 100644
--- a/dune/tectonic/localfriction.hh
+++ b/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 &,
diff --git a/dune/tectonic/minimisation.hh b/dune/tectonic/minimisation.hh
index 9bae48be..ad079b02 100644
--- a/dune/tectonic/minimisation.hh
+++ b/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;
--
GitLab