update utilities-list.txt

Added missing object files back to lagrangian library. Otherwise solvers
won't compile. If this removal was intentional, it was not complete.
Reverting the change for now.

.gitignore

@@ -0,0 +1,8 @@
+*.o
+*.d
+*.a
+*.dep
+log_*
+log.*
+*~
+

README

@@ -25,16 +25,17 @@ License
     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).
+    This code provides models and solvers 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 
+CFDEM(R) coupling provides an open source parallel coupled CFD-DEM framework 
+combining the strengths of LIGGGHTS(R) DEM code and the Open Source 
+CFD package OpenFOAM(R)(*). The CFDEM(R)coupling 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 
+code LIGGGHTS(R). 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.
@@ -54,7 +55,7 @@ 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
+- "doc" directory including the documentation of CFDEM(R)coupling
 - "tutorials" directory including basic tutorial cases showing the functionality
 
 
@@ -64,18 +65,17 @@ 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) -
+CFDEM(R)coupling stands for Computational Fluid Dynamics (CFD) -
 Discrete Element Method (DEM) coupling.
 
-CFDEMcoupling is an open-source code, distributed freely under the terms of the 
+CFDEM(R)coupling is an open-source code, distributed freely under the terms of the 
 GNU Public License (GPL).
 
-Core development of CFDEMcoupling is done by 
+Core development of CFDEM(R)coupling 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.
+(*) "OpenFOAM(R)"_of is a registered trade mark of OpenCFD Limited, a wholly owned subsidiary of the ESI Group.
+This offering is not approved or endorsed by OpenCFD Limited, the producer of the OpenFOAM software and owner of the OPENFOAM® and OpenCFD® trade marks.
 \*---------------------------------------------------------------------------*/

applications/.gitignore

@@ -0,0 +1,7 @@
+*.o
+*.d
+*.a
+*.dep
+log_*
+log.*
+*~

applications/solvers/cfdemSolverIB/Make/options

@@ -11,9 +11,8 @@ EXE_INC = \
     -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/fvOptions/lnInclude
-
 EXE_LIBS = \
-    -L$(FOAM_USER_LIBBIN)\
+    -L$(CFDEM_LIB_DIR)\
     -lincompressibleRASModels \
     -lincompressibleLESModels \
     -lincompressibleTransportModels \
@@ -22,3 +21,4 @@ EXE_LIBS = \
     -ldynamicMesh \
     -lfvOptions \
     -l$(CFDEM_LIB_NAME)
+

applications/solvers/cfdemSolverIB/cfdemSolverIB.C

@@ -72,7 +72,7 @@ int main(int argc, char *argv[])
     #include "initContinuityErrs.H"
 
     #if defined(version22)
-        #include "createFvOptions.H"
+    #include "createFvOptions.H"
     #endif
 
     // create cfdemCloud
@@ -128,10 +128,12 @@ int main(int argc, char *argv[])
             for (int corr=0; corr<nCorr; corr++)
             {
                 volScalarField rUA = 1.0/UEqn.A();
+                surfaceScalarField rUAf(fvc::interpolate(rUA));
 
                 U = rUA*UEqn.H();
                 #ifdef version23
-                phi = (fvc::interpolate(U) & mesh.Sf()); // there is a new version in 23x
+                phi = (fvc::interpolate(U) & mesh.Sf())
+                    + rUAf*fvc::ddtCorr(U, phi);
                 #else
                 phi = (fvc::interpolate(U) & mesh.Sf())
                     + fvc::ddtPhiCorr(rUA, U, phi);

applications/solvers/cfdemSolverPiso/Make/options

@@ -7,7 +7,7 @@ EXE_INC = \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
 
 EXE_LIBS = \
-    -L$(FOAM_USER_LIBBIN)\
+    -L$(CFDEM_LIB_DIR)\
     -lincompressibleRASModels \
     -lincompressibleLESModels \
     -lincompressibleTransportModels \

applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C

@@ -42,6 +42,7 @@ Description
 #include "implicitCouple.H"
 #include "clockModel.H"
 #include "smoothingModel.H"
+#include "forceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,11 +60,10 @@ int main(int argc, char *argv[])
     #include "checkModelType.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
+    Info<< "\nStarting time loop\n" << endl;
     while (runTime.loop())
     {
-        Info<< "\nStarting time loop\n" << endl;
-            particleCloud.clockM().start(1,"Global");
+        particleCloud.clockM().start(1,"Global");
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
@@ -72,13 +72,24 @@ int main(int argc, char *argv[])
 
         // do particle stuff
         particleCloud.clockM().start(2,"Coupling");
-        particleCloud.evolve(voidfraction,Us,U);
+        bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
+
+        if(hasEvolved)
+        {
+            particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
+        }
     
         Info << "update Ksl.internalField()" << endl;
         Ksl = particleCloud.momCoupleM(0).impMomSource();
-        particleCloud.smoothingM().smoothen(Ksl);
         Ksl.correctBoundaryConditions();
 
+       //Force Checks
+       vector fTotal(0,0,0);
+       vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField()));
+       reduce(fImpTotal, sumOp<vector>());
+       Info << "TotalForceExp: " << fTotal << endl;
+       Info << "TotalForceImp: " << fImpTotal << endl;
+
         #include "solverDebugInfo.H"
         particleCloud.clockM().stop("Coupling");
 
@@ -91,23 +102,19 @@ int main(int argc, char *argv[])
                 // Momentum predictor
                 fvVectorMatrix UEqn
                 (
-                    fvm::ddt(voidfraction,U) + fvm::Sp(fvc::ddt(voidfraction),U)
-                  + fvm::div(phi,U) + fvm::Sp(fvc::div(phi),U)
+                    fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
+                  + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
 //                + turbulence->divDevReff(U)
                   + particleCloud.divVoidfractionTau(U, voidfraction)
                   ==
                   - fvm::Sp(Ksl/rho,U)
                 );
 
-                if (modelType=="B")
-                 UEqn == - fvc::grad(p) + Ksl/rho*Us;
-                else
-                    UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us;
-
                 UEqn.relax();
-
-                if (momentumPredictor)
-                    solve(UEqn);
+                if (momentumPredictor && (modelType=="B" || modelType=="Bfull"))
+                    solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
+                else if (momentumPredictor)
+                    solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
 
                 // --- PISO loop
 
@@ -120,11 +127,17 @@ int main(int argc, char *argv[])
 
                     surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
                     volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
+                    surfaceScalarField rUAfvoidfraction("(voidfraction2|A(U)F)", fvc::interpolate(rUAvoidfraction));
 
                     U = rUA*UEqn.H();
 
-                    phi = (fvc::interpolate(U*voidfraction) & mesh.Sf() );
-                     //+ fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
+                    #ifdef version23
+                    phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
+                        + rUAfvoidfraction*fvc::ddtCorr(U, phi);
+                    #else
+                    phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
+                        + fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
+                    #endif
                     surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
                     surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
 
@@ -157,13 +170,14 @@ int main(int argc, char *argv[])
                         if (nonOrth == nNonOrthCorr)
                         {
                             phiGes -= pEqn.flux();
+                            phi = phiGes;
                         }
 
                     } // end non-orthogonal corrector loop
 
                     #include "continuityErrorPhiPU.H"
 
-                    if (modelType=="B")
+                    if (modelType=="B" || modelType=="Bfull")
                         U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
                     else
                         U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
@@ -191,7 +205,7 @@ int main(int argc, char *argv[])
     }
 
     Info<< "End\n" << endl;
-    
+
     return 0;
 }
 

applications/solvers/cfdemSolverPisoMS/Make/options

@@ -8,7 +8,7 @@ EXE_INC = \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
 
 EXE_LIBS = \
-    -L$(FOAM_USER_LIBBIN)\
+    -L$(CFDEM_LIB_DIR)\
     -lincompressibleRASModels \
     -lincompressibleLESModels \
     -lincompressibleTransportModels \

applications/solvers/cfdemSolverPisoScalar/Make/options

@@ -7,7 +7,7 @@ EXE_INC = \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
 
 EXE_LIBS = \
-    -L$(FOAM_USER_LIBBIN)\
+    -L$(CFDEM_LIB_DIR)\
     -lincompressibleRASModels \
     -lincompressibleLESModels \
     -lincompressibleTransportModels \

applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C

@@ -30,7 +30,7 @@ Application
 Description
     Transient solver for incompressible flow.
     Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
-    The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6, 
+    The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6,
     where additional functionality for CFD-DEM coupling is added.
 \*---------------------------------------------------------------------------*/
 
@@ -40,30 +40,26 @@ Description
 
 #include "cfdemCloud.H"
 #include "implicitCouple.H"
-#include "forceModel.H"
 #include "smoothingModel.H"
+#include "forceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
     #include "setRootCase.H"
-
     #include "createTime.H"
     #include "createMesh.H"
     #include "createFields.H"
-
     #include "initContinuityErrs.H"
 
     // create cfdemCloud
     #include "readGravitationalAcceleration.H"
     cfdemCloud particleCloud(mesh);
-
     #include "checkModelType.H"
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
     Info<< "\nStarting time loop\n" << endl;
-
     while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
@@ -72,11 +68,15 @@ int main(int argc, char *argv[])
         #include "CourantNo.H"
 
         // do particle stuff
-        Info << "- evolve()" << endl;
-        particleCloud.evolve(voidfraction,Us,U);
+        bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
 
-        Ksl.internalField() = particleCloud.momCoupleM(0).impMomSource();
-        particleCloud.smoothingM().smoothen(Ksl);
+        if(hasEvolved)
+        {
+            particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
+        }
+    
+        Info << "update Ksl.internalField()" << endl;
+        Ksl = particleCloud.momCoupleM(0).impMomSource();
         Ksl.correctBoundaryConditions();
 
 
@@ -87,107 +87,116 @@ int main(int argc, char *argv[])
         Tsource.correctBoundaryConditions();
 
         // solve scalar transport equation
-        phi = fvc::interpolate(U*voidfraction) & mesh.Sf();
-
-        solve
+        fvScalarMatrix TEqn
         (
-            fvm::ddt(voidfraction,T)
-          + fvm::div(phi, T)
-          - fvm::laplacian(DT*voidfraction, T)
-          ==
-            Tsource
+           fvm::ddt(voidfraction,T) - fvm::Sp(fvc::ddt(voidfraction),T)
+         + fvm::div(phi, T) - fvm::Sp(fvc::div(phi),T)
+         - fvm::laplacian(DT*voidfraction, T)
+         ==
+           Tsource
         );
+        TEqn.relax();
+        TEqn.solve();
 
-        // Pressure-velocity PISO corrector
+        if(particleCloud.solveFlow())
         {
-            // Momentum predictor
-            fvVectorMatrix UEqn
-            (
-                fvm::ddt(voidfraction,U)
-              + fvm::div(phi, U)
-              + turbulence->divDevReff(U)
-             == 
-              - fvm::Sp(Ksl/rho,U)
-            );
-
-            UEqn.relax();
-
-            if (momentumPredictor)
+            // Pressure-velocity PISO corrector
             {
-                //solve UEqn
-                if (modelType=="B")
+                // Momentum predictor
+                fvVectorMatrix UEqn
+                (
+                    fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
+                  + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
+//                + turbulence->divDevReff(U)
+                  + particleCloud.divVoidfractionTau(U, voidfraction)
+                  ==
+                  - fvm::Sp(Ksl/rho,U)
+                );
+
+                UEqn.relax();
+                if (momentumPredictor && (modelType=="B" || modelType=="Bfull"))
                     solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
-                else
+                else if (momentumPredictor)
                     solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
+
+                // --- PISO loop
+
+                //for (int corr=0; corr<nCorr; corr++)
+                int nCorrSoph = nCorr + 5 * pow((1-particleCloud.dataExchangeM().timeStepFraction()),1);
+
+                for (int corr=0; corr<nCorrSoph; corr++)
+                {
+                    volScalarField rUA = 1.0/UEqn.A();
+
+                    surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
+                    volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
+                    surfaceScalarField rUAfvoidfraction("(voidfraction2|A(U)F)", fvc::interpolate(rUAvoidfraction));
+
+                    U = rUA*UEqn.H();
+
+                    #ifdef version23
+                    phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
+                        + rUAfvoidfraction*fvc::ddtCorr(U, phi);
+                    #else
+                    phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
+                        + fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
+                    #endif
+                    surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
+                    surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
+
+                    if (modelType=="A")
+                        rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
+
+                    // Non-orthogonal pressure corrector loop
+                    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+                    {
+                        // Pressure corrector
+                        fvScalarMatrix pEqn
+                        (
+                            fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + particleCloud.ddtVoidfraction()
+                        );
+                        pEqn.setReference(pRefCell, pRefValue);
+
+                        if
+                        (
+                            corr == nCorr-1
+                         && nonOrth == nNonOrthCorr
+                        )
+                        {
+                            pEqn.solve(mesh.solver("pFinal"));
+                        }
+                        else
+                        {
+                            pEqn.solve();
+                        }
+
+                        if (nonOrth == nNonOrthCorr)
+                        {
+                            phiGes -= pEqn.flux();
+                            phi = phiGes;
+                        }
+
+                    } // end non-orthogonal corrector loop
+
+                    #include "continuityErrorPhiPU.H"
+
+                    if (modelType=="B" || modelType=="Bfull")
+                        U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
+                    else
+                        U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
+
+                    U.correctBoundaryConditions();
+
+                } // end piso loop
             }
 
-            // --- PISO loop
-
-            //for (int corr=0; corr<nCorr; corr++)
-            int nCorrSoph = nCorr + 5 * pow((1-particleCloud.dataExchangeM().timeStepFraction()),1);
-            for (int corr=0; corr<nCorrSoph; corr++)
-            {
-                volScalarField rUA = 1.0/UEqn.A();
-                surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
-                volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
-
-                U = rUA*UEqn.H();
-
-                #ifdef version23
-                phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() );
-                #else
-                phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
-                      + fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
-                #endif
-                surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
-                surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
-                if (modelType=="A")
-                    rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
-
-                // Non-orthogonal pressure corrector loop
-                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-                {
-                    // Pressure corrector
-                    fvScalarMatrix pEqn
-                    (
-                        fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + particleCloud.ddtVoidfraction()
-                    );
-                    pEqn.setReference(pRefCell, pRefValue);
-
-                    if
-                    (
-                        corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
-                    )
-                    {
-                        pEqn.solve(mesh.solver("pFinal"));
-                    }
-                    else
-                    {
-                        pEqn.solve();
-                    }
-
-                    if (nonOrth == nNonOrthCorr)
-                    {
-                        phiGes -= pEqn.flux();
-                    }
-
-                } // end non-orthogonal corrector loop
-
-                #include "continuityErrorPhiPU.H"
-
-                if (modelType=="B")
-                    U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-                else
-                    U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-
-                U.correctBoundaryConditions();
-
-            } // end piso loop
+            turbulence->correct();
+        }// end solveFlow
+        else
+        {
+            Info << "skipping flow solution." << endl;
         }
 
-        turbulence->correct();
-
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/utilities/cfdemPostproc/Make/options

@@ -8,7 +8,7 @@ EXE_INC = \
 
 
 EXE_LIBS = \
-    -L$(FOAM_USER_LIBBIN)\
+    -L$(CFDEM_LIB_DIR)\
     -lincompressibleRASModels \
     -lincompressibleLESModels \
     -lincompressibleTransportModels \

applications/utilities/cfdemPostproc/cfdemPostproc.C

@@ -71,6 +71,7 @@ int main(int argc, char *argv[])
     double **voidfractions_;
     double **particleWeights_;
     double **particleVolumes_;
+    double **particleV_;
     double **cellIDs_;
     
     particleCloud.dataExchangeM().allocateArray(positions_,0.,3);
@@ -80,6 +81,7 @@ int main(int argc, char *argv[])
     particleCloud.dataExchangeM().allocateArray(voidfractions_,0.,1);
     particleCloud.dataExchangeM().allocateArray(particleWeights_,0.,1);
     particleCloud.dataExchangeM().allocateArray(particleVolumes_,0.,1);
+    particleCloud.dataExchangeM().allocateArray(particleV_,0.,1);
     particleCloud.get_cellIDs(cellIDs_);  // get ref to cellIDs
     //particleCloud.dataExchangeM().allocateArray(cellIDs_,0.,1);
     
@@ -105,7 +107,7 @@ int main(int argc, char *argv[])
         particleCloud.locateM().findCell(NULL,positions_,cellIDs_,particleCloud.numberOfParticles());
         particleCloud.setPos(positions_);
 
-        particleCloud.voidFractionM().setvoidFraction(NULL,voidfractions_,particleWeights_,particleVolumes_);
+        particleCloud.voidFractionM().setvoidFraction(NULL,voidfractions_,particleWeights_,particleVolumes_,particleV_);
 
         voidfraction.internalField() = particleCloud.voidFractionM().voidFractionInterp();
         voidfraction.correctBoundaryConditions();
@@ -135,6 +137,7 @@ int main(int argc, char *argv[])
     particleCloud.dataExchangeM().destroy(voidfractions_,1);
     particleCloud.dataExchangeM().destroy(particleWeights_,1);
     particleCloud.dataExchangeM().destroy(particleVolumes_,1);
+    particleCloud.dataExchangeM().destroy(particleV_,1);
     //particleCloud.dataExchangeM().destroy(cellIDs_); // destroyed in cloud
 
     Info<< "End\n" << endl;

doc/.gitignore

@@ -0,0 +1,7 @@
+*.o
+*.d
+*.a
+*.dep
+log_*
+log.*
+*~

doc/CFDEMcoupling_Manual.html

@@ -164,7 +164,7 @@ In order to get the latest code version, please use the git repository at http:/
 </P>
 <PRE>modelType 
 </PRE>
-<P>"modelType" refers to the formulation of the equations to be solved. Choose "A" or "B", according to Zhou et al. (2010): "Discrete particle simulation of particle-fluid flow: model formulations and their applicability", JFM. "A" requires the use of the force models gradPForce and viscForce, whereas "B" requires the force model "Archimedes".
+<P>"modelType" refers to the formulation of the equations to be solved. Choose "A", "B" or "Bfull", according to Zhou et al. (2010): "Discrete particle simulation of particle-fluid flow: model formulations and their applicability", JFM. "A" requires the use of the force models gradPForce and viscForce, whereas "B" requires the force model "Archimedes". "Bfull" refers to model type I.
 </P>
 <PRE>couplingInterval 
 </PRE>
@@ -212,30 +212,31 @@ listing below of styles within certain commands.
 <TR ALIGN="center"><TD ><A HREF = "clockModel_noClock.html">clockModel_noClock</A></TD><TD ><A HREF = "clockModel_standardClock.html">clockModel_standardClock</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "dataExchangeModel.html">dataExchangeModel</A></TD><TD ><A HREF = "dataExchangeModel_noDataExchange.html">dataExchangeModel_noDataExchange</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "dataExchangeModel_oneWayVTK.html">dataExchangeModel_oneWayVTK</A></TD><TD ><A HREF = "dataExchangeModel_twoWayFiles.html">dataExchangeModel_twoWayFiles</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "dataExchangeModel_twoWayM2M.html">dataExchangeModel_twoWayM2M</A></TD><TD ><A HREF = "dataExchangeModel_twoWayMPI.html">dataExchangeModel_twoWayMPI</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "dataExchangeModel_twoWayMPI.html">dataExchangeModel_twoWayMPI</A></TD><TD ><A HREF = "dataExchangeModel_twoWayMany2Many.html">dataExchangeModel_twoWayMany2Many</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "forceModel.html">forceModel</A></TD><TD ><A HREF = "forceModel_Archimedes.html">forceModel_Archimedes</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "forceModel_ArchimedesIB.html">forceModel_ArchimedesIB</A></TD><TD ><A HREF = "forceModel_DiFeliceDrag.html">forceModel_DiFeliceDrag</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "forceModel_GidaspowDrag.html">forceModel_GidaspowDrag</A></TD><TD ><A HREF = "forceModel_KochHillDrag.html">forceModel_KochHillDrag</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "forceModel_LaEuScalarTemp.html">forceModel_LaEuScalarTemp</A></TD><TD ><A HREF = "forceModel_MeiLift.html">forceModel_MeiLift</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "forceModel_SchillerNaumannDrag.html">forceModel_SchillerNaumannDrag</A></TD><TD ><A HREF = "forceModel_ShirgaonkarIB.html">forceModel_ShirgaonkarIB</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "forceModel_gradPForce.html">forceModel_gradPForce</A></TD><TD ><A HREF = "forceModel_noDrag.html">forceModel_noDrag</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "forceModel_particleCellVolume.html">forceModel_particleCellVolume</A></TD><TD ><A HREF = "forceModel_virtualMassForce.html">forceModel_virtualMassForce</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "forceModel_viscForce.html">forceModel_viscForce</A></TD><TD ><A HREF = "liggghtsCommandModel.html">liggghtsCommandModel</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "forceModel_fieldStore.html">forceModel_fieldStore</A></TD><TD ><A HREF = "forceModel_gradPForce.html">forceModel_gradPForce</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "forceModel_noDrag.html">forceModel_noDrag</A></TD><TD ><A HREF = "forceModel_particleCellVolume.html">forceModel_particleCellVolume</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "forceModel_virtualMassForce.html">forceModel_virtualMassForce</A></TD><TD ><A HREF = "forceModel_viscForce.html">forceModel_viscForce</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "forceSubModel.html">forceSubModel</A></TD><TD ><A HREF = "forceSubModel_ImEx.html">forceSubModel_ImEx</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "forceSubModel_ImExCorr.html">forceSubModel_ImExCorr</A></TD><TD ><A HREF = "liggghtsCommandModel.html">liggghtsCommandModel</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "liggghtsCommandModel_execute.html">liggghtsCommandModel_execute</A></TD><TD ><A HREF = "liggghtsCommandModel_readLiggghtsData.html">liggghtsCommandModel_readLiggghtsData</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "liggghtsCommandModel_runLiggghts.html">liggghtsCommandModel_runLiggghts</A></TD><TD ><A HREF = "liggghtsCommandModel_writeLiggghts.html">liggghtsCommandModel_writeLiggghts</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "locateModel.html">locateModel</A></TD><TD ><A HREF = "locateModel_engineSearch.html">locateModel_engineSearch</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "locateModel_engineSearchIB.html">locateModel_engineSearchIB</A></TD><TD ><A HREF = "locateModel_standardSearch.html">locateModel_standardSearch</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "locateModel_turboEngineSearchM2M.html">locateModel_turboEngineM2MSearch</A></TD><TD ><A HREF = "locateModel_turboEngineSearch.html">locateModel_turboEngineSearch</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "meshMotionModel.html">meshMotionModel</A></TD><TD ><A HREF = "meshMotionModel_noMeshMotion.html">meshMotionModel_noMeshMotion</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "momCoupleModel.html">momCoupleModel</A></TD><TD ><A HREF = "momCoupleModel_explicitCouple.html">momCoupleModel_explicitCouple</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "momCoupleModel_implicitCouple.html">momCoupleModel_implicitCouple</A></TD><TD ><A HREF = "momCoupleModel_noCouple.html">momCoupleModel_noCouple</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "probeModel.html">probeModel</A></TD><TD ><A HREF = "probeModel_noProbe.html">probeModel_noProbe</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "regionModel.html">regionModel</A></TD><TD ><A HREF = "regionModel_allRegion.html">regionModel_allRegion</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "smoothingModel.html">smoothingModel</A></TD><TD ><A HREF = "smoothingModel_constDiffSmoothing.html">smoothingModel_constDiffSmoothing</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "smoothingModel_noSmoothing.html">smoothingModel_noSmoothing</A></TD><TD ><A HREF = "voidFractionModel.html">voidfractionModel</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "voidFractionModel_GaussVoidFraction.html">voidfractionModel_GaussVoidFraction</A></TD><TD ><A HREF = "voidFractionModel_IBVoidFraction.html">voidfractionModel_IBVoidFraction</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "voidFractionModel_bigParticleVoidFraction.html">voidfractionModel_bigParticleVoidFraction</A></TD><TD ><A HREF = "voidFractionModel_centreVoidFraction.html">voidfractionModel_centreVoidFraction</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "voidFractionModel_dividedVoidFraction.html">voidfractionModel_dividedVoidFraction</A> 
