diff --git a/configure.ac b/configure.ac
index 30bd4d5d1e623c84dc2ffa8e05d1ebd135155639..f1c6b26683e419d75a683c127f39e19aa0cc6edf 100644
--- a/configure.ac
+++ b/configure.ac
@@ -12,7 +12,6 @@ DUNE_CHECK_ALL
 AC_CONFIG_FILES([
   Makefile
   src/Makefile
-  src/tests/Makefile
   dune/Makefile
   dune/tectonic/Makefile
   m4/Makefile
diff --git a/dune/tectonic/nicefunction.hh b/dune/tectonic/frictionpotential.hh
similarity index 85%
rename from dune/tectonic/nicefunction.hh
rename to dune/tectonic/frictionpotential.hh
index 7943a2c71e454bf3a5671042e6ca2e8134dcd2fe..d01659d10ff8a84c26f23e48a44f8bc546ed3d49 100644
--- a/dune/tectonic/nicefunction.hh
+++ b/dune/tectonic/frictionpotential.hh
@@ -12,9 +12,9 @@
 #include "frictiondata.hh"
 
 namespace Dune {
-class NiceFunction {
+class FrictionPotentialWrapper {
 protected:
-  NiceFunction(std::vector<double> const &_kinks) : kinks(_kinks) {}
+  FrictionPotentialWrapper(std::vector<double> const &_kinks) : kinks(_kinks) {}
 
 private:
   std::vector<double> const kinks;
@@ -22,9 +22,9 @@ class NiceFunction {
 public:
   std::vector<double> const &get_kinks() const { return kinks; }
 
-  NiceFunction() : kinks() {}
+  FrictionPotentialWrapper() : kinks() {}
 
-  virtual ~NiceFunction() {}
+  virtual ~FrictionPotentialWrapper() {}
 
   double virtual leftDifferential(double s) const = 0;
   double virtual rightDifferential(double s) const = 0;
@@ -45,10 +45,10 @@ class NiceFunction {
 //        = 0                       otherwise
 // with V_m = V_0 exp(-K/a),
 // i.e.   K = -a log(V_m / V_0) = mu_0 + b log(V_0 / L) + b alpha
-class RuinaFunction : public NiceFunction {
+class FrictionPotential : public FrictionPotentialWrapper {
 public:
-  RuinaFunction(double coefficient, FrictionData const &fd, double state)
-      : NiceFunction(),
+  FrictionPotential(double coefficient, FrictionData const &fd, double state)
+      : FrictionPotentialWrapper(),
         coefficientProduct(coefficient * fd.a * fd.normalStress),
         // state is assumed to be logarithmic
         V_m(fd.V0 *
@@ -110,9 +110,9 @@ class RuinaFunction : public NiceFunction {
   double const V_m;
 };
 
-class TrivialFunction : public NiceFunction {
+class TrivialFunction : public FrictionPotentialWrapper {
 public:
-  TrivialFunction() : NiceFunction() {}
+  TrivialFunction() : FrictionPotentialWrapper() {}
 
   double virtual evaluate(double) const { return 0; }
 
diff --git a/dune/tectonic/globalnonlinearity.hh b/dune/tectonic/globalnonlinearity.hh
index 2c600b0b4d45d410385790ae702f8ef481354f16..38a7e4d8420fa34f2754e1e184002bff81eaad95 100644
--- a/dune/tectonic/globalnonlinearity.hh
+++ b/dune/tectonic/globalnonlinearity.hh
@@ -8,7 +8,7 @@
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/matrixindexset.hh>
 
-#include "nicefunction.hh"
+#include "frictionpotential.hh"
 #include "localnonlinearity.hh"
 
 namespace Dune {
diff --git a/dune/tectonic/globalruinanonlinearity.hh b/dune/tectonic/globalruinanonlinearity.hh
index bf7565f02c7d6acacda20fb35394d5afafb10694..0d9350c0a9b59cf20b587d2f0dda49133bee0eec 100644
--- a/dune/tectonic/globalruinanonlinearity.hh
+++ b/dune/tectonic/globalruinanonlinearity.hh
@@ -11,7 +11,7 @@
 
 #include "globalnonlinearity.hh"
 #include "localnonlinearity.hh"
-#include "nicefunction.hh"
+#include "frictionpotential.hh"
 
 namespace Dune {
 template <class MatrixType, class VectorType>
@@ -32,7 +32,7 @@ class GlobalRuinaNonlinearity
       restrictions[i] = nodalIntegrals[i] == 0
                             ? trivialNonlinearity
                             : make_shared<LocalNonlinearity<dim> const>(
-                                  make_shared<RuinaFunction const>(
+                                  make_shared<FrictionPotential const>(
                                       nodalIntegrals[i], fd, state[i]));
     }
   }
diff --git a/dune/tectonic/localnonlinearity.hh b/dune/tectonic/localnonlinearity.hh
index ede2c21bac936451a199713d4ff0ed338e5e9627..eef910c7cc93bba932a64f1922d83e3e13d09f0b 100644
--- a/dune/tectonic/localnonlinearity.hh
+++ b/dune/tectonic/localnonlinearity.hh
@@ -10,7 +10,7 @@
 
 #include <dune/fufem/interval.hh>
 
-#include "nicefunction.hh"
+#include "frictionpotential.hh"
 
 namespace Dune {
 template <int dimension> class LocalNonlinearity {
@@ -18,7 +18,8 @@ template <int dimension> class LocalNonlinearity {
   using VectorType = FieldVector<double, dimension>;
   using MatrixType = FieldMatrix<double, dimension, dimension>;
 
-  LocalNonlinearity(shared_ptr<NiceFunction const> func) : func(func) {}
+  LocalNonlinearity(shared_ptr<FrictionPotentialWrapper const> func)
+      : func(func) {}
 
   std::vector<double> const &get_kinks() const { return func->get_kinks(); }
 
@@ -149,7 +150,7 @@ template <int dimension> class LocalNonlinearity {
   }
 
 private:
-  shared_ptr<NiceFunction const> const func;
+  shared_ptr<FrictionPotentialWrapper const> const func;
 };
 }
 #endif
diff --git a/src/Makefile.am b/src/Makefile.am
index a5112384f496bb7f8c794aee0b34a92751992192..eb5768776b6d8bbdddeb43ed29556297098b7dcc 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,5 +1,3 @@
-SUBDIRS = tests
-
 bin_PROGRAMS = \
 	one-body-sample-2D \
 	one-body-sample-3D
diff --git a/src/tests/Makefile.am b/src/tests/Makefile.am
deleted file mode 100644
index f203d6dd0371acaa0a49a310282659dba38b7e68..0000000000000000000000000000000000000000
--- a/src/tests/Makefile.am
+++ /dev/null
@@ -1,67 +0,0 @@
-check_PROGRAMS = \
-	test-circle-1 \
-	test-circle-10 \
-	test-gradient-horrible \
-	test-gradient-horrible-logarithmic \
-	test-gradient-identity \
-	test-gradient-kinks \
-	test-gradient-parabola \
-	test-gradient-sample \
-	test-gradient-sample-3d \
-	test-gradient-sample-nonsmooth \
-	test-gradient-sample-steep \
-	test-gradient-sample-steep2 \
-	test-gradient-sample-verysteep \
-	test-gradient-sample2 \
-	test-gradient-trivial \
-	test-minimise2
-
-test_circle_1_SOURCES                      = test-circle.cc
-test_circle_1_CPPFLAGS                     = $(AM_CPPFLAGS) -DDUNE_TECTONIC_TEST_CIRCLE_SCALE=1
-test_circle_10_SOURCES                     = test-circle.cc
-test_circle_10_CPPFLAGS                    = $(AM_CPPFLAGS) -DDUNE_TECTONIC_TEST_CIRCLE_SCALE=10
-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
-test_gradient_kinks_SOURCES                = test-gradient-kinks.cc
-test_gradient_parabola_SOURCES             = test-gradient-parabola.cc
-test_gradient_sample_SOURCES               = test-gradient-sample.cc
-test_gradient_sample_3d_SOURCES            = test-gradient-sample-3d.cc
-test_gradient_sample_nonsmooth_SOURCES     = test-gradient-sample-nonsmooth.cc
-test_gradient_sample_steep_SOURCES         = test-gradient-sample-steep.cc
-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
-
-TESTS = $(check_PROGRAMS)
-
-AM_CXXFLAGS = \
-	-Wall \
-	-Wextra \
-	-Wno-c++11-compat \
-	-Wno-c++11-extensions \
-	-Wno-deprecated-declarations \
-	-Wno-empty-body \
-	-Wno-mismatched-tags \
-	-Wno-missing-declarations \
-	-Wno-overloaded-virtual \
-	-Wno-reorder \
-	-Wno-sign-compare \
-	-Wno-tautological-compare \
-	-Wno-type-limits \
-	-Wno-unused-parameter \
-	-Wno-unused-variable \
-	-Wno-unneeded-internal-declaration
-AM_CPPFLAGS = \
-	-DDUNE_FMatrix_WITH_CHECKING \
-	$(DUNE_CPPFLAGS) \
-	-I$(top_srcdir)
-
-# The libraries have to be given in reverse order (most basic libraries
-# last).
-LDADD = \
-	$(DUNE_LDFLAGS) $(DUNE_LIBS)
-AM_LDFLAGS = \
-	$(DUNE_LDFLAGS)
diff --git a/src/tests/test-circle.cc b/src/tests/test-circle.cc
deleted file mode 100644
index 97d19abfaa8ea452dab4d6c46c520a086a9caf7e..0000000000000000000000000000000000000000
--- a/src/tests/test-circle.cc
+++ /dev/null
@@ -1,73 +0,0 @@
-/* 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;
-  using Functional = Dune::EllipticEnergy<dim>;
-  using SmallMatrix = Functional::SmallMatrix;
-  using SmallVector = Functional::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);
-  }
-}
diff --git a/src/tests/test-gradient-horrible-logarithmic.cc b/src/tests/test-gradient-horrible-logarithmic.cc
deleted file mode 100644
index bd2ab10630f3d9695f860f35d760c8cd3f895027..0000000000000000000000000000000000000000
--- a/src/tests/test-gradient-horrible-logarithmic.cc
+++ /dev/null
@@ -1,49 +0,0 @@
-/* Checks if the algorithm converges regardless of where it starts */
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <cassert>
-
-#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;
-  using Functional = Dune::EllipticEnergy<dim>;
-  using SmallMatrix = Functional::SmallMatrix;
-  using SmallVector = Functional::SmallVector;
-
-  SmallMatrix const A = { { 3, 1.5 }, { 1.5, 4 } };
-  SmallVector const b = { 1, 2 };
-
-  auto const f = Dune::make_shared<Dune::HorribleFunctionLogarithmic const>();
-  auto const phi = Dune::make_shared<Functional::NonlinearityType const>(f);
-  Functional const J(A, b, phi);
-
-  double ret1;
-  {
-    SmallVector start = { 17, 34 };
-    ret1 = functionTester(J, start, 6);
-  }
-
-  double ret2;
-  {
-    // Something random
-    SmallVector start = { 279, -96 };
-    ret2 = functionTester(J, start, 12);
-  }
-  assert(std::abs(ret1 - ret2) < 1e-5);
-
-  double ret3;
-  {
-    SmallVector start = { 0, 0 };
-    ret3 = functionTester(J, start, 4);
-  }
-  assert(std::abs(ret1 - ret3) < 1e-5);
-}
diff --git a/src/tests/test-gradient-horrible.cc b/src/tests/test-gradient-horrible.cc
deleted file mode 100644
index 8e912b7982c5f2b3045b7464c7e240df598e4a91..0000000000000000000000000000000000000000
--- a/src/tests/test-gradient-horrible.cc
+++ /dev/null
@@ -1,49 +0,0 @@
-/* Checks if the algorithm converges regardless of where it starts */
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <cassert>
-
-#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;
-  using Functional = Dune::EllipticEnergy<dim>;
-  using SmallMatrix = Functional::SmallMatrix;
-  using SmallVector = Functional::SmallVector;
-
-  SmallMatrix const A = { { 3, 1.5 }, { 1.5, 4 } };
-  SmallVector const b = { 1, 2 };
-
-  auto const f = Dune::make_shared<Dune::HorribleFunction const>();
-  auto const phi = Dune::make_shared<Functional::NonlinearityType const>(f);
-  Functional const J(A, b, phi);
-
-  double ret1;
-  {
-    SmallVector start = { 17, 34 };
-    ret1 = functionTester(J, start, 7);
-  }
-
-  double ret2;
-  {
-    // Something random
-    SmallVector start = { 279, -96 };
-    ret2 = functionTester(J, start, 8);
-  }
-  assert(std::abs(ret1 - ret2) < 1e-5);
-
-  double ret3;
-  {
-    SmallVector start = { 0, 0 };
-    ret3 = functionTester(J, start, 4);
-  }
-  assert(std::abs(ret1 - ret3) < 1e-5);
-}
diff --git a/src/tests/test-gradient-identity.cc b/src/tests/test-gradient-identity.cc
deleted file mode 100644
index 906ee5324a67282efa8a729a32f01c89d8a85534..0000000000000000000000000000000000000000
--- a/src/tests/test-gradient-identity.cc
+++ /dev/null
@@ -1,49 +0,0 @@
-/* Checks if the algorithm converges regardless of where it starts */
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <cassert>
-
-#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;
-  using Functional = Dune::EllipticEnergy<dim>;
-  using SmallMatrix = Functional::SmallMatrix;
-  using SmallVector = Functional::SmallVector;
-
-  SmallMatrix const A = { { 3, 1.5 }, { 1.5, 4 } };
-  SmallVector const b = { 1, 2 };
-
-  auto const f = Dune::make_shared<Dune::LinearFunction const>(1);
-  auto const phi = Dune::make_shared<Functional::NonlinearityType const>(f);
-  Functional const J(A, b, phi);
-
-  double ret1;
-  {
-    SmallVector start = { 17, 34 };
-    ret1 = functionTester(J, start, 6);
-  }
-
-  double ret2;
-  {
-    // Something random
-    SmallVector start = { 279, -96 };
-    ret2 = functionTester(J, start, 10);
-  }
-  assert(std::abs(ret1 - ret2) < 1e-5);
-
-  double ret3;
-  {
-    SmallVector start = { 0, 0 };
-    ret3 = functionTester(J, start, 3);
-  }
-  assert(std::abs(ret1 - ret3) < 1e-5);
-}
diff --git a/src/tests/test-gradient-kinks.cc b/src/tests/test-gradient-kinks.cc
deleted file mode 100644
index 857cd07de2e29818ff77d1c312520e63a3f572da..0000000000000000000000000000000000000000
--- a/src/tests/test-gradient-kinks.cc
+++ /dev/null
@@ -1,57 +0,0 @@
-/* Checks if the algorithm converges regardless of where it starts */
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <cassert>
-
-#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;
-  using Functional = Dune::EllipticEnergy<dim>;
-  using SmallMatrix = Functional::SmallMatrix;
-  using SmallVector = Functional::SmallVector;
-
-  SmallMatrix const A = { { 4, 1.5 }, { 1.5, 3 } };
-  SmallVector const b = { -8, 14 };
-
-  auto const f = Dune::make_shared<Dune::ThreeKinkFunction const>();
-  auto const phi = Dune::make_shared<Functional::NonlinearityType const>(f);
-  Functional const J(A, b, phi);
-
-  size_t runs = 5;
-
-  double ret1;
-  {
-    SmallVector start = { 17, 34 };
-    ret1 = functionTester(J, start, runs);
-  }
-
-  double ret2;
-  {
-    SmallVector start = { 1, 0 };
-    ret2 = functionTester(J, start, runs);
-  }
-  assert(std::abs(ret1 - ret2) < 1e-5);
-
-  double ret3;
-  {
-    SmallVector start = { 0, 1 };
-    ret3 = functionTester(J, start, runs);
-  }
-  assert(std::abs(ret1 - ret3) < 1e-5);
-
-  double ret4;
-  {
-    SmallVector start = { 0, 0 };
-    ret4 = functionTester(J, start, runs);
-  }
-  assert(std::abs(ret1 - ret4) < 1e-5);
-}
diff --git a/src/tests/test-gradient-method-helper.hh b/src/tests/test-gradient-method-helper.hh
deleted file mode 100644
index 310e60d82a92cddf0744e481c0e3ff67fce0f4e7..0000000000000000000000000000000000000000
--- a/src/tests/test-gradient-method-helper.hh
+++ /dev/null
@@ -1,56 +0,0 @@
-#ifndef TEST_GRADIENT_METHOD_HELPER_HH
-#define TEST_GRADIENT_METHOD_HELPER_HH
-
-#include <boost/format.hpp>
-
-#include <dune/tnnmg/problem-classes/bisection.hh>
-
-#include <dune/tectonic/ellipticenergy.hh>
-#include <dune/tectonic/minimisation.hh>
-#include <dune/tectonic/minimisation2.hh>
-
-template <int dim>
-double functionTester(Dune::EllipticEnergy<dim> J,
-                      typename Dune::EllipticEnergy<dim>::SmallVector &start,
-                      size_t runs) {
-  Bisection const bisection(
-      0.0,   // acceptError: Stop if the search interval has
-             //              become smaller than this number
-      1.0,   // acceptFactor: ?
-      1e-12, // requiredResidual: ?
-      true,  // fastQuadratic
-      0);    // safety: acceptance factor for inexact minimization
-  double const original = J(start);
-  Dune::minimise(J, start, runs, bisection);
-  double const final = J(start);
-  std::cout << boost::format("%8g -> %e (%e)") % original % final % start
-            << std::endl;
-  return final;
-}
-
-template <int dim>
-double functionTester2(Dune::EllipticEnergy<dim> J,
-                       typename Dune::EllipticEnergy<dim>::SmallVector &start,
-                       size_t runs) {
-  Bisection const bisection(
-      0.0,   // acceptError: Stop if the search interval has
-             //              become smaller than this number
-      1.0,   // acceptFactor: ?
-      1e-12, // requiredResidual: ?
-      true,  // fastQuadratic
-      0);    // safety: acceptance factor for inexact minimization
-  Dune::minimise2(J, start, runs, bisection);
-  double const final = J(start);
-  std::cout << boost::format("         -> %e (%e)") % final % start
-            << std::endl;
-  return final;
-}
-
-template <int dim>
-double two_distance(typename Dune::EllipticEnergy<dim>::SmallVector const &x,
-                    typename Dune::EllipticEnergy<dim>::SmallVector const &y) {
-  typename Dune::EllipticEnergy<dim>::SmallVector tmp = x;
-  tmp -= y;
-  return tmp.two_norm();
-}
-#endif
diff --git a/src/tests/test-gradient-method-nicefunction.hh b/src/tests/test-gradient-method-nicefunction.hh
deleted file mode 100644
index 4eaba7c081597d75f96b93c7383f708c095817d4..0000000000000000000000000000000000000000
--- a/src/tests/test-gradient-method-nicefunction.hh
+++ /dev/null
@@ -1,176 +0,0 @@
-#ifndef TEST_GRADIENT_METHOD_NICEFUNCTION_HH
-#define TEST_GRADIENT_METHOD_NICEFUNCTION_HH
-
-#include <cassert>
-
-#include <dune/tectonic/nicefunction.hh>
-
-namespace Dune {
-class MyFunction : public NiceFunction {
-protected:
-  MyFunction(std::vector<double> const &_kinks) : NiceFunction(_kinks) {}
-
-public:
-  MyFunction() : NiceFunction() {}
-
-  double virtual second_deriv(double) const {
-    assert(false);
-    return 0;
-  }
-
-  double virtual regularity(double) const {
-    assert(false);
-    return 0;
-  }
-};
-
-class Parabola : public MyFunction {
-  double virtual evaluate(double x) const { return x * x; }
-
-  double virtual leftDifferential(double s) const { return 2 * s; }
-
-  double virtual rightDifferential(double s) const { return 2 * s; }
-
-  double virtual second_deriv(double s) const { return 2; }
-
-  double virtual regularity(double s) const { return 2; }
-
-  bool virtual smoothesNorm() const { return true; }
-};
-
-class LinearFunction : public MyFunction {
-public:
-  LinearFunction(double a) : MyFunction(), coefficient(a) {}
-
-  double virtual evaluate(double x) const { return coefficient * x; }
-
-  double virtual leftDifferential(double s) const { return coefficient; }
-
-  double virtual rightDifferential(double s) const { return coefficient; }
-
-  double virtual second_deriv(double) const { return 0; }
-
-  double virtual regularity(double s) const { return 0; }
-
-private:
-  double const coefficient;
-};
-
-template <int slope> class SampleFunction : public MyFunction {
-public:
-  SampleFunction() : MyFunction({ 1 }) {}
-
-  double virtual evaluate(double x) const {
-    return (x < 1) ? x : (slope * (x - 1) + 1);
-  }
-
-  double virtual leftDifferential(double s) const {
-    return (s <= 1) ? 1 : slope;
-  }
-
-  double virtual rightDifferential(double s) const {
-    return (s < 1) ? 1 : slope;
-  }
-};
-
-// slope in [n-1,n] is n
-class HorribleFunction : public MyFunction {
-  // TODO: Handle kinks
-
-public:
-  double virtual evaluate(double x) const {
-    double const fl = floor(x);
-    double const sum = fl * (fl + 1) / 2;
-    return sum + (fl + 1) * (x - fl);
-  }
-
-  double virtual leftDifferential(double x) const {
-    double const fl = floor(x);
-    if (x - fl < 1e-14)
-      return fl;
-    else
-      return fl + 1;
-  }
-
-  double virtual rightDifferential(double x) const {
-    double const c = ceil(x);
-    if (c - x < 1e-14)
-      return c + 1;
-    else
-      return c;
-  }
-};
-
-// slope in [n-1,n] is log(n+1)
-class HorribleFunctionLogarithmic : public MyFunction {
-public:
-  // TODO: Handle kinks
-
-  double virtual evaluate(double x) const {
-    double y = 0;
-    size_t const fl = floor(x);
-    for (size_t i = 1; i <= fl;)
-      y += std::log(
-          ++i); // factorials grow to fast so we compute this incrementally
-
-    return y + std::log(fl + 2) * (x - fl);
-  }
-
-  double virtual leftDifferential(double x) const {
-    double const fl = floor(x);
-    if (x - fl < 1e-14)
-      return std::log(fl + 1);
-    else
-      return std::log(fl + 2);
-  }
-
-  double virtual rightDifferential(double x) const {
-    double const c = ceil(x);
-    if (c - x < 1e-14)
-      return std::log(c + 2);
-    else
-      return std::log(c + 1);
-  }
-};
-
-class ThreeKinkFunction : public MyFunction {
-public:
-  ThreeKinkFunction() : MyFunction({ 5, 10, 15 }) {}
-
-  double virtual evaluate(double x) const {
-    double acc = 0;
-    double last_kink = 0;
-    int i;
-    auto const &kinks = this->get_kinks();
-    for (i = 0; i < kinks.size(); ++i) {
-      if (x <= kinks[i])
-        break;
-
-      acc += (i + 1) * (kinks[i] - last_kink);
-      last_kink = kinks[i];
-    }
-    return acc + (i + 1) * (x - last_kink);
-  }
-
-  double virtual leftDifferential(double x) const {
-    if (x <= 5)
-      return 1;
-    if (x <= 10)
-      return 2;
-    if (x <= 15)
-      return 3;
-    return 4;
-  }
-
-  double virtual rightDifferential(double x) const {
-    if (x < 5)
-      return 1;
-    if (x < 10)
-      return 2;
-    if (x < 15)
-      return 3;
-    return 4;
-  }
-};
-}
-#endif
diff --git a/src/tests/test-gradient-parabola.cc b/src/tests/test-gradient-parabola.cc
deleted file mode 100644
index 23f4594eb55de3db3c784acbd05f5c4461da88ef..0000000000000000000000000000000000000000
--- a/src/tests/test-gradient-parabola.cc
+++ /dev/null
@@ -1,70 +0,0 @@
-/* Checks if the descent direction is computed correctly using the
-   analytic solution; also checks if the algorithm converges
-   to the right solution regardless of where it starts */
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <cassert>
-
-#include <dune/common/shared_ptr.hh>
-
-#include <dune/tectonic/ellipticenergy.hh>
-
-#include "test-gradient-method-helper.hh"
-#include "test-gradient-method-nicefunction.hh"
-
-int main() {
-  int const dim = 2;
-  using Functional = Dune::EllipticEnergy<dim>;
-  using SmallMatrix = Functional::SmallMatrix;
-  using SmallVector = Functional::SmallVector;
-
-  SmallMatrix const A = { { 3, 1.5 }, { 1.5, 4 } };
-  SmallVector const b = { 1, 2 };
-
-  // |x|^2 as the nonlinearity is the same as having no nonlinearity
-  // but twice the identity matrix added to A. In other words, we're
-  // solving A + 2*id = b
-  auto const f = Dune::make_shared<Dune::Parabola const>();
-  auto const phi = Dune::make_shared<Functional::NonlinearityType const>(f);
-  Functional const J(A, b, phi);
-
-  SmallVector solution = { 4.0 / 37.0, 34.0 / 111.0 }; // Analytic solution
-  SmallMatrix M = A;
-  M[0][0] += 2;
-  M[1][1] += 2;
-
-  double ret1;
-  {
-    SmallVector start = b;
-    start *= 17;
-
-    SmallVector analytic_descent = b;
-    M.mmv(start, analytic_descent);
-
-    SmallVector numerical_descent;
-    J.descentDirection(start, numerical_descent);
-    assert(two_distance<dim>(numerical_descent, analytic_descent) < 1e-10);
-
-    ret1 = functionTester(J, start, 6);
-    assert(two_distance<dim>(start, solution) < 1e-6);
-  }
-
-  double ret2;
-  {
-    // Something random
-    SmallVector start = { 279, -96 };
-    SmallVector analytic_descent = b;
-    M.mmv(start, analytic_descent);
-
-    SmallVector numerical_descent;
-    J.descentDirection(start, numerical_descent);
-    assert(two_distance<dim>(numerical_descent, analytic_descent) < 1e-10);
-
-    ret2 = functionTester(J, start, 15);
-    assert(two_distance<dim>(start, solution) < 1e-6);
-  }
-  assert(std::abs(ret1 - ret2) < 1e-11);
-}
diff --git a/src/tests/test-gradient-sample-3d.cc b/src/tests/test-gradient-sample-3d.cc
deleted file mode 100644
index 7b92869524233f9b9db716164626725be5f5242a..0000000000000000000000000000000000000000
--- a/src/tests/test-gradient-sample-3d.cc
+++ /dev/null
@@ -1,49 +0,0 @@
-/* Checks if the algorithm converges regardless of where it starts */
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <cassert>
-
-#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 = 3;
-  using Functional = Dune::EllipticEnergy<dim>;
-  using SmallMatrix = Functional::SmallMatrix;
-  using SmallVector = Functional::SmallVector;
-
-  SmallMatrix const A = { { 3, 1.5, 0 }, { 1.5, 4, 1.5 }, { 0, 1.5, 5 } };
-  SmallVector const b = { 1, 2, 3 };
-
-  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);
-
-  double ret1;
-  {
-    SmallVector start = { 17, 34, 51 };
-    ret1 = functionTester(J, start, 9);
-  }
-
-  double ret2;
-  {
-    // Something random
-    SmallVector start = { 279, -96, -15 };
-    ret2 = functionTester(J, start, 15);
-  }
-  assert(std::abs(ret1 - ret2) < 1e-5);
-
-  double ret3;
-  {
-    SmallVector start = { 0, 0, 0 };
-    ret3 = functionTester(J, start, 5);
-  }
-  assert(std::abs(ret1 - ret3) < 1e-5);
-}
diff --git a/src/tests/test-gradient-sample-nonsmooth.cc b/src/tests/test-gradient-sample-nonsmooth.cc
deleted file mode 100644
index e8ac1f28ad7a71f275e5aac9ceeb3e8ecf710b0b..0000000000000000000000000000000000000000
--- a/src/tests/test-gradient-sample-nonsmooth.cc
+++ /dev/null
@@ -1,70 +0,0 @@
-/* Checks if the descent direction is computed correctly
-   using the analytic solution */
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <cassert>
-
-#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;
-  using Functional = Dune::EllipticEnergy<dim>;
-  using SmallMatrix = Functional::SmallMatrix;
-  using SmallVector = Functional::SmallVector;
-
-  SmallMatrix const A = { { 3, 1.5 }, { 1.5, 4 } };
-  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);
-
-  SmallVector start;
-  /*
-    for x = b/|b|:
-
-    j_(x)
-    = Ax - b + 1/|x| x
-    = 1/sqrt(5) (6, 9.5) - (1, 2) + 1/sqrt(5) (1, 2)
-    = (7/sqrt(5) - 1, 11.5/sqrt(5) - 2)
-
-    j+(x)
-    = Ax - b + 2/|x| x
-    = 1/sqrt(5) (6, 9.5) - (1, 2) + 2/sqrt(5) (1, 2)
-    = (8/sqrt(5) - 1, 13.5/sqrt(5) - 2)
-  */
-  {
-    start = b;
-    start /= (b.two_norm() + 1e-12); // Make sure the norm is below 1;
-    assert(start.two_norm() < 1);
-
-    SmallVector numerical_descent;
-    J.descentDirection(start, numerical_descent);
-    SmallVector analytic_descent = { -(7 / std::sqrt(5) - 1),
-                                     -(11.5 / std::sqrt(5) - 2) };
-    assert(two_distance<dim>(analytic_descent, numerical_descent) < 1e-10);
-
-    functionTester(J, start, 6);
-  }
-  {
-    start = b;
-    start /= (b.two_norm() - 1e-12); // Make sure the norm is above 1;
-    assert(start.two_norm() > 1);
-
-    SmallVector numerical_descent;
-    J.descentDirection(start, numerical_descent);
-    SmallVector analytic_descent = { -(8 / std::sqrt(5) - 1),
-                                     -(13.5 / std::sqrt(5) - 2) };
-    assert(two_distance<dim>(analytic_descent, numerical_descent) < 1e-10);
-
-    functionTester(J, start, 6);
-  }
-}
diff --git a/src/tests/test-gradient-sample-steep.cc b/src/tests/test-gradient-sample-steep.cc
deleted file mode 100644
index 47b87d55fc921bf1e4ba2db4168e52667de45e3b..0000000000000000000000000000000000000000
--- a/src/tests/test-gradient-sample-steep.cc
+++ /dev/null
@@ -1,58 +0,0 @@
-/* Checks if the algorithm converges regardless of where it starts */
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <cassert>
-
-#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;
-  using Functional = Dune::EllipticEnergy<dim>;
-  using SmallMatrix = Functional::SmallMatrix;
-  using SmallVector = Functional::SmallVector;
-
-  SmallMatrix const A = { { 1, 0 }, { 0, 1 } };
-  SmallVector const b = { 1, 2.5 };
-
-  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);
-
-  double ret1;
-  {
-    SmallVector start = { 0, 1 };
-    ret1 = functionTester(J, start, 2);
-    SmallVector minimum;
-    functionTester2(J, minimum, 2);
-    assert(two_distance<dim>(start, minimum) < 1e-5);
-  }
-
-  double ret2;
-  {
-    // Something random
-    SmallVector start = { 279, -96 };
-    ret2 = functionTester(J, start, 3);
-    SmallVector minimum;
-    functionTester2(J, minimum, 3);
-    assert(two_distance<dim>(start, minimum) < 1e-5);
-  }
-  assert(std::abs(ret1 - ret2) < 1e-5);
-
-  double ret3;
-  {
-    SmallVector start = { 0, 0 };
-    ret3 = functionTester(J, start, 1);
-    SmallVector minimum;
-    functionTester2(J, minimum, 1);
-    assert(two_distance<dim>(start, minimum) < 1e-5);
-  }
-  assert(std::abs(ret1 - ret3) < 1e-5);
-}
diff --git a/src/tests/test-gradient-sample-steep2.cc b/src/tests/test-gradient-sample-steep2.cc
deleted file mode 100644
index 47b87d55fc921bf1e4ba2db4168e52667de45e3b..0000000000000000000000000000000000000000
--- a/src/tests/test-gradient-sample-steep2.cc
+++ /dev/null
@@ -1,58 +0,0 @@
-/* Checks if the algorithm converges regardless of where it starts */
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <cassert>
-
-#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;
-  using Functional = Dune::EllipticEnergy<dim>;
-  using SmallMatrix = Functional::SmallMatrix;
-  using SmallVector = Functional::SmallVector;
-
-  SmallMatrix const A = { { 1, 0 }, { 0, 1 } };
-  SmallVector const b = { 1, 2.5 };
-
-  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);
-
-  double ret1;
-  {
-    SmallVector start = { 0, 1 };
-    ret1 = functionTester(J, start, 2);
-    SmallVector minimum;
-    functionTester2(J, minimum, 2);
-    assert(two_distance<dim>(start, minimum) < 1e-5);
-  }
-
-  double ret2;
-  {
-    // Something random
-    SmallVector start = { 279, -96 };
-    ret2 = functionTester(J, start, 3);
-    SmallVector minimum;
-    functionTester2(J, minimum, 3);
-    assert(two_distance<dim>(start, minimum) < 1e-5);
-  }
-  assert(std::abs(ret1 - ret2) < 1e-5);
-
-  double ret3;
-  {
-    SmallVector start = { 0, 0 };
-    ret3 = functionTester(J, start, 1);
-    SmallVector minimum;
-    functionTester2(J, minimum, 1);
-    assert(two_distance<dim>(start, minimum) < 1e-5);
-  }
-  assert(std::abs(ret1 - ret3) < 1e-5);
-}
diff --git a/src/tests/test-gradient-sample-verysteep.cc b/src/tests/test-gradient-sample-verysteep.cc
deleted file mode 100644
index 06af4323182c427df2b4821b71ee4b1d3cd9bb7f..0000000000000000000000000000000000000000
--- a/src/tests/test-gradient-sample-verysteep.cc
+++ /dev/null
@@ -1,58 +0,0 @@
-/* Checks if the algorithm converges regardless of where it starts */
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <cassert>
-
-#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;
-  using Functional = Dune::EllipticEnergy<dim>;
-  using SmallMatrix = Functional::SmallMatrix;
-  using SmallVector = Functional::SmallVector;
-
-  SmallMatrix const A = { { 1, 0 }, { 0, 1 } };
-  SmallVector const b = { 1, 2.5 };
-
-  auto const f = Dune::make_shared<Dune::SampleFunction<100> const>();
-  auto const phi = Dune::make_shared<Functional::NonlinearityType const>(f);
-  Functional const J(A, b, phi);
-
-  double ret1;
-  {
-    SmallVector start = { 0, 1 };
-    ret1 = functionTester(J, start, 2);
-    SmallVector minimum;
-    functionTester2(J, minimum, 2);
-    assert(two_distance<dim>(start, minimum) < 1e-5);
-  }
-
-  double ret2;
-  {
-    // Something random
-    SmallVector start = { 279, -96 };
-    ret2 = functionTester(J, start, 4);
-    SmallVector minimum;
-    functionTester2(J, minimum, 4);
-    assert(two_distance<dim>(start, minimum) < 1e-5);
-  }
-  assert(std::abs(ret1 - ret2) < 1e-8);
-
-  double ret3;
-  {
-    SmallVector start = { 0, 0 };
-    ret3 = functionTester(J, start, 1);
-    SmallVector minimum;
-    functionTester2(J, minimum, 1);
-    assert(two_distance<dim>(start, minimum) < 1e-5);
-  }
-  assert(std::abs(ret1 - ret3) < 1e-5);
-}
diff --git a/src/tests/test-gradient-sample.cc b/src/tests/test-gradient-sample.cc
deleted file mode 100644
index 938c82bc9555b9c33bbc409228eaee093c9d1fe3..0000000000000000000000000000000000000000
--- a/src/tests/test-gradient-sample.cc
+++ /dev/null
@@ -1,72 +0,0 @@
-/* Checks if the descent direction is computed correctly using the
-   analytic solution; also checks if the algorithm converges
-   regardless of where it starts */
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <cassert>
-
-#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;
-  using Functional = Dune::EllipticEnergy<dim>;
-  using SmallMatrix = Functional::SmallMatrix;
-  using SmallVector = Functional::SmallVector;
-
-  SmallMatrix const A = { { 3, 1.5 }, { 1.5, 4 } };
-  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);
-
-  /*
-    j(x)
-    = Ax - b + 2/|x| x
-    = 17*(6, 9.5) - (1, 2) + 2/sqrt(5) (1, 2)
-    = (102 - 1 + 2/sqrt(5), 161.5 - 2 + 4/sqrt(5))
-  */
-  SmallVector analytic_descent = { -(102 - 1 + 2 / std::sqrt(5)),
-                                   -(161.5 - 2 + 4 / std::sqrt(5)) };
-  SmallVector numerical_descent;
-  J.descentDirection(SmallVector({ 17, 34 }), numerical_descent);
-  assert(two_distance<dim>(analytic_descent, numerical_descent) < 1e-10);
-
-  double ret1;
-  {
-    SmallVector start = { 17, 34 };
-    ret1 = functionTester(J, start, 8);
-    SmallVector minimum;
-    functionTester2(J, minimum, 8);
-    assert(two_distance<dim>(start, minimum) < 1e-5);
-  }
-
-  double ret2;
-  {
-    // Something random
-    SmallVector start = { 279, -96 };
-    ret2 = functionTester(J, start, 10);
-    SmallVector minimum;
-    functionTester2(J, minimum, 10);
-    assert(two_distance<dim>(start, minimum) < 1e-5);
-  }
-  assert(std::abs(ret1 - ret2) < 1e-5);
-
-  double ret3;
-  {
-    SmallVector start = { 0, 0 };
-    ret3 = functionTester(J, start, 3);
-    SmallVector minimum;
-    functionTester2(J, minimum, 3);
-    assert(two_distance<dim>(start, minimum) < 1e-5);
-  }
-  assert(std::abs(ret1 - ret3) < 1e-5);
-}
diff --git a/src/tests/test-gradient-sample2.cc b/src/tests/test-gradient-sample2.cc
deleted file mode 100644
index 8b09f463ec82ee505fae46d081c81c56626bedce..0000000000000000000000000000000000000000
--- a/src/tests/test-gradient-sample2.cc
+++ /dev/null
@@ -1,51 +0,0 @@
-/* Checks if the algorithm converges regardless of where it starts;
-   in particular, that converging in one step does not cause us to
-   diverge again */
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <cassert>
-
-#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;
-  using Functional = Dune::EllipticEnergy<dim>;
-  using SmallMatrix = Functional::SmallMatrix;
-  using SmallVector = Functional::SmallVector;
-
-  SmallMatrix const A = { { 1, 0 }, { 0, 1 } };
-  SmallVector const b = { 1, 1 };
-
-  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);
-
-  double ret1;
-  {
-    SmallVector start = b;
-    ret1 = functionTester(J, start, 2);
-  }
-
-  double ret2;
-  {
-    // Something random
-    SmallVector start = { 279, -96 };
-    ret2 = functionTester(J, start, 7);
-  }
-  assert(std::abs(ret1 - ret2) < 1e-5);
-
-  double ret3;
-  {
-    SmallVector start = { 0, 0 };
-    ret3 = functionTester(J, start, 2);
-  }
-  assert(std::abs(ret1 - ret3) < 1e-5);
-}
diff --git a/src/tests/test-gradient-trivial.cc b/src/tests/test-gradient-trivial.cc
deleted file mode 100644
index 044d0c0dff6c79d76b83fed91249120592840fb9..0000000000000000000000000000000000000000
--- a/src/tests/test-gradient-trivial.cc
+++ /dev/null
@@ -1,62 +0,0 @@
-/* Checks if the descent direction is computed correctly using the
-   analytic solution; also checks if the algorithm converges
-   to the right solution regardless of where it starts */
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <cassert>
-
-#include <dune/common/shared_ptr.hh>
-
-#include <dune/tectonic/ellipticenergy.hh>
-
-#include "test-gradient-method-helper.hh"
-
-int main() {
-  int const dim = 2;
-  using Functional = Dune::EllipticEnergy<dim>;
-  using SmallMatrix = Functional::SmallMatrix;
-  using SmallVector = Functional::SmallVector;
-
-  SmallMatrix const A = { { 3, 1.5 }, { 1.5, 4 } };
-  SmallVector const b = { 1, 2 };
-
-  auto const f = Dune::make_shared<Dune::TrivialFunction const>();
-  auto const phi = Dune::make_shared<Functional::NonlinearityType const>(f);
-  Functional const J(A, b, phi);
-
-  SmallVector solution = { 4.0 / 39.0, 6.0 / 13.0 }; // Analytic solution
-
-  double ret1;
-  {
-    SmallVector start = b;
-    start *= 17;
-
-    SmallVector analytic_descent = b;
-    A.mmv(start, analytic_descent);
-    SmallVector numerical_descent;
-    J.descentDirection(start, numerical_descent);
-    assert(two_distance<dim>(numerical_descent, analytic_descent) < 1e-10);
-
-    ret1 = functionTester(J, start, 6);
-    assert(two_distance<dim>(start, solution) < 1e-5);
-  }
-
-  double ret2;
-  {
-    // Something random
-    SmallVector start = { 279, -96 };
-
-    SmallVector analytic_descent = b;
-    A.mmv(start, analytic_descent);
-    SmallVector numerical_descent;
-    J.descentDirection(start, numerical_descent);
-    assert(two_distance<dim>(numerical_descent, analytic_descent) < 1e-10);
-
-    ret2 = functionTester(J, start, 25);
-    assert(two_distance<dim>(start, solution) < 1e-6);
-  }
-  assert(std::abs(ret1 - ret2) < 1e-5);
-}
diff --git a/src/tests/test-minimise2.cc b/src/tests/test-minimise2.cc
deleted file mode 100644
index a84886cb270475a35b3cfdb1282fc53869b81242..0000000000000000000000000000000000000000
--- a/src/tests/test-minimise2.cc
+++ /dev/null
@@ -1,58 +0,0 @@
-/* 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;
-  using Functional = Dune::EllipticEnergy<dim>;
-  using SmallMatrix = Functional::SmallMatrix;
-  using SmallVector = Functional::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);
-  }
-}
diff --git a/src/tests/test_circle_1.m b/src/tests/test_circle_1.m
deleted file mode 100644
index a14afb42ed3fdab8189c6790863092dc277402f8..0000000000000000000000000000000000000000
--- a/src/tests/test_circle_1.m
+++ /dev/null
@@ -1,38 +0,0 @@
-clear all;
-
-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 = 1.0;
-
-t = -5:0.1:5;
-x = scale*cos(t);
-y = scale*sin(t);
-z = arrayfun(@(vec1, vec2) J([vec1; vec2]), x, y);
-
-hold on
-plot3(x,y,z,'k');
-
-title 'scale = 1'
-
-# Minimum taken from test-circle.cc with scale=1
-
-plot3(-1.487905e-01, +9.888687e-01, J([-1.487905e-01; +9.888687e-01]), 'r+');
-
-x = scale * (-1:.05:1);
-y = scale * (-1:.05:1);
-[X, Y] = meshgrid(x,y);
-
-for i=1:length(y)
-  val(i,:) = arrayfun(@(vec1, vec2) J([vec1; vec2]), X(i,:), Y(i,:));
-end
-
-surf(x, y, val)
-hold off;
diff --git a/src/tests/test_circle_10.m b/src/tests/test_circle_10.m
deleted file mode 100644
index 22224cc5b7c1c7ee28e460f56f43431e9f659ac0..0000000000000000000000000000000000000000
--- a/src/tests/test_circle_10.m
+++ /dev/null
@@ -1,40 +0,0 @@
-clear all;
-
-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);
-
-hold on
-plot3(x,y,z,'k');
-
-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])])
-
-x = scale * (-1:.05:1);
-y = scale * (-1:.05:1);
-[X, Y] = meshgrid(x,y);
-
-for i=1:length(y)
-  val(i,:) = arrayfun(@(vec1, vec2) J([vec1; vec2]), X(i,:), Y(i,:));
-end
-
-surf(x, y, val)
-hold off;