diff --git a/source/Applications/DLR-F16-Porous/f16.cpp b/source/Applications/DLR-F16-Porous/f16.cpp
index 5f08df3e1f70d64b123ce3b801004a949aa09ef1..8062acf7c2ca183123f4b03c3522ee83d1af1ea1 100644
--- a/source/Applications/DLR-F16-Porous/f16.cpp
+++ b/source/Applications/DLR-F16-Porous/f16.cpp
@@ -485,6 +485,26 @@ void run(string configname)
          UBLOG(logINFO, "PID = "<<myid<<" Physical Memory currently used by current process: "<<Utilities::getPhysMemUsedByMe()/1073741824.0<<" GB");
       }
 
+      if (myid==0)
+      {
+         UBLOG(logINFO, "Parameters:");
+         UBLOG(logINFO, "* Re                  = "<<Re);
+         UBLOG(logINFO, "* Ma                  = "<<Ma);
+         UBLOG(logINFO, "* velocity (uReal)    = "<<uReal<<" m/s");
+         UBLOG(logINFO, "* viscosity (nuReal)  = "<<nuReal<<" m^2/s");
+         UBLOG(logINFO, "* chord length (lReal)= "<<lReal<<" m");
+         UBLOG(logINFO, "* velocity LB (uLB)   = "<<uLB);
+         UBLOG(logINFO, "* viscosity LB (nuLB) = "<<nuLB);
+         UBLOG(logINFO, "* chord length (l_LB) = "<<lLB<<" dx_base");
+         UBLOG(logINFO, "* dx_base             = "<<deltaXcoarse<<" m");
+         UBLOG(logINFO, "* dx_refine           = "<<deltaXfine<<" m");
+         UBLOG(logINFO, "* blocknx             = "<<blockNx[0]<<"x"<<blockNx[1]<<"x"<<blockNx[2]);
+         UBLOG(logINFO, "* refineDistance      = "<<refineDistance);
+         UBLOG(logINFO, "* number of levels    = "<<refineLevel+1);
+         UBLOG(logINFO, "* number of threads   = "<<numOfThreads);
+         UBLOG(logINFO, "* number of processes = "<<comm->getNumberOfProcesses());
+         UBLOG(logINFO, "* path = "<<pathOut);
+       }
 
       if (newStart)
       {
@@ -505,28 +525,7 @@ void run(string configname)
 
          if (myid==0)
          {
-            UBLOG(logINFO, "Parameters:");
-            UBLOG(logINFO, "* Re                  = "<<Re);
-            UBLOG(logINFO, "* Ma                  = "<<Ma);
-            UBLOG(logINFO, "* velocity (uReal)    = "<<uReal<<" m/s");
-            UBLOG(logINFO, "* viscosity (nuReal)  = "<<nuReal<<" m^2/s");
-            UBLOG(logINFO, "* chord length (lReal)= "<<lReal<<" m");
-            UBLOG(logINFO, "* velocity LB (uLB)   = "<<uLB);
-            UBLOG(logINFO, "* viscosity LB (nuLB) = "<<nuLB);
-            UBLOG(logINFO, "* chord length (l_LB) = "<<lLB<<" dx_base");
-            UBLOG(logINFO, "* dx_base             = "<<deltaXcoarse<<" m");
-            UBLOG(logINFO, "* dx_refine           = "<<deltaXfine<<" m");
-            UBLOG(logINFO, "* blocknx             = "<<blockNx[0]<<"x"<<blockNx[1]<<"x"<<blockNx[2]);
-            UBLOG(logINFO, "* refineDistance      = "<<refineDistance);
-            UBLOG(logINFO, "* number of levels    = "<<refineLevel+1);
-            UBLOG(logINFO, "* number of threads   = "<<numOfThreads);
-            UBLOG(logINFO, "* number of processes = "<<comm->getNumberOfProcesses());
-            UBLOG(logINFO, "* path = "<<pathOut);
-            UBLOG(logINFO, "Preprozess - start");
-         }
-
-         if (myid==0)
-         {
+            UBLOG(logINFO, "Preprocessing - start");
             UBLOG(logINFO, "PID = "<<myid<<" Point 3");
             UBLOG(logINFO, "PID = "<<myid<<" Physical Memory currently used by current process: "<<Utilities::getPhysMemUsedByMe()/1073741824.0<<" GB");
          }
@@ -997,6 +996,51 @@ void run(string configname)
       }
       else
       {
+         //if (myid==0) UBLOG(logINFO, "writeGeoObject");
+         //0.075 m / 0.25c
+         //SPtr<IntegrateValuesHelper> mic8(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -0.075+4.25*deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -0.3+5.25*deltaXcoarse));
+
+         //
+
+         //double dist = 0.0744;
+
+         //GbCuboid3DPtr mic0(new GbCuboid3D(0.3+deltaXfine, 0.015, 0.000517+0.00037+7.0*deltaXfine, 0.3+2.0*deltaXfine, 0.015+deltaXfine, 0.000517+0.00037+8.0*deltaXfine));
+         //if (myid==0) GbSystem3D::writeGeoObject(mic0.get(), pathOut+"/geo/mic0new", WbWriterVtkXmlBinary::getInstance());
+
+         //GbCuboid3DPtr mic1(new GbCuboid3D(0.3, 0.015, -dist-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist));
+         //if (myid==0) GbSystem3D::writeGeoObject(mic1.get(), pathOut+"/geo/mic1new", WbWriterVtkXmlBinary::getInstance());
+
+         //GbCuboid3DPtr mic2(new GbCuboid3D(0.3, 0.015, -dist*2.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*2.0));
+         //if (myid==0) GbSystem3D::writeGeoObject(mic2.get(), pathOut+"/geo/mic2new", WbWriterVtkXmlBinary::getInstance());
+
+         //GbCuboid3DPtr mic3(new GbCuboid3D(0.3, 0.015, -dist*3.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*3.0));
+         //if (myid==0) GbSystem3D::writeGeoObject(mic3.get(), pathOut+"/geo/mic3new", WbWriterVtkXmlBinary::getInstance());
+
+         //GbCuboid3DPtr mic4(new GbCuboid3D(0.3, 0.015, -dist*4.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*4.0));
+         //if (myid==0) GbSystem3D::writeGeoObject(mic4.get(), pathOut+"/geo/mic4new", WbWriterVtkXmlBinary::getInstance());
+
+         //GbCuboid3DPtr mic5(new GbCuboid3D(0.3, 0.015, -dist*5.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*5.0));
+         //if (myid==0) GbSystem3D::writeGeoObject(mic5.get(), pathOut+"/geo/mic5new", WbWriterVtkXmlBinary::getInstance());
+
+         //GbCuboid3DPtr mic6(new GbCuboid3D(0.3, 0.015, -dist*6.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*6.0));
+         //if (myid==0) GbSystem3D::writeGeoObject(mic6.get(), pathOut+"/geo/mic6new", WbWriterVtkXmlBinary::getInstance());
+
+         //GbCuboid3DPtr mic7(new GbCuboid3D(0.3, 0.015, -dist*7.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*7.0));
+         //if (myid==0) GbSystem3D::writeGeoObject(mic7.get(), pathOut+"/geo/mic7new", WbWriterVtkXmlBinary::getInstance());
+
+         //GbCuboid3DPtr mic8(new GbCuboid3D(0.3, 0.015, -dist*8.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*8.0));
+         //if (myid==0) GbSystem3D::writeGeoObject(mic8.get(), pathOut+"/geo/mic8new", WbWriterVtkXmlBinary::getInstance());
+
+         //GbCuboid3DPtr geoAddWallP(new GbCuboid3D(0.269984375, g_minX2-blockLength, 0.0016, 0.2702421875, g_maxX2+blockLength, 0.0076));
+         //if (myid==0) GbSystem3D::writeGeoObject(geoAddWallP.get(), pathOut+"/geo/geoAddWallP", WbWriterVtkXmlASCII::getInstance());
+         //SPtr<D3Q27Interactor> addWallPIntr = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoAddWallP, grid, noSlipBCAdapter, Interactor3D::SOLID));
+
+         //SetBcBlocksBlockVisitor setBcVisitor(addWallPIntr);
+         //grid->accept(setBcVisitor);
+         //addWallPIntr->initInteractor();
+         //return;
+
+
          //restartCoProcessor->restart((int)restartStep);
          migCoProcessor->restart((int)restartStep);
          grid->setTimeStep(restartStep);
@@ -1107,48 +1151,99 @@ void run(string configname)
          TimeAveragedValuesCoProcessor::Density | TimeAveragedValuesCoProcessor::Velocity | TimeAveragedValuesCoProcessor::Fluctuations));
       tav->setWithGhostLayer(true);
 