+<TR ALIGN="center"><TD ><A HREF = "locateModel_turboEngineSearch.html">locateModel_turboEngineSearch</A></TD><TD ><A HREF = "meshMotionModel.html">meshMotionModel</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "meshMotionModel_noMeshMotion.html">meshMotionModel_noMeshMotion</A></TD><TD ><A HREF = "momCoupleModel.html">momCoupleModel</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "momCoupleModel_explicitCouple.html">momCoupleModel_explicitCouple</A></TD><TD ><A HREF = "momCoupleModel_implicitCouple.html">momCoupleModel_implicitCouple</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "momCoupleModel_noCouple.html">momCoupleModel_noCouple</A></TD><TD ><A HREF = "probeModel.html">probeModel</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "probeModel_noProbe.html">probeModel_noProbe</A></TD><TD ><A HREF = "regionModel.html">regionModel</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "regionModel_allRegion.html">regionModel_allRegion</A></TD><TD ><A HREF = "smoothingModel.html">smoothingModel</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "smoothingModel_constDiffSmoothing.html">smoothingModel_constDiffSmoothing</A></TD><TD ><A HREF = "smoothingModel_noSmoothing.html">smoothingModel_noSmoothing</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "voidFractionModel.html">voidfractionModel</A></TD><TD ><A HREF = "voidFractionModel_GaussVoidFraction.html">voidfractionModel_GaussVoidFraction</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "voidFractionModel_IBVoidFraction.html">voidfractionModel_IBVoidFraction</A></TD><TD ><A HREF = "voidFractionModel_bigParticleVoidFraction.html">voidfractionModel_bigParticleVoidFraction</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "voidFractionModel_centreVoidFraction.html">voidfractionModel_centreVoidFraction</A></TD><TD ><A HREF = "voidFractionModel_dividedVoidFraction.html">voidfractionModel_dividedVoidFraction</A> 
 </TD></TR></TABLE></DIV>
 
 </HTML>

doc/CFDEMcoupling_Manual.txt

@@ -133,7 +133,7 @@ Reasonable example settings for the "couplingProperties" dictionary are given in
 
 modelType :pre
 
-"modelType" refers to the formulation of the equations to be solved. Choose "A" or "B", according to Zhou et al. (2010): "Discrete particle simulation of particle-fluid flow: model formulations and their applicability", JFM. "A" requires the use of the force models gradPForce and viscForce, whereas "B" requires the force model "Archimedes".
+"modelType" refers to the formulation of the equations to be solved. Choose "A", "B" or "Bfull", according to Zhou et al. (2010): "Discrete particle simulation of particle-fluid flow: model formulations and their applicability", JFM. "A" requires the use of the force models gradPForce and viscForce, whereas "B" requires the force model "Archimedes". "Bfull" refers to model type I.
 
 couplingInterval :pre
 
@@ -221,6 +221,10 @@ Reasonable example settings for the "liggghtsCommands" dictionary are given in t
 
 
 
+
+
+
+
 
 
 
@@ -250,8 +254,8 @@ listing below of styles within certain commands.
 "dataExchangeModel_noDataExchange"_dataExchangeModel_noDataExchange.html,
 "dataExchangeModel_oneWayVTK"_dataExchangeModel_oneWayVTK.html,
 "dataExchangeModel_twoWayFiles"_dataExchangeModel_twoWayFiles.html,
-"dataExchangeModel_twoWayM2M"_dataExchangeModel_twoWayM2M.html,
 "dataExchangeModel_twoWayMPI"_dataExchangeModel_twoWayMPI.html,
+"dataExchangeModel_twoWayMany2Many"_dataExchangeModel_twoWayMany2Many.html,
 "forceModel"_forceModel.html,
 "forceModel_Archimedes"_forceModel_Archimedes.html,
 "forceModel_ArchimedesIB"_forceModel_ArchimedesIB.html,
@@ -262,11 +266,15 @@ listing below of styles within certain commands.
 "forceModel_MeiLift"_forceModel_MeiLift.html,
 "forceModel_SchillerNaumannDrag"_forceModel_SchillerNaumannDrag.html,
 "forceModel_ShirgaonkarIB"_forceModel_ShirgaonkarIB.html,
+"forceModel_fieldStore"_forceModel_fieldStore.html,
 "forceModel_gradPForce"_forceModel_gradPForce.html,
 "forceModel_noDrag"_forceModel_noDrag.html,
 "forceModel_particleCellVolume"_forceModel_particleCellVolume.html,
 "forceModel_virtualMassForce"_forceModel_virtualMassForce.html,
 "forceModel_viscForce"_forceModel_viscForce.html,
+"forceSubModel"_forceSubModel.html,
+"forceSubModel_ImEx"_forceSubModel_ImEx.html,
+"forceSubModel_ImExCorr"_forceSubModel_ImExCorr.html,
 "liggghtsCommandModel"_liggghtsCommandModel.html,
 "liggghtsCommandModel_execute"_liggghtsCommandModel_execute.html,
 "liggghtsCommandModel_readLiggghtsData"_liggghtsCommandModel_readLiggghtsData.html,
@@ -276,7 +284,6 @@ listing below of styles within certain commands.
 "locateModel_engineSearch"_locateModel_engineSearch.html,
 "locateModel_engineSearchIB"_locateModel_engineSearchIB.html,
 "locateModel_standardSearch"_locateModel_standardSearch.html,
-"locateModel_turboEngineM2MSearch"_locateModel_turboEngineSearchM2M.html,
 "locateModel_turboEngineSearch"_locateModel_turboEngineSearch.html,
 "meshMotionModel"_meshMotionModel.html,
 "meshMotionModel_noMeshMotion"_meshMotionModel_noMeshMotion.html,

doc/cfdemSolverPisoScalar.html

@@ -11,12 +11,14 @@
 </H3>
 <P><B>Description:</B>
 </P>
-<P>"cfdemSolverPisoScalar" is a coupled CFD-DEM solver using CFDEMcoupling, an open source parallel coupled CFD-DEM framework. Based on pisoFoam(R)(*), a finite volume based solver for turbulent Navier-Stokes equations applying PISO algorithm, "cfdemSolverPisoScalar" has additional functionality for a coupling to the DEM code "LIGGGHTS" as well as a scalar transport equation. The volume averaged Navier-Stokes Equations are solved accounting for momentum exchange and volume displacement of discrete particles whose trajectories are calculated in the DEM code LIGGGHTS. The scalar transport equation is coupled to scalar properties of the particle phase, thus convective heat transfer in a fluid granular system can be modeled with "cfdemSolverPisoScalar".
+<P>"cfdemSolverPisoScalar" is a coupled CFD-DEM solver using CFDEMcoupling, an open source parallel coupled CFD-DEM framework. Based on pisoFoam(R)(*), a finite volume based solver for turbulent Navier-Stokes equations applying PISO algorithm, "cfdemSolverPisoScalar" has additional functionality for a coupling to the DEM code "LIGGGHTS" as well as a scalar transport equation. The volume averaged Navier-Stokes Equations are solved accounting for momentum exchange and volume displacement of discrete particles, whose trajectories are calculated in the DEM code LIGGGHTS. The scalar transport equation is coupled to scalar properties of the particle phase, thus convective heat transfer in a fluid granular system can be modeled with "cfdemSolverPisoScalar". 
 </P>
 <P>see:
 </P>
 <P>GONIVA, C., KLOSS, C., HAGER,A. and PIRKER, S. (2010): "An Open Source CFD-DEM Perspective", Proc. of OpenFOAM Workshop, Göteborg, June 22.-24.
 </P>
+<P>The heat transfer equation is implemented according to Nield & Bejan (2013), Convection in Porous Media, DOI 10.1007/978-1-4614-5541-7_2, Springer
+</P>
 <HR>
 
 <P>(*) This offering is not approved or endorsed by OpenCFD Limited, the producer of the OpenFOAM software and owner of the OPENFOAM®  and OpenCFD®  trade marks. OPENFOAM® is a registered trade mark of OpenCFD Limited, a wholly owned subsidiary of the ESI Group.

doc/cfdemSolverPisoScalar.txt

@@ -9,12 +9,15 @@ cfdemSolverPisoScalar command :h3
 
 [Description:]
 
-"cfdemSolverPisoScalar" is a coupled CFD-DEM solver using CFDEMcoupling, an open source parallel coupled CFD-DEM framework. Based on pisoFoam(R)(*), a finite volume based solver for turbulent Navier-Stokes equations applying PISO algorithm, "cfdemSolverPisoScalar" has additional functionality for a coupling to the DEM code "LIGGGHTS" as well as a scalar transport equation. The volume averaged Navier-Stokes Equations are solved accounting for momentum exchange and volume displacement of discrete particles whose trajectories are calculated in the DEM code LIGGGHTS. The scalar transport equation is coupled to scalar properties of the particle phase, thus convective heat transfer in a fluid granular system can be modeled with "cfdemSolverPisoScalar".
+"cfdemSolverPisoScalar" is a coupled CFD-DEM solver using CFDEMcoupling, an open source parallel coupled CFD-DEM framework. Based on pisoFoam(R)(*), a finite volume based solver for turbulent Navier-Stokes equations applying PISO algorithm, "cfdemSolverPisoScalar" has additional functionality for a coupling to the DEM code "LIGGGHTS" as well as a scalar transport equation. The volume averaged Navier-Stokes Equations are solved accounting for momentum exchange and volume displacement of discrete particles, whose trajectories are calculated in the DEM code LIGGGHTS. The scalar transport equation is coupled to scalar properties of the particle phase, thus convective heat transfer in a fluid granular system can be modeled with "cfdemSolverPisoScalar". 
 
 see:
 
 GONIVA, C., KLOSS, C., HAGER,A. and PIRKER, S. (2010): "An Open Source CFD-DEM Perspective", Proc. of OpenFOAM Workshop, Göteborg, June 22.-24.
 
+The heat transfer equation is implemented according to Nield & Bejan (2013), Convection in Porous Media, DOI 10.1007/978-1-4614-5541-7_2, Springer
+
+
 
 :line
 (*) This offering is not approved or endorsed by OpenCFD Limited, the producer of the OpenFOAM software and owner of the OPENFOAM®  and OpenCFD®  trade marks. OPENFOAM® is a registered trade mark of OpenCFD Limited, a wholly owned subsidiary of the ESI Group.

doc/clockModel_standardClock.html

