diff --git a/README b/README
deleted file mode 100644
index fd2ca756..00000000
--- a/README
+++ /dev/null
@@ -1,81 +0,0 @@
-/*---------------------------------------------------------------------------*\
-    CFDEMcoupling - Open Source CFD-DEM coupling
-
-    CFDEMcoupling is part of the CFDEMproject
-    www.cfdem.com
-                                Christoph Goniva, christoph.goniva@cfdem.com
-                                Copyright 2009-2012 JKU Linz
-                                Copyright 2012-     DCS Computing GmbH, Linz
--------------------------------------------------------------------------------
-License
-    This file is part of CFDEMcoupling.
-
-    CFDEMcoupling is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 3 of the License, or (at your
-    option) any later version.
-
-    CFDEMcoupling is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with CFDEMcoupling; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Description
-    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
-    and OpenFOAM. Note: this code is not part of OpenFOAM (see DISCLAIMER).
-\*---------------------------------------------------------------------------*/
-
-
-CFDEM coupling provides an open source parallel coupled CFD-DEM framework 
-combining the strengths of LIGGGHTS DEM code and the Open Source 
-CFD package OpenFOAM(R)(*). The CFDEMcoupling toolbox allows to expand 
-standard CFD solvers of OpenFOAM(R)(*) to include a coupling to the DEM 
-code LIGGGHTS. In this toolbox the particle representation within the 
-CFD solver is organized by "cloud" classes. Key functionalities are organised 
-in sub-models (e.g. force models, data exchange models, etc.) which can easily 
-be selected and combined by dictionary settings.
-
-The coupled solvers run fully parallel on distributed-memory clusters. 
-
-Features are:
-
-- its modular approach allows users to easily implement new models
-- its MPI parallelization enables to use it for large scale problems
-- the "forum"_lws on CFD-DEM gives the possibility to exchange with other 
-  users / developers
-- the use of GIT allows to easily update to the latest version
-- basic documentation is provided
-
-The file structure:
-
-- "src" directory including the source files of the coupling toolbox and models
-- "applications" directory including the solver files for coupled CFD-DEM simulations
-- "doc" directory including the documentation of CFDEMcoupling
-- "tutorials" directory including basic tutorial cases showing the functionality
-
-
-
-Details on installation are given on the "www.cfdem.com"
-
-The functionality of this CFD-DEM framwork is described via "tutorial cases" showing 
-how to use different solvers and models.
-
-CFDEMcoupling stands for Computational Fluid Dynamics (CFD) -
-Discrete Element Method (DEM) coupling.
-
-CFDEMcoupling is an open-source code, distributed freely under the terms of the 
-GNU Public License (GPL).
-
-Core development of CFDEMcoupling is done by 
-Christoph Goniva and Christoph Kloss, both at DCS Computing GmbH, 2012
-
-
-\*---------------------------------------------------------------------------*/
-(*) "OpenFOAM(R)"_of is a registered trade mark of the ESI Group.
-This offering is not affiliated, approved or endorsed by ESI Group, 
-the producer of the OpenFOAM® software and owner of the OpenFOAM® trade mark.
-\*---------------------------------------------------------------------------*/
diff --git a/doc/CFDEMcoupling_Manual.pdf b/doc/CFDEMcoupling_Manual.pdf
index 934f51e0..ba0dd6d1 100644
Binary files a/doc/CFDEMcoupling_Manual.pdf and b/doc/CFDEMcoupling_Manual.pdf differ
diff --git a/doc/githubAccess_public.pdf b/doc/githubAccess_public.pdf
index faf9d5ee..892d2ad2 100644
Binary files a/doc/githubAccess_public.pdf and b/doc/githubAccess_public.pdf differ
diff --git a/src/lagrangian/cfdemParticle/Make/files b/src/lagrangian/cfdemParticle/Make/files
index 80745e03..49170aa4 100644
--- a/src/lagrangian/cfdemParticle/Make/files
+++ b/src/lagrangian/cfdemParticle/Make/files
@@ -6,7 +6,6 @@ voidFractionModels = subModels/voidFractionModel
 locateModels = subModels/locateModel
 meshMotionModels = subModels/meshMotionModel
 momCoupleModels = subModels/momCoupleModel
-regionModels = subModels/regionModel
 dataExchangeModels = subModels/dataExchangeModel
 averagingModels = subModels/averagingModel
 clockModels = subModels/clockModel
@@ -22,20 +21,33 @@ $(forceModels)/forceModel/newForceModel.C
 $(forceModels)/noDrag/noDrag.C
 $(forceModels)/checkCouplingInterval/checkCouplingInterval.C
 $(forceModels)/DiFeliceDrag/DiFeliceDrag.C
+$(forceModels)/DiFeliceDragNLift/DiFeliceDragNLift.C
 $(forceModels)/GidaspowDrag/GidaspowDrag.C
 $(forceModels)/SchillerNaumannDrag/SchillerNaumannDrag.C
 $(forceModels)/Archimedes/Archimedes.C
 $(forceModels)/ArchimedesIB/ArchimedesIB.C
 $(forceModels)/interface/interface.C
 $(forceModels)/ShirgaonkarIB/ShirgaonkarIB.C
+$(forceModels)/fieldTimeAverage/fieldTimeAverage.C
+$(forceModels)/fieldBound/fieldBound.C
+$(forceModels)/volWeightedAverage/volWeightedAverage.C
+$(forceModels)/totalMomentumExchange/totalMomentumExchange.C
 $(forceModels)/KochHillDrag/KochHillDrag.C
+$(forceModels)/BeetstraDrag/multiphaseFlowBasic/multiphaseFlowBasic.C
+$(forceModels)/BeetstraDrag/BeetstraDrag.C
+$(forceModels)/LaEuScalarLiquid/LaEuScalarLiquid.C
 $(forceModels)/LaEuScalarTemp/LaEuScalarTemp.C
+$(forceModels)/LaEuScalarDust/LaEuScalarDust.C
 $(forceModels)/virtualMassForce/virtualMassForce.C
 $(forceModels)/gradPForce/gradPForce.C
+$(forceModels)/gradULiftForce/gradULiftForce.C
 $(forceModels)/viscForce/viscForce.C
 $(forceModels)/MeiLift/MeiLift.C
-$(forceModels)/fieldTimeAverage/fieldTimeAverage.C
-$(forceModels)/volWeightedAverage/volWeightedAverage.C
+$(forceModels)/KochHillDragNLift/KochHillDragNLift.C
+$(forceModels)/solidsPressureForce/solidsPressureForce.C
+$(forceModels)/periodicPressure/periodicPressure.C
+$(forceModels)/periodicPressureControl/periodicPressureControl.C
+$(forceModels)/averageSlipVel/averageSlipVel.C
 
 $(forceModelsMS)/forceModelMS/forceModelMS.C
 $(forceModelsMS)/forceModelMS/newForceModelMS.C