-      SPtr<IntegrateValuesHelper> mic1(new IntegrateValuesHelper(grid, comm, 0.3-deltaXfine, 0.015, 0.0005, 0.3, 0.015+deltaXfine, 0.0005+deltaXfine));
-      if (myid==0) GbSystem3D::writeGeoObject(mic1->getBoundingBox().get(), pathOut+"/geo/mic1", WbWriterVtkXmlBinary::getInstance());
-      SPtr<UbScheduler> stepMV(new UbScheduler(1, 0, 1000000));
-      SPtr<TimeseriesCoProcessor> tsp1(new TimeseriesCoProcessor(grid, stepMV, mic1, pathOut+"/mic/mic1", comm));
+      //SPtr<IntegrateValuesHelper> mic1(new IntegrateValuesHelper(grid, comm, 0.3-deltaXfine, 0.015, 0.0005, 0.3, 0.015+deltaXfine, 0.0005+deltaXfine));
+      //if (myid==0) GbSystem3D::writeGeoObject(mic1->getBoundingBox().get(), pathOut+"/geo/mic1", WbWriterVtkXmlBinary::getInstance());
+      //SPtr<UbScheduler> stepMV(new UbScheduler(1, 0, 1000000));
+      //SPtr<TimeseriesCoProcessor> tsp1(new TimeseriesCoProcessor(grid, stepMV, mic1, pathOut+"/mic/mic1", comm));
+
+      //SPtr<IntegrateValuesHelper> mic2(new IntegrateValuesHelper(grid, comm, 0.3+deltaXfine, 0.015, 0.001685, 0.3, 0.015+deltaXfine, 0.001685+deltaXfine));
+      //if (myid==0) GbSystem3D::writeGeoObject(mic2->getBoundingBox().get(), pathOut+"/geo/mic2", WbWriterVtkXmlBinary::getInstance());
+      //SPtr<TimeseriesCoProcessor> tsp2(new TimeseriesCoProcessor(grid, stepMV, mic2, pathOut+"/mic/mic2", comm));
 
-      SPtr<IntegrateValuesHelper> mic2(new IntegrateValuesHelper(grid, comm, 0.3+deltaXfine, 0.015, 0.001685, 0.3, 0.015+deltaXfine, 0.001685+deltaXfine));
-      if (myid==0) GbSystem3D::writeGeoObject(mic2->getBoundingBox().get(), pathOut+"/geo/mic2", WbWriterVtkXmlBinary::getInstance());
-      SPtr<TimeseriesCoProcessor> tsp2(new TimeseriesCoProcessor(grid, stepMV, mic2, pathOut+"/mic/mic2", comm));
+      //SPtr<IntegrateValuesHelper> mic3(new IntegrateValuesHelper(grid, comm, 0.3-deltaXcoarse, 0.015, -0.46+4.25*deltaXcoarse, 0.3, 0.015+deltaXcoarse, -0.46+5.25*deltaXcoarse));
+      //if (myid==0) GbSystem3D::writeGeoObject(mic3->getBoundingBox().get(), pathOut+"/geo/mic3", WbWriterVtkXmlBinary::getInstance());
+      //SPtr<TimeseriesCoProcessor> tsp3(new TimeseriesCoProcessor(grid, stepMV, mic3, pathOut+"/mic/mic3", comm));
 
-      SPtr<IntegrateValuesHelper> mic3(new IntegrateValuesHelper(grid, comm, 0.3-deltaXcoarse, 0.015, -0.46+4.25*deltaXcoarse, 0.3, 0.015+deltaXcoarse, -0.46+5.25*deltaXcoarse));
-      if (myid==0) GbSystem3D::writeGeoObject(mic3->getBoundingBox().get(), pathOut+"/geo/mic3", WbWriterVtkXmlBinary::getInstance());
-      SPtr<TimeseriesCoProcessor> tsp3(new TimeseriesCoProcessor(grid, stepMV, mic3, pathOut+"/mic/mic3", comm));
+      //SPtr<IntegrateValuesHelper> mic4(new IntegrateValuesHelper(grid, comm, 0.3-deltaXcoarse, 0.015, 0.46-5.25*deltaXcoarse, 0.3, 0.015+deltaXcoarse, 0.46-4.25*deltaXcoarse));
+      //if (myid==0) GbSystem3D::writeGeoObject(mic4->getBoundingBox().get(), pathOut+"/geo/mic4", WbWriterVtkXmlBinary::getInstance());
+      //SPtr<TimeseriesCoProcessor> tsp4(new TimeseriesCoProcessor(grid, stepMV, mic4, pathOut+"/mic/mic4", comm));
 
