diff --git a/dune/phasefield/archive/acfs_ibvp.hh b/dune/phasefield/archive/acfs_ibvp.hh
deleted file mode 100644
index 7b61d9d6892c22d9d1b5eb213aa6c4249f8ce225..0000000000000000000000000000000000000000
--- a/dune/phasefield/archive/acfs_ibvp.hh
+++ /dev/null
@@ -1,280 +0,0 @@
-
-#pragma once
-//#include "../myincludes.hh"
-#include <cmath>
-
-template<typename GV, typename RF, std::size_t dim_range>
-class ZeroVectorFunction :
-  public Dune::PDELab::AnalyticGridFunctionBase<
-  Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim_range>,
-  ZeroVectorFunction<GV,RF,dim_range> >,
-  public Dune::PDELab::InstationaryFunctionDefaults
-{
-public:
-  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim_range> Traits;
-  typedef Dune::PDELab::AnalyticGridFunctionBase<Traits, ZeroVectorFunction> BaseT;
-
-  typedef typename Traits::DomainType DomainType;
-  typedef typename Traits::RangeType RangeType;
-
-  ZeroVectorFunction(const GV & gv) : BaseT(gv) {}
-
-  inline void evaluateGlobal(const DomainType & x, RangeType & y) const
-  {
-    y=0;
-  }
-};
-
-template<typename GV, typename RF>
-class ZeroScalarFunction
-  : public ZeroVectorFunction<GV,RF,1>
-{
-public:
-
-  ZeroScalarFunction(const GV & gv) : ZeroVectorFunction<GV,RF,1>(gv) {}
-
-};
-
-
-
-template<typename GV, typename RF, int dim>
-class PhiZero :
-  public Dune::PDELab::AnalyticGridFunctionBase<Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim>, PhiZero<GV,RF,dim> >
-{
-  const RF mean;
-  const RF maxnoise;
-  const RF eps;
-
-public:
-  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim> Traits;
-  typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiZero<GV,RF,dim> > BaseT;
-  typedef typename Traits::DomainType DomainType;
-  typedef typename Traits::RangeType RangeType;
-
-  PhiZero(const GV & gv, const RF& mean_, const RF& maxnoise_, const RF& eps_) :
-     BaseT(gv), mean(mean_), maxnoise(maxnoise_), eps(eps_) {}
-
-  inline void evaluateGlobal(const DomainType & x, RangeType & y) const
-  {
-    y = 0.0;
-    //Random Noise
-    //y[0] = mean+ maxnoise*(-1 + 2*((RF)rand())/RAND_MAX);
-    //Deterministic "pseudeorandom"
-    //y[0] = mean+ maxnoise* sin(sqrt(3)*100000*x[0]+sqrt(2))*cos(sqrt(2)*1000000*x[1]+sqrt(5))
-
-    //Bubble
-    //RF radius = sqrt( (x[0]-0.5)*(x[0]-0.5) + (x[1]-0.5)*(x[1]-0.5) );
-    //y[0] = 0.5*(1+tanh(1/eps/sqrt(2)*(radius-0.2)));
-
-    //Stripes
-    // RF d = 0;
-    // if(x[0]<3./8)
-    //   d = x[0]-1./4;
-    // else if(x[0] < 5./8)
-    //   d = 1./8 - (x[0]-3./8);
-    // else
-    //   d = -1./8 + (x[0]-5./8);
-    // y[0] = 0.5*(1+tanh(1/eps/sqrt(2)*d));
-
-    //Left/Right
-    //y[0] = 0.5*(1+tanh(1/eps/sqrt(2)*(x[0]-0.75)));
-
-    //Top/Bottom (Rayleigh-Taylor Instability)
-    y[0] = 0.5*(1-tanh(1/eps/sqrt(2)*(x[1]-0.75-0.01*cos(x[0]*2*3.1415))));
-  }
-};
-
-template<typename GV, typename RF, int dim>
-class ThreephasePhiZero :
-  public Dune::PDELab::AnalyticGridFunctionBase<Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim>, ThreephasePhiZero<GV,RF,dim> >
-{
-  const RF mean;
-  const RF maxnoise;
-  const RF eps;
-  const int phaseNumber;
-
-public:
-  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim> Traits;
-  typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,ThreephasePhiZero<GV,RF,dim> > BaseT;
-  typedef typename Traits::DomainType DomainType;
-  typedef typename Traits::RangeType RangeType;
-
-  ThreephasePhiZero(const GV & gv, const int& num_, const RF& mean_, const RF& maxnoise_, const RF& eps_) :
-     BaseT(gv), mean(mean_), maxnoise(maxnoise_), eps(eps_), phaseNumber(num_) {}
-
-  inline void evaluateGlobal(const DomainType & x, RangeType & y) const
-  {
-    y = 0.0;
-
-    //Stripes
-    // if(phaseNumber == 0)
-    // {
-    //   RF d = (x[0]-3./4);
-    //   y[0] = 0.5*(1+tanh(1/eps/sqrt(2)*d));
-    // }
-    // if(phaseNumber == 1)
-    // {
-    //   RF d = -(x[0]-1./4);
-    //   y[0] = 0.5*(1+tanh(1/eps/sqrt(2)*d));
-    // }
-
-    //Contact point
-    RF d1 = x[0]-0.5;
-    RF yx = 0.5*(1+tanh(1/eps/sqrt(2)*d1));
-
-    RF d2 = x[1]-0.75;
-    if(phaseNumber == 0)
-      d2 = -d2;
-    y[0]=0.5*(1+tanh(1/eps/sqrt(2)*d2))*yx;
-
-  }
-};
-
-template<typename GV, typename RF, int dim>
-class UZero :
-  public Dune::PDELab::AnalyticGridFunctionBase<Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim>, UZero<GV,RF,dim> >
-{
-  const RF mean;
-  const RF maxnoise;
-
-public:
-  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim> Traits;
-  typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,UZero<GV,RF,dim> > BaseT;
-  typedef typename Traits::DomainType DomainType;
-  typedef typename Traits::RangeType RangeType;
-
-  UZero(const GV & gv, const RF& mean_, const RF& maxnoise_) :
-     BaseT(gv), mean(mean_), maxnoise(maxnoise_) {}
-
-  inline void evaluateGlobal(const DomainType & x, RangeType & y) const
-  {
-    y[0] = mean;
-    //if (x[0]> 1-1e-5)
-    //{
-    //  y[0]=mean-0.3;
-    //}
-  }
-};
-
-
-
-class U_BC : public Dune::PDELab::DirichletConstraintsParameters
-{
-public:
-  template<typename I>
-  bool isDirichlet(const I & intersection, const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
-  {
-    Dune::FieldVector<typename I::ctype, I::coorddimension>
-      xg = intersection.geometry().global( coord );
-    if (xg[0]> 1-1e-5)
-    {
-      return true;
-    }
-    return false;
-  }
-};
-
-
-class SinglePointPressureBC : public Dune::PDELab::DirichletConstraintsParameters
-{
-public:
-  template<typename I> bool isDirichlet(const I & intersection,
-                                        const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
-  {
-    Dune::FieldVector<typename I::ctype, I::coorddimension>
-      xg = intersection.geometry().global( coord );
-    //std::cout<< "x0: " << xg[0] << "    x1: " << xg[1] << std::endl;
-    if (xg[0]< 1e-2)
-    {
-      //std::cout<< "Found Pressure Zero Point";
-      //return true;
-    }
-    return false;
-  }
-};
-
-template<typename GV, typename RF, int dim>
-class IniVelocity :
-  public Dune::PDELab::AnalyticGridFunctionBase<
-  Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim>,
-  IniVelocity<GV,RF,dim> >
-{
-  const RF meanflow_;
-  const RF timemax;
-  RF time;
-
-
-public:
-  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim> Traits;
-  typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,IniVelocity<GV,RF,dim> > BaseT;
-
-  typedef typename Traits::DomainType DomainType;
-  typedef typename Traits::RangeType RangeType;
-
-  IniVelocity(const GV & gv, const RF& meanflow, const RF& tm) :
-    BaseT(gv)
-    , meanflow_(meanflow), timemax(tm)
-  { time=0.0;}
-
-  template <typename T>
-  void setTime(T t){
-    time = t;
-  }
-
-  inline void evaluateGlobal(const DomainType & x, RangeType & y) const
-  {
-    y = 0.0;
-    RF radius = sqrt( (x[0]-0.5)*(x[0]-0.5) + (x[1]-0.5)*(x[1]-0.5) );
-    RF r = meanflow_*(0.5-radius);
-    if(r< 0)
-    {
-      r=0;
-    }
-
-    //y[0] = -r*(x[1]-0.5) * 1e-4;     //Rotating flow
-    //y[1] = r*(x[0]-0.5)  * 1e-4;
-
-    //y[1]=x[0];              //Shear flow
-    //y[0]=0;
-
-    //y[0]=0;                   //No flow
-    //y[1]=0;
-    RF tscale=fmin(time/timemax,1);
-
-    y[0]=0;                   //Shear flow boundary
-    if(x[0]>0.9999)
-    {
-      //y[1]=4*meanflow_*x[1]*(1-x[1])*tscale;
-      if (x[1]>0 && x[1]<0.75)
-      {
-        y[1]=64/9*meanflow_*x[1]*(0.75-x[1])*tscale;
-      }
-    }
-  }
-};
-
-
-template<typename GV, typename RF, int dim>
-class TUPressure :
-  public Dune::PDELab::AnalyticGridFunctionBase<
-  Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim>,
-  TUPressure<GV,RF,dim> >
-{
-public:
-  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim> Traits;
-  typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,TUPressure<GV,RF,dim> > BaseT;
-
-  typedef typename Traits::DomainType DomainType;
-  typedef typename Traits::RangeType RangeType;
-
-  TUPressure(const GV & gv) :
-    BaseT(gv)
-  { }
-
-  inline void evaluateGlobal(const DomainType & x, RangeType & y) const
-  {
-    y = 0.0;
-    //y[0]=0.1*x[0];
-  }
-  template <typename T>  void setTime(T t){}
-};
diff --git a/dune/phasefield/archive/acfs_problem.hh b/dune/phasefield/archive/acfs_problem.hh
deleted file mode 100644
index 1d6a836b29fc459926155efde62931074af8296d..0000000000000000000000000000000000000000
--- a/dune/phasefield/archive/acfs_problem.hh
+++ /dev/null
@@ -1,160 +0,0 @@
-
-#pragma once
-//#include "../myincludes.hh"
-#include <cmath>
-
-template<typename Number>
-class Problem
-{
-  Number t;
-
-  const Number TurningPointLeft  = 0.5 - 1/sqrt(12);
-  const Number TurningPointRight = 0.5 + 1/sqrt(12);
-public:
-  typedef Number value_type;
-  Number eps, delta;
-  Number D;
-  Number dt;
-  Number uAst;
-  Number u_regularization;
-  Number mu;
-
-  Number rho_0;
-
-  Number rho_1;
-  Number rho_2;
-
-  Number M;
-  Number M_reduced;
-  Number fscaling;
-  Number SurfaceTensionScaling;
-
-  //Number rho (const Number& x) {return rho_0;}
-  Number rho (const Number& x){return x*rho_2 + (1-x)*rho_1;}
-
-  Number Sigma1,Sigma2,Sigma3,SigmaT;
-
-
-  // Constructor takes parameter
-  Problem (const Number& dt_, const Dune::ParameterTree& ptree)
-           : t(0.0)
-  {
-    eps = ptree.get("problem.eps",(Number)0.1);
-    delta = ptree.get("problem.delta",(Number)0.1);
-    M = ptree.get("problem.Mobility",(Number)0.1);
-    M_reduced = ptree.get("problem.MobilityReduced",(Number)0.1);
-    u_regularization = ptree.get("problem.URegularization",(Number)0.01);
-
-    D = ptree.get("physics.D",(Number)0.1);
-    uAst = ptree.get("physics.uAst",(Number)0.1);
-    rho_0 = ptree.get("physics.rho",(Number)1);
-    mu = ptree.get("physics.mu",(Number)0.001);
-    fscaling = ptree.get("phsyics.fscaling",(Number)1);
-    SurfaceTensionScaling = ptree.get("physics.SurfaceTensionScaling",(Number)1);
-    rho_1 = ptree.get("3p.rho1",(Number)1);
-    rho_2 = ptree.get("3p.rho2",(Number)1);
-
-    Sigma1 = ptree.get("3p.Sigma1",(Number)1);
-    Sigma2 = ptree.get("3p.Sigma2",(Number)1);
-    Sigma3 = ptree.get("3p.Sigma3",(Number)1);
-    SigmaT = 1/(1/Sigma1 + 1/Sigma2 + 1/Sigma3);
-
-    dt=dt_;
-  }
-
-  Number Mob (const Number& u) const
-  {
-    return M; //M_reduced + (M-M_reduced)*q(u);
-  }
-
-  //Double Well
-  Number W (const Number& u) const
-  {
-    return u*u*(1-u)*(1-u);
-  }
-  // Double Well Derivative(/Secant)
-  Number dW (const Number& x, const Number& y) const
-  {
-    return (x+y)*(1 +x*x +y*y) -2*(x*x +x*y +y*y);
-  }
-
-  // Double Well Derivative
-  Number dW (const Number& x) const
-  {
-    return 2*x*(1-x)*(1-x) - 2*x*x*(1-x);
-  }
-
-  // Double Well Derivative of Convex Part
-  Number dWConvex (const Number& x) const
-  {
-    if(x<TurningPointLeft)
-      return dW(x)-dW(TurningPointLeft);
-    if(x>TurningPointRight)
-      return dW(x)-dW(TurningPointRight);
-    return 0;
-  }
-
-  // Double Well Derivative of Convcarve Part
-  Number dWConcarve (const Number& x) const
-  {
-    if(x<TurningPointLeft)
-      return dW(TurningPointLeft);
-    if(x>TurningPointRight)
-      return dW(TurningPointRight);
-    return dW(x);
-  }
-
-  Number dWdpf1 (const Number& pf1, const Number& pf2, const Number& pf3) const
-  {
-    return (Sigma1-SigmaT)*dW(pf1)-SigmaT*dW(pf2)-SigmaT*dW(pf3);
-  }
-
-  Number dWdpf2 (const Number& pf1, const Number& pf2, const Number& pf3) const
-  {
-    return (Sigma2-SigmaT)*dW(pf2)-SigmaT*dW(pf1)-SigmaT*dW(pf3);
-  }
-
-  Number dWdpf3 (const Number& pf1, const Number& pf2, const Number& pf3) const
-  {
-    return (Sigma3-SigmaT)*dW(pf3)-SigmaT*dW(pf1)-SigmaT*dW(pf2);
-  }
-
-  Number f (const Number& u) const
-  {
-    if(u<0)
-      return fscaling*(-1.0+1000*u);
-    else
-      return fscaling*(u*u-1);
-  }
-
-  Number q (const Number& pf) const
-  {
-    if (pf<0 || pf>1)
-      return 0;
-    return pf*(1-pf);
-  }
-
-  Number q (const Number& pf1, const Number& pf2, const Number& pf3) const
-  {
-    if (pf1<0 || pf3<0)
-      return 0;
-    return pf1*pf3;
-  }
-
-  Number Diffusion(const Number& pf_c)
-  {
-    return D*fmax(1e-3,(pf_c-0.05)/0.95);
-  }
-
-  Number VelDissipation(const Number& pf_f)
-  {
-    return 10*(1-pf_f)*(1-pf_f);
-  }
-
-  //! Set time in instationary case
-  void setTime (Number t_)
-  {
-    t = t_;
-  }
-
-};
diff --git a/dune/phasefield/archive/ff_parameters.hh b/dune/phasefield/archive/ff_parameters.hh
deleted file mode 100644
index 1d6a836b29fc459926155efde62931074af8296d..0000000000000000000000000000000000000000
--- a/dune/phasefield/archive/ff_parameters.hh
+++ /dev/null
@@ -1,160 +0,0 @@
-
-#pragma once
-//#include "../myincludes.hh"
-#include <cmath>
-
-template<typename Number>
-class Problem
-{
-  Number t;
-
-  const Number TurningPointLeft  = 0.5 - 1/sqrt(12);
-  const Number TurningPointRight = 0.5 + 1/sqrt(12);
-public:
-  typedef Number value_type;
-  Number eps, delta;
-  Number D;
-  Number dt;
-  Number uAst;
-  Number u_regularization;
-  Number mu;
-
-  Number rho_0;
-
-  Number rho_1;
-  Number rho_2;
-
-  Number M;
-  Number M_reduced;
-  Number fscaling;
-  Number SurfaceTensionScaling;
-
-  //Number rho (const Number& x) {return rho_0;}
-  Number rho (const Number& x){return x*rho_2 + (1-x)*rho_1;}
-
-  Number Sigma1,Sigma2,Sigma3,SigmaT;
-
-
-  // Constructor takes parameter
-  Problem (const Number& dt_, const Dune::ParameterTree& ptree)
-           : t(0.0)
-  {
-    eps = ptree.get("problem.eps",(Number)0.1);
-    delta = ptree.get("problem.delta",(Number)0.1);
-    M = ptree.get("problem.Mobility",(Number)0.1);
-    M_reduced = ptree.get("problem.MobilityReduced",(Number)0.1);
-    u_regularization = ptree.get("problem.URegularization",(Number)0.01);
-
-    D = ptree.get("physics.D",(Number)0.1);
-    uAst = ptree.get("physics.uAst",(Number)0.1);
-    rho_0 = ptree.get("physics.rho",(Number)1);
-    mu = ptree.get("physics.mu",(Number)0.001);
-    fscaling = ptree.get("phsyics.fscaling",(Number)1);
-    SurfaceTensionScaling = ptree.get("physics.SurfaceTensionScaling",(Number)1);
-    rho_1 = ptree.get("3p.rho1",(Number)1);
-    rho_2 = ptree.get("3p.rho2",(Number)1);
-
-    Sigma1 = ptree.get("3p.Sigma1",(Number)1);
-    Sigma2 = ptree.get("3p.Sigma2",(Number)1);
-    Sigma3 = ptree.get("3p.Sigma3",(Number)1);
-    SigmaT = 1/(1/Sigma1 + 1/Sigma2 + 1/Sigma3);
-
-    dt=dt_;
-  }
-
-  Number Mob (const Number& u) const
-  {
-    return M; //M_reduced + (M-M_reduced)*q(u);
-  }
-
-  //Double Well
-  Number W (const Number& u) const
-  {
-    return u*u*(1-u)*(1-u);
-  }
-  // Double Well Derivative(/Secant)
-  Number dW (const Number& x, const Number& y) const
-  {
-    return (x+y)*(1 +x*x +y*y) -2*(x*x +x*y +y*y);
-  }
-
-  // Double Well Derivative
-  Number dW (const Number& x) const
-  {
-    return 2*x*(1-x)*(1-x) - 2*x*x*(1-x);
-  }
-
-  // Double Well Derivative of Convex Part
-  Number dWConvex (const Number& x) const
-  {
-    if(x<TurningPointLeft)
-      return dW(x)-dW(TurningPointLeft);
-    if(x>TurningPointRight)
-      return dW(x)-dW(TurningPointRight);
-    return 0;
-  }
-
-  // Double Well Derivative of Convcarve Part
-  Number dWConcarve (const Number& x) const
-  {
-    if(x<TurningPointLeft)
-      return dW(TurningPointLeft);
-    if(x>TurningPointRight)
-      return dW(TurningPointRight);
-    return dW(x);
-  }
-
-  Number dWdpf1 (const Number& pf1, const Number& pf2, const Number& pf3) const
-  {
-    return (Sigma1-SigmaT)*dW(pf1)-SigmaT*dW(pf2)-SigmaT*dW(pf3);
-  }
-
-  Number dWdpf2 (const Number& pf1, const Number& pf2, const Number& pf3) const
-  {
-    return (Sigma2-SigmaT)*dW(pf2)-SigmaT*dW(pf1)-SigmaT*dW(pf3);
-  }
-
-  Number dWdpf3 (const Number& pf1, const Number& pf2, const Number& pf3) const
-  {
-    return (Sigma3-SigmaT)*dW(pf3)-SigmaT*dW(pf1)-SigmaT*dW(pf2);
-  }
-
-  Number f (const Number& u) const
-  {
-    if(u<0)
-      return fscaling*(-1.0+1000*u);
-    else
-      return fscaling*(u*u-1);
-  }
-
-  Number q (const Number& pf) const
-  {
-    if (pf<0 || pf>1)
-      return 0;
-    return pf*(1-pf);
-  }
-
-  Number q (const Number& pf1, const Number& pf2, const Number& pf3) const
-  {
-    if (pf1<0 || pf3<0)
-      return 0;
-    return pf1*pf3;
-  }
-
-  Number Diffusion(const Number& pf_c)
-  {
-    return D*fmax(1e-3,(pf_c-0.05)/0.95);
-  }
-
-  Number VelDissipation(const Number& pf_f)
-  {
-    return 10*(1-pf_f)*(1-pf_f);
-  }
-
-  //! Set time in instationary case
-  void setTime (Number t_)
-  {
-    t = t_;
-  }
-
-};
diff --git a/dune/phasefield/archive/tmp___myincludes.hh b/dune/phasefield/archive/tmp___myincludes.hh
deleted file mode 100644
index 86fe4121616bbfa13ebf9ecde4742476c0dfff24..0000000000000000000000000000000000000000
--- a/dune/phasefield/archive/tmp___myincludes.hh
+++ /dev/null
@@ -1,59 +0,0 @@
-#pragma once
-
-
-// C includes
-#include<sys/stat.h>
-// C++ includes
-#include<math.h>
-#include<iostream>
-// dune-common includes
-#include<dune/common/parallel/mpihelper.hh>
-#include<dune/common/parametertreeparser.hh>
-#include<dune/common/timer.hh>
-// dune-geometry includes
-#include<dune/geometry/referenceelements.hh>
-#include<dune/geometry/quadraturerules.hh>
-// dune-grid includes
-#include<dune/grid/onedgrid.hh>
-#include<dune/grid/yaspgrid.hh>
-#include<dune/grid/utility/structuredgridfactory.hh>
-#include<dune/grid/io/file/vtk.hh>
-#include<dune/grid/io/file/gmshreader.hh>
-#if HAVE_UG
-#include<dune/grid/uggrid.hh>
-#endif
-#if HAVE_DUNE_ALUGRID
-#include<dune/alugrid/grid.hh>
-#include<dune/alugrid/dgf.hh>
-#include<dune/grid/io/file/dgfparser/dgfparser.hh>
-#endif
-// dune-istl included by pdelab
-// dune-pdelab includes
-#include<dune/pdelab/common/function.hh>
-#include<dune/pdelab/common/vtkexport.hh>
-#include<dune/pdelab/finiteelementmap/pkfem.hh>
-#include<dune/pdelab/finiteelementmap/qkfem.hh>
-#include<dune/pdelab/finiteelementmap/p0fem.hh>
-#include<dune/pdelab/constraints/common/constraints.hh>
-#include<dune/pdelab/constraints/common/constraintsparameters.hh>
-#include<dune/pdelab/constraints/conforming.hh>
-#include<dune/pdelab/function/callableadapter.hh>
-#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
-#include<dune/pdelab/gridfunctionspace/subspace.hh>
-#include<dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
-#include<dune/pdelab/gridfunctionspace/interpolate.hh>
-#include<dune/pdelab/gridfunctionspace/vtk.hh>
-#include<dune/pdelab/gridoperator/gridoperator.hh>
-#include<dune/pdelab/gridoperator/onestep.hh>
-#include<dune/pdelab/localoperator/defaultimp.hh>
-#include<dune/pdelab/localoperator/pattern.hh>
-#include<dune/pdelab/localoperator/flags.hh>
-#include<dune/pdelab/localoperator/variablefactories.hh>
-#include<dune/pdelab/backend/istl.hh>
-#include<dune/pdelab/stationary/linearproblem.hh>
-#include<dune/pdelab/instationary/onestep.hh>
-#include<dune/pdelab/newton/newton.hh>
-#include<dune/pdelab/adaptivity/adaptivity.hh>
-
-//#include <dune/xt/grid/grids.hh>
-//#include <dune/xt/grid/view/periodic.hh>