diff --git a/dune/tectonic/ellipticenergy.hh b/dune/tectonic/ellipticenergy.hh index 36d5ec3cc05459b9f4edc4466f8ccd310b4a9ecc..d078f7c7f0795f2e3528b851e1bb97b97ffd4139 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 02b5b2e2e726148d21fad7b152ff86e34dd16e0f..5edabfca11b22d25d913ef060f83a5021f56ac2f 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 9bae48be8e2f0189db270f91c27ce1f9077842ed..ad079b027c187c60e4567603527e6bea8c748e08 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;