diff --git a/dune/phasefield/porenetwork/pn_initial.hh b/dune/phasefield/porenetwork/pn_initial.hh
index be66d9ed1dfeeb6afb767b9205ef84dede559ede..f070f93c875e9ff07080c6d665fbc6fff00e1518 100644
--- a/dune/phasefield/porenetwork/pn_initial.hh
+++ b/dune/phasefield/porenetwork/pn_initial.hh
@@ -2,6 +2,8 @@
 
 #define INITIAL_SHAPE_SQUARE 0
 #define INITIAL_SHAPE_CIRCLE 1
+#define INITIAL_SHAPE_ELLIPSE 2
+#define INITIAL_SHAPE_TWOCIRCLES 3
 
 template<typename GV, typename RF, typename Params>
 class PhiInitial_Quartersquare :
@@ -41,6 +43,28 @@ public:
       RF rad = sqrt((x[0]-1)*(x[0]-1) + (x[1]-1)*(x[1]-1));
       y = param.dw.eval_shape(((1-r)-rad)*1/eps);
     }
+    else if(myshape == INITIAL_SHAPE_ELLIPSE)
+    {
+      RF rad = sqrt((x[0]-1)/1.2*(x[0]-1)/1.2 + (x[1]-1)/0.8*(x[1]-1)/0.8);
+      y = param.dw.eval_shape(((1-r)-rad)*1/eps);
+    }
+    else if(myshape == INITIAL_SHAPE_TWOCIRCLES)
+    {
+      float radius1 = 1.75; //radius circle 1
+      float radius2 = 1.75; //radius circle 2
+      float x01 = 0.9;      //xo circle 1
+      float x11 = 1.5;      //x1 circle 1
+      float x02 = 1.5;      //x0 circle 2
+      float x12 = 1.2;      //x1 circle 2
+      RF rad = sqrt((x[0]-x01)*(x[0]-x01) + (x[1]-x11)*(x[1]-x11));
+      RF y1 = param.dw.eval_shape(((1-r)-(rad)/radius1)*1/eps);
+      RF rad2 = sqrt((x[0]-x02)*(x[0]-x02) + (x[1]-x12)*(x[1]-x12));
+      RF y2 = param.dw.eval_shape(((1-r)-(rad2)/radius2)*1/eps);
+      y = y1 * y2;
+      //RF rad1 = sqrt((x[0]-2)*(x[0]-2) + (x[1]-1)*(x[1]-1));
+      //RF rad2 = sqrt((x[0]-1)*(x[0]-1) + (x[1]-2)*(x[1]-2));
+      //y = param.dw.eval_shape(((1-r)-rad2)*1/eps) * param.dw.eval_shape(((1-r)-rad1)*1/eps);
+    }
     else
     {
       std::cout << "WARNING: pn_initial.hh cannot find shape" << myshape << std::endl;
diff --git a/src/pntest/CMakeLists.txt b/src/pntest/CMakeLists.txt
index 05b01592c924f23d1f16e2d14b31b887687cc08b..67c0524d28bc37515a1e6a83030ca9477daeb300 100644
--- a/src/pntest/CMakeLists.txt
+++ b/src/pntest/CMakeLists.txt
@@ -1,3 +1,3 @@
-#add_executable(pntest pntest.cc)
-add_executable(pn2ptest pn2ptest.cc)
+add_executable(pntest pntest.cc)
+#add_executable(pn2ptest pn2ptest.cc)
 dune_symlink_to_source_files(FILES grids pn1p.ini pn2p.ini)
diff --git a/src/pntest/pntest.cc b/src/pntest/pntest.cc
index 7c239ff6ee81a17c69f473b2b79be1af573294a5..5ae546edcbc342260f11ea480d3e50aa05793190 100644
--- a/src/pntest/pntest.cc
+++ b/src/pntest/pntest.cc
@@ -3,6 +3,7 @@
 #endif
 
 #include<iostream>
+#include<fstream>
 // dune includes
 #include<dune/common/parallel/mpihelper.hh>
 #include<dune/common/parametertreeparser.hh>
@@ -55,27 +56,37 @@ int main(int argc, char** argv)
 
     //Set initial Shape
     //param.initial.shape = INITIAL_SHAPE_SQUARE;
-    param.initial.shape = INITIAL_SHAPE_CIRCLE;
+    //param.initial.shape = INITIAL_SHAPE_CIRCLE;
+    //param.initial.shape = INITIAL_SHAPE_ELLIPSE;
+    param.initial.shape = INITIAL_SHAPE_TWOCIRCLES;
 
     //Create Simulation
     Pn_1PFlow<Parameters> pn_1pFlow(GridFilename,OutputFilename,param, 10e-6); //Last Parameter is initial pore Radius in m
 
-
     //Example of usecase for Pn_1PFlow
     std::cout << "Kf = " << pn_1pFlow.getKf() << " phiSolid = " << pn_1pFlow.getPhiSolid() << std::endl;
     std::cout << "Pore Area = " << pn_1pFlow.getPoreArea() << " m^2" << std::endl;
     std::cout << "Pore Area Fraction (compared to inital values) = " << pn_1pFlow.getPoreAreaFraction() << std::endl;
-
     pn_1pFlow.writeVTK(0.0);
+    std::string filename("test.csv");
+    std::ofstream file_out;
+    file_out.open(filename, std::ios_base::app);
+    int i = 0;
+    while(1 - pn_1pFlow.getPhiSolid() > 0.15) {
+    i++;
     //pn_1pFlow.grow(2e-6);  //Growth is measured in m
     pn_1pFlow.growArea((10e-6)*(10e-6)*0.1);   // Growth of Area in m^2
-    pn_1pFlow.writeVTK(1.0);
+    //pn_1pFlow.writeVTK(1.0);
+    pn_1pFlow.writeVTK(i);
     std::cout << "Kf = " << pn_1pFlow.getKf() << " phiSolid = " << pn_1pFlow.getPhiSolid() << std::endl;
     std::cout << "Total Flux = - Kf/viscosity * gradient p" << std::endl;
     std::cout << "Kf has unit m^4" << std::endl;
 
     std::cout << "Pore Area = " << pn_1pFlow.getPoreArea()  << " m^2" << std::endl;
     std::cout << "Pore Area Fraction (compared to inital values) = " << pn_1pFlow.getPoreAreaFraction() << std::endl;
+    file_out << pn_1pFlow.getKf() << "," << 1 - pn_1pFlow.getPhiSolid() << std::endl;
+
+  }
 
   }
   catch (Dune::Exception &e){