@@ -65,21 +77,18 @@ $(locateModels)/turboEngineSearch/turboEngineSearch.C
 $(locateModels)/turboEngineSearchM2M/turboEngineSearchM2M.C
 $(locateModels)/engineSearchIB/engineSearchIB.C
 
-
 $(meshMotionModels)/meshMotionModel/meshMotionModel.C
 $(meshMotionModels)/meshMotionModel/newMeshMotionModel.C
 $(meshMotionModels)/noMeshMotion/noMeshMotion.C
+$(meshMotionModels)/DEMdrivenMeshMotion/DEMdrivenMeshMotion.C
 
 $(momCoupleModels)/momCoupleModel/momCoupleModel.C
 $(momCoupleModels)/momCoupleModel/newMomCoupleModel.C
 $(momCoupleModels)/explicitCouple/explicitCouple.C
+$(momCoupleModels)/explicitCoupleSource/explicitCoupleSource.C
 $(momCoupleModels)/implicitCouple/implicitCouple.C
 $(momCoupleModels)/noCouple/noCouple.C
 
-$(regionModels)/regionModel/regionModel.C
-$(regionModels)/regionModel/newRegionModel.C
-$(regionModels)/allRegion/allRegion.C
-
 $(dataExchangeModels)/dataExchangeModel/dataExchangeModel.C
 $(dataExchangeModels)/dataExchangeModel/newDataExchangeModel.C
 $(dataExchangeModels)/oneWayVTK/oneWayVTK.C
@@ -110,4 +119,4 @@ $(smoothingModels)/smoothingModel/newSmoothingModel.C
 $(smoothingModels)/noSmoothing/noSmoothing.C
 $(smoothingModels)/constDiffSmoothing/constDiffSmoothing.C
 
