diff --git a/dune/tectonic/ellipticenergy.hh b/dune/tectonic/ellipticenergy.hh
index d64b8ab123d864462ac1c2c8ae318334e27c8922..97e44277e435b3f5463494da1eb0b44597615450 100644
--- a/dune/tectonic/ellipticenergy.hh
+++ b/dune/tectonic/ellipticenergy.hh
@@ -52,6 +52,24 @@ template <int dim> class EllipticEnergy {
}
// returns false if the direction is tangential
+ void descentDirectionNew(SmallVector const &x, SmallVector &ret) const {
+ SmallVector pg;
+ upperGradient(x, pg);
+ SmallVector mg;
+ lowerGradient(x, mg);
+ double const pgx = pg * x;
+ double const mgx = mg * x;
+ if (pgx >= 0 && mgx >= 0) {
+ ret = pg;
+ ret *= -1;
+ } else if (pgx <= 0 && mgx <= 0) {
+ ret = mg;
+ ret *= -1;
+ } else {
+ assert(false);
+ }
+ }
+
bool descentDirection(SmallVector const &x, SmallVector &ret) const {
// Check the squared norm rather than each component because
// complementaryProjection() divides by it
@@ -234,5 +252,65 @@ void minimise(Functional const &J, typename Functional::SmallVector &x,
}
}
}
+
+// NOTE: only works for the 2D case
+template <class Functional>
+void minimise2(Functional const &J, typename Functional::SmallVector &x,
+ size_t steps, Bisection const &bisection) {
+ typedef typename Functional::SmallVector SmallVector;
+ minimisationInitialiser(J, bisection, x);
+
+ { // Make a single step where we already know we're not differentiable
+ SmallVector descDir;
+ if (x == SmallVector(0))
+ J.descentAtZero(descDir);
+ else
+ J.descentDirectionNew(x, descDir);
+ descentMinimisation(J, x, descDir, bisection);
+ }
+
+ // From here on, we're in a C1 region and stay there.
+ for (size_t i = 1; i < steps; ++i) {
+ SmallVector descDir;
+ J.descentDirectionNew(x, descDir);
+ descentMinimisation(J, x, descDir, bisection);
+ }
+}
+
+template <class Functional>
+void minimisationInitialiser(Functional const &J, Bisection const &bisection,
+ typename Functional::SmallVector &startingPoint) {
+ typedef typename Functional::SmallVector SmallVector;
+
+ std::vector<double> const kinks = { 5, 10,
+ 15 }; // FIXME: We happen to know these
+
+ SmallVector x_old(0);
+ double Jx_old = J(x_old);
+
+ for (auto &kink : kinks) {
+ SmallVector x1 = { kink, 0 }; // Random vector that has the right norm
+ SmallVector const descDir1 = { x1[1], -x1[0] };
+ tangentialMinimisation(J, x1, descDir1, bisection);
+ double const Jx1 = J(x1);
+
+ SmallVector x2 = { -x1[0], -x1[1] };
+ SmallVector const descDir2 = { x2[1], -x2[0] };
+ tangentialMinimisation(J, x2, descDir2, bisection);
+ double const Jx2 = J(x2);
+
+ double const Jx = (Jx2 < Jx1 ? Jx2 : Jx1);
+ SmallVector const x = (Jx2 < Jx1 ? x2 : x1);
+
+ if (Jx >= Jx_old) {
+ startingPoint = x_old;
+ return;
+ }
+
+ Jx_old = Jx;
+ x_old = x;
+ }
+ startingPoint = x_old;
+}
}
#endif
diff --git a/src/Makefile.am b/src/Makefile.am
index 974967c314653e3b70eda30d7f46284513e3db2a..2467aafbdfe3640ba86b5db998b0fa8ee7df0d13 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -13,7 +13,8 @@ check_PROGRAMS = \
test-gradient-sample-steep2 \
test-gradient-sample-verysteep \
test-gradient-sample2 \
- test-gradient-trivial
+ test-gradient-trivial \
+ test-minimise2
test_circle_1_SOURCES = test-circle.cc
test_circle_1_CPPFLAGS = $(AM_CPPFLAGS) -DDUNE_TECTONIC_TEST_CIRCLE_SCALE=1
@@ -32,6 +33,7 @@ test_gradient_sample_steep2_SOURCES = test-gradient-sample-steep2.cc
test_gradient_sample_verysteep_SOURCES = test-gradient-sample-verysteep.cc
test_gradient_sample2_SOURCES = test-gradient-sample2.cc
test_gradient_trivial_SOURCES = test-gradient-trivial.cc
+test_minimise2_SOURCES = test-minimise2.cc
bin_PROGRAMS = \
one-body-sample-2D \
diff --git a/src/test-minimise2.cc b/src/test-minimise2.cc
new file mode 100644
index 0000000000000000000000000000000000000000..1646bc0ce53733cace53a0cad91256996f5788ba
--- /dev/null
+++ b/src/test-minimise2.cc
@@ -0,0 +1,58 @@
+/* Check that minimise2 always produces better results than minimise
+ when the latter starts from 0 */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <cassert>
+
+#include <boost/format.hpp>
+
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/tectonic/ellipticenergy.hh>
+
+#include "test-gradient-method-nicefunction.hh"
+#include "test-gradient-method-helper.hh"
+
+int main() {
+ int const dim = 2;
+ typedef Dune::EllipticEnergy<dim> Functional;
+ typedef Functional::SmallMatrix SmallMatrix;
+ typedef Functional::SmallVector SmallVector;
+
+ SmallMatrix const A = { { 4, 1.5 }, { 1.5, 3 } };
+
+ std::vector<SmallVector> const bs = {
+ { 1, 2 }, { -8, -14 }, { -16, -28 }, { -24, -42 }, { -32, -56 },
+ };
+
+ auto const f = Dune::make_shared<Dune::ThreeKinkFunction const>();
+ auto const phi = Dune::make_shared<Functional::NonlinearityType const>(f);
+
+ std::vector<Functional> const Js = { Functional(A, bs[0], phi),
+ Functional(A, bs[1], phi),
+ Functional(A, bs[2], phi),
+ Functional(A, bs[3], phi),
+ Functional(A, bs[4], phi) };
+
+ Bisection const bisection(0.0, 1.0, 1e-12, false, 0);
+
+ size_t const runs = 5;
+
+ SmallVector x;
+
+ for (auto const &J : Js) {
+ x = { 0, 0 };
+ Dune::minimise(J, x, runs, bisection);
+ double const xmin = J(x);
+
+ x = { 0, 0 };
+ Dune::minimise2(J, x, runs, bisection);
+ double const xmin2 = J(x);
+ std::cout << xmin << std::endl;
+ std::cout << xmin2 << std::endl << std::endl;
+ assert(xmin2 <= xmin);
+ }
+}