diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index cf54f59515f90a58c1e1eb171445be79326b9f85..4ff180ec69020b192786d2596ded80cdf4547c3a 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -5,7 +5,6 @@ set(TESTS
   frameinvariancetest
   globalgfetestfunctionbasistest
   harmonicenergytest
-  interillustration
   localgeodesicfefunctiontest
   localgeodesicfestiffnesstest
   localgfetestfunctiontest
diff --git a/test/interillustration.cc b/test/interillustration.cc
deleted file mode 100644
index 0e05bd186c82ec14fe0cac4087cbf456333fcf27..0000000000000000000000000000000000000000
--- a/test/interillustration.cc
+++ /dev/null
@@ -1,214 +0,0 @@
-#include "config.h"
-
-#include <dune/grid/uggrid.hh>
-#include <dune/grid/geometrygrid.hh>
-#include <dune/grid/io/file/vtk.hh>
-
-#include <dune/localfunctions/lagrange/p1.hh>
-#include <dune/localfunctions/lagrange/p2.hh>
-
-#include <dune/functions/functionspacebases/lagrangebasis.hh>
-#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
-
-#include <dune/gfe/unitvector.hh>
-#include <dune/gfe/localgeodesicfefunction.hh>
-#include <dune/gfe/localprojectedfefunction.hh>
-#include <dune/gfe/localtangentfefunction.hh>
-
-
-using namespace Dune;
-
-/** \brief Encapsulates the grid deformation for the GeometryGrid class */
-template <class HostGridView>
-class DeformationFunction
-: public Dune :: DiscreteCoordFunction< double, 3, DeformationFunction<HostGridView> >
-{
-  typedef DeformationFunction<HostGridView> This;
-  typedef Dune :: DiscreteCoordFunction< double, 3, This > Base;
-
-  static const int dim = HostGridView::dimension;
-
-public:
-
-  DeformationFunction(const HostGridView& gridView,
-                      const std::vector<FieldVector<double,3> >& deformedPosition)
-  : gridView_(gridView),
-  deformedPosition_(deformedPosition)
-  {}
-
-  void evaluate (const typename HostGridView::template Codim<dim>::Entity& hostEntity, unsigned int corner,
-                 Dune::FieldVector<double,3> &y ) const
-  {
-    const typename HostGridView::IndexSet& indexSet = gridView_.indexSet();
-    int idx = indexSet.index(hostEntity);
-    y = deformedPosition_[idx];
-  }
-
-  void evaluate (const typename HostGridView::template Codim<0>::Entity& hostEntity, unsigned int corner,
-                 Dune::FieldVector<double,3> &y ) const
-  {
-    const typename HostGridView::IndexSet& indexSet = gridView_.indexSet();
-    int idx = indexSet.subIndex(hostEntity, corner,dim);
-    y = deformedPosition_[idx];
-  }
-
-private:
-
-  HostGridView gridView_;
-
-  const std::vector<FieldVector<double,3> > deformedPosition_;
-
-};
-
-/** \brief
- *
- * \param partialDerivative The dof with respect to which we are computing the variation
- */
-template <class TargetSpace, class LocalFEFunctionType>
-void interpolate(const LocalFEFunctionType& localGeodesicFEFunction,
-                 int partialDerivative,
-                const std::string& nameSuffix)
-{
-  static const int dim = 2;
-
-  typedef UGGrid<dim> GridType;
-
-  GridFactory<GridType> factory;
-
-  factory.insertVertex({0,0});
-  factory.insertVertex({1,0});
-  factory.insertVertex({0,1});
-
-  factory.insertElement(GeometryTypes::triangle, {0,1,2});
-
-  shared_ptr<GridType> grid = factory.createGrid();
-
-  grid->globalRefine(3);
-
-  typedef GridType::LeafGridView GridView;
-  GridView gridView = grid->leafGridView();
-
-  ///////////////////////////////////////////////////////////////
-  //   Sample the interpolating function on the grid
-  ///////////////////////////////////////////////////////////////
-
-  std::vector<typename TargetSpace::CoordinateType> samples(gridView.size(dim));
-
-  for (auto v=gridView.begin<dim>(); v!=gridView.end<dim>(); ++v)
-    samples[gridView.indexSet().index(*v)] = localGeodesicFEFunction.evaluate(v->geometry().corner(0)).globalCoordinates();
-
-  // Sample a variation vector field
-  std::vector<typename TargetSpace::EmbeddedTangentVector> variation0(gridView.size(dim));
-  std::vector<typename TargetSpace::EmbeddedTangentVector> variation1(gridView.size(dim));
-  for (const auto& v : vertices(gridView))
-  {
-    FieldMatrix<double,3,3> derivative;
-    localGeodesicFEFunction.evaluateDerivativeOfValueWRTCoefficient(v.geometry().corner(0),
-                                                                    partialDerivative,  // select the Lagrange node
-                                                                    derivative);
-
-    Dune::FieldMatrix<double,2,3> basis = localGeodesicFEFunction.coefficient(partialDerivative).orthonormalFrame();
-
-    derivative.mtv(basis[0], variation0[gridView.indexSet().index(v)]);
-    derivative.mtv(basis[1], variation1[gridView.indexSet().index(v)]);
-  }
-
-  // sample a checkerboard pattern for nicer pictures
-  uint pattern = 8;
-  std::vector<double> colors(gridView.size(0));
-  for (auto e=gridView.begin<0>(); e!=gridView.end<0>(); ++e)
-  {
-    FieldVector<double,dim> center = e->geometry().center();
-
-    uint i = pattern * center[0];
-    uint j = pattern * center[1];
-
-    FieldVector<double,dim> local;
-    local[0] = (center[0] - (i/double(pattern)))*pattern;
-    local[1] = (center[1] - (j/double(pattern)))*pattern;
-
-    colors[gridView.indexSet().index(*e)] = (local[0] + local[1] <= 1);
-  }
-
-  ///////////////////////////////////////////////////////////////
-  //   Write a grid with the interpolated positions
-  ///////////////////////////////////////////////////////////////
-
-  typedef Dune::GeometryGrid<GridType,DeformationFunction<typename GridType::LeafGridView> > DeformedGridType;
-
-  DeformationFunction<GridView> deformationFunction(gridView, samples);
-
-  // stupid, can't instantiate deformedGrid with a const grid
-  DeformedGridType deformedGrid(const_cast<GridType&>(*grid), deformationFunction);
-
-  typedef Functions::LagrangeBasis<typename DeformedGridType::LeafGridView, 1> FEBasis;
-  FEBasis feBasis(deformedGrid.leafGridView());
-
-  auto variation0Function = Functions::makeDiscreteGlobalBasisFunction<typename TargetSpace::EmbeddedTangentVector>(feBasis, variation0);
-  auto variation1Function = Functions::makeDiscreteGlobalBasisFunction<typename TargetSpace::EmbeddedTangentVector>(feBasis, variation1);
-
-  Dune::VTKWriter<typename DeformedGridType::LeafGridView> vtkWriter(deformedGrid.leafGridView());
-  vtkWriter.addCellData(colors, "colors");
-
-  vtkWriter.addVertexData(variation0Function, VTK::FieldInfo("variation 0", VTK::FieldInfo::Type::scalar, variation0[0].size()));
-  vtkWriter.addVertexData(variation1Function, VTK::FieldInfo("variation 1", VTK::FieldInfo::Type::scalar, variation1[0].size()));
-
-  vtkWriter.write("sphere-patch-" + nameSuffix);
-}
-
-int main(int argc, char* argv[]) try
-{
-  static const int dim = 2;
-  typedef UnitVector<double,3> TargetSpace;
-
-  //////////////////////////////////////////////////////////////
-  //   Set up a first-order interpolation function on the sphere
-  //////////////////////////////////////////////////////////////
-
-  std::vector<TargetSpace> coefficients(3);
-
-  coefficients[0] = TargetSpace(FieldVector<double,3>({1,0,0}));
-  coefficients[1] = TargetSpace(FieldVector<double,3>({0,1,0}));
-  coefficients[2] = TargetSpace(FieldVector<double,3>({0,0,1}));
-
-  typedef LocalGeodesicFEFunction<dim, double, P1LocalFiniteElement<double,double,dim>, TargetSpace> P1LocalGFEFunctionType;
-  //typedef GFE::LocalProjectedFEFunction<dim, double, LocalFiniteElement, TargetSpace> LocalProjectedFEFunctionType;
-  //typedef GFE::LocalTangentFEFunction<dim, double, LocalFiniteElement, TargetSpace> LocalTangentFEFunctionType;
-
-  P1LocalFiniteElement<double,double,dim> localFiniteElement;
-  P1LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,coefficients);
-  interpolate<TargetSpace, P1LocalGFEFunctionType>(localGeodesicFEFunction, 0, "riemannian-p1");
-
-  //interpolate<TargetSpace,LocalProjectedFEFunctionType>(coefficients, "projected");
-  //interpolate<TargetSpace,LocalTangentFEFunctionType>(coefficients, "tangent");
-
-  //////////////////////////////////////////////////////////////
-  //   Set up a second-order interpolation function on the sphere
-  //////////////////////////////////////////////////////////////
-
-  coefficients.resize(6);
-
-  coefficients[0] = TargetSpace(FieldVector<double,3>({1,0,0}));
-  coefficients[1] = TargetSpace(FieldVector<double,3>({0.5,0.5,0.15}));
-  coefficients[2] = TargetSpace(FieldVector<double,3>({0,1,0}));
-  coefficients[3] = TargetSpace(FieldVector<double,3>({0.5,0.15,0.5}));
-  coefficients[4] = TargetSpace(FieldVector<double,3>({0.15,0.5,0.5}));
-  coefficients[5] = TargetSpace(FieldVector<double,3>({0,0,1}));
-
-  typedef LocalGeodesicFEFunction<dim, double, P2LocalFiniteElement<double,double,dim>, TargetSpace> P2LocalGFEFunctionType;
-  //typedef GFE::LocalProjectedFEFunction<dim, double, LocalFiniteElement, TargetSpace> LocalProjectedFEFunctionType;
-  //typedef GFE::LocalTangentFEFunction<dim, double, LocalFiniteElement, TargetSpace> LocalTangentFEFunctionType;
-
-  P2LocalFiniteElement<double,double,dim> p2LocalFiniteElement;
-  P2LocalGFEFunctionType p2LocalGeodesicFEFunction(p2LocalFiniteElement,coefficients);
-  interpolate<TargetSpace, P2LocalGFEFunctionType>(p2LocalGeodesicFEFunction, 0, "riemannian-p2-vertex");
-  interpolate<TargetSpace, P2LocalGFEFunctionType>(p2LocalGeodesicFEFunction, 1, "riemannian-p2-edge");
-
-  //interpolate<TargetSpace,LocalProjectedFEFunctionType>(coefficients, "projected");
-  //interpolate<TargetSpace,LocalTangentFEFunctionType>(coefficients, "tangent");
-
-} catch (Exception& e) {
-
-    std::cout << e.what() << std::endl;
-
-}