test-circle.cc

/* Assures that a circle never has more than two minima and that the
   bisection takes us to one in a single step */

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#ifndef DUNE_TECTONIC_TEST_CIRCLE_SCALE
#error DUNE_TECTONIC_TEST_CIRCLE_SCALE cannot be unset
#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 } };
  SmallVector const b = { 1, 2 };

  auto const f = Dune::make_shared<Dune::SampleFunction<2> const>();
  auto const phi = Dune::make_shared<Functional::NonlinearityType const>(f);
  Functional const J(A, b, phi);

  Bisection const bisection(0.0, 1.0, 1e-12, false, 0);

  double scale = DUNE_TECTONIC_TEST_CIRCLE_SCALE;

  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 = { scale * std::sin(radii[i]), scale * std::cos(radii[i]) };
    SmallVector descDir = { x[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;
  for (size_t i = 0; i < intervals; ++i) {
    double const alpha = i / (double)intervals * 2 * M_PI;
    SmallVector x = { scale * std::sin(alpha), scale * std::cos(alpha) };
    SmallVector descDir = { x[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);
  }
}