diff --git a/dune/tectonic/minimisation.hh b/dune/tectonic/minimisation.hh index 79f1d4686f07ed15b9cf84a2eff18bfe8798f902..5407fc3f3e0211bb0c2373c58d6b93772da2fa22 100644 --- a/dune/tectonic/minimisation.hh +++ b/dune/tectonic/minimisation.hh @@ -117,63 +117,5 @@ 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 = J.get_kinks(); - - 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(0); - x2.axpy(-1, x1); - 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) - break; - - Jx_old = Jx; - x_old = x; - } - startingPoint = x_old; -} } #endif diff --git a/dune/tectonic/minimisation2.hh b/dune/tectonic/minimisation2.hh new file mode 100644 index 0000000000000000000000000000000000000000..c78062c9fd4e548bbb1eeb43391adb8f144b24d8 --- /dev/null +++ b/dune/tectonic/minimisation2.hh @@ -0,0 +1,70 @@ +#ifndef MINIMISATION2_HH +#define MINIMISATION2_HH + +#include <dune/common/fmatrix.hh> +#include <dune/common/fvector.hh> + +#include <dune/tnnmg/problem-classes/bisection.hh> + +#include "minimisation.hh" + +namespace Dune { +// 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 = J.get_kinks(); + + 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(0); + x2.axpy(-1, x1); + 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) + break; + + Jx_old = Jx; + x_old = x; + } + startingPoint = x_old; +} +} +#endif diff --git a/src/test-gradient-method-helper.hh b/src/test-gradient-method-helper.hh index a9479ebb0dbbec8af0b6068bba1194a767cad2d0..2bb22ad9e587a4c589fbecfeed93b153e23bec8c 100644 --- a/src/test-gradient-method-helper.hh +++ b/src/test-gradient-method-helper.hh @@ -7,6 +7,7 @@ #include <dune/tectonic/ellipticenergy.hh> #include <dune/tectonic/minimisation.hh> +#include <dune/tectonic/minimisation2.hh> template <int dim> double functionTester(Dune::EllipticEnergy<dim> J,