diff --git a/dune/phasefield/base_parameters.hh b/dune/phasefield/base_parameters.hh
index e35943de17b28b1dffe3f0cf727f0c5dfaf20f59..940f3c86bd2a8a348e7288307f3e305a85f09b75 100644
--- a/dune/phasefield/base_parameters.hh
+++ b/dune/phasefield/base_parameters.hh
@@ -78,6 +78,51 @@ public:
   }
 };
 
+template<typename Number>
+class Reaction_Calcite_Precipitation // Calcite precipitation CO3^2- + Ca^2+ -> CaCO3
+                                     // unit of concentration in mol[CO3] /l (equilibrium at around 0.14*10^-3 mol/l)
+                                     // we only know the reaction is happening faster than the timescale of seconds
+                                     // As we assume the concentration of Ca^2+ to be big (around 0.16 mol/l)
+                                     // we get a linear reaction kinetic
+{
+public:
+  const Number reactionRate;
+  const Number equilibriumConcentration;
+  const Number bulkRate;
+  Reaction_Calcite_Precipitation(const Dune::ParameterTree & param) :
+    reactionRate(param.get("ReactionRate",(Number)10)),
+    equilibriumConcentration(param.get("uEquilibrium",(Number)0.00014)),
+    bulkRate(param.get("BulkRate",(Number)0.00001))
+  {
+  }
+
+  Number eval( const Number& u) const
+  {
+    return reactionRate*(u-equilibriumConcentration);
+  }
+
+  Number eval_bulkreaction() const
+  {
+    return bulkRate;
+  }
+};
+
+template<typename Number>
+class Reaction_None
+{
+public:
+  Reaction_None(const Dune::ParameterTree & param)
+  {
+  }
+
+  Number eval( const Number& u) const
+  {
+    return 0;
+  }
+
+  Number eval_bulkreaction() const { return 0;}
+};
+
 //______________________________________________________________________________________________________
 
 template<typename Number>
@@ -240,6 +285,19 @@ public:
   {
   }
 };
+
+template <typename Number, typename DW, typename Mobility, typename Reaction, typename VelDissipation>
+class Params_summerschool : public Params_fs_r<Number, DW, Mobility, Reaction, VelDissipation>
+{
+public:
+  Number CellHeight;
+
+  Params_summerschool(const Dune::ParameterTree& ptree) :
+    Params_fs_r<Number, DW, Mobility, Reaction, VelDissipation>(ptree),
+    CellHeight(ptree.get("domain.CellHeight",0.1))
+  {
+  }
+};
 //
 // template<typename Number>
 // class Problem
diff --git a/dune/phasefield/boundary_fn.hh b/dune/phasefield/boundary_fn.hh
index becbe5700cac6dd7195195038018a4e902d28463..9ad4ac4d7dea1a0dfcfad4aee92ba3fb69a7a2e5 100644
--- a/dune/phasefield/boundary_fn.hh
+++ b/dune/phasefield/boundary_fn.hh
@@ -5,6 +5,10 @@
 #include <dune/pdelab/constraints/common/constraintsparameters.hh>
 
 #define GRID_EPS 1e-7