@@ -23,7 +23,7 @@
 </P>
 <P>The "standardClock" model is a basic clockModel model which measures the run time between every ".start(int arrayPos,string name)" and ".stop(string name)" statement placed in the code. If a ".start(name)" is called more than once (e.g. in a loop) the accumulated times are calculated.  After the simulation has finished, the data is stored in $caseDir/CFD/clockData/$startTime/*.txt .
 Since the measurements are stored in an array, it is necessary to put a variable <I>arrayPos</I> (type integer) at the start command. Those do not need to be in ascending order and positions may be omitted. The standard size of this array is 30 and can be changed at the initialization of the standardClock class. If <I>arrayPos</I> is out of bounds, the array size will be doubled. The stop command does not need <I>arrayPos</I>, since the class remembers the positions. The string name is intended for easier evaluation afterwards an may be omitted like ".start(int arrayPos)" and ".stop()". The command ".stop(string name)" is a safety feature, because if the name is not equal to the started name, output will be produced for information.
-After the case ran you may use the matPlot.py script located in $CFDEM_UT_DIR/vizClock/ to produce a graphical output of your measurements. The usage is like 'python < matPlot.py' and you have to be in the directory of the desired time step, where there is a file called "timeEvalFull.txt", which contains averaged and maximum data with respect to the number of processes.
+After the case ran you may use the matPlot.py script located in $CFDEM_UT_DIR/vizClock/ to produce a graphical output of your measurements. The usage is like 'python < matPlot.py' and you have to be in the directory of the desired time step, where there is a file called "timeEvalFull.txt", which contains averaged and maximum data with respect to the number of processes. There is an alias called "vizClock" to run this python routine for visualizing the data.
 </P>
 <P><B>Restrictions:</B> none.
 </P>

doc/clockModel_standardClock.txt

@@ -21,7 +21,7 @@ clockModel standardClock; :pre
 
 The "standardClock" model is a basic clockModel model which measures the run time between every ".start(int arrayPos,string name)" and ".stop(string name)" statement placed in the code. If a ".start(name)" is called more than once (e.g. in a loop) the accumulated times are calculated.  After the simulation has finished, the data is stored in $caseDir/CFD/clockData/$startTime/*.txt .
 Since the measurements are stored in an array, it is necessary to put a variable {arrayPos} (type integer) at the start command. Those do not need to be in ascending order and positions may be omitted. The standard size of this array is 30 and can be changed at the initialization of the standardClock class. If {arrayPos} is out of bounds, the array size will be doubled. The stop command does not need {arrayPos}, since the class remembers the positions. The string name is intended for easier evaluation afterwards an may be omitted like ".start(int arrayPos)" and ".stop()". The command ".stop(string name)" is a safety feature, because if the name is not equal to the started name, output will be produced for information.
-After the case ran you may use the matPlot.py script located in $CFDEM_UT_DIR/vizClock/ to produce a graphical output of your measurements. The usage is like 'python < matPlot.py' and you have to be in the directory of the desired time step, where there is a file called "timeEvalFull.txt", which contains averaged and maximum data with respect to the number of processes.
+After the case ran you may use the matPlot.py script located in $CFDEM_UT_DIR/vizClock/ to produce a graphical output of your measurements. The usage is like 'python < matPlot.py' and you have to be in the directory of the desired time step, where there is a file called "timeEvalFull.txt", which contains averaged and maximum data with respect to the number of processes. There is an alias called "vizClock" to run this python routine for visualizing the data.
 
 [Restrictions:] none.
 

doc/dataExchangeModel_twoWayM2M.txt

@@ -1,42 +0,0 @@
-"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c
-
-:link(lws,http://www.cfdem.com)
-:link(lc,CFDEMcoupling_Manual.html#comm)
-
-:line
-
-dataExchangeModel_twoWayM2M command :h3
-
-[Syntax:]
-
-Defined in couplingProperties dictionary.
-
-dataExchangeModel twoWayM2M;
-twoWayM2MProps
-\{
-    liggghtsPath "path";
-\}; :pre
-
-{path} = path to the DEM simulation input file :ulb,l
-:ule
-
-[Examples:]
-
-dataExchangeModel twoWayM2M;
-twoWayM2MProps
-\{
-    liggghtsPath "../DEM/in.liggghts_init";
-\} :pre
-
-[Description:]
-
-The data exchange model performs the data exchange between the DEM code and the CFD code. The twoWayM2M model is a model that can exchange particle properties from DEM to CFD and from CFD to DEM. Data is exchanged via MPI technique using the many to many mapping scheme. The DEM run is executed by the coupling model, via a liggghtsCommandModel object.
-
-[Restrictions:]
-
-Should be used in combination with the turboEngineSearchM2M locate model to achieve best performance!
-
-[Related commands:]
-
-"dataExchangeModel"_dataExchangeModel.html
-

doc/dataExchangeModel_twoWayMany2Many.html

@@ -7,14 +7,14 @@
 
 <HR>
 
-<H3>dataExchangeModel_twoWayM2M command 
+<H3>dataExchangeModel_twoWayMany2Many command 
 </H3>
 <P><B>Syntax:</B>
 </P>
 <P>Defined in couplingProperties dictionary.
 </P>
-<PRE>dataExchangeModel twoWayM2M;
-twoWayM2MProps
+<PRE>dataExchangeModel twoWayMany2Many;
+twoWayMany2ManyProps
 {
     liggghtsPath "path";
 }; 
@@ -25,19 +25,19 @@ twoWayM2MProps
 </UL>
 <P><B>Examples:</B>
 </P>
-<PRE>dataExchangeModel twoWayM2M;
-twoWayM2MProps
+<PRE>dataExchangeModel twoWayMany2Many;
+twoWayMany2ManyProps
 {
     liggghtsPath "../DEM/in.liggghts_init";
 } 
 </PRE>
 <P><B>Description:</B>
 </P>
-<P>The data exchange model performs the data exchange between the DEM code and the CFD code. The twoWayM2M model is a model that can exchange particle properties from DEM to CFD and from CFD to DEM. Data is exchanged via MPI technique using the many to many mapping scheme. The DEM run is executed by the coupling model, via a liggghtsCommandModel object.
+<P>The data exchange model performs the data exchange between the DEM code and the CFD code. The twoWayMany2Many model is a model that can exchange particle properties from DEM to CFD and from CFD to DEM. Data is exchanged via MPI technique using the many to many mapping scheme. The DEM run is executed by the coupling model, via a liggghtsCommandModel object.
 </P>
 <P><B>Restrictions:</B>
 </P>
-<P>Should be used in combination with the turboEngineSearchM2M locate model to achieve best performance!
+<P>Must be used in combination with the engineSearchMany2Many locate model!
 </P>
 <P><B>Related commands:</B>
 </P>

doc/dataExchangeModel_twoWayMany2Many.txt

@@ -0,0 +1,42 @@
+"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c
+
+:link(lws,http://www.cfdem.com)
+:link(lc,CFDEMcoupling_Manual.html#comm)
+
+:line
+
+dataExchangeModel_twoWayMany2Many command :h3
+
+[Syntax:]
+
+Defined in couplingProperties dictionary.
+
+dataExchangeModel twoWayMany2Many;
+twoWayMany2ManyProps
+\{
+    liggghtsPath "path";
+\}; :pre
+
+{path} = path to the DEM simulation input file :ulb,l
+:ule
+
+[Examples:]
+
+dataExchangeModel twoWayMany2Many;
+twoWayMany2ManyProps
+\{
+    liggghtsPath "../DEM/in.liggghts_init";
+\} :pre
+
+[Description:]
+
+The data exchange model performs the data exchange between the DEM code and the CFD code. The twoWayMany2Many model is a model that can exchange particle properties from DEM to CFD and from CFD to DEM. Data is exchanged via MPI technique using the many to many mapping scheme. The DEM run is executed by the coupling model, via a liggghtsCommandModel object.
+
+[Restrictions:]
+
+Must be used in combination with the engineSearchMany2Many locate model!
+
+[Related commands:]
+
+"dataExchangeModel"_dataExchangeModel.html
+

doc/forceModel.html

@@ -33,7 +33,7 @@
 </P>
 <P><B>Description:</B>
 </P>
-<P>The force model performs the calculation of forces (e.g. fluid-particle interaction forces) acting on each DEM particle. All force models selected are executed sequentially and the forces on the particles are superposed.
+<P>The force model performs the calculation of forces (e.g. fluid-particle interaction forces) acting on each DEM particle. All force models selected are executed sequentially and the forces on the particles are superposed. If the fluid density field is needed, by default a field named "rho" will be used. Via the forceSubModel an alternative field can be chosen.
 </P>
 <P><B>Restrictions:</B>
 </P>

doc/forceModel.txt

@@ -31,7 +31,7 @@ Note: This examples list might not be complete - please look for other models (f
 
 [Description:]
 
-The force model performs the calculation of forces (e.g. fluid-particle interaction forces) acting on each DEM particle. All force models selected are executed sequentially and the forces on the particles are superposed.
+The force model performs the calculation of forces (e.g. fluid-particle interaction forces) acting on each DEM particle. All force models selected are executed sequentially and the forces on the particles are superposed. If the fluid density field is needed, by default a field named "rho" will be used. Via the forceSubModel an alternative field can be chosen.
 
 [Restrictions:]
 

doc/forceModel_Archimedes.html

@@ -19,13 +19,10 @@
 );
 ArchimedesProps
 {
-    densityFieldName "density";
     gravityFieldName "gravity";
 }; 
 </PRE>
-<UL><LI><I>density</I> = name of the finite volume density field 
-
-<LI><I>gravity</I> = name of the finite volume gravity field 
+<UL><LI><I>gravity</I> = name of the finite volume gravity field 
 
 
 </UL>
@@ -37,7 +34,6 @@ ArchimedesProps
 );
 ArchimedesProps
 {
-    densityFieldName "rho";
     gravityFieldName "g";
 } 
 </PRE>

doc/forceModel_Archimedes.txt

@@ -17,12 +17,10 @@ forceModels
 );
 ArchimedesProps
 \{
-    densityFieldName "density";
     gravityFieldName "gravity";
 \}; :pre
 
-{density} = name of the finite volume density field :ulb,l
-{gravity} = name of the finite volume gravity field :l
+{gravity} = name of the finite volume gravity field :ulb,l
 :ule
 
 [Examples:]
@@ -33,7 +31,6 @@ forceModels
 );
 ArchimedesProps
 \{
-    densityFieldName "rho";
     gravityFieldName "g";
 \} :pre
 

doc/forceModel_ArchimedesIB.html

@@ -19,14 +19,11 @@
 );
 ArchimedesIBProps
 {
-    densityFieldName "density";
     gravityFieldName "gravity";
     voidfractionFieldName "voidfraction";
 }; 
 </PRE>
-<UL><LI><I>density</I> = name of the finite volume density field 
-
-<LI><I>gravity</I> = name of the finite volume gravity field 
+<UL><LI><I>gravity</I> = name of the finite volume gravity field 
 
 <LI><I>voidfraction</I> = name of the finite volume voidfraction field 
 
@@ -40,7 +37,6 @@ ArchimedesIBProps
 );
 ArchimedesIBProps
 {
-    densityFieldName "rho";
     gravityFieldName "g";
     voidfractionFieldName "voidfractionNext";
 } 

doc/forceModel_ArchimedesIB.txt

@@ -17,13 +17,11 @@ forceModels
 );
 ArchimedesIBProps
 \{
-    densityFieldName "density";
     gravityFieldName "gravity";
     voidfractionFieldName "voidfraction";
 \}; :pre
 
-{density} = name of the finite volume density field :ulb,l
-{gravity} = name of the finite volume gravity field :l
+{gravity} = name of the finite volume gravity field :ulb,l
 {voidfraction} = name of the finite volume voidfraction field :l
 :ule
 
@@ -35,7 +33,6 @@ forceModels
 );
 ArchimedesIBProps
 \{
-    densityFieldName "rho";
     gravityFieldName "g";
     voidfractionFieldName "voidfractionNext";
 \} :pre

doc/forceModel_DiFeliceDrag.html

@@ -20,15 +20,12 @@
 DiFeliceDragProps
 {
     velFieldName "U";
-    densityFieldName "density";
-    interpolation;
+    interpolation switch1;
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
 
-<LI><I>density</I> = name of the finite volume gravity field 
-
-<LI><I>interpolation</I> = flag to use interpolated voidfraction and velocity values (normally off) 
+<LI><I>switch1</I> = flag to use interpolated voidfraction and velocity values (normally off) 
 
 
 </UL>
@@ -41,8 +38,7 @@ DiFeliceDragProps
 DiFeliceDragProps
 {
     velFieldName "U";
-    densityFieldName "rho";
-    interpolation;
+    interpolation true;
 } 
 </PRE>
 <P><B>Description:</B>

doc/forceModel_DiFeliceDrag.txt

@@ -18,13 +18,11 @@ forceModels
 DiFeliceDragProps
 \{
     velFieldName "U";
-    densityFieldName "density";
-    interpolation;
+    interpolation switch1;
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
-{density} = name of the finite volume gravity field :l
-{interpolation} = flag to use interpolated voidfraction and velocity values (normally off) :l
+{switch1} = flag to use interpolated voidfraction and velocity values (normally off) :l
 :ule
 
 [Examples:]
@@ -36,8 +34,7 @@ forceModels
 DiFeliceDragProps
 \{
     velFieldName "U";
-    densityFieldName "rho";
-    interpolation;
+    interpolation true;
 \} :pre
 
 [Description:]

doc/forceModel_GidaspowDrag.html

@@ -20,24 +20,24 @@
 GidaspowDragProps
 {
     velFieldName "U";
-    densityFieldName "density";
     voidfractionFieldName "voidfraction";
+    granVelFieldName "Us";
     phi "scalar";
-    interpolation;
-    implDEM;
+    interpolation switch1;
+    implForceDEM switch2;
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
 
-<LI><I>density</I> = name of the finite volume gravity field 
-
 <LI><I>voidfraction</I> = name of the finite volume voidfraction field 
 
+<LI><I>Us</I> = name of the finite volume cell averaged particle velocity field 
+
 <LI><I>phi</I> = drag correction factor (in doubt 1) 
 
-<LI><I>interpolation</I> = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values 
+<LI><I>switch1</I> = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values 
 
-<I>implDEM</I> = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
+<I>switch2</I> = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
 
 </UL>
 <P><B>Examples:</B>
@@ -49,8 +49,8 @@ GidaspowDragProps
 GidaspowDragProps
 {
     velFieldName "U";
-    densityFieldName "rho";
     voidfractionFieldName "voidfraction";
+    granVelFieldName "Us";
 } 
 </PRE>
 <P><B>Description:</B>

doc/forceModel_GidaspowDrag.txt

@@ -18,19 +18,19 @@ forceModels
 GidaspowDragProps
 \{
     velFieldName "U";
-    densityFieldName "density";
     voidfractionFieldName "voidfraction";
+    granVelFieldName "Us";
     phi "scalar";
-    interpolation;
-    implDEM;
+    interpolation switch1;
+    implForceDEM switch2;
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
-{density} = name of the finite volume gravity field :l
 {voidfraction} = name of the finite volume voidfraction field :l
+{Us} = name of the finite volume cell averaged particle velocity field :l
 {phi} = drag correction factor (in doubt 1) :l
-{interpolation} = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values :l
-{implDEM} = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
+{switch1} = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values :l
+{switch2} = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
 :ule
 
 [Examples:]
@@ -42,8 +42,8 @@ forceModels
 GidaspowDragProps
 \{
     velFieldName "U";
-    densityFieldName "rho";
     voidfractionFieldName "voidfraction";
+    granVelFieldName "Us";
 \} :pre
 
 [Description:]

doc/forceModel_KochHillDrag.html

@@ -20,21 +20,18 @@
 KochHillDragProps
 {
     velFieldName "U";
-    densityFieldName "density";
     voidfractionFieldName "voidfraction"; 
-    interpolation;
-    implDEM;
+    interpolation "bool1";
+    implForceDEM "bool2";
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
 
-<LI><I>density</I> = name of the finite volume gravity field 
-
 <LI><I>voidfraction</I> = name of the finite volume voidfraction field 
 
-<LI><I>interpolation</I> = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values 
+<LI><I>bool1</I> = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values 
 
-<I>implDEM</I> = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
+<I>bool2</I> = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
 
 </UL>
 <P><B>Examples:</B>
@@ -46,7 +43,6 @@ KochHillDragProps
 KochHillDragProps
 {
     velFieldName "U";
-    densityFieldName "rho";
     voidfractionFieldName "voidfraction"; 
 } 
 </PRE>

doc/forceModel_KochHillDrag.txt

@@ -18,17 +18,15 @@ forceModels
 KochHillDragProps
 \{
     velFieldName "U";
-    densityFieldName "density";
     voidfractionFieldName "voidfraction"; 
-    interpolation;
-    implDEM;
+    interpolation "bool1";
+    implForceDEM "bool2";
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
-{density} = name of the finite volume gravity field :l
 {voidfraction} = name of the finite volume voidfraction field :l
-{interpolation} = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values :l
-{implDEM} = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
+{bool1} = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values :l
+{bool2} = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
 :ule
 
 [Examples:]
@@ -40,7 +38,6 @@ forceModels
 KochHillDragProps
 \{
     velFieldName "U";
-    densityFieldName "rho";
     voidfractionFieldName "voidfraction"; 
 \} :pre
 

doc/forceModel_LaEuScalarTemp.html

@@ -21,23 +21,19 @@ LaEuScalarTempProps
 {
     velFieldName "U";
     tempFieldName "T";
-    tempSourceFieldName "Tsource";
     voidfractionFieldName "voidfraction";
     partTempName "Temp";
     partHeatFluxName "convectiveHeatFlux";
     lambda value;
     Cp value1;
-    densityFieldName "density";
-    interpolation;
-    verbose;
+    interpolation "switch1";
+    verbose "switch2";
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
 
 <LI><I>T</I> = name of the finite volume scalar temperature field 
 
-<LI><I>Tsource</I> = name of the finite volume scalar temperature source field 
-
 <LI><I>voidfraction</I> = name of the finite volume voidfraction field 
 
 <LI><I>Temp</I> = name of the DEM data representing the particles temperature 
@@ -48,11 +44,9 @@ LaEuScalarTempProps
 
 <LI><I>value1</I> = fluid specific heat capacity [W*s/(kg*K)] 
 
-<LI><I>density</I> = name of the finite volume fluid density field 
+<LI><I>switch1</I> = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values 
 
-<LI><I>interpolation</I> = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values 
-
-<LI><I>verbose</I> = (normally off) for verbose run 
+<LI><I>switch2</I> = (normally off) for verbose run 
 
 
 </UL>
@@ -66,13 +60,11 @@ LaEuScalarTempProps
 {
     velFieldName "U";
     tempFieldName "T";
-    tempSourceFieldName "Tsource";
     voidfractionFieldName "voidfraction";
     partTempName "Temp";
     partHeatFluxName "convectiveHeatFlux";
     lambda 0.0256;
     Cp 1007;
-    densityFieldName "rho";
 } 
 </PRE>
 <P><B>Description:</B>
@@ -81,7 +73,7 @@ LaEuScalarTempProps
 </P>
 <P><B>Restrictions:</B>
 </P>
-<P>Goes only with cfdemSolverScalar.
+<P>Goes only with cfdemSolverScalar. The force model has to be the second (!!!) model in the forces list.
 </P>
 <P><B>Related commands:</B>
 </P>

doc/forceModel_LaEuScalarTemp.txt

@@ -19,28 +19,24 @@ LaEuScalarTempProps
 \{
     velFieldName "U";
     tempFieldName "T";
-    tempSourceFieldName "Tsource";
     voidfractionFieldName "voidfraction";
     partTempName "Temp";
     partHeatFluxName "convectiveHeatFlux";
     lambda value;
     Cp value1;
-    densityFieldName "density";
-    interpolation;
-    verbose;
+    interpolation "switch1";
+    verbose "switch2";
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
 {T} = name of the finite volume scalar temperature field :l
-{Tsource} = name of the finite volume scalar temperature source field :l
 {voidfraction} = name of the finite volume voidfraction field :l
 {Temp} = name of the DEM data representing the particles temperature :l
 {convectiveHeatFlux} = name of the DEM data representing the particle-fluid convective heat flux :l
 {value} = fluid thermal conductivity \[W/(m*K)\] :l
 {value1} = fluid specific heat capacity \[W*s/(kg*K)\] :l
-{density} = name of the finite volume fluid density field :l
-{interpolation} = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values :l
-{verbose} = (normally off) for verbose run :l
+{switch1} = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values :l
+{switch2} = (normally off) for verbose run :l
 :ule
 
 [Examples:]
@@ -53,13 +49,11 @@ LaEuScalarTempProps
 \{
     velFieldName "U";
     tempFieldName "T";
-    tempSourceFieldName "Tsource";
     voidfractionFieldName "voidfraction";
     partTempName "Temp";
     partHeatFluxName "convectiveHeatFlux";
     lambda 0.0256;
     Cp 1007;
-    densityFieldName "rho";
 \} :pre
 
 [Description:]
@@ -68,7 +62,7 @@ This "forceModel" does not influence the particles or the fluid flow! Using the
 
 [Restrictions:]
 
-Goes only with cfdemSolverScalar.
+Goes only with cfdemSolverScalar. The force model has to be the second (!!!) model in the forces list.
 
 [Related commands:]
 

doc/forceModel_MeiLift.html

@@ -20,21 +20,18 @@
 MeiLiftProps
 {
     velFieldName "U";
-    densityFieldName "density";
     useSecondOrderTerms;
-    interpolation;
-    verbose;
+    interpolation "switch1";
+    verbose "switch2";
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
 
-<LI><I>density</I> = name of the finite volume fluid density field 
-
 <LI><I>useSecondOrderTerms</I> = switch to activate second order terms in the lift force model 
 
-<LI><I>interpolation</I> = switch to activate tri-linear interpolation of the flow quantities at the particle position 
+<LI><I>switch1</I> = switch to activate tri-linear interpolation of the flow quantities at the particle position 
 
-<LI><I>verbose</I> = switch to activate the report of per-particle quantities to the screen 
+<LI><I>switch2</I> = switch to activate the report of per-particle quantities to the screen 
 
 
 </UL>
@@ -47,10 +44,9 @@ MeiLiftProps
 MeiLiftProps
 {
     velFieldName "U";
-    densityFieldName "rho";
     useSecondOrderTerms;
-    interpolation;
-    verbose;
+    interpolation true;
+    verbose true;
 } 
 </PRE>
 <P><B>Description:</B>

doc/forceModel_MeiLift.txt

@@ -18,17 +18,15 @@ forceModels
 MeiLiftProps
 \{
     velFieldName "U";
-    densityFieldName "density";
     useSecondOrderTerms;
-    interpolation;
-    verbose;
+    interpolation "switch1";
+    verbose "switch2";
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
-{density} = name of the finite volume fluid density field :l
 {useSecondOrderTerms} = switch to activate second order terms in the lift force model :l
-{interpolation} = switch to activate tri-linear interpolation of the flow quantities at the particle position :l
-{verbose} = switch to activate the report of per-particle quantities to the screen :l
+{switch1} = switch to activate tri-linear interpolation of the flow quantities at the particle position :l
+{switch2} = switch to activate the report of per-particle quantities to the screen :l
 :ule
 
 [Examples:]
@@ -40,10 +38,9 @@ forceModels
 MeiLiftProps
 \{
     velFieldName "U";
-    densityFieldName "rho";
     useSecondOrderTerms;
-    interpolation;
-    verbose;
+    interpolation true;
+    verbose true;
 \} :pre
 
 [Description:]

doc/forceModel_SchillerNaumannDrag.html

@@ -20,13 +20,10 @@
 SchillerNaumannDragProps
 {
     velFieldName "U";
-    densityFieldName "density";
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
 
-<LI><I>density</I> = name of the finite volume gravity field 
-
 
 </UL>
 <P><B>Examples:</B>
@@ -38,7 +35,6 @@ SchillerNaumannDragProps
 SchillerNaumannDragProps
 {
     velFieldName "U";
-    densityFieldName "rho";
 } 
 </PRE>
 <P><B>Description:</B>

doc/forceModel_SchillerNaumannDrag.txt

@@ -18,11 +18,9 @@ forceModels
 SchillerNaumannDragProps
 \{
     velFieldName "U";
-    densityFieldName "density";
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
-{density} = name of the finite volume gravity field :l
 :ule
 
 [Examples:]
@@ -34,7 +32,6 @@ forceModels
 SchillerNaumannDragProps
 \{
     velFieldName "U";
-    densityFieldName "rho";
 \} :pre
 
 [Description:]

doc/forceModel_ShirgaonkarIB.html

@@ -20,14 +20,11 @@
 ShirgaonkarIBProps
 {
     velFieldName "U";
-    densityFieldName "density";
     pressureFieldName "pressure";
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
 
-<LI><I>density</I> = name of the finite volume density field 
-
 <LI><I>pressure</I> = name of the finite volume pressure field 
 
 
@@ -41,7 +38,6 @@ ShirgaonkarIBProps
 ShirgaonkarIBProps
 {
     velFieldName "U";
-    densityFieldName "rho";
     pressureFieldName "p";
 } 
 </PRE>

doc/forceModel_ShirgaonkarIB.txt

@@ -18,12 +18,10 @@ forceModels
 ShirgaonkarIBProps
 \{
     velFieldName "U";
-    densityFieldName "density";
     pressureFieldName "pressure";
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
-{density} = name of the finite volume density field :l
 {pressure} = name of the finite volume pressure field :l
 :ule
 
@@ -36,7 +34,6 @@ forceModels
 ShirgaonkarIBProps
 \{
     velFieldName "U";
-    densityFieldName "rho";
     pressureFieldName "p";
 \} :pre
 

doc/forceModel_fieldStore.html

@@ -0,0 +1,68 @@
+<HTML>
+<CENTER><A HREF = "http://www.cfdem.com">CFDEMproject WWW Site</A> - <A HREF = "CFDEMcoupling_Manual.html#comm">CFDEM Commands</A> 
+</CENTER>
+
+
+
+
+<HR>
+
+<H3>forceModel_fieldStore command 
+</H3>
+<P><B>Syntax:</B>
+</P>
+<P>Defined in couplingProperties dictionary.
+</P>
+<PRE>forceModels
+(
+    fieldStore
+);
+fieldStoreProps
+{
+    scalarFieldNames
+    (
+        "scalarField"
+    );
+    vectorFieldNames
+    (
+        "vectorField"
+    );
+}; 
+</PRE>
+<UL><LI><I>scalarField</I> = names of the finite volume scalar fields to be stored 
+
+<LI><I>vectorField</I> = names of the finite volume vector fields to be stored 
+
+
+</UL>
+<P><B>Examples:</B>
+</P>
+<PRE>forceModels
+(
+    fieldStore
+);
+fieldStoreProps
+{
+    scalarFieldNames
+    (
+        "voidfraction"
+    );
+    vectorFieldNames
+    (
+        "U"
+    );
+} 
+</PRE>
+<P><B>Description:</B>
+</P>
+<P>This "forceModel" does not influence the particles or the flow - it is a tool to store a scalar/vector field! This is especially useful if you use a boundary condition which cannot interpreted correctly in your postporcessor (e.g. paraview).
+</P>
+<P><B>Restrictions:</B>
+</P>
+<P>none.
+</P>
+<P><B>Related commands:</B>
+</P>
+<P><A HREF = "forceModel.html">forceModel</A>
+</P>
+</HTML>

doc/forceModel_fieldStore.txt

@@ -0,0 +1,63 @@
+"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c
+
+:link(lws,http://www.cfdem.com)
+:link(lc,CFDEMcoupling_Manual.html#comm)
+
+:line
+
+forceModel_fieldStore command :h3
+
+[Syntax:]
+
+Defined in couplingProperties dictionary.
+
+forceModels
+(
+    fieldStore
+);
+fieldStoreProps
+\{
+    scalarFieldNames
+    (
+        "scalarField"
+    );
+    vectorFieldNames
+    (
+        "vectorField"
+    );
+\}; :pre
+
+{scalarField} = names of the finite volume scalar fields to be stored :ulb,l
+{vectorField} = names of the finite volume vector fields to be stored :l
+:ule
+
+[Examples:]
+
+forceModels
+(
+    fieldStore
+);
+fieldStoreProps
+\{
+    scalarFieldNames
+    (
+        "voidfraction"
+    );
+    vectorFieldNames
+    (
+        "U"
+    );
+\} :pre
+
+[Description:]
+
+This "forceModel" does not influence the particles or the flow - it is a tool to store a scalar/vector field! This is especially useful if you use a boundary condition which cannot interpreted correctly in your postporcessor (e.g. paraview).
+
+[Restrictions:]
+
+none.
+
+[Related commands:]
+
+"forceModel"_forceModel.html
+

doc/forceModel_gradPForce.html

@@ -20,18 +20,15 @@
 gradPForceProps
 {
     pFieldName "pressure";
-    densityFieldName "density";
     velocityFieldName "U";
-    interpolation;
+    interpolation switch1;
 }; 
 </PRE>
 <UL><LI><I>pressure</I> = name of the finite volume fluid pressure field 
 
-<LI><I>density</I> = name of the finite volume gravity field 
-
 <LI><I>U</I> = name of the finite volume fluid velocity field 
 
-<LI><I>interpolation</I> = flag to use interpolated pressure values (normally off) 
+<LI><I>switch1</I> = flag to use interpolated pressure values (normally off) 
 
 
 </UL>
@@ -44,9 +41,8 @@ gradPForceProps
 gradPForceProps
 {
     pFieldName "p";
-    densityFieldName "rho";
     velocityFieldName "U";
-    interpolation;
+    interpolation true;
 } 
 </PRE>
 <P><B>Description:</B>

doc/forceModel_gradPForce.txt

@@ -18,15 +18,13 @@ forceModels
 gradPForceProps
 \{
     pFieldName "pressure";
-    densityFieldName "density";
     velocityFieldName "U";
-    interpolation;
+    interpolation switch1;
 \}; :pre
 
 {pressure} = name of the finite volume fluid pressure field :ulb,l
-{density} = name of the finite volume gravity field :l
 {U} = name of the finite volume fluid velocity field :l
-{interpolation} = flag to use interpolated pressure values (normally off) :l
+{switch1} = flag to use interpolated pressure values (normally off) :l
 :ule
 
 [Examples:]
@@ -38,9 +36,8 @@ forceModels
 gradPForceProps
 \{
     pFieldName "p";
-    densityFieldName "rho";
     velocityFieldName "U";
-    interpolation;
+    interpolation true;
 \} :pre
 
 [Description:]

doc/forceModel_virtualMassForce.html

@@ -20,13 +20,10 @@
 virtualMassForceProps
 {
     velFieldName "U";
-    densityFieldName "density";
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
 
-<LI><I>density</I> = name of the finite volume fluid density field 
-
 
 </UL>
 <P><B>Examples:</B>
@@ -38,7 +35,6 @@ virtualMassForceProps
 virtualMassForceProps
 {
     velFieldName "U";
-    densityFieldName "rho";
 } 
 </PRE>
 <P><B>Description:</B>

doc/forceModel_virtualMassForce.txt

@@ -18,11 +18,9 @@ forceModels
 virtualMassForceProps
 \{
     velFieldName "U";
-    densityFieldName "density";
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
-{density} = name of the finite volume fluid density field :l
 :ule
 
 [Examples:]
@@ -34,7 +32,6 @@ forceModels
 virtualMassForceProps
 \{
     velFieldName "U";
-    densityFieldName "rho";
 \} :pre
 
 [Description:]

doc/forceModel_viscForce.html

@@ -20,15 +20,12 @@
 viscForceProps
 {
     velocityFieldName "U";
-    densityFieldName "density";
-    interpolation;
+    interpolation "switch";
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
 
-<LI><I>density</I> = name of the finite volume gravity field 
-
-<LI><I>interpolation</I> = flag to use interpolated stress values (normally off) 
+<LI><I>switch</I> = flag to use interpolated stress values (normally off) 
 
 
 </UL>
@@ -41,7 +38,6 @@ viscForceProps
 viscForceProps
 {
     velocityFieldName "U";
-    densityFieldName "density";
 } 
 </PRE>
 <P><B>Description:</B>

doc/forceModel_viscForce.txt

@@ -18,13 +18,11 @@ forceModels
 viscForceProps
 \{
     velocityFieldName "U";
-    densityFieldName "density";
-    interpolation;
+    interpolation "switch";
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
-{density} = name of the finite volume gravity field :l
-{interpolation} = flag to use interpolated stress values (normally off) :l
+{switch} = flag to use interpolated stress values (normally off) :l
 :ule
 
 [Examples:]
@@ -36,7 +34,6 @@ forceModels
 viscForceProps
 \{
     velocityFieldName "U";
-    densityFieldName "density";
 \} :pre
 
 [Description:]

doc/forceSubModel.html

@@ -0,0 +1,49 @@
+<HTML>
+<CENTER><A HREF = "http://www.cfdem.com">CFDEMproject WWW Site</A> - <A HREF = "CFDEMcoupling_Manual.html#comm">CFDEM Commands</A> 
+</CENTER>
+
+
+
+
+<HR>
+
+<H3>forceSubModel command 
+</H3>
+<P><B>Syntax:</B>
+</P>
+<P>Defined in couplingProperties sub-dictionary of the force model in use. If no force sub-model is applied ImEx is used as default. If the keyword "forceSubModels" is provided, a choice of sub model is demanded.
+</P>
+<PRE>forceSubModels
+(
+    model_x
+    model_y
+); 
+</PRE>
+<UL><LI>model = name of force sub-model to be applied 
+</UL>
+<P><B>Examples:</B>
+</P>
+<PRE>forceSubModels
+(
+    ImEx
+); 
+</PRE>
+<P>Note: This examples list might not be complete - please look for other models (forceSubModel_XY) in this documentation.
+</P>
+<P><B>Description:</B>
+</P>
+<P>The force sub model is designed to hold the settings a force model can have. For now it handles the treatExplicit, treatDEM and implDEM option. 
+</P>
+<P><B>Restrictions:</B>
+</P>
+<P>None.
+</P>
+<P><B>Related commands:</B>
+</P>
+<P><A HREF = "forceSubModel_ImEx.html">ImEx</A>
+</P>
+<P>Note: This examples list may be incomplete - please look for other models (forceSubModel_XY) in this documentation.
+</P>
+<P><B>Default:</B> none.
+</P>
+</HTML>

doc/forceSubModel.txt

@@ -0,0 +1,45 @@
+"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c
+
+:link(lws,http://www.cfdem.com)
+:link(lc,CFDEMcoupling_Manual.html#comm)
+
+:line
+
+forceSubModel command :h3
+
+[Syntax:]
+
+Defined in couplingProperties sub-dictionary of the force model in use. If no force sub-model is applied ImEx is used as default. If the keyword "forceSubModels" is provided, a choice of sub model is demanded.
+
+forceSubModels
+(
+    model_x
+    model_y
+); :pre
+
+model = name of force sub-model to be applied :ul
+
+[Examples:]
+
+forceSubModels
+(
+    ImEx
+); :pre
+
+Note: This examples list might not be complete - please look for other models (forceSubModel_XY) in this documentation.
+
+[Description:]
+
+The force sub model is designed to hold the settings a force model can have. For now it handles the treatExplicit, treatDEM and implDEM option. 
+
+[Restrictions:]
+
+None.
+
+[Related commands:]
+
+"ImEx"_forceSubModel_ImEx.html
+
+Note: This examples list may be incomplete - please look for other models (forceSubModel_XY) in this documentation.
+
+[Default:] none.

doc/forceSubModel_ImEx.html

@@ -0,0 +1,45 @@
+<HTML>
+<CENTER><A HREF = "http://www.cfdem.com">CFDEMproject WWW Site</A> - <A HREF = "CFDEMcoupling_Manual.html#comm">CFDEM Commands</A> 
+</CENTER>
+
+
+
+
+<HR>
+
+<H3>forceSubModel_ImEx command 
+</H3>
+<P><B>Syntax:</B>
+</P>
+<P>Defined in couplingProperties sub-dictionary of the force model in use.
+</P>
+<P>forceSubModels
+(
+    ImEx;
+);
+</P>
+<P>treatExplicit true;  // optional for some force models.
+treatDEM true;       // optional for some force models.
+implDEM true;        // optional for some force models.
+</P>
+<P><B>Examples:</B>
+</P>
+<P>forceSubModels
+(
+    ImEx;
+);
+treatExplicit true;  // optional for some force models.
+</P>
+<P><B>Description:</B>
+</P>
+<P> If no force sub-model is applied ImEx is used as default. If the keyword "forceSubModels" is provided, a choice of sub model is demanded. Depending on the force model different keywords are read and can therefrore be set (see the log file). If the keyword is provided, its value is used.
+</P>
+<P><B>Restrictions:</B>
+</P>
+<P>none.
+</P>
+<P><B>Related commands:</B>
+</P>
+<P><A HREF = "forceSubModel.html">forceSubModel</A>
+</P>
+</HTML>

doc/forceSubModel_ImEx.txt

@@ -0,0 +1,42 @@
+"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c
+
+:link(lws,http://www.cfdem.com)
+:link(lc,CFDEMcoupling_Manual.html#comm)
+
+:line
+
+forceSubModel_ImEx command :h3
+
+[Syntax:]
+
+Defined in couplingProperties sub-dictionary of the force model in use.
+
+forceSubModels
+(
+    ImEx;
+);
+
+treatExplicit true;  // optional for some force models.
+treatDEM true;       // optional for some force models.
+implDEM true;        // optional for some force models.
+
+[Examples:]
+
+forceSubModels
+(
+    ImEx;
+);
+treatExplicit true;  // optional for some force models.
+
+[Description:]
+
+ If no force sub-model is applied ImEx is used as default. If the keyword "forceSubModels" is provided, a choice of sub model is demanded. Depending on the force model different keywords are read and can therefrore be set (see the log file). If the keyword is provided, its value is used.
+
+[Restrictions:]
+
+none.
+
+[Related commands:]
+
+"forceSubModel"_forceSubModel.html
+

doc/forceSubModel_ImExCorr.html

@@ -0,0 +1,46 @@
+<HTML>
+<CENTER><A HREF = "http://www.cfdem.com">CFDEMproject WWW Site</A> - <A HREF = "CFDEMcoupling_Manual.html#comm">CFDEM Commands</A> 
+</CENTER>
+
+
+
+
+<HR>
+
+<H3>forceSubModel_ImExCorr command 
+</H3>
+<P><B>Syntax:</B>
+</P>
+<P>Defined in couplingProperties sub-dictionary of the force model in use.
+</P>
+<P>forceSubModels
+(
+    ImExCorr;
+);
+</P>
+<P>treatExplicit true;  // optional for some force models.
+treatDEM true;       // optional for some force models.
+implDEM true;        // optional for some force models.
+explicitInterpCorr true; // optional for some force models.
+</P>
+<P><B>Examples:</B>
+</P>
+<P>forceSubModels
+(
+    ImExCorr;
+);
+treatExplicit true;  // optional for some force models.
+</P>
+<P><B>Description:</B>
+</P>
+<P> Same as ImEx, but it additionally reads "explicitInterpCorr" to correct the error steming from interpolation of Ufluid and averaging of Uparticles.
+</P>
+<P><B>Restrictions:</B>
+</P>
+<P>none.
+</P>
+<P><B>Related commands:</B>
+</P>
+<P><A HREF = "forceSubModel.html">forceSubModel</A>
+</P>
+</HTML>

doc/forceSubModel_ImExCorr.txt

@@ -0,0 +1,43 @@
+"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c
+
+:link(lws,http://www.cfdem.com)
+:link(lc,CFDEMcoupling_Manual.html#comm)
+
+:line
+
+forceSubModel_ImExCorr command :h3
+
+[Syntax:]
+
+Defined in couplingProperties sub-dictionary of the force model in use.
+
+forceSubModels
+(
+    ImExCorr;
+);
+
+treatExplicit true;  // optional for some force models.
+treatDEM true;       // optional for some force models.
+implDEM true;        // optional for some force models.
+explicitInterpCorr true; // optional for some force models.
+
+[Examples:]
+
+forceSubModels
+(
+    ImExCorr;
+);
+treatExplicit true;  // optional for some force models.
+
+[Description:]
+
+ Same as ImEx, but it additionally reads "explicitInterpCorr" to correct the error steming from interpolation of Ufluid and averaging of Uparticles.
+
+[Restrictions:]
+
+none.
+
+[Related commands:]
+
+"forceSubModel"_forceSubModel.html
+

doc/liggghtsCommandModel_writeLiggghts.html

@@ -26,13 +26,12 @@ writeLiggghtsProps
     verbose;
 } 
 </PRE>
-<UL><LI><I>switch1</I> = switch (choose on/off) to select if only last step is stored or every write step. 
+<UL><LI><I>switch1</I> = switch (choose on/off) to select if only last step is stored or every write step (default on). 
 
-<LI><I>name</I> = name of the restart file to be written in /$caseDir/DEM/  default default "liggghts.restartCFDEM" 
+<LI><I>name</I> = name of the restart file to be written in /$caseDir/DEM/  default (default "liggghts.restartCFDEM") 
 
-<LI><I>switch2</I> = switch (choose on/off) to select if only one restart file $name or many files $name_$timeStamp are written 
-
-<LI><I>verbose</I> = (normally off) for verbose run 
+<LI><I>switch2</I> = switch (choose on/off) to select if only one restart file $name or many files $name_$timeStamp are written (default off):l
+<I>verbose</I> = (default off) for verbose run 
 
 
 </UL>

doc/liggghtsCommandModel_writeLiggghts.txt

@@ -24,10 +24,10 @@ writeLiggghtsProps
     verbose;
 \} :pre
 
-{switch1} = switch (choose on/off) to select if only last step is stored or every write step. :ulb,l
-{name} = name of the restart file to be written in /$caseDir/DEM/  default default "liggghts.restartCFDEM" :l
-{switch2} = switch (choose on/off) to select if only one restart file $name or many files $name_$timeStamp are written :l
-{verbose} = (normally off) for verbose run :l
+{switch1} = switch (choose on/off) to select if only last step is stored or every write step (default on). :ulb,l
+{name} = name of the restart file to be written in /$caseDir/DEM/  default (default "liggghts.restartCFDEM") :l
+{switch2} = switch (choose on/off) to select if only one restart file $name or many files $name_$timeStamp are written (default off):l
+{verbose} = (default off) for verbose run :l
 :ule
 
 [Examples:]

doc/locateModel_turboEngineSearchM2M.html

@@ -1,58 +0,0 @@
-<HTML>
-<CENTER><A HREF = "http://www.cfdem.com">CFDEMproject WWW Site</A> - <A HREF = "CFDEMcoupling_Manual.html#comm">CFDEM Commands</A> 
-</CENTER>
-
-
-
-
-<HR>
-
-<H3>locateModel_turboEngineM2MSearch command 
-</H3>
-<P><B>Syntax:</B>
-</P>
-<P>Defined in couplingProperties dictionary.
-</P>
-<PRE>locateModel turboEngineM2M;
-turboEngineM2MProps
-{
-    turboEngineProps
-    {
-        treeSearch switch1;
-    }
-} 
-</PRE>
-<UL><LI><I>switch1</I> = switch to use tree search algorithm 
-
-
-</UL>
-<P><B>Examples:</B>
-</P>
-<PRE>locateModel turboEngineM2M;
-turboEngineM2MProps
-{
-    turboEngineProps
-    {
-        treeSearch true;
-    }
-} 
-</PRE>
-<P><B>Description:</B>
-</P>
-<P>The locateModel "turboEngineM2M" locates the CFD cell and cellID corresponding to a given position. The algorithm is improved compared to engine search to show better parallel performance.
-</P>
-<P>The turboEngineM2MSearch locate Model can be used with different settings to use different algorithms:
-</P>
-<UL><LI>faceDecomp false; treeSearch false;     will execute some geometric (linear) search using the last known cellID (recommended) 
-
-<LI>faceDecomp false; treeSearch true;      will use a recursive tree structure to find the cell. 
-
-
-</UL>
-<P><B>Restrictions:</B> This model can only be used with many to many data exchange model!
-</P>
-<P><B>Related commands:</B>
-</P>
-<P><A HREF = "locateModel.html">locateModel</A>
-</P>
-</HTML>

doc/locateModel_turboEngineSearchM2M.txt

@@ -1,52 +0,0 @@
-"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c
-
-:link(lws,http://www.cfdem.com)
-:link(lc,CFDEMcoupling_Manual.html#comm)
-
-:line
-
-locateModel_turboEngineM2MSearch command :h3
-
-[Syntax:]
-
-Defined in couplingProperties dictionary.
-
-locateModel turboEngineM2M;
-turboEngineM2MProps
-\{
-    turboEngineProps
-    \{
-        treeSearch switch1;
-    \}
-\} :pre
-
-{switch1} = switch to use tree search algorithm :ulb,l
-:ule
-
-[Examples:]
-
-locateModel turboEngineM2M;
-turboEngineM2MProps
-\{
-    turboEngineProps
-    \{
-        treeSearch true;
-    \}
-\} :pre
-
-[Description:]
-
-The locateModel "turboEngineM2M" locates the CFD cell and cellID corresponding to a given position. The algorithm is improved compared to engine search to show better parallel performance.
-
-The turboEngineM2MSearch locate Model can be used with different settings to use different algorithms:
-
-faceDecomp false; treeSearch false;     will execute some geometric (linear) search using the last known cellID (recommended) :ulb,l
-faceDecomp false; treeSearch true;      will use a recursive tree structure to find the cell. :l
-:ule
-
-[Restrictions:] This model can only be used with many to many data exchange model!
-
-[Related commands:]
-
-"locateModel"_locateModel.html
-

doc/momCoupleModel.html

@@ -33,7 +33,8 @@
 </P>
 <P>Note that the variable "imExSplitFactor" can be set in the couplingProperties in order to treat implicitly defined forces (in the implementation of the force model) as explicit ones. "imExSplitFactor 1.0;" is set by default, meaning that all implicit forces will be considered implicitly, whereas "imExSplitFactor 0.0;" would mean that implicitly defined forces will be treated in an explicit fashion.
 </P>
-<P><B>Description:</B>
+<P>Note that the switch "treatVoidCellsAsExplicitForce true;" can be set in the couplingProperties in order to change the treatment of cells which are void of particles. This is only relevant if (i) smoothing is used, and (ii) implicit force coupling is performed. By default, the particle veloctiy field (Us) will be smoothed to obtain a meaningful reference quantity for the implicit force coupling. In case "treatVoidCellsAsExplicitForce true;" is set, however, Us will not be smoothed and implicit forces (after the smoothing has been performed) in cells void of particles be treated as explicit ones. This avoids the problem of defining Us in cells that are void of particles, but for which an implicit coupling force is obtained in the smoothing process.
+<B>Description:</B>
 </P>
 <P>The momCoupleModel is the base class for momentum exchange between DEM and CFD simulation.
 </P>

doc/momCoupleModel.txt

@@ -31,6 +31,7 @@ Forces can be coupled in an implicit way to the fluid solver (i.e., when solving
 
 Note that the variable "imExSplitFactor" can be set in the couplingProperties in order to treat implicitly defined forces (in the implementation of the force model) as explicit ones. "imExSplitFactor 1.0;" is set by default, meaning that all implicit forces will be considered implicitly, whereas "imExSplitFactor 0.0;" would mean that implicitly defined forces will be treated in an explicit fashion.
 
+Note that the switch "treatVoidCellsAsExplicitForce true;" can be set in the couplingProperties in order to change the treatment of cells which are void of particles. This is only relevant if (i) smoothing is used, and (ii) implicit force coupling is performed. By default, the particle veloctiy field (Us) will be smoothed to obtain a meaningful reference quantity for the implicit force coupling. In case "treatVoidCellsAsExplicitForce true;" is set, however, Us will not be smoothed and implicit forces (after the smoothing has been performed) in cells void of particles be treated as explicit ones. This avoids the problem of defining Us in cells that are void of particles, but for which an implicit coupling force is obtained in the smoothing process.
 [Description:]
 
 The momCoupleModel is the base class for momentum exchange between DEM and CFD simulation.

src/.gitignore

@@ -0,0 +1,7 @@
+*.o
+*.d
+*.a
+*.dep
+log_*
+log.*
+*~

src/finiteVolume/Make/files

@@ -1,5 +1,8 @@
 fvPatchFields = fields/fvPatchFields
 derivedFvPatchFields = $(fvPatchFields)/derived
+
+general = cfdTools/general
+
 $(derivedFvPatchFields)/uniformFixedValueVoidfraction/uniformFixedValueVoidfractionFvPatchFields.C
 $(derivedFvPatchFields)/uniformFixedValueTube/uniformFixedValueTubeFvPatchFields.C
 

src/lagrangian/cfdemParticle/Make/files

@@ -1,11 +1,14 @@
 cfdemCloud = cfdemCloud
+cfdTools = cfdTools
 forceModels = subModels/forceModel
+forceSubModels = subModels/forceModel/forceSubModels
 forceModelsMS = subModels/forceModelMS
 IOModels = subModels/IOModel
 voidFractionModels = subModels/voidFractionModel
 locateModels = subModels/locateModel
 meshMotionModels = subModels/meshMotionModel
 momCoupleModels = subModels/momCoupleModel
+regionModels = subModels/regionModel
 dataExchangeModels = subModels/dataExchangeModel
 averagingModels = subModels/averagingModel
 clockModels = subModels/clockModel
@@ -14,68 +17,52 @@ smoothingModels = subModels/smoothingModel
 probeModels = subModels/probeModel
 
 $(cfdemCloud)/cfdemCloud.C
-derived/cfdemCloudBiDisperse/cfdemCloudBiDisperse.C
 derived/cfdemCloudIB/cfdemCloudIB.C
 derived/cfdemCloudMS/cfdemCloudMS.C
 
+$(cfdTools)/global.C
+$(cfdTools)/newGlobal.C
+
 $(forceModels)/forceModel/forceModel.C
 $(forceModels)/forceModel/newForceModel.C
 $(forceModels)/noDrag/noDrag.C
 $(forceModels)/checkCouplingInterval/checkCouplingInterval.C
 $(forceModels)/DiFeliceDrag/DiFeliceDrag.C
-$(forceModels)/DiFeliceDragNLift/DiFeliceDragNLift.C
+$(forceModels)/fieldStore/fieldStore.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)/interfaceParticleProbe/interfaceParticleProbe.C
-$(forceModels)/fieldStore/fieldStore.C
-$(forceModels)/fieldTimeAverage/fieldTimeAverage.C
-$(forceModels)/fieldBound/fieldBound.C
-$(forceModels)/volWeightedAverage/volWeightedAverage.C
-$(forceModels)/totalMomentumExchange/totalMomentumExchange.C
 $(forceModels)/KochHillDrag/KochHillDrag.C
 $(forceModels)/KochHillRWDrag/KochHillRWDrag.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)/HollowayDrag/HollowayDrag.C
 $(forceModels)/viscForce/viscForce.C
 $(forceModels)/MeiLift/MeiLift.C
-$(forceModels)/melting/melting.C
-$(forceModels)/KochHillDragNLift/KochHillDragNLift.C
-$(forceModels)/stokesSpheroidDrag/stokesSpheroidDrag.C
-$(forceModels)/solidsPressureForce/solidsPressureForce.C
-$(forceModels)/periodicPressure/periodicPressure.C
-$(forceModels)/periodicPressureControl/periodicPressureControl.C
-$(forceModels)/averageSlipVel/averageSlipVel.C
 $(forceModels)/particleCellVolume/particleCellVolume.C
+$(forceModels)/fieldTimeAverage/fieldTimeAverage.C
+$(forceModels)/volWeightedAverage/volWeightedAverage.C
 
 $(forceModelsMS)/forceModelMS/forceModelMS.C
 $(forceModelsMS)/forceModelMS/newForceModelMS.C
 $(forceModelsMS)/DiFeliceDragMS/DiFeliceDragMS.C
-$(forceModelsMS)/GidaspowDragMS/GidaspowDragMS.C
-$(forceModelsMS)/noDragMS/noDragMS.C
+
+$(forceSubModels)/forceSubModel/newForceSubModel.C
+$(forceSubModels)/forceSubModel/forceSubModel.C
+$(forceSubModels)/ImEx/ImEx.C
 
 $(probeModels)/probeModel/probeModel.C
 $(probeModels)/probeModel/newProbeModel.C
 $(probeModels)/noProbe/noProbe.C
 $(probeModels)/particleProbe/particleProbe.C
-/*$(probeModels)/interfaceParticleProbe/interfaceParticleProbe.C*/
 
 $(IOModels)/IOModel/IOModel.C
 $(IOModels)/IOModel/newIOModel.C
 $(IOModels)/noIO/noIO.C
 $(IOModels)/basicIO/basicIO.C