-LIB = $(FOAM_USER_LIBBIN)/lib$(CFDEM_LIB_NAME)
+LIB = $(CFDEM_LIB_DIR)/lib$(CFDEM_LIB_NAME)
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/Allrun.sh b/tutorials/cfdemSolverPiso/ErgunTestCG/Allrun.sh
new file mode 100755
index 00000000..0de01eac
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/Allrun.sh
@@ -0,0 +1,22 @@
+#!/bin/bash
+
+#===================================================================#
+# allrun script for testcase as part of test routine 
+# run settlingTest
+# Christoph Goniva - Sept. 2010
+#===================================================================#
+
+#- define variables
+casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")"
+
+# check if mesh was built
+if [ -d "$casePath/CFD/constant/polyMesh/boundary" ]; then
+    echo "mesh was built before - using old mesh"
+else
+    echo "mesh needs to be built"
+    cd $casePath/CFD
+    blockMesh
+fi
+
+#- run parallel CFD-DEM in new terminal
+gnome-terminal --title='cfdemSolverPiso ErgunTestCG CFD'  -e "bash $casePath/parCFDDEMrun.sh" 
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/Ksl b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/Ksl
new file mode 100644
index 00000000..f44ecabe
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/Ksl
@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      Ksl;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -3 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    wall
+    {
+        type            zeroGradient;
+        //type            fixedValue;
+        //value           uniform 0;
+    }
+    inlet
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/U b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/U
new file mode 100644
index 00000000..bab71f9e
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/U
@@ -0,0 +1,89 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+
+    wall
+    {
+        //type            fixedValue;
+        //value           uniform (0 0 0);
+        type            slip;
+    }
+
+    inlet
+    {
+        /*type            flowRateInletVelocity;
+        flowRate        0.001;
+        value           uniform (0 0 0);*/
+
+        /*type            fixedValue;
+        value           uniform (0 0 0.0001);*/
+
+        //type            zeroGradient;
+
+/*        type            groovyBC;
+        variables       "Uend=vector(0,0,0.02);tEnd=0.1;";
+        valueExpression "((time() < tEnd) ? Uend/tEnd*time():Uend)";
+        value           uniform (0 0 0);*/
+
+        /*// 2.0.x, 1.6,ext
+        type            timeVaryingUniformFixedValue;
+        fileName        "steps_0p1s";
+        outOfBounds     clamp;
+        value           uniform (0 0 0);*/
+
+        // 2.1.x
+        type            uniformFixedValueVoidfraction;
+        uniformValue     table
+        (
+            (0.000 (0 0 0.002))
+            (0.020 (0 0 0.002))
+            (0.021 (0 0 0.004))
+            (0.040 (0 0 0.004))
+            (0.041 (0 0 0.006))
+            (0.060 (0 0 0.006))
+            (0.061 (0 0 0.008))
+            (0.080 (0 0 0.008))
+            (0.081 (0 0 0.010))
+            (0.100 (0 0 0.010))
+            (0.101 (0 0 0.012))
+            (0.120 (0 0 0.012))
+            (0.121 (0 0 0.014))
+            (0.140 (0 0 0.014))
+            (0.141 (0 0 0.016))
+            (0.160 (0 0 0.016))
+            (0.161 (0 0 0.018))
+            (0.180 (0 0 0.018))
+            (0.181 (0 0 0.020))
+            (0.200 (0 0 0.020))
+        );
+    }
+
+    outlet
+    {
+        /*type            fluxCorrectedVelocity; //inletOutlet;
+        value           uniform (0 0 0);
+        inletValue      uniform (0 0 0);*/
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/Us b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/Us
new file mode 100644
index 00000000..1af9a050
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/Us
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      Us;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    wall
+    {
+        type            zeroGradient;
+        //type            fixedValue;
+        //value           uniform (1 0 0);
+    }
+    inlet
+    {
+        type            zeroGradient;
+        //type            fixedValue;
+        //value           uniform (1 0 0);;
+    }
+    outlet
+    {
+        type            zeroGradient;
+        //type            fixedValue;
+        //value           uniform (1 0 0);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/epsilon b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/epsilon
new file mode 100644
index 00000000..37ffab5e
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/epsilon
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 0.003;
+
+boundaryField
+{
+    wall
+    {
+        type            epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0.003;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 0.003;
+        value           uniform 0.003;
+    }
+    inlet
+    {
+        /*type            turbulentMixingLengthDissipationRateInlet;
+        mixingLength    0.005;
+        value           uniform 0.003;*/
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/k b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/k
new file mode 100644
index 00000000..bbf04afb
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/k
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.001;
+
+boundaryField
+{
+    wall
+    {
+        type            kqRWallFunction;
+        value           uniform 0.001;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 0.001;
+        value           uniform 0.001;
+    }
+    inlet
+    {
+        /*type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.05;
+        value           uniform 0.001;*/
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/nut b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/nut
new file mode 100644
index 00000000..dcd36c89
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/nut
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    wall
+    {
+        type            nutWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/p b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/p
new file mode 100644
index 00000000..d111e6be
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/p
@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 1.0e5;
+
+boundaryField
+{
+
+    wall
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+        type            zeroGradient;
+        //type            fixedValue;
+        //value           uniform 100000;
+    }
+    outlet
+    {
+        //type            zeroGradient;
+        type            fixedValue;
+        value           uniform 100000;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/rho b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/rho
new file mode 100644
index 00000000..2c9e8cd0
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/rho
@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      rho;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -3 0 0 0 0 0];
+
+internalField   uniform 30;
+
+boundaryField
+{
+
+    wall
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/voidfraction b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/voidfraction
new file mode 100644
index 00000000..8298f8c4
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/voidfraction
@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      voidfraction;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    wall
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+    inlet
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/RASProperties b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/RASProperties
new file mode 100644
index 00000000..daec1829
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/RASProperties
@@ -0,0 +1,25 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      RASProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RASModel        laminar;
+
+turbulence      off;
+
+printCoeffs     on;
+
+
+// ************************************************************************* //
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties
new file mode 100644
index 00000000..c93374b1
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties
@@ -0,0 +1,177 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.4                                   |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "";
+    local           "";
+
+    class           dictionary;
+    object          couplingProperties;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//===========================================================================//
+// sub-models & settings
+
+modelType "A"; // A or B
+
+couplingInterval 100;
+
+voidFractionModel divided;
+
+locateModel engine;
+
+meshMotionModel noMeshMotion;
+
+regionModel allRegion;
+
+IOModel "basicIO";
+
+dataExchangeModel twoWayMPI;//twoWayFiles;//oneWayVTK;//
+
+averagingModel dense;//dilute;//
+
+clockModel off;//standardClock;//
+
+forceModels
+(
+    //GidaspowDrag
+    //DiFeliceDrag
+    KochHillDrag
+    gradPForce
+    viscForce
+    //Archimedes
+    //volWeightedAverage
+    //totalMomentumExchange
+);
+
+momCoupleModels
+(
+    implicitCouple
+);
+
+turbulenceModelType "RASProperties";//"LESProperties";//
+
+//===========================================================================//
+// sub-model properties
+
+implicitCoupleProps
+{
+    velFieldName "U";
+    granVelFieldName "Us";
+    voidfractionFieldName "voidfraction";
+}
+
+ArchimedesProps
+{
+    densityFieldName "rho";
+    gravityFieldName "g";
+}
+gradPForceProps
+{
+    pFieldName "p";
+    densityFieldName "rho";
+    voidfractionFieldName "voidfraction";
+    velocityFieldName "U";
+    //interpolation;
+}
+
+viscForceProps
+{
+    velocityFieldName "U";
+    densityFieldName "rho";
+    interpolation;
+}
+volWeightedAverageProps
+{
+    scalarFieldNames
+    (
+        voidfraction
+    );
+    vectorFieldNames
+    (
+    );
+    upperThreshold 0.999;
+    lowerThreshold 0;
+    verbose;
+}
+totalMomentumExchangeProps
+{
+    implicitMomExFieldName "Ksl";
+    explicitMomExFieldName "none";
+    fluidVelFieldName "U";
+    granVelFieldName "Us";
+    densityFieldName "rho";
+}
+GidaspowDragProps
+{
+    velFieldName "U";
+    densityFieldName "rho";
+}
+DiFeliceDragProps
+{
+    velFieldName "U";
+    densityFieldName "rho";
+    voidfractionFieldName "voidfraction";
+}
+
+KochHillDragProps
+{
+    velFieldName "U";
+    densityFieldName "rho";
+    voidfractionFieldName "voidfraction";
+    //verbose;
+}
+
+virtualMassForceProps
+{
+    velFieldName "U";
+    densityFieldName "rho";
+}
+
+oneWayVTKProps
+{
+    couplingFilename "vtk_out%4.4d.vtk";
+    maxNumberOfParticles 30000;
+}
+
+twoWayFilesProps
+{
+    maxNumberOfParticles 10100;
+}
+
+centreProps
+{
+    alphaMin 0.10;
+}
+
+engineProps
+{
+    treeSearch true;
+}
+
+dividedProps
+{
+    alphaMin 0.1;
+    scaleUpVol 1.0;
+}
+
+twoWayMPIProps
+{
+    maxNumberOfParticles 10100;
+    liggghtsPath "../DEM/in.liggghts_resume";
+}
+// ************************************************************************* //
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/g b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/g
new file mode 100644
index 00000000..317bdd50
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/g
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           ( 0 0 -9.81 );
+
+
+// ************************************************************************* //
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/liggghtsCommands b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/liggghtsCommands
new file mode 100644
index 00000000..ee883115
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/liggghtsCommands
@@ -0,0 +1,30 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.4                                   |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "";
+    local           "";
+
+    class           dictionary;
+    object          liggghtsCommands;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+liggghtsCommandModels
+(
+   runLiggghts
+);
+// ************************************************************************* //
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/polyMesh/blockMeshDict b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/polyMesh/blockMeshDict
new file mode 100644
index 00000000..5e0fe81c
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/polyMesh/blockMeshDict
@@ -0,0 +1,154 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.0                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+version  2.0;
+format   ascii;
+class       dictionary;
+object      blockMeshDict;
+}
+// ************************************
+
+
+
+
+
+
+   meshGenApp blockMesh;
+   convertToMeters 1;
+
+    //64 mm column diameter
+    //15 cm length
+   
+   
+   
+    //Width of middle square section
+   
+   
+   
+   
+    //how many cells in the square section
+    //how many cells from square section to perimeter
+    // how many cells from top to bottom
+
+   vertices
+   (
+    ( 0.0069 0.0069 0.0  ) // Vertex fiveoclocksqb = 0 
+    (-0.0069 0.0069 0.0  ) // Vertex sevenoclocksqb = 1 
+    (-0.0069 -0.0069 0.0 ) // Vertex elevenoclocksqb = 2 
+    ( 0.0069 -0.0069 0.0 ) // Vertex oneoclocksqb = 3 
+   
+    ( 0.00975807358913172 0.00975807357161699 0.0 ) // Vertex fiveoclockcb = 4 
+    (-0.00975807358913172 0.00975807357161699 0.0) // Vertex sevenoclockcb = 5 
+    (-0.00975807358913172 -0.00975807357161699 0.0) // Vertex elevenoclockcb = 6 
+    ( 0.00975807358913172 -0.00975807357161699 0.0) // Vertex oneoclockcb = 7 
+
+    ( 0.0069 0.0069 0.0553) // Vertex fiveoclocksqt = 8 
+    (-0.0069 0.0069 0.0553) // Vertex sevenoclocksqt = 9 
+    (-0.0069 -0.0069 0.0553) // Vertex elevenoclocksqt = 10 
+    ( 0.0069 -0.0069 0.0553) // Vertex oneoclocksqt = 11 
+   
+    ( 0.00975807358913172 0.00975807357161699 0.0553) // Vertex fiveoclockct = 12 
+    (-0.00975807358913172 0.00975807357161699 0.0553) // Vertex sevenoclockct = 13 
+    (-0.00975807358913172 -0.00975807357161699 0.0553) // Vertex elevenoclockct = 14 
+    ( 0.00975807358913172 -0.00975807357161699 0.0553) // Vertex oneoclockct = 15 
+   );				
+
+   blocks
+   (
+    //square block
+    hex (
+       2 3 0 1   
+       10 11 8 9   
+       )
+    (8 8 24)
+    simpleGrading (1 1 1)
+
+    //slice1
+    hex (
+       1 0 4 5
+       9 8 12 13
+       )
+    (8 4 24)
+    simpleGrading (1 1 1)
+
+    //slice2
+    hex (
+       6 2 1  5 
+       14 10 9   13 
+       )
+   ( 4 8 24)
+simpleGrading (1 1 1)
+
+   //slice3
+   hex (
+         6 7 3 2   
+         14 15 11 10   
+       )
+   (8 4 24)
+simpleGrading (1 1 1)
+
+   //slice4
+   hex (
+         3 7 4 0  
+         11 15 12 8  
+       )
+   (4 8 24)
+simpleGrading (1 1 1)
+
+   );
+
+
+   //create the quarter circles
+   edges
+   (
+    arc 4 5 (0.0 0.0138 0.0 )
+    arc 5 6 (-0.0138 0.0 0.0)
+    arc 6 7 (0.0 -0.0138 0.0 )
+    arc 7 4 (0.0138 0.0 0.0)
+
+    arc 12 13 (0.0 0.0138 0.0553 )
+    arc 13 14 (-0.0138 0.0 0.0553 )
+    arc 14 15 (0.0 -0.0138 0.0553 )
+    arc 15 12 (0.0138 0.0 0.0553 )
+
+   );
+
+   patches
+   (
+    patch inlet
+    (
+     (0 3 2 1)
+     (0 4 7 3)
+     (4 0 1 5)
+     (1 2 6 5)
+     (3 7 6 2)
+    )
+
+    patch outlet
+    (
+     (8 11 10 9)
+     (8 12 15 11)
+     (12 8 9 13)
+     (9 10 14 13)
+     (11 15 14 10)
+    )
+
+    wall wall
+    (
+     (5 4 12 13)
+     (5 13 14 6)
+     (6 14 15 7)
+     (7 15 12 4)
+    )
+
+);
+
+mergePatchPairs
+(
+);
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/transportProperties b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/transportProperties
new file mode 100644
index 00000000..374975c9
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/transportProperties
@@ -0,0 +1,38 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+transportModel  Newtonian;
+
+nu              nu [ 0 2 -1 0 0 0 0 ] 1.5e-04;
+
+CrossPowerLawCoeffs
+{
+    nu0             nu0 [ 0 2 -1 0 0 0 0 ] 1e-06;
+    nuInf           nuInf [ 0 2 -1 0 0 0 0 ] 1e-06;
+    m               m [ 0 0 1 0 0 0 0 ] 1;
+    n               n [ 0 0 0 0 0 0 0 ] 1;
+}
+
+BirdCarreauCoeffs
+{
+    nu0             nu0 [ 0 2 -1 0 0 0 0 ] 1e-06;
+    nuInf           nuInf [ 0 2 -1 0 0 0 0 ] 1e-06;
+    k               k [ 0 0 1 0 0 0 0 ] 0;
+    n               n [ 0 0 0 0 0 0 0 ] 1;
+}
+// ************************************************************************* //
+
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/turbulenceProperties b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/turbulenceProperties
new file mode 100644
index 00000000..f6753662
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/turbulenceProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  RASModel;
+
+
+// ************************************************************************* //
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/octave/loaddata.m b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/octave/loaddata.m
new file mode 100644
index 00000000..9897140c
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/octave/loaddata.m
@@ -0,0 +1,69 @@
+function [data,colname]=loaddata(filename,columns,headerlines)
+
+% data = loaddata(filename,columns,headerlines)
+% data = loaddata(filename,0,headerlines): automatc detection of 
+%           columns by word count in first headerline or first data line
+%           (if headerline==0), separated by blanks or tabs
+%           columns=-1: write colum assignment code lines
+%              works only if headerline>0 and colnames separated by 
+%              not more than 1 blank!
+%           columns=-2: columns separated by 2 blanks (fluent/scheme auswertung-instat.scm)
+% [data, colname] = loaddata(filename,0,headerlines) 
+%           colname = cell string array of column names from headerline
+
+fprintf(1,'loading %s ...  ', filename);
+f=fopen(filename,'r');
+if f==-1
+  fprintf(1,'\n*** error: could not open "%s" ...\n', filename);
+  data=[];
+else
+  for i=1:headerlines
+    if i==1 s=fgets(f); else fgets(f); end
+  end
+  if headerlines==0
+    s=fgets(f); frewind(f);
+  end
+  if columns<=0 % & headerlines>=1
+    fprintf(1,'\n');
+    fprintf(1,'  %s',s);
+    pos = findstr(sprintf('\t'),s); % trennzeichenpositionen in zeile suchen
+    if length(pos)>0
+      endpos=length(s);
+      %if isspace(s(endpos-1)) endpos=endpos-1; end
+      pos = [0 pos(1,:) endpos]; % anfangs- und endposition hinzufügen
+    elseif columns==-2 % spalteneinträge durch 2 blanks getrennt (fluent/scheme)
+      pos=[0];
+      for i=1:length(s)-1
+        if isspace(s(i)) & isspace(s(i+1))
+          pos(end+1)=i+1;
+        end
+      end
+      pos(end+1)=length(s);      
+    else % spalteneinträge durch blanks getrennt     
+      %pos = findstr(' ',s); % trennzeichenpositionen in zeile suchen
+      pos=[];
+      word=0;
+      for i=1:length(s)
+        if isspace(s(i))
+          word=0;
+        elseif ~word
+          word=1;
+          pos(end+1)=i-1;
+        end
+      end
+      pos(end+1)=length(s);
+    end
+    if headerlines>=1
+      for i=1:length(pos)-1 % alle spalten
+        colname{i}=s(pos(i)+1:pos(i+1)-1);
+        if columns==-1 fprintf(' = data(strmatch(''%s'',colname),:); %% column %d\n', colname{i},i); end
+        if columns==-2 fprintf('%% column %d: %s\n', i, colname{i}); end
+      end   
+    end
+    columns = length(pos)-1;
+    fprintf(1,'  total: %d columns ',columns);
+  end
+  data=fscanf(f,'%f',[columns,inf]);
+  fclose(f);
+  fprintf(1,'done.\n');
+end
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/octave/scaling.m b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/octave/scaling.m
new file mode 100644
index 00000000..df1bda10
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/octave/scaling.m
@@ -0,0 +1,22 @@
+close all;
+clear;
+clc;
+
+%====================================%
+% setting up the scaling factor
+%====================================%
+
+%-constants
+g=9.81
+
+%- particle props
+dp=0.0008
+rhoP=200
+Vp=dp^3*pi/6
+np=19531
+
+%- scaling
+dParcel=dp*2
+VParcel=dParcel^3*pi/6
+nParcel=np*Vp/VParcel
+dragForceScale=dParcel/dp
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/octave/totalPressureDrop.m b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/octave/totalPressureDrop.m
new file mode 100644
index 00000000..2014d778
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/octave/totalPressureDrop.m
@@ -0,0 +1,96 @@
+close all;
+clear;
+clc;
+
+%====================================%
+% simulation data 1
+%====================================%
+rhoG = 30			% density in kg/m3
+%path = '../probes/0/p';
+path = '../probes/0/p';
+columns=22;
+headerlines=4;
+data = loaddata(path,columns,headerlines);
+data=transpose(data);
+[x,y]=size(data)
+dp_sim = (data(:,2)-data(:,y))*rhoG;
+t_sim = data(:,1);
+%fprintf('final pressureDrop of sim = %f Pa\n',dp_sim(length(dp_sim)) )
+
+%====================================%
+% analytical calculation
+%====================================%
+
+%===================
+% Ergun Equation
+%===================
+fprintf('\ncalc Ergun eqn:\n')
+dp = 0.0008			% particle diameter
+phip = 1			% sphericity
+epsilon = 0.451335              % void fraction
+Ustart = 0.002
+Uend = 0.02
+timeStepSize = 0.001;            % time interval of pressure data
+Tstart = 0;
+Tend = t_sim(length(t_sim));
+deltaU=(Uend-Ustart)/((Tend-Tstart)/timeStepSize);
+U = Ustart+deltaU:deltaU:Uend;  % velocity over time
+Ua = U / epsilon;		% physical velocity
+L = 0.0156			% length of bed
+nuG = 1.5*10^-4			% kinemat Visk in m2/s
+muG = nuG*rhoG			% dynam visc in Pa s
+
+dpErgun= L * (
+                150*((1-epsilon)^2/epsilon^3)*((muG.*U)/(phip*dp)^2) 
+              +1.75*((1-epsilon)/epsilon^3)*((rhoG.*U.^2)/(phip*dp))
+        );
+
+fprintf('NOTE: this pressure is divided by density (according to CFD solver)\n')
+fprintf('so the result does not depend on density\n')
+
+%fprintf('final pressure drop (Ergun eqn)= %f Pa\n',dpErgun)
+
+%==================================
+% min fluidization velocity in m/s
+%==================================
+rhoP = 7000                      % particle density in kg/m3
+g = 9.81                        % gravity m/s2
+Umf = dp^2*(rhoP-rhoG)*g/(150*muG)*(epsilon^3*phip^2)/(1-epsilon)
+ReMF = Umf*dp*rhoG/muG 			% must be <20 !!!
+%Umf = sqrt(phip*dp^2/1.75*(rhoP-rhoG)/rhoG*g*epsilon^3) % Re>1000
+
+dpUmf= L * (
+                150*((1-epsilon)^2/epsilon^3)*((muG.*Umf)/(phip*dp)^2) 
+              +1.75*((1-epsilon)/epsilon^3)*((rhoG.*Umf.^2)/(phip*dp))
+        );
+
+%====================================%
+% plot data
+%====================================%
+length(U)
+length(dp_sim)
+figure(2)
+plot(U,dp_sim)
+title("Ergun pressure drop vs. simulation")
+a=strcat("analytical (Ergun), Umf=",num2str(Umf),", dpUmf=",num2str(dpUmf));
+legend(a,"simulation")
+xlabel("velocity in [m/s]")
+ylabel("pressure drop [Pa]")
+axis([0,Uend,0,dpErgun(length(dpErgun))])
+
+figure(1)
+plot(U,dpErgun,U,dp_sim,[Umf,Uend],dpUmf*ones(1,2))
+title("Ergun pressure drop vs. simulation")
+a=strcat("analytical (Ergun), Umf=",num2str(Umf),", dpUmf=",num2str(dpUmf));
+legend(a,"simulation")
+xlabel("velocity in [m/s]")
+ylabel("pressure drop [Pa]")
+axis([0,Uend,0,dpErgun(length(dpErgun))])
+
+%print('cfdemSolverPiso_settlingTest.eps','-deps2')
+print -color "cfdemSolverPiso_ErgunTestMPI.eps"
+replot;
+
+
+
+
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/steps_0p1s b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/steps_0p1s
new file mode 100644
index 00000000..0b9befc7
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/steps_0p1s
@@ -0,0 +1,22 @@
+(
+    (0.000 (0 0 0.002))
+    (0.010 (0 0 0.002))
+    (0.011 (0 0 0.004))
+    (0.020 (0 0 0.004))
+    (0.021 (0 0 0.006))
+    (0.030 (0 0 0.006))
+    (0.031 (0 0 0.008))
+    (0.040 (0 0 0.008))
+    (0.041 (0 0 0.010))
+    (0.050 (0 0 0.010))
+    (0.051 (0 0 0.012))
+    (0.060 (0 0 0.012))
+    (0.061 (0 0 0.014))
+    (0.070 (0 0 0.014))
+    (0.071 (0 0 0.016))
+    (0.080 (0 0 0.016))
+    (0.081 (0 0 0.018))
+    (0.090 (0 0 0.018))
+    (0.091 (0 0 0.020))
+    (0.100 (0 0 0.020))
+)
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/controlDict b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/controlDict
new file mode 100644
index 00000000..8a1c1733
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/controlDict
@@ -0,0 +1,168 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     pisoFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         0.2;//0.01;
+
+deltaT          0.001;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.005;//0.01;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  no;
+
+maxCo 0.1;
+
+//libs ( "libgroovyBC.so" );
+libs ( "libfiniteVolumeCFDEM.so" );
+
+functions
+(
+
+    probes
+    {
+        type        probes;
+        // Where to load it from
+        functionObjectLibs ( "libsampling.so" );  
+        // Name of the directory for probe data
+        name        probes;
+        probeLocations
+        (
+            (0 0 0.0001)
+            (0 0 0.0026)
+            (0 0 0.0051)
+            (0 0 0.0076)
+            (0 0 0.0101)
+            (0 0 0.0126)
+            (0 0 0.0151)
+            (0 0 0.0176)
+            (0 0 0.0201)
+            (0 0 0.0226)
+            (0 0 0.0251)
+            (0 0 0.0276)
+            (0 0 0.0301)
+            (0 0 0.0326)
+            (0 0 0.0351)
+            (0 0 0.0375)
+            (0 0 0.0401)
+            (0 0 0.0426)
+            (0 0 0.0451)
+            (0 0 0.0476)
+            (0 0 0.0529)
+        );
+
+        // Fields to be probed
+        fields ( p U voidfraction volAverage_voidfraction);
+
+        // Write at same frequency as fields
+        outputControl   timeStep;//outputTime;
+        outputInterval  1;
+    }
+
+    probes_inletFlux
+    {
+        type            faceSource;
+        functionObjectLibs ("libfieldFunctionObjects.so");
+        enabled         true;
+        outputControl   timeStep;   //outputTime;
+        log             true;       // log to screen?
+        valueOutput     true;       // Write values at run-time output times?
+        source          patch;   // Type of face source:
+        sourceName      inlet;         // faceZone name, see below
+        operation       sum;
+        fields
+        (
+            phi
+        );
+    }
+
+    probes_outletFlux
+    {
+        type            faceSource;
+        functionObjectLibs ("libfieldFunctionObjects.so");
+        enabled         true;
+        outputControl   timeStep;   //outputTime;
+        log             true;       // log to screen?
+        valueOutput     true;       // Write values at run-time output times?
+        source          patch;   // Type of face source:
+        sourceName      outlet;         // faceZone name, see below
+        operation       sum;
+        fields
+        (
+            phi
+        );
+    }
+
+    probes_wallFlux
+    {
+        type            faceSource;
+        functionObjectLibs ("libfieldFunctionObjects.so");
+        enabled         true;
+        outputControl   timeStep;   //outputTime;
+        log             true;       // log to screen?
+        valueOutput     true;       // Write values at run-time output times?
+        source          patch;   // Type of face source:
+        sourceName      wall;         // faceZone name, see below
+        operation       sum;
+        fields
+        (
+            phi
+        );
+    }
+   /*pressureDrop
+   {
+     type patchAverage;
+     functionObjectLibs
+     (
+         "libsimpleFunctionObjects.so"
+     );
+     verbose true;
+     patches
+     (
+         inlet
+         outlet
+     );
+     fields
+     (
+         p
+     );
+     factor  1;
+   }*/
+);
+// ************************************************************************* //
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/controlDict.foam b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/controlDict.foam
new file mode 100644
index 00000000..7493598b
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/controlDict.foam
@@ -0,0 +1,117 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     pisoFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         0.1;//0.01;
+
+deltaT          0.001;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.001;//0.01;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  no;
+
+maxCo 0.1;
+
+libs ("libOpenFOAM.so" "libgroovyBC.so");
+
+functions
+(
+
+    probes
+    {
+        type        probes;
+        // Where to load it from
+        functionObjectLibs ( "libsampling.so" );  
+        // Name of the directory for probe data
+        name        probes;
+        probeLocations
+        (
+            (0 0 0.0001)
+            (0 0 0.0026)
+            (0 0 0.0051)
+            (0 0 0.0076)
+            (0 0 0.0101)
+            (0 0 0.0126)
+            (0 0 0.0151)
+            (0 0 0.0176)
+            (0 0 0.0201)
+            (0 0 0.0226)
+            (0 0 0.0251)
+            (0 0 0.0276)
+            (0 0 0.0301)
+            (0 0 0.0326)
+            (0 0 0.0351)
+            (0 0 0.0375)
+            (0 0 0.0401)
+            (0 0 0.0426)
+            (0 0 0.0451)
+            (0 0 0.0476)
+            (0 0 0.0529)
+        );
+
+        // Fields to be probed
+        fields ( p U voidfraction volAverage_voidfraction);
+
+        // Write at same frequency as fields
+        outputControl   timeStep;//outputTime;
+        outputInterval  1;
+    }
+
+   /*pressureDrop
+   {
+     type patchAverage;
+     functionObjectLibs
+     (
+         "libsimpleFunctionObjects.so"
+     );
+     verbose true;
+     patches
+     (
+         inlet
+         outlet
+     );
+     fields
+     (
+         p
+     );
+     factor  1;
+   }*/
+);
+// ************************************************************************* //
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/decomposeParDict b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/decomposeParDict
new file mode 100644
index 00000000..bccec657
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/decomposeParDict
@@ -0,0 +1,86 @@
+/*-------------------------------*- C++ -*---------------------------------*\
+|    =========                                                              |
+|    \\      /     OpenFOAM                                                 |
+|     \\    /                                                               |
+|      \\  /       The Open Source CFD Toolbox                              |
+|       \\/                                        http://www.OpenFOAM.org  |
+\*-------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    note        "mesh decomposition control dictionary";
+    location    "system";
+    object      decomposeParDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains  4;
+
+//- Keep owner and neighbour on same processor for faces in zones:
+// preserveFaceZones (heater solid1 solid3);
+
+//method          scotch;
+// method          hierarchical;
+ method          simple;
+// method          metis;
+// method          manual;
+
+simpleCoeffs
+{
+    n           (1 2 2);
+    delta       0.001;
+}
+
+hierarchicalCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+    order       xyz;
+}
+
+metisCoeffs
+{
+ /*
+    processorWeights
+    (
+        1
+        1
+        1
+        1
+    );
+  */
+}
+
+scotchCoeffs
+{
+    //processorWeights
+    //(
+    //    1
+    //    1
+    //    1
+    //    1
+    //);
+    //writeGraph  true;
+    //strategy "b";
+}
+
+manualCoeffs
+{
+    dataFile    "decompositionData";
+}
+
+
+//// Is the case distributed
+//distributed     yes;
+//// Per slave (so nProcs-1 entries) the directory above the case.
+//roots           
+//(
+//    "/tmp"
+//    "/tmp"
+//);
+
+
+// ************************************************************************* //
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/funkySetFieldsDict b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/funkySetFieldsDict
new file mode 100755
index 00000000..39414dd5
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/funkySetFieldsDict
@@ -0,0 +1,33 @@
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "/.automount/werner/Werner/bgschaid/bgschaid-foamStuff/DepositionTests";
+    case            "neutralTest";
+    instance        "system";
+    local           "";
+
+    class           dictionary;
+    object          funkySetFieldsDict;
+}
+
+expressions
+(
+	init
+	{
+		field voidfraction;
+		expression "1";
+                valuePatches 2(inlet outlet); // set value 1 at those patches
+	}
+	setGamma
+	{
+		field voidfraction;
+		expression "0.43";
+		//condition "pow(pow(pos().x-0.2,2) + pow(pos().y-0,2)+ pow(pos().z-0,2),0.5) < 0.04";
+                condition "pos().z < 0.01467 && pos().z > 0.0";
+                keepPatches true;
+                //valuePatches 1(inlet); // set value 1 at those patches
+	}
+	
+);
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/fvSchemes b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/fvSchemes
new file mode 100644
index 00000000..fac85f68
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/fvSchemes
@@ -0,0 +1,76 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(p)         Gauss linear;
+    grad(U)         Gauss linear;
+}
+
+divSchemes
+{
+    default         Gauss linear;
+    div(phi,U)      Gauss limitedLinearV 1;
+    div(phi,k)      Gauss limitedLinear 1;
+    div(phi,epsilon) Gauss limitedLinear 1;
+    div(phi,R)      Gauss limitedLinear 1;
+    div(R)          Gauss linear;
+    div(phi,nuTilda) Gauss limitedLinear 1;
+    div((viscousTerm*dev(grad(U).T()))) Gauss linear;
+    div((nu*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(grad(U).T()))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+    laplacian(viscousTerm,U) Gauss linear corrected;
+    laplacian(nu,U) Gauss linear corrected;
+    laplacian(nuEff,U) Gauss linear corrected;
+    laplacian((1|A(U)),p) Gauss linear corrected;
+    laplacian((voidfraction2|A(U)),p) Gauss linear corrected;
+    laplacian(DkEff,k) Gauss linear corrected;
+    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
+    laplacian(DREff,R) Gauss linear corrected;
+    laplacian(DnuTildaEff,nuTilda) Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+    interpolate(U)  linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p               ;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/fvSolution b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/fvSolution
new file mode 100644
index 00000000..955f5b14
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/fvSolution
@@ -0,0 +1,86 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-06;
+        relTol          0.1;
+    }
+
+    pFinal
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-06;
+        relTol          0;
+    }
+
+    U
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0;
+    }
+
+    k
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0;
+    }
+
+    epsilon
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0;
+    }
+
+    R
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0;
+    }
+
+    nuTilda
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0;
+    }
+}
+
+PISO
+{
+    nCorrectors     4;
+    nNonOrthogonalCorrectors 1;
+    pRefCell        0;
+    pRefValue       0;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/in.liggghts_init b/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/in.liggghts_init
new file mode 100755
index 00000000..4f82b3f3
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/in.liggghts_init
@@ -0,0 +1,78 @@
+#######################################################
+# tutorial case showing coarse graining functionality #
+# Christoph Goniva  08.02.2013                        #
+#######################################################
+echo both
+
+#######################################################
+# variables                                           #
+#######################################################
+variable cg equal 2
+variable nPorg equal 19531
+variable nPscal equal ${nPorg}/(${cg}*${cg}*${cg})
+#######################################################
+
+atom_style	granular
+atom_modify	map array
+communicate	single vel yes
+coarsegraining ${cg}
+
+boundary	m m m
+newton		off
+
+units		si
+
+region		reg block -0.015 0.015 -0.015 0.015 -0.001 0.0554 units box
+create_box	1 reg
+
+neighbor	0.001 bin
+neigh_modify	delay 0
+
+
+#Material properties required for new pair styles
+
+fix 		m1 all property/global youngsModulus peratomtype 5.e6
+fix 		m2 all property/global poissonsRatio peratomtype 0.45
+fix 		m3 all property/global coefficientRestitution peratomtypepair 1 0.3
+fix 		m4 all property/global coefficientFriction peratomtypepair 1 0.5
+fix 		m5 all property/global characteristicVelocity scalar 1.
+
+#pair style
+pair_style 	gran/hooke/history
+pair_coeff	* *
+
+#timestep, gravity
+timestep	0.000005
+fix		gravi all gravity 9.81 vector 0.0 0.0 -1.0
+
+#walls
+fix			zwalls1 all wall/gran/hooke/history     primitive type 1  zplane 0.0
+fix			zwalls2 all wall/gran/hooke/history     primitive type 1  zplane 0.0553
+fix         cylwalls all wall/gran/hooke/history    primitive type 1 zcylinder 0.01385 0. 0. tangential_damping off
+
+#particle distributions and insertion
+region		bc cylinder z 0.0 0.0 0.012 0. 0.055 units box
+fix		    pts1 all particletemplate/sphere 1 atom_type 1 density constant 200 radius constant 0.0004
+fix		    pdd1 all particledistribution/discrete 1.  1 pts1 1.0
+
+fix		ins all insert/pack seed 100001 distributiontemplate pdd1 vel constant 0. 0. -1. insert_every once overlapcheck yes all_in yes particles_in_region ${nPscal} region bc
+
+#apply nve integration to all particles that are inserted as single particles
+fix		integr all nve/sphere
+
+#screen output
+compute		1 all erotate/sphere
+thermo_style	custom step atoms ke c_1 vol
+thermo		1000
+thermo_modify	lost ignore norm no
+compute_modify	thermo_temp dynamic yes
+
+#insert the first particles so that dump is not empty
+run		1
+dump		dmp all custom 5000 post/dump.liggghts_init id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius 
+#force : f_couple_cfd[0] f_couple_cfd[1] f_couple_cfd[2]
+#node : f_couple_cfd[6]
+#cell id : f_couple_cfd[7]
+
+run		10000 upto
+write_restart 	liggghts.restart_coarseGraining_${cg}
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/in.liggghts_resume b/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/in.liggghts_resume
new file mode 100644
index 00000000..b64765fc
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/in.liggghts_resume
@@ -0,0 +1,86 @@
+#######################################################
+# tutorial case showing coarse graining functionality #
+# Christoph Goniva  08.02.2013                        #
+#######################################################
+echo both
+
+#######################################################
+# variables                                           #
+#######################################################
+variable cg equal 2
+variable nPorg equal 19531
+variable nPscal equal ${nPorg}/(${cg}*${cg}*${cg})
+#######################################################
+
+atom_style	granular
+atom_modify	map array
+communicate	single vel yes
+coarsegraining ${cg}
+
+boundary	m m m
+newton		off
+
+units		si
+processors	1 2 2
+
+#read the restart file
+read_restart 	../DEM/liggghts.restart_coarseGraining_${cg}
+
+#do not do this here, the simulation box is in the restart file!
+#region		reg block -0.015 0.015 -0.015 0.015 -0.001 0.0554 units box
+#create_box	1 reg
+
+neighbor	0.0005 bin
+neigh_modify	delay 0
+
+#Material properties required for new pair styles
+
+fix 		m1 all property/global youngsModulus peratomtype 5.e6
+fix 		m2 all property/global poissonsRatio peratomtype 0.45
+fix 		m3 all property/global coefficientRestitution peratomtypepair 1 0.3
+fix 		m4 all property/global coefficientFriction peratomtypepair 1 0.5
+fix 		m5 all property/global characteristicVelocity scalar 1.
+
+#pair style
+pair_style 	gran/hooke/history  #Hertzian without cohesion
+pair_coeff	* *
+
+#timestep, gravity
+timestep	0.00001
+fix		gravi all gravity 9.81 vector 0.0 0.0 -1.0
+
+#walls
+fix			zwalls1 all wall/gran/hooke/history     primitive type 1  zplane 0.0
+fix			zwalls2 all wall/gran/hooke/history     primitive type 1  zplane 0.0553
+fix         cylwalls all wall/gran/hooke/history    primitive type 1 zcylinder 0.01385 0. 0. tangential_damping off
+
+# change the particles density
+set             group all density 7000
+
+#cfd coupling
+fix		cfd all couple/cfd couple_every 100 mpi
+fix		cfd2 all couple/cfd/force
+
+#apply nve integration to all particles that are inserted as single particles
+fix		integr all nve/sphere
+
+#center of mass
+compute         centerOfMass all com
+
+#compute total dragforce
+#compute		dragtotal all reduce sum f_dragforce[1] f_dragforce[2] f_dragforce[3]
+
+#screen output
+compute		1 all erotate/sphere
+thermo_style	custom step atoms ke c_1 vol c_centerOfMass[3] #c_dragtotal[1] c_dragtotal[2] c_dragtotal[3]
+thermo		10
+thermo_modify	lost ignore norm no
+compute_modify	thermo_temp dynamic yes
+
+#insert the first particles so that dump is not empty
+dump 		myDump all stl 1 post/dump_*.stl
+#run		1
+dump		dmp all custom 5000 ../DEM/post/dump*.liggghts_restart id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius 
+undump		myDump
+
+run		1
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/liggghts.restart b/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/liggghts.restart
new file mode 100644
index 00000000..a15c95f0
Binary files /dev/null and b/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/liggghts.restart differ
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/post/dummy b/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/post/dummy
new file mode 100644
index 00000000..ebf63dfd
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/post/dummy
@@ -0,0 +1 @@
+dummyfile
diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/parCFDDEMrun.sh b/tutorials/cfdemSolverPiso/ErgunTestCG/parCFDDEMrun.sh
new file mode 100644
index 00000000..4bc31e5a
--- /dev/null
+++ b/tutorials/cfdemSolverPiso/ErgunTestCG/parCFDDEMrun.sh
@@ -0,0 +1,101 @@
+#!/bin/bash
+
+#===================================================================#
+# allrun script for testcase as part of test routine 
+# run settlingTest CFD part
+# Christoph Goniva - Feb. 2011
+#===================================================================#
+
+#- source CFDEM env vars
+. ~/.bashrc
+
+#- include functions
+source $CFDEM_SRC_DIR/etc/functions.sh
+
+#--------------------------------------------------------------------------------#
+#- define variables
+casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")"
+logpath=$casePath
+headerText="run_parallel_cfdemSolverPiso_ErgunTestMPI_CFDDEM"
+logfileName="log_$headerText"
+solverName="cfdemSolverPiso"
+nrProcs="4"
+machineFileName="none"   # yourMachinefileName | none
+debugMode="off"          # on | off | strict
+testHarnessPath="$CFDEM_TEST_HARNESS_PATH"
+runOctave="true"
+postproc="false"
+#--------------------------------------------------------------------------------#
+
+#- call function to run a parallel CFD-DEM case
+parCFDDEMrun $logpath $logfileName $casePath $headerText $solverName $nrProcs $machineFileName $debugMode
+
+
+if [ $runOctave == "true" ]
+    then
+        #------------------------------#
+        #  octave
+
+        #- change path
+        cd octave
+
+        #- rmove old graph
+        rm cfdemSolverPiso_ErgunTestMPI.eps
+
+        #- run octave
+        octave totalPressureDrop.m
+
+        #- show plot 
+        evince cfdemSolverPiso_ErgunTestMPI.eps
+
+        #- copy log file to test harness
+        cp ../../$logfileName $testHarnessPath
+        cp cfdemSolverPiso_ErgunTestMPI.eps $testHarnessPath
+fi
+
+if [ $postproc == "true" ]
+  then
+
+    #- keep terminal open (if started in new terminal)
+    echo "simulation finisehd? ...press enter to proceed"
+    read
+
+    #- get VTK data from liggghts dump file
+    cd $casePath/DEM/post
+    python -i $CFDEM_LPP_DIR/lpp.py dump*.liggghts_restart
+
+    #- get VTK data from CFD sim
+    cd $casePath/CFD
+    foamToVTK                                                   #- serial run of foamToVTK
+    #source $CFDEM_SRC_DIR/etc/functions.sh                       #- include functions
+    #pseudoParallelRun "foamToVTK" $nrPostProcProcessors          #- pseudo parallel run of foamToVTK
+
+    #- start paraview
+    paraview
+
+    #- keep terminal open (if started in new terminal)
+    echo "...press enter to clean up case"
+    echo "press Ctr+C to keep data"
+    read
+
+fi
+
+#- clean up case
+echo "deleting data at: $casePath :\n"
+rm -r $casePath/CFD/0.*
+rm -r $casePath/CFD/log.*
+rm -r $casePath/CFD/octave/octave-core
+rm -r $casePath/CFD/VTK
+rm -r $casePath/CFD/processor*
+rm -r $casePath/CFD/couplingFiles/*
+rm -r $casePath/DEM/post/*
+rm -r $casePath/DEM/log.*
+rm -r $casePath/DEM/liggghts.restartCFDEM*
+rm -r $casePath/CFD/probes*
+rm -r $casePath/CFD/particles
+rm -r $casePath/CFD/lagrangian
+rm -r $casePath/CFD/clockData
+echo "done"
+
+#- preserve post directory
+echo "dummyfile" >> $casePath/DEM/post/dummy