diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 986d7a311fd1ff764814cfc56d26daa3b9ca4ccd..99814cfa14d45c1211ee887aa324873de28747bc 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -15,3 +15,4 @@ add_subdirectory(paper_thinstrip_full)
 add_subdirectory(paper_thinstrip_mixed)
 add_subdirectory(summerschool_inflow)
 add_subdirectory(pntest)
+add_subdirectory(example_2p_cahnhilliard)
diff --git a/src/example_2p_cahnhilliard/2p_cahnhilliard_equation.cc b/src/example_2p_cahnhilliard/2p_cahnhilliard_equation.cc
new file mode 100644
index 0000000000000000000000000000000000000000..109c68ed00a91145e30d08d6bde1d199530b4e21
--- /dev/null
+++ b/src/example_2p_cahnhilliard/2p_cahnhilliard_equation.cc
@@ -0,0 +1,83 @@
+#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_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)
+{
+  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("ch2p.ini",ptree);
+    ptreeparser.readOptions(argc,argv,ptree);
+
+    const int refinement = ptree.get<int>("domain.level");
+
+    typedef Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming> Grid;
+    typedef Grid::ctype DF;
+
+    std::vector<int> boundary_index_map;
+    std::vector<int> element_index_map;
+
+    std::string filename = ptree.get("domain.filename", "square.msh");
+
+    Dune::GridFactory<Grid> factory;
+    Dune::GmshReader<Grid>::read(factory,filename,true,true);
+    std::shared_ptr<Grid> gridp(factory.createGrid());
+
+    typedef Grid::LeafGridView GV;
+    GV gv = gridp->leafGridView();
+
+    using ES = Dune::PDELab::AllEntitySet<GV>;
+    ES es(gv);
+
+    const int l = 2;
+    typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,l> PHI_FEM;
+
+    PHI_FEM phiFem(es);
+
+    driver_2p(*gridp, gv, es, 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;
+  }
+}
diff --git a/src/example_2p_cahnhilliard/2p_driver.hh b/src/example_2p_cahnhilliard/2p_driver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a5437f87da3591de7260bfbcf06e2c8fd309d639
--- /dev/null
+++ b/src/example_2p_cahnhilliard/2p_driver.hh
@@ -0,0 +1,189 @@
+#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_2p_cahnhilliard.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/ff_chns.hh>
+#include <dune/phasefield/initial_fn.hh>
+#include <dune/phasefield/debug/vector_debug.hh>
+
+
+
+//_________________________________________________________________________________________________________________________
+
+template<typename Grid, typename GV, typename ES, typename PHI_FEM>
+void driver_2p(Grid& grid, GV& gv, ES& es, const PHI_FEM& phiFem, Dune::ParameterTree& ptree, Dune::MPIHelper& helper)
+{
+  const int dim = GV::dimension;
+  typedef double RF;
+
+  // We use the Paramter File for the Cahn--Hilliard Navier--Stokes Model,
+  //as it includes all relevant parameters of the Cahn--Hilliard Model
+  typedef Params_ff< RF, DoubleWell_poly<RF> ,  Mobility_Const<RF> > Parameters;
+  Parameters param(ptree);
+
+  TimestepManager timestepManager(ptree.sub("Time"),param.time,2);
+
+  //_________________________________________________________________________________________________________________________
+  // Make grid function spaces
+
+  // Again, we use the Cahn--Hilliard Navier--Stokes Model, in operator splitting form the Cahn--Hilliard grid function space
+  // is exactly the space we need.
+  typedef GFS_ff_chns<RF, ES, PHI_FEM, PHI_FEM, PHI_FEM> GFS;
+  GFS gfs(es,phiFem,phiFem,phiFem);
+  typedef typename GFS::NS NS_GFS;
+  typedef typename GFS::CH CH_GFS;
+
+
+  //_________________________________________________________________________________________________________________________
+  //Constraints
+
+  typedef Bdry_ff_chns<2, Dune::PDELab::DirichletConstraintsParameters,         // v_x, v_y
+              Dune::PDELab::NoDirichletConstraintsParameters,                   // p
+              Dune::PDELab::NoDirichletConstraintsParameters,                   // phi
+              Dune::PDELab::NoDirichletConstraintsParameters>                   // mu
+              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_ff_chns< ZeroInitial<GV,RF,Parameters>,                      //v
+                      ZeroInitial<GV,RF,Parameters>,                       //p
+                      PhiInitial_SpinodalDecomposition<GV,RF,Parameters>,  //phi
+                      ZeroInitial<GV,RF,Parameters>>  Initial;             //mu
+  Initial initial(gv,param);
+  Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
+
+  CH_V chVOld = chV;
+
+  //_________________________________________________________________________________________________________________________
+  // LocalOperators
+
+
+  typedef FEM_2P_CahnHilliard<Parameters, CH_GFS, CH_V > CH_LOP;
+  CH_LOP chLop(param, gfs.ch, chVOld);
+
+  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_ff_chns<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);
+
+
+  // Refine Intial Values
+  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> RIND0;
+        RIND0 rInd0(gv,gfs.ch,ptree.sub("Refinement"));
+        rInd0.calculateIndicator(gv,chV,false);
+        rInd0.markGrid(grid,2);
+        rInd0.refineGrid(grid,gv,chV,nsV,gfs.ch,gfs.ns,1);
+        //rInd0.myRefineGrid(grid,gv,chV,nsV,gfs.ch,gfs.ns,1);   // if used with PeriodicGridView
+        Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
+  }
+
+  gfs.updateConstraintsContainer(bdry);
+  interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
+
+
+  //_________________________________________________________________________________________________________________________
+  printVector(chV,10,"chV");
+  printVector(nsV,10,"nsV");
+
+
+  //____________________________________________________________________________________________________________________
+  chVOld = chV;
+  nsVOld = nsV;
+  while (!timestepManager.isFinished())
+    {
+
+      // Grid Operators
+      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,600,200,3,1);
+
+      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.doOnlyCHTimestep(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_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND;
+      RIND rInd(gv,gfs.ch,ptree.sub("Refinement"));
+      rInd.calculateIndicator(gv,chV,false);
+      rInd.markGrid(grid,2);
+      rInd.refineGrid(grid,chV,chVOld,gfs.ch,2);
+      //rInd.myRefineGrid(grid,gv,chV,chVOld,nsV,nsVOld,gfs.ch,gfs.ns,1);  //For PeriodicGridView
+
+      // recompute constraints
+      gfs.updateConstraintsContainer(bdry);
+      interpolateToBdry(gfs,initial.iniCH,chV);
+      interpolateToBdry(gfs,initial.iniCH,chVOld);
+    }
+}
diff --git a/src/example_2p_cahnhilliard/CMakeLists.txt b/src/example_2p_cahnhilliard/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3b66bd759c7a4918bbb69b6a48b48a5f3885cc19
--- /dev/null
+++ b/src/example_2p_cahnhilliard/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_executable(2p_cahnhilliard_equation 2p_cahnhilliard_equation.cc)
+dune_symlink_to_source_files(FILES grids ch2p.ini)
diff --git a/src/example_2p_cahnhilliard/chns.ini b/src/example_2p_cahnhilliard/chns.ini
new file mode 100644
index 0000000000000000000000000000000000000000..bcab18aa9cde7def748616cf18426941f7e0f3a7
--- /dev/null
+++ b/src/example_2p_cahnhilliard/chns.ini
@@ -0,0 +1,75 @@
+[Physics]
+uAst = 1
+D = 0.02
+rho1 = 1.0
+rho2 = 1.0
+mu=0.005
+ReactionRate=0.5
+#ReactionLimiter = 1000
+SurfaceTensionScaling=0.1
+VelDissipationRate = 1000
+#VelDissipationMax = 0.9
+CurvatureAlpha=0.001
+
+[domain]
+level = 0
+filename = grids/tper1.msh
+len = 1
+periodic = 01
+#YScalingFactor = 0.5
+YScaling = 1 #0.125
+
+[output]
+filename = output1
+
+[Phasefield]
+eps=0.03   #0.03
+delta=0.03 #0.03
+DoubleWellDelta=0.05 #0.05
+Mobility=0.04
+Sigma1=1
+Sigma2=1
+Sigma3=1
+
+[Time]
+tMax = 1
+dt = 0.0015
+dtMax = 0.01
+dtMin = 0.00015
+#saveintervall = 0.5
+#dtIncreaseFactor=1.2
+#dtDecreaseFactor=0.5
+
+[Initial]
+# FlowVyy= 10# pConst/(Domain_width * mu)
+pConst=0.05 #0.1
+uConst=0.5 #0.5
+uInflow=0.3 #0.3
+
+[RefinementInitial]
+etaRefine  = 0.2
+etaCoarsen = 0.18
+maxRefineLevel = 4
+minRefineLevel = -1
+
+[Refinement]
+etaRefine  = 0.2
+etaCoarsen = 0.1
+maxRefineLevel = 4
+minRefineLevel = -1
+
+[NS_Newton]
+ReassembleThreshold = 0.0
+LineSearchMaxIterations = 30
+MaxIterations = 20
+AbsoluteLimit = 1e-9
+Reduction = 1e-7
+VerbosityLevel = 2
+
+[Phasefield_Newton]
+ReassembleThreshold = 0.0
+LineSearchMaxIterations = 30
+MaxIterations = 100
+AbsoluteLimit = 1e-9
+Reduction = 1e-7
+VerbosityLevel = 2
diff --git a/src/example_2p_cahnhilliard/grids/periodic.dgf b/src/example_2p_cahnhilliard/grids/periodic.dgf
new file mode 100644
index 0000000000000000000000000000000000000000..b9e8e57f555a4985c5b054b7d6d33db87213c04a
--- /dev/null
+++ b/src/example_2p_cahnhilliard/grids/periodic.dgf
@@ -0,0 +1,23 @@
+DGF
+
+INTERVAL
+0 0            % lower left corner
+1 1            % upper right corner
+8 8            % number of cells in each direction
+#
+
+PERIODICFACETRANSFORMATION
+1 0, 0 1 + 1 0  % make periodic in x-direction
+1 0, 0 1 + 0 1  % make periodic in y-direction
+#
+
+BOUNDARYDOMAIN
+DEFAULT 1
+#
+
+GRIDPARAMETER
+OVERLAP 0
+#
+
+#
+ 
diff --git a/src/example_2p_cahnhilliard/grids/xperiodic.dgf b/src/example_2p_cahnhilliard/grids/xperiodic.dgf
new file mode 100644
index 0000000000000000000000000000000000000000..19762d4b0d7ed6d6db6eee27cd4d1bc4282f1cc0
--- /dev/null
+++ b/src/example_2p_cahnhilliard/grids/xperiodic.dgf
@@ -0,0 +1,22 @@
+DGF
+
+INTERVAL
+0 0            % lower left corner
+1 1            % upper right corner
+8 8            % number of cells in each direction
+#
+
+PERIODICFACETRANSFORMATION
+1 0, 0 1 + 1 0  % make periodic in x-direction
+#
+
+BOUNDARYDOMAIN
+DEFAULT 1
+#
+
+GRIDPARAMETER
+OVERLAP 0
+#
+
+#
+