Make sure circles do not have more than two minima

src/Makefile.am

 check_PROGRAMS = \
+	test-circle \
 	test-gradient-horrible \
 	test-gradient-horrible-logarithmic \
 	test-gradient-identity \
@@ -12,6 +13,7 @@ check_PROGRAMS = \
 	test-gradient-sample2 \
 	test-gradient-trivial
 
+test_circle_SOURCES                        = test-circle.cc
 test_gradient_horrible_SOURCES             = test-gradient-horrible.cc
 test_gradient_horrible_logarithmic_SOURCES = test-gradient-horrible-logarithmic.cc
 test_gradient_identity_SOURCES             = test-gradient-identity.cc

src/test-circle.cc

+/* Assures that a circle never has more than two minima and that the
+   bisection takes us to one after no more than <iterations> steps */
+
+#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 A;
+  A[0][0] = 4;
+  A[0][1] = A[1][0] = 1.5;
+  A[1][1] = 3;
+  SmallVector b;
+  b[0] = 1;
+  b[1] = 2;
+
+  auto const f = Dune::make_shared<Dune::SampleFunction<2> const>();
+  auto const phi = Dune::make_shared<Functional::NonlinearityType const>(f);
+  Functional J(A, b, phi);
+
+  Bisection const bisection(0.0, 1.0, 1e-12, false, 0);
+
+  double scale = 10;
+
+  std::vector<SmallVector> minima(4);
+  std::vector<double> radii = { M_PI / 4.0,     3 * M_PI / 4.0,
+                                5 * M_PI / 4.0, 7 * M_PI / 4.0 };
+
+  boost::format const formatter("J([%+e]) = %g");
+
+  for (size_t i = 0; i < radii.size(); ++i) {
+    SmallVector x;
+    x[0] = scale * std::sin(radii[i]);
+    x[1] = scale * std::cos(radii[i]);
+    SmallVector descDir;
+    descDir[0] = x[1];
+    descDir[1] = -x[0];
+    tangentialMinimisation(J, x, descDir, bisection);
+    minima[i] = x;
+    std::cout << boost::format(formatter) % x % J(x) << std::endl;
+  }
+
+  double const intervals = 10000;
+  double const iterations = 2;
+  for (size_t i = 0; i < intervals; ++i) {
+    double alpha = i / (double)intervals * 2 * M_PI;
+    SmallVector x;
+    x[0] = scale * std::sin(alpha);
+    x[1] = scale * std::cos(alpha);
+    SmallVector descDir;
+
+    for (int i = 0; i < iterations; ++i) {
+      descDir[0] = x[1];
+      descDir[1] = -x[0];
+      tangentialMinimisation(J, x, descDir, bisection);
+    }
+
+    bool minimum_hit(false);
+    for (auto const &minimum : minima) {
+      if (two_distance<dim>(x, minimum) < 1e-12) {
+        minimum_hit = true;
+        break;
+      }
+    }
+    if (!minimum_hit)
+      std::cout << std::endl << boost::format(formatter) % x % J(x)
+                << std::endl;
+    assert(minimum_hit);
+  }
+}

src/test_circle_10.m

+if exist('graphics_toolkit','file')
+  graphics_toolkit('fltk')
+end
+
+
+A   = [4 1.5; 1.5 3]; b=[1; 2];
+f   = @(x) x + (x > 1) .* (x - 1);
+phi = @(x) f(norm(x,2));
+J   = @(x) .5*dot(A*x,x) - dot(b,x) + phi(x);
+
+scale = 10.0;
+
+t = -5:0.1:5;
+x = scale*cos(t);
+y = scale*sin(t);
+z = arrayfun(@(vec1, vec2) J([vec1; vec2]), x, y);
+
+plot3(x,y,z);
+
+title 'scale = 10'
+
+# Minima taken from test-circle.cc
+
+line([-5.3444 6.36022],
+     [8.45206 -7.71671],
+     [J([-5.3444; 8.45206]) J([6.36022; -7.71671])])