+#define GRID_WIDTH 6      // 1.5 for summer school //2 for contactpointtest //6 for full summer school
+
+#define GRID_BOTTOM 0   //0 For summer school   //0.35 for contactpointtest
+#define GRID_TOP 0.125  //0.125 for summer school      //1 for contactpointtest
 
 class InflowOutflow_v : public Dune::PDELab::DirichletConstraintsParameters
 {
@@ -14,7 +18,7 @@ public:
   {
     Dune::FieldVector<typename I::ctype, I::coorddimension>
       xg = intersection.geometry().global( coord );
-    if (xg[0]> 1-GRID_EPS && xg[1]>-1+GRID_EPS && xg[1]<1-GRID_EPS)
+    if ( (xg[0]> GRID_WIDTH-GRID_EPS  /*|| xg[0] < GRID_EPS*/) && xg[1]>GRID_BOTTOM+GRID_EPS && xg[1]<GRID_TOP-GRID_EPS)
     {
       return false;
     }
@@ -32,7 +36,7 @@ public:
   {
     Dune::FieldVector<typename I::ctype, I::coorddimension>
       xg = intersection.geometry().global( coord );
-    if (xg[0]> 1-GRID_EPS && xg[1]>-1+GRID_EPS && xg[1]<1-GRID_EPS)
+    if ((xg[0]> GRID_WIDTH-GRID_EPS /*|| xg[0] < GRID_EPS*/) && xg[1]>GRID_BOTTOM+GRID_EPS && xg[1]<GRID_TOP-GRID_EPS)
     {
       return true;
     }
@@ -64,7 +68,7 @@ public:
   {
     Dune::FieldVector<typename I::ctype, I::coorddimension>
       xg = intersection.geometry().global( coord );
-    if (xg[0]< GRID_EPS || xg[0] >1-GRID_EPS)
+    if ( (xg[0]< GRID_EPS || xg[0] >GRID_WIDTH-GRID_EPS) && xg[1]>GRID_BOTTOM+GRID_EPS && xg[1]<GRID_TOP-GRID_EPS)
     {
       return true;
     }
diff --git a/dune/phasefield/chns_base.hh b/dune/phasefield/chns_base.hh
index 38d762b67d704b9a62a8d0f9b0e037f78c2473f0..cf3c9b02e98cba20f66b77274c6aaa5c88e61bfa 100644
--- a/dune/phasefield/chns_base.hh
+++ b/dune/phasefield/chns_base.hh
@@ -40,6 +40,7 @@ public:
 
   double lastsave;
   double saveintervall;
+  double saveintervallstart;
 
   int verb;
 
@@ -56,6 +57,7 @@ public:
     dtDecreaseFactor = param.get("dtDecreaseFactor",(double)0.5);
 
     saveintervall = param.get("saveintervall",(double)0.0);
+    saveintervallstart = param.get("saveintervallstart",(double)0.0);
 
     time.t=param.get("tIni",(double)0.0);
     lastsave = time.t - saveintervall;
@@ -99,7 +101,7 @@ public:
       if (verb > 0)   std::cout << "New Timestep "<< time.dt << std::endl;
     }
 
-    if( time.t>lastsave+saveintervall )
+    if( time.t < saveintervallstart || time.t>lastsave+saveintervall )
     {
         lastsave = time.t;
         vtkwriter.write(time.t,Dune::VTK::ascii);
diff --git a/dune/phasefield/ffs_chns_r.hh b/dune/phasefield/ffs_chns_r.hh
index 87ca6eae63a569762a14d54f7a0f64158b9c8526..f529752d90d42de4025e5175e44a924c474c105b 100644
--- a/dune/phasefield/ffs_chns_r.hh
+++ b/dune/phasefield/ffs_chns_r.hh
@@ -197,34 +197,36 @@ public:
 
 //_________________________________________________________________________________________________________________________
 
-template<const int dim, typename BdryV_Scalar, typename BdryP,
+template<typename BdryV1, typename BdryV2, typename BdryP,
           typename BdryPhi, typename BdryMu, typename BdryPhi2, typename BdryMu2, typename BdryU>
-class Bdry_ffs_chns_r
+class Bdry_ffs_chns_r_compositeV_2d
 {
 public:
-  BdryV_Scalar bdryV_Scalar;
+  BdryV1 bdryV1;
+  BdryV2 bdryV2;
   BdryP bdryP;
   BdryPhi bdryPhi;
   BdryMu bdryMu;
   BdryPhi2 bdryPhi2;
   BdryMu2 bdryMu2;
   BdryU bdryU;
-  typedef Dune::PDELab::PowerConstraintsParameters<BdryV_Scalar,dim>    BdryV;
+  typedef Dune::PDELab::CompositeConstraintsParameters<BdryV1,BdryV2>    BdryV;
   typedef Dune::PDELab::CompositeConstraintsParameters<BdryV, BdryP> BdryNS;
   typedef Dune::PDELab::CompositeConstraintsParameters<BdryPhi, BdryMu, BdryPhi2, BdryMu2, BdryU> BdryCH;
   BdryV bdryV;
   BdryNS bdryNS;
   BdryCH bdryCH;
 
-  Bdry_ffs_chns_r() :
-      bdryV_Scalar(),
+  Bdry_ffs_chns_r_compositeV_2d() :
+      bdryV1(),
+      bdryV2(),
       bdryP(),
       bdryPhi(),
       bdryMu(),
       bdryPhi2(),
       bdryMu2(),
       bdryU(),
-      bdryV(bdryV_Scalar),
+      bdryV(bdryV1,bdryV2),
       bdryNS(bdryV,bdryP),
       bdryCH(bdryPhi,bdryMu,bdryPhi2,bdryMu2,bdryU)
   {
diff --git a/dune/phasefield/initial_fn.hh b/dune/phasefield/initial_fn.hh
index a9522f69f7b841857d95e4fee56890023c3ea0e7..835b6b5411dca0c1bc29d53868f96983b1fe82da 100644
--- a/dune/phasefield/initial_fn.hh
+++ b/dune/phasefield/initial_fn.hh
@@ -104,7 +104,6 @@ public:
 
   inline void evaluateGlobal(const  typename Traits::DomainType & x, typename Traits::RangeType & y) const
   {
-    y = 0.5*(1-tanh(1/eps/sqrt(2)*(x[1]-0.75-0.01*cos(x[0]*2*3.1415))));
     y = 0.0;
 
     RF d1 = x[0]-0.5;
@@ -119,6 +118,38 @@ public:
 
 //______________________________________________________________________________________________________________________
 
+template<typename GV, typename RF, typename Params, const int Phase>
+class PhiInitial_3pFlowTest :
+  public Dune::PDELab::AnalyticGridFunctionBase<
+    Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_3pFlowTest<GV,RF,Params,Phase> >
+{
+  RF eps;
+
+public:
+  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
+
+  PhiInitial_3pFlowTest(const GV & gv, const Params & params) :
+     Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiInitial_3pFlowTest<GV,RF,Params,Phase> > (gv)
+     {
+       eps = params.eps;
+     }
+
+  inline void evaluateGlobal(const  typename Traits::DomainType & x, typename Traits::RangeType & y) const
+  {
+    y = 0.0;
+
+    RF d1 = x[0]-0.3;
+
+
+    RF d2 = x[1]-0.25;
+    RF yy = 0.5*(1+tanh(1/eps/sqrt(2)*d2));
+    if(Phase == 0)
+      d1 = -d1;
+    y=0.5*(1+tanh(1/eps/sqrt(2)*d1))*yy;
+  }
+};
+//______________________________________________________________________________________________________________________
+
 template<typename GV, typename RF, typename Params>
 class PhiInitial_Sphere :
   public Dune::PDELab::AnalyticGridFunctionBase<
@@ -151,6 +182,8 @@ class PhiInitial_Rectangle :
   public Dune::PDELab::AnalyticGridFunctionBase<
     Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_Rectangle<GV,RF,Params> >
 {
+  RF xcenter;
+  RF ycenter;
   RF eps;
   RF radius;
 
@@ -162,12 +195,14 @@ public:
      {
        eps = params.eps;
        radius = params.initial.ptree.get("Radius",0.2);
+       xcenter = params.initial.ptree.get("xCenter",0.5);
+       ycenter = params.initial.ptree.get("yCenter",0.5);
      }
 
   inline void evaluateGlobal(const  typename Traits::DomainType & x, typename Traits::RangeType & y) const
   {
-    RF dist1 = std::abs(x[0]-0.5);
-    RF dist2 = std::abs(x[1]-0.5);
+    RF dist1 = std::abs(x[0]-xcenter);
+    RF dist2 = std::abs(x[1]-ycenter);
     RF m1 = 0.5*(1+tanh(1/eps/sqrt(2)*(dist1-radius)));
     RF m2 = 0.5*(1+tanh(1/eps/sqrt(2)*(dist2-radius)));
     y[0] = 1-(1-m1)*(1-m2);
@@ -176,6 +211,109 @@ public:
 
 //______________________________________________________________________________________________________________________
 
+template<typename GV, typename RF, typename Params>
+class PhiInitial_2Precip :
+  public Dune::PDELab::AnalyticGridFunctionBase<
+    Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_2Precip<GV,RF,Params> >
+{
+  RF eps;
+  RF x1center;
+  RF y1center;
+  RF radius1;
+  RF x2center;
+  RF y2center;
+  RF radius2;
+  RF x3center;
+  RF y3center;
+  RF radius3;
+
+
+public:
+  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
+
+  PhiInitial_2Precip(const GV & gv, const Params & params) :
+     Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiInitial_2Precip<GV,RF,Params> > (gv)
+     {
+       eps = params.eps;
+       radius1 = params.initial.ptree.get("Radius1",0.2);
+       x1center = params.initial.ptree.get("x1Center",0.5);
+       y1center = params.initial.ptree.get("y1Center",0.5);
+       radius2 = params.initial.ptree.get("Radius2",0.2);
+       x2center = params.initial.ptree.get("x2Center",0.5);
+       y2center = params.initial.ptree.get("y2Center",0.5);
+       radius3 = params.initial.ptree.get("Radius3",0.2);
+       x3center = params.initial.ptree.get("x3Center",0.5);
+       y3center = params.initial.ptree.get("y3Center",0.5);
+     }
+
+  inline void evaluateGlobal(const  typename Traits::DomainType & x, typename Traits::RangeType & y) const
+  {
+    RF dist1 = std::sqrt( (x[0]-x1center)*(x[0]-x1center) + (x[1]-y1center)*(x[1]-y1center) );
+    RF y1 = 0.5*(1+tanh(1/eps/sqrt(2)*(dist1-radius1)));
+
+    RF dist2 = std::sqrt( (x[0]-x2center)*(x[0]-x2center) + (x[1]-y2center)*(x[1]-y2center) );
+    RF y2 = 0.5*(1+tanh(1/eps/sqrt(2)*(dist2-radius2)));
+
+    RF dist3 = std::sqrt( (x[0]-x3center)*(x[0]-x3center) + (x[1]-y3center)*(x[1]-y3center) );
+    RF y3 = 0.5*(1+tanh(1/eps/sqrt(2)*(dist3-radius3)));
+
+    y[0] = y1*y2*y3;
+  }
+};
+
+template<typename GV, typename RF, typename Params>
+class PhiInitial_SummerSchoolFull :
+  public Dune::PDELab::AnalyticGridFunctionBase<
+    Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_SummerSchoolFull<GV,RF,Params> >
+{
+  RF eps;
+  RF radius;
+  std::vector<RF> xValues;
+  std::vector<RF> yValues;
+
+public:
+  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
+
+  PhiInitial_SummerSchoolFull(const GV & gv, const Params & params) :
+     Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiInitial_SummerSchoolFull<GV,RF,Params> > (gv)
+     {
+       eps = params.eps;
+       radius = params.initial.ptree.get("Radius",0.2);
+
+       xValues.push_back(0.560); yValues.push_back(0.070);
+       xValues.push_back(0.831); yValues.push_back(0.029);
+       xValues.push_back(1.353); yValues.push_back(0.008);
+       xValues.push_back(1.298); yValues.push_back(0.105);
+       xValues.push_back(1.387); yValues.push_back(0.108);
+       xValues.push_back(1.959); yValues.push_back(0.014);
+       xValues.push_back(2.397); yValues.push_back(-0.001);
+       xValues.push_back(3.014); yValues.push_back(0.074);
+       xValues.push_back(3.703); yValues.push_back(0.244);
+       xValues.push_back(3.808); yValues.push_back(-0.015);
+       xValues.push_back(4.190); yValues.push_back(0.071);
+       xValues.push_back(4.816); yValues.push_back(0.026);
+
+     }
+
+  inline void evaluateGlobal(const  typename Traits::DomainType & x, typename Traits::RangeType & y) const
+  {
+    RF dist;
+    RF temp;
+    RF yc = 1;
+    for(int i=0; i<xValues.size(); i++)
+    {
+      dist = std::sqrt( (x[0]-xValues[i])*(x[0]-xValues[i]) + (x[1]-yValues[i])*(x[1]-yValues[i]) );
+      temp = 0.5*(1+tanh(1/eps/sqrt(2)*(dist-radius)));
+      yc = yc*temp;
+    }
+    y[0] = yc;
+  }
+};
+
+//______________________________________________________________________________________________________________________
+
+#include <dune/phasefield/boundary_fn.hh>
+
 template<typename GV, typename RF, typename Params, int dim>
 class VelocityInitial_Inflow :
   public Dune::PDELab::AnalyticGridFunctionBase<
@@ -207,9 +345,10 @@ public:
   inline void evaluateGlobal(const DomainType & x, RangeType & y) const
   {
     y = 0.0;
-    if(x[0]<1e-7) //Inflow boundary
+    if(x[0]<GRID_EPS) //Inflow boundary
     {
-      y[0] = std::max((RF)0,3/2 * meanflow * x[1]*(1-x[1]));
+      RF factor = meanflow * 6/( (GRID_TOP-GRID_BOTTOM)*(GRID_TOP-GRID_BOTTOM) );
+      y[0] = std::max((RF)0, factor * (x[1]-GRID_BOTTOM)*(GRID_TOP-x[1]));
     }
   }
 };
@@ -242,3 +381,27 @@ public:
     }
   }
 };
+
+template<typename GV, typename RF, typename Params>
+class PInitial_Inflow :
+  public Dune::PDELab::AnalyticGridFunctionBase<
+    Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PInitial_Inflow<GV,RF,Params> >
+{
+
+public:
+  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
+
+  PInitial_Inflow(const GV & gv, const Params & params) :
+     Dune::PDELab::AnalyticGridFunctionBase<Traits,PInitial_Inflow<GV,RF,Params> > (gv)
+     {
+     }
+
+  inline void evaluateGlobal(const  typename Traits::DomainType & x, typename Traits::RangeType & y) const
+  {
+    y = 0;
+    if(x[0]<1e-7)
+    {
+      y=20;
+    }
+  }
+};
diff --git a/dune/phasefield/localoperator/fem_base_chns.hh b/dune/phasefield/localoperator/fem_base_chns.hh
index 3826fbcd173f7d495b98604bb96b34da15c0766c..33ece2387c105ec27708d439664f36f156fa706b 100644
--- a/dune/phasefield/localoperator/fem_base_chns.hh
+++ b/dune/phasefield/localoperator/fem_base_chns.hh
@@ -163,6 +163,87 @@ public:
   }
 };
 
+
+//_________________________________________________________________________________________________________________________
+
+template<typename RF, const int dim, typename CHLFS>
+class CHVariables_Threephase_Extension : public CHVariables<RF,dim,CHLFS>
+{
+public:
+  const typename CHLFS::template Child<2>::Type& pf2space;
+  const typename CHLFS::template Child<3>::Type& xi2space;
+
+  RF pf2;
+  Dune::FieldVector<RF,dim> gradpf2;
+  RF xi2;
+  Dune::FieldVector<RF,dim> gradxi2;
+
+  template<typename Pos, typename Cache, typename LV, typename S>
+  CHVariables_Threephase_Extension(const Cache& cache,
+                            const Pos& position, const CHLFS& chLfs, const LV& lv, S& s) :
+    CHVariables<RF,dim,CHLFS>(cache, position, chLfs, lv, s),
+    pf2space(child(chLfs,Dune::Indices::_2)),
+    xi2space(child(chLfs,Dune::Indices::_3)),
+    pf2(0.0),
+    gradpf2(0.0),
+    xi2(0.0),
+    gradxi2(0.0)
+  {
+    for (size_t i=0; i<this->pfspace.size(); i++)
+    {
+      pf2 += lv[pf2space.localIndex(i)]*this->basis.phi[i];
+      gradpf2.axpy(lv[pf2space.localIndex(i)],this->basis.gradphi[i][0]);
+      xi2 += lv[xi2space.localIndex(i)]*this->basis.phi[i];
+      gradxi2.axpy(lv[xi2space.localIndex(i)],this->basis.gradphi[i][0]);
+      }
+  }
+};
+
+
+//_________________________________________________________________________________________________________________________
+
+template<typename RF, const int dim, typename CHLFS>
+class CHVariables_Threephase_U_PFOld_Extension : public CHVariables_Threephase_Extension<RF,dim,CHLFS>
+{
+public:
+  const typename CHLFS::template Child<4>::Type& uspace;
+
+  RF pfOld;
+  Dune::FieldVector<RF,dim> gradpfOld;
+  RF pf2Old;
+  Dune::FieldVector<RF,dim> gradpf2Old;
+  RF u;
+  Dune::FieldVector<RF,dim> gradu;
+  RF uOld;
+
+
+  template<typename Pos, typename Cache, typename LV, typename S>
+  CHVariables_Threephase_U_PFOld_Extension(const Cache& cache,
+                            const Pos& position, const CHLFS& chLfs, const LV& lv, const LV& lvOld, S& s) :
+    CHVariables_Threephase_Extension<RF,dim,CHLFS>(cache, position, chLfs, lv, s),
+    uspace(child(chLfs,Dune::Indices::_4)),
+    pfOld(0.0),
+    gradpfOld(0.0),
+    pf2Old(0.0),
+    gradpf2Old(0.0),
+    u(0.0),
+    gradu(0.0),
+    uOld(0.0)
+  {
+    for (size_t i=0; i<this->pfspace.size(); i++)
+    {
+      u += lv[uspace.localIndex(i)]*this->basis.phi[i];
+      gradu.axpy(lv[uspace.localIndex(i)],this->basis.gradphi[i][0]);
+      uOld += lvOld [uspace.localIndex(i)]*this->basis.phi[i];
+      pfOld += lvOld[this->pfspace.localIndex(i)]*this->basis.phi[i];
+      gradpfOld.axpy(lvOld[this->pfspace.localIndex(i)],this->basis.gradphi[i][0]);
+      pf2Old += lvOld[this->pf2space.localIndex(i)]*this->basis.phi[i];
+      gradpf2Old.axpy(lvOld[this->pf2space.localIndex(i)],this->basis.gradphi[i][0]);
+    }
+  }
+};
+
+
 //_________________________________________________________________________________________________________________________
 
 template<typename RF, const int dim, typename NSLFS>
diff --git a/dune/phasefield/localoperator/fem_ffs_chns_navierstokes.hh b/dune/phasefield/localoperator/fem_ffs_chns_navierstokes.hh
index f8a5ecd5d0d97b5828346a611c9f507944850da0..c2e668d2f2f863e44e49320c35e23d7d3a00877c 100644
--- a/dune/phasefield/localoperator/fem_ffs_chns_navierstokes.hh
+++ b/dune/phasefield/localoperator/fem_ffs_chns_navierstokes.hh
@@ -8,6 +8,9 @@
 #include<dune/pdelab/localoperator/variablefactories.hh>
 #include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
 
+#include<dune/phasefield/localoperator/fem_base_chns.hh>
+
+
 
 
 template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
@@ -19,36 +22,20 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
     {
     private:
       typedef typename CHContainer::field_type RF;
+      typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity   ENTITY;
 
-      typedef typename CHGFS::template Child<0>::Type FEM;
-      typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasis;
-      Dune::PDELab::LocalBasisCache<LocalBasis> cache;
+      CHLocalBasisCache<CHGFS> chCache;
+      VLocalBasisCache<NSGFS> vCache;
+      PLocalBasisCache<NSGFS> pCache;
 
-      typedef typename NSGFS::template Child<0>::Type::template Child<0>::Type VFEM;
-      typedef typename VFEM::Traits::FiniteElementType::Traits::LocalBasisType VLocalBasis;
-      Dune::PDELab::LocalBasisCache<VLocalBasis> vcache;
+      typedef GFS_LocalViewer<CHGFS, CHContainer, ENTITY, RF> CHLView;
+      mutable CHLView chLocal;
+      //mutable CHLView chOldLocal;
 
-      Param& param; // parameter functions
+      typedef GFS_LocalViewer<NSGFS, NSContainer, ENTITY, RF> NSLView;
+      mutable NSLView nsOldLocal;
 
-      //FOR EVALUATING the new and old cahn-Hilliard values
-      typedef Dune::PDELab::LocalFunctionSpace<CHGFS> CHLFS;
-      mutable CHLFS chLfs;
-      typedef Dune::PDELab::LFSIndexCache<CHLFS> CHLFSCache;
-      mutable CHLFSCache chLfsCache;
-      typedef typename CHContainer::template ConstLocalView<CHLFSCache> CHView;
-      mutable CHView chView;
-      mutable CHView chOldView;
-      mutable std::vector<typename LocalBasis::Traits::RangeFieldType> chLocal;
-      mutable std::vector<typename LocalBasis::Traits::RangeFieldType> chOldLocal;
-
-      //FOR EVALUATING the old navier-Stokes values
-      typedef Dune::PDELab::LocalFunctionSpace<NSGFS> NSLFS;
-      mutable NSLFS nsLfs;
-      typedef Dune::PDELab::LFSIndexCache<NSLFS> NSLFSCache;
-      mutable NSLFSCache nsLfsCache;
-      typedef typename NSContainer::template ConstLocalView<NSLFSCache> NSView;
-      mutable NSView nsView;
-      mutable std::vector<typename VLocalBasis::Traits::RangeFieldType> nsOldLocal;
+      Param& param; // parameter functions
 
       //FOR DETECTING ENTITY changes
       typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity   ENT;
@@ -68,10 +55,10 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
 
       FEM_FFS_CHNS_NavierStokes (Param& param_, CHGFS& chGfs_, CHContainer& chV_, CHContainer& chVOld_,
                          NSGFS& nsGfs_, NSContainer& nsVOld_) :
-        param(param_),
-        chLfs(chGfs_), chLfsCache(chLfs), chView(chV_), chOldView(chVOld_),
-        chLocal(chGfs_.maxLocalSize()), chOldLocal(chGfs_.maxLocalSize()),
-        nsLfs(nsGfs_), nsLfsCache(nsLfs), nsView(nsVOld_), nsOldLocal(nsGfs_.maxLocalSize())
+        chLocal(chGfs_,chV_),
+        //chOldLocal(chGfs_,chVOld_),
+        nsOldLocal(nsGfs_,nsVOld_),
+        param(param_)
       {
       }
 
@@ -79,132 +66,43 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
       template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
       void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
       {
-        //Left hand side, terms depending on the solution x
-        using namespace Dune::Indices;
 
         //Quadrature Rule
         const int dim = EG::Entity::dimension;
-        const int order = 2; //TODO: Choose more educated
+        const int order = 5; //TODO: Choose more educated
         auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
-        //const auto& myordering = lfsv.gridFunctionSpace().ordering();
-
-        const auto& pf1space = child(chLfs,_0);
-        const auto& xi1space = child(chLfs,_1);
-        const auto& pf2space = child(chLfs,_2);
-        const auto& xi2space = child(chLfs,_3);
-        const auto& femspace = child(chLfs,_0);
-
-        const auto& vspace = child(lfsv,_0);
-        const auto& pspace = child(lfsv,_1);
-        const auto& vfemspace = child(vspace,_0);
-
-        if(eg.entity() != OldEnt)
-        {
-          OldEnt = eg.entity();
 
-          chLfs.bind(eg.entity());
-          chLfsCache.update();
-
-          chView.bind(chLfsCache);
-          chView.read(chLocal);
-          chView.unbind();
-
-          chOldView.bind(chLfsCache);
-          chOldView.read(chOldLocal);
-          chOldView.unbind();
-
-          nsLfs.bind(eg.entity());
-          nsLfsCache.update();
-
-          nsView.bind(nsLfsCache);
-          nsView.read_sub_container(vspace,nsOldLocal);
-          nsView.unbind();
-        }
+        chLocal.DoEvaluation(eg);
+        //chOldLocal.DoEvaluation(eg);
+        nsOldLocal.DoSubEvaluation(eg,child(nsOldLocal.lfs,Dune::Indices::_0));
 
         for(const auto& ip : rule)
         {
-          // evaluate shape functions and  gradient of shape functions
-          auto& phi = cache.evaluateFunction(ip.position(), femspace.finiteElement().localBasis());
-          auto& gradphihat = cache.evaluateJacobian(ip.position(), femspace.finiteElement().localBasis());
-
-          auto& vphi = vcache.evaluateFunction(ip.position(), vfemspace.finiteElement().localBasis());
-          auto& gradvphihat = vcache.evaluateJacobian(ip.position(), vfemspace.finiteElement().localBasis());
-
-          // transform gradients of shape functions to real element
           const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
-          auto gradphi = makeJacobianContainer(femspace);
-          auto gradvphi = makeJacobianContainer(vfemspace);
-          for (size_t i=0; i<femspace.size(); i++)
-          {
-            S.mv(gradphihat[i][0],gradphi[i][0]);
-          }
-          for (size_t i=0; i<vfemspace.size(); i++)
-          {
-            S.mv(gradvphihat[i][0],gradvphi[i][0]);
-          }
-
-          //compute current values of the unknowns
-          Dune::FieldVector<RF,dim> v(0.0);
-          Dune::FieldVector<RF,dim> vOld(0.0);
-          Dune::FieldMatrix<RF,dim,dim> jacv(0.0);
-          RF p=0.0;
-          RF pf1=0.0;
-          RF pf2=0.0;
-          RF xi1=0.0;
-          RF xi2=0.0;
-          Dune::FieldVector<RF,dim> gradpf1(0.0);
-          Dune::FieldVector<RF,dim> gradpf2(0.0);
-          Dune::FieldVector<RF,dim> gradxi1(0.0);
-          Dune::FieldVector<RF,dim> gradxi2(0.0);
-
-          for (size_t i=0; i<femspace.size(); i++)
-          {
-
-            p += x(pspace,i)*phi[i];
+          auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
 
-            pf1 += chLocal[pf1space.localIndex(i)]*phi[i];
-            gradpf1.axpy(chLocal[pf1space.localIndex(i)],gradphi[i][0]);
-            pf2 += chLocal[pf2space.localIndex(i)]*phi[i];
-            gradpf2.axpy(chLocal[pf2space.localIndex(i)],gradphi[i][0]);
+          CHVariables_Threephase_Extension<RF,dim,decltype(chLocal.lfs)> ch(chCache, ip.position(), chLocal.lfs, chLocal.lv, S);
+          NSVariables_P_VOld_Extension<RF,dim,LFSV> ns(vCache, pCache, ip.position(), lfsv, x.base(), nsOldLocal.lv, S);
 
-            xi1 += chLocal[xi1space.localIndex(i)]*phi[i];
-            gradxi1.axpy(chLocal[xi1space.localIndex(i)],gradphi[i][0]);
-            xi2 += chLocal[xi2space.localIndex(i)]*phi[i];
-            gradxi2.axpy(chLocal[xi2space.localIndex(i)],gradphi[i][0]);
-          }
-          for (size_t i=0; i<vfemspace.size(); i++)
-          {
-            for(size_t k=0; k<dim; k++)
-            {
-              v[k] += x(child(vspace,k),i)*vphi[i];
-              vOld[k] += nsOldLocal[child(vspace,k).localIndex(i)]*vphi[i];
-              jacv[k].axpy(x(child(vspace,k),i),gradvphi[i][0]);
-            }
-          }
-          //std::cout << "v:" << v << std::endl;
-          //std::cout << "vOld:" << vOld << std::endl;
-          //std::cout << "jacv:" << jacv << std::endl;
-          //std::cout << "p:" << p << std::endl;
 
-          RF pf3=1-pf1-pf2;
-          RF pf_f=pf1+pf2+2*param.delta*pf3;
+          RF pf3=1-ch.pf-ch.pf2;
+          RF pf_f=ch.pf+ch.pf2+2*param.delta*pf3;
           Dune::FieldVector<RF,dim> gradpf_f(0.0);
-          gradpf_f.axpy((1-2*param.delta),gradpf1+gradpf2);
+          gradpf_f.axpy((1-2*param.delta),ch.gradpf+ch.gradpf2);
 
-          RF rho_f=param.phys.rho_f(pf1,pf2);
+          RF rho_f=param.phys.rho_f(ch.pf,ch.pf2);
           RF rho_f_reg = rho_f + (param.phys.rho_1+param.phys.rho_2)*param.delta;
 
-          auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
-          for (size_t i=0; i<femspace.size(); i++)
+          for (size_t i=0; i<ns.pspace.size(); i++)
           {
             for(int k=0; k<dim; k++)
             {
               // integrate: (div v)*phi
-              //std::cout << -1.0* jacv[k][k]*phi[i]*factor << std::endl;
-              r.accumulate(pspace,i, -1.0* pf_f*(v*gradphi[i][0])*factor);
+              //r.accumulate(ns.pspace,i, -1.0* pf_f*(v*gradphi[i][0])*factor);            Weak form
+              r.accumulate(ns.pspace,i, -1.0 * (pf_f*ns.div_v  + (gradpf_f*ns.v) )* ns.pbasis.phi[i] * factor);
             }
           }
-          for (size_t i=0; i<vfemspace.size(); i++)
+          for (size_t i=0; i<ns.vfemspace.size(); i++)
           {
             for(int k=0; k<dim; k++)
             {
@@ -215,18 +113,18 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
 
               //RF lambda = 3;
               //const auto v_nabla_v = v * jacv[k];
-              RF gradxi3k = -param.tw.Sigma3*(gradxi1[k]/param.tw.Sigma1 + gradxi2[k]/param.tw.Sigma2);
+              RF gradxi3k = -param.tw.Sigma3*(ch.gradxi[k]/param.tw.Sigma1 + ch.gradxi2[k]/param.tw.Sigma2);
               RF S = 0;
               //S += xi1*(gradpf1[k]*pf_f - pf1 * gradpf_f[k])/pf_f;
               //S += xi1*(gradpf2[k]*pf_f - pf2 * gradpf_f[k])/pf_f;
               //S += 2*param.delta*(gradxi3k-gradxi2[k]-gradxi1[k]);
-              S += pf1*gradxi1[k] + pf2*gradxi2[k];
-              S = S*param.SurfaceTensionScaling*vphi[i];
+              S += ch.pf*ch.gradxi[k] + ch.pf2*ch.gradxi2[k];
+              S = S*param.SurfaceTensionScaling*ns.vbasis.phi[i];
 
-              r.accumulate(child(vspace,k),i, (   rho_f_reg*(v[k]-vOld[k])*vphi[i]/param.time.dt          //TimeEvolution
-                                                - p * (pf_f*gradvphi[i][0][k] + vphi[i]*gradpf_f[k])            //Pressure
-                                                + param.phys.mu * (jacv[k]*gradvphi[i][0]) //+vphi[i]*(jacv[k]*gradpf)   )  //Viscosity
-                                                + param.vDis.eval(pf_f)*v[k]*vphi[i]              //Velocity Dissipation in Solid
+              r.accumulate(child(ns.vspace,k),i, (   rho_f_reg*(ns.v[k]-ns.vOld[k])*ns.vbasis.phi[i]/param.time.dt          //TimeEvolution
+                                                - ns.p * (pf_f*ns.vbasis.gradphi[i][0][k] + ns.vbasis.phi[i]*gradpf_f[k])            //Pressure
+                                                + param.phys.mu * (ns.jacv[k]*ns.vbasis.gradphi[i][0]) //+vphi[i]*(jacv[k]*gradpf)   )  //Viscosity
+                                                + param.vDis.eval(pf_f)*ns.v[k]*ns.vbasis.phi[i]              //Velocity Dissipation in Solid
                                                 //- rho*g[k]*vphi[i]                //Gravitation
                                                 + S                                //Surface Tension
                                               )*factor);
diff --git a/dune/phasefield/localoperator/fem_ffs_chns_r_cahnhilliard.hh b/dune/phasefield/localoperator/fem_ffs_chns_r_cahnhilliard.hh
index 6b4a53fe7ad6af6686280a7ff045729497507ec8..f5914b70e88deed9d8ede98ea9c9f05e0521a8bb 100644
--- a/dune/phasefield/localoperator/fem_ffs_chns_r_cahnhilliard.hh
+++ b/dune/phasefield/localoperator/fem_ffs_chns_r_cahnhilliard.hh
@@ -8,6 +8,8 @@
 #include<dune/pdelab/localoperator/variablefactories.hh>
 #include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
 
+#include<dune/phasefield/localoperator/fem_base_chns.hh>
+
 //Parameter class (for functions), FEM space and RangeField
 template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
 class FEM_FFS_CHNS_R_CahnHilliard :
@@ -18,38 +20,19 @@ class FEM_FFS_CHNS_R_CahnHilliard :
 {
 private:
   typedef typename CHContainer::field_type RF;
+  typedef typename CHGFS::Traits::GridViewType::template Codim<0>::Entity   ENTITY;
 
-  typedef typename CHGFS::template Child<0>::Type FEM;
-  typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasis;
-  Dune::PDELab::LocalBasisCache<LocalBasis> cache;
+  CHLocalBasisCache<CHGFS> chCache;
+  VLocalBasisCache<NSGFS> vCache;
 
-  typedef typename NSGFS::template Child<0>::Type::template Child<0>::Type VFEM;
-  typedef typename VFEM::Traits::FiniteElementType::Traits::LocalBasisType VLocalBasis;
-  Dune::PDELab::LocalBasisCache<VLocalBasis> vcache;
+  typedef GFS_LocalViewer<CHGFS, CHContainer, ENTITY, RF> CHLView;
+  mutable CHLView chOldLocal;
 
-  Param& param; // parameter functions
+  typedef GFS_LocalViewer<NSGFS, NSContainer, ENTITY, RF> NSLView;
+  mutable NSLView nsLocal;
 
-  //FOR EVALUATING the old cahn-Hilliard values
-  typedef Dune::PDELab::LocalFunctionSpace<CHGFS> CHLFS;
-  mutable CHLFS chLfs;
-  typedef Dune::PDELab::LFSIndexCache<CHLFS> CHLFSCache;
-  mutable CHLFSCache chLfsCache;
-  typedef typename CHContainer::template ConstLocalView<CHLFSCache> CHView;
-  mutable CHView chView;
-  mutable std::vector<typename LocalBasis::Traits::RangeFieldType> chOldLocal;
-
-  //FOR EVALUATING the old navier-Stokes values
-  typedef Dune::PDELab::LocalFunctionSpace<NSGFS> NSLFS;
-  mutable NSLFS nsLfs;
-  typedef Dune::PDELab::LFSIndexCache<NSLFS> NSLFSCache;
-  mutable NSLFSCache nsLfsCache;
-  typedef typename NSContainer::template ConstLocalView<NSLFSCache> NSView;
-  mutable NSView nsView;
-  mutable std::vector<typename VLocalBasis::Traits::RangeFieldType> nsOldLocal;
-
-  //FOR DETECTING ENTITY changes
-  typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity   ENT;
-  mutable ENT OldEnt;
+
+  Param& param; // parameter functions
 
 public:
   enum {doPatternVolume=true};
@@ -57,10 +40,10 @@ public:
   enum {doAlphaVolume=true};
 
   //constructor
-  FEM_FFS_CHNS_R_CahnHilliard(Param& param_, CHGFS& chGfs_, CHContainer& chVOld_, NSGFS& nsGfs_, NSContainer& nsVOld_) :
-    param(param_),
-    chLfs(chGfs_), chLfsCache(chLfs), chView(chVOld_), chOldLocal(chGfs_.maxLocalSize()),
-    nsLfs(nsGfs_), nsLfsCache(nsLfs), nsView(nsVOld_), nsOldLocal(nsGfs_.maxLocalSize())
+  FEM_FFS_CHNS_R_CahnHilliard(Param& param_, CHGFS& chGfs_, CHContainer& chVOld_, NSGFS& nsGfs_, NSContainer& nsV) :
+    chOldLocal(chGfs_,chVOld_),
+    nsLocal(nsGfs_,nsV),
+    param(param_)
   {
   }
 
@@ -68,156 +51,76 @@ public:
   template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
   void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
   {
-    //Left hand side, terms depending on the solution x
-    using namespace Dune::Indices;
-
     //Quadrature Rule
     const int dim = EG::Entity::dimension;
     const int order = 2; //TODO: Choose more educated
     auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
-    //const auto& myordering = lfsv.gridFunctionSpace().ordering();
-
-    const auto& pf1space = child(lfsv,_0);
-    const auto& xi1space = child(lfsv,_1);
-    const auto& pf2space = child(lfsv,_2);
-    const auto& xi2space = child(lfsv,_3);
-    const auto& uspace = child(lfsv,_4);
 
-    const auto& femspace = child(lfsv,_0);
-
-    const auto& vspace = child(nsLfs,_0);
-    const auto& vfemspace = child(vspace,_0);
-
-    if(eg.entity() != OldEnt)
-    {
-      OldEnt = eg.entity();
-
-      chLfs.bind(eg.entity());
-      chLfsCache.update();
-      chView.bind(chLfsCache);
-      chView.read(chOldLocal);
-      chView.unbind();
-
-      nsLfs.bind(eg.entity());
-      nsLfsCache.update();
-      nsView.bind(nsLfsCache);
-
-      nsView.read_sub_container(vspace,nsOldLocal);
-      nsView.unbind();
-    }
+    chOldLocal.DoEvaluation(eg);
+    //chOldLocal.DoEvaluation(eg);
+    nsLocal.DoSubEvaluation(eg,child(nsLocal.lfs,Dune::Indices::_0));
 
     for(const auto& ip : rule)
     {
-      // evaluate shape functions and  gradient of shape functions
-      auto& phi = cache.evaluateFunction(ip.position(), femspace.finiteElement().localBasis());
-      auto& gradphihat = cache.evaluateJacobian(ip.position(), femspace.finiteElement().localBasis());
-
-      auto& vphi = vcache.evaluateFunction(ip.position(), vfemspace.finiteElement().localBasis());
-
-      // transform gradients of shape functions to real element
       const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
-      auto gradphi = makeJacobianContainer(femspace);
-      for (size_t i=0; i<femspace.size(); i++)
-        S.mv(gradphihat[i][0],gradphi[i][0]);
-
-      //compute current values of the unknowns
-      Dune::FieldVector<RF,dim> v(0.0);
-      RF pf1=0.0;
-      RF pf1Old=0.0;
-      RF xi1=0.0;
-
-      RF pf2=0.0;
-      RF pf2Old=0.0;
-      RF xi2=0.0;
-
-      RF u=0.0;
-      RF uOld=0.0;
-
-      Dune::FieldVector<RF,dim> gradpf1(0.0);
-      Dune::FieldVector<RF,dim> gradxi1(0.0);
-      Dune::FieldVector<RF,dim> gradpf2(0.0);
-      Dune::FieldVector<RF,dim> gradxi2(0.0);
-
-      Dune::FieldVector<RF,dim> gradu(0.0);
-
-      for (size_t i=0; i<femspace.size(); i++)
-      {
-
-        gradpf1.axpy(x(pf1space,i),gradphi[i][0]);
-        pf1 += x(pf1space,i)*phi[i];
-        gradxi1.axpy(x(xi1space,i),gradphi[i][0]);
-        xi1 += x(xi1space,i)*phi[i];
+      auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
 
-        gradpf2.axpy(x(pf2space,i),gradphi[i][0]);
-        pf2 += x(pf2space,i)*phi[i];
-        gradxi2.axpy(x(xi2space,i),gradphi[i][0]);
-        xi2 += x(xi2space,i)*phi[i];
+      CHVariables_Threephase_U_PFOld_Extension<RF,dim,LFSV> ch(chCache, ip.position(), lfsv, x.base(),chOldLocal.lv, S);
+      NSVariables<RF,dim,decltype(nsLocal.lfs)> ns(vCache, ip.position(), nsLocal.lfs, nsLocal.lv, S);
 
-        gradu.axpy(x(uspace,i),gradphi[i][0]);
-        u += x(uspace,i)*phi[i];
+      RF pf3 = 1-ch.pf-ch.pf2;
+      RF xi3 = -param.tw.Sigma3*(ch.xi/param.tw.Sigma1 + ch.xi2/param.tw.Sigma2);
 
-        pf1Old += chOldLocal[pf1space.localIndex(i)]*phi[i];
-        pf2Old += chOldLocal[pf2space.localIndex(i)]*phi[i];
-        uOld += chOldLocal[uspace.localIndex(i)]*phi[i];
-      }
-      for (size_t i=0; i<vfemspace.size(); i++)
-      {
-        for(size_t k=0; k<dim; k++)
-        {
-          v[k] += nsOldLocal[child(vspace,k).localIndex(i)]*vphi[i];
-        }
-      }
-      RF pf3 = 1-pf1-pf2;
-      RF xi3 = -param.tw.Sigma3*(xi1/param.tw.Sigma1 + xi2/param.tw.Sigma2);
-
-      RF pf_c = pf1;
-      RF pf_c_Old = pf1Old;
+      RF pf_c = ch.pf;
+      RF pf_c_Old = ch.pfOld;
       RF pf_c_reg = pf_c+param.delta;
       RF pf_c_reg_Old = pf_c_Old + param.delta;
 
-      auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
-      for (size_t i=0; i<femspace.size(); i++)
+      for (size_t i=0; i<ch.pfspace.size(); i++)
       {
-        RF R1 = -param.tw.eval_q(pf1,pf2,pf3)/param.eps*
-                  (param.reaction.eval(u)+param.phys.curvatureAlpha*(xi1-xi3))*phi[i];
-        RF J1 = param.mobility.eval(pf1,pf2,pf3)*param.eps/param.tw.Sigma1 * (gradxi1*gradphi[i][0]);
+        RF R1 = 0;//-param.tw.eval_q(ch.pf,ch.pf2,pf3)/param.eps*
+                  //(param.reaction.eval(ch.u)+param.phys.curvatureAlpha*(ch.xi-xi3))*phi[i];
+        RF J1 = param.mobility.eval(ch.pf,ch.pf2,pf3)*param.eps/param.tw.Sigma1 * (ch.gradxi*ch.basis.gradphi[i][0]);
 
       // integrate: (pf-pfold)*phi_i/delta_t + gradpf * grad phi_i + ...
-        r.accumulate(pf1space,i,factor*(
-          (pf1-pf1Old)*phi[i]/param.time.dt
+        r.accumulate(ch.pfspace,i,factor*(
+          (ch.pf-ch.pfOld)*ch.basis.phi[i]/param.time.dt
           +J1
-          -pf1*(v*gradphi[i][0])
+          +(ch.pf*ns.div_v + (ch.gradpf*ns.v))*ch.basis.phi[i]
+          //-pf1*(v*gradphi[i][0])   weak transport
           -R1
         ));
 
-        r.accumulate(pf2space,i,factor*(
-          (pf2-pf2Old)*phi[i]/param.time.dt
-          +param.mobility.eval(pf1,pf2,pf3)*param.eps/param.tw.Sigma2 * (gradxi2*gradphi[i][0])
-          -pf2*(v*gradphi[i][0])
+        r.accumulate(ch.pf2space,i,factor*(
+          (ch.pf2-ch.pf2Old)*ch.basis.phi[i]/param.time.dt
+          +param.mobility.eval(ch.pf,ch.pf2,pf3)*param.eps/param.tw.Sigma2 * (ch.gradxi2*ch.basis.gradphi[i][0])
+          +(ch.pf2*ns.div_v + (ch.gradpf2*ns.v))*ch.basis.phi[i]
+          //-pf2*(v*gradphi[i][0])
         ));
 
-        r.accumulate(xi1space,i,factor*(
-          xi1*phi[i]
-          - 1/param.eps* param.tw.eval_dpf1(pf1,pf2,pf3) *phi[i]
-          - param.eps * param.tw.Sigma1* (gradpf1*gradphi[i][0])
+        r.accumulate(ch.xispace,i,factor*(
+          ch.xi*ch.basis.phi[i]
+          - 1/param.eps* param.tw.eval_dpf1(ch.pf,ch.pf2,pf3) *ch.basis.phi[i]
+          - param.eps * param.tw.Sigma1* (ch.gradpf*ch.basis.gradphi[i][0])
         ));
 
-        r.accumulate(xi2space,i,factor*(
-          xi2*phi[i]
-          - 1/param.eps* param.tw.eval_dpf2(pf1,pf2,pf3) *phi[i]
-          - param.eps * param.tw.Sigma2* (gradpf2*gradphi[i][0])
+        r.accumulate(ch.xi2space,i,factor*(
+          ch.xi2*ch.basis.phi[i]
+          - 1/param.eps* param.tw.eval_dpf2(ch.pf,ch.pf2,pf3) *ch.basis.phi[i]
+          - param.eps * param.tw.Sigma2* (ch.gradpf2*ch.basis.gradphi[i][0])
         ));
 
       // integrate (pf-pfold)*(u-uast)*phi_i/delta_t + pfold*(u-uold)*phi_i/delta_t + D*pf*gradu*grad phi_i = 0
-        r.accumulate(uspace,i,factor*(
+        r.accumulate(ch.uspace,i,factor*(
 
           //(u-param.uAst)*phi[i]/param.dt //For Testing
 
-          (pf_c_reg-pf_c_reg_Old)*(u-param.phys.uAst)*phi[i]/param.time.dt
-          + pf_c_reg_Old*(u-uOld)*phi[i]/param.time.dt
-          + param.phys.D*pf_c_reg_Old*(gradu*gradphi[i][0])
-          + J1*(u-param.phys.uAst)
-          - pf_c*(u-param.phys.uAst)*(v*gradphi[i][0])
+          (pf_c_reg-pf_c_reg_Old)*(ch.u-param.phys.uAst)*ch.basis.phi[i]/param.time.dt
+          + pf_c_reg_Old*(ch.u-ch.uOld)*ch.basis.phi[i]/param.time.dt
+          + param.phys.D*pf_c_reg_Old*(ch.gradu*ch.basis.gradphi[i][0])
+          + J1*(ch.u-param.phys.uAst)
+          //- pf_c*(u-param.phys.uAst)*(v*gradphi[i][0])
+          + (pf_c*(ch.u-param.phys.uAst)*ns.div_v + (ch.u-param.phys.uAst)*(ch.gradpf*ns.v) + pf_c*(ch.gradu*ns.v) )*ch.basis.phi[i]
           ));
       }
     }
diff --git a/dune/phasefield/localoperator/fem_fs_chns_r_cahnhilliard.hh b/dune/phasefield/localoperator/fem_fs_chns_r_cahnhilliard.hh
index 7ac9d0dcbc498e2a99c8fd6e4547de624edbc52a..3c4a2dac985ee072b0d657b6b0cb1841fddea668 100644
--- a/dune/phasefield/localoperator/fem_fs_chns_r_cahnhilliard.hh
+++ b/dune/phasefield/localoperator/fem_fs_chns_r_cahnhilliard.hh
@@ -56,7 +56,7 @@ public:
   {
     //Quadrature Rule
     const int dim = EG::Entity::dimension;
-    const int order = 2; //TODO: Choose more educated
+    const int order = 5; //TODO: Choose more educated
     auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
 
     chOldLocal.DoEvaluation(eg);
@@ -72,13 +72,18 @@ public:
       NSVariables<RF,dim,decltype(nsLocal.lfs)> ns(vCache, ip.position(), nsLocal.lfs, nsLocal.lv, S);
 
       RF pf_c = ch.pf;
+      RF pf_f = ch.pf + 2*param.delta*(1-ch.pf);
+      Dune::FieldVector<RF,dim> gradpf_f(0.0);
+      gradpf_f.axpy((1-2*param.delta),ch.gradpf);
+
       RF pf_c_Old = ch.pfOld;
       RF pf_c_reg = pf_c+param.delta;
       RF pf_c_reg_Old = pf_c_Old + param.delta;
 
       for (size_t i=0; i<ch.pfspace.size(); i++)
       {
-        RF R1 = -std::max(param.dw.eval_q(ch.pf) - 0.05, 0.0)/param.eps*
+        RF Rspeed = param.dw.eval_q((ch.pfOld-0.05)/0.7)/0.7;
+        RF R1 = -Rspeed/param.eps*
                   (param.reaction.eval(ch.u)+param.phys.curvatureAlpha*ch.xi)*ch.basis.phi[i]; //TODO Scale with sigma13
         RF J1 = param.mobility.eval(ch.pf)*param.eps * (ch.gradxi*ch.basis.gradphi[i][0]);
 
@@ -106,7 +111,9 @@ public:
           + pf_c_reg_Old*(ch.u-ch.uOld)*ch.basis.phi[i]/param.time.dt
           + param.phys.D*pf_c_reg_Old*(ch.gradu*ch.basis.gradphi[i][0])
           + J1*(ch.u-param.phys.uAst)
-          + (pf_c*(ch.u-param.phys.uAst)*ns.div_v + (ch.u-param.phys.uAst)*(ch.gradpf*ns.v) + pf_c*(ch.gradu*ns.v) )*ch.basis.phi[i]
+          //+ (pf_c*(ch.u-param.phys.uAst)*ns.div_v + (ch.u-param.phys.uAst)*(ch.gradpf*ns.v) + pf_c*(ch.gradu*ns.v) )*ch.basis.phi[i]
+          + ( pf_f*(ns.v*ch.gradu) )*ch.basis.phi[i]
+          - pf_c_reg*param.reaction.eval_bulkreaction()*ch.basis.phi[i]//Bulk reaction
           ));
       }
     }
diff --git a/dune/phasefield/localoperator/fem_fs_chns_r_instatstokes.hh b/dune/phasefield/localoperator/fem_fs_chns_r_instatstokes.hh
index d4551beb76f90c383eec4aeeb8b84207b09fd0ea..e991d24fbba3d0c50f7e9a172cf98b66e0599001 100644
--- a/dune/phasefield/localoperator/fem_fs_chns_r_instatstokes.hh
+++ b/dune/phasefield/localoperator/fem_fs_chns_r_instatstokes.hh
@@ -71,7 +71,7 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
       {
         //Quadrature Rule
         const int dim = EG::Entity::dimension;
-        const int order = 4; //TODO: Choose more educated
+        const int order = 5; //TODO: Choose more educated
         auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
 
         chLocal.DoEvaluation(eg);
@@ -107,8 +107,8 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
             {
               //SIMPLIFY S=0, TODO: Implement this correctly
               r.accumulate(child(ns.vspace,k),i, (   rho_f_reg*(ns.v[k]-ns.vOld[k])*ns.vbasis.phi[i]/param.time.dt          //TimeEvolution
-                                                - ns.p * (pf_f*ns.vbasis.gradphi[i][0][k] + ns.vbasis.phi[i]*gradpf_f[k])  //Pressure in weak form
-                                                //+ pf_f * (gradp[k] * vphi[i])                           //Pressure in strong form
+                                                //- ns.p * (pf_f*ns.vbasis.gradphi[i][0][k] + ns.vbasis.phi[i]*gradpf_f[k])  //Pressure in weak form
+                                                + pf_f * (ns.gradp[k] * ns.vbasis.phi[i])                           //Pressure in strong form
                                                 + mu_reg * (ns.jacv[k]*ns.vbasis.gradphi[i][0])        //Viscosity
                                                 + param.vDis.eval(pf_f)*ns.v[k]*ns.vbasis.phi[i]              //Velocity Dissipation in Solid
                                                 //- rho*g[k]*vphi[i]                              //Gravitation
diff --git a/dune/phasefield/localoperator/fem_fs_chns_r_thinplateflow.hh b/dune/phasefield/localoperator/fem_fs_chns_r_thinplateflow.hh
new file mode 100644
index 0000000000000000000000000000000000000000..43d85cf860530a899871c16ac44b920cc782cba8
--- /dev/null
+++ b/dune/phasefield/localoperator/fem_fs_chns_r_thinplateflow.hh
@@ -0,0 +1,127 @@
+#pragma once
+
+#include<dune/pdelab/common/quadraturerules.hh>
+#include<dune/pdelab/finiteelement/localbasiscache.hh>
+//#include<dune/pdelab/localoperator/idefault.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/gridfunctionspace/gridfunctionspace.hh>
+
+#include<dune/phasefield/localoperator/fem_base_chns.hh>
+
+
+
+template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
+    class FEM_FS_CHNS_R_ThinPlateFlow :
+      public Dune::PDELab::FullVolumePattern,
+      public Dune::PDELab::LocalOperatorDefaultFlags,
+      public Dune::PDELab::NumericalJacobianVolume<FEM_FS_CHNS_R_ThinPlateFlow<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
+      public Dune::PDELab::NumericalJacobianApplyVolume<FEM_FS_CHNS_R_ThinPlateFlow<Param, CHGFS, CHContainer, NSGFS, NSContainer> >
+    {
+    private:
+      typedef typename CHContainer::field_type RF;
+      typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity   ENTITY;
+
+      CHLocalBasisCache<CHGFS> chCache;
+      VLocalBasisCache<NSGFS> vCache;
+      PLocalBasisCache<NSGFS> pCache;
+
+      typedef GFS_LocalViewer<CHGFS, CHContainer, ENTITY, RF> CHLView;
+      mutable CHLView chLocal;
+      //mutable CHLView chOldLocal;
+
+      typedef GFS_LocalViewer<NSGFS, NSContainer, ENTITY, RF> NSLView;
+      mutable NSLView nsOldLocal;
+
+      Param& param; // parameter functions
+
+
+    public:
+      //! Boundary condition indicator type
+    //  using BC = StokesBoundaryCondition; // ??
+
+
+//      static const bool navier = P::assemble_navier;
+//      static const bool full_tensor = P::assemble_full_tensor;
+
+      // pattern assembly flags
+      enum { doPatternVolume = true };
+
+      // residual assembly flags
+      enum { doAlphaVolume = true };
+      //enum { doLambdaVolume = true };
+
+      //TODO enum { doLambdaBoundary = true };
+
+
+      FEM_FS_CHNS_R_ThinPlateFlow (Param& param_, CHGFS& chGfs_, CHContainer& chV_,
+                         NSGFS& nsGfs_, NSContainer& nsVOld_) :
+        chLocal(chGfs_,chV_),
+        //chOldLocal(chGfs_,chVOld_),
+        nsOldLocal(nsGfs_,nsVOld_),
+        param(param_)
+      {
+      }
+
+
+      template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
+      void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
+      {
+        //Quadrature Rule
+        const int dim = EG::Entity::dimension;
+        const int order = 5; //TODO: Choose more educated
+        auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
+
+        chLocal.DoEvaluation(eg);
+        //chOldLocal.DoEvaluation(eg);
+        nsOldLocal.DoSubEvaluation(eg,child(nsOldLocal.lfs,Dune::Indices::_0));
+
+        for(const auto& ip : rule)
+        {
+          const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
+          auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
+
+          CHVariables<RF,dim,decltype(chLocal.lfs)> ch(chCache, ip.position(), chLocal.lfs, chLocal.lv, S);
+          NSVariables_P_VOld_Extension<RF,dim,LFSV> ns(vCache, pCache, ip.position(), lfsv, x.base(), nsOldLocal.lv, S);
+
+          RF rho_f=param.phys.rho_1*ch.pf;
+          RF rho_f_reg= rho_f + param.delta;
+
+          RF pf_f = ch.pf + 2*param.delta*(1-ch.pf);
+          Dune::FieldVector<RF,dim> gradpf_f(0.0);
+          gradpf_f.axpy((1-2*param.delta),ch.gradpf);
+
+          RF mu_reg = pf_f*param.phys.mu;
+
+          for (size_t i=0; i<ns.pspace.size(); i++)
+          {
+            //r.accumulate(pspace,i, -1.0* pf_f*(v*gradphi[i][0])*factor);
+            r.accumulate(ns.pspace,i, -1.0 * (pf_f*ns.div_v  + (gradpf_f*ns.v) )* ns.pbasis.phi[i] * factor);
+          }
+
+          for (size_t i=0; i<ns.vfemspace.size(); i++)
+          {
+            for(int k=0; k<dim; k++)
+            {
+              //SIMPLIFY S=0, TODO: Implement this correctly
+              r.accumulate(child(ns.vspace,k),i, (   rho_f_reg*(ns.v[k]-ns.vOld[k])*ns.vbasis.phi[i]/param.time.dt          //TimeEvolution
+                                                //- ns.p * (pf_f*ns.vbasis.gradphi[i][0][k] + ns.vbasis.phi[i]*gradpf_f[k])  //Pressure in weak form
+                                                + pf_f * (ns.gradp[k] * ns.vbasis.phi[i])                           //Pressure in strong form
+                                                + mu_reg * (ns.jacv[k]*ns.vbasis.gradphi[i][0])        //Viscosity
+                                                + param.vDis.eval(pf_f)*ns.v[k]*ns.vbasis.phi[i]              //Velocity Dissipation in Solid
+                                                //- rho*g[k]*vphi[i]                              //Gravitation
+                                                //+ S                                             //Surface Tension
+                                                //+ Reaction TODO
+                                                + mu_reg * 8 / (param.CellHeight * param.CellHeight) * ns.v[k]*ns.vbasis.phi[i] //Drag from bottom and top plate
+                                              )*factor);
+            // integrate (grad u)^T * grad phi_i                                                  //Transpose for Viscosity
+            for(int dd=0; dd<dim; ++dd)
+                r.accumulate(child(ns.vspace,k),i, param.phys.mu * (ns.jacv[dd][k] * (ns.vbasis.gradphi[i][0][dd])) * factor);
+            }
+          }
+        }
+      }
+
+    };
diff --git a/dune/phasefield/localoperator/fem_fs_r_only.hh b/dune/phasefield/localoperator/fem_fs_r_only.hh
new file mode 100644
index 0000000000000000000000000000000000000000..902946f1460be96f6eaaff87f46a9169bffeca37
--- /dev/null
+++ b/dune/phasefield/localoperator/fem_fs_r_only.hh
@@ -0,0 +1,83 @@
+#pragma once
+
+#include<dune/pdelab/common/quadraturerules.hh>
+#include<dune/pdelab/finiteelement/localbasiscache.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/gridfunctionspace/gridfunctionspace.hh>
+
+#include<dune/phasefield/localoperator/fem_base_chns.hh>
+
+
+//Parameter class (for functions), FEM space and RangeField
+template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
+class FEM_FS_R_Only :
+  public Dune::PDELab::FullVolumePattern,
+  public Dune::PDELab::LocalOperatorDefaultFlags,
+  public Dune::PDELab::NumericalJacobianVolume<FEM_FS_R_Only<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
+  public Dune::PDELab::NumericalJacobianApplyVolume<FEM_FS_R_Only<Param, CHGFS, CHContainer, NSGFS, NSContainer> >
+{
+private:
+  typedef typename CHContainer::field_type RF;
+  typedef typename CHGFS::Traits::GridViewType::template Codim<0>::Entity   ENTITY;
+
+  CHLocalBasisCache<CHGFS> chCache;
+
+  typedef GFS_LocalViewer<CHGFS, CHContainer, ENTITY, RF> CHLView;
+  mutable CHLView chOldLocal;
+
+  Param& param; // parameter functions
+
+public:
+  enum {doPatternVolume=true};
+  enum {doAlphaVolume=true};
+
+  //constructor
+  FEM_FS_R_Only(Param& param_, CHGFS& chGfs_, CHContainer& chVOld_, NSGFS& nsGfs_, NSContainer& nsV) :
+    chOldLocal(chGfs_,chVOld_),
+    param(param_)
+  {
+  }
+
+
+  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
+  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
+  {
+    //Quadrature Rule
+    const int dim = EG::Entity::dimension;
+    const int order = 5; //TODO: Choose more educated
+    auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
+
+    chOldLocal.DoEvaluation(eg);
+
+    for(const auto& ip : rule)
+    {
+      const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
+      auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
+
+      CHVariables_U_PFOld_Extension<RF,dim,LFSV> ch(chCache, ip.position(), lfsv, x.base(),chOldLocal.lv, S);
+
+      for (size_t i=0; i<ch.pfspace.size(); i++)
+      {
+        RF Rspeed = param.dw.eval_q((ch.pfOld-0.05)/0.7)/0.7;
+        RF R1 = -Rspeed/param.eps*
+                  (param.reaction.eval(ch.uOld))*ch.basis.phi[i];
+
+        r.accumulate(ch.pfspace,i,factor*(
+          (ch.pf-ch.pfOld)*ch.basis.phi[i]/param.time.dt
+          -R1
+        ));
+
+        r.accumulate(ch.xispace,i,factor*(
+          (ch.xi-0)*ch.basis.phi[i]
+        ));
+
+        r.accumulate(ch.uspace,i,factor*(
+          (ch.u-ch.uOld)*ch.basis.phi[i]
+        ));
+      }
+    }
+  }
+};
diff --git a/dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh b/dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh
index 4b79d18f9c7d9f75d187f4abb2949bf957914f70..f87547532e537daaaab63a69d33754913ccd4eee 100644
--- a/dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh
+++ b/dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh
@@ -13,13 +13,14 @@ class ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
 
 public:
   explicit ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES (const GO& grid_operator,
-    unsigned maxiter_ = 5000, unsigned steps_ = 5, int verbose_ = 1)
+    unsigned maxiter_ = 5000, unsigned restartiter_ = 200, unsigned steps_ = 5, int verbose_ = 1)
     : _grid_operator(grid_operator)
     , gfs(grid_operator.trialGridFunctionSpace())
     , phelper(gfs,verbose_)
     , maxiter(maxiter_)
     , steps(steps_)
     , verbose(verbose_)
+    , restartiter(restartiter_)
   {}
 
   template<class Vector>
@@ -62,7 +63,7 @@ public:
 #endif
     int verb=0;
     if (gfs.gridView().comm().rank()==0) verb=verbose;
-    Solver<VectorType, VectorType, VectorType> solver(oop,psp,parsmoother,reduction,100,maxiter,verb);
+    Solver<VectorType, VectorType, VectorType> solver(oop,psp,parsmoother,reduction,restartiter,maxiter,verb);
     Dune::InverseOperatorResult stat;
     //make r consistent
     if (gfs.gridView().comm().size()>1){
@@ -93,4 +94,5 @@ private:
   unsigned maxiter;
   unsigned steps;
   int verbose;
+  unsigned restartiter;
 };
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index fb18a3807bda0ce703f5d5974ae90aac6f881ba0..9ca989cd621852505aefdb396d9ecb751913f389 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -3,3 +3,5 @@ target_link_dune_default_libraries("dune-phasefield")
 add_subdirectory(test_2f0s)
 add_subdirectory(test_3pchns_r)
 add_subdirectory(test_2pchns_r)
+add_subdirectory(test_contactpointdynamics)
+add_subdirectory(summerschool)
diff --git a/src/rectangle_periodic.geo b/src/rectangle_periodic.geo
new file mode 100644
index 0000000000000000000000000000000000000000..c427e578ff5be0fdc905ca27eca3f92a2b3d34c2
--- /dev/null
+++ b/src/rectangle_periodic.geo
@@ -0,0 +1,22 @@
+// mesh width associated with points
+lc = 0.02;
+
+Point(1) = {0, 0, 0, lc};
+Point(2) = {2, 0, 0, lc};
+Point(3) = {2, 1, 0, lc};
+Point(4) = {0, 1, 0, lc};
+
+Line(1) = {2,1};
+Line(2) = {3,2};
+
+Line(3) = {3,4};
+Line(4) = {4,1};
+
+Periodic Line {1} = {3};
+Periodic Line {2} = {4};
+
+Line Loop(100) = {-1,-2,3,4};  
+Plane Surface(200) = {100};
+Physical Surface(300) = {200};
+
+Mesh.Algorithm = 5;
diff --git a/src/summerschool.geo b/src/summerschool.geo
new file mode 100644
index 0000000000000000000000000000000000000000..0c68dd4346d74082ee6a009c781f9b46c72314ac
--- /dev/null
+++ b/src/summerschool.geo
@@ -0,0 +1,35 @@
+// mesh width associated with points
+lc = 0.013;
+
+R = 1/4;  // Radius of pore
+h = 1/8;  // Height of inlet
+ix = 1/2; // Inlet length
+dx = R - Sqrt(R*R - h*h/4); // Additional x-distance to the edge of the circle 
+
+Point(1) = {0,          0,  0, lc};
+Point(2) = {ix+dx,      0,  0, lc};
+Point(3) = {ix+R+R-dx,  0,  0, lc};
+Point(4) = {ix+R+R+ix,  0,  0, lc};
+Point(5) = {ix+R+R+ix,  h,  0, lc};
+Point(6) = {ix+R+R-dx,  h,  0, lc};
+Point(7) = {ix+dx,      h,  0, lc};
+Point(8) = {0,          h,  0, lc};
+
+Point(9) = {ix+R, h/2, 0, lc};
+
+Line(1) = {1,2};
+Circle(2) = {2,9,3};
+Line(3) = {3,4};
+Line(4) = {4,5};
+Line(5) = {5,6};
+Circle(6) = {6,9,7};
+Line(7) = {7,8};
+Line(8) = {1,8};
+
+Periodic Line {4} = {8};
+
+Line Loop(100) = {1,2,3,4,5,6,7,-8};  
+Plane Surface(200) = {100};
+Physical Surface(300) = {200};
+
+Mesh.Algorithm = 5;
diff --git a/src/summerschool/2p_driver.hh b/src/summerschool/2p_driver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3916f9eae2ad206811a7c61e27ee9bb5a3f3937e
--- /dev/null
+++ b/src/summerschool/2p_driver.hh
@@ -0,0 +1,273 @@
+#pragma once
+
+// C includes
+#include<sys/stat.h>
+// C++ includes
+#include<iostream>
+// Dune includes
+#include<dune/grid/io/file/vtk.hh>
+#include<dune/pdelab/newton/newton.hh>
+
+#include <dune/phasefield/localoperator/fem_fs_chns_r_cahnhilliard.hh>
+#include <dune/phasefield/localoperator/fem_fs_chns_r_thinplateflow.hh>
+#include <dune/phasefield/localoperator/fem_fs_r_only.hh>
+#include <dune/phasefield/localoperator/pf_diff_estimator.hh>
+
+#include <dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh>
+
+#include <dune/phasefield/chns_base.hh>
+#include <dune/phasefield/fs_chns_r.hh>
+#include <dune/phasefield/initial_fn.hh>
+#include <dune/phasefield/boundary_fn.hh>
+#include <dune/phasefield/debug/vector_debug.hh>
+
+#include "summerschool_timestepmanager.hh"
+
+
+
+//_________________________________________________________________________________________________________________________
+
+template<typename Grid, typename ES, typename V_FEM, typename P_FEM, typename PHI_FEM>
+void driver_fs (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem, const PHI_FEM& phiFem, Dune::ParameterTree& ptree, Dune::MPIHelper& helper)
+{
+
+  typedef typename Grid::LeafGridView GV;
+  //std::bitset<2> periodic_directions(std::string("01"));
+  //typedef Dune::XT::Grid::PeriodicGridView<LGV> GV;
+  GV gv = grid.leafGridView();
+  const int dim = GV::dimension;
+  typedef double RF;
+
+
+  typedef Params_summerschool< RF,
+                        DoubleWell_limited< RF, DoubleWell_poly<RF>, Limiter_FakeDiverging<RF> >,//DoubleWell_poly<RF>,
+                        Mobility_Const<RF>,
+                        Reaction_Calcite_Precipitation<RF>,
+                        VelDissipation_Quadratic<RF> > Parameters;
+  Parameters param(ptree);
+  SummerschoolTimestepManager timestepManager(ptree.sub("Time"),param.time,2);
+
+
+//_________________________________________________________________________________________________________________________
+// Make grid function spaces
+
+
+  typedef GFS_fs_chns_r<RF, ES, V_FEM, P_FEM, PHI_FEM, PHI_FEM> GFS;
+  GFS gfs(es,vFem,pFem,phiFem,phiFem);
+  typedef typename GFS::NS NS_GFS;
+  typedef typename GFS::CH CH_GFS;
+
+
+  //_________________________________________________________________________________________________________________________
+  //Constraints
+
+  typedef Bdry_fs_chns_r_compositeV_2d<InflowOutflow_v,
+              Dune::PDELab::DirichletConstraintsParameters,
+              InflowOutflow_p,
+              InflowOutflow_phi,
+              Dune::PDELab::NoDirichletConstraintsParameters,
+              InflowOutflow_u> Bdry;
+  Bdry bdry;
+  gfs.updateConstraintsContainer(bdry);
+
+  //_________________________________________________________________________________________________________________________
+  //Coefficient vectors
+
+  using CH_V = Dune::PDELab::Backend::Vector<CH_GFS, RF>;
+  using NS_V = Dune::PDELab::Backend::Vector<NS_GFS, RF>;
+  CH_V chV(gfs.ch);
+  NS_V nsV(gfs.ns);
+
+  //_________________________________________________________________________________________________________________________
+  //Initial Values
+
+    typedef Ini_fs_chns_r< VelocityInitial_Inflow<GV,RF,Parameters,dim>,
+                      PInitial_Inflow<GV,RF,Parameters>,
+                      PhiInitial_SummerSchoolFull<GV,RF,Parameters>,
+                      ZeroInitial<GV,RF,Parameters>,
+                      UInitial_Inflow<GV,RF,Parameters> >  Initial;
+    Initial initial(gv,param);
+    Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
+    Dune::PDELab::interpolate(initial.iniNS , gfs.ns, nsV);
+
+  //____________________________________________________________________________________________________________________
+  // Refine Initial Data
+
+  for(int i=0; i<ptree.get("Refinement.maxRefineLevel",2); i++)
+  {
+    // mark elements for refinement
+    typedef RefinementIndicator<GV,RF,PF_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND;
+    RIND rInd(gv,gfs.ch,ptree.sub("Refinement"));
+    rInd.calculateIndicator(gv,chV,false);
+    rInd.markGrid(grid,2);
+
+    // do refinement
+    auto chT = transferSolutions(gfs.ch,4,chV);
+    auto nsT = transferSolutions(gfs.ns,4,nsV);
+    Dune::PDELab::adaptGrid(grid, chT,nsT);
+
+    Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
+    Dune::PDELab::interpolate(initial.iniNS , gfs.ns, nsV);
+  }
+
+
+  gfs.updateConstraintsContainer(bdry);
+  interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
+
+  CH_V chVOld = chV;
+  NS_V nsVOld = nsV;
+
+  //_________________________________________________________________________________________________________________________
+  // LocalOperators
+
+  typedef FEM_FS_CHNS_R_ThinPlateFlow<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> NS_LOP;
+  NS_LOP nsLop(param, gfs.ch, chV, gfs.ns, nsVOld);
+
+  typedef FEM_FS_CHNS_R_CahnHilliard<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> CH_LOP;
+  CH_LOP chLop(param, gfs.ch, chVOld, gfs.ns, nsV);
+
+  typedef FEM_FS_R_Only<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> REACTION_LOP;
+  REACTION_LOP reactionLop(param, gfs.ch, chVOld, gfs.ns, nsV);
+
+  typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
+  MBE mbe(10);// Maximal number of nonzeroes per row can be cross-checked by patternStatistics().
+
+
+  //_________________________________________________________________________________________________________________________
+  // Construct discrete Grid Functions
+  typedef DGF_fs_chns_r<GFS, NS_V, CH_V> DGF;
+  DGF dgf(gfs, nsV, chV);
+
+  //_________________________________________________________________________________________________________________________
+  // prepare VTK writer
+  std::string filename=ptree.get("output.filename","output");
+  struct stat st;
+  if( stat( filename.c_str(), &st ) != 0 )
+  {
+    int stat = 0;
+    stat = mkdir(filename.c_str(),S_IRWXU|S_IRWXG|S_IRWXO);
+    if( stat != 0 && stat != -1)
+      std::cout << "Error: Cannot create directory "  << filename << std::endl;
+  }
+  auto stationaryvtkwriter = std::make_shared<Dune::SubsamplingVTKWriter<GV> >(gv,1);
+  Dune::VTKSequenceWriter<GV> vtkwriter(stationaryvtkwriter,filename.c_str(),filename.c_str(),"");
+  dgf.addToVTKWriter(vtkwriter);
+
+
+
+  vtkwriter.write(param.time.t,Dune::VTK::ascii);
+
+  //_________________________________________________________________________________________________________________________
+  printVector(chV,10,"chV");
+  printVector(nsV,10,"nsV");
+
+
+  //____________________________________________________________________________________________________________________
+  //int i_it=0;
+
+
+  while (!timestepManager.isFinished())
+    {
+
+      // Grid Operators
+      typedef Dune::PDELab::GridOperator<NS_GFS,NS_GFS,NS_LOP,MBE,RF,RF,RF,typename GFS::NS_CC,typename GFS::NS_CC> NS_GO;
+      NS_GO nsGo(gfs.ns,gfs.nsCC,gfs.ns,gfs.nsCC,nsLop,mbe);
+
+      typedef Dune::PDELab::GridOperator<CH_GFS,CH_GFS,CH_LOP,MBE,RF,RF,RF,typename GFS::CH_CC,typename GFS::CH_CC> CH_GO;
+      CH_GO chGo(gfs.ch,gfs.chCC,gfs.ch,gfs.chCC,chLop,mbe);
+
+      typedef Dune::PDELab::GridOperator<CH_GFS,CH_GFS,REACTION_LOP,MBE,RF,RF,RF,typename GFS::CH_CC,typename GFS::CH_CC> REACTION_GO;
+      REACTION_GO reactionGo(gfs.ch,gfs.chCC,gfs.ch,gfs.chCC,reactionLop,mbe);
+
+      // Linear solver
+      typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
+                  <CH_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> CH_LS; CH_LS ls_ch(chGo,1000,200,3,1);
+
+      typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
+                  <NS_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> NS_LS; NS_LS ls_ns(nsGo,4000,400,3,1);
+
+      typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
+                  <REACTION_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> REACTION_LS; REACTION_LS ls_reaction(reactionGo,500,100,3,1);
+
+      // Nonlinear solver
+      typedef Dune::PDELab::Newton<NS_GO,NS_LS,NS_V> PDESolver1;
+      PDESolver1 ns_newton(nsGo,nsV,ls_ns);
+      ns_newton.setParameters(ptree.sub("NS_Newton"));
+
+      typedef Dune::PDELab::Newton<CH_GO,CH_LS,CH_V> PDESolver2;
+      PDESolver2 ch_newton(chGo,chV,ls_ch);
+      ch_newton.setParameters(ptree.sub("Phasefield_Newton"));
+
+      typedef Dune::PDELab::Newton<REACTION_GO,REACTION_LS,CH_V> PDESolver3;
+      PDESolver3 reaction_newton(reactionGo,chV,ls_reaction);
+      reaction_newton.setParameters(ptree.sub("Phasefield_Newton"));
+
+      //_________________________________________________________________________________________________________________________
+      // do time step
+
+
+      interpolateToBdry(gfs,initial.iniNS,initial.iniCH, nsV, chV);
+
+      int result;
+      if((param.time.t > ptree.get("Time.reactionOnlyStart",(double)5) )
+          && (timestepManager.lastTimestepReactionOnly == false)
+          && (timestepManager.dtMax >= param.time.dt)
+        )
+      {
+        result = timestepManager.doReactionOnlyTimestep(reaction_newton,vtkwriter);
+      }
+      else
+      {
+        result = timestepManager.doTimestep(ns_newton,ch_newton,vtkwriter);
+        if(result > 0)
+        {
+          timestepManager.lastTimestepReactionOnly = false;
+        }
+      }
+
+      if(result <0)
+      {
+        return;
+      }
+      else if(result == 0)
+      {
+        nsV = nsVOld;
+        chV = chVOld;
+        continue;
+      }
+
+      chVOld = chV;
+      nsVOld = nsV;
+
+      //_________________________________________________________________________________________________________________________
+
+      // mark elements for refinement
+      typedef RefinementIndicator<GV,RF,PF_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND;
+      RIND rInd(gv,gfs.ch,ptree.sub("Refinement"));
+      rInd.calculateIndicator(gv,chV,false);
+      rInd.markGrid(grid,2);
+
+      printVector(chV,3,"chV");
+      printVector(nsV,3,"nsV");
+
+      // do refinement
+      std::cout << "adapt grid and solution...   ";
+      auto chTransfer = transferSolutions(gfs.ch,4,chV,chVOld);
+      auto nsTransfer = transferSolutions(gfs.ns,4,nsV,nsVOld);
+      Dune::PDELab::adaptGrid(grid, chTransfer,nsTransfer);
+      std::cout << "   ...adapt ended" << std::endl;
+
+      printVector(chV,3,"chV");
+      printVector(nsV,3,"nsV");
+
+      //vtkwriter.write(timestepManager.t+1e-6,Dune::VTK::ascii);
+
+      // recompute constraints
+      gfs.updateConstraintsContainer(bdry);
+      interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
+      interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsVOld,chVOld);
+
+
+      //vtkwriter.write(timestepManager.t+2e-6,Dune::VTK::ascii);
+    }
+}
diff --git a/src/summerschool/CMakeLists.txt b/src/summerschool/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..c08acc1adbb08518a80770c55d707e9785a50045
--- /dev/null
+++ b/src/summerschool/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_executable(summerschool summerschool.cc)
+dune_symlink_to_source_files(FILES grids ss_calcite.ini)
diff --git a/src/summerschool/ss_calcite.ini b/src/summerschool/ss_calcite.ini
new file mode 100644
index 0000000000000000000000000000000000000000..afa6379272077f55723e30caebc04917a0858e1c
--- /dev/null
+++ b/src/summerschool/ss_calcite.ini
@@ -0,0 +1,71 @@
+[Physics]
+uAst = 27.1	       #Molar density of CaCO3, in mol/l
+uEquilibrium = 0.00014 #Solubility of CaCO3 in water, in mol/l
+ReactionRate = 10
+BulkRate = 0.00001796 #Urease reaction rate, depending un urea concentration and experimentally matched, in mol/(l*s)
+D = 0.000804 #of CO3^2- in Water see Zeebe,2011, in mm^2/s
+rho1 = 0.997 #Water, in 10^-6 kg/mm^3
+rho2 = 2.71  #Kalcite, in 10^-6 kg/mm^3
+mu=1.01      #Viscosity of Water, in mPa*s = 10^-6 kg/mm
+VelDissipationRate = 1000000
+CurvatureAlpha=0
+
+[domain]
+filename = grids/summerschool_full_coarse.msh
+CellHeight=0.085
+
+
+[Phasefield]
+eps=0.002
+delta=0.012
+Mobility=0.0001#0.001 * 0.1
+DoubleWellDelta=0.1
+DoubleWellGamma=0.08
+
+[Time]
+tMax = 100000 #25000 * 4
+dt = 0.0015
+dtMax = 0.1
+dtMin = 0.00015
+saveintervallstart = 100
+saveintervall = 49.99
+reactionOnlyDt = 10
+reactionOnlyStart = 2
+#dtIncreaseFactor=1.2
+#dtDecreaseFactor=0.5
+
+[Initial]
+MeanInflow=2 #mm/s (in mean-> *1.5 for max velocity)
+uConst=0.00020
+uInflow=0.00050 #TODO: Make a parameter study
+Radius=0.02
+
+Radius1=0.015
+x1Center=0.8
+y1Center=-0.1
+Radius2=0.015
+x2Center=0.6
+y2Center=0.11
+Radius3=0.015
+x3Center=0.8
+y3Center=0.225
+
+[Refinement]
+etaRefine  = 0.30
+etaCoarsen = 0.15
+maxRefineLevel = 7
+minRefineLevel = 0
+
+[NS_Newton]
+ReassembleThreshold = 0.0
+LineSearchMaxIterations = 30
+MaxIterations = 30
+AbsoluteLimit = 1e-9
+Reduction = 1e-7
+VerbosityLevel = 1
+
+[Phasefield_Newton]
+ReassembleThreshold = 0.0
+LineSearchMaxIterations = 30
+MaxIterations = 30
+VerbosityLevel = 1
diff --git a/src/summerschool/summerschool.cc b/src/summerschool/summerschool.cc
new file mode 100644
index 0000000000000000000000000000000000000000..1072646d95f72287498e82e190aa9009f226d27a
--- /dev/null
+++ b/src/summerschool/summerschool.cc
@@ -0,0 +1,99 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include<iostream>
+// dune includes
+#include<dune/common/parallel/mpihelper.hh>
+#include<dune/common/parametertreeparser.hh>
+#include<dune/common/timer.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
+// pdelab includes
+#include<dune/pdelab/finiteelementmap/pkfem.hh>
+
+#include "2p_driver.hh"
+
+//===============================================================
+// Main program with grid setup
+//===============================================================
+
+
+
+
+int main(int argc, char** argv)
+{
+  #if !HAVE_SUPERLU
+    std::cout << "Error: These examples work only if SuperLU is available." << std::endl;
+    exit(1);
+  #endif
+
+  #if HAVE_UG
+  try{
+    // Maybe initialize Mpi
+    Dune::MPIHelper&
+      helper = Dune::MPIHelper::instance(argc, argv);
+    if(Dune::MPIHelper::isFake)
+      std::cout<< "This is a sequential program." << std::endl;
+    else
+      std::cout << "Parallel code run on "
+                << helper.size() << " process(es)" << std::endl;
+
+    // open ini file
+    Dune::ParameterTree ptree;
+    Dune::ParameterTreeParser ptreeparser;
+    ptreeparser.readINITree("ss_calcite.ini",ptree);
+    ptreeparser.readOptions(argc,argv,ptree);
+
+    typedef Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming> Grid;
+    //typedef Dune::UGGrid<2> Grid;
+    typedef Grid::ctype DF;
+
+    std::vector<int> boundary_index_map;
+    std::vector<int> element_index_map;
+
+    std::string filename = ptree.get("domain.filename",
+                                     "turbtube2d.msh");
+    Dune::GridFactory<Grid> factory;
+    if (helper.rank()==0)
+      Dune::GmshReader<Grid>::read(factory,filename,boundary_index_map,element_index_map,true,false);
+    std::shared_ptr<Grid> gridp(factory.createGrid());
+
+    typedef Grid::LeafGridView GV;
+    GV gv = gridp->leafGridView();
+
+    std::cout << "OK" << std::endl;
+
+    using ES = Dune::PDELab::NonOverlappingEntitySet<GV>;
+    ES es(gv);
+
+    const int k = 2;
+    typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,k> V_FEM;
+    typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,k-1> P_FEM;
+    const int l = 2;
+    typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,l> PHI_FEM;
+
+    V_FEM vFem(es);
+    P_FEM pFem(es);
+    PHI_FEM phiFem(es);
+
+    driver_fs(*gridp, es,vFem,pFem,phiFem, ptree, helper);
+
+
+  }
+  catch (Dune::Exception &e){
+    std::cerr << "Dune reported error: " << e << std::endl;
+    return 1;
+  }
+  catch (...){
+    std::cerr << "Unknown exception thrown!" << std::endl;
+    return 1;
+  }
+  #endif
+}
diff --git a/src/summerschool/summerschool_timestepmanager.hh b/src/summerschool/summerschool_timestepmanager.hh
new file mode 100644
index 0000000000000000000000000000000000000000..cb751f7b501c06a7788b4a3bfa96d4aa723e8e84
--- /dev/null
+++ b/src/summerschool/summerschool_timestepmanager.hh
@@ -0,0 +1,59 @@
+#pragma once
+
+#include<dune/common/parametertreeparser.hh>
+
+#include<dune/pdelab/finiteelementmap/p0fem.hh>
+#include <dune/pdelab/adaptivity/adaptivity.hh>
+#include<dune/pdelab/gridoperator/gridoperator.hh>
+#include<dune/pdelab/backend/istl.hh>
+
+#include <dune/phasefield/base_parameters.hh>
+#include <dune/phasefield/chns_base.hh>
+
+//_________________________________________________________________________________________________________________________
+
+class SummerschoolTimestepManager :
+public TimestepManager
+{
+public:
+  double reactionOnlyDt;
+  bool lastTimestepReactionOnly;
+  SummerschoolTimestepManager(const Dune::ParameterTree & param, Params_Time& time_, int verbosity) :
+  TimestepManager(param, time_, verbosity)
+  {
+    reactionOnlyDt = param.get("reactionOnlyDt",(double)30);
+    lastTimestepReactionOnly = false;
+  }
+
+  template<typename REACTION_NEWTON, typename VTKWRITER>
+  int doReactionOnlyTimestep(REACTION_NEWTON& reaction_newton, VTKWRITER& vtkwriter)
+  {
+    try
+    {
+      if (verb > 1)   std::cout << "= Begin newtonapply REACTION:" << std::endl;
+
+      double tmpdt = time.dt;
+      time.dt = reactionOnlyDt;
+      reaction_newton.apply();
+      time.dt = tmpdt;
+    }
+    catch(...)
+    {
+      std::cout << "===============================" << std::endl;
+      std::cout << "Problem with Newton Solver, aborting :(" << std::endl;
+      std::cout << "===============================" << std::endl;
+      return -1;
+    }
+
+    lastTimestepReactionOnly = true;
+    time.t += reactionOnlyDt;
+
+    if( time.t < saveintervallstart || time.t>lastsave+2*saveintervall ) // Write reaction steps only if there are no other steps inbetween
+    {
+        lastsave = time.t;
+        vtkwriter.write(time.t,Dune::VTK::ascii);
+    }
+    return 1;
+  }
+
+};
diff --git a/src/summerschool_full.geo b/src/summerschool_full.geo
new file mode 100644
index 0000000000000000000000000000000000000000..383a89c605829c1846ae833dc1e576904c17279a
--- /dev/null
+++ b/src/summerschool_full.geo
@@ -0,0 +1,80 @@
+// mesh width associated with points
+lc = 0.015;
+
+R = 1/4;  // Radius of pore
+h = 1/8;  // Height of inlet
+ix = 1/2; // Inlet length
+dx = R - Sqrt(R*R - h*h/4); // Additional x-distance to the edge of the circle 
+Wp = ix*2+R*2; // Unit cell length
+
+Point(1) = {0,          0,  0, lc};
+Point(2) = {ix+dx,      0,  0, lc};
+Point(3) = {ix+R+R-dx,  0,  0, lc};
+
+Point(12) = {Wp+ix+dx,      0,  0, lc};
+Point(13) = {Wp+ix+R+R-dx,  0,  0, lc};
+
+Point(22) = {Wp+Wp+ix+dx,      0,  0, lc};
+Point(23) = {Wp+Wp+ix+R+R-dx,  0,  0, lc};
+
+Point(32) = {Wp+Wp+Wp+ix+dx,      0,  0, lc};
+Point(33) = {Wp+Wp+Wp+ix+R+R-dx,  0,  0, lc};
+Point(34) = {Wp+Wp+Wp+ix+R+R+ix,  0,  0, lc};
+
+Point(35) = {Wp+Wp+Wp+ix+R+R+ix,  h,  0, lc};
+Point(36) = {Wp+Wp+Wp+ix+R+R-dx,  h,  0, lc};
+Point(37) = {Wp+Wp+Wp+ix+dx,      h,  0, lc};
+
+Point(26) = {Wp+Wp+ix+R+R-dx,  h,  0, lc};
+Point(27) = {Wp+Wp+ix+dx,      h,  0, lc};
+
+Point(16) = {Wp+ix+R+R-dx,  h,  0, lc};
+Point(17) = {Wp+ix+dx,      h,  0, lc};
+
+Point(6) = {ix+R+R-dx,  h,  0, lc};
+Point(7) = {ix+dx,      h,  0, lc};
+Point(8) = {0,          h,  0, lc};
+
+Point(9) = {ix+R, h/2, 0, lc};
+Point(19) = {Wp+ix+R, h/2, 0, lc};
+Point(29) = {Wp+Wp+ix+R, h/2, 0, lc};
+Point(39) = {Wp+Wp+Wp+ix+R, h/2, 0, lc};
+
+
+Line(1) = {1,2};
+Circle(2) = {2,9,3};
+Line(3) = {3,12};
+
+Circle(12) = {12,19,13};
+Line(13) = {13,22};
+
+Circle(22) = {22,29,23};
+Line(23) = {23,32};
+
+Circle(32) = {32,39,33};
+Line(33) = {33,34};
+
+Line(4) = {34,35};
+
+Line(35) = {35,36};
+Circle(36) = {36,39,37};
+Line(37) = {37,26};
+
+Circle(26) = {26,29,27};
+Line(27) = {27,16};
+
+Circle(16) = {16,19,17};
+Line(17) = {17,6};
+
+Circle(6) = {6,9,7};
+Line(7) = {7,8};
+
+Line(8) = {1,8};
+
+Periodic Line {4} = {8};
+
+Line Loop(100) = {1,2,3,12,13,22,23,32,33,4,35,36,37,26,27,16,17,6,7,-8};  
+Plane Surface(200) = {100};
+Physical Surface(300) = {200};
+
+Mesh.Algorithm = 5;
diff --git a/src/test_2pchns_r/2p_driver.hh b/src/test_2pchns_r/2p_driver.hh
index c636cd8d830765e4599fa7e0e196ed7af4a374c9..3c22e70039e9306a0ef0096ceb870dd452ebe2f8 100644
--- a/src/test_2pchns_r/2p_driver.hh
+++ b/src/test_2pchns_r/2p_driver.hh
@@ -79,7 +79,7 @@ void driver_fs (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem, c
   //Initial Values
 
     typedef Ini_fs_chns_r< VelocityInitial_Inflow<GV,RF,Parameters,dim>,
-                      ZeroInitial<GV,RF,Parameters>,
+                      PInitial_Inflow<GV,RF,Parameters>,
                       PhiInitial_Rectangle<GV,RF,Parameters>,
                       ZeroInitial<GV,RF,Parameters>,
                       UInitial_Inflow<GV,RF,Parameters> >  Initial;
@@ -172,10 +172,10 @@ void driver_fs (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem, c
 
       // Linear solver
       typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
-                  <CH_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> CH_LS; CH_LS ls_ch(chGo,500,3,1);
+                  <CH_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> CH_LS; CH_LS ls_ch(chGo,500,100,3,1);
 
       typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
-                  <NS_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> NS_LS; NS_LS ls_ns(nsGo,500,3,1);
+                  <NS_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> NS_LS; NS_LS ls_ns(nsGo,1000,200,3,1);
 
       // Nonlinear solver
       typedef Dune::PDELab::Newton<NS_GO,NS_LS,NS_V> PDESolver1;
diff --git a/src/test_2pchns_r/2pchns_r.ini b/src/test_2pchns_r/2pchns_r.ini
index c9bb158fb539ad75c78c2e0895ea22239fa10a76..c6bbba6794b9abb2f88ff87c20b9a3fa23c87ed7 100644
--- a/src/test_2pchns_r/2pchns_r.ini
+++ b/src/test_2pchns_r/2pchns_r.ini
@@ -6,7 +6,7 @@ rho2 = 1.0
 mu=0.01
 ReactionRate=0.1
 #ReactionLimiter = 1000
-VelDissipationRate = 100
+VelDissipationRate = 1000
 CurvatureAlpha=0.001
 
 [domain]
@@ -16,12 +16,13 @@ level = 0
 #filename = grids/square_7k_del.msh
 #filename = grids/rectangle_3k_per.msh
 #filename = grids/rectangle_53k_per.msh
-filename = grids/square_2k_del_per.msh
+#filename = grids/square_2k_del_per.msh
+filename = grids/summerschool.msh
 
 [Phasefield]
-eps=0.005
+eps=0.002
 delta=0.01
-Mobility=0.005
+Mobility=0.0005
 DoubleWellDelta=0.005
 DoubleWellGamma=0.0045
 
@@ -35,10 +36,12 @@ dtMin = 0.0015
 #dtDecreaseFactor=0.5
 
 [Initial]
-MeanInflow=1
+MeanInflow=10000
 uConst=1.0
 uInflow=1.3
-Radius=0.15
+Radius=0.05
+xCenter=0.75
+yCenter=0.1
 
 [Refinement]
 etaRefine  = 0.2
diff --git a/src/test_contactpointdynamics/3p_driver.hh b/src/test_contactpointdynamics/3p_driver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a8c8440b2dcf83b8a9a29940831923f1cb1a2579
--- /dev/null
+++ b/src/test_contactpointdynamics/3p_driver.hh
@@ -0,0 +1,218 @@
+#pragma once
+
+// C includes
+#include<sys/stat.h>
+// C++ includes
+#include<iostream>
+// Dune includes
+#include<dune/grid/io/file/vtk.hh>
+#include<dune/pdelab/newton/newton.hh>
+
+#include <dune/phasefield/localoperator/fem_ffs_chns_r_cahnhilliard.hh>
+#include <dune/phasefield/localoperator/fem_ffs_chns_navierstokes.hh>
+#include <dune/phasefield/localoperator/pf_3p_diff_estimator.hh>
+
+#include <dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh>
+
+#include <dune/phasefield/chns_base.hh>
+#include <dune/phasefield/ffs_chns_r.hh>
+#include <dune/phasefield/initial_fn.hh>
+#include <dune/phasefield/debug/vector_debug.hh>
+
+
+
+//_________________________________________________________________________________________________________________________
+
+template<typename Grid, typename ES, typename V_FEM, typename P_FEM, typename PHI_FEM>
+void driver_3p (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem, const PHI_FEM& phiFem, Dune::ParameterTree& ptree, Dune::MPIHelper& helper)
+{
+
+  typedef typename Grid::LeafGridView GV;
+  //std::bitset<2> periodic_directions(std::string("01"));
+  //typedef Dune::XT::Grid::PeriodicGridView<LGV> GV;
+  GV gv = grid.leafGridView();
+  const int dim = GV::dimension;
+  typedef double RF;
+
+
+  typedef Params_ffs_r< RF,
+                        TripleWell<RF,DoubleWell_limited< RF, DoubleWell_poly<RF>, Limiter_FakeDiverging<RF> > >,
+                        Mobility_Const<RF>,
+                        Reaction_None<RF>,
+                        VelDissipation_Quadratic<RF> > Parameters;
+  Parameters param(ptree);
+  TimestepManager timestepManager(ptree.sub("Time"),param.time,2);
+
+
+//_________________________________________________________________________________________________________________________
+// Make grid function spaces
+
+
+  typedef GFS_ffs_chns_r<RF, ES, V_FEM, P_FEM, PHI_FEM, PHI_FEM> GFS;
+  GFS gfs(es,vFem,pFem,phiFem,phiFem);
+  typedef typename GFS::NS NS_GFS;
+  typedef typename GFS::CH CH_GFS;
+
+
+  //_________________________________________________________________________________________________________________________
+  //Constraints
+
+  typedef Bdry_ffs_chns_r_compositeV_2d<InflowOutflow_v,            // v_x
+              Dune::PDELab::DirichletConstraintsParameters,         // v_y
+              InflowOutflow_p,                                      // p
+              InflowOutflow_phi,                                    // phi1
+              Dune::PDELab::NoDirichletConstraintsParameters,       // mu1
+              InflowOutflow_phi,                                    // phi2
+              Dune::PDELab::NoDirichletConstraintsParameters,       // mu2
+              InflowOutflow_u>                                      // u
+              Bdry;
+  Bdry bdry;
+  gfs.updateConstraintsContainer(bdry);
+
+  //_________________________________________________________________________________________________________________________
+  //Coefficient vectors
+
+  using CH_V = Dune::PDELab::Backend::Vector<CH_GFS, RF>;
+  using NS_V = Dune::PDELab::Backend::Vector<NS_GFS, RF>;
+  CH_V chV(gfs.ch);
+  NS_V nsV(gfs.ns);
+  CH_V chVOld = chV;
+  NS_V nsVOld = nsV;
+
+  //_________________________________________________________________________________________________________________________
+  //Initial Values
+
+    typedef Ini_ffs_chns_r< VelocityInitial_Inflow<GV,RF,Parameters,dim>,
+                      ZeroInitial<GV,RF,Parameters>,
+                      PhiInitial_3pFlowTest<GV,RF,Parameters,0>,
+                      ZeroInitial<GV,RF,Parameters>,
+                      PhiInitial_3pFlowTest<GV,RF,Parameters,1>,
+                      ZeroInitial<GV,RF,Parameters>,
+                      ZeroInitial<GV,RF,Parameters> >  Initial;
+    Initial initial(gv,param);
+    Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
+    Dune::PDELab::interpolate(initial.iniNS , gfs.ns, nsV);
+
+  //_________________________________________________________________________________________________________________________
+  // LocalOperators
+
+  typedef FEM_FFS_CHNS_NavierStokes<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> NS_LOP;
+  NS_LOP nsLop(param, gfs.ch, chV, chVOld, gfs.ns, nsVOld);
+
+  typedef FEM_FFS_CHNS_R_CahnHilliard<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> CH_LOP;
+  CH_LOP chLop(param, gfs.ch, chVOld, gfs.ns, nsV);
+
+  typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
+  MBE mbe(10);// Maximal number of nonzeroes per row can be cross-checked by patternStatistics().
+
+
+  //_________________________________________________________________________________________________________________________
+  // Construct discrete Grid Functions
+  typedef DGF_ffs_chns_r<GFS, NS_V, CH_V> DGF;
+  DGF dgf(gfs, nsV, chV);
+
+  //_________________________________________________________________________________________________________________________
+  // prepare VTK writer
+  std::string filename=ptree.get("output.filename","output");
+  struct stat st;
+  if( stat( filename.c_str(), &st ) != 0 )
+  {
+    int stat = 0;
+    stat = mkdir(filename.c_str(),S_IRWXU|S_IRWXG|S_IRWXO);
+    if( stat != 0 && stat != -1)
+      std::cout << "Error: Cannot create directory "  << filename << std::endl;
+  }
+  auto stationaryvtkwriter = std::make_shared<Dune::SubsamplingVTKWriter<GV> >(gv,0);
+  Dune::VTKSequenceWriter<GV> vtkwriter(stationaryvtkwriter,filename.c_str(),filename.c_str(),"");
+  dgf.addToVTKWriter(vtkwriter);
+
+  vtkwriter.write(param.time.t,Dune::VTK::ascii);
+
+//_________________________________________________________________________________________________________________________
+  printVector(chV,10,"chV");
+  printVector(nsV,10,"nsV");
+
+
+  //____________________________________________________________________________________________________________________
+  //int i_it=0;
+  chVOld = chV;
+  nsVOld = nsV;
+  while (!timestepManager.isFinished())
+    {
+
+      // Grid Operators
+      typedef Dune::PDELab::GridOperator<NS_GFS,NS_GFS,NS_LOP,MBE,RF,RF,RF,typename GFS::NS_CC,typename GFS::NS_CC> NS_GO;
+      NS_GO nsGo(gfs.ns,gfs.nsCC,gfs.ns,gfs.nsCC,nsLop,mbe);
+
+      typedef Dune::PDELab::GridOperator<CH_GFS,CH_GFS,CH_LOP,MBE,RF,RF,RF,typename GFS::CH_CC,typename GFS::CH_CC> CH_GO;
+      CH_GO chGo(gfs.ch,gfs.chCC,gfs.ch,gfs.chCC,chLop,mbe);
+
+      // Linear solver
+      typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
+                  <CH_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> CH_LS; CH_LS ls_ch(chGo,500,100,3,1);
+
+      typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
+                  <NS_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> NS_LS; NS_LS ls_ns(nsGo,1000,200,3,1);
+
+      // Nonlinear solver
+      typedef Dune::PDELab::Newton<NS_GO,NS_LS,NS_V> PDESolver1;
+      PDESolver1 ns_newton(nsGo,nsV,ls_ns);
+      ns_newton.setParameters(ptree.sub("NS_Newton"));
+
+      typedef Dune::PDELab::Newton<CH_GO,CH_LS,CH_V> PDESolver2;
+      PDESolver2 ch_newton(chGo,chV,ls_ch);
+      ch_newton.setParameters(ptree.sub("Phasefield_Newton"));
+
+      //_________________________________________________________________________________________________________________________
+      // do time step
+
+
+      interpolateToBdry(gfs,initial.iniNS,initial.iniCH, nsV, chV);
+
+      int result = timestepManager.doTimestep(ns_newton,ch_newton,vtkwriter);
+      if(result <0)
+      {
+        return;
+      }
+      else if(result == 0)
+      {
+        nsV = nsVOld;
+        chV = chVOld;
+        continue;
+      }
+
+      chVOld = chV;
+      nsVOld = nsV;
+
+      //_________________________________________________________________________________________________________________________
+
+      // mark elements for refinement
+      typedef RefinementIndicator<GV,RF,PF_3P_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND;
+      RIND rInd(gv,gfs.ch,ptree.sub("Refinement"));
+      rInd.calculateIndicator(gv,chV,false);
+      rInd.markGrid(grid,2);
+
+      printVector(chV,3,"chV");
+      printVector(nsV,3,"nsV");
+
+      // do refinement
+      std::cout << "adapt grid and solution" << std::endl;
+      auto chTransfer = transferSolutions(gfs.ch,4,chV,chVOld);
+      auto nsTransfer = transferSolutions(gfs.ns,4,nsV,nsVOld);
+      Dune::PDELab::adaptGrid(grid, chTransfer,nsTransfer);
+      std::cout << "adapt grid and solution ended" << std::endl;
+
+      printVector(chV,3,"chV");
+      printVector(nsV,3,"nsV");
+
+      //vtkwriter.write(timestepManager.t+1e-6,Dune::VTK::ascii);
+
+      // recompute constraints
+      gfs.updateConstraintsContainer(bdry);
+      interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
+      interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsVOld,chVOld);
+
+
+      //vtkwriter.write(timestepManager.t+2e-6,Dune::VTK::ascii);
+    }
+}
diff --git a/src/test_contactpointdynamics/CMakeLists.txt b/src/test_contactpointdynamics/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..828a797599b3a215214638ea659c666f5391cd38
--- /dev/null
+++ b/src/test_contactpointdynamics/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_executable(cp_equation cp_equation.cc)
+dune_symlink_to_source_files(FILES grids chns.ini)
diff --git a/src/test_contactpointdynamics/chns.ini b/src/test_contactpointdynamics/chns.ini
new file mode 100644
index 0000000000000000000000000000000000000000..86d599d5909bccf1d7db22c1daa8a45cdfd878b0
--- /dev/null
+++ b/src/test_contactpointdynamics/chns.ini
@@ -0,0 +1,56 @@
+[Physics]
+uAst = 2
+D = 0.01
+rho1 = 1.0
+rho2 = 1.0
+mu=0.005
+ReactionRate=0
+#ReactionLimiter = 1000
+SurfaceTensionScaling=0.1
+VelDissipationRate = 1000
+CurvatureAlpha=0.001
+
+[domain]
+level = 0
+filename = grids/rectangle_periodic.msh
+
+[Phasefield]
+eps=0.03
+delta=0.03
+DoubleWellDelta=0.5
+Mobility=0.03
+Sigma1=1
+Sigma2=1
+Sigma3=1
+
+[Time]
+tMax = 30
+dt = 0.015
+#dtMax = 0.015
+dtMin = 0.0015
+#saveintervall = 0.5
+#dtIncreaseFactor=1.2
+#dtDecreaseFactor=0.5
+
+[Initial]
+MeanInflow=10
+
+[Refinement]
+etaRefine  = 0.2
+etaCoarsen = 0.1
+maxRefineLevel = 3
+minRefineLevel = -1
+
+[NS_Newton]
+ReassembleThreshold = 0.0
+LineSearchMaxIterations = 30
+MaxIterations = 30
+AbsoluteLimit = 1e-9
+Reduction = 1e-7
+VerbosityLevel = 1
+
+[Phasefield_Newton]
+ReassembleThreshold = 0.0
+LineSearchMaxIterations = 30
+MaxIterations = 30
+VerbosityLevel = 1
diff --git a/src/test_contactpointdynamics/cp_equation.cc b/src/test_contactpointdynamics/cp_equation.cc
new file mode 100644
index 0000000000000000000000000000000000000000..b17ed0b8d64f87cecb1664e0df1bfdd0807687d6
--- /dev/null
+++ b/src/test_contactpointdynamics/cp_equation.cc
@@ -0,0 +1,116 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include<iostream>
+// dune includes
+#include<dune/common/parallel/mpihelper.hh>
+#include<dune/common/parametertreeparser.hh>
+#include<dune/common/timer.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
+// pdelab includes
+#include<dune/pdelab/finiteelementmap/pkfem.hh>
+
+#include "3p_driver.hh"
+
+//===============================================================
+// Main program with grid setup
+//===============================================================
+
+
+
+
+int main(int argc, char** argv)
+{
+  #if !HAVE_SUPERLU
+    std::cout << "Error: These examples work only if SuperLU is available." << std::endl;
+    exit(1);
+  #endif
+
+  #if HAVE_UG
+  try{
+    // Maybe initialize Mpi
+    Dune::MPIHelper&
+      helper = Dune::MPIHelper::instance(argc, argv);
+    if(Dune::MPIHelper::isFake)
+      std::cout<< "This is a sequential program." << std::endl;
+    else
+      std::cout << "Parallel code run on "
+                << helper.size() << " process(es)" << std::endl;
+
+    // open ini file
+    Dune::ParameterTree ptree;
+    Dune::ParameterTreeParser ptreeparser;
+    ptreeparser.readINITree("chns.ini",ptree);
+    ptreeparser.readOptions(argc,argv,ptree);
+
+    // read ini file
+//    const int dim = ptree.get<int>("grid.dim");
+    const int refinement = ptree.get<int>("domain.level");
+
+
+    if (true/*dim==2*/)
+        {
+
+          typedef Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming> Grid;
+          //typedef Dune::UGGrid<2> Grid;
+          typedef Grid::ctype DF;
+
+          std::vector<int> boundary_index_map;
+          std::vector<int> element_index_map;
+
+          std::string filename = ptree.get("domain.filename",
+                                           "turbtube2d.msh");
+          Dune::GridFactory<Grid> factory;
+          if (helper.rank()==0)
+            Dune::GmshReader<Grid>::read(factory,filename,boundary_index_map,element_index_map,true,false);
+          std::shared_ptr<Grid> gridp(factory.createGrid());
+          //gridp->globalRefine(refinement);
+          //gridp->loadBalance();
+          //std::cout << " after load balance /" << helper.rank() << "/ " << gridp->size(0) << std::endl;
+
+          //typedef Grid::LeafGridView LGV;
+
+          //std::bitset<2> periodic_directions(std::string("01"));
+          //typedef Dune::XT::Grid::PeriodicGridView<LGV> GV;
+          //GV gv(gridp->leafGridView(),periodic_directions);
+
+          typedef Grid::LeafGridView GV;
+          GV gv = gridp->leafGridView();
+
+          std::cout << "OK" << std::endl;
+
+          using ES = Dune::PDELab::NonOverlappingEntitySet<GV>;
+          ES es(gv);
+
+          const int k = 2;
+          typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,k> V_FEM;
+          typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,k-1> P_FEM;
+          const int l = 1;
+          typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,l> PHI_FEM;
+
+          V_FEM vFem(es);
+          P_FEM pFem(es);
+          PHI_FEM phiFem(es);
+
+          driver_3p(*gridp, es,vFem,pFem,phiFem, ptree, helper);
+          }
+
+        }
+  catch (Dune::Exception &e){
+    std::cerr << "Dune reported error: " << e << std::endl;
+    return 1;
+  }
+  catch (...){
+    std::cerr << "Unknown exception thrown!" << std::endl;
+    return 1;
+  }
+  #endif
+}