-$(IOModels)/tempIO/tempIO.C
-$(IOModels)/colorIO/colorIO.C
 $(IOModels)/trackIO/trackIO.C
 $(IOModels)/sophIO/sophIO.C
 
@@ -83,50 +70,46 @@ $(voidFractionModels)/voidFractionModel/voidFractionModel.C
 $(voidFractionModels)/voidFractionModel/newVoidFractionModel.C
 $(voidFractionModels)/centreVoidFraction/centreVoidFraction.C
 $(voidFractionModels)/dividedVoidFraction/dividedVoidFraction.C
-$(voidFractionModels)/dividedVoidFractionBiDi/dividedVoidFractionBiDi.C
 $(voidFractionModels)/dividedVoidFractionMS/dividedVoidFractionMS.C
 $(voidFractionModels)/bigParticleVoidFraction/bigParticleVoidFraction.C
 $(voidFractionModels)/GaussVoidFraction/GaussVoidFraction.C
 $(voidFractionModels)/IBVoidFraction/IBVoidFraction.C
-$(voidFractionModels)/weightedNeigbhorsVoidFraction/weightedNeigbhorsVoidFraction.C
 
 $(locateModels)/locateModel/locateModel.C
 $(locateModels)/locateModel/newLocateModel.C
 $(locateModels)/standardSearch/standardSearch.C
 $(locateModels)/engineSearch/engineSearch.C