-      SPtr<IntegrateValuesHelper> mic4(new IntegrateValuesHelper(grid, comm, 0.3-deltaXcoarse, 0.015, 0.46-5.25*deltaXcoarse, 0.3, 0.015+deltaXcoarse, 0.46-4.25*deltaXcoarse));
-      if (myid==0) GbSystem3D::writeGeoObject(mic4->getBoundingBox().get(), pathOut+"/geo/mic4", WbWriterVtkXmlBinary::getInstance());
-      SPtr<TimeseriesCoProcessor> tsp4(new TimeseriesCoProcessor(grid, stepMV, mic4, pathOut+"/mic/mic4", comm));
+      //SPtr<IntegrateValuesHelper> mic5(new IntegrateValuesHelper(grid, comm, 0.3+deltaXfine, 0.015, 0.000517+0.00037+7.0*deltaXfine, 0.3+2.0*deltaXfine, 0.015+deltaXfine, 0.000517+0.00037+8.0*deltaXfine));
+      //if (myid==0) GbSystem3D::writeGeoObject(mic5->getBoundingBox().get(), pathOut+"/geo/mic5", WbWriterVtkXmlBinary::getInstance());
+      //SPtr<TimeseriesCoProcessor> tsp5(new TimeseriesCoProcessor(grid, stepMV, mic5, pathOut+"/mic/mic5", comm));
 
-      SPtr<IntegrateValuesHelper> mic5(new IntegrateValuesHelper(grid, comm, 0.3+deltaXfine, 0.015, 0.000517+0.00037+7.0*deltaXfine, 0.3+2.0*deltaXfine, 0.015+deltaXfine, 0.000517+0.00037+8.0*deltaXfine));
-      if (myid==0) GbSystem3D::writeGeoObject(mic5->getBoundingBox().get(), pathOut+"/geo/mic5", WbWriterVtkXmlBinary::getInstance());
-      SPtr<TimeseriesCoProcessor> tsp5(new TimeseriesCoProcessor(grid, stepMV, mic5, pathOut+"/mic/mic5", comm));
+      ////0.46 m / 1.5c
+      //SPtr<IntegrateValuesHelper> mic6(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -0.4599-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -0.4599));
+      //if (myid==0) GbSystem3D::writeGeoObject(mic6->getBoundingBox().get(), pathOut+"/geo/mic6", WbWriterVtkXmlBinary::getInstance());
+      //SPtr<TimeseriesCoProcessor> tsp6(new TimeseriesCoProcessor(grid, stepMV, mic6, pathOut+"/mic/mic6", comm));
 
-      SPtr<IntegrateValuesHelper> mic6(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -0.46+4.25*deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -0.46+5.25*deltaXcoarse));
-      if (myid==0) GbSystem3D::writeGeoObject(mic6->getBoundingBox().get(), pathOut+"/geo/mic6", WbWriterVtkXmlBinary::getInstance());
-      SPtr<TimeseriesCoProcessor> tsp6(new TimeseriesCoProcessor(grid, stepMV, mic6, pathOut+"/mic/mic6", comm));
+      ////0.3 m / 1.0c
+      //SPtr<IntegrateValuesHelper> mic7(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -0.299, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -0.299+deltaXcoarse));
+      //if (myid==0) GbSystem3D::writeGeoObject(mic7->getBoundingBox().get(), pathOut+"/geo/mic7", WbWriterVtkXmlBinary::getInstance());
+      //SPtr<TimeseriesCoProcessor> tsp7(new TimeseriesCoProcessor(grid, stepMV, mic7, pathOut+"/mic/mic7", comm));
+
+      ////0.075 m / 0.25c
+      //SPtr<IntegrateValuesHelper> mic8(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -0.0744-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -0.0744));
+      //if (myid==0) GbSystem3D::writeGeoObject(mic8->getBoundingBox().get(), pathOut+"/geo/mic8", WbWriterVtkXmlBinary::getInstance());
+      //SPtr<TimeseriesCoProcessor> tsp8(new TimeseriesCoProcessor(grid, stepMV, mic8, pathOut+"/mic/mic8", comm));
+
+      double dist = 0.0744; //0.25c
+      SPtr<UbScheduler> stepMV(new UbScheduler(1, 0, 1000000));
+      //0.0
+      SPtr<IntegrateValuesHelper> mic0new(new IntegrateValuesHelper(grid, comm, 0.3+deltaXfine, 0.015, 0.000517+0.00037+7.0*deltaXfine, 0.3+2.0*deltaXfine, 0.015+deltaXfine, 0.000517+0.00037+8.0*deltaXfine));
+      if (myid==0) GbSystem3D::writeGeoObject(mic0new->getBoundingBox().get(), pathOut+"/geo/mic0new", WbWriterVtkXmlBinary::getInstance());
+      SPtr<TimeseriesCoProcessor> tsp0(new TimeseriesCoProcessor(grid, stepMV, mic0new, pathOut+"/mic/mic0new", comm));
+      //0.25c
+      SPtr<IntegrateValuesHelper> mic1new(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -dist-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist));
+      if (myid==0) GbSystem3D::writeGeoObject(mic1new->getBoundingBox().get(), pathOut+"/geo/mic1new", WbWriterVtkXmlBinary::getInstance());
+      SPtr<TimeseriesCoProcessor> tsp1(new TimeseriesCoProcessor(grid, stepMV, mic1new, pathOut+"/mic/mic1new", comm));
+      //0.5c
+      SPtr<IntegrateValuesHelper> mic2new(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -dist*2.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*2.0));
+      if (myid==0) GbSystem3D::writeGeoObject(mic2new->getBoundingBox().get(), pathOut+"/geo/mic2new", WbWriterVtkXmlBinary::getInstance());
+      SPtr<TimeseriesCoProcessor> tsp2(new TimeseriesCoProcessor(grid, stepMV, mic2new, pathOut+"/mic/mic2new", comm));
+      //0.75c
+      SPtr<IntegrateValuesHelper> mic3new(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -dist*3.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*3.0));
+      if (myid==0) GbSystem3D::writeGeoObject(mic3new->getBoundingBox().get(), pathOut+"/geo/mic3new", WbWriterVtkXmlBinary::getInstance());
+      SPtr<TimeseriesCoProcessor> tsp3(new TimeseriesCoProcessor(grid, stepMV, mic3new, pathOut+"/mic/mic3new", comm));
+      //1.0c
+      SPtr<IntegrateValuesHelper> mic4new(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -dist*4.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*4.0));
+      if (myid==0) GbSystem3D::writeGeoObject(mic4new->getBoundingBox().get(), pathOut+"/geo/mic4new", WbWriterVtkXmlBinary::getInstance());
+      SPtr<TimeseriesCoProcessor> tsp4(new TimeseriesCoProcessor(grid, stepMV, mic4new, pathOut+"/mic/mic4new", comm));
+      //1.25c
+      SPtr<IntegrateValuesHelper> mic5new(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -dist*5.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*5.0));
+      if (myid==0) GbSystem3D::writeGeoObject(mic5new->getBoundingBox().get(), pathOut+"/geo/mic5new", WbWriterVtkXmlBinary::getInstance());
+      SPtr<TimeseriesCoProcessor> tsp5(new TimeseriesCoProcessor(grid, stepMV, mic5new, pathOut+"/mic/mic5new", comm));
+      //1.5c
+      SPtr<IntegrateValuesHelper> mic6new(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -dist*6.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*6.0));
+      if (myid==0) GbSystem3D::writeGeoObject(mic6new->getBoundingBox().get(), pathOut+"/geo/mic6new", WbWriterVtkXmlBinary::getInstance());
+      SPtr<TimeseriesCoProcessor> tsp6(new TimeseriesCoProcessor(grid, stepMV, mic6new, pathOut+"/mic/mic6new", comm));
+      //1.75c
+      SPtr<IntegrateValuesHelper> mic7new(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -dist*7.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*7.0));
+      if (myid==0) GbSystem3D::writeGeoObject(mic7new->getBoundingBox().get(), pathOut+"/geo/mic7new", WbWriterVtkXmlBinary::getInstance());
+      SPtr<TimeseriesCoProcessor> tsp7(new TimeseriesCoProcessor(grid, stepMV, mic7new, pathOut+"/mic/mic7new", comm));
+      //2.0c
+      SPtr<IntegrateValuesHelper> mic8new(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -dist*8.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*8.0));
+      if (myid==0) GbSystem3D::writeGeoObject(mic8new->getBoundingBox().get(), pathOut+"/geo/mic8new", WbWriterVtkXmlBinary::getInstance());
+      SPtr<TimeseriesCoProcessor> tsp8(new TimeseriesCoProcessor(grid, stepMV, mic8new, pathOut+"/mic/mic8new", comm));
 
       omp_set_num_threads(numOfThreads);
       SPtr<UbScheduler> stepGhostLayer(new UbScheduler(1));
       SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endTime));
       calculator->addCoProcessor(nupsCoProcessor);
