diff --git a/dune/tectonic/ellipticenergy.hh b/dune/tectonic/ellipticenergy.hh index 89b403e22259520b433f635a6fbcc1dbd3e29b4a..39aa3f0faa5b72d5206eb74472ee21b5dc61c351 100644 --- a/dune/tectonic/ellipticenergy.hh +++ b/dune/tectonic/ellipticenergy.hh @@ -31,27 +31,32 @@ 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 + SmallVector d; + smoothGradient(zero, d); + d *= -1; + + Interval<double> D; + phi->directionalSubDiff(zero, d, D); + double const nonlinearDecline = D[1]; + double const smoothDecline = -(d * d); + double const combinedDecline = smoothDecline + nonlinearDecline; + + if (combinedDecline < 0) { + ret = d; + } else { + ret = 0; + } + } + // returns false if the direction is tangential 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) { - // If there is a direction of descent, this is it - SmallVector d; - smoothGradient(x, d); - d *= -1; - - Interval<double> D; - phi->directionalSubDiff(x, d, D); - double const nonlinearDecline = D[1]; - double const smoothDecline = -(d * d); - double const combinedDecline = smoothDecline + nonlinearDecline; - - if (combinedDecline < 0) { - ret = d; - } else { - ret = 0; - } + descentAtZero(ret); return true; }