-$(locateModels)/engineSearchMany2Many/engineSearchMany2Many.C
 $(locateModels)/turboEngineSearch/turboEngineSearch.C
-$(locateModels)/turboEngineSearchM2M/turboEngineSearchM2M.C
+$(locateModels)/engineSearchMany2Many/engineSearchMany2Many.C
 $(locateModels)/engineSearchIB/engineSearchIB.C
-$(locateModels)/hyperEngineSearch/hyperEngineSearch.C
-$(locateModels)/ijkSearch/ijkSearch.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
 $(dataExchangeModels)/twoWayFiles/twoWayFiles.C
 $(dataExchangeModels)/noDataExchange/noDataExchange.C
 $(dataExchangeModels)/twoWayMPI/twoWayMPI.C
-$(dataExchangeModels)/twoWayM2M/twoWayM2M.C
 $(dataExchangeModels)/twoWayMany2Many/twoWayMany2Many.C
 
 $(averagingModels)/averagingModel/averagingModel.C
 $(averagingModels)/averagingModel/newAveragingModel.C
 $(averagingModels)/dilute/dilute.C
 $(averagingModels)/dense/dense.C
-$(averagingModels)/denseBiDi/denseBiDi.C
 
 $(clockModels)/clockModel/clockModel.C
 $(clockModels)/clockModel/newClockModel.C
@@ -135,7 +118,6 @@ $(clockModels)/noClock/noClock.C
 
 $(liggghtsCommandModels)/liggghtsCommandModel/liggghtsCommandModel.C
 $(liggghtsCommandModels)/liggghtsCommandModel/newLiggghtsCommandModel.C
-$(liggghtsCommandModels)/colorParticles/colorParticles.C
 $(liggghtsCommandModels)/execute/execute.C
 $(liggghtsCommandModels)/runLiggghts/runLiggghts.C
 $(liggghtsCommandModels)/writeLiggghts/writeLiggghts.C
@@ -145,6 +127,5 @@ $(smoothingModels)/smoothingModel/smoothingModel.C
 $(smoothingModels)/smoothingModel/newSmoothingModel.C
 $(smoothingModels)/noSmoothing/noSmoothing.C
 $(smoothingModels)/constDiffSmoothing/constDiffSmoothing.C
-$(smoothingModels)/localPSizeDiffSmoothing/localPSizeDiffSmoothing.C
 
 LIB = $(CFDEM_LIB_DIR)/lib$(CFDEM_LIB_NAME)

src/lagrangian/cfdemParticle/Make/options

@@ -1,6 +1,9 @@
 sinclude $(GENERAL_RULES)/mplib$(WM_MPLIB)
 sinclude $(RULES)/mplib$(WM_MPLIB)
 
+GIT_VERSION := $(shell git describe --abbrev=4 --dirty --always --tags)
+PFLAGS+= -DGITVERSION=\"$(GIT_VERSION)\"
+
 
 EXE_INC = \
      $(PFLAGS) \
@@ -14,12 +17,11 @@ EXE_INC = \
     -I$(LIB_SRC)/OpenFOAM/containers/HashTables/labelHashSet \
     -I$(CFDEM_LIGGGHTS_SRC_DIR) \
     -I$(CFDEM_M2MLIB_PATH) \
-    -I$(CFDEM_Many2ManyLIB_PATH) \
     -I$(CFDEM_SRC_DIR)/cfdTools \
 
 LIB_LIBS = \
      $(PLIBS) \
-    -L$(FOAM_USER_LIBBIN) \
+    -L$(CFDEM_LIB_DIR) \
     -lfiniteVolume \
     -lincompressibleRASModels \
     -lincompressibleLESModels \
@@ -27,12 +29,6 @@ LIB_LIBS = \
     -llagrangian \
     -lmpi_cxx \
     -L$(CFDEM_LIGGGHTS_SRC_DIR) \
-    -l$(CFDEM_LIGGGHTS_LIB_NAME) \
-    -L$(CFDEM_M2MLIB_PATH) \
-    -lcouple \
+    -Wl,--whole-archive -l$(CFDEM_LIGGGHTS_LIB_NAME) -Wl,--no-whole-archive \
     -L$(CFDEM_Many2ManyLIB_PATH) \
-    -lcoupleMany2Many \
-
-/* add -I$(CFDEM_POEMSLIB_PATH) \ to EXE_INC */
-/*    -L$(CFDEM_POEMSLIB_PATH) \ */
-/*    -lpoems   */
+    -lcoupleMany2Many

src/lagrangian/cfdemParticle/cfdTools/checkModelType.H

@@ -4,8 +4,7 @@
     word modelType = particleCloud.modelType();
 
     //Warning << "model type not being checked" << endl;