-      //calculator->addCoProcessor(restartCoProcessor);
-      calculator->addCoProcessor(migCoProcessor);
-      calculator->addCoProcessor(writeMQSelectCoProcessor);
-      calculator->addCoProcessor(writeMQCoProcessor);
+      calculator->addCoProcessor(tsp0);
       calculator->addCoProcessor(tsp1);
       calculator->addCoProcessor(tsp2);
       calculator->addCoProcessor(tsp3);
       calculator->addCoProcessor(tsp4);
       calculator->addCoProcessor(tsp5);
       calculator->addCoProcessor(tsp6);
+      calculator->addCoProcessor(tsp7);
+      calculator->addCoProcessor(tsp8);
+      calculator->addCoProcessor(restartCoProcessor);
+      calculator->addCoProcessor(writeMQSelectCoProcessor);
+      calculator->addCoProcessor(writeMQCoProcessor);
       calculator->addCoProcessor(tav);
 
-
       if (myid==0) UBLOG(logINFO, "Simulation-start");
       calculator->calculate();
       if (myid==0) UBLOG(logINFO, "Simulation-end");
diff --git a/source/Applications/DLR-F16-Solid/f16.cpp b/source/Applications/DLR-F16-Solid/f16.cpp
index 27b7bc62350c7568f95092d778efeaafbc01e7dc..177fa424682675f6967a063cad1981cb0c98bbeb 100644
--- a/source/Applications/DLR-F16-Solid/f16.cpp
+++ b/source/Applications/DLR-F16-Solid/f16.cpp
@@ -196,6 +196,26 @@ void run(string configname)
          UBLOG(logINFO, "PID = "<<myid<<" Physical Memory currently used by current process: "<<Utilities::getPhysMemUsedByMe()/1073741824.0<<" GB");
       }
 
+      if (myid==0)
+      {
+         UBLOG(logINFO, "Parameters:");
+         UBLOG(logINFO, "* Re                  = "<<Re);
+         UBLOG(logINFO, "* Ma                  = "<<Ma);
+         UBLOG(logINFO, "* velocity (uReal)    = "<<uReal<<" m/s");
+         UBLOG(logINFO, "* viscosity (nuReal)  = "<<nuReal<<" m^2/s");
+         UBLOG(logINFO, "* chord length (lReal)= "<<lReal<<" m");
+         UBLOG(logINFO, "* velocity LB (uLB)   = "<<uLB);
+         UBLOG(logINFO, "* viscosity LB (nuLB) = "<<nuLB);
+         UBLOG(logINFO, "* chord length (l_LB) = "<<lLB<<" dx_base");
+         UBLOG(logINFO, "* dx_base             = "<<deltaXcoarse<<" m");
+         UBLOG(logINFO, "* dx_refine           = "<<deltaXfine<<" m");
+         UBLOG(logINFO, "* blocknx             = "<<blockNx[0]<<"x"<<blockNx[1]<<"x"<<blockNx[2]);
+         UBLOG(logINFO, "* refineDistance      = "<<refineDistance);
+         UBLOG(logINFO, "* number of levels    = "<<refineLevel+1);
+         UBLOG(logINFO, "* number of threads   = "<<numOfThreads);
+         UBLOG(logINFO, "* number of processes = "<<comm->getNumberOfProcesses());
+         UBLOG(logINFO, "* path = "<<pathOut);
+      }
 
       if (newStart)
       {
@@ -216,24 +236,7 @@ void run(string configname)
 
          if (myid==0)
          {
-            UBLOG(logINFO, "Parameters:");
-            UBLOG(logINFO, "* Re                  = "<<Re);
-            UBLOG(logINFO, "* Ma                  = "<<Ma);
-            UBLOG(logINFO, "* velocity (uReal)    = "<<uReal<<" m/s");
-            UBLOG(logINFO, "* viscosity (nuReal)  = "<<nuReal<<" m^2/s");
-            UBLOG(logINFO, "* chord length (lReal)= "<<lReal<<" m");
-            UBLOG(logINFO, "* velocity LB (uLB)   = "<<uLB);
-            UBLOG(logINFO, "* viscosity LB (nuLB) = "<<nuLB);
-            UBLOG(logINFO, "* chord length (l_LB) = "<<lLB<<" dx_base");
-            UBLOG(logINFO, "* dx_base             = "<<deltaXcoarse<<" m");
-            UBLOG(logINFO, "* dx_refine           = "<<deltaXfine<<" m");
-            UBLOG(logINFO, "* blocknx             = "<<blockNx[0]<<"x"<<blockNx[1]<<"x"<<blockNx[2]);
-            UBLOG(logINFO, "* refineDistance      = "<<refineDistance);
-            UBLOG(logINFO, "* number of levels    = "<<refineLevel+1);
-            UBLOG(logINFO, "* number of threads   = "<<numOfThreads);
-            UBLOG(logINFO, "* number of processes = "<<comm->getNumberOfProcesses());
-            UBLOG(logINFO, "* path = "<<pathOut);
-            UBLOG(logINFO, "Preprozess - start");
+            UBLOG(logINFO, "Preprocessing - start");
          }
 
          
@@ -714,49 +717,98 @@ void run(string configname)
          TimeAveragedValuesCoProcessor::Density | TimeAveragedValuesCoProcessor::Velocity | TimeAveragedValuesCoProcessor::Fluctuations));
       tav->setWithGhostLayer(true);
 
-      SPtr<IntegrateValuesHelper> mic1(new IntegrateValuesHelper(grid, comm, 0.3-deltaXfine, 0.015, 0.0005, 0.3, 0.015+deltaXfine, 0.0005+deltaXfine));
-      if (myid==0) GbSystem3D::writeGeoObject(mic1->getBoundingBox().get(), pathOut+"/geo/mic1", WbWriterVtkXmlBinary::getInstance());
-      SPtr<UbScheduler> stepMV(new UbScheduler(1, 0, 1000000));
-      SPtr<TimeseriesCoProcessor> tsp1(new TimeseriesCoProcessor(grid, stepMV, mic1, pathOut+"/mic/mic1", comm));
+      //SPtr<IntegrateValuesHelper> mic1(new IntegrateValuesHelper(grid, comm, 0.3-deltaXfine, 0.015, 0.0005, 0.3, 0.015+deltaXfine, 0.0005+deltaXfine));
+      //if (myid==0) GbSystem3D::writeGeoObject(mic1->getBoundingBox().get(), pathOut+"/geo/mic1", WbWriterVtkXmlBinary::getInstance());
+      //SPtr<UbScheduler> stepMV(new UbScheduler(1, 0, 1000000));
+      //SPtr<TimeseriesCoProcessor> tsp1(new TimeseriesCoProcessor(grid, stepMV, mic1, pathOut+"/mic/mic1", comm));
 
-      SPtr<IntegrateValuesHelper> mic2(new IntegrateValuesHelper(grid, comm, 0.3+deltaXfine, 0.015, 0.001685, 0.3, 0.015+deltaXfine, 0.001685+deltaXfine));
-      if (myid==0) GbSystem3D::writeGeoObject(mic2->getBoundingBox().get(), pathOut+"/geo/mic2", WbWriterVtkXmlBinary::getInstance());
-      SPtr<TimeseriesCoProcessor> tsp2(new TimeseriesCoProcessor(grid, stepMV, mic2, pathOut+"/mic/mic2", comm));
+      //SPtr<IntegrateValuesHelper> mic2(new IntegrateValuesHelper(grid, comm, 0.3+deltaXfine, 0.015, 0.001685, 0.3, 0.015+deltaXfine, 0.001685+deltaXfine));
+      //if (myid==0) GbSystem3D::writeGeoObject(mic2->getBoundingBox().get(), pathOut+"/geo/mic2", WbWriterVtkXmlBinary::getInstance());
+      //SPtr<TimeseriesCoProcessor> tsp2(new TimeseriesCoProcessor(grid, stepMV, mic2, pathOut+"/mic/mic2", comm));
 
-      SPtr<IntegrateValuesHelper> mic3(new IntegrateValuesHelper(grid, comm, 0.3-deltaXcoarse, 0.015, -0.46+4.25*deltaXcoarse, 0.3, 0.015+deltaXcoarse, -0.46+5.25*deltaXcoarse));
-      if (myid==0) GbSystem3D::writeGeoObject(mic3->getBoundingBox().get(), pathOut+"/geo/mic3", WbWriterVtkXmlBinary::getInstance());
-      SPtr<TimeseriesCoProcessor> tsp3(new TimeseriesCoProcessor(grid, stepMV, mic3, pathOut+"/mic/mic3", comm));
+      //SPtr<IntegrateValuesHelper> mic3(new IntegrateValuesHelper(grid, comm, 0.3-deltaXcoarse, 0.015, -0.46+4.25*deltaXcoarse, 0.3, 0.015+deltaXcoarse, -0.46+5.25*deltaXcoarse));
+      //if (myid==0) GbSystem3D::writeGeoObject(mic3->getBoundingBox().get(), pathOut+"/geo/mic3", WbWriterVtkXmlBinary::getInstance());
+      //SPtr<TimeseriesCoProcessor> tsp3(new TimeseriesCoProcessor(grid, stepMV, mic3, pathOut+"/mic/mic3", comm));
 
-      SPtr<IntegrateValuesHelper> mic4(new IntegrateValuesHelper(grid, comm, 0.3-deltaXcoarse, 0.015, 0.46-5.25*deltaXcoarse, 0.3, 0.015+deltaXcoarse, 0.46-4.25*deltaXcoarse));
-      if (myid==0) GbSystem3D::writeGeoObject(mic4->getBoundingBox().get(), pathOut+"/geo/mic4", WbWriterVtkXmlBinary::getInstance());
-      SPtr<TimeseriesCoProcessor> tsp4(new TimeseriesCoProcessor(grid, stepMV, mic4, pathOut+"/mic/mic4", comm));
+      //SPtr<IntegrateValuesHelper> mic4(new IntegrateValuesHelper(grid, comm, 0.3-deltaXcoarse, 0.015, 0.46-5.25*deltaXcoarse, 0.3, 0.015+deltaXcoarse, 0.46-4.25*deltaXcoarse));
+      //if (myid==0) GbSystem3D::writeGeoObject(mic4->getBoundingBox().get(), pathOut+"/geo/mic4", WbWriterVtkXmlBinary::getInstance());
+      //SPtr<TimeseriesCoProcessor> tsp4(new TimeseriesCoProcessor(grid, stepMV, mic4, pathOut+"/mic/mic4", comm));
 
-      SPtr<IntegrateValuesHelper> mic5(new IntegrateValuesHelper(grid, comm, 0.3+deltaXfine, 0.015, 0.000517+0.00037+7.0*deltaXfine, 0.3+2.0*deltaXfine, 0.015+deltaXfine, 0.000517+0.00037+8.0*deltaXfine));
-      if (myid==0) GbSystem3D::writeGeoObject(mic5->getBoundingBox().get(), pathOut+"/geo/mic5", WbWriterVtkXmlBinary::getInstance());
-      SPtr<TimeseriesCoProcessor> tsp5(new TimeseriesCoProcessor(grid, stepMV, mic5, pathOut+"/mic/mic5", comm));
+      //SPtr<IntegrateValuesHelper> mic5(new IntegrateValuesHelper(grid, comm, 0.3+deltaXfine, 0.015, 0.000517+0.00037+7.0*deltaXfine, 0.3+2.0*deltaXfine, 0.015+deltaXfine, 0.000517+0.00037+8.0*deltaXfine));
+      //if (myid==0) GbSystem3D::writeGeoObject(mic5->getBoundingBox().get(), pathOut+"/geo/mic5", WbWriterVtkXmlBinary::getInstance());
+      //SPtr<TimeseriesCoProcessor> tsp5(new TimeseriesCoProcessor(grid, stepMV, mic5, pathOut+"/mic/mic5", comm));
 