-
-    if (modelType=="B"){
+    if (modelType=="Bfull"){
         Info << "solving volume averaged Navier Stokes equations of type B\n"<< endl;
 
         // check if Archimedes is used
@@ -18,6 +17,41 @@
         if(!found)
             FatalError <<"Archimedes model not found!\n" << abort(FatalError);
 
+        // check if gradPForce is used
+        found=false;
+        forAll(particleCloud.forceModels(),i)
+        {
+            if(particleCloud.forceModels()[i]=="gradPForce")
+                found=true;
+        }
+        if(!found)
+            FatalError <<"gradPForce model not found!\n" << abort(FatalError);
+
+        // check if viscForce is used
+        found=false;
+        forAll(particleCloud.forceModels(),i)
+        {
+            if(particleCloud.forceModels()[i]=="viscForce")
+                found=true;
+        }
+        if(!found)
+            FatalError <<"viscForce model not found!\n" << abort(FatalError);
+
+    }else if(modelType=="B"){
+        Info << "solving volume averaged Navier Stokes equations of type B\n"<< endl;
+
+        // check if Archimedes is used
+        bool found=false;
+        forAll(particleCloud.forceModels(),i)
+        {
+            if(particleCloud.forceModels()[i]=="Archimedes")
+                found=true;
+        }
+        if(!found)
+            FatalError <<"Archimedes model not found!\n" << abort(FatalError);
+
+
+
         // check if gradP and viscForce are used
         found=false;
         forAll(particleCloud.forceModels(),i)

src/lagrangian/cfdemParticle/cfdTools/global.C

@@ -0,0 +1,80 @@
+/*---------------------------------------------------------------------------*\
+    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(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "global.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(global, 0);
+
+defineRunTimeSelectionTable(global, dictionary);
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+void global::info()
+{
+    Info << "\nYou are currently using:" << endl;
+    Info << "OF version: " << FOAMversion << endl;
+    Info << "OF build: " << FOAMbuild << endl;
+    Info << "CFDEM build: " << CFDEMversion << "\n" << endl;
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+global::global
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    dict_(dict),
+    particleCloud_(sm),
+    CFDEMversion(GITVERSION)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+global::~global()
+{}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/cfdTools/global.H

@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+    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(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
+
+Class
+    global
+
+SourceFiles
+    global.Cver
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef global_H
+#define global_H
+
+#include "fvCFD.H"
+#include "cfdemCloud.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class global Declaration
+\*---------------------------------------------------------------------------*/
+
+class global
+{
+
+protected:
+
+    // Protected data
+    const dictionary& dict_;
+
+    cfdemCloud& particleCloud_;
+
+    const char* const CFDEMversion;
+
+
+    // Protected member functions
+
+public:
+
+    //- Runtime type information
+    TypeName("global");
+
+    // Declare runtime constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            global,
+            dictionary,
+            (
+                const dictionary& dict,
+                cfdemCloud& sm
+            ),
+            (dict,sm)
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        global
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+
+    // Destructor
+
+        virtual ~global();
+
+
+    // Selector
+
+        static autoPtr<global> New
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+    // Member Function
+
+    void info();
+
+    // Access
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/cfdTools/mathExtra.H

@@ -0,0 +1,265 @@
+/* ----------------------------------------------------------------------
+    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(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
+    
+    Copyright of this contribution:
+    Copyright 2014-     TU Graz, IPPT
+------------------------------------------------------------------------- */
+
+#ifndef CFDEM_MATH_EXTRA_H
+#define CFDEM_MATH_EXTRA_H
+
+#include <vector>
+//#include "math.h"
+#include "stdio.h"
+#include "string.h"
+#include "error.h"
+#include "ctype.h"
+
+#define TOLERANCE_ORTHO 1e-10
+
+namespace MathExtra 
+{
+
+//  inline void     outerProduct(double *vec1, double *vec2, double **m);
+//  inline double  spheroidGeometry(int index, double bi, double ai);
+
+
+
+//--------------------------------------------------------------------
+//   Outer Product of two vectors
+inline void outerProduct(double *vec1, double *vec2, double **m)
+{
+  int i, j;
+  //debug output
+  for( i = 0; i < 3; ++i )
+//    printf("OUTER PRODUCT: Input: vec1 element %d = %g", i, vec1[i]);
+  for( i = 0; i < 3; ++i )
+//    printf("OUTER PRODUCT: Input: vec2 element %d=%g", i, vec2[i]);
+  
+  //calculation
+  for( i = 0; i < 3; ++i )
+    for( j = 0; j < 3; ++j )
+    {
+      m[i][j] = vec1[i] * vec2[j];
+      printf("OUTER PRODUCT: Result: m[%d][%d]=%g", i, j, m[i][j]);
+    }
+    
+    
+}
+//--------------------------------------------------------------------
+// Compute the major, minor axis and eccentricity parameters of a prolate spheroid
+inline bool spheroidGeometry(double radius, double aspectRatio,  //inputs
+                                 double& ai,    double& bi,              //outputs
+                                 double& ei,    double& Le              //outputs
+                                ) //
+{
+    //INPUT
+    //  radius      ...volume-equivalent radius of the spheroid
+    //  aspectRatio ...major/minor aspect ratio
+    
+    //OUTPUT
+    //  ai  ... 
+    //  bi  ... 
+    //  ei  ... 
+    //  Le  ... 
+  
+    if(radius<=0.0)     //avoid troubles in case radius is 0 or negative
+        return false;
+
+ 	ai = radius * std::pow(aspectRatio*aspectRatio,0.33333333333333333333333);
+    bi = ai / aspectRatio;
+ 	ei = std::sqrt( 
+ 	                1.0
+                   - 1.0 / (aspectRatio*aspectRatio)
+ 	              );
+ 	Le = std::log(
+ 		           (1.0+ei)
+ 		          /(1.0-ei)
+                 );
+
+	return true; 		    
+}
+
+//--------------------------------------------------------------------
+// Compute the major, minor axis and eccentricity parameters of a prolate spheroid
+inline double Pi()
+{
+    return 3.1415926535897932384626433832795;
+}
+
+
+//--------------------------------------------------------------------
+// Compute the major, minor axis and eccentricity parameters of a prolate spheroid
+inline bool spheroidGeometry2(double radius, double aspectRatio,  //inputs
+                                 double& ai,     double& bi,           //outputs
+                                 double& XAe,    double& YAe,         //outputs
+                                 double& XCe,    double& YCe,         //outputs
+                                 double& YHe
+                                ) //
+{
+    //INPUT
+    //  radius      ...volume-equivalent radius of the spheroid
+    //  aspectRatio ...major/minor aspect ratio
+    
+    //OUTPUT
+    //  XAe  ...Eccentricity dependet parameter
+    //  YAe  ...Eccentricity dependet parameter
+    //  XCe  ...Eccentricity dependet parameter
+    //  XCe  ...Eccentricity dependet parameter
+    //  YCe  ...Eccentricity dependet parameter
+    //  YHe  ...Eccentricity dependet parameter
+  
+
+    double ei(0.0), Le(0.0);  
+    bool   result = 
+          spheroidGeometry(radius, aspectRatio,  //inputs
+                                 ai,     bi,              //outputs
+                                 ei,    Le              //outputs
+                              );
+    if(!result)
+        return false;
+        
+    XAe= 2.6666666666666666666666667
+         *ei*ei*ei
+         /(-2.0*ei+(1.0+ei*ei)*Le); 
+    YAe= 5.333333333333333333333333333
+        *ei*ei*ei
+        /(2.0*ei+(3*ei*ei-1.0)*Le);
+    XCe= 1.333333333333333333333333333
+        *ei*ei*ei
+        *(1.0-ei*ei)
+        /(2.0*ei-(1.0-ei*ei)*Le);
+    YCe= 1.3333333333333333333333
+        *ei*ei*ei
+        *(2.0-ei*ei)
+        /(-2.0*ei+(1.0+ei*ei)*Le);
+    YHe= 1.3333333333333333333333
+        *ei*ei*ei*ei*ei
+        /(-2.0*ei+(1.0+ei*ei)*Le);
+
+	return true; 		    
+
+}
+
+//--------------------------------------------------------------------
+// zeroize a 3x3x3 tensor
+inline void zeroize333(double tensor[3][3][3] )
+{
+    for(int iX=0; iX<3; iX++)
+        for(int iY=0; iY<3; iY++)
+            for(int iZ=0; iZ<3; iZ++)
+                tensor[iX][iY][iZ] = 0.0;
+}
+
+//--------------------------------------------------------------------
+// zeroize a 3x3 tensor
+inline void zeroize33(double tensor[3][3] )
+{
+    for(int iX=0; iX<3; iX++)
+        for(int iY=0; iY<3; iY++)
+                tensor[iX][iY] = 0.0;
+}
+
+//--------------------------------------------------------------------
+// multiply a 3x3x3 tensor with a scalar
+inline void multiply333(double scalar, double tensor[3][3][3] )
+{
+    for(int iX=0; iX<3; iX++)
+        for(int iY=0; iY<3; iY++)
+            for(int iZ=0; iZ<3; iZ++)
+                tensor[iX][iY][iZ] *= scalar;
+}
+//--------------------------------------------------------------------
+// Compute a dot and dyadic product of with a vector
+inline void permutationTensor(double tensor[3][3][3] )
+{
+    zeroize333(tensor);
+    tensor[0][1][2] = 1.0; 
+    tensor[1][2][0] = 1.0; 
+    tensor[2][0][1] = 1.0; 
+    tensor[0][2][1] =-1.0; 
+    tensor[2][1][0] =-1.0; 
+    tensor[1][0][2] =-1.0; 
+}
+
+
+
+//--------------------------------------------------------------------
+// Compute a dot product of the permutation tensor and 
+// then a dyadic product of with a vector
+inline bool permutationDotDyadic( 
+                                 double vector[3], 
+                                 double  tensor[3][3][3]      
+                                )
+{
+    //Generate permutation tensor
+    double permutationT[3][3][3];
+    permutationTensor(permutationT);
+    
+    //Step 1: compute dot prodcut of permutation tensor
+    double permutationDotProd[3][3];
+    zeroize33(permutationDotProd);
+
+    for(int iX=0; iX<3; iX++)
+        for(int iY=0; iY<3; iY++)
+            for(int iZ=0; iZ<3; iZ++)
+                permutationDotProd[iX][iY] += permutationT[iX][iY][iZ]
+                                            * vector[iZ];
+
+    for(int iX=0; iX<3; iX++)
+        for(int iY=0; iY<3; iY++)
+            for(int iZ=0; iZ<3; iZ++)
+                tensor[iX][iY][iZ] = permutationDotProd[iX][iY]
+                                   * vector[iZ]; 
+    return true;
+
+}
+
+//--------------------------------------------------------------------
+// Compute a dot and dyadic product of with a vector
+inline bool doubleDotTensor333Tensor33(double tensor333[3][3][3], 
+                                       double tensor33[3][3],
+                                       double result[3]
+                                      )
+{
+    result[0]=0.0;result[1]=0.0;result[2]=0.0;
+    
+    for(int iX=0; iX<3; iX++)
+        for(int iY=0; iY<3; iY++)
+            for(int iZ=0; iZ<3; iZ++)
+                result[iX] += tensor333[iX][iY][iZ] * tensor33[iY][iZ];
+
+    return true;
+}
+
+
+}; //end of namespace
+
+#endif

src/lagrangian/cfdemParticle/cfdTools/newGlobal.C

@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+    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(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "global.H"
+#include "dilute.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+autoPtr<global> global::New
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+{
+    word globalType
+    (
+        dict.lookup("global")
+    );
+
+    Info<< "Selecting global "
+         << globalType << endl;
+
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(globalType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalError
+            << "global::New(const dictionary&, const spray&) : "
+            << endl
+            << "    unknown globalType type "
+            << globalType
+            << ", constructor not in hash table" << endl << endl
+            << "    Valid global types are :"
+            << endl;
+        Info<< dictionaryConstructorTablePtr_->toc()
+            << abort(FatalError);
+    }
+
+    return autoPtr<global>(cstrIter()(dict,sm));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/cfdTools/versionInfo.H

@@ -1,9 +1,44 @@
-word CFDEMversion="cfdem-2.6.4";
-word compatibleLIGGGHTSversion="3.0.1";
-word OFversion="2.2.x-commit-61b850bc107bdd60bbf1bf9a6417b9faf701d128";
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
 
-Info << "\nCFDEMcoupling version: " << CFDEMversion << "\n" <<  endl;
-Info << "\n, compatible to LIGGGHTS version: " << compatibleLIGGGHTSversion << "\n" <<  endl;
-Info << "\n, compatible to OF version: " << OFversion << "\n" <<  endl;
+    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(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef versionInfo_H
+#define versionInfo_H
+
+word CFDEMversion="cfdem-2.8.2";
+word compatibleLIGGGHTSversion="3.0.6";
+word OFversion="2.3.x-commit-4d6f4a3115ff76ec4154c580eb041bc95ba4ec09";
+
+Info << "\nCFDEMcoupling version: " << CFDEMversion <<  endl;
+Info << ", compatible to LIGGGHTS version: " << compatibleLIGGGHTSversion <<  endl;
+Info << ", compatible to OF version and build: " << OFversion <<  endl;
+
+#endif

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C

@@ -31,6 +31,7 @@ Description
 
 #include "fileName.H"
 #include "cfdemCloud.H"
+#include "global.H"
 #include "forceModel.H"
 #include "locateModel.H"
 #include "momCoupleModel.H"
@@ -80,6 +81,7 @@ Foam::cfdemCloud::cfdemCloud
     positions_(NULL),
     velocities_(NULL),
     fluidVel_(NULL),
+    fAcc_(NULL),
     impForces_(NULL),
     expForces_(NULL),
     DEMForces_(NULL),
@@ -89,6 +91,7 @@ Foam::cfdemCloud::cfdemCloud
     cellIDs_(NULL),
     particleWeights_(NULL),
     particleVolumes_(NULL),
+    particleV_(NULL),
     numberOfParticles_(0),
     numberOfParticlesChanged_(false),
     arraysReallocated_(false),
@@ -99,7 +102,9 @@ Foam::cfdemCloud::cfdemCloud
     cg_(1.),
     cgOK_(true),
     impDEMdrag_(false),
+    impDEMdragAcc_(false),
     imExSplitFactor_(1.0),
+    treatVoidCellsAsExplicitForce_(false),
     useDDTvoidfraction_(false),
     ddtVoidfraction_
     (   
@@ -117,7 +122,7 @@ Foam::cfdemCloud::cfdemCloud
     turbulence_
     (
         #if defined(version21) || defined(version16ext)
-            #ifdef comp
+            #ifdef compre
                 mesh.lookupObject<compressible::turbulenceModel>
             #else
                 mesh.lookupObject<incompressible::turbulenceModel>
@@ -213,12 +218,16 @@ Foam::cfdemCloud::cfdemCloud
     )
 {
     #include "versionInfo.H"
+    global buildInfo(couplingProperties_,*this);
+    buildInfo.info();
 
     Info << "If BC are important, please provide volScalarFields -imp/expParticleForces-" << endl;
     if (couplingProperties_.found("solveFlow"))
         solveFlow_=Switch(couplingProperties_.lookup("solveFlow"));
     if (couplingProperties_.found("imExSplitFactor"))
         imExSplitFactor_ = readScalar(couplingProperties_.lookup("imExSplitFactor"));
+    if (couplingProperties_.found("treatVoidCellsAsExplicitForce"))
+        treatVoidCellsAsExplicitForce_ = readBool(couplingProperties_.lookup("treatVoidCellsAsExplicitForce"));
     if (couplingProperties_.found("verbose")) verbose_=true;
     if (couplingProperties_.found("ignore")) ignore_=true;
     if (turbulenceModelType_=="LESProperties")
@@ -276,6 +285,7 @@ Foam::cfdemCloud::~cfdemCloud()
     dataExchangeM().destroy(positions_,3);
     dataExchangeM().destroy(velocities_,3);
     dataExchangeM().destroy(fluidVel_,3);
+    dataExchangeM().destroy(fAcc_,3);
     dataExchangeM().destroy(impForces_,3);
     dataExchangeM().destroy(expForces_,3);
     dataExchangeM().destroy(DEMForces_,3);
@@ -285,6 +295,7 @@ Foam::cfdemCloud::~cfdemCloud()
     dataExchangeM().destroy(cellIDs_,1);
     dataExchangeM().destroy(particleWeights_,1);
     dataExchangeM().destroy(particleVolumes_,1);
+    dataExchangeM().destroy(particleV_,1);
 }
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
 void Foam::cfdemCloud::getDEMdata()
@@ -292,6 +303,9 @@ void Foam::cfdemCloud::getDEMdata()
     dataExchangeM().getData("radius","scalar-atom",radii_);
     dataExchangeM().getData("x","vector-atom",positions_);
     dataExchangeM().getData("v","vector-atom",velocities_);
+
+    if(impDEMdragAcc_)
+        dataExchangeM().getData("dragAcc","vector-atom",fAcc_); // array is used twice - might be necessary to clean it first
 }
 
 void Foam::cfdemCloud::giveDEMdata()
@@ -302,6 +316,7 @@ void Foam::cfdemCloud::giveDEMdata()
 
         if(impDEMdrag_)
         {
+            if(verbose_) Info << "sending Ksl and uf" << endl;
             dataExchangeM().giveData("Ksl","scalar-atom",Cds_);
             dataExchangeM().giveData("uf","vector-atom",fluidVel_);
         }
@@ -391,11 +406,6 @@ label Foam::cfdemCloud::particleCell(int index)
     return cellI;
 }
 
-double Foam::cfdemCloud::d(int index)
-{
-    return 2*radii()[index][0];
-}
-
 vector Foam::cfdemCloud::position(int index)
 {
     vector pos;
@@ -427,12 +437,6 @@ int Foam::cfdemCloud::nrForceModels()
     return forceModels_.size();
 }
 
-scalar Foam::cfdemCloud::radius(int index)
-{
-    scalar r = radii()[index][0];
-    return r;
-}
-
 scalar Foam::cfdemCloud::voidfraction(int index)
 {
     return voidfractions()[index][0];
@@ -452,6 +456,16 @@ label Foam::cfdemCloud::liggghtsCommandModelIndex(word name)
     return index;
 }
 
+std::vector<double*>* Foam::cfdemCloud::getVprobe()
+{
+ return probeModel_->getVprobe();
+}
+
+std::vector<double>* Foam::cfdemCloud::getSprobe()
+{
+ return probeModel_->getSprobe();
+}
+
 // * * * * * * * * * * * * * * * WRITE  * * * * * * * * * * * * * //
 
 // * * *   write cfdemCloud internal data   * * * //
@@ -466,6 +480,7 @@ bool Foam::cfdemCloud::evolve
     numberOfParticlesChanged_ = false;
     arraysReallocated_=false;
     bool doCouple=false;
+    probeModel_->clearProbes();
 
     if(!ignore())
     {
@@ -507,7 +522,7 @@ bool Foam::cfdemCloud::evolve
             // set void fraction field
             clockM().start(19,"setvoidFraction");
             if(verbose_) Info << "- setvoidFraction()" << endl;
-            voidFractionM().setvoidFraction(NULL,voidfractions_,particleWeights_,particleVolumes_);
+            voidFractionM().setvoidFraction(NULL,voidfractions_,particleWeights_,particleVolumes_,particleV_);
             if(verbose_) Info << "setvoidFraction done." << endl;
             clockM().stop("setvoidFraction");
 
@@ -519,7 +534,12 @@ bool Foam::cfdemCloud::evolve
             //Smoothen "next" fields            
             smoothingM().dSmoothing();
             smoothingM().smoothen(voidFractionM().voidFractionNext());
-            smoothingM().smoothenReferenceField(averagingM().UsNext());
+
+            //only smoothen if we use implicit force coupling in cells void of particles
+            //because we need unsmoothened Us field to detect cells for explicit 
+            //force coupling
+            if(!treatVoidCellsAsExplicitForce())
+                smoothingM().smoothenReferenceField(averagingM().UsNext());
             
             clockM().stop("setVectorAverage");
         }
@@ -595,15 +615,40 @@ bool Foam::cfdemCloud::reAllocArrays() const
         dataExchangeM().allocateArray(positions_,0.,3);
         dataExchangeM().allocateArray(velocities_,0.,3);
         dataExchangeM().allocateArray(fluidVel_,0.,3);
+        dataExchangeM().allocateArray(fAcc_,0.,3);
         dataExchangeM().allocateArray(impForces_,0.,3);
         dataExchangeM().allocateArray(expForces_,0.,3);
         dataExchangeM().allocateArray(DEMForces_,0.,3);
         dataExchangeM().allocateArray(Cds_,0.,1);
         dataExchangeM().allocateArray(radii_,0.,1);
         dataExchangeM().allocateArray(voidfractions_,1.,voidFractionM().maxCellsPerParticle());
-        dataExchangeM().allocateArray(cellIDs_,0.,voidFractionM().maxCellsPerParticle());
+        dataExchangeM().allocateArray(cellIDs_,-1.,voidFractionM().maxCellsPerParticle());
         dataExchangeM().allocateArray(particleWeights_,0.,voidFractionM().maxCellsPerParticle());
         dataExchangeM().allocateArray(particleVolumes_,0.,voidFractionM().maxCellsPerParticle());
+        dataExchangeM().allocateArray(particleV_,0.,1);
+        arraysReallocated_ = true;
+        return true;
+    }
+    return false;
+}
+
+bool Foam::cfdemCloud::reAllocArrays(int nP, bool forceRealloc) const
+{
+    if( (numberOfParticlesChanged_ && !arraysReallocated_) || forceRealloc)
+    {
+        // get arrays of new length
+        dataExchangeM().allocateArray(positions_,0.,3,nP);
+        dataExchangeM().allocateArray(velocities_,0.,3,nP);
+        dataExchangeM().allocateArray(fluidVel_,0.,3,nP);
+        dataExchangeM().allocateArray(impForces_,0.,3,nP);
+        dataExchangeM().allocateArray(expForces_,0.,3,nP);
+        dataExchangeM().allocateArray(DEMForces_,0.,3,nP);
+        dataExchangeM().allocateArray(Cds_,0.,1,nP);
+        dataExchangeM().allocateArray(radii_,0.,1,nP);
+        dataExchangeM().allocateArray(voidfractions_,1.,voidFractionM().maxCellsPerParticle(),nP);
+        dataExchangeM().allocateArray(cellIDs_,0.,voidFractionM().maxCellsPerParticle(),nP);
+        dataExchangeM().allocateArray(particleWeights_,0.,voidFractionM().maxCellsPerParticle(),nP);
+        dataExchangeM().allocateArray(particleVolumes_,0.,voidFractionM().maxCellsPerParticle(),nP);
         arraysReallocated_ = true;
         return true;
     }
@@ -615,7 +660,7 @@ tmp<fvVectorMatrix> cfdemCloud::divVoidfractionTau(volVectorField& U,volScalarFi
     return
     (
       - fvm::laplacian(voidfractionNuEff(voidfraction), U)
-      - fvc::div(voidfractionNuEff(voidfraction)*dev(fvc::grad(U)().T()))
+      - fvc::div(voidfractionNuEff(voidfraction)*dev2(fvc::grad(U)().T()))
     );
 }
 
@@ -648,11 +693,11 @@ void cfdemCloud::calcDdtVoidfraction(volScalarField& voidfraction) const
 
 tmp<volScalarField> cfdemCloud::voidfractionNuEff(volScalarField& voidfraction) const
 {
-    if (modelType_=="B")
+    if (modelType_=="B" || modelType_=="Bfull")
     {
         return tmp<volScalarField>
         (
-            #ifdef comp
+            #ifdef compre
                 new volScalarField("viscousTerm", (turbulence_.mut() + turbulence_.mu()))
             #else
                 new volScalarField("viscousTerm", (turbulence_.nut() + turbulence_.nu()))
@@ -663,7 +708,7 @@ tmp<volScalarField> cfdemCloud::voidfractionNuEff(volScalarField& voidfraction)
     {
         return tmp<volScalarField>
         (
-            #ifdef comp
+            #ifdef compre
                 new volScalarField("viscousTerm", voidfraction*(turbulence_.mut() + turbulence_.mu()))
             #else
                 new volScalarField("viscousTerm", voidfraction*(turbulence_.nut() + turbulence_.nu()))

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H

@@ -44,6 +44,7 @@ SourceFiles
 
 // choose version
 #include "OFversion.H"
+#include <vector>
 
 #include "fvCFD.H"
 #include "IFstream.H"
@@ -102,6 +103,8 @@ protected:
 
     mutable double **fluidVel_;
 
+    mutable double **fAcc_;
+
     mutable double **impForces_;
 
     mutable double **expForces_;
@@ -120,6 +123,8 @@ protected:
 
     mutable double **particleVolumes_;
 
+    mutable double **particleV_;
+
     int numberOfParticles_;
 
     bool numberOfParticlesChanged_;
@@ -140,14 +145,18 @@ protected:
 
     bool impDEMdrag_;
 
+    bool impDEMdragAcc_;
+
     mutable scalar imExSplitFactor_;
 
+    mutable bool treatVoidCellsAsExplicitForce_; //will treat the coupling force in cells with no Us data explicitly
+
     bool useDDTvoidfraction_;
 
     mutable volScalarField ddtVoidfraction_;
 
     #if defined(version21) || defined(version16ext)
-        #ifdef comp
+        #ifdef compre
             const compressible::turbulenceModel& turbulence_;
         #else
             const incompressible::turbulenceModel& turbulence_;
@@ -200,6 +209,7 @@ public:
     friend class dataExchangeModel;
     friend class voidFractionModel;
     friend class forceModel;
+    friend class forceSubModel;
 
 // Constructors
 
@@ -223,8 +233,6 @@ public:
 
        label particleCell(int);
 
-       double d(int);
-
        vector position(int);
 
        vector velocity(int);
@@ -235,8 +243,6 @@ public:
 
        virtual int nrForceModels();
 
-       scalar radius(int);
-
        scalar voidfraction(int);
 
        label liggghtsCommandModelIndex(word);
@@ -247,8 +253,12 @@ public:
 
        inline const bool& impDEMdrag() const;
 
+       inline const bool& impDEMdragAcc() const;
+
        inline const scalar& imExSplitFactor() const;
 
+       inline const bool& treatVoidCellsAsExplicitForce() const;
+
        inline const bool& ignore() const;
 
        inline const fvMesh& mesh() const;
@@ -265,6 +275,8 @@ public:
 
        inline double ** fluidVels() const;
 
+       inline double ** fAccs() const;
+
        inline double ** impForces() const;
 
        inline double ** expForces() const;
@@ -285,7 +297,13 @@ public:
 
        inline double ** particleWeights() const;
 
-       inline double ** particleVolumes() const;
+       virtual inline label body(int);
+
+       virtual inline double particleVolume(int);
+
+       inline scalar radius(int);
+
+       virtual inline double d(int);
 
        virtual inline double dMin() {return -1;}
        virtual inline double dMax() {return -1;}
@@ -297,10 +315,18 @@ public:
        virtual inline int ** particleTypes() const {return NULL;};
        virtual label particleType(label index) const {return -1;};
 
+       //access to the particle's rotation and torque data
+       virtual inline double ** DEMTorques() const {return NULL;};
+       virtual inline double ** omegaArray() const {return NULL;};
+       virtual vector omega(int) const {return Foam::vector(0,0,0);};
+
        //access to the particles' orientation information
        virtual inline double ** exArray() const {return NULL;};
        virtual vector ex(int) const {return Foam::vector(0,0,0);};
 
+       //Detector if SRF module is enable or not
+       virtual inline bool SRFOn(){return false;}
+
        inline int numberOfParticles() const;
 
        inline bool numberOfParticlesChanged() const;
@@ -336,7 +362,7 @@ public:
        inline autoPtr<liggghtsCommandModel>* liggghtsCommand() const;
 
        #if defined(version21) || defined(version16ext)
-            #ifdef comp
+            #ifdef compre
                 inline const compressible::turbulenceModel& turbulence() const;
             #else
                 inline const incompressible::turbulenceModel& turbulence() const;
@@ -352,6 +378,9 @@ public:
 
         virtual bool reAllocArrays() const;
 
+        virtual bool reAllocArrays(int nP, bool forceRealloc) const;  //force number of particles during reallocation
+
+
     // IO
         void writeScalarFieldToTerminal(double**&);
 
@@ -369,6 +398,10 @@ public:
         tmp<volScalarField> voidfractionNuEff(volScalarField&) const;
 
         void resetArray(double**&,int,int,double resetVal=0.);
+
+        std::vector<double*>* getVprobe();
+         
+        std::vector<double>* getSprobe();
 };
 
 

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H

@@ -55,11 +55,21 @@ inline const bool& cfdemCloud::impDEMdrag() const
     return impDEMdrag_;
 };
 
+inline const bool& cfdemCloud::impDEMdragAcc() const
+{ 
+    return impDEMdragAcc_;
+};
+
 inline const scalar& cfdemCloud::imExSplitFactor() const
 {
 	    return imExSplitFactor_;
 };
 
+inline const bool& cfdemCloud::treatVoidCellsAsExplicitForce() const
+{
+    return treatVoidCellsAsExplicitForce_;
+}
+
 inline const scalar& cfdemCloud::cg() const
 { 
     return cg_;
@@ -105,6 +115,11 @@ inline double ** cfdemCloud::fluidVels() const
     return fluidVel_;
 }
 
+inline double ** cfdemCloud::fAccs() const
+{
+    return fAcc_;
+}
+
 inline double ** cfdemCloud::impForces() const
 {
     return impForces_;
@@ -162,9 +177,24 @@ inline double ** cfdemCloud::particleWeights() const
     return particleWeights_;
 }
 
-inline double ** cfdemCloud::particleVolumes() const
+inline label Foam::cfdemCloud::body(int index)
 {
-    return particleVolumes_;
+    return index;
+}
+
+inline double cfdemCloud::particleVolume(int index)
+{
+    return particleV_[index][0];
+}
+
+inline scalar cfdemCloud::radius(int index)
+{
+    return radii_[index][0];
+}
+
+inline double cfdemCloud::d(int index)
+{
+    return 2*radii_[index][0];
 }
 
 inline int cfdemCloud::numberOfParticles() const
@@ -254,7 +284,7 @@ inline autoPtr<liggghtsCommandModel>* cfdemCloud::liggghtsCommand() const
 }
 
 #if defined(version21) || defined(version16ext)
-    #ifdef comp
+    #ifdef compre
         inline const compressible::turbulenceModel& cfdemCloud::turbulence() const
     #else
         inline const incompressible::turbulenceModel& cfdemCloud::turbulence() const

src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.C

@@ -121,7 +121,7 @@ bool Foam::cfdemCloudIB::evolve()
 
           // set void fraction field
           if(verbose_) Info << "- setvoidFraction()" << endl;
-          voidFractionM().setvoidFraction(NULL,voidfractions_,particleWeights_,particleVolumes_);
+          voidFractionM().setvoidFraction(NULL,voidfractions_,particleWeights_,particleVolumes_,particleV_);
           if(verbose_) Info << "setvoidFraction done." << endl;
         }
         
@@ -190,6 +190,7 @@ void Foam::cfdemCloudIB::calcVelocityCorrection
             }
         //}
     }
+    U.correctBoundaryConditions();
 
     // make field divergence free - set reference value in case it is needed
     fvScalarMatrix phiIBEqn

src/lagrangian/cfdemParticle/derived/cfdemCloudMS/cfdemCloudMS.C

@@ -56,21 +56,25 @@ cfdemCloudMS::cfdemCloudMS
     cellIDsCM_(NULL),
     bodies_(NULL),
     nrigids_(NULL),
-    exCM_(NULL),
-    eyCM_(NULL),
-    ezCM_(NULL),
-    VclumpCM_(NULL),
-    SclumpCM_(NULL),
-    scalingCM_(NULL),
-    typeCM_(NULL),
-    Cclump_ex_(NULL),
-    Cclump_ey_(NULL),
+    clumpType_(NULL),
+    nClumpTypes_(1),
+    clumpVol_(NULL),
+    clumpDH_(NULL),
+    clumpWeights_(NULL),
+    //exCM_(NULL),
+    //eyCM_(NULL),
+    //ezCM_(NULL),
+    //SclumpCM_(NULL),
+    //scalingCM_(NULL),
+    //Cclump_ex_(NULL),
+    //Cclump_ey_(NULL),
     impForcesCM_(NULL),
     expForcesCM_(NULL),
     DEMForcesCM_(NULL),
-    particleWeightsCM_(NULL),
     numberOfClumps_(-1),
     numberOfClumpsChanged_(false),
+    manDHdev_(false),
+    dHbyV_(scalarList(0)),
     useforcePerClump_(false),
     forceModels_(couplingProperties_.lookup("forceModelsMS"))
 {
@@ -84,6 +88,7 @@ cfdemCloudMS::cfdemCloudMS
             forceModels_[i]
         );
     }
+
 }
 
 
@@ -91,24 +96,25 @@ cfdemCloudMS::cfdemCloudMS
 
 cfdemCloudMS::~cfdemCloudMS()
 {
-    delete positionsCM_;
-    delete velocitiesCM_;
-    delete cellIDsCM_;
-    delete bodies_;
-    delete nrigids_;
-    delete exCM_;
-    delete eyCM_;
-    delete ezCM_;
-    delete VclumpCM_;
-    delete SclumpCM_;
-    delete scalingCM_;
-    delete typeCM_;
-    delete Cclump_ex_;
-    delete Cclump_ey_;
-    delete impForcesCM_;
-    delete expForcesCM_;
-    delete DEMForcesCM_;
-    delete particleWeightsCM_;
+    dataExchangeM().destroy(positionsCM_,3);
+    dataExchangeM().destroy(velocitiesCM_,3);
+    dataExchangeM().destroy(cellIDsCM_,1);
+    dataExchangeM().destroy(bodies_,1);
+    dataExchangeM().destroy(nrigids_,1);
+    dataExchangeM().destroy(clumpType_,1);
+    dataExchangeM().destroy(clumpVol_,1);
+    dataExchangeM().destroy(clumpDH_,1);
+    dataExchangeM().destroy(clumpWeights_,1);
+    //delete exCM_;
+    //delete eyCM_;
+    //delete ezCM_;
+    //delete SclumpCM_;
+    //delete scalingCM_;
+    //delete Cclump_ex_;
+    //delete Cclump_ey_;
+    dataExchangeM().destroy(impForcesCM_,3);
+    dataExchangeM().destroy(expForcesCM_,3);
+    dataExchangeM().destroy(DEMForcesCM_,3);
 }
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -117,32 +123,50 @@ cfdemCloudMS::~cfdemCloudMS()
 void cfdemCloudMS::getDEMdata()
 {
     cfdemCloud::getDEMdata();
-    // update NClumpTypes in data exch. model
-    //dataExchangeM().checkNClumpTypes();
 
     dataExchangeM().getData("xcm","vector-multisphere",positionsCM_);   // position of the centre of mass
     dataExchangeM().getData("vcm","vector-multisphere",velocitiesCM_);  // velocity of the centre of mass
     dataExchangeM().getData("body","scalar-atom",bodies_);              // clump-particle connex
-
     dataExchangeM().getData("nrigid","scalar-multisphere",nrigids_);    // # particles in clump
-    dataExchangeM().getData("ex_space","vector-multisphere",exCM_);     // axis of inertia
-    dataExchangeM().getData("ey_space","vector-multisphere",eyCM_);     // axis of inertia
-    dataExchangeM().getData("ez_space","vector-multisphere",ezCM_);     // axis of inertia
+    dataExchangeM().getData("clumptype","scalar-multisphere",clumpType_);   // type of the clump
 
-//    dataExchangeM().getScalarData("Vclump",VclumpCM_);   // Volume of the clump
-//    dataExchangeM().getScalarData("Sclump",SclumpCM_);   // surface area of the clump
+    nClumpTypes_=dataExchangeM().getNumberOfTypes();                        // nr of clump types
+    double* typeVol_=dataExchangeM().getTypeVol();                          // volume of the clump
 
-//    dataExchangeM().getScalarData("scaling",scalingCM_); // scaling of the clump
-//    dataExchangeM().getScalarData("typeCM",typeCM_);    // type of the clump
+    //- save clump volume and mass
+    double **typeDH(NULL);
+    dataExchangeM().allocateArray(typeDH,-1,1,nClumpTypes()+1);
+    if(manDHdev_) // use manually defined dH
+    {
+        for(int k = 1;k <= nClumpTypes(); k++)
+            typeDH[k][0]=dHbyV_[k-1]*typeVol_[k];
+    }
+    else // calc dH from volAeqivalent shpere
+    {  
+        for(int k = 1;k <= nClumpTypes(); k++)
+            typeDH[k][0]=pow(typeVol_[k]*1.9099,1./3.); // 6/pi=1.9099 // calc a hydraulic diameter as d of vol equal sphere
+    }
 
-//    dataExchangeM().getScalarData("Cclump_ex",Cclump_ex_);   // cross section of the clump in ex normal direction
-//    dataExchangeM().getScalarData("Cclump_ey",Cclump_ey_);   // cross section of the clump in ey normal direction
+    int ct(0);
+    for(int ind = 0;ind < numberOfClumps(); ind++)
+    {
+        ct=clumpType()[0][ind];
+        clumpVol_[ind][0] = typeVol_[ct];      
+        clumpDH_[ind][0]=typeDH[ct][0];
+        //Info << "ct=" << ct << endl;
+        //Info << "clumpVol()[ind][0]=" << clumpVol()[ind][0] << endl;
+        //Info << "clumpDH()[ind][0]=" << clumpDH()[ind][0] << endl;
+    }
+    dataExchangeM().destroy(typeDH,1);
+    // --
 
-    // calc Saequi, surface area of the vol aequivalent sphere
-      // Saequi=pi*pow(6*VclumpCM_/pi,2./3.);
-
-    // calc Caequi, cross section of the vol aequivalent sphere
-      // Caequi=Saequi/4;
+    //dataExchangeM().getData("ex_space","vector-multisphere",exCM_);     // axis of inertia
+    //dataExchangeM().getData("ey_space","vector-multisphere",eyCM_);     // axis of inertia
+    //dataExchangeM().getData("ez_space","vector-multisphere",ezCM_);     // axis of inertia
+    //dataExchangeM().getScalarData("Sclump",SclumpCM_);   // surface area of the clump
+    //dataExchangeM().getScalarData("scaling",scalingCM_); // scaling of the clump
+    //dataExchangeM().getScalarData("Cclump_ex",Cclump_ex_);   // cross section of the clump in ex normal direction
+    //dataExchangeM().getScalarData("Cclump_ey",Cclump_ey_);   // cross section of the clump in ey normal direction
 }
 
 void Foam::cfdemCloudMS::giveDEMdata()
@@ -174,27 +198,26 @@ bool cfdemCloudMS::reAllocArrays() const
 {
     if(cfdemCloud::reAllocArrays())
     {
-        int nClumpTypes = dataExchangeM().nClumpTypes();
-
         // get arrays of new length
         dataExchangeM().allocateArray(positionsCM_,0,3,"nbodies");
         dataExchangeM().allocateArray(velocitiesCM_,0,3,"nbodies");
-        dataExchangeM().allocateArray(cellIDsCM_,0,1,"nbodies");
-        dataExchangeM().allocateArray(bodies_,0,3);
-        dataExchangeM().allocateArray(nrigids_,0,1,"nbodies");
-        dataExchangeM().allocateArray(exCM_,0,3,"nbodies");
-        dataExchangeM().allocateArray(eyCM_,0,3,"nbodies");
-        dataExchangeM().allocateArray(ezCM_,0,3,"nbodies");
-        dataExchangeM().allocateArray(VclumpCM_,0,3,nClumpTypes);
-        dataExchangeM().allocateArray(SclumpCM_,0,3,nClumpTypes);
-        dataExchangeM().allocateArray(scalingCM_,0,3,"nbodies");
-        dataExchangeM().allocateArray(typeCM_,0,3,"nbodies");
-        dataExchangeM().allocateArray(Cclump_ex_,0,3,nClumpTypes);
-        dataExchangeM().allocateArray(Cclump_ey_,0,3,nClumpTypes);
+        dataExchangeM().allocateArray(cellIDsCM_,-1,1,"nbodies");
+        dataExchangeM().allocateArray(bodies_,0,1);
+        dataExchangeM().allocateArray(nrigids_,0,1,"nbodies");      
+        dataExchangeM().allocateArray(clumpType_,0,1,"nbodies");
+        dataExchangeM().allocateArray(clumpVol_,0,1,"nbodies");
+        dataExchangeM().allocateArray(clumpDH_,1.,1,"nbodies");
+        dataExchangeM().allocateArray(clumpWeights_,1,1,"nbodies");
+        //dataExchangeM().allocateArray(exCM_,0,3,"nbodies");
+        //dataExchangeM().allocateArray(eyCM_,0,3,"nbodies");
+        //dataExchangeM().allocateArray(ezCM_,0,3,"nbodies");
+        //dataExchangeM().allocateArray(SclumpCM_,0,3,nClumpTypes);
+        //dataExchangeM().allocateArray(scalingCM_,0,3,"nbodies");
+        //dataExchangeM().allocateArray(Cclump_ex_,0,3,nClumpTypes);
+        //dataExchangeM().allocateArray(Cclump_ey_,0,3,nClumpTypes);
         dataExchangeM().allocateArray(impForcesCM_,0,3,"nbodies");
         dataExchangeM().allocateArray(expForcesCM_,0,3,"nbodies");
         dataExchangeM().allocateArray(DEMForcesCM_,0,3,"nbodies");
-        dataExchangeM().allocateArray(particleWeightsCM_,1,1,"nbodies");        // filed is never changed-correct only for centre
         return true;
     }
     return false;
@@ -203,7 +226,6 @@ bool cfdemCloudMS::reAllocArrays() const
 void Foam::cfdemCloudMS::setNumberOfParticles(int nP)
 {
     cfdemCloud::setNumberOfParticles(nP);
-
     int nC = dataExchangeM().getNumberOfClumps();
 
     if(nC != numberOfClumps())
@@ -211,6 +233,9 @@ void Foam::cfdemCloudMS::setNumberOfParticles(int nP)
         numberOfClumpsChanged_ = true;
         numberOfClumps_ = nC;
     }
+
+    // in case last particle has left an ma-_tag_ms is not up to date
+    numberOfClumps_ = min(numberOfParticles(),numberOfClumps_);
 }
 
 void Foam::cfdemCloudMS::findCells()
@@ -246,14 +271,14 @@ void Foam::cfdemCloudMS::setParticleForceField()
         (
         forceM(0).impParticleForces(),
         impForcesCM_,
-        particleWeightsCM_,
+        clumpWeights_,
         numberOfClumps()
         );
         averagingM().setVectorSumSimple
         (
         forceM(0).expParticleForces(),
         expForcesCM_,
-        particleWeightsCM_,
+        clumpWeights_,
         numberOfClumps()
         );
     }
@@ -262,45 +287,11 @@ void Foam::cfdemCloudMS::setParticleForceField()
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 //             PUBLIC MEMBER FUNCTIONS
 
-vector Foam::cfdemCloudMS::positionCM(int index)
-{
-    vector pos;
-    for(int i=0;i<3;i++) pos[i] = positionsCM()[index][i];
-    return pos;
-}
-
-vector Foam::cfdemCloudMS::velocityCM(int index)
-{
-    vector vel;
-    for(int i=0;i<3;i++) vel[i] = velocitiesCM()[index][i];
-    return vel;
-}
-
-label Foam::cfdemCloudMS::cellIDCM(int index)
-{
-    return cellIDsCM_[index][0];
-}
-
-label Foam::cfdemCloudMS::body(int index)
-{
-    return bodies_[0][index]-1;
-}
-
-label Foam::cfdemCloudMS::nrigid(int index)
-{
-    return nrigids_[0][index];
-//    return nrigids_[0][0];
-}
-
 const forceModel& Foam::cfdemCloudMS::forceM(int i)
 {
     return forceModel_[i];
 }
 
-int Foam::cfdemCloudMS::nrForceModels()
-{
-    return forceModels_.size();
-}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/lagrangian/cfdemParticle/derived/cfdemCloudMS/cfdemCloudMS.H

@@ -59,37 +59,36 @@ class cfdemCloudMS
     public cfdemCloud
 {
 private:
-
     mutable double **positionsCM_;
     mutable double **velocitiesCM_;
     mutable double **cellIDsCM_;
     mutable int **bodies_;
     mutable int **nrigids_;
 
-    mutable double **exCM_;
-    mutable double **eyCM_;
-    mutable double **ezCM_;
+    mutable int **clumpType_;
+    mutable int nClumpTypes_;
+    mutable double **clumpVol_;
+    mutable double **clumpDH_;
+    mutable double **clumpWeights_;
 
-    mutable double **VclumpCM_;
-    mutable double **SclumpCM_;
-
-    mutable double **scalingCM_;
-    mutable double **typeCM_;
-
-    mutable double **Cclump_ex_;
-    mutable double **Cclump_ey_;
+    //mutable double **exCM_;
+    //mutable double **eyCM_;
+    //mutable double **ezCM_;
+    //mutable double **SclumpCM_;
+    //mutable double **scalingCM_;
+    //mutable double **Cclump_ex_;
+    //mutable double **Cclump_ey_;
 
     mutable double **impForcesCM_;
     mutable double **expForcesCM_;
     mutable double **DEMForcesCM_;
 
-    mutable double **particleWeightsCM_;
-
     int numberOfClumps_;
     bool numberOfClumpsChanged_;
+    Switch manDHdev_;
+    scalarList dHbyV_;
 
     bool useforcePerClump_;
-
     const wordList forceModels_;
     autoPtr<class forceModelMS>* forceModel_;
 
@@ -119,21 +118,32 @@ public:
 
     // Member Functions
       // Access
+       inline label body(int);
 
-       vector positionCM(int);
+       inline double particleVolume(int);
 
-       vector velocityCM(int);
+       inline vector positionCM(int);
 
-       label cellIDCM(int);
+       inline vector velocityCM(int);
 
-       label body(int);
+       inline label cellIDCM(int);
 
-       label nrigid(int);
+       inline label nrigid(int);
+
+       inline int nrForceModels();
 
        inline double **& positionsCM() const;
 
        inline double **& velocitiesCM() const;
 
+       inline int **& clumpType() const;
+
+       inline int nClumpTypes() const;
+
+       inline double **& clumpVol() const;
+
+       inline double **& clumpDH() const;
+
        inline double **& impForcesCM() const;
 
        inline double **& expForcesCM() const;
@@ -144,12 +154,16 @@ public:
 
        inline bool numberOfClumpsChanged() const;
 
+       inline void setManDHdev(Switch);
+
+       inline void setDHbyV(scalarList);
+
+       inline scalarList& dHbyV();
+
        inline bool useforcePerClump() const;
 
        const forceModel& forceM(int);
 
-       int nrForceModels();
-
        bool evolve(volScalarField&,volVectorField&,volVectorField&);
 };
 

src/lagrangian/cfdemParticle/derived/cfdemCloudMS/cfdemCloudMSI.H

@@ -34,6 +34,55 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+inline label Foam::cfdemCloudMS::body(int index)
+{
+    return bodies_[0][index]-1;
+}
+
+inline double cfdemCloudMS::particleVolume(int index)
+{
+    int ind = body(index); // particle to clump ID
+
+    double Vp(SMALL);
+    if(ind >= 0 && ind <numberOfClumps())
+        Vp=clumpVol_[ind][0];
+    else
+        Warning << "ind=" << ind << endl;
+    
+    int nR(nrigid(ind));
+    if(nR>0) Vp/=nR;
+    return Vp;
+}
+
+inline vector Foam::cfdemCloudMS::positionCM(int index)
+{
+    vector pos;
+    for(int i=0;i<3;i++) pos[i] = positionsCM()[index][i];
+    return pos;
+}
+
+inline vector Foam::cfdemCloudMS::velocityCM(int index)
+{
+    vector vel;
+    for(int i=0;i<3;i++) vel[i] = velocitiesCM()[index][i];
+    return vel;
+}
+
+inline label Foam::cfdemCloudMS::cellIDCM(int index)
+{
+    return cellIDsCM_[index][0];
+}
+
+inline label Foam::cfdemCloudMS::nrigid(int index)
+{
+    return nrigids_[0][index];
+}
+
+inline int Foam::cfdemCloudMS::nrForceModels()
+{
+    return forceModels_.size();
+}
+
 inline double **& cfdemCloudMS::positionsCM() const
 {
     return positionsCM_;
@@ -44,6 +93,26 @@ inline double **& cfdemCloudMS::velocitiesCM() const
     return velocitiesCM_;
 }
 
+inline int cfdemCloudMS::nClumpTypes() const
+{
+    return nClumpTypes_;
+}
+
+inline int **& cfdemCloudMS::clumpType() const
+{
+    return clumpType_;
+}
+
+inline double **& cfdemCloudMS::clumpVol() const
+{
+    return clumpVol_;
+}
+
+inline double **& cfdemCloudMS::clumpDH() const
+{
+    return clumpDH_;
+}
+
 inline double **& cfdemCloudMS::impForcesCM() const
 {
     return impForcesCM_;
@@ -69,6 +138,21 @@ inline bool cfdemCloudMS::numberOfClumpsChanged() const
     return numberOfClumpsChanged_;
 }
 
+inline void cfdemCloudMS::setManDHdev(Switch v)
+{
+    manDHdev_ = v;
+}
+
+inline void cfdemCloudMS::setDHbyV(scalarList v)
+{
+    dHbyV_ = v;
+}
+
+inline scalarList& cfdemCloudMS::dHbyV()
+{
+    return dHbyV_;
+}
+
 inline bool cfdemCloudMS::useforcePerClump() const
 {
     return useforcePerClump_;

src/lagrangian/cfdemParticle/etc/OFversion/OFversion.H

@@ -1,11 +1,11 @@
-//#define version23  // currently being tested
-#define version22  // currently being used
+#define version23  // currently being tested
+//#define version22  // currently being used
 //#define version21
 //#define version16ext
 //#define version15
 
 // choose comp/incomp
-//#define comp   // if comp is on - you must use Make/options_comp!
+//#define compre   // if comp is on - you must use Make/options_comp!
 
 //define multi sphere
 //#define multisphere

src/lagrangian/cfdemParticle/etc/bashrc

@@ -37,12 +37,25 @@
 #- export environment variables (adapt to your paths)
 #------------------------------------------------------------------------------
 
+#check if default lammps lib path should be used
+if [[ $CFDEM_LAMMPS_LIB_DIR == "" ]]; then
+    export CFDEM_LAMMPS_LIB_DIR=$CFDEM_LIGGGHTS_SRC_DIR/../lib/
+else
+    echo "using already defined CFDEM_LAMMPS_LIB_DIR=$CFDEM_LAMMPS_LIB_DIR ."
+fi
+
 #- LIGGGHTS lib name
 export CFDEM_LIGGGHTS_LIB_NAME=lmp_$CFDEM_LIGGGHTS_MAKEFILE_NAME
 
 #- CFDEM lib name
 export CFDEM_LIB_NAME=lagrangianCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
 
+#- CFDEM compressible lib name
+export CFDEM_LIB_COMP_NAME=lagrangianCFDEMcomp-$CFDEM_VERSION-$WM_PROJECT_VERSION
+
+#-----------------------------------------------------
+# additional libraries
+
 #- LMP Many2Many lib path and makefile
 export CFDEM_Many2ManyLIB_PATH=$CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library
 export CFDEM_Many2ManyLIB_MAKEFILENAME=$CFDEM_LIGGGHTS_MAKEFILE_NAME
@@ -52,9 +65,17 @@ export CFDEM_M2MLIB_PATH=$CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataE
 export CFDEM_M2MLIB_MAKEFILENAME=$CFDEM_LIGGGHTS_MAKEFILE_NAME
 
 #- LMP POEMS lib path and makefile
-export CFDEM_POEMSLIB_PATH=$CFDEM_LIGGGHTS_SRC_DIR/../lib/poems
+export CFDEM_POEMSLIB_PATH=$CFDEM_LAMMPS_LIB_DIR/poems
 export CFDEM_POEMSLIB_MAKEFILENAME=g++
 
+#- LMP ASPHERE lib path and makefile
+export CFDEM_ASPHERELIB_PATH=$CFDEM_LAMMPS_LIB_DIR/poems
+export CFDEM_ASPHERELIB_MAKEFILENAME=g++
+
+#-C3PO library
+export C3PO_SRC_DIR=$CFDEM_SRC_DIR/c3po
+#-----------------------------------------------------
+
 #- path to test harness
 export CFDEM_TEST_HARNESS_PATH=$CFDEM_PROJECT_USER_DIR/log/logFilesCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
 
@@ -100,6 +121,9 @@ alias cfdemUt='cd $CFDEM_UT_DIR'
 #- shortcut to run path
 alias cfdemRun='cd $CFDEM_PROJECT_USER_DIR/run'
 
+#- shortcut to user solver path
+alias cfdemUsrSol='cd $CFDEM_PROJECT_USER_DIR/applications/solvers'
+
 #- shortcut to documentation path
 alias cfdemDoc='cd $CFDEM_DOC_DIR'
 
@@ -142,6 +166,9 @@ alias cfdemCompCFDEMuti='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compil
 #- shortcut to test basic tutorials
 alias cfdemTestTUT='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/testTutorials.sh'
 
+#- shortcut to visualize the clock model data
+alias vizClock='python $CFDEM_UT_DIR/vizClock/matPlot.py'
+
 #- recursive touch of current directory
 alias touchRec='find ./* -exec touch {} \;'
 
@@ -157,6 +184,11 @@ export -f cfdemLiggghtsPar
 cfdemGrep() { grep -rl "$1" ./* | xargs gedit; }
 export -f cfdemGrep
 
+#- shortcut lo list files in a directory
+#cfdemListFiles() { find $1 | sed s:""$1"":: > listOfFiles.txt; }   #leave out the dir iteslf in list
+cfdemListFiles() { find $1 > listOfFiles.txt; }                     #keep the dir in list
+export -f cfdemListFiles
+
 # check if the run directory exists
 if [ -d "$CFDEM_PROJECT_USER_DIR" ]; then
     :
@@ -164,11 +196,11 @@ else
         echo "make new dirs $CFDEM_PROJECT_USER_DIR ? (y/n)"
         read YN
         if [ $YN == "y" ]; then
-            cd $CFDEM_PROJECT_DIR
-            mkdir $LOGNAME-$CFDEM_VERSION-$WM_PROJECT_VERSION
-            cd $LOGNAME-$CFDEM_VERSION-$WM_PROJECT_VERSION
+            mkdir -p $CFDEM_PROJECT_USER_DIR
+            cd $CFDEM_PROJECT_USER_DIR
             mkdir run
             mkdir -p log/logFilesCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
+            mkdir -p applications/solvers
         else
             echo "aborted by user."
             #exit

src/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_sol.sh

@@ -51,15 +51,16 @@ else
     logpath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")/$logDir"
 
     ##number of solvers compiled at a time
-    nsteps=$WM_NCOMPPROCS
-    echo "do compilation on $nsteps procs"
-    nchunk=`echo $njobs/$nsteps+1 | bc`
-    if [[ $WM_NCOMPPROCS == "" ]]; then
-        echo "do compilation in serial"
+
+    if [[ $WM_NCOMPPROCS == "" ]] || [ $WM_NCOMPPROCS -eq 1 ]; then
         nsteps=1
-        nchunk=1
-    else
-        echo "do compilation on $nsteps procs in $nchunk chunks"      
+        let nchunk=$njobs+1 # +1, to wait for the last compilation too
+        echo "do compilation in serial"
+    else    
+        nsteps=$WM_NCOMPPROCS
+        nchunk=`echo $njobs/$nsteps+1 | bc`
+        echo "do compilation on $nsteps procs in $nchunk chunks" 
+        let nchunk++ # +1, to wait for the last compilation too     
     fi
 
     counter=0
@@ -108,6 +109,8 @@ else
                 let counter++
             fi
         done
+
+        sleep 1 # wait a second until compilation starts
     done
 
     echo "compilation done."

src/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_src.sh

@@ -20,48 +20,48 @@ mkdir -p $logDir
 #================================================================================#
 # compile src
 #================================================================================#
-    whitelist="$CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/library-list.txt"
-    echo ""
-    echo "Please provide the libraries to be compiled in the $CWD/$whitelist file."
+whitelist="$CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/library-list.txt"
+echo ""
+echo "Please provide the libraries to be compiled in the $CWD/$whitelist file."
 
-    if [ ! -f "$CWD/$whitelist" ];then
-        echo "$whitelist does not exist in $CWD. Nothing will be done."
-        NLINES=0
-        COUNT=0
-    else
-        NLINES=`wc -l < $CWD/$whitelist`
-        COUNT=0
-    fi
+if [ ! -f "$CWD/$whitelist" ];then
+    echo "$whitelist does not exist in $CWD. Nothing will be done."
+    NLINES=0
+    COUNT=0
+else
+    NLINES=`wc -l < $CWD/$whitelist`
+    COUNT=0
+fi
 
-    while [ $COUNT -lt $NLINES ]
-    do
-            let COUNT++  
-            LINE=`head -n $COUNT $CWD/$whitelist | tail -1`
-  
-            # white lines
-            if [[ "$LINE" == "" ]]; then
-                echo "compile $LINE"
-                continue
-            # comments
-            elif [[ "$LINE" == \#* ]]; then
-                continue
-             # paths
-            elif [[ "$LINE" == */dir ]]; then
-                echo "will change path..."
-                LINE=$(echo "${LINE%????}")
-                path="$CFDEM_SRC_DIR/$LINE"
-                cd $path
-                #continue
-            fi
-            #--------------------------------------------------------------------------------#
-            #- define variables
-            logpath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")/$logDir"
-            logfileName="log_compileCFDEMcoupling_"$(basename $LINE)""
-            casePath="$path"
-            headerText="$logfileName""-$NOW"
-            #--------------------------------------------------------------------------------#
-            # remove old log file
-            rm "$logpath/$logfileName"*
-            compileLib $logpath $logfileName $casePath $headerText
-    done
+while [ $COUNT -lt $NLINES ]
+do
+        let COUNT++  
+        LINE=`head -n $COUNT $CWD/$whitelist | tail -1`
+
+        # white lines
+        if [[ "$LINE" == "" ]]; then
+            echo "compile $LINE"
+            continue
+        # comments
+        elif [[ "$LINE" == \#* ]]; then
+            continue
+         # paths
+        elif [[ "$LINE" == */dir ]]; then
+            echo "will change path..."
+            LINE=$(echo "${LINE%????}")
+            path="$CFDEM_SRC_DIR/$LINE"
+            cd $path
+            #continue
+        fi
+        #--------------------------------------------------------------------------------#
+        #- define variables
+        logpath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")/$logDir"
+        logfileName="log_compileCFDEMcoupling_"$(basename $LINE)""
+        casePath="$path"
+        headerText="$logfileName""-$NOW"
+        #--------------------------------------------------------------------------------#
+        # remove old log file
+        rm "$logpath/$logfileName"*
+        compileLib $logpath $logfileName $casePath $headerText
+done
 

src/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_uti.sh

@@ -1,36 +1,115 @@
 #!/bin/bash
 
 #===================================================================#
-# compile routine for CFDEMcoupling solvers, part of CFDEMproject 
+# compile routine for CFDEMcoupling utilities, part of CFDEMproject 
 # Christoph Goniva - May. 2012, DCS Computing GmbH
 #===================================================================#
 
+whitelist="utilities-list.txt"
+
 #- include functions
 source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
-
-NOW="$(date +"%Y-%m-%d-%H:%M")"
 logDir="log"
-
-
 cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
 mkdir -p $logDir
 
-#================================================================================#
-# compile utilities
-#================================================================================#
+CWD="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")"
+NOW="$(date +"%Y-%m-%d-%H:%M")"
 
+echo ""
+echo "This routine will compile the utilities specified in utilities-list.txt"
+echo ""
+#echo "Are the variables CFDEM_UT_DIR=$CFDEM_UT_DIR"
+#echo "and CFDEM_SRC_DIR=$CFDEM_SRC_DIR/lagrangian/cfdemParticle correct? (y/n)"
+#read YN
+#if [ "$YN" != "y" ];then
+#  	echo "Aborted by user."
+#  	exit 1
+#fi
+
+echo ""
+echo "Please provide the utilities to be compiled in the $CWD/$whitelist file."
+echo "structure:"
+echo "path  to provide the path relative to CFDEM_UT_DIR"
+echo ""
+echo "example:"
+echo "cfdemPostproc/dir"
+echo ""
+
+if [ ! -f "$CWD/$whitelist" ];then
+    echo "$whitelist does not exist in $CWD"
+else
+    njobs=`wc -l < $CWD/$whitelist` 
+    echo ""
+    echo "running compilation in pseudo-parallel mode of $njobs utilities"
 
-for utName in "cfdemPostproc"
-do
     #--------------------------------------------------------------------------------#
-    #- define variables
     logpath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")/$logDir"
-    logfileName="log_compileCFDEMcoupling""_$utName" 
-    casePath="$CFDEM_UT_DIR/$utName"
-    headerText="$logfileName""_$utName""-$NOW"
-    #--------------------------------------------------------------------------------#
-    compileSolver $logpath $logfileName $casePath $headerText
-done
 
-echo "Note: the list of utilities compiled might be incomplete."
-echo "please check $CFDEM_UT_DIR for more utilities available"
+    ##number of utilities compiled at a time
+
+    if [[ $WM_NCOMPPROCS == "" ]] || [ $WM_NCOMPPROCS -eq 1 ]; then
+        nsteps=1
+        let nchunk=$njobs+1 # +1, to wait for the last compilation too
+        echo "do compilation in serial"
+    else    
+        nsteps=$WM_NCOMPPROCS
+        nchunk=`echo $njobs/$nsteps+1 | bc`
+        echo "do compilation on $nsteps procs in $nchunk chunks" 
+        let nchunk++ # +1, to wait for the last compilation too     
+    fi
+
+    counter=0
+    for i in `seq $nchunk`
+    do
+
+        #wait until prev. compilation is finished
+        echo "waiting..."
+        until [ `ps -C make | wc -l` -eq 1 ]; 
+        do 
+            sleep 2
+        done
+
+        for j in `seq $nsteps`
+        do
+            let solNr=($i-1)*$nsteps+$j
+            LINE=`head -n $solNr $CWD/$whitelist | tail -1`
+
+            # white lines
+            if [[ "$LINE" == "" ]]; then
+                continue
+            # comments
+            elif [[ "$LINE" == \#* ]]; then
+                continue
+            # paths
+            elif [[ "$LINE" == */dir ]]; then
+            #echo "change path"
+                LINE=$(echo "${LINE%????}")
+                path="$CFDEM_UT_DIR/$LINE"
+                #cd $path
+                let solNr++
+            fi
+
+            if [[ "$counter" -lt "$njobs" ]]; then
+                #--------------------------------------------------------------------------------#
+                #- define variables
+                #logpath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")/$logDir"
+                logfileName="log_compileCFDEMcoupling""_$LINE" 
+                casePath="$CFDEM_UT_DIR/$LINE"
+                headerText="$logfileName""_$LINE""-$NOW"
+                parallel="true"
+                #--------------------------------------------------------------------------------#
+
+                #echo "compiling $LINE"
+                compileSolver $logpath $logfileName $casePath $headerText $parallel
+                let counter++
+            fi
+        done
+
+        sleep 1 # wait a second until compilation starts
+    done
+
+    echo "compilation done."
+fi
+
+

src/lagrangian/cfdemParticle/etc/compileLIGGGHTS.sh

@@ -14,6 +14,12 @@ logDir="log"
 cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
 mkdir -p $logDir
 
+#================================================================================#
+# copy LIGGGHTS patch files if available
+#================================================================================#
+echo "copying patch files for LIGGGHTS if available"
+cp $CFDEM_SRC_DIR/LIGGGHTSpatch/* $CFDEM_LIGGGHTS_SRC_DIR
+
 #--------------------------------------------------------------------------------#
 #- define variables
 logpath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")/$logDir"

src/lagrangian/cfdemParticle/etc/compileLIGGGHTS_lib.sh

@@ -33,6 +33,8 @@ mkdir -p $logDir
         COUNT=0
     fi
 
+    logpath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")/$logDir"
+
     while [ $COUNT -lt $NLINES ]
     do
             let COUNT++  
@@ -57,7 +59,7 @@ mkdir -p $logDir
 
             #--------------------------------------------------------------------------------#
             #- define variables
-            logpath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")/$logDir"
+            #logpath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")/$logDir"
             logfileName="log_compile$LINE""lib"
             headerText="$logfileName""-$NOW"
             libVarMakefileName="CFDEM_$LINE""LIB_MAKEFILENAME"

src/lagrangian/cfdemParticle/etc/cshrc

@@ -43,6 +43,9 @@ setenv CFDEM_LIGGGHTS_LIB_NAME lmp_$CFDEM_LIGGGHTS_MAKEFILE_NAME
 #- CFDEM lib name
 setenv CFDEM_LIB_NAME lagrangianCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
 
+#- CFDEM compressible lib name
+setenv CFDEM_LIB_COMP_NAME lagrangianCFDEMcomp-$CFDEM_VERSION-$WM_PROJECT_VERSION
+
 #- LMP Many2Many lib path and makefile
 setenv CFDEM_Many2ManyLIB_PATH $CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library
 setenv CFDEM_Many2ManyLIB_MAKEFILENAME $CFDEM_LIGGGHTS_MAKEFILE_NAME
@@ -98,6 +101,9 @@ alias cfdemUt 'cd $CFDEM_UT_DIR'
 #- shortcut to run path
 alias cfdemRun 'cd $CFDEM_PROJECT_USER_DIR/run'
 
+#- shortcut to user solver path
+alias cfdemUsrSol 'cd $CFDEM_PROJECT_USER_DIR/applications/solvers'
+
 #- shortcut to documentation path
 alias cfdemDoc 'cd $CFDEM_DOC_DIR'
 
@@ -140,6 +146,9 @@ alias cfdemCompCFDEMuti 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compil
 #- shortcut to test basic tutorials
 alias cfdemTestTUT 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/testTutorials.sh'
 
+#- shortcut to visualize the clock model data
+alias vizClock 'python $CFDEM_UT_DIR/vizClock/matPlot.py'
+
 #- shortcut to run liggghts in serial
 alias cfdemLiggghts '$CFDEM_LIGGGHTS_SRC_DIR/lmp_$CFDEM_LIGGGHTS_MAKEFILE_NAME'
 
@@ -157,20 +166,20 @@ if ( -d "$CFDEM_PROJECT_USER_DIR" ) then
     :
 else
     if ( -d "$CFDEM_PROJECT_DIR" ) then
-        cd $CFDEM_PROJECT_DIR/..
-        echo "make new dirs $CFDEM_PROJECT_DIR/../$LOGNAME-$CFDEM_VERSION-$WM_PROJECT_VERSION ? (y/n)"
+        echo "make new dirs $CFDEM_PROJECT_USER_DIR ? (y/n)"
 	set YN=$<
         if ( $YN == "y" ) then
-            mkdir $LOGNAME-$CFDEM_VERSION-$WM_PROJECT_VERSION
-            cd $LOGNAME-$CFDEM_VERSION-$WM_PROJECT_VERSION
+            mkdir -p $CFDEM_PROJECT_USER_DIR
+            cd $CFDEM_PROJECT_USER_DIR
             mkdir run
             mkdir -p log/logFilesCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
+            mkdir -p applications/solvers
         else
             echo "aborted by user."
             exit
         endif
     else
-        echo "error in CFDEMcoupling's bashrc."
+        echo "error in CFDEMcoupling's cshrc."
         exit
     endif
 endif

src/lagrangian/cfdemParticle/etc/functions.sh

@@ -71,6 +71,21 @@ compileLib()
 
     #- wclean and wmake
     #if [ $doClean != "noClean" ]; then
+        # check library to compile is compressible
+        str=$casePath
+        i=$((${#str}-4))
+        ending=${str:$i:4}
+        if [[ $ending == "Comp" ]]; then
+                echo "Compiling a compressible library - so doing an rmdepall of incomp library first."
+                echo "Please make sure to have the compressible libraries first in the library-list.txt!"
+                cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle
+                echo "changing to $PWD"
+                rmdepall 2>&1 | tee -a $logpath/$logfileName
+                cd $casePath
+                echo "changing to $PWD"
+            else
+                echo "Compiling a incompressible library."
+        fi
         rmdepall 2>&1 | tee -a $logpath/$logfileName
         wclean 2>&1 | tee -a $logpath/$logfileName
     #fi
@@ -139,6 +154,7 @@ compileLIGGGHTS()
     logpath="$1"
     logfileName="$2"
     headerText="$3"
+    clean="$4"
     #--------------------------------------------------------------------------------#
 
     #- clean up old log file
@@ -157,9 +173,14 @@ compileLIGGGHTS()
     echo 2>&1 | tee -a $logpath/$logfileName
 
     #- wclean and wmake
-    rm $CFDEM_LIGGGHTS_SRC_DIR/"lmp_"$CFDEM_LIGGGHTS_MAKEFILE_NAME
-    rm $CFDEM_LIGGGHTS_SRC_DIR/"lib"$CFDEM_LIGGGHTS_LIB_NAME".a"
-    make clean-all 2>&1 | tee -a $logpath/$logfileName
+    if [[ $clean == "false" ]]; then
+        echo "not cleaning LIGGGHTS"
+    else
+        rm $CFDEM_LIGGGHTS_SRC_DIR/"lmp_"$CFDEM_LIGGGHTS_MAKEFILE_NAME
+        rm $CFDEM_LIGGGHTS_SRC_DIR/"lib"$CFDEM_LIGGGHTS_LIB_NAME".a"
+        make clean-all 2>&1 | tee -a $logpath/$logfileName
+        echo "cleaning LIGGGHTS"
+    fi
     if [[ $WM_NCOMPPROCS == "" ]]; then
         echo "compiling LIGGGHTS on one CPU"
         make $CFDEM_LIGGGHTS_MAKEFILE_NAME 2>&1 | tee -a $logpath/$logfileName
@@ -189,7 +210,13 @@ compileLMPlib()
     rm $logpath/$logfileName
 
     #- change path
-    cd $libraryPath
+    if [ -d "$libraryPath" ]; then
+        cd $libraryPath
+    else
+        echo ""
+        echo "lib path $libraryPath does not exist - check settings in .../etc/bashrc."
+        read
+    fi
 
     #- header
     echo 2>&1 | tee -a $logpath/$logfileName
@@ -524,7 +551,7 @@ parCFDrun()
 
     #- run applictaion
     if [ $machineFileName == "none" ]; then
-        mpirun -np $nrProcs $solverName 2>&1 | tee -a $logpath/$logfileName
+        mpirun -np $nrProcs $solverName -parallel 2>&1 | tee -a $logpath/$logfileName
     else
         mpirun -machinefile $machineFileName -np $nrProcs $debugMode $solverName -parallel 2>&1 | tee -a $logpath/$logfileName
     fi
@@ -550,6 +577,7 @@ parCFDDEMrun()
     machineFileName="$7"
     debugMode="$8"
     reconstuctCase="$9"
+    cleanCase="$10"
     #--------------------------------------------------------------------------------#
 
     if [ $debugMode == "on" ]; then

src/lagrangian/cfdemParticle/etc/library-liggghts-list.txt

@@ -1,4 +1,4 @@
 #syntax: makefileName/dir
 #note: dir is not a path, just a keyword here
 ###############################################
-M2M/dir
+Many2Many/dir

src/lagrangian/cfdemParticle/etc/library-list.txt

@@ -1,11 +1,2 @@
 lagrangian/cfdemParticle/dir
-
-#====================================================='
-#- RADL
-fvOptions/dir
-cylPorousMedia/dir
-
-#====================================================='
-#- other
 finiteVolume/dir
-

src/lagrangian/cfdemParticle/etc/solver-list.txt

@@ -1,13 +1,4 @@
+cfdemSolverPisoMS/dir
 cfdemSolverPiso/dir
 cfdemSolverIB/dir
 cfdemSolverPisoScalar/dir
-cfdemSolverPimpleImEx/dir
-cfdemSolverIBInterLubrication/dir
-cfdemSolverIBScalar/dir
-cfdemSolverInterDyM/dir
-cfdemSolverInterDyMPC/dir
-cfdemSolverBubble/dir
-cfdemSolverPisoMS/dir
-cfdemSolverPimpleDyM_22x/dir
-cfdemSolverPimpleDyMMS_22x/dir
-cfdemSolverPimpleDyMScalar_22x/dir

src/lagrangian/cfdemParticle/etc/tutorial-list.txt

@@ -7,51 +7,15 @@
 #===================================================================#
 
 cfdemSolverPiso/settlingTestMPI/dir
+
 cfdemSolverPiso/ErgunTestMPI/dir
+
 cfdemSolverPiso/ErgunTestMPI_cgs/dir
+
 cfdemSolverPiso/ErgunTestMPI_restart/dir
+
 cfdemSolverIB/twoSpheresGlowinskiMPI/dir
+
 cfdemSolverPisoScalar/packedBedTemp/dir
 
-#===================================================================#
-# RADL
-cfdemSolverPimpleImEx/settlingTestMPI/dir
-cfdemSolverPimpleImEx/ErgunTestMPI/dir
-#cfdemSolverPimpleImEx/crossFlow/dir
-#cfdemSolverIB/periodicCase/dir
-#cfdemSolverIB/cfdemIBPeriodicCubicalBox_fullyPeriodic/dir
-#cfdemSolverIBInterLubrication/twoCoatedParticlesRelMotion_smallTest/dir
-#cfdemSolverIBScalar/cfdemIBPeriodicCubicalBoxScalar/dir
-
-#===================================================================#
-# NesteJacobs
-#Projects/Neste/cfdemSolverBubble/3pFBreactor/dir
-#Projects/Neste/cfdemSolverInterDyM/3pFBreactor/dir
-
-#===================================================================#
-# not in release:
-
-#cfdemSolverPiso/settlingTestBigParticleMPI/dir
-#cfdemSolverPiso/ErgunTestCG/dir
-#cfdemSolverPiso/ErgunTestM2M/dir
-#cfdemSolverPiso/HopperEmptying/dir
-
-cfdemSolverPimpleDyM/ErgunTestMPI/dir
-
-#cfdemSolverPisoMS/settlingTestMPI/dir
-#cfdemSolverPisoMS/ErgunTestMPI/dir
-
-#cfdemSolverInterDyM/twoPhaseSettlingTest/dir
-#cfdemSolverInterDyM/ErgunTestMPI/dir
-#cfdemSolverInterDyM/granularPiston/dir
-#cfdemSolverInterDyM/sugarNcoffee/dir
-
-#cfdemSolverBubble/ErgunTestMPI_pureLiquid/dir
-
-#- these examples are already designed for 2.3.x
-#cfdemSolverInterDyMPC/sugarNcoffee/dir
-#cfdemSolverInterDyMPC/granularPiston/dir
-#cfdemSolverInterDyMPC/meltingPot/dir
-
-
-
+cfdemSolverPiso/ErgunTestCG/dir

src/lagrangian/cfdemParticle/etc/utilities-list.txt

@@ -0,0 +1 @@
+cfdemPostproc/dir

src/lagrangian/cfdemParticle/subModels/IOModel/IOModel/IOModel.C

@@ -50,6 +50,15 @@ int IOModel::dumpDEMdata() const
     return -1;
 }
 
+bool IOModel::dumpNow() const
+{
+    //bool dmp(false);
+    //if (time_.value()+SMALL > time_.endTime().value()-time_.deltaT().value() || time_.outputTime())
+    //    dmp=true;
+
+    return time_.outputTime();
+}
+
 fileName IOModel::createTimeDir(fileName path) const
 {
     fileName timeDirPath(path/time_.timeName());

src/lagrangian/cfdemParticle/subModels/IOModel/IOModel/IOModel.H

@@ -113,6 +113,8 @@ public:
 
     virtual int dumpDEMdata() const;
 
+    bool dumpNow() const;
+
     fileName createTimeDir(fileName) const;
 
     fileName createLagrangianDir(fileName) const;

src/lagrangian/cfdemParticle/subModels/IOModel/basicIO/basicIO.C

@@ -84,7 +84,7 @@ basicIO::~basicIO()
 
 int basicIO::dumpDEMdata() const
 {
-    if (time_.outputTime())
+    if (dumpNow())
     {
         // make time directory
         if (parOutput_) lagPath_=buildFilePath(dirName_);

src/lagrangian/cfdemParticle/subModels/IOModel/sophIO/sophIO.C

@@ -79,7 +79,7 @@ int sophIO::dumpDEMdata() const
 {
     int npProcs(-1);
 
-    if (time_.outputTime())
+    if (dumpNow())
     {
         npProcs=basicIO::dumpDEMdata();
 

src/lagrangian/cfdemParticle/subModels/IOModel/trackIO/trackIO.C

@@ -77,7 +77,7 @@ int trackIO::dumpDEMdata() const
 {
     int npProcs(-1);
 
-    if (time_.outputTime())
+    if (dumpNow())
     {
         npProcs = sophIO::dumpDEMdata();