-      //SPtr<IntegrateValuesHelper> mic6(new IntegrateValuesHelper(grid, comm, 0.3, 0.015,  -0.46+4.25*deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -0.46+5.25*deltaXcoarse));
+      ////0.46 m / 1.5c
+      //SPtr<IntegrateValuesHelper> mic6(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -0.4599-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -0.4599));
       //if (myid==0) GbSystem3D::writeGeoObject(mic6->getBoundingBox().get(), pathOut+"/geo/mic6", WbWriterVtkXmlBinary::getInstance());
       //SPtr<TimeseriesCoProcessor> tsp6(new TimeseriesCoProcessor(grid, stepMV, mic6, pathOut+"/mic/mic6", comm));
 
-      //SPtr<IntegrateValuesHelper> mic7(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -0.3+4.25*deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -0.3+5.25*deltaXcoarse));
+      ////0.3 m / 1.0c
+      //SPtr<IntegrateValuesHelper> mic7(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -0.299, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -0.299+deltaXcoarse));
       //if (myid==0) GbSystem3D::writeGeoObject(mic7->getBoundingBox().get(), pathOut+"/geo/mic7", WbWriterVtkXmlBinary::getInstance());
       //SPtr<TimeseriesCoProcessor> tsp7(new TimeseriesCoProcessor(grid, stepMV, mic7, pathOut+"/mic/mic7", comm));
 
+      ////0.075 m / 0.25c
+      //SPtr<IntegrateValuesHelper> mic8(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -0.0744-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -0.0744));
+      //if (myid==0) GbSystem3D::writeGeoObject(mic8->getBoundingBox().get(), pathOut+"/geo/mic8", WbWriterVtkXmlBinary::getInstance());
+      //SPtr<TimeseriesCoProcessor> tsp8(new TimeseriesCoProcessor(grid, stepMV, mic8, pathOut+"/mic/mic8", comm));
+
+      double dist = 0.0744; //0.25c
+      SPtr<UbScheduler> stepMV(new UbScheduler(1, 0, 1000000));
+      //0.0
+      SPtr<IntegrateValuesHelper> mic0new(new IntegrateValuesHelper(grid, comm, 0.3+deltaXfine, 0.015, 0.000517+0.00037+7.0*deltaXfine, 0.3+2.0*deltaXfine, 0.015+deltaXfine, 0.000517+0.00037+8.0*deltaXfine));
+      if (myid==0) GbSystem3D::writeGeoObject(mic0new->getBoundingBox().get(), pathOut+"/geo/mic0new", WbWriterVtkXmlBinary::getInstance());
+      SPtr<TimeseriesCoProcessor> tsp0(new TimeseriesCoProcessor(grid, stepMV, mic0new, pathOut+"/mic/mic0new", comm));
+      //0.25c
+      SPtr<IntegrateValuesHelper> mic1new(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -dist-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist));
+      if (myid==0) GbSystem3D::writeGeoObject(mic1new->getBoundingBox().get(), pathOut+"/geo/mic1new", WbWriterVtkXmlBinary::getInstance());
+      SPtr<TimeseriesCoProcessor> tsp1(new TimeseriesCoProcessor(grid, stepMV, mic1new, pathOut+"/mic/mic1new", comm));
+      //0.5c
+      SPtr<IntegrateValuesHelper> mic2new(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -dist*2.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*2.0));
+      if (myid==0) GbSystem3D::writeGeoObject(mic2new->getBoundingBox().get(), pathOut+"/geo/mic2new", WbWriterVtkXmlBinary::getInstance());
+      SPtr<TimeseriesCoProcessor> tsp2(new TimeseriesCoProcessor(grid, stepMV, mic2new, pathOut+"/mic/mic2new", comm));
+      //0.75c
+      SPtr<IntegrateValuesHelper> mic3new(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -dist*3.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*3.0));
+      if (myid==0) GbSystem3D::writeGeoObject(mic3new->getBoundingBox().get(), pathOut+"/geo/mic3new", WbWriterVtkXmlBinary::getInstance());
+      SPtr<TimeseriesCoProcessor> tsp3(new TimeseriesCoProcessor(grid, stepMV, mic3new, pathOut+"/mic/mic3new", comm));
+      //1.0c
+      SPtr<IntegrateValuesHelper> mic4new(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -dist*4.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*4.0));
+      if (myid==0) GbSystem3D::writeGeoObject(mic4new->getBoundingBox().get(), pathOut+"/geo/mic4new", WbWriterVtkXmlBinary::getInstance());
+      SPtr<TimeseriesCoProcessor> tsp4(new TimeseriesCoProcessor(grid, stepMV, mic4new, pathOut+"/mic/mic4new", comm));
+      //1.25c
+      SPtr<IntegrateValuesHelper> mic5new(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -dist*5.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*5.0));
+      if (myid==0) GbSystem3D::writeGeoObject(mic5new->getBoundingBox().get(), pathOut+"/geo/mic5new", WbWriterVtkXmlBinary::getInstance());
+      SPtr<TimeseriesCoProcessor> tsp5(new TimeseriesCoProcessor(grid, stepMV, mic5new, pathOut+"/mic/mic5new", comm));
+      //1.5c
+      SPtr<IntegrateValuesHelper> mic6new(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -dist*6.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*6.0));
+      if (myid==0) GbSystem3D::writeGeoObject(mic6new->getBoundingBox().get(), pathOut+"/geo/mic6new", WbWriterVtkXmlBinary::getInstance());
+      SPtr<TimeseriesCoProcessor> tsp6(new TimeseriesCoProcessor(grid, stepMV, mic6new, pathOut+"/mic/mic6new", comm));
+      //1.75c
+      SPtr<IntegrateValuesHelper> mic7new(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -dist*7.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*7.0));
+      if (myid==0) GbSystem3D::writeGeoObject(mic7new->getBoundingBox().get(), pathOut+"/geo/mic7new", WbWriterVtkXmlBinary::getInstance());
+      SPtr<TimeseriesCoProcessor> tsp7(new TimeseriesCoProcessor(grid, stepMV, mic7new, pathOut+"/mic/mic7new", comm));
+      //2.0c
+      SPtr<IntegrateValuesHelper> mic8new(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -dist*8.0-deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -dist*8.0));
+      if (myid==0) GbSystem3D::writeGeoObject(mic8new->getBoundingBox().get(), pathOut+"/geo/mic8new", WbWriterVtkXmlBinary::getInstance());
+      SPtr<TimeseriesCoProcessor> tsp8(new TimeseriesCoProcessor(grid, stepMV, mic8new, pathOut+"/mic/mic8new", comm));
+
       omp_set_num_threads(numOfThreads);
       SPtr<UbScheduler> stepGhostLayer(new UbScheduler(1));
       SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endTime));
       calculator->addCoProcessor(nupsCoProcessor);
+      calculator->addCoProcessor(tsp0);
+      calculator->addCoProcessor(tsp1);
+      calculator->addCoProcessor(tsp2);
+      calculator->addCoProcessor(tsp3);
+      calculator->addCoProcessor(tsp4);
+      calculator->addCoProcessor(tsp5);
+      calculator->addCoProcessor(tsp6);
+      calculator->addCoProcessor(tsp7);
+      calculator->addCoProcessor(tsp8);
       calculator->addCoProcessor(restartCoProcessor);
-      //calculator->addCoProcessor(writeMQSelectCoProcessor);
+      calculator->addCoProcessor(writeMQSelectCoProcessor);
       calculator->addCoProcessor(writeMQCoProcessor);
-      //calculator->addCoProcessor(tsp1);
-      //calculator->addCoProcessor(tsp2);
-      //calculator->addCoProcessor(tsp3);
-      //calculator->addCoProcessor(tsp4);
-      //calculator->addCoProcessor(tsp5);
-      //calculator->addCoProcessor(tsp6);
-      //calculator->addCoProcessor(tav);
+      calculator->addCoProcessor(tav);
 
 
       if (myid==0) UBLOG(logINFO, "Simulation-start");
diff --git a/source/Applications/pChannel/configBombadilpChannel.cfg b/source/Applications/pChannel/configBombadilpChannel.cfg
index c94f7669ad153820f6a88bdcf841446dc2d90d66..c9cda307ad4ef83f94633c9a51a50914fe9ebca2 100644
--- a/source/Applications/pChannel/configBombadilpChannel.cfg
+++ b/source/Applications/pChannel/configBombadilpChannel.cfg
@@ -58,8 +58,10 @@ pmL = 1e-3 1e-3 1e-3
 #grid
 refineLevel = 0
 #deltaXfine  = 10e-6
-deltaXfine  = 20e-6
-#deltaXfine  = 80e-6 #level 2
+
+#deltaXfine  = 20e-6
+deltaXfine  = 80e-6 #level 2
+
 
 blocknx = 10 10 10
 
diff --git a/source/Applications/pChannel/pChannel.cpp b/source/Applications/pChannel/pChannel.cpp
index e7fa6abcf13f93748f0619eabacd5658fab55e8e..f698d8d26bc2963464a0b30f709fbd9737e5fa96 100644
--- a/source/Applications/pChannel/pChannel.cpp
+++ b/source/Applications/pChannel/pChannel.cpp
@@ -23,6 +23,147 @@ double rangeRandom1()
 
 using namespace std;
 
+std::vector<int> x1Nbr;
+std::vector<int> x2Nbr;
+std::vector<int> x3Nbr;
+
+std::vector<int> x1NbrTemp;
+std::vector<int> x2NbrTemp;
+std::vector<int> x3NbrTemp;
+
+//void findSolidNeighbor(SPtr<GbVoxelMatrix3D> voxelMatrix, int x1, int x2, int x3)
+//{
+//   for (int k3 = -1; k3<=1; k3++)
+//   {
+//      for (int k2 = -1; k2<=1; k2++)
+//      {
+//         for (int k1 = -1; k1<=1; k1++)
+//         {
+//            int j1 = x1+k1;
+//            int j2 = x2+k2;
+//            int j3 = x3+k3;
+//            if (j1>=0&&j1<nodesX1 && j2>=0&&j2<nodesX2 && j3>=0&&j3<nodesX3)
+//            {
+//               if ((*voxelMatrix)(j1, j2, j3)==GbVoxelMatrix3D::FLUID)
+//               {
+//                  if (flagMatrix(j1, j2, j3)==0)
+//                  {
+//                     voxelMatrixTemp(j1, j2, j3) = GbVoxelMatrix3D::SOLID;
+//                     flagMatrix(j1, j2, j3) = 1;
+//                     x1NbrTemp.push_back(j1);
+//                     x2NbrTemp.push_back(j2);
+//                     x3NbrTemp.push_back(j3);
+//                  }
+//               }
+//            }
+//         }
+//      }
+//   }
+//}
+
+void changePorosity(SPtr<GbVoxelMatrix3D> sample, vector<int> pmNX)
+{
+   int minX1 = 0;
+   int minX2 = 0;
+   int minX3 = 0;
+   int maxX1 = pmNX[0];
+   int maxX2 = pmNX[1];
+   int maxX3 = pmNX[2];
+   sample->calculateNumberOfSolidAndFluid();
+   double nSolid = sample->getNumberOfSolid();
+   double nFluid = sample->getNumberOfFluid();
+   double porosityStart = nFluid/(nSolid+nFluid);
+   double porosityEnd = 0.5;
+   double porosityStep = (porosityEnd-porosityStart)/(double)maxX3;
+   double totallSliceVoxel = maxX1*maxX2;
+   vector<int> fluidThreshold;
+
+   SPtr<GbVoxelMatrix3D> sampleTemp = SPtr<GbVoxelMatrix3D>(sample->clone());
+
+   int count=1;
+
+   for (int ix3=minX3; ix3<maxX3; ix3++)
+   {
+      int cFluid = 0;
+      for (int ix2=minX2; ix2<maxX2; ix2++)
+      {
+         for (int ix1=minX1; ix1<maxX1; ix1++)
+         {
+            if ((*sample)(ix1, ix2, ix3) == GbVoxelMatrix3D::FLUID)
+            {
+               cFluid++;
+            }
+         }
+      }
+      double slicePorosity = (double)cFluid/totallSliceVoxel;
+      double porosityPercent = (porosityStep*(double)count)/slicePorosity;
+      fluidThreshold.push_back((totallSliceVoxel-(double)cFluid)*porosityPercent);
+      count++;
+   }
+   int solidCount = 0;
+   count=0;
+
+   for (int ix3=minX3; ix3<maxX3; ix3++)
+   {
+      //while (fluidThreshold[count] > 0)
+      //{
+     // int fTh = fluidThreshold[count];
+
+         int solidCount = 0;
+         for (int ix2=minX2; ix2<maxX2; ix2++)
+         {
+            for (int ix1=minX1; ix1<maxX1; ix1++)
+            {
+               if ((*sample)(ix1, ix2, ix3) == GbVoxelMatrix3D::SOLID)
+               {
+                  bool flagSolid = true;
+                  for (int k3 = -1; k3<=1; k3++)
+                  {
+                     for (int k2 = -1; k2<=1; k2++)
+                     {
+                        for (int k1 = -1; k1<=1; k1++)
+                        {
+                           int j1 = ix1+k1;
+                           int j2 = ix2+k2;
+                           int j3 = ix3;//+k3;
+                           if (j1>=0&&j1<maxX1 && j2>=0&&j2<maxX2 && j3>=0&&j3<maxX3)
+                           {
+                              if ((*sample)(j1, j2, j3) == GbVoxelMatrix3D::FLUID)
+                              {
+                                 flagSolid = flagSolid && false;
+                              }
+                           }
+                        }
+                     }
+                  }
+                  if (!flagSolid)
+                  {
+                     (*sample)(ix1, ix2, ix3)=GbVoxelMatrix3D::FLUID;
+                      fluidThreshold[count]--;
+                     solidCount++;
+                  }
+                  if ( fluidThreshold[count] == 0)
+                  {
+                     ix1=maxX1;
+                     ix2=maxX2;
+                  }
+               }
+            }
+         }
+         UBLOG(logINFO, "count = " << count);
+         UBLOG(logINFO, "fTh = " <<  fluidThreshold[count]);
+         UBLOG(logINFO, "solidCount = " << solidCount);
+         //sample = sampleTemp;
+         //sampleTemp = SPtr<GbVoxelMatrix3D>(sample->clone());
+         
+        count++;     
+      
+       }
+      
+   //}
+   //sampleTemp->writeToLegacyVTKBinary("d:/temp/ChannelFlow/geo/sampleTemp.vtk");
+}
+
 //////////////////////////////////////////////////////////////////////////
 void run(string configname)
 {
@@ -56,7 +197,7 @@ void run(string configname)
       bool            thinWall          = config.getValue<bool>("thinWall");
       double          Re                = config.getValue<double>("Re");
       double          channelHigh       = config.getValue<double>("channelHigh");
-      bool            changeQs          = config.getValue<bool>("changeQs"); 
+      bool            changeQs          = config.getValue<bool>("changeQs");
       double          timeAvStart       = config.getValue<double>("timeAvStart");
       double          timeAvStop        = config.getValue<double>("timeAvStop");
       bool            averaging         = config.getValue<bool>("averaging");
@@ -127,12 +268,12 @@ void run(string configname)
       bcProc = SPtr<BCProcessor>(new BCProcessor());
 
       SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new IncompressibleCumulantLBMKernel());
-      
+
       mu::Parser fctForcingX1;
       fctForcingX1.SetExpr("Fx1");
       fctForcingX1.DefineConst("Fx1", 1.0e-6);
       kernel->setWithForcing(true);
-      
+
       kernel->setBCProcessor(bcProc);
 
       //////////////////////////////////////////////////////////////////////////
@@ -187,7 +328,7 @@ void run(string configname)
       {
          if (myid == 0) UBLOG(logINFO, "new start...");
 
-         
+
 
          grid->setPeriodicX1(true);
          grid->setPeriodicX2(true);
@@ -297,7 +438,7 @@ void run(string configname)
          intHelper.setBC();
 
          ////porous media
-         if(true)
+         if (true)
          {
             string samplePathname = pathGeo + "/" + sampleFilename;
 
@@ -314,6 +455,11 @@ void run(string configname)
             sample->setVoxelMatrixDelta(voxelDeltaX[0], voxelDeltaX[1], voxelDeltaX[2]);
             sample->setVoxelMatrixMininum(origin[0], origin[1], origin[2]);
 
+            changePorosity(sample, pmNX);
+
+            sample->writeToLegacyVTKBinary(pathOut+"/geo/sample.vtk");
+            return;
+
             int bounceBackOption = 1;
             bool vxFile = false;
             int i = 0;
@@ -421,7 +567,7 @@ void run(string configname)
 
          if (myid == 0) UBLOG(logINFO, "Restart - end");
       }
-     
+
       SPtr<UbScheduler> nupsSch(new UbScheduler(nupsStep[0], nupsStep[1], nupsStep[2]));
       std::shared_ptr<NUPSCounterCoProcessor> nupsCoProcessor(new NUPSCounterCoProcessor(grid, nupsSch, numOfThreads, comm));
 
@@ -450,15 +596,15 @@ void run(string configname)
       //DecrViscSch->addSchedule(10, 0, 1000);
       //DecreaseViscosityCoProcessor decrViscPPPtr(grid, DecrViscSch, &decrViscFunc, comm);
 
-	  //if (changeQs)
-	  //{
-		 // double z1 = pmL[2];
-		 // SPtr<IntegrateValuesHelper> intValHelp2(new IntegrateValuesHelper(grid, comm,
-			//  coord[0], coord[1], z1 - deltaXfine,
-			//  coord[3], coord[4], z1 + deltaXfine));
-		 // if (myid == 0) GbSystem3D::writeGeoObject(intValHelp2->getBoundingBox().get(), pathOut + "/geo/intValHelp2", WbWriterVtkXmlBinary::getInstance());
-		 // Utilities::ChangeRandomQs(intValHelp2);
-	  //}
+     //if (changeQs)
+     //{
+       // double z1 = pmL[2];
+       // SPtr<IntegrateValuesHelper> intValHelp2(new IntegrateValuesHelper(grid, comm,
+         //  coord[0], coord[1], z1 - deltaXfine,
+         //  coord[3], coord[4], z1 + deltaXfine));
+       // if (myid == 0) GbSystem3D::writeGeoObject(intValHelp2->getBoundingBox().get(), pathOut + "/geo/intValHelp2", WbWriterVtkXmlBinary::getInstance());
+       // Utilities::ChangeRandomQs(intValHelp2);
+     //}
 
       std::vector<double> levelCoords;
       std::vector<int> levels;
@@ -489,15 +635,15 @@ void run(string configname)
       //SPtr<CoProcessor> tav(new TimeAveragedValuesCoProcessor(grid, pathOut, WbWriterVtkXmlBinary::getInstance(), tavSch, comm,
       //   TimeAveragedValuesCoProcessor::Velocity | TimeAveragedValuesCoProcessor::Fluctuations | TimeAveragedValuesCoProcessor::Triplecorrelations,
       //   levels, levelCoords, bounds));
-      
-      
+
+
       //create line time series
-      SPtr<UbScheduler> tpcSch(new UbScheduler(1,1,3));
+      SPtr<UbScheduler> tpcSch(new UbScheduler(1, 1, 3));
       //GbPoint3DPtr p1(new GbPoint3D(0.0,0.005,0.01));
       //GbPoint3DPtr p2(new GbPoint3D(0.064,0.005,0.01));
       //SPtr<GbLine3D> line(new GbLine3D(p1.get(),p2.get()));
-      SPtr<GbLine3D> line(new GbLine3D(new GbPoint3D(0.0,0.005,0.01),new GbPoint3D(0.064,0.005,0.01)));
-      LineTimeSeriesCoProcessor lineTs(grid, tpcSch,pathOut+"/TimeSeries/line1.csv",line, 0,comm);
+      SPtr<GbLine3D> line(new GbLine3D(new GbPoint3D(0.0, 0.005, 0.01), new GbPoint3D(0.064, 0.005, 0.01)));
+      LineTimeSeriesCoProcessor lineTs(grid, tpcSch, pathOut+"/TimeSeries/line1.csv", line, 0, comm);
       if (myid==0) lineTs.writeLine(pathOut+"/geo/line1");
 
       if (myid == 0)
@@ -509,7 +655,7 @@ void run(string configname)
 
       omp_set_num_threads(numOfThreads);
       SPtr<UbScheduler> stepGhostLayer(new UbScheduler(1));
-      SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endTime));
+      SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, (int)endTime));
       calculator->addCoProcessor(nupsCoProcessor);
       calculator->addCoProcessor(AdjForcCoProcessor);
       calculator->addCoProcessor(migCoProcessor);