Merge pull request #30 from ParticulateFlow/release

Release 17.08

applications/solvers/cfdemSolverIB/Make/options

@@ -1,6 +1,7 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \

applications/solvers/cfdemSolverPiso/Make/options

@@ -1,6 +1,7 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \

applications/solvers/cfdemSolverPisoMS/Make/options

@@ -1,6 +1,7 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \

applications/solvers/cfdemSolverPisoScalar/Make/options

@@ -1,6 +1,7 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \

applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C

@@ -25,7 +25,7 @@ License
     along with CFDEMcoupling.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    cfdemSolverPiso
+    cfdemSolverPisoScalar
 
 Description
     Transient solver for incompressible flow.

applications/solvers/cfdemSolverRhoPimple/Make/options

@@ -4,6 +4,7 @@ PFLAGS+= -Dcompre
 
 EXE_INC = \
      $(PFLAGS) \
+    -I$(CFDEM_OFVERSION_DIR) \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \

applications/utilities/cfdemPostproc/Make/options

@@ -1,6 +1,7 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \

etc/bashrc

@@ -0,0 +1,259 @@
+#----------------------------------*-sh-*--------------------------------------
+# CFDEMcoupling
+#
+# Copyright 2009-2012 JKU Linz
+# Copyright 2012-2015 DCS Computing GmbH, Linz
+# Copyright 2015-     JKU Linz
+#------------------------------------------------------------------------------
+#
+# File
+#     etc/bashrc
+#
+# Description
+#     Startup file for CFDEMcoupling
+#     Sourced from ~/.profile or ~/.bashrc
+#
+# you can test the correctness using cfdemSystemTest.sh
+#------------------------------------------------------------------------------
+
+export CFDEM_PROJECT=CFDEM
+export CFDEM_VERSION=17.08
+
+################################################################################
+# USER EDITABLE PART: Changes made here may be lost with the next upgrade
+#
+# Please set to the appropriate path if the default is not correct.
+#
+# activate compatible OpenFOAM version
+. $HOME/OpenFOAM/OpenFOAM-4.x/etc/bashrc
+#
+# CFDEMcoupling
+export CFDEM_INST_DIR=$HOME/$CFDEM_PROJECT
+export CFDEM_PROJECT_DIR=$CFDEM_INST_DIR/CFDEMcoupling
+export CFDEM_PROJECT_USER_DIR=$CFDEM_INST_DIR/$USER-$WM_PROJECT_VERSION
+#
+# LIGGGHTS
+export CFDEM_LIGGGHTS_INST_DIR=$CFDEM_INST_DIR/LIGGGHTS
+#
+# LPP installation path
+export CFDEM_LPP_INST_DIR=$CFDEM_INST_DIR/LPP
+#
+# Path to additional libraries
+export CFDEM_ADD_LIBS_DIR=
+#
+# END OF (NORMAL) USER EDITABLE PART
+################################################################################
+
+#- The old dirs to be cleaned from the environment variables
+cfdemOldDirs="$CFDEM_LIGGGHTS_BIN_DIR $CFDEM_APP_DIR $CFDEM_LIB_DIR"
+
+export CFDEM_LIGGGHTS_SRC_DIR=$CFDEM_LIGGGHTS_INST_DIR/src
+export CFDEM_LIGGGHTS_BIN_DIR=$CFDEM_LIGGGHTS_INST_DIR/src-build
+
+#------------------------------------------------------------------------------
+# Source initialization functions
+#------------------------------------------------------------------------------
+. $WM_PROJECT_DIR/etc/config.sh/functions
+
+#------------------------------------------------------------------------------
+# Evaluate command-line parameters
+# these can be used to set/unset values
+#------------------------------------------------------------------------------
+_foamEval $@
+
+#------------------------------------------------------------------------------
+
+export CFDEM_SRC_DIR=$CFDEM_PROJECT_DIR/src
+export CFDEM_SOLVER_DIR=$CFDEM_PROJECT_DIR/applications/solvers
+export CFDEM_UT_DIR=$CFDEM_PROJECT_DIR/applications/utilities
+export CFDEM_DOC_DIR=$CFDEM_PROJECT_DIR/doc
+export CFDEM_TUT_DIR=$CFDEM_PROJECT_DIR/tutorials
+export CFDEM_LPP_DIR=$CFDEM_LPP_INST_DIR/src
+
+#------------------------------------------------------------------------------
+
+#- 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
+
+#- Check if additional libraries should be compiled together with solvers
+if [[ $CFDEM_ADD_LIBS_DIR == "" ]]; then
+    export CFDEM_ADD_LIBS_DIR=$CFDEM_PROJECT_DIR/etc
+else
+    echo "using CFDEM_ADD_LIBS_DIR=$CFDEM_ADD_LIBS_DIR defined by user."
+fi
+
+#------------------------------------------------------------------------------
+
+#- LMP Many2Many lib path and makefile
+export CFDEM_Many2ManyLIB_PATH=$CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library
+export CFDEM_Many2ManyLIB_MAKEFILENAME=fedora_fpic
+
+#------------------------------------------------------------------------------
+
+#- path to test harness
+export CFDEM_TEST_HARNESS_PATH=$CFDEM_PROJECT_USER_DIR/log/logFilesCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
+
+#- path to libraries
+export CFDEM_LIB_DIR=$CFDEM_PROJECT_DIR/platforms/$WM_OPTIONS/lib
+
+#- path to apps
+export CFDEM_APP_DIR=$CFDEM_PROJECT_DIR/platforms/$WM_OPTIONS/bin
+
+#- create directories
+mkdir -p $CFDEM_LIB_DIR
+mkdir -p $CFDEM_APP_DIR
+
+#- path to OF version flag file
+export CFDEM_OFVERSION_DIR=$CFDEM_PROJECT_DIR/etc/OFversion
+
+#------------------------------------------------------------------------------
+# Clean standard environment variables (PATH, LD_LIBRARY_PATH)
+#------------------------------------------------------------------------------
+foamClean=$WM_PROJECT_DIR/bin/foamCleanPath
+
+#- Clean PATH
+cleaned=`$foamClean "$PATH" "$cfdemOldDirs"` && PATH="$cleaned"
+
+#- Clean LD_LIBRARY_PATH
+cleaned=`$foamClean "$LD_LIBRARY_PATH" "$cfdemOldDirs"` \
+    && LD_LIBRARY_PATH="$cleaned"
+
+export PATH LD_LIBRARY_PATH
+
+#- add binary directories to $PATH
+_foamAddPath $CFDEM_APP_DIR:$CFDEM_LIGGGHTS_BIN_DIR
+_foamAddLib $CFDEM_LIB_DIR
+
+#------------------------------------------------------------------------------
+# Cleanup environment:
+#------------------------------------------------------------------------------
+unset cleaned foamClean cfdemOldDirs
+
+#------------------------------------------------------------------------------
+# Unload initialization functions:
+#------------------------------------------------------------------------------
+. $WM_PROJECT_DIR/etc/config.sh/functions
+
+#------------------------------------------------------------------------------
+# Settings for lpp postproc tool
+#------------------------------------------------------------------------------
+
+#- nr of procs for lpp tool
+export CFDEM_LPP_NPROCS=4
+
+#- nr of procs for lpp tool
+export CFDEM_LPP_CHUNKSIZE=1
+
+#- shortcut to run lpp
+alias lpp='python -i $CFDEM_LPP_DIR/lpp.py --cpunum $CFDEM_LPP_NPROCS --chunksize $CFDEM_LPP_CHUNKSIZE'
+
+#------------------------------------------------------------------------------
+# aliases for easy navigation (no changes necessary)
+#------------------------------------------------------------------------------
+
+#- shortcut to cfdem path
+alias cfdem='cd $CFDEM_PROJECT_DIR'
+
+#- shortcut to src path
+alias cfdemSrc='cd $CFDEM_SRC_DIR'
+
+#- shortcut to tutorial path
+alias cfdemTut='cd $CFDEM_TUT_DIR'
+
+#- shortcut to solver path
+alias cfdemSol='cd $CFDEM_SOLVER_DIR'
+
+#- shortcut to utilities path
+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'
+
+#- shortcut to open the doxygen with firefox
+alias cfdemDox='firefox $CFDEM_DOC_DIR/doxygen/html/index.html'
+
+#- shortcut to LIGGGHTS path
+alias cfdemLIG='cd $CFDEM_LIGGGHTS_SRC_DIR'
+
+#- shortcut to system test
+alias cfdemSysTest='bash $CFDEM_PROJECT_DIR/etc/cfdemSystemTest.sh'
+
+#- shortcut to pull LIGGGHTS
+alias cfdemPullLIG='bash $CFDEM_PROJECT_DIR/etc/pullLIGGGHTS.sh'
+
+#- shortcut to pull CFDEMcoupling
+alias cfdemPullCFDEMcoupling='bash $CFDEM_PROJECT_DIR/etc/pullCFDEMcoupling.sh'
+
+#- shortcut to clean CFDEM
+alias cfdemCleanCFDEM='bash $CFDEM_PROJECT_DIR/etc/cleanCFDEMcoupling.sh'
+
+#- shortcut to compile LIGGGHTS + sublibraries
+alias cfdemCompLIG='bash $CFDEM_PROJECT_DIR/etc/compileLIGGGHTS.sh'
+
+#- shortcut to compile CFDEMcoupling +LIGGGHTS
+alias cfdemCompCFDEMall='bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_all.sh'
+
+#- shortcut to compile CFDEMcoupling (src+solvers)
+alias cfdemCompCFDEM='bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling.sh'
+
+#- shortcut to compile CFDEMcoupling src
+alias cfdemCompCFDEMsrc='bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_src.sh'
+
+#- shortcut to compile CFDEMcoupling solvers
+alias cfdemCompCFDEMsol='bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_sol.sh'
+
+#- shortcut to compile CFDEMcoupling utilities
+alias cfdemCompCFDEMuti='bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_uti.sh'
+
+#- shortcut to test basic tutorials
+alias cfdemTestTUT='bash $CFDEM_PROJECT_DIR/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 {} \;'
+
+#- shortcut to run liggghts in serial
+cfdemLiggghts() { $CFDEM_LIGGGHTS_BIN_DIR/liggghts < $1; }
+export -f cfdemLiggghts
+
+#- shortcut to run liggghts in parallel
+cfdemLiggghtsPar() { mpirun -np $2 $CFDEM_LIGGGHTS_BIN_DIR/liggghts < $1; }
+export -f cfdemLiggghtsPar
+
+#- shortcut to open files including a pattern
+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 user directory exists
+if [ -d "$CFDEM_PROJECT_USER_DIR" ]; then
+    :
+else
+        echo "make new dirs $CFDEM_PROJECT_USER_DIR ? (y/n)"
+        read YN
+        if [ $YN == "y" ]; then
+            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
+        fi
+fi

etc/cfdemSystemTest.sh

@@ -1,12 +1,12 @@
 #!/bin/bash
 
 #===================================================================#
-# sytsem settings test routine for cfdem project 
+# sytsem settings test routine for cfdem project
 # Christoph Goniva - May. 2011, DCS Computing GmbH
 #===================================================================#
 
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 
 #- show gcc settings
 checkGPP="true"
@@ -33,18 +33,17 @@ checkDirComment "$CFDEM_SRC_DIR" '$CFDEM_SRC_DIR' "yes"
 checkDirComment "$CFDEM_SOLVER_DIR" '$CFDEM_SOLVER_DIR' "yes"
 checkDirComment "$CFDEM_TUT_DIR" '$CFDEM_TUT_DIR' "yes"
 checkDirComment "$CFDEM_LIGGGHTS_SRC_DIR" '$CFDEM_LIGGGHTS_SRC_DIR' "yes"
+checkEnvComment "$CFDEM_LIGGGHTS_BIN_DIR" '$CFDEM_LIGGGHTS_BIN_DIR' "yes"
 checkDirComment "$CFDEM_LPP_DIR" '$CFDEM_LPP_DIR' "yes"
 checkDirComment "$CFDEM_ADD_LIBS_DIR" '$CFDEM_ADD_LIBS_DIR' "yes"
-checkDirComment "$CFDEM_PIZZA_DIR" '$CFDEM_PIZZA_DIR' "no"
 checkDirComment "$CFDEM_TEST_HARNESS_PATH" '$CFDEM_TEST_HARNESS_PATH' "no"
-#checkDirComment "$C3PO_SRC_DIR" '$C3PO_SRC_DIR' "no"
 echo ""
 
 echo "library names"
-echo '$CFDEM_LIGGGHTS_LIB_NAME = '"$CFDEM_LIGGGHTS_LIB_NAME"
 echo '$CFDEM_LIB_NAME = '"$CFDEM_LIB_NAME"
-echo '$LD_LIBRARY_PATH  = '"$LD_LIBRARY_PATH"
-echo '$WM_NCOMPPROCS  = '"$WM_NCOMPPROCS"
+echo '$PATH = '"$PATH"
+echo '$LD_LIBRARY_PATH = '"$LD_LIBRARY_PATH"
+echo '$WM_NCOMPPROCS = '"$WM_NCOMPPROCS"
 
 echo "*******************"
 
@@ -78,6 +77,6 @@ if [ $checkAddOn == "true" ]
     if [ $(checkDir $filePath) == "true" ]; then
         source $filePath/etc/$packageName"SystemTest.sh"
     else
-        echo "$packageName does not exist." 
+        echo "$packageName does not exist."
     fi
 fi

etc/cleanCFDEMcoupling.sh

@@ -6,7 +6,7 @@
 #===================================================================#
 
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 
 #================================================================================#
 # clean src remove object files

etc/compileCFDEMcoupling.sh

@@ -6,26 +6,26 @@
 #===================================================================#
 
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 
 NOW="$(date +"%Y-%m-%d-%H:%M")"
 logDir="log"
 
 
-cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
+cd $CFDEM_PROJECT_DIR/etc
 mkdir -p $logDir
 
 #================================================================================#
 # compile src
 #================================================================================#
-bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_src.sh
+bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_src.sh
 
 #================================================================================#
 # compile solvers
 #================================================================================#
-bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_sol.sh
+bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_sol.sh
 
 #================================================================================#
 # compile utilities
 #================================================================================#
-bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_uti.sh
+bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_uti.sh

etc/compileCFDEMcoupling_all.sh

@@ -7,26 +7,26 @@
 #===================================================================#
 
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 
 NOW="$(date +"%Y-%m-%d-%H:%M")"
 logDir="log"
 
 
-cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
+cd $CFDEM_PROJECT_DIR/etc
 mkdir -p $logDir
 
 #================================================================================#
 # compile LIGGGHTS src
 #================================================================================#
-bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileLIGGGHTS.sh
+bash $CFDEM_PROJECT_DIR/etc/compileLIGGGHTS.sh
 
 #================================================================================#
 # compile LIGGGHTS libraries
 #================================================================================#
-bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileLIGGGHTS_lib.sh
+bash $CFDEM_PROJECT_DIR/etc/compileLIGGGHTS_lib.sh
 
 #================================================================================#
 # compile CFDEMcoupling
 #================================================================================#
-bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling.sh
+bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling.sh

etc/compileCFDEMcoupling_sol.sh

@@ -8,9 +8,9 @@
 whitelist="solver-list.txt"
 
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 logDir="log"
-cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
+cd $CFDEM_PROJECT_DIR/etc
 mkdir -p $logDir
 
 #- remove old success/fail logs

etc/compileCFDEMcoupling_src.sh

@@ -7,20 +7,20 @@
 #===================================================================#
  
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 
 NOW="$(date +"%Y-%m-%d-%H:%M")"
 logDir="log"
 
 
-cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
+cd $CFDEM_PROJECT_DIR/etc
 mkdir -p $logDir
 
 
 #================================================================================#
 # compile src
 #================================================================================#
-whitelist="$CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/library-list.txt"
+whitelist="$CFDEM_PROJECT_DIR/etc/library-list.txt"
 echo ""
 echo "Please provide the libraries to be compiled in the $CWD/$whitelist file."
 

etc/compileCFDEMcoupling_uti.sh

@@ -8,9 +8,9 @@
 whitelist="utilities-list.txt"
 
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 logDir="log"
-cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
+cd $CFDEM_PROJECT_DIR/etc
 mkdir -p $logDir
 
 CWD="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")"

etc/compileLIGGGHTS.sh

@@ -6,12 +6,12 @@
 #===================================================================
 
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 
 NOW="$(date +"%Y-%m-%d-%H:%M")"
 logDir="log"
 
-cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
+cd $CFDEM_PROJECT_DIR/etc
 mkdir -p $logDir
 
 #================================================================================#
@@ -38,4 +38,4 @@ compileLIGGGHTS $logpath $logfileName $headerText
 #================================================================================#
 # compile LIGGGHTS libraries
 #================================================================================#
-bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileLIGGGHTS_lib.sh
+bash $CFDEM_PROJECT_DIR/etc/compileLIGGGHTS_lib.sh

etc/compileLIGGGHTS_lib.sh

@@ -6,19 +6,19 @@
 #===================================================================#
  
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 
 NOW="$(date +"%Y-%m-%d-%H:%M")"
 logDir="log"
 
 
-cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
+cd $CFDEM_PROJECT_DIR/etc
 mkdir -p $logDir
 
 #================================================================================#
 # compile src
 #================================================================================#
-    whitelist="$CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/library-liggghts-list.txt"
+    whitelist="$CFDEM_PROJECT_DIR/etc/library-liggghts-list.txt"
     echo ""
     echo "Compiling sub-libraries of LIGGGHTS now..."
     echo "Please provide the libraries to be compiled in the $CWD/$whitelist file."

etc/cshrc

@@ -0,0 +1,289 @@
+#----------------------------------*-sh-*--------------------------------------
+# CFDEMcoupling
+#
+# Copyright 2009-2012 JKU Linz
+# Copyright 2012-2015 DCS Computing GmbH, Linz
+# Copyright 2015-     JKU Linz
+#------------------------------------------------------------------------------
+#
+# File
+#     etc/cshrc
+#
+# Description
+#     Startup file for CFDEMcoupling
+#     Sourced from ~/.cshrc
+#------------------------------------------------------------------------------
+
+setenv CFDEM_PROJECT CFDEM
+setenv CFDEM_VERSION 17.08
+
+################################################################################
+# USER EDITABLE PART: Changes made here may be lost with the next upgrade
+#
+# Please set to the appropriate path if the default is not correct.
+#
+# activate compatible OpenFOAM version
+. $HOME/OpenFOAM/OpenFOAM-4.x/etc/cshrc
+#
+# CFDEMcoupling
+setenv CFDEM_INST_DIR $HOME/$CFDEM_PROJECT
+setenv CFDEM_PROJECT_DIR $CFDEM_INST_DIR/CFDEMcoupling
+setenv CFDEM_PROJECT_USER_DIR $CFDEM_INST_DIR/$USER-$WM_PROJECT_VERSION
+#
+# LIGGGHTS
+setenv CFDEM_LIGGGHTS_INST_DIR $CFDEM_INST_DIR/LIGGGHTS
+#
+# LPP installation path
+setenv CFDEM_LPP_INST_DIR $CFDEM_INST_DIR/LPP
+#
+# Path to additional libraries
+setenv CFDEM_ADD_LIBS_DIR
+#
+# END OF (NORMAL) USER EDITABLE PART
+################################################################################
+
+#- The old dirs to be cleaned from the environment variables
+set cfdemOldDirs=
+if ( $?CFDEM_LIGGGHTS_BIN_DIR ) then
+   set cfdemOldDirs="$cfdemOldDirs $CFDEM_LIGGGHTS_BIN_DIR"
+endif
+if ( $?CFDEM_APP_DIR ) then
+   set cfdemOldDirs="$cfdemOldDirs $CFDEM_APP_DIR"
+endif
+if ( $?CFDEM_LIB_DIR ) then
+   set cfdemOldDirs="$cfdemOldDirs $CFDEM_LIB_DIR"
+endif
+
+setenv CFDEM_LIGGGHTS_SRC_DIR $CFDEM_LIGGGHTS_INST_DIR/src
+setenv CFDEM_LIGGGHTS_BIN_DIR $CFDEM_LIGGGHTS_INST_DIR/src-build
+
+#------------------------------------------------------------------------------
+# Source files, possibly with some verbosity
+#------------------------------------------------------------------------------
+alias _foamSource 'if ($?FOAM_VERBOSE && $?prompt) echo "Sourcing: \!*"; if (\!* != "") source \!*'
+
+#------------------------------------------------------------------------------
+# Evaluate command-line parameters
+# these can be used to set/unset values
+#------------------------------------------------------------------------------
+while ( $#argv > 0 )
+    switch ($argv[1])
+    case -*:
+        # stray option (not meant for us here) -> get out
+        break
+        breaksw
+    case *=:
+        # name=       -> unsetenv name
+        if ($?FOAM_VERBOSE && $?prompt) echo "unsetenv $argv[1]:s/=//"
+        eval "unsetenv $argv[1]:s/=//"
+        breaksw
+    case *=*:
+        # name=value  -> setenv name value
+        if ($?FOAM_VERBOSE && $?prompt) echo "setenv $argv[1]:s/=/ /"
+        eval "setenv $argv[1]:s/=/ /"
+        breaksw
+    default:
+        # filename: source it
+        if ( -f "$1" ) then
+            _foamSource "$1"
+        else
+            _foamSource `$WM_PROJECT_DIR/bin/foamEtcFile -silent "$1"`
+        endif
+        breaksw
+    endsw
+    shift
+end
+
+#------------------------------------------------------------------------------
+
+setenv CFDEM_SRC_DIR $CFDEM_PROJECT_DIR/src
+setenv CFDEM_SOLVER_DIR $CFDEM_PROJECT_DIR/applications/solvers
+setenv CFDEM_UT_DIR $CFDEM_PROJECT_DIR/applications/utilities
+setenv CFDEM_DOC_DIR $CFDEM_PROJECT_DIR/doc
+setenv CFDEM_TUT_DIR $CFDEM_PROJECT_DIR/tutorials
+setenv CFDEM_LPP_DIR $CFDEM_LPP_INST_DIR/src
+
+#------------------------------------------------------------------------------
+
+#- 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
+
+#- Check if additional libraries should be compiled together with solvers
+if ( ! ($?CFDEM_ADD_LIBS_DIR) ) then
+    setenv CFDEM_ADD_LIBS_DIR $CFDEM_PROJECT_DIR/etc
+else
+    echo "using CFDEM_ADD_LIBS_DIR=$CFDEM_ADD_LIBS_DIR defined by user."
+endif
+
+#------------------------------------------------------------------------------
+
+#- LMP Many2Many lib path and makefile
+setenv CFDEM_Many2ManyLIB_PATH $CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library
+setenv CFDEM_Many2ManyLIB_MAKEFILENAME fedora_fpic
+
+#------------------------------------------------------------------------------
+
+#- path to test harness
+setenv CFDEM_TEST_HARNESS_PATH $CFDEM_PROJECT_USER_DIR/log/logFilesCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
+
+#- path to libraries
+setenv CFDEM_LIB_DIR $CFDEM_PROJECT_DIR/platforms/$WM_OPTIONS/lib
+
+#- path to apps
+setenv CFDEM_APP_DIR $CFDEM_PROJECT_DIR/platforms/$WM_OPTIONS/bin
+
+#- create directories
+mkdir -p $CFDEM_LIB_DIR
+mkdir -p $CFDEM_APP_DIR
+
+#- path to OF version flag file
+setenv CFDEM_OFVERSION_DIR $CFDEM_PROJECT_DIR/etc/OFversion
+
+#------------------------------------------------------------------------------
+# Clean standard environment variables (PATH, LD_LIBRARY_PATH)
+#------------------------------------------------------------------------------
+set foamClean=$WM_PROJECT_DIR/bin/foamCleanPath
+
+#- prevent local variables from shadowing setenv variables
+unset PATH LD_LIBRARY_PATH
+
+if (! $?LD_LIBRARY_PATH ) setenv LD_LIBRARY_PATH ''
+
+#- Clean PATH
+set cleaned=`$foamClean "$PATH" "$cfdemOldDirs"`
+if ( $status == 0 ) setenv PATH $cleaned
+
+#- Clean LD_LIBRARY_PATH
+set cleaned=`$foamClean "$LD_LIBRARY_PATH" "$cfdemOldDirs"`
+if ( $status == 0 ) setenv LD_LIBRARY_PATH $cleaned
+
+#- add binary directories to $PATH
+_foamAddPath $CFDEM_APP_DIR:$CFDEM_LIGGGHTS_BIN_DIR
+_foamAddLib $CFDEM_LIB_DIR
+
+#------------------------------------------------------------------------------
+# Cleanup environment:
+#------------------------------------------------------------------------------
+unset cleaned foamClean cfdemOldDirs
+unalias _foamSource
+
+#------------------------------------------------------------------------------
+#- settings for lpp postproc tool
+#------------------------------------------------------------------------------
+
+#- nr of procs for lpp tool
+setenv CFDEM_LPP_NPROCS 4
+
+#- nr of procs for lpp tool
+setenv CFDEM_LPP_CHUNKSIZE 1
+
+#- shortcut to run lpp
+alias lpp 'python -i $CFDEM_LPP_DIR/lpp.py --cpunum $CFDEM_LPP_NPROCS --chunksize $CFDEM_LPP_CHUNKSIZE \!:1'
+
+#------------------------------------------------------------------------------
+# aliases for easy navigation (no changes necessary)
+#------------------------------------------------------------------------------
+
+#- shortcut to cfdem path
+alias cfdem 'cd $CFDEM_PROJECT_DIR'
+
+#- shortcut to src path
+alias cfdemSrc 'cd $CFDEM_SRC_DIR'
+
+#- shortcut to tutorial path
+alias cfdemTut 'cd $CFDEM_TUT_DIR'
+
+#- shortcut to solver path
+alias cfdemSol 'cd $CFDEM_SOLVER_DIR'
+
+#- shortcut to utilities path
+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'
+
+#- shortcut to open the doxygen with firefox
+alias cfdemDox 'firefox $CFDEM_DOC_DIR/doxygen/html/index.html'
+
+#- shortcut to LIGGGHTS path
+alias cfdemLIG 'cd $CFDEM_LIGGGHTS_SRC_DIR'
+
+#- shortcut to system test
+alias cfdemSysTest 'bash $CFDEM_PROJECT_DIR/etc/cfdemSystemTest.sh'
+
+#- shortcut to pull LIGGGHTS
+alias cfdemPullLIG 'bash $CFDEM_PROJECT_DIR/etc/pullLIGGGHTS.sh'
+
+#- shortcut to pull CFDEMcoupling
+alias cfdemPullCFDEMcoupling 'bash $CFDEM_PROJECT_DIR/etc/pullCFDEMcoupling.sh'
+
+#- shortcut to clean CFDEM
+alias cfdemCleanCFDEM 'bash $CFDEM_PROJECT_DIR/etc/cleanCFDEMcoupling.sh'
+
+#- shortcut to compile LIGGGHTS + sublibraries
+alias cfdemCompLIG 'bash $CFDEM_PROJECT_DIR/etc/compileLIGGGHTS.sh'
+
+#- shortcut to compile CFDEMcoupling +LIGGGHTS
+alias cfdemCompCFDEMall 'bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_all.sh'
+
+#- shortcut to compile CFDEMcoupling (src+solvers)
+alias cfdemCompCFDEM 'bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling.sh'
+
+#- shortcut to compile CFDEMcoupling src
+alias cfdemCompCFDEMsrc 'bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_src.sh'
+
+#- shortcut to compile CFDEMcoupling solvers
+alias cfdemCompCFDEMsol 'bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_sol.sh'
+
+#- shortcut to compile CFDEMcoupling utilities
+alias cfdemCompCFDEMuti 'bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_uti.sh'
+
+#- shortcut to test basic tutorials
+alias cfdemTestTUT 'bash $CFDEM_PROJECT_DIR/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_BIN_DIR/liggghts'
+
+#- shortcut to run liggghts in parallel (no fcts in csh...)
+alias cfdemLiggghtsPar 'echo "mpirun -np xx -machinefile mynodes $CFDEM_LIGGGHTS_BIN_DIR/liggghts < in.liggghts_init"'
+
+#- shortcut to run liggghts in parallel
+# unfortunately no functions available in csh
+
+#- shortcut to run lpp
+alias lpp 'python -i $CFDEM_LPP_DIR/lpp.py \!:1'
+
+#- check if the user directory exists
+if ( -d "$CFDEM_PROJECT_USER_DIR" ) then
+    :
+else
+    if ( -d "$CFDEM_PROJECT_DIR" ) then
+        echo "make new dirs $CFDEM_PROJECT_USER_DIR ? (y/n)"
+        set YN=$<
+        if ( $YN == "y" ) then
+            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 cshrc."
+        exit
+    endif
+endif

etc/functions.sh

@@ -249,7 +249,7 @@ cleanCFDEM()
 
     #**********************************************
     #cleaning libraries
-    whitelist="$CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/library-list.txt"
+    whitelist="$CFDEM_PROJECT_DIR/etc/library-list.txt"
     echo ""
     echo "Please provide the libraries to be cleaned in the $CWD/$whitelist file."
 
@@ -303,7 +303,7 @@ cleanCFDEM()
 
     #**********************************************
     #cleaning solvers
-    whitelist="$CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/solver-list.txt"
+    whitelist="$CFDEM_PROJECT_DIR/etc/solver-list.txt"
     echo ""
     echo "Please provide the solvers to be cleaned in the $CWD/$whitelist file."
 

etc/pullLIGGGHTS.sh

@@ -6,12 +6,12 @@
 #===================================================================
 
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 
 NOW="$(date +"%Y-%m-%d-%H:%M")"
 logDir="log"
 
-cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
+cd $CFDEM_PROJECT_DIR/etc
 mkdir -p $logDir
 
 #--------------------------------------------------------------------------------#

src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.C

@@ -195,6 +195,7 @@ template<class Type>
 void uniformFixedValueTubeFvPatchField<Type>::write(Ostream& os) const
 {
     fvPatchField<Type>::write(os);
+    this->writeEntry("value", os);
     uniformValue_->writeData(os);
     os.writeKeyword("tubeLength") << tubeLength_ << token::END_STATEMENT << nl;
     os.writeKeyword("tubeDiameter") << tubeDiameter_ << token::END_STATEMENT << nl;

src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueVoidfraction/uniformFixedValueVoidfractionFvPatchField.C

@@ -145,6 +145,7 @@ template<class Type>
 void uniformFixedValueVoidfractionFvPatchField<Type>::write(Ostream& os) const
 {
     fvPatchField<Type>::write(os);
+    this->writeEntry("value", os);
     uniformValue_->writeData(os);
 }
 

src/lagrangian/cfdemParticle/Make/files

@@ -111,6 +111,7 @@ $(voidFractionModels)/dividedVoidFractionMS/dividedVoidFractionMS.C
 $(voidFractionModels)/bigParticleVoidFraction/bigParticleVoidFraction.C
 $(voidFractionModels)/GaussVoidFraction/GaussVoidFraction.C
 $(voidFractionModels)/IBVoidFraction/IBVoidFraction.C
+$(voidFractionModels)/trilinearVoidFraction/trilinearVoidFraction.C
 
 $(locateModels)/locateModel/locateModel.C
 $(locateModels)/locateModel/newLocateModel.C

src/lagrangian/cfdemParticle/Make/options

@@ -9,6 +9,7 @@ include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 EXE_INC = \
      $(PFLAGS) \
      $(PINC) \
+    -I$(CFDEM_OFVERSION_DIR) \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \

src/lagrangian/cfdemParticle/cfdTools/debugInfo.H

@@ -1,7 +1,7 @@
-
+{
   //========================================================================//
-  scalar countCell=0;              // number of cells touched by particles
-  int points=0;                    // number of particles and sub-points
+  label countCell=0;               // number of cells touched by particles
+  label points=0;                  // number of particles and sub-points
   scalar totalParticleWeights=0;   // total weight of all particles and sub-points
   vector totalForce_array(0,0,0);  // total force on particles based on particle array
   vector totalForce_field(0,0,0);  // forceField of forceM(), used to calc Ksl
@@ -93,3 +93,4 @@
   Info <<"meanR_array = "<< meanR_array << endl;
   Info <<"=============================================================================" << endl;
   Info << endl;
+}

src/lagrangian/cfdemParticle/cfdTools/mathExtra.H

@@ -27,7 +27,7 @@ License
 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
 ------------------------------------------------------------------------- */
@@ -36,15 +36,37 @@ Description
 #define CFDEM_MATH_EXTRA_H
 
 #include <vector>
-//#include "math.h"
 #include <stdio.h>
 #include <string.h>
 #include <error.h>
 #include <ctype.h>
+#include "mathematicalConstants.H"
 
 #define TOLERANCE_ORTHO 1e-10
 
-namespace MathExtra 
+namespace Foam
+{
+namespace constant
+{
+namespace mathematical
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    const scalar piByFour(0.25*pi);
+    const scalar fourPiByThree(4.0*pi/3.0);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace mathematical
+} // End namespace constant
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+
+namespace MathExtra
 {
 
 //  inline void     outerProduct(double *vec1, double *vec2, double **m);
@@ -81,13 +103,13 @@ inline bool spheroidGeometry(double radius, double aspectRatio, //inputs
     //INPUT
     //  radius      ...volume-equivalent radius of the spheroid
     //  aspectRatio ...major/minor aspect ratio
-    
+
     //OUTPUT
-    //  ai  ... 
-    //  bi  ... 
-    //  ei  ... 
-    //  Le  ... 
-  
+    //  ai  ...
+    //  bi  ...
+    //  ei  ...
+    //  Le  ...
+
     if(radius<=0.0)     //avoid troubles in case radius is 0 or negative
         return false;
 
@@ -119,7 +141,7 @@ inline bool spheroidGeometry2(double radius,  double aspectRatio, //inputs
     //INPUT
     //  radius      ...volume-equivalent radius of the spheroid
     //  aspectRatio ...major/minor aspect ratio
-    
+
     //OUTPUT
     //  XAe  ...Eccentricity dependet parameter
     //  YAe  ...Eccentricity dependet parameter
@@ -212,7 +234,7 @@ inline bool permutationDotDyadic(
     //Generate permutation tensor
     double permutationT[3][3][3];
     permutationTensor(permutationT);
-    
+
     //Step 1: compute dot prodcut of permutation tensor
     double permutationDotProd[3][3];
     zeroize33(permutationDotProd);
@@ -227,7 +249,7 @@ inline bool permutationDotDyadic(
         for(int iY=0; iY<3; iY++)
             for(int iZ=0; iZ<3; iZ++)
                 tensor[iX][iY][iZ] = permutationDotProd[iX][iY]
-                                   * vector[iZ]; 
+                                   * vector[iZ];
     return true;
 
 }
@@ -240,7 +262,7 @@ inline bool doubleDotTensor333Tensor33(double tensor333[3][3][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++)

src/lagrangian/cfdemParticle/cfdTools/setupProbeModel.H

@@ -1,3 +1,2 @@
- //set probeModel parameters for this force model
- particleCloud_.probeM().setOutputFile();
- particleCloud_.probeM().setCounter();
+//set probeModel parameters for this force model
+if (probeIt_) { particleCloud_.probeM().setOutputFile(typeName+".logDat"); }

src/lagrangian/cfdemParticle/cfdTools/versionInfo.H

@@ -34,8 +34,8 @@ Description
 #ifndef versionInfo_H
 #define versionInfo_H
 
-word CFDEMversion="PFM 17.02";
-word compatibleLIGGGHTSversion="PFM 17.02";
+word CFDEMversion="PFM 17.08";
+word compatibleLIGGGHTSversion="PFM 17.08";
 word OFversion="4.x";
 
 Info << "\nCFDEMcoupling version: " << CFDEMversion <<  endl;

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C

@@ -81,6 +81,7 @@ cfdemCloud::cfdemCloud
     solveFlow_(true),
     verbose_(false),
     ignore_(false),
+    allowCFDsubTimestep_(true),
     limitDEMForces_(false),
     modelType_(couplingProperties_.lookup("modelType")),
     positions_(NULL),
@@ -126,6 +127,7 @@ cfdemCloud::cfdemCloud
         mesh,
         dimensionedScalar("zero", dimensionSet(0,0,-1,0,0), 0)  // 1/s
     ),
+    checkPeriodicCells_(false),
     turbulence_
     (
         mesh.lookupObject<turbulenceModel>
@@ -196,7 +198,7 @@ cfdemCloud::cfdemCloud
         clockModel::New
         (
             couplingProperties_,
-            *this
+            mesh.time()
         )
     ),
     smoothingModel_
@@ -225,6 +227,14 @@ cfdemCloud::cfdemCloud
         solveFlow_=Switch(couplingProperties_.lookup("solveFlow"));
     if (couplingProperties_.found("imExSplitFactor"))
         imExSplitFactor_ = readScalar(couplingProperties_.lookup("imExSplitFactor"));
+
+    if(imExSplitFactor_ > 1.0)
+        FatalError << "You have set imExSplitFactor > 1 in your couplingProperties. Must be <= 1."
+                   << abort(FatalError);
+    if(imExSplitFactor_ < 0.0)
+        FatalError << "You have set imExSplitFactor < 0 in your couplingProperties. Must be >= 0."
+                   << abort(FatalError);
+
     if (couplingProperties_.found("treatVoidCellsAsExplicitForce"))
         treatVoidCellsAsExplicitForce_ = readBool(couplingProperties_.lookup("treatVoidCellsAsExplicitForce"));
     if (couplingProperties_.found("verbose")) verbose_=true;
@@ -242,17 +252,6 @@ cfdemCloud::cfdemCloud
     else
         Info << "ignoring ddt(voidfraction)" << endl;
 
-    forceModel_ = new autoPtr<forceModel>[nrForceModels()];
-    for (int i=0;i<nrForceModels();i++)
-    {
-        forceModel_[i] = forceModel::New
-        (
-            couplingProperties_,
-            *this,
-            forceModels_[i]
-        );
-    }
-
     momCoupleModel_ = new autoPtr<momCoupleModel>[momCoupleModels_.size()];
     for (int i=0;i<momCoupleModels_.size();i++)
     {
@@ -264,6 +263,17 @@ cfdemCloud::cfdemCloud
         );
     }
 
+    forceModel_ = new autoPtr<forceModel>[nrForceModels()];
+    for (int i=0;i<nrForceModels();i++)
+    {
+        forceModel_[i] = forceModel::New
+        (
+            couplingProperties_,
+            *this,
+            forceModels_[i]
+        );
+    }
+
     // run liggghts commands from cfdem
     liggghtsCommand_ = new autoPtr<liggghtsCommandModel>[liggghtsCommandModelList_.size()];
     for (int i=0;i<liggghtsCommandModelList_.size();i++)
@@ -289,7 +299,44 @@ cfdemCloud::cfdemCloud
     }
 
     dataExchangeM().setCG();
-    if (!cgOK_ && cg_ > 1) FatalError<< "at least one of your models is not fit for cg !!!"<< abort(FatalError);
+    Switch cgWarnOnly_(couplingProperties_.lookupOrDefault<Switch>("cgWarnOnly", true));
+    if (!cgOK_ && cg_ > 1)
+    {
+        if (cgWarnOnly_)
+            Warning << "at least one of your models is not fit for cg !!!" << endl;
+        else
+            FatalError << "at least one of your models is not fit for cg !!!" << abort(FatalError);
+    }
+
+    // check if simulation is a fully periodic box
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    int nPatchesCyclic(0);
+    int nPatchesNonCyclic(0);
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+        if (isA<cyclicPolyPatch>(pp) || isA<cyclicAMIPolyPatch>(pp))
+            ++nPatchesCyclic;
+        else if (!isA<processorPolyPatch>(pp))
+            ++nPatchesNonCyclic;
+    }
+
+    if (nPatchesNonCyclic == 0)
+    {
+        checkPeriodicCells_ = true;
+    }
+
+    //hard set checkperiodic cells if wished
+    if(this->couplingProperties().found("checkPeriodicCells"))
+    {
+        checkPeriodicCells_ = couplingProperties().lookupOrDefault<Switch>("checkPeriodicCells", checkPeriodicCells_);
+    }
+
+    if (nPatchesCyclic > 0 && nPatchesNonCyclic > 0)
+    {
+        if (verbose_) Info << "nPatchesNonCyclic=" << nPatchesNonCyclic << ", nPatchesCyclic=" << nPatchesCyclic << endl;
+        Warning << "Periodic handing is disabled because the domain is not fully periodic!\n" << endl;
+    }
 }
 
 // * * * * * * * * * * * * * * * * Destructors  * * * * * * * * * * * * * * //
@@ -412,12 +459,15 @@ void cfdemCloud::setVectorAverages()
     );
     if(verbose_) Info << "setVectorAverage done." << endl;
 }
+
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
+
 void cfdemCloud::checkCG(bool ok)
 {
     if(!cgOK_) return;
     if(!ok) cgOK_ = ok;
 }
+
 void cfdemCloud::setPos(double**& pos)
 {
     for(int index = 0;index <  numberOfParticles(); ++index)
@@ -427,40 +477,32 @@ void cfdemCloud::setPos(double**& pos)
         }
     }
 }
+
 // * * * * * * * * * * * * * * * ACCESS  * * * * * * * * * * * * * //
 
-label cfdemCloud::particleCell(int index)
+label cfdemCloud::particleCell(int index) const
 {
-    label cellI = cellIDs()[index][0];
-    return cellI;
+    return cellIDs()[index][0];
 }
 
-vector cfdemCloud::position(int index)
+vector cfdemCloud::position(int index) const
 {
-    vector pos;
-    for(int i=0;i<3;i++) pos[i] = positions()[index][i];
-    return pos;
+    return vector(positions()[index][0],positions()[index][1],positions()[index][2]);
 }
 
-vector cfdemCloud::velocity(int index)
+vector cfdemCloud::velocity(int index) const
 {
-    vector vel;
-    for(int i=0;i<3;i++) vel[i] = velocities()[index][i];
-    return vel;
+    return vector(velocities()[index][0],velocities()[index][1],velocities()[index][2]);
 }
 
-vector cfdemCloud::expForce(int index)
+vector cfdemCloud::expForce(int index) const
 {
-    vector force;
-    for(int i=0;i<3;i++) force[i] = DEMForces()[index][i];
-    return force;
+    return vector(DEMForces()[index][0],DEMForces()[index][1],DEMForces()[index][2]);
 }
 
-vector cfdemCloud::fluidVel(int index)
+vector cfdemCloud::fluidVel(int index) const
 {
-    vector vel;
-    for(int i=0;i<3;i++) vel[i] = fluidVels()[index][i];
-    return vel;
+    return vector(fluidVels()[index][0],fluidVels()[index][1],fluidVels()[index][2]);
 }
 
 const forceModel& cfdemCloud::forceM(int i)
@@ -468,43 +510,31 @@ const forceModel& cfdemCloud::forceM(int i)
     return forceModel_[i];
 }
 
-int cfdemCloud::nrForceModels()
+label cfdemCloud::nrForceModels() const
 {
     return forceModels_.size();
 }
 
-int cfdemCloud::nrMomCoupleModels()
+label cfdemCloud::nrMomCoupleModels() const
 {
     return momCoupleModels_.size();
 }
 
-scalar cfdemCloud::voidfraction(int index)
+scalar cfdemCloud::voidfraction(int index) const
 {
     return voidfractions()[index][0];
 }
 
-label cfdemCloud::liggghtsCommandModelIndex(word name)
+label cfdemCloud::liggghtsCommandModelIndex(word name) const
 {
-    int index=-1;
     forAll(liggghtsCommandModelList_,i)
     {
         if(liggghtsCommand()[i]().name() == name)
         {
-            index = i;
-            break;
+            return i;
         }
     }
-    return index;
-}
-
-std::vector< std::vector<double*> >* cfdemCloud::getVprobe()
-{
- return probeModel_->getVprobe();
-}
-
-std::vector< std::vector<double> >* cfdemCloud::getSprobe()
-{
- return probeModel_->getSprobe();
+    return -1;
 }
 
 // * * * * * * * * * * * * * * * WRITE  * * * * * * * * * * * * * //
@@ -590,6 +620,10 @@ bool cfdemCloud::evolve
         //      IMPLICIT FORCE CONTRIBUTION AND SOLVER USE EXACTLY THE SAME AVERAGED
         //      QUANTITIES AT THE GRID!
         Info << "\n timeStepFraction() = " << dataExchangeM().timeStepFraction() << endl;
+        if(dataExchangeM().timeStepFraction() > 1.0000001)
+        {
+            FatalError << "cfdemCloud::dataExchangeM().timeStepFraction()>1: Do not do this, since dangerous. This might be due to the fact that you used a adjustable CFD time step. Please use a fixed CFD time step." << abort(FatalError);
+        }
         clockM().start(24,"interpolateEulerFields");
 
         // update voidFractionField
@@ -675,29 +709,6 @@ bool cfdemCloud::reAllocArrays()
     return false;
 }
 
-bool cfdemCloud::reAllocArrays(int nP, bool forceRealloc)
-{
-    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_,-1,voidFractionM().maxCellsPerParticle(),nP);
-        dataExchangeM().allocateArray(particleWeights_,0.,voidFractionM().maxCellsPerParticle(),nP);
-        dataExchangeM().allocateArray(particleVolumes_,0.,voidFractionM().maxCellsPerParticle(),nP);
-        arraysReallocated_ = true;
-        return true;
-    }
-    return false;
-}
-
 tmp<fvVectorMatrix> cfdemCloud::divVoidfractionTau(volVectorField& U,volScalarField& voidfraction) const
 {
     return

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H

@@ -49,7 +49,7 @@ SourceFiles
 #include "fvCFD.H"
 #include "IFstream.H"
 
-#include <turbulenceModel.H>
+#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -92,6 +92,8 @@ protected:
 
     bool ignore_;
 
+    bool allowCFDsubTimestep_;
+
     bool limitDEMForces_;
 
     scalar maxDEMForce_;
@@ -160,6 +162,8 @@ protected:
 
     mutable volScalarField ddtVoidfraction_;
 
+    mutable Switch checkPeriodicCells_;
+
     const turbulenceModel& turbulence_;
 
     autoPtr<forceModel>* forceModel_;
@@ -224,45 +228,49 @@ public:
 // public Member Functions
 
     // Access
+       bool allowCFDsubTimestep() { return allowCFDsubTimestep_; }
+
+       void setAllowCFDsubTimestep(bool b) { allowCFDsubTimestep_ = b; }
+
        void checkCG(bool);
 
        void setPos(double **&);
 
-       word modelType(){ return modelType_; }
+       const word& modelType() const { return modelType_; }
 
-       label particleCell(int);
+       label particleCell(int) const;
 
-       vector position(int);
+       vector position(int) const;
 
-       vector velocity(int);
+       vector velocity(int) const;
 
-       vector expForce(int);
+       vector expForce(int) const;
 
-       vector fluidVel(int);
+       vector fluidVel(int) const;
 
        virtual const forceModel& forceM(int);
 
-       virtual int nrForceModels();
+       virtual label nrForceModels() const;
 
-       virtual int nrMomCoupleModels();
+       virtual label nrMomCoupleModels() const;
 
-       scalar voidfraction(int);
+       scalar voidfraction(int) const;
 
-       label liggghtsCommandModelIndex(word);
+       label liggghtsCommandModelIndex(word) const;
 
-       inline void setCG(double) const;
+       inline void setCG(double);
 
-       inline const scalar& cg() const;
+       inline scalar cg() const;
 
-       inline const bool& impDEMdrag() const;
+       inline bool impDEMdrag() const;
 
-       inline const bool& impDEMdragAcc() const;
+       inline bool impDEMdragAcc() const;
 
-       inline const scalar& imExSplitFactor() const;
+       inline scalar imExSplitFactor() const;
 
-       inline const bool& treatVoidCellsAsExplicitForce() const;
+       inline bool treatVoidCellsAsExplicitForce() const;
 
-       inline const bool& ignore() const;
+       inline bool ignore() const;
 
        inline const fvMesh& mesh() const;
 
@@ -300,13 +308,13 @@ public:
 
        inline double ** particleWeights() const;
 
-       virtual inline label body(int);
+       virtual inline label body(int) const;
 
-       virtual inline double particleVolume(int);
+       virtual inline double particleVolume(int) const;
 
-       inline scalar radius(int);
+       inline scalar radius(int) const;
 
-       virtual inline double d(int);
+       virtual inline double d(int) const;
 
        inline scalar d32(bool recalc=true);
        virtual inline double dMin() {return -1;}
@@ -322,11 +330,11 @@ public:
        //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);}
+       virtual vector omega(int) const {return 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);}
+       virtual vector ex(int) const {return vector(0,0,0);}
 
        //Detector if SRF module is enable or not
        virtual inline bool SRFOn(){return false;}
@@ -339,7 +347,7 @@ public:
 
        inline bool arraysReallocated() const;
 
-       inline const wordList& forceModels();
+       inline const wordList& forceModels() const;
 
        inline const voidFractionModel& voidFractionM() const;
 
@@ -347,11 +355,11 @@ public:
 
        inline const momCoupleModel& momCoupleM(int) const;
 
-       inline const dataExchangeModel& dataExchangeM() const;
+       inline dataExchangeModel& dataExchangeM();
 
        inline const IOModel& IOM() const;
 
-       inline const probeModel& probeM() const;
+       inline probeModel& probeM();
 
        inline const averagingModel& averagingM() const;
 
@@ -376,13 +384,10 @@ public:
 
         virtual bool reAllocArrays();
 
-        virtual bool reAllocArrays(int nP, bool forceRealloc);  //force number of particles during reallocation
-
-
     // IO
-        void writeScalarFieldToTerminal(double**&);
+        void writeScalarFieldToTerminal(double**&) const;
 
-        void writeVectorFieldToTerminal(double**&);
+        void writeVectorFieldToTerminal(double**&) const;
 
     // functions
         tmp<fvVectorMatrix> divVoidfractionTau(volVectorField& ,volScalarField&) const;
@@ -397,11 +402,9 @@ public:
 
         void resetArray(double**&,int,int,double resetVal=0.);
 
-        std::vector< std::vector<double*> >* getVprobe();
-
-        std::vector< std::vector<double> >* getSprobe();
-
         void otherForces(volVectorField&);
+
+        bool checkPeriodicCells() { return checkPeriodicCells_; }
 };
 
 

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H

@@ -44,38 +44,38 @@ namespace Foam
 {
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-inline void cfdemCloud::setCG(double cg) const
+inline void cfdemCloud::setCG(double cg)
 {
     cg_ = cg;
     Info << "cg is set to: " << cg_ << endl;
 }
 
-inline const bool& cfdemCloud::impDEMdrag() const
+inline bool cfdemCloud::impDEMdrag() const
 {
     return impDEMdrag_;
 }
 
-inline const bool& cfdemCloud::impDEMdragAcc() const
+inline bool cfdemCloud::impDEMdragAcc() const
 {
     return impDEMdragAcc_;
 }
 
-inline const scalar& cfdemCloud::imExSplitFactor() const
+inline scalar cfdemCloud::imExSplitFactor() const
 {
     return imExSplitFactor_;
 }
 
-inline const bool& cfdemCloud::treatVoidCellsAsExplicitForce() const
+inline bool cfdemCloud::treatVoidCellsAsExplicitForce() const
 {
     return treatVoidCellsAsExplicitForce_;
 }
 
-inline const scalar& cfdemCloud::cg() const
+inline scalar cfdemCloud::cg() const
 {
     return cg_;
 }
 
-inline const bool& cfdemCloud::ignore() const
+inline bool cfdemCloud::ignore() const
 {
     return ignore_;
 }
@@ -177,24 +177,24 @@ inline double ** cfdemCloud::particleWeights() const
     return particleWeights_;
 }
 
-inline label cfdemCloud::body(int index)
+inline label cfdemCloud::body(int index) const
 {
     return index;
 }
 
-inline double cfdemCloud::particleVolume(int index)
+inline double cfdemCloud::particleVolume(int index) const
 {
     return particleV_[index][0];
 }
 
-inline scalar cfdemCloud::radius(int index)
+inline scalar cfdemCloud::radius(int index) const
 {
     return radii_[index][0];
 }
 
-inline double cfdemCloud::d(int index)
+inline double cfdemCloud::d(int index) const
 {
-    return 2*radii_[index][0];
+    return 2.*radii_[index][0];
 }
 
 inline double cfdemCloud::d32(bool recalc)
@@ -237,7 +237,7 @@ inline bool cfdemCloud::arraysReallocated() const
     return arraysReallocated_;
 }
 
-inline const wordList& cfdemCloud::forceModels()
+inline const wordList& cfdemCloud::forceModels() const
 {
     return forceModels_;
 }
@@ -252,9 +252,9 @@ inline const momCoupleModel& cfdemCloud::momCoupleM(int i) const
     return momCoupleModel_[i];
 }
 
-inline const dataExchangeModel& cfdemCloud::dataExchangeM() const
+inline dataExchangeModel& cfdemCloud::dataExchangeM()
 {
-    return dataExchangeModel_;
+    return dataExchangeModel_();
 }
 
 inline const IOModel& cfdemCloud::IOM() const
@@ -262,9 +262,9 @@ inline const IOModel& cfdemCloud::IOM() const
     return IOModel_;
 }
 
-inline const probeModel& cfdemCloud::probeM() const
+inline probeModel& cfdemCloud::probeM()
 {
-    return probeModel_;
+    return probeModel_();
 }
 
 inline const voidFractionModel& cfdemCloud::voidFractionM() const

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudIO.C

@@ -48,7 +48,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * IO * * * * * * * * * * * * * //
 
-void cfdemCloud::writeScalarFieldToTerminal(double**& array)
+void cfdemCloud::writeScalarFieldToTerminal(double**& array) const
 {
     // init double array
     for (int i=0; i<numberOfParticles(); i++)
@@ -57,7 +57,7 @@ void cfdemCloud::writeScalarFieldToTerminal(double**& array)
     }
 }
 
-void cfdemCloud::writeVectorFieldToTerminal(double**& array)
+void cfdemCloud::writeVectorFieldToTerminal(double**& array) const
 {
     // init double array
     for (int i=0; i<numberOfParticles(); i++)

src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergy.C

@@ -92,12 +92,12 @@ void cfdemCloudEnergy::speciesExecute()
 }
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-int cfdemCloudEnergy::nrEnergyModels()
+label cfdemCloudEnergy::nrEnergyModels() const
 {
     return energyModels_.size();
 }
 
-bool& cfdemCloudEnergy::implicitEnergyModel()
+bool cfdemCloudEnergy::implicitEnergyModel() const
 {
     return implicitEnergyModel_;
 }

src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergy.H

@@ -53,21 +53,21 @@ class cfdemCloudEnergy
 protected:
 
     const wordList energyModels_;
-    
+
     bool implicitEnergyModel_;
-    
+
     const wordList chemistryModels_;
 
     autoPtr<energyModel>* energyModel_;
-    
+
     autoPtr<thermCondModel> thermCondModel_;
-    
+
     autoPtr<chemistryModel> chemistryModel_;
-    
+
     void calcEnergyContributions();
 
     void speciesExecute();
-    
+
 public:
 
     friend class energyModel;
@@ -89,23 +89,23 @@ public:
 
 
        const energyModel& energyM(int);
-       
+
        const thermCondModel& thermCondM();
-       
+
        const chemistryModel& chemistryM();
 
-       int nrEnergyModels();
+       label nrEnergyModels() const;
+
+       inline const wordList& energyModels() const;
+
+       bool implicitEnergyModel() const;
 
-       inline const wordList& energyModels();
-       
-       bool& implicitEnergyModel();
-       
        void energyContributions(volScalarField&);
-       
+
        void energyCoefficients(volScalarField&);
-       
+
        bool evolve(volScalarField&,volVectorField&,volVectorField&);
-       
+
        void postFlow();
 
 };

src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergyI.H

@@ -25,7 +25,7 @@ namespace Foam
 {
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-inline const wordList& cfdemCloudEnergy::energyModels()
+inline const wordList& cfdemCloudEnergy::energyModels() const
 {
     return energyModels_;
 }

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

@@ -91,7 +91,7 @@ bool cfdemCloudIB::reAllocArrays()
     if(cfdemCloud::reAllocArrays())
     {
         // get arrays of new length
-        dataExchangeM().allocateArray(angularVelocities_,0,3);
+        dataExchangeM().allocateArray(angularVelocities_,0.,3);
         return true;
     }
     return false;
@@ -224,11 +224,9 @@ void cfdemCloudIB::calcVelocityCorrection
 
 }
 
-vector cfdemCloudIB::angularVelocity(int index)
+vector cfdemCloudIB::angularVelocity(int index) const
 {
-    vector vel;
-    for(int i=0;i<3;i++) vel[i] = angularVelocities_[index][i];
-    return vel;
+    return vector(angularVelocities_[index][0],angularVelocities_[index][1],angularVelocities_[index][2]);
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.H

@@ -90,7 +90,7 @@ public:
         void calcVelocityCorrection(volScalarField&,volVectorField&,volScalarField&,volScalarField&); // this could be moved to an IB mom couple model
 
       // Access
-        vector angularVelocity(int);
+        vector angularVelocity(int) const;
 
         inline double ** angularVelocities() const
         {

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

@@ -128,14 +128,14 @@ void cfdemCloudMS::getDEMdata()
 
     //- save clump volume and mass
     double **typeDH(NULL);
-    dataExchangeM().allocateArray(typeDH,-1,1,nClumpTypes()+1);
+    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
     }
@@ -144,7 +144,7 @@ void cfdemCloudMS::getDEMdata()
     for(int ind = 0;ind < numberOfClumps(); ind++)
     {
         ct=clumpType()[0][ind];
-        clumpVol_[ind][0] = typeVol_[ct];      
+        clumpVol_[ind][0] = typeVol_[ct];
         clumpDH_[ind][0]=typeDH[ct][0];
         //Info << "ct=" << ct << endl;
         //Info << "clumpVol()[ind][0]=" << clumpVol()[ind][0] << endl;
@@ -192,18 +192,18 @@ bool cfdemCloudMS::reAllocArrays()
     if(cfdemCloud::reAllocArrays())
     {
         // get arrays of new length
-        dataExchangeM().allocateArray(positionsCM_,0,3,"nbodies");
-        dataExchangeM().allocateArray(velocitiesCM_,0,3,"nbodies");
+        dataExchangeM().allocateArray(positionsCM_,0.,3,"nbodies");
+        dataExchangeM().allocateArray(velocitiesCM_,0.,3,"nbodies");
         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(clumpVol_,0.,1,"nbodies");
         dataExchangeM().allocateArray(clumpDH_,1.,1,"nbodies");
-        dataExchangeM().allocateArray(clumpWeights_,1,1,"nbodies");
-        dataExchangeM().allocateArray(impForcesCM_,0,3,"nbodies");
-        dataExchangeM().allocateArray(expForcesCM_,0,3,"nbodies");
-        dataExchangeM().allocateArray(DEMForcesCM_,0,3,"nbodies");
+        dataExchangeM().allocateArray(clumpWeights_,1.,1,"nbodies");
+        dataExchangeM().allocateArray(impForcesCM_,0.,3,"nbodies");
+        dataExchangeM().allocateArray(expForcesCM_,0.,3,"nbodies");
+        dataExchangeM().allocateArray(DEMForcesCM_,0.,3,"nbodies");
         return true;
     }
     return false;

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

@@ -117,9 +117,9 @@ public:
 
     // Member Functions
       // Access
-       inline label body(int);
+       inline label body(int) const;
 
-       inline double particleVolume(int);
+       inline double particleVolume(int) const;
 
        inline vector positionCM(int);
 
@@ -127,9 +127,9 @@ public:
 
        inline label cellIDCM(int);
 
-       inline label nrigid(int);
+       inline label nrigid(int) const;
 
-       inline int nrForceModels();
+       inline label nrForceModels() const;
 
        inline double **& positionsCM() const;
 

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

@@ -36,12 +36,12 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline label cfdemCloudMS::body(int index)
+inline label cfdemCloudMS::body(int index) const
 {
     return bodies_[0][index]-1;
 }
 
-inline double cfdemCloudMS::particleVolume(int index)
+inline double cfdemCloudMS::particleVolume(int index) const
 {
     int ind = body(index); // particle to clump ID
 
@@ -50,7 +50,7 @@ inline double cfdemCloudMS::particleVolume(int index)
         Vp=clumpVol_[ind][0];
     else
         Warning << "ind=" << ind << endl;
-    
+
     int nR(nrigid(ind));
     if(nR>0) Vp/=nR;
     return Vp;
@@ -75,12 +75,12 @@ inline label cfdemCloudMS::cellIDCM(int index)
     return cellIDsCM_[index][0];
 }
 
-inline label cfdemCloudMS::nrigid(int index)
+inline label cfdemCloudMS::nrigid(int index) const
 {
     return nrigids_[0][index];
 }
 
-inline int cfdemCloudMS::nrForceModels()
+inline label cfdemCloudMS::nrForceModels() const
 {
     return forceModels_.size();
 }

src/lagrangian/cfdemParticle/etc/bashrc

@@ -1,216 +0,0 @@
-#----------------------------------*-sh-*--------------------------------------
-# CFDEMcoupling
-# Christoph Goniva
-# June 2012
-#------------------------------------------------------------------------------
-#
-# Script
-#     etc/bashrc
-#
-# Description
-#     Startup file for cfdem exporting environment variables
-#     Sourced  ~/.bashrc
-#
-#- adapt this and add to ./bashrc 
-#- you can test the correctness using cfdemSystemTest.sh
-
-##================================================#
-##- source cfdem env vars
-#export CFDEM_VERSION=PUBLIC
-#export CFDEM_PROJECT_DIR=$HOME/CFDEM/CFDEMcoupling-$CFDEM_VERSION-$WM_PROJECT_VERSION
-#export CFDEM_SRC_DIR=$CFDEM_PROJECT_DIR/src
-#export CFDEM_SOLVER_DIR=$CFDEM_PROJECT_DIR/applications/solvers
-#export CFDEM_DOC_DIR=$CFDEM_PROJECT_DIR/doc
-#export CFDEM_UT_DIR=$CFDEM_PROJECT_DIR/applications/utilities
-#export CFDEM_TUT_DIR=$CFDEM_PROJECT_DIR/tutorials
-#export CFDEM_PROJECT_USER_DIR=$HOME/CFDEM/$LOGNAME-$CFDEM_VERSION-$WM_PROJECT_VERSION
-#export CFDEM_bashrc=$CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/bashrc
-#export CFDEM_LIGGGHTS_SRC_DIR=$HOME/LIGGGHTS/LIGGGHTS-PUBLIC/src
-#export CFDEM_LIGGGHTS_BIN_DIR=$HOME/LIGGGHTS/LIGGGHTS-PUBLIC/src-build
-#export CFDEM_LIGGGHTS_MAKEFILE_NAME=fedora_fpic # obsolete with CMake
-#export CFDEM_LPP_DIR=$HOME/LIGGGHTS/mylpp/src
-#export CFDEM_PIZZA_DIR=$HOME/LIGGGHTS/PIZZA/gran_pizza_17Aug10/src
-#. $CFDEM_bashrc
-#================================================#
-#------------------------------------------------------------------------------
-
-#------------------------------------------------------------------------------
-#- 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 CFDEM_LAMMPS_LIB_DIR=$CFDEM_LAMMPS_LIB_DIR defined by user."
-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
-
-#check if additional libraries should be compiled together with solvers
-if [[ $CFDEM_ADD_LIBS_DIR == "" ]]; then
-    export CFDEM_ADD_LIBS_DIR=$CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
-else
-    echo "using CFDEM_ADD_LIBS_DIR=$CFDEM_ADD_LIBS_DIR defined by user."
-fi
-
-#-----------------------------------------------------
-# 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
-
-#- LMP M2M lib path and makefile
-export CFDEM_M2MLIB_PATH=$CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayM2M/library
-export CFDEM_M2MLIB_MAKEFILENAME=$CFDEM_LIGGGHTS_MAKEFILE_NAME
-
-#- LMP POEMS lib path and makefile
-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
-
-#- path to libraries
-export CFDEM_LIB_DIR=$FOAM_USER_LIBBIN
-
-#- path to apps
-export CFDEM_APP_DIR=$FOAM_USER_APPBIN
-
-#- path to OF version flag file
-export CFDEM_OFVERSION_DIR=$CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/OFversion
-#------------------------------------------------------------------------------
-#- settings for lpp postproc tool
-#------------------------------------------------------------------------------
-
-#- nr of procs for lpp tool
-export CFDEM_LPP_NPROCS=4
-
-#- nr of procs for lpp tool
-export CFDEM_LPP_CHUNKSIZE=1
-
-#- shortcut to run lpp
-alias lpp='python -i $CFDEM_LPP_DIR/lpp.py --cpunum $CFDEM_LPP_NPROCS --chunksize $CFDEM_LPP_CHUNKSIZE'
-
-#------------------------------------------------------------------------------
-#- aliases for easy navigation (no changes necessary)
-#------------------------------------------------------------------------------
-#- shortcut to cfdem path
-alias cfdem='cd $CFDEM_PROJECT_DIR'
-
-#- shortcut to src path
-alias cfdemSrc='cd $CFDEM_SRC_DIR'
-
-#- shortcut to tutorial path
-alias cfdemTut='cd $CFDEM_TUT_DIR'
-
-#- shortcut to solver path
-alias cfdemSol='cd $CFDEM_SOLVER_DIR'
-
-#- shortcut to utilities path
-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'
-
-#- shortcut to open the doxygen with firefox
-alias cfdemDox='firefox $CFDEM_DOC_DIR/doxygen/html/index.html'
-
-#- shortcut to LIGGGHTS path
-alias cfdemLIG='cd $CFDEM_LIGGGHTS_SRC_DIR'
-
-#- shortcut to system test
-alias cfdemSysTest='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/cfdemSystemTest.sh'
-
-#- shortcut to pull LIGGGHTS
-alias cfdemPullLIG='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/pullLIGGGHTS.sh'
-
-#- shortcut to pull CFDEMcoupling
-alias cfdemPullCFDEMcoupling='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/pullCFDEMcoupling.sh'
-
-#- shortcut to clean CFDEM
-alias cfdemCleanCFDEM='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/cleanCFDEMcoupling.sh'
-
-#- shortcut to compile LIGGGHTS + sublibraries
-alias cfdemCompLIG='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileLIGGGHTS.sh'
-
-#- shortcut to compile CFDEMcoupling +LIGGGHTS
-alias cfdemCompCFDEMall='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_all.sh'
-
-#- shortcut to compile CFDEMcoupling (src+solvers)
-alias cfdemCompCFDEM='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling.sh'
-
-#- shortcut to compile CFDEMcoupling src
-alias cfdemCompCFDEMsrc='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_src.sh'
-
-#- shortcut to compile CFDEMcoupling solvers
-alias cfdemCompCFDEMsol='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_sol.sh'
-
-#- shortcut to compile CFDEMcoupling utilities
-alias cfdemCompCFDEMuti='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_uti.sh'
-
-#- 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 {} \;'
-
-#- shortcut to run liggghts in serial
-cfdemLiggghts() { $CFDEM_LIGGGHTS_SRC_DIR/lmp_$CFDEM_LIGGGHTS_MAKEFILE_NAME < $1; }
-export -f cfdemLiggghts
-
-#- shortcut to run liggghts in parallel
-cfdemLiggghtsPar() { mpirun -np $2 $CFDEM_LIGGGHTS_SRC_DIR/lmp_$CFDEM_LIGGGHTS_MAKEFILE_NAME < $1; }
-export -f cfdemLiggghtsPar
-
-#- shortcut to open files including a pattern
-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
-    :
-else
-        echo "make new dirs $CFDEM_PROJECT_USER_DIR ? (y/n)"
-        read YN
-        if [ $YN == "y" ]; then
-            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
-        fi
-fi

src/lagrangian/cfdemParticle/etc/cshrc

@@ -1,202 +0,0 @@
-#----------------------------------*-sh-*--------------------------------------
-# CFDEMcoupling
-# Christoph Goniva
-# June 2012
-#------------------------------------------------------------------------------
-#
-# Script
-#     etc/cshrc
-#
-# Description
-#     Startup file for cfdem exporting environment variables
-#     Sourced  ~/.cshrc
-#
-#- adapt this and add to ./cshrc 
-#- you can test the correctness using cfdemSystemTest.sh
-
-##================================================#
-##- source cfdem env vars
-#setenv CFDEM_VERSION PUBLIC
-#setenv CFDEM_PROJECT_DIR $HOME/CFDEM/CFDEMcoupling-$CFDEM_VERSION-$WM_PROJECT_VERSION
-#setenv CFDEM_SRC_DIR $CFDEM_PROJECT_DIR/src/lagrangian/cfdemParticle
-#setenv CFDEM_SOLVER_DIR $CFDEM_PROJECT_DIR/applications/solvers
-#setenv CFDEM_DOC_DIR $CFDEM_PROJECT_DIR/doc
-#setenv CFDEM_UT_DIR $CFDEM_PROJECT_DIR/applications/utilities
-#setenv CFDEM_TUT_DIR $CFDEM_PROJECT_DIR/tutorials
-#setenv CFDEM_PROJECT_USER_DIR $HOME/CFDEM/$LOGNAME-$CFDEM_VERSION-$WM_PROJECT_VERSION
-#setenv CFDEM_bashrc $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/cshrc
-#setenv CFDEM_LIGGGHTS_SRC_DIR $HOME/LIGGGHTS/LIGGGHTS-PUBLIC/src
-#setenv CFDEM_LIGGGHTS_MAKEFILE_NAME fedora_fpic
-#setenv CFDEM_LPP_DIR $HOME/LIGGGHTS/mylpp/src
-#setenv CFDEM_PIZZA_DIR $HOME/LIGGGHTS/PIZZA/gran_pizza_17Aug10/src
-#source $CFDEM_bashrc
-#================================================#
-#------------------------------------------------------------------------------
-
-#------------------------------------------------------------------------------
-#- export environment variables (adapt to your paths)
-#------------------------------------------------------------------------------
-
-#check if default lammps lib path should be used
-if ( ! ($?CFDEM_LAMMPS_LIB_DIR) ) then
-    setenv CFDEM_LAMMPS_LIB_DIR $CFDEM_LIGGGHTS_SRC_DIR"/../lib/"
-else
-    echo "using CFDEM_LAMMPS_LIB_DIR=$CFDEM_LAMMPS_LIB_DIR defined by user."
-endif
-
-#- LIGGGHTS lib name
-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
-
-#check if additional libraries should be compiled together with solvers
-if ( ! ($?CFDEM_ADD_LIBS_DIR) ) then
-    setenv CFDEM_ADD_LIBS_DIR $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
-else
-    echo "using CFDEM_ADD_LIBS_DIR=$CFDEM_ADD_LIBS_DIR defined by user."
-endif
-
-#-----------------------------------------------------
-# additional libraries
-
-#- 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
-
-#- LMP M2M lib path
-setenv CFDEM_M2MLIB_PATH $CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayM2M/library
-setenv CFDEM_M2MLIB_MAKEFILENAME $CFDEM_LIGGGHTS_MAKEFILE_NAME
-
-#- LMP POEMS lib path
-setenv CFDEM_POEMSLIB_PATH $CFDEM_LIGGGHTS_SRC_DIR/../lib/poems
-setenv CFDEM_POEMSLIB_MAKEFILENAME g++
-
-#- path to test harness
-setenv CFDEM_TEST_HARNESS_PATH $CFDEM_PROJECT_USER_DIR/log/logFilesCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
-
-#- path to libraries
-setenv CFDEM_LIB_DIR $FOAM_USER_LIBBIN
-
-#- path to apps
-setenv CFDEM_APP_DIR $FOAM_USER_APPBIN
-
-#------------------------------------------------------------------------------
-#- settings for lpp postproc tool
-#------------------------------------------------------------------------------
-
-#- nr of procs for lpp tool
-setenv CFDEM_LPP_NPROCS 4
-
-#- nr of procs for lpp tool
-setenv CFDEM_LPP_CHUNKSIZE 1
-
-#- shortcut to run lpp
-alias lpp 'python -i $CFDEM_LPP_DIR/lpp.py --cpunum $CFDEM_LPP_NPROCS --chunksize $CFDEM_LPP_CHUNKSIZE \!:1'
-
-#------------------------------------------------------------------------------
-#- aliases for easy navigation (no changes necessary)
-#------------------------------------------------------------------------------
-#- shortcut to cfdem path
-alias cfdem 'cd $CFDEM_PROJECT_DIR'
-
-#- shortcut to src path
-alias cfdemSrc 'cd $CFDEM_SRC_DIR'
-
-#- shortcut to tutorial path
-alias cfdemTut 'cd $CFDEM_TUT_DIR'
-
-#- shortcut to solver path
-alias cfdemSol 'cd $CFDEM_SOLVER_DIR'
-
-#- shortcut to utilities path
-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'
-
-#- shortcut to open the doxygen with firefox
-alias cfdemDox 'firefox $CFDEM_DOC_DIR/doxygen/html/index.html'
-
-#- shortcut to LIGGGHTS path
-alias cfdemLIG 'cd $CFDEM_LIGGGHTS_SRC_DIR'
-
-#- shortcut to system test
-alias cfdemSysTest 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/cfdemSystemTest.sh'
-
-#- shortcut to pull LIGGGHTS
-alias cfdemPullLIG 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/pullLIGGGHTS.sh'
-
-#- shortcut to pull CFDEMcoupling
-alias cfdemPullCFDEMcoupling 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/pullCFDEMcoupling.sh'
-
-#- shortcut to clean CFDEM
-alias cfdemCleanCFDEM 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/cleanCFDEMcoupling.sh'
-
-#- shortcut to compile LIGGGHTS + sublibraries
-alias cfdemCompLIG 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileLIGGGHTS.sh'
-
-#- shortcut to compile CFDEMcoupling +LIGGGHTS
-alias cfdemCompCFDEMall 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_all.sh'
-
-#- shortcut to compile CFDEMcoupling (src+solvers)
-alias cfdemCompCFDEM 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling.sh'
-
-#- shortcut to compile CFDEMcoupling src
-alias cfdemCompCFDEMsrc 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_src.sh'
-
-#- shortcut to compile CFDEMcoupling solvers
-alias cfdemCompCFDEMsol 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_sol.sh'
-
-#- shortcut to compile CFDEMcoupling utilities
-alias cfdemCompCFDEMuti 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_uti.sh'
-
-#- 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'
-
-#- shortcut to run liggghts in parallel (no fcts in csh...)
-alias cfdemLiggghtsPar 'echo "mpirun -np xx -machinefile mynodes $CFDEM_LIGGGHTS_SRC_DIR/lmp_$CFDEM_LIGGGHTS_MAKEFILE_NAME < in.liggghts_init"'
-
-#- shortcut to run liggghts in parallel
-# unfortunately no functions available in csh
-
-#- shortcut to run lpp
-alias lpp 'python -i $CFDEM_LPP_DIR/lpp.py \!:1'
-
-# check if the directory exists
-if ( -d "$CFDEM_PROJECT_USER_DIR" ) then
-    :
-else
-    if ( -d "$CFDEM_PROJECT_DIR" ) then
-        echo "make new dirs $CFDEM_PROJECT_USER_DIR ? (y/n)"
-	set YN=$<
-        if ( $YN == "y" ) then
-            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 cshrc."
-        exit
-    endif
-endif

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

@@ -78,55 +78,71 @@ fileName IOModel::createLagrangianDir(fileName path) const
 fileName IOModel::buildFilePath(word dirName) const
 {
     // create file structure
-	fileName path("");
+    fileName path("");
     if(parOutput_)
     {
-    	path=fileName(particleCloud_.mesh().time().path()/particleCloud_.mesh().time().timeName()/dirName/"particleCloud");
-    	mkDir(path,0777);
-    } else
+        path = fileName(particleCloud_.mesh().time().path()/particleCloud_.mesh().time().timeName()/dirName/"particleCloud");
+        mkDir(path,0777);
+    }
+    else
     {
-		path=fileName("."/dirName);
-    	mkDir(path,0777);
-    	mkDir(fileName(path/"constant"),0777);
-    	OFstream* stubFile = new OFstream(fileName(path/"particles.foam"));
-    	delete stubFile;
-        }
+        path = fileName("."/dirName);
+        mkDir(path,0777);
+        mkDir(path/"constant",0777);
+        OFstream stubFile(path/"particles.foam");
+    }
     return path;
 }
 
-void IOModel::streamDataToPath(fileName path, double** array,int nPProc,word name,word type,word className,word finaliser) const
+void IOModel::streamDataToPath(fileName path, double** array,int nPProc,word name,word type,word className) const
 {
-    vector vec;
-    OFstream* fileStream = new OFstream(fileName(path/name));
-    *fileStream << "FoamFile\n";
-    *fileStream << "{version 2.0; format ascii;class "<< className << "; location 0;object  "<< name <<";}\n";
-    *fileStream << nPProc <<"\n";
+    OFstream fileStream(path/name);
 
+    fileStream
+        << "FoamFile\n{\n"
+        << "    version     " << fileStream.version() << ";\n"
+        << "    format      " << fileStream.format() << ";\n"
+        << "    class       " << className << ";\n"
+        << "    location    " << 0 << ";\n"
+        << "    object      " << name << ";\n"
+        << "}" << nl;
 
-    if(type!="origProcId")*fileStream << "(\n";
-    else if(type=="origProcId")
+    fileStream << nPProc <<"\n";
+
+    if (type == "origProcId")
     {
-        if(nPProc>0) *fileStream <<"{0}"<< "\n";
-        else *fileStream <<"{}"<< "\n";
+        if (nPProc > 0) fileStream << "{0}" << "\n";
+        else fileStream << "{}" << "\n";
+        return;
     }
 
-    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    fileStream << token::BEGIN_LIST << nl;
+
+    int ** cellIDs = particleCloud_.cellIDs();
+    for (int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
-        if (particleCloud_.cellIDs()[index][0] > -1) // particle Found
+        if (cellIDs[index][0] > -1) // particle Found
         {
-            if (type=="scalar"){
-                *fileStream << array[index][0] << " \n";
-            }else if (type=="position" || type=="vector"){
-                for(int i=0;i<3;i++) vec[i] = array[index][i];
-                *fileStream <<"( "<< vec[0] <<" "<<vec[1]<<" "<<vec[2]<<" ) "<< finaliser << " \n";
-            }else if (type=="label"){
-                *fileStream << index << finaliser << " \n";
+            if (type == "scalar")
+            {
+                fileStream << array[index][0] << " \n";
+            }
+            else if (type == "position")
+            {
+                fileStream << "( "<< array[index][0] << " " << array[index][1] << " " << array[index][2] << " ) " << cellIDs[index][0] << " \n";
+            }
+            else if (type == "vector")
+            {
+                fileStream << "( "<< array[index][0] << " " << array[index][1] << " " << array[index][2] << " ) " << " \n";
+            }
+            else if (type == "label")
+            {
+                fileStream << index << " \n";
             }
         }
     }
 
-    if(type!="origProcId")*fileStream << ")\n";
-    delete fileStream;
+    fileStream << token::END_LIST << nl;
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -143,12 +159,12 @@ IOModel::IOModel
     time_(sm.mesh().time()),
     parOutput_(true)
 {
-	if (
+    if (
             particleCloud_.dataExchangeM().myType()=="oneWayVTK" ||
             dict_.found("serialOutput")
        )
     {
-        parOutput_=false;
+        parOutput_ = false;
         Warning << "IO model is in serial write mode, only data on proc 0 is written" << endl;
     }
 }

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

@@ -121,7 +121,7 @@ public:
 
     fileName buildFilePath(word) const;
 
-    void streamDataToPath(fileName,double**,int,word,word type,word className,word finaliser) const;
+    void streamDataToPath(fileName,double**,int,word,word type,word className) const;
 };
 
 

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

@@ -92,19 +92,19 @@ int basicIO::dumpDEMdata() const
         {
             Info << "createTimeDir(path_), path="<<path_ << endl;
             Info << "lagPath_=createTimeDir(fileName(lagPath_/lagrangian)), lagPath="<<path_ << endl;
-        	lagPath_=createTimeDir(path_);
-        	lagPath_=createTimeDir(fileName(lagPath_/"lagrangian"));
+            lagPath_ = createTimeDir(path_);
+            lagPath_ = createTimeDir(fileName(lagPath_/"lagrangian"));
         }
         // calc the number of particles on proc
         int count(0);
         for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
             if (particleCloud_.cellIDs()[index][0] > -1) count++;
-        nPProc_=count;
-        
+        nPProc_ = count;
+
         // stream data to file
-        streamDataToPath(lagPath_, particleCloud_.positions(),nPProc_,"positions","vector","Cloud<passiveParticle>","0");
-        streamDataToPath(lagPath_, particleCloud_.velocities(),nPProc_,"v","vector","vectorField","");
-        streamDataToPath(lagPath_, particleCloud_.radii(),nPProc_,"r","scalar","scalarField","");
+        streamDataToPath(lagPath_, particleCloud_.positions(),nPProc_,"positions","position","Cloud<passiveParticle>");
+        streamDataToPath(lagPath_, particleCloud_.velocities(),nPProc_,"v","vector","vectorField");
+        streamDataToPath(lagPath_, particleCloud_.radii(),nPProc_,"r","scalar","scalarField");
     }
     return nPProc_;
 }

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

@@ -84,10 +84,10 @@ int sophIO::dumpDEMdata() const
         npProcs=basicIO::dumpDEMdata();
 
         // stream data to file
-        streamDataToPath(lagPath_, particleCloud_.voidfractions(),npProcs,"voidfractions","scalar","scalarField","");
-        streamDataToPath(lagPath_, particleCloud_.impForces(),npProcs,"impForces","vector","vectorField","");
-        streamDataToPath(lagPath_, particleCloud_.expForces(),npProcs,"expForces","vector","vectorField","");
-        streamDataToPath(lagPath_, particleCloud_.DEMForces(),npProcs,"DEMForces","vector","vectorField","");
+        streamDataToPath(lagPath_, particleCloud_.voidfractions(),npProcs,"voidfractions","scalar","scalarField");
+        streamDataToPath(lagPath_, particleCloud_.impForces(),npProcs,"impForces","vector","vectorField");
+        streamDataToPath(lagPath_, particleCloud_.expForces(),npProcs,"expForces","vector","vectorField");
+        streamDataToPath(lagPath_, particleCloud_.DEMForces(),npProcs,"DEMForces","vector","vectorField");
     }
     return npProcs;
 }

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

@@ -82,8 +82,8 @@ int trackIO::dumpDEMdata() const
         npProcs = sophIO::dumpDEMdata();
 
         // stream data to file
-        streamDataToPath(lagPath_, particleCloud_.velocities(),npProcs,"origId","label","labelField","");
-        streamDataToPath(lagPath_, particleCloud_.velocities(),npProcs,"origProcId","origProcId","labelField","");
+        streamDataToPath(lagPath_, particleCloud_.velocities(),npProcs,"origId","label","labelField");
+        streamDataToPath(lagPath_, particleCloud_.velocities(),npProcs,"origProcId","origProcId","labelField");
     }
 
     return npProcs;

src/lagrangian/cfdemParticle/subModels/averagingModel/averagingModel/averagingModel.C

@@ -262,7 +262,7 @@ void averagingModel::setDSauter
     for(int index=0; index< particleCloud_.numberOfParticles(); index++)
     {
         if(myParticleType!=0) //in case a particle type is specified, only consider particles of the right type
-            if(myParticleType != particleCloud_.particleType(index)) continue; 
+            if(myParticleType != particleCloud_.particleType(index)) continue;
 
         radius         = particleCloud_.radii()[index][0] / scale_; //the primary particle diameter
         radiusPow2 = radius*radius;
@@ -313,7 +313,7 @@ void averagingModel::resetWeightFields() const
 }
 
 
-void Foam::averagingModel::undoWeightFields(double**const& mask) const
+void averagingModel::undoWeightFields(double**const& mask) const
 {
     for(int index=0; index< particleCloud_.numberOfParticles(); index++)
     {
@@ -326,41 +326,14 @@ void Foam::averagingModel::undoWeightFields(double**const& mask) const
     }
 }
 
-tmp<volVectorField> Foam::averagingModel::UsInterp() const
+tmp<volVectorField> averagingModel::UsInterp() const
 {
-    tmp<volVectorField> tsource
+    const scalar tsf = particleCloud_.dataExchangeM().timeStepFraction();
+
+    return tmp<volVectorField>
     (
-        new volVectorField
-        (
-            IOobject
-            (
-                "Us_averagingModel",
-                particleCloud_.mesh().time().timeName(),
-                particleCloud_.mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            particleCloud_.mesh(),
-            dimensionedVector
-            (
-                "zero",
-                dimensionSet(0, 1, -1, 0, 0),
-                vector::zero
-            )
-        )
+        new volVectorField("Us_averagingModel", (1. - tsf) * UsPrev_ + tsf * UsNext_)
     );
-
-    if (particleCloud_.dataExchangeM().couplingStep() > 1)
-    {
-        tsource.ref() = (1 - particleCloud_.dataExchangeM().timeStepFraction()) * UsPrev_
-                    + particleCloud_.dataExchangeM().timeStepFraction() * UsNext_;
-    }
-    else
-    {
-        tsource.ref() = UsNext_;
-    }
-
-    return tsource;
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C

@@ -147,11 +147,11 @@ species::species
 
 species::~species()
 {
-    delete partTemp_;
-    delete partRho_;
+    particleCloud_.dataExchangeM().destroy(partTemp_,1);
+    particleCloud_.dataExchangeM().destroy(partRho_,1);
 
-    for (int i=0; i<speciesNames_.size();i++) delete [] concentrations_[i];
-    for (int i=0; i<speciesNames_.size();i++) delete [] changeOfSpeciesMass_[i];
+    for (int i=0; i<speciesNames_.size();i++) particleCloud_.dataExchangeM().destroy(concentrations_[i],1);
+    for (int i=0; i<speciesNames_.size();i++) particleCloud_.dataExchangeM().destroy(changeOfSpeciesMass_[i],1);
 
 }
 
@@ -176,10 +176,10 @@ void species::allocateMyArrays() const
 
 void species::execute()
 {
-  // realloc the arrays
+    // realloc the arrays
     allocateMyArrays();
 
-  // get Y_i, T, rho at particle positions, fill arrays with them and push to LIGGGHTS
+    // get Y_i, T, rho at particle positions, fill arrays with them and push to LIGGGHTS
 
     label  cellI=0;
     scalar Tfluid(0);
@@ -194,7 +194,7 @@ void species::execute()
     for (int index=0; index<particleCloud_.numberOfParticles(); index++)
     {
         cellI=particleCloud_.cellIDs()[index][0];
-        if (cellI >=0)
+        if (cellI >= 0)
         {
             if(interpolation_)
             {
@@ -213,12 +213,12 @@ void species::execute()
             partRho_[index][0]=rhofluid;
             for (int i=0; i<speciesNames_.size();i++)
             {
-	        Yfluid_[i] = Y_[i][cellI];
-            concentrations_[i][index][0]=Yfluid_[i];
+                Yfluid_[i] = Y_[i][cellI];
+                concentrations_[i][index][0] = Yfluid_[i];
             }
         }
 
-        if(particleCloud_.verbose() && index >=0 && index < 2)
+        if(particleCloud_.verbose() && index >= 0 && index < 2)
         {
             /*for(int i =0; i<speciesNames_.size();i++)
             {

src/lagrangian/cfdemParticle/subModels/clockModel/clockModel/clockModel.C

@@ -33,7 +33,7 @@ Description
 #include <mpi.h>
 #include "clockModel.H"
 #include <unistd.h>
-
+#include <time.h>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,19 +48,18 @@ defineRunTimeSelectionTable(clockModel, dictionary);
 
 // * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * //
 
-void Foam::clockModel::start(int pos) const
+void clockModel::start(int pos) const
 {
     start(pos,"");
-    return;
 }
 
-void Foam::clockModel::start(int pos, const std::string& ident) const
+void clockModel::start(int pos, const std::string& ident) const
 {
-    if (particleCloud_.mesh().time().value() > startTime_)
+    if (time_.value() > startTime_)
     {
         if (pos >= n_) // alternatively one fixed size?
         {
-            n_ = 2*n_;
+            n_ += 32;
             deltaT_.resize(n_,0);
             identifier_.resize(n_,"");
             nOfRuns_.resize(n_,0);
@@ -73,14 +72,13 @@ void Foam::clockModel::start(int pos, const std::string& ident) const
         parent_[pos]=curParent_;
         curParent_ = pos;
         nOfRuns_[pos] += 1;
-        deltaT_[pos]-=std::clock();
+        deltaT_[pos] -= std::clock();
     }
-    return;
 }
 
-void Foam::clockModel::stop() const
+void clockModel::stop() const
 {
-    if (particleCloud_.mesh().time().value() > startTime_)
+    if (time_.value() > startTime_)
     {
         deltaT_[curParent_] += std::clock();
         curLev_ -= 1;
@@ -93,17 +91,16 @@ void Foam::clockModel::stop() const
             curParent_ = -1;
         }
     }
-    return;
 }
 
-void Foam::clockModel::stop(const std::string& ident) const
+void clockModel::stop(const std::string& ident) const
 {
-    if (particleCloud_.mesh().time().value() > startTime_)
+    if (time_.value() > startTime_)
     {
         deltaT_[curParent_] += std::clock();
         if (curParent_ > 0 && identifier_[curParent_].compare(ident) != 0)
         {
-            Pout<<"Warning: stop identifier did not equal start identifier! "<<ident<<" & "<<identifier_[curParent_]<<nl;
+            Pout << "Warning: stop identifier did not equal start identifier! " << ident << " & " << identifier_[curParent_] << nl;
         }
         curLev_ -= 1;
         if (curParent_ >= 0)
@@ -115,10 +112,9 @@ void Foam::clockModel::stop(const std::string& ident) const
             curParent_ = -1;
         }
     }
-    return;
 }
 
-std::string Foam::clockModel::eval() const
+std::string clockModel::eval() const
 {
     std::ostringstream strs("Measurements in CPU-seconds:\n");
     strs << "Name\tdeltaT\tnOfRuns\tlevel\tparentNr\tparentName\n";
@@ -149,7 +145,7 @@ std::string Foam::clockModel::eval() const
     return strs.str();
 }
 
-void Foam::clockModel::evalFile() const
+void clockModel::evalFile() const
 {
     std::ofstream outFile;
     std::string fileName(path_/"timeEval.txt");
@@ -159,10 +155,10 @@ void Foam::clockModel::evalFile() const
     outFile.close();
 }
 
-void Foam::clockModel::evalPar() const
+void clockModel::evalPar() const
 {
     int myrank, numprocs;
-    MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
+    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
     MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
 
     std::ofstream outFile;
@@ -227,29 +223,27 @@ void Foam::clockModel::evalPar() const
         outFile << strs.str();
         outFile.close();
     }
-
-    return;
 }
 
 
-void Foam::clockModel::initElems()
+void clockModel::initElems()
 {
     //init elems
-    for (int i = 0;i < n_; i++)
+    for (int i = 0; i < n_; ++i)
     {
         deltaT_[i] = 0;
-        identifier_[i] = "";
+        identifier_[i].clear();
         nOfRuns_[i] = 0;
         level_[i] = -1;
         parent_[i] = -2;
     }
 }
 
-std::vector<int> Foam::clockModel::calcShift() const
+std::vector<int> clockModel::calcShift() const
 {
     std::vector<int> shifts(n_, 0);
 
-    for (int i=1; i<n_; i++)
+    for (int i=1; i<n_; ++i)
     {
         if (parent_[i] == -2)
         {
@@ -263,7 +257,7 @@ std::vector<int> Foam::clockModel::calcShift() const
     return shifts;
 }
 
-void Foam::clockModel::normHist() const
+void clockModel::normHist() const
 {
     int myrank, numprocs;
     MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
@@ -300,28 +294,28 @@ void Foam::clockModel::normHist() const
     Info << "===========================" << endl;
 
     getRAMUsage();
-    return;
 }
 
-void Foam::clockModel::plotHist(double buffIn,const std::string& identifier,int numprocs,int myrank) const
+void clockModel::plotHist(double buffIn,const std::string& identifier,int numprocs,int myrank) const
 {
     double* globalTime_all = NULL;
+
     if (myrank == 0) globalTime_all = new double[numprocs];
     MPI_Gather(&buffIn, 1, MPI_DOUBLE, globalTime_all, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
 
     if (myrank == 0)
         for (int j=0; j<numprocs; j++)
-            printf("%4f  ",globalTime_all[j]);
+            printf("%4f  ", globalTime_all[j]);
 
     Info << "\t" << identifier << endl;
 
     delete [] globalTime_all;
 }
 
-void Foam::clockModel::Hist() const
+void clockModel::Hist() const
 {
-    int myrank=-10;
-    MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
+    int myrank = -1;
+    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
 
     //Global = 1 / Coupling = 2 / LIGGGHTS = 3 /Flow = 26
 
@@ -333,11 +327,9 @@ void Foam::clockModel::Hist() const
     Pout << "[" << myrank << "]: " << "Coupling - LIGGGHTS" << " " << ((deltaT_[2]-deltaT_[3])/CLOCKS_PER_SEC) << '\n';
     //Flow = 26
     Pout << "[" << myrank << "]: " << identifier_[26] << " " << (deltaT_[26]/CLOCKS_PER_SEC) << '\n';
-
-    return;
 }
 
-void Foam::clockModel::getRAMUsage() const
+void clockModel::getRAMUsage() const
 {
     int myrank, numprocs;
     MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
@@ -360,6 +352,7 @@ void Foam::clockModel::getRAMUsage() const
     int SwapMem = 0;
     int temp = 0;
     strs.str("");
+
     if (inFile.is_open())	//search in File smaps for Rss and Swap entries
     {
         while (inFile.good())
@@ -382,15 +375,16 @@ void Foam::clockModel::getRAMUsage() const
             }
         }
     }
-    double SwapMB = static_cast<double>(SwapMem)/1024.0; //kB -> MB
-    double RssMB = static_cast<double>(RssMem)/1024.0;
+    double SwapMB = SwapMem/1024.0; //kB -> MB
+    double RssMB = RssMem/1024.0;
 
     inFile.close();
 
     // set up communication between Procs and plot Stuff
     Info << " RAM USAGE HISTOGRAM in MB" << endl;
     plotHist(RssMB,"RSS memory used",numprocs,myrank);
-    if (SwapMem > 0) 
+
+    if (SwapMem > 0)
     {
         plotHist(SwapMB,"WARNING: Swap",numprocs,myrank);
     }
@@ -398,26 +392,23 @@ void Foam::clockModel::getRAMUsage() const
 
     //Pout << "SWAP Memory used: " << SwapMem <<"MB\n";
     //Pout << "Rss Memory used: " << RssMem <<"MB\n";
-
-    return;
 }
 // * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * //
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
-Foam::clockModel::clockModel
+clockModel::clockModel
 (
     const dictionary& dict,
-    cfdemCloud& sm
+    const Time& time
 )
 :
     dict_(dict),
-    particleCloud_(sm),
+    time_(time),
     path_("clockData"),
-    startTime_(sm.mesh().time().startTime().value()+sm.mesh().time().deltaT().value()+SMALL),  // delay start of measurement by deltaT
-    //startTime_(0),                                //no delay
-    n_(30),
+    startTime_(time.startTime().value()+time.deltaT().value()+SMALL),  // delay start of measurement by deltaT
+    n_(32),
     deltaT_(n_),
     identifier_(n_),
     nOfRuns_(n_),
@@ -432,7 +423,7 @@ Foam::clockModel::clockModel
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-Foam::clockModel::~clockModel()
+clockModel::~clockModel()
 {}
 
 

src/lagrangian/cfdemParticle/subModels/clockModel/clockModel/clockModel.H

@@ -42,10 +42,8 @@ SourceFiles
 #define START(x) start(__COUNTER__,x)
 
 #include "fvCFD.H"
-#include "cfdemCloud.H"
-#include "dataExchangeModel.H"
-
 #include <vector>
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -61,7 +59,7 @@ protected:
 
     // Protected data
     const dictionary& dict_;
-    cfdemCloud& particleCloud_;
+    const Time& time_;
     fileName path_;
     scalar startTime_;
 
@@ -88,9 +86,9 @@ public:
             dictionary,
             (
                 const dictionary& dict,
-                cfdemCloud& sm
+                const Time& time
             ),
-            (dict,sm)
+            (dict,time)
         );
 
 
@@ -100,7 +98,7 @@ public:
         clockModel
         (
             const dictionary& dict,
-            cfdemCloud& sm
+            const Time& time
         );
 
 
@@ -114,7 +112,7 @@ public:
         static autoPtr<clockModel> New
         (
             const dictionary& dict,
-            cfdemCloud& sm
+            const Time& time
         );
 
 
@@ -129,9 +127,9 @@ public:
         virtual void evalPar() const;
         void initElems();
         std::vector<int> calcShift() const;	                //detects empty indices in vector, when times are evaluated
-        void Hist() const;                                  //calc Histogram
-        virtual void normHist() const;                      //calc normalized Histogram
-        void plotHist(double,const std::string&,int,int) const;    //plot histogramm to terminal
+        void Hist() const;                                      //calc Histogram
+        virtual void normHist() const;                          //calc normalized Histogram
+        void plotHist(double,const std::string&,int,int) const; //plot histogramm to terminal
         void getRAMUsage() const;
 
 };

src/lagrangian/cfdemParticle/subModels/clockModel/clockModel/newClockModel.C

@@ -44,7 +44,7 @@ namespace Foam
 autoPtr<clockModel> clockModel::New
 (
     const dictionary& dict,
-    cfdemCloud& sm
+    const Time& time
 )
 {
     word clockModelType
@@ -52,7 +52,7 @@ autoPtr<clockModel> clockModel::New
         dict.lookup("clockModel")
     );
 
-    Info<< "Selecting clockModel "
+    Info << "Selecting clockModel "
          << clockModelType << endl;
 
 
@@ -73,7 +73,7 @@ autoPtr<clockModel> clockModel::New
             << abort(FatalError);
     }
 
-    return autoPtr<clockModel>(cstrIter()(dict,sm));
+    return autoPtr<clockModel>(cstrIter()(dict,time));
 }
 
 

src/lagrangian/cfdemParticle/subModels/clockModel/noClock/noClock.C

@@ -56,10 +56,10 @@ addToRunTimeSelectionTable
 noClock::noClock
 (
     const dictionary& dict,
-    cfdemCloud& sm
+    const Time& time
 )
 :
-    clockModel(dict,sm)
+  clockModel(dict,time)
 {
     initElems();
 }

src/lagrangian/cfdemParticle/subModels/clockModel/noClock/noClock.H

@@ -47,7 +47,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class noDrag Declaration
+                           Class noClock Declaration
 \*---------------------------------------------------------------------------*/
 
 class noClock
@@ -67,7 +67,7 @@ public:
         noClock
         (
             const dictionary& dict,
-            cfdemCloud& sm
+            const Time& time
         );
 
     // Destructor

src/lagrangian/cfdemParticle/subModels/clockModel/standardClock/standardClock.C

@@ -30,7 +30,6 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "error.H"
-#include "IOModel.H"
 #include "standardClock.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -57,12 +56,14 @@ addToRunTimeSelectionTable
 standardClock::standardClock
 (
     const dictionary& dict,
-    cfdemCloud& sm
+    const Time& time
 )
 :
-    clockModel(dict,sm)
+  clockModel(dict,time)
 {
-    path_=particleCloud_.IOM().createTimeDir(path_);
+    path_ = path_/time_.timeName();
+    mkDir(path_,0777);
+
     initElems();
 }
 

src/lagrangian/cfdemParticle/subModels/clockModel/standardClock/standardClock.H

@@ -47,7 +47,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class noDrag Declaration
+                           Class standardClock Declaration
 \*---------------------------------------------------------------------------*/
 
 class standardClock
@@ -67,7 +67,7 @@ public:
         standardClock
         (
             const dictionary& dict,
-            cfdemCloud& sm
+            const Time& time
         );
 
     // Destructor

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.C

@@ -45,7 +45,7 @@ defineRunTimeSelectionTable(dataExchangeModel, dictionary);
 
 // * * * * * * * * * * * * * * protected Member Functions  * * * * * * * * * * * * * //
 
-void Foam::dataExchangeModel::setNumberOfParticles(int numberOfParticles) const
+void dataExchangeModel::setNumberOfParticles(int numberOfParticles) const
 {
     particleCloud_.setNumberOfParticles(numberOfParticles);
 }
@@ -55,7 +55,7 @@ void Foam::dataExchangeModel::setNumberOfParticles(int numberOfParticles) const
 //====
 // double **
 
-void Foam::dataExchangeModel::allocateArray
+void dataExchangeModel::allocateArray
 (
     double**& array,
     double initVal,
@@ -77,7 +77,7 @@ void Foam::dataExchangeModel::allocateArray
     }
 }
 
-void Foam::dataExchangeModel::allocateArray
+void dataExchangeModel::allocateArray
 (
     double**& array,
     double initVal,
@@ -92,7 +92,7 @@ void Foam::dataExchangeModel::allocateArray
     allocateArray(array,initVal,width,len);
 }
 
-void Foam::dataExchangeModel::destroy(double** array,int /*len*/) const
+void dataExchangeModel::destroy(double** array,int /*len*/) const
 {
     if (array == NULL) return;
 
@@ -103,7 +103,7 @@ void Foam::dataExchangeModel::destroy(double** array,int /*len*/) const
 
 //====
 // int **
-void Foam::dataExchangeModel::allocateArray
+void dataExchangeModel::allocateArray
 (
     int**& array,
     int initVal,
@@ -125,7 +125,7 @@ void Foam::dataExchangeModel::allocateArray
     }
 }
 
-void Foam::dataExchangeModel::allocateArray
+void dataExchangeModel::allocateArray
 (
     int**& array,
     int initVal,
@@ -140,7 +140,7 @@ void Foam::dataExchangeModel::allocateArray
     allocateArray(array,initVal,width,len);
 }
 
-void Foam::dataExchangeModel::destroy(int** array,int /*len*/) const
+void dataExchangeModel::destroy(int** array,int /*len*/) const
 {
     if (array == NULL) return;
 
@@ -152,7 +152,7 @@ void Foam::dataExchangeModel::destroy(int** array,int /*len*/) const
 
 //====
 // int *
-void Foam::dataExchangeModel::allocateArray
+void dataExchangeModel::allocateArray
 (
     int*& array,
     int initVal,
@@ -166,16 +166,15 @@ void Foam::dataExchangeModel::allocateArray
         array[i] = initVal;
 }
 
-void Foam::dataExchangeModel::destroy(int* array) const
+void dataExchangeModel::destroy(int* array) const
 {
-    if (array == NULL) return;
     delete [] array;
 }
 //====
 
 //====
 // double *
-void Foam::dataExchangeModel::allocateArray
+void dataExchangeModel::allocateArray
 (
     double*& array,
     double initVal,
@@ -189,16 +188,14 @@ void Foam::dataExchangeModel::allocateArray
         array[i] = initVal;
 }
 
-void Foam::dataExchangeModel::destroy(double* array) const
+void dataExchangeModel::destroy(double* array) const
 {
-    if (array == NULL) return;
-
     delete [] array;
 }
 //====
 
 
-bool Foam::dataExchangeModel::couple(int i) const
+bool dataExchangeModel::couple(int i) const
 {
     bool coupleNow = false;
     if (doCoupleNow())
@@ -209,36 +206,31 @@ bool Foam::dataExchangeModel::couple(int i) const
     return coupleNow;
 }
 
-scalar Foam::dataExchangeModel::timeStepFraction() const
+scalar dataExchangeModel::timeStepFraction() const
 {
     //return fraction between previous coupling TS and actual TS
-    //scalar DEMtime = DEMts_ * couplingInterval_;
-    //scalar frac = ( ( particleCloud_.mesh().time().value()-particleCloud_.mesh().time().startTime().value() ) - (couplingStep_) * DEMtime) / DEMtime; //Chr 05.03.2013
-    scalar frac = ( particleCloud_.mesh().time().value()-particleCloud_.mesh().time().startTime().value() - couplingStep_ * couplingTime() ) / couplingTime();
-    if (frac < 1e-4) frac = 1.;
-
-    return frac;
+    return ( particleCloud_.mesh().time().value()-particleCloud_.mesh().time().startTime().value() - (couplingStep_-1) * couplingTime() ) / couplingTime();
 }
 
-int Foam::dataExchangeModel::getNumberOfParticles() const
+int dataExchangeModel::getNumberOfParticles() const
 {
     Warning << "ask for nr of particles - which is not supported for this dataExchange model" << endl;
     return -1;
 }
 
-int Foam::dataExchangeModel::getNumberOfClumps() const
+int dataExchangeModel::getNumberOfClumps() const
 {
     Warning << "ask for nr of clumps - which is not supported for this dataExchange model" << endl;
     return -1;
 }
 
-int Foam::dataExchangeModel::getNumberOfTypes() const
+int dataExchangeModel::getNumberOfTypes() const
 {
     Warning << "ask for nr of types - which is not supported for this dataExchange model" << endl;
     return -1;
 }
 
-double* Foam::dataExchangeModel::getTypeVol() const
+double* dataExchangeModel::getTypeVol() const
 {
     Warning << "ask for type volume - which is not supported for this dataExchange model" << endl;
     return NULL;

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H

@@ -198,22 +198,29 @@ public:
 
         inline void checkTSsize() const
         {
-            if(particleCloud_.mesh().time().deltaT().value() >  couplingInterval_ * DEMts_ + SMALL)
+            if (particleCloud_.mesh().time().deltaT().value() >  couplingInterval_ * DEMts_ + SMALL)
             {
                 Info << "particleCloud_.mesh().time().deltaT().value() = " << particleCloud_.mesh().time().deltaT().value() << endl;
                 Info << "couplingInterval_ = " << couplingInterval_ << endl;
                 Info << "DEMts_ = " << DEMts_ << endl;
-                FatalError<<"\nError - TS bigger than coupling interval!\n"<< abort(FatalError);
+                FatalError << "\nError - CFD time-step bigger than coupling time (= DEM time step * coupling interval)!\n" << abort(FatalError);
             }
+            if (std::fabs((round(couplingTime()/particleCloud_.mesh().time().deltaT().value())*particleCloud_.mesh().time().deltaT().value())-couplingTime()) > SMALL)
+            {
+                Info << "particleCloud_.mesh().time().deltaT().value() = " << particleCloud_.mesh().time().deltaT().value() << endl;
+                Info << "couplingInterval_ = " << couplingInterval_ << endl;
+                Info << "DEMts_ = " << DEMts_ << endl;
+                Warning << "\nWarning - Coupling time (= DEM time step * coupling interval) is not a multiple of  CFD time-step!\n" << endl;
+            }
+            if (!particleCloud_.allowCFDsubTimestep())
+                if (particleCloud_.mesh().time().deltaT().value() < couplingInterval_ * DEMts_ + SMALL)
+                    FatalError << "\nYour models require: CFD time-step = coupling interval (= DEM time step * coupling interval)! \n" << abort(FatalError);
+
+            // warn if sub-TS
+            if (particleCloud_.mesh().time().deltaT().value() < couplingTime() - SMALL)
+                Warning << "You are using sub-time-steps (i.e. CFD TS < coupling time)! Check your settings properly." << endl;
         }
 
-        /*inline bool checkExactTiming() const
-        {
-            return false;
-        }*/
-
-        //void checkNClumpTypes() const {};
-
         inline void readDEMtsfromDict(dictionary& propsDict)
         {
             DEMts_ = readScalar(propsDict.lookup("DEMts"));
@@ -222,10 +229,8 @@ public:
 
         inline bool doCoupleNow() const
         {
-            if (particleCloud_.mesh().time().value()-particleCloud_.mesh().time().startTime().value()
-             - ((1+couplingStep_)*(DEMts_*couplingInterval_))
-             + 1e-10 > 0)                                         // Chr 27.03.2013 : first coupling after DEMts_*couplingInterval_
-//             > particleCloud_.mesh().time().deltaT().value()/2) // Chr 27.03.2013
+            if (particleCloud_.mesh().time().value()-particleCloud_.mesh().time().startTime().value()-SMALL
+                > couplingStep_*DEMts_*couplingInterval_)
             {
                 return true;
             }
@@ -245,16 +250,16 @@ public:
             for (int i=0;i<n;i++)
                 for (int j=0;j<3;j++)
                     particleCloud_.positions_[i][j]=pos[i*3+j];
-        };
+        }
         inline void setCellIDs(label n,int* ID) const
         {
             for (int i=0;i<n;i++)
                 particleCloud_.cellIDs_[i][0]=ID[i];
-        };
+        }
 
         virtual word myType() const=0;
 
-        virtual void setCG() const { Warning << "setCG() not executed correctly!" << endl; }
+        virtual void setCG() { Warning << "setCG() not executed correctly!" << endl; }
 };
 
 

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.C

@@ -127,7 +127,7 @@ void twoWayMPI::giveData
 
 //============
 // double **
-void Foam::twoWayMPI::allocateArray
+void twoWayMPI::allocateArray
 (
     double**& array,
     double initVal,
@@ -138,7 +138,7 @@ void Foam::twoWayMPI::allocateArray
     allocate_external_double(array, width, length, initVal, lmp);
 }
 
-void Foam::twoWayMPI::allocateArray
+void twoWayMPI::allocateArray
 (
     double**& array,
     double initVal,
@@ -149,7 +149,7 @@ void Foam::twoWayMPI::allocateArray
     allocate_external_double(array, width, length, initVal, lmp);
 }
 
-void Foam::twoWayMPI::destroy(double** array,int /*len*/) const
+void twoWayMPI::destroy(double** array,int /*len*/) const
 {
     if (array == NULL) return;
 
@@ -160,7 +160,7 @@ void Foam::twoWayMPI::destroy(double** array,int /*len*/) const
 
 //============
 // int **
-void Foam::twoWayMPI::allocateArray
+void twoWayMPI::allocateArray
 (
     int**& array,
     int initVal,
@@ -171,7 +171,7 @@ void Foam::twoWayMPI::allocateArray
     allocate_external_int(array, width, length, initVal, lmp);
 }
 
-void Foam::twoWayMPI::allocateArray
+void twoWayMPI::allocateArray
 (
     int**& array,
     int initVal,
@@ -182,7 +182,7 @@ void Foam::twoWayMPI::allocateArray
     allocate_external_int(array, width, length, initVal, lmp);
 }
 
-void Foam::twoWayMPI::destroy(int** array,int /*len*/) const
+void twoWayMPI::destroy(int** array,int /*len*/) const
 {
     if (array == NULL) return;
 
@@ -192,21 +192,19 @@ void Foam::twoWayMPI::destroy(int** array,int /*len*/) const
 }
 //============
 // int *
-void Foam::twoWayMPI::destroy(int* array) const
+void twoWayMPI::destroy(int* array) const
 {
-    if (array == NULL) return;
     free(array);
 }
 //============
 // double *
-void Foam::twoWayMPI::destroy(double* array) const
+void twoWayMPI::destroy(double* array) const
 {
-    if (array == NULL) return;
     free(array);
 }
 //============
 
-bool Foam::twoWayMPI::couple(int i) const
+bool twoWayMPI::couple(int i) const
 {
     bool coupleNow = false;
     if (i==0)
@@ -354,12 +352,12 @@ bool Foam::twoWayMPI::couple(int i) const
     return coupleNow;
 }
 
-int Foam::twoWayMPI::getNumberOfParticles() const
+int twoWayMPI::getNumberOfParticles() const
 {
     return liggghts_get_maxtag(lmp);
 }
 
-int Foam::twoWayMPI::getNumberOfClumps() const
+int twoWayMPI::getNumberOfClumps() const
 {
     #ifdef multisphere
         return liggghts_get_maxtag_ms(lmp);
@@ -369,7 +367,7 @@ int Foam::twoWayMPI::getNumberOfClumps() const
     return -1;
 }
 
-int Foam::twoWayMPI::getNumberOfTypes() const
+int twoWayMPI::getNumberOfTypes() const
 {
     #ifdef multisphere
         return liggghts_get_ntypes_ms(lmp);
@@ -378,7 +376,7 @@ int Foam::twoWayMPI::getNumberOfTypes() const
     return -1;
 }
 
-double* Foam::twoWayMPI::getTypeVol() const
+double* twoWayMPI::getTypeVol() const
 {
     #ifdef multisphere
         return liggghts_get_vclump_ms(lmp);

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.H

@@ -160,7 +160,7 @@ public:
 
         word myType() const { return typeName; }
 
-        void setCG() const { particleCloud_.setCG(lmp->force->cg()); }
+        void setCG() { particleCloud_.setCG(lmp->force->cg()); }
 };
 
 

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/twoWayMany2Many.C

@@ -62,7 +62,7 @@ twoWayMany2Many::twoWayMany2Many
 :
     dataExchangeModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
-    pbm_(sm.mesh().boundaryMesh()), 
+    pbm_(sm.mesh().boundaryMesh()),
     pData_(sm.mesh().globalData()),
     procPatches_(pData_.processorPatches()),
     procPatchIndices_(pData_.processorPatchIndices()),
@@ -280,7 +280,7 @@ void twoWayMany2Many::giveData
 
 //============
 // double **
-void Foam::twoWayMany2Many::allocateArray
+void twoWayMany2Many::allocateArray
 (
     double**& array,
     double initVal,
@@ -295,7 +295,7 @@ void Foam::twoWayMany2Many::allocateArray
             array[i][j] = initVal;
 }
 
-void Foam::twoWayMany2Many::allocateArray
+void twoWayMany2Many::allocateArray
 (
     double**& array,
     double initVal,
@@ -310,14 +310,14 @@ void Foam::twoWayMany2Many::allocateArray
             array[i][j] = initVal;
 }
 
-void inline Foam::twoWayMany2Many::destroy(double** array,int len) const
+void inline twoWayMany2Many::destroy(double** array,int len) const
 {
     lmp->memory->destroy(array);
 }
 
 //============
 // int **
-void Foam::twoWayMany2Many::allocateArray
+void twoWayMany2Many::allocateArray
 (
     int**& array,
     int initVal,
@@ -332,7 +332,7 @@ void Foam::twoWayMany2Many::allocateArray
             array[i][j] = initVal;
 }
 
-void Foam::twoWayMany2Many::allocateArray
+void twoWayMany2Many::allocateArray
 (
     int**& array,
     int initVal,
@@ -347,14 +347,14 @@ void Foam::twoWayMany2Many::allocateArray
             array[i][j] = initVal;
 }
 
-void inline Foam::twoWayMany2Many::destroy(int** array,int len) const
+void inline twoWayMany2Many::destroy(int** array,int len) const
 {
     lmp->memory->destroy(array);
 }
 
 //============
 // double *
-void Foam::twoWayMany2Many::allocateArray(double*& array, double initVal, int length) const
+void twoWayMany2Many::allocateArray(double*& array, double initVal, int length) const
 {
     int len = max(length,1);
     lmp->memory->grow(array, len, "m2m:dbl*");
@@ -362,14 +362,14 @@ void Foam::twoWayMany2Many::allocateArray(double*& array, double initVal, int le
         array[i] = initVal;
 }
 
-void inline Foam::twoWayMany2Many::destroy(double* array) const
+void inline twoWayMany2Many::destroy(double* array) const
 {
     lmp->memory->destroy(array);
 }
 
 //==============
 // int *
-void Foam::twoWayMany2Many::allocateArray(int*& array, int initVal, int length) const
+void twoWayMany2Many::allocateArray(int*& array, int initVal, int length) const
 {
     int len = max(length,1);
     lmp->memory->grow(array, len, "m2m:int*");
@@ -377,14 +377,14 @@ void Foam::twoWayMany2Many::allocateArray(int*& array, int initVal, int length)
         array[i] = initVal;
 }
 
-void inline Foam::twoWayMany2Many::destroy(int* array) const
+void inline twoWayMany2Many::destroy(int* array) const
 {
     lmp->memory->destroy(array);
 }
 //==============
 
 
-bool Foam::twoWayMany2Many::couple(int i) const
+bool twoWayMany2Many::couple(int i) const
 {
     bool coupleNow = false;
     if (i==0)
@@ -443,19 +443,19 @@ bool Foam::twoWayMany2Many::couple(int i) const
     return coupleNow;
 }
 
-int Foam::twoWayMany2Many::getNumberOfParticles() const
+int twoWayMany2Many::getNumberOfParticles() const
 {
     return liggghts_get_maxtag(lmp);
 }
 
-int Foam::twoWayMany2Many::getNumberOfClumps() const
+int twoWayMany2Many::getNumberOfClumps() const
 {
     Warning << "Foam::twoWayMany2Many::getNumberOfClumps() - changes necessary here" << endl;
     //return liggghts_get_maxtag_ms(lmp);
     return 1;
 }
 
-void Foam::twoWayMany2Many::syncIDs() const
+void twoWayMany2Many::syncIDs() const
 {
     particleCloud_.clockM().start(5,"recv_DEM_ids");
 
@@ -575,7 +575,7 @@ void Foam::twoWayMany2Many::syncIDs() const
     particleCloud_.clockM().stop("setup_Comm");
 }
 
-void Foam::twoWayMany2Many::locateParticle(int* id_lammpsSync, bool id_lammps_alloc_flag) const
+void twoWayMany2Many::locateParticle(int* id_lammpsSync, bool id_lammps_alloc_flag) const
 {
 #if defined(version21)
 
@@ -669,7 +669,7 @@ void Foam::twoWayMany2Many::locateParticle(int* id_lammpsSync, bool id_lammps_al
     }
     particleCloud_.clockM().stop("locate_Stage1");
     particleCloud_.clockM().start(8,"locate_Stage2");
- 
+
     PstreamBuffers pBufs(Pstream::nonBlocking);
 
     forAll(particleTransferID, i)
@@ -802,7 +802,7 @@ void Foam::twoWayMany2Many::locateParticle(int* id_lammpsSync, bool id_lammps_al
     if (firstRun_)
     {
         int* id_foam_nowhere_all;
-        Foam::dataExchangeModel::allocateArray(id_foam_nowhere_all,1,nlocal_foam_lostAll);
+        dataExchangeModel::allocateArray(id_foam_nowhere_all,1,nlocal_foam_lostAll);
         MPI_Allreduce(id_foamLostAll, id_foam_nowhere_all, nlocal_foam_lostAll, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
 
         int i = 0;
@@ -820,7 +820,7 @@ void Foam::twoWayMany2Many::locateParticle(int* id_lammpsSync, bool id_lammps_al
 
                         for (int k=0;k<3;k++)
                             id_lammpsVec_[j*3+k] = id_lammpsVec_[(nlocal_lammps_-1)*3+k];
-              
+
                         nlocal_lammps_ -= 1;
                         break;
                     }
@@ -828,7 +828,7 @@ void Foam::twoWayMany2Many::locateParticle(int* id_lammpsSync, bool id_lammps_al
             }
             i++;
         }
-        Foam::dataExchangeModel::destroy(id_foam_nowhere_all);
+        dataExchangeModel::destroy(id_foam_nowhere_all);
         id_foam_nowhere_all = NULL;
         if (id_lammps_alloc_flag) destroy(id_lammps_);
         id_lammps_ = NULL;

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/twoWayMany2Many.H

@@ -114,7 +114,7 @@ private:
     mutable int *cellID_foam_;
     mutable double *pos_foam_;
 
-    const polyBoundaryMesh& pbm_; 
+    const polyBoundaryMesh& pbm_;
     const globalMeshData& pData_;
     const labelList& procPatches_;
     const labelList& procPatchIndices_;
@@ -209,7 +209,7 @@ public:
         void syncIDs() const;
         void locateParticle(int*, bool) const;
         word myType() const { return typeName; }
-        void setCG() const { particleCloud_.setCG(lmp->force->cg()); }
+        void setCG() { particleCloud_.setCG(lmp->force->cg()); }
 };
 
 

src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C

@@ -134,19 +134,19 @@ heatTransferGunn::heatTransferGunn
     if (calcPartTempField_)
     {
         if (propsDict_.found("partRefTemp"))
-	    partRefTemp_.value()=readScalar(propsDict_.lookup ("partRefTemp"));
-	partTempField_.writeOpt() = IOobject::AUTO_WRITE;
-	partRelTempField_.writeOpt() = IOobject::AUTO_WRITE;
-	partTempField_.write();
-	partRelTempField_.write();
-	Info <<  "Particle temperature field activated." << endl;
+            partRefTemp_.value()=readScalar(propsDict_.lookup ("partRefTemp"));
+        partTempField_.writeOpt() = IOobject::AUTO_WRITE;
+        partRelTempField_.writeOpt() = IOobject::AUTO_WRITE;
+        partTempField_.write();
+        partRelTempField_.write();
+        Info <<  "Particle temperature field activated." << endl;
     }
     if (verbose_)
     {
         ReField_.writeOpt() = IOobject::AUTO_WRITE;
-	NuField_.writeOpt() = IOobject::AUTO_WRITE;
-	ReField_.write();
-	NuField_.write();
+        NuField_.writeOpt() = IOobject::AUTO_WRITE;
+        ReField_.write();
+        NuField_.write();
     }
 }
 
@@ -155,10 +155,10 @@ heatTransferGunn::heatTransferGunn
 
 heatTransferGunn::~heatTransferGunn()
 {
-    delete partTemp_;
-    delete partHeatFlux_;
-    delete partRe_;
-    delete partNu_;
+    particleCloud_.dataExchangeM().destroy(partTemp_,1);
+    particleCloud_.dataExchangeM().destroy(partHeatFlux_,1);
+    particleCloud_.dataExchangeM().destroy(partRe_,1);
+    particleCloud_.dataExchangeM().destroy(partNu_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
@@ -168,7 +168,7 @@ void heatTransferGunn::allocateMyArrays() const
     double initVal=0.0;
     particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1);  // field/initVal/with/lenghtFromLigghts
     particleCloud_.dataExchangeM().allocateArray(partHeatFlux_,initVal,1);
-    
+
     if(verbose_)
     {
         particleCloud_.dataExchangeM().allocateArray(partRe_,initVal,1);
@@ -188,10 +188,10 @@ void heatTransferGunn::calcEnergyContribution()
 
     // get DEM data
     particleCloud_.dataExchangeM().getData(partTempName_,"scalar-atom",partTemp_);
-    
+
     if(calcPartTempField_)
-    {       
-	partTempField_.primitiveFieldRef() = 0.0;
+    {
+        partTempField_.primitiveFieldRef() = 0.0;
         particleCloud_.averagingM().resetWeightFields();
         particleCloud_.averagingM().setScalarAverage
         (
@@ -202,10 +202,10 @@ void heatTransferGunn::calcEnergyContribution()
             NULL
         );
 
-	volScalarField sumTp (particleCloud_.averagingM().UsWeightField() * partTempField_);
-	dimensionedScalar aveTemp("aveTemp",dimensionSet(0,0,0,1,0,0,0), gSum(sumTp) / particleCloud_.numberOfParticles());
-	partRelTempField_ = (partTempField_ - aveTemp) / (aveTemp - partRefTemp_);
-	Info << "heatTransferGunn: average part. temp = " << aveTemp.value() << endl;
+        volScalarField sumTp (particleCloud_.averagingM().UsWeightField() * partTempField_);
+        dimensionedScalar aveTemp("aveTemp",dimensionSet(0,0,0,1,0,0,0), gSum(sumTp) / particleCloud_.numberOfParticles());
+        partRelTempField_ = (partTempField_ - aveTemp) / (aveTemp - partRefTemp_);
+        Info << "heatTransferGunn: average part. temp = " << aveTemp.value() << endl;
     }
 
     #ifdef compre
@@ -213,7 +213,7 @@ void heatTransferGunn::calcEnergyContribution()
     #else
        const volScalarField mufField = particleCloud_.turbulence().nu()*rho_;
     #endif
-    
+
 
     // calc La based heat flux
     scalar voidfraction(1);
@@ -251,9 +251,9 @@ void heatTransferGunn::calcEnergyContribution()
                     Ufluid = U_[cellI];
                     Tfluid = tempField_[cellI];
                 }
-                
+
                 if (voidfraction < 0.01)
-		    voidfraction = 0.01;
+                    voidfraction = 0.01;
 
                 // calc relative velocity
                 Us = particleCloud_.velocity(index);
@@ -264,7 +264,7 @@ void heatTransferGunn::calcEnergyContribution()
                 Pr = max(SMALL, Cp_ * muf / kf0_);
 
                 Nup = Nusselt(voidfraction, Rep, Pr);
-                
+
 
                 scalar h = kf0_ * Nup / ds;
                 scalar As = ds * ds * M_PI; // surface area of sphere
@@ -272,14 +272,14 @@ void heatTransferGunn::calcEnergyContribution()
                 // calc convective heat flux [W]
                 heatFlux(index, h, As, Tfluid);
                 heatFluxCoeff(index, h, As);
-		
-		if(verbose_)
-		{
-		    partRe_[index][0] = Rep;
-		    partNu_[index][0] = Nup;
-		}
-		
-		if(particleCloud_.verbose() && index >=0 && index <2)
+
+                if(verbose_)
+                {
+                    partRe_[index][0] = Rep;
+                    partNu_[index][0] = Nup;
+                }
+
+                if(particleCloud_.verbose() && index >=0 && index <2)
                 {
                     Info << "partHeatFlux = " << partHeatFlux_[index][0] << endl;
                     Info << "magUr = " << magUr << endl;
@@ -308,27 +308,27 @@ void heatTransferGunn::calcEnergyContribution()
     if(verbose_)
     {
         ReField_.primitiveFieldRef() = 0.0;
-	NuField_.primitiveFieldRef() = 0.0;
-	particleCloud_.averagingM().resetWeightFields();
+        NuField_.primitiveFieldRef() = 0.0;
+        particleCloud_.averagingM().resetWeightFields();
         particleCloud_.averagingM().setScalarAverage
         (
             ReField_,
             partRe_,
             particleCloud_.particleWeights(),
-	    particleCloud_.averagingM().UsWeightField(),
+            particleCloud_.averagingM().UsWeightField(),
             NULL
         );
-	particleCloud_.averagingM().resetWeightFields();
+        particleCloud_.averagingM().resetWeightFields();
         particleCloud_.averagingM().setScalarAverage
         (
             NuField_,
             partNu_,
             particleCloud_.particleWeights(),
-	    particleCloud_.averagingM().UsWeightField(),
+            particleCloud_.averagingM().UsWeightField(),
             NULL
         );
     }
-    
+
     // limit source term
     forAll(QPartFluid_,cellI)
     {
@@ -336,13 +336,13 @@ void heatTransferGunn::calcEnergyContribution()
 
         if(mag(EuFieldInCell) > maxSource_ )
         {
-	     Pout << "limiting source term\n"  << endl  ;
+             Pout << "limiting source term\n"  << endl  ;
              QPartFluid_[cellI] = sign(EuFieldInCell) * maxSource_;
         }
     }
-    
+
     QPartFluid_.correctBoundaryConditions();
-    
+
     giveData(0);
 
 }
@@ -354,13 +354,10 @@ void heatTransferGunn::addEnergyContribution(volScalarField& Qsource) const
 
 scalar heatTransferGunn::Nusselt(scalar voidfraction, scalar Rep, scalar Pr) const
 {
-    scalar Nup(0.0);
-    Nup = (7 - 10 * voidfraction + 5 * voidfraction * voidfraction) *
+    return (7 - 10 * voidfraction + 5 * voidfraction * voidfraction) *
                         (1 + 0.7 * Foam::pow(Rep,0.2) * Foam::pow(Pr,0.33)) +
-			(1.33 - 2.4 * voidfraction + 1.2 * voidfraction * voidfraction) *
-			Foam::pow(Rep,0.7) * Foam::pow(Pr,0.33);
-                        
-    return Nup;
+                        (1.33 - 2.4 * voidfraction + 1.2 * voidfraction * voidfraction) *
+                        Foam::pow(Rep,0.7) * Foam::pow(Pr,0.33);
 }
 
 void heatTransferGunn::heatFlux(label index, scalar h, scalar As, scalar Tfluid)
@@ -370,7 +367,7 @@ void heatTransferGunn::heatFlux(label index, scalar h, scalar As, scalar Tfluid)
 
 void heatTransferGunn::heatFluxCoeff(label index, scalar h, scalar As)
 {
-   //no heat transfer coefficient in explicit model   
+    //no heat transfer coefficient in explicit model
 }
 
 void heatTransferGunn::giveData(int call)
@@ -379,7 +376,7 @@ void heatTransferGunn::giveData(int call)
     {
         Info << "total convective particle-fluid heat flux [W] (Eulerian) = " << gSum(QPartFluid_*1.0*QPartFluid_.mesh().V()) << endl;
 
-        particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);    
+        particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);
     }
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnImplicit/heatTransferGunnImplicit.C

@@ -58,7 +58,7 @@ heatTransferGunnImplicit::heatTransferGunnImplicit
     partHeatFluxCoeff_(NULL)
 {
     allocateMyArrays();
-    
+
     // no limiting necessary for implicit heat transfer
     maxSource_ = 1e30;
 }
@@ -68,7 +68,7 @@ heatTransferGunnImplicit::heatTransferGunnImplicit
 
 heatTransferGunnImplicit::~heatTransferGunnImplicit()
 {
-    delete partHeatFluxCoeff_;
+    particleCloud_.dataExchangeM().destroy(partHeatFluxCoeff_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
@@ -87,7 +87,7 @@ void heatTransferGunnImplicit::calcEnergyContribution()
     heatTransferGunn::calcEnergyContribution();
 
     QPartFluidCoeff_.primitiveFieldRef() = 0.0;
-    
+
     particleCloud_.averagingM().setScalarSum
     (
         QPartFluidCoeff_,
@@ -97,7 +97,7 @@ void heatTransferGunnImplicit::calcEnergyContribution()
     );
 
     QPartFluidCoeff_.primitiveFieldRef() /= -QPartFluidCoeff_.mesh().V();
-    
+
 //    QPartFluidCoeff_.correctBoundaryConditions();
 
 }
@@ -123,7 +123,7 @@ void heatTransferGunnImplicit::giveData(int call)
     {
         //Info << "total convective particle-fluid heat flux [W] (Eulerian) = " << gSum(QPartFluid_*1.0*QPartFluid_.mesh().V()) << endl;
 
-        particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);    
+        particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);
     }
 }
 
@@ -133,7 +133,7 @@ void heatTransferGunnImplicit::postFlow()
     scalar Tfluid(0.0);
     scalar Tpart(0.0);
     interpolationCellPoint<scalar> TInterpolator_(tempField_);
-        
+
     for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
     {
             cellI = particleCloud_.cellIDs()[index][0];
@@ -146,12 +146,12 @@ void heatTransferGunnImplicit::postFlow()
                 }
                 else
                     Tfluid = tempField_[cellI];
-                
+
                 Tpart = partTemp_[index][0];
                 partHeatFlux_[index][0] = (Tfluid - Tpart) * partHeatFluxCoeff_[index][0];
             }
     }
-                
+
     giveData(1);
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C

@@ -77,7 +77,7 @@ reactionHeat::reactionHeat
 
 reactionHeat::~reactionHeat()
 {
-    delete reactionHeat_;
+    particleCloud_.dataExchangeM().destroy(reactionHeat_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //

src/lagrangian/cfdemParticle/subModels/forceModel/Archimedes/Archimedes.C

@@ -68,10 +68,10 @@ Archimedes::Archimedes
 {
 
     //Append the field names to be probed
-    particleCloud_.probeM().initialize(typeName, "archimedesF.logDat");
+    particleCloud_.probeM().initialize(typeName, typeName+".logDat");
     particleCloud_.probeM().vectorFields_.append("archimedesForce");  //first entry must the be the force
     particleCloud_.probeM().scalarFields_.append("Vp");
-    particleCloud_.probeM().writeHeader();  
+    particleCloud_.probeM().writeHeader();
 
 
     if (propsDict_.found("twoDimensional"))
@@ -84,18 +84,20 @@ Archimedes::Archimedes
     setForceSubModels(propsDict_);
 
     // define switches which can be read from dict
-    forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(1,true); // activate treatForceDEM switch
+    forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
+    forceSubM(0).setSwitchesList(SW_TREAT_FORCE_DEM,true); // activate treatForceDEM switch
     forceSubM(0).readSwitches();
 
-    if (modelType_=="A" || modelType_=="Bfull"){
-        if(!forceSubM(0).switches()[1]) // treatDEM != true
+    if (modelType_=="A" || modelType_=="Bfull")
+    {
+        if(!forceSubM(0).switches()[SW_TREAT_FORCE_DEM]) // treatDEM != true
         {
             Warning << "Usually model type A and Bfull need Archimedes only on DEM side only (treatForceDEM=true)! are you sure about your settings?" << endl;
         }
     }
-    if (modelType_=="B"){
-        if(forceSubM(0).switches()[1]) // treatDEM = true
+    else if (modelType_=="B")
+    {
+        if(forceSubM(0).switches()[SW_TREAT_FORCE_DEM]) // treatDEM = true
         {
             Warning << "Usually model type B needs Archimedes on CFD and DEM side (treatForceDEM=false)! are you sure about your settings?" << endl;
         }
@@ -119,12 +121,12 @@ void Archimedes::setForce() const
 
     #include "setupProbeModel.H"
 
-    for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
+    for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
         //if(mask[index][0])
         //{
             label cellI = particleCloud_.cellIDs()[index][0];
-            force=vector::zero;
+            force = vector::zero;
 
             if (cellI > -1) // particle Found
             {
@@ -133,7 +135,9 @@ void Archimedes::setForce() const
                     scalar r = particleCloud_.radius(index);
                     force = -g_.value()*forceSubM(0).rhoField()[cellI]*r*r*M_PI; // circle area
                     Warning << "Archimedes::setForce() : this functionality is not tested!" << endl;
-                }else{
+                }
+                else
+                {
                     force = -g_.value()*forceSubM(0).rhoField()[cellI]*particleCloud_.particleVolume(index);
                 }
 

src/lagrangian/cfdemParticle/subModels/forceModel/ArchimedesIB/ArchimedesIB.C

@@ -69,9 +69,9 @@ ArchimedesIB::ArchimedesIB
     g_(sm.mesh().lookupObject<uniformDimensionedVectorField> (gravityFieldName_))
 {
     //Append the field names to be probed
-    particleCloud_.probeM().initialize(typeName, "archimedesIBF.logDat");
+    particleCloud_.probeM().initialize(typeName, typeName+".logDat");
     particleCloud_.probeM().vectorFields_.append("archimedesIBForce");  //first entry must the be the force
-    particleCloud_.probeM().writeHeader();  
+    particleCloud_.probeM().writeHeader();
 
     if (propsDict_.found("twoDimensional"))
     {
@@ -83,12 +83,12 @@ ArchimedesIB::ArchimedesIB
     setForceSubModels(propsDict_);
 
     // define switches which can be read from dict
-    forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch
+    forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
 
     // read those switches defined above, if provided in dict
     forceSubM(0).readSwitches();
 
-    forceSubM(0).setSwitches(1,true); // treatDEM = true
+    forceSubM(0).setSwitches(SW_TREAT_FORCE_DEM,true); // treatDEM = true
     Info << "accounting for Archimedes only on DEM side!" << endl;
 
     particleCloud_.checkCG(true);
@@ -109,7 +109,7 @@ void ArchimedesIB::setForce() const
 
     #include "setupProbeModel.H"
 
-    for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
+    for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
         //if(mask[index][0])
         //{
@@ -120,8 +120,8 @@ void ArchimedesIB::setForce() const
                 if (cellI > -1) // particle Found
                 {
                     //force += -g_.value()*forceSubM(0).rhoField()[cellI]*forceSubM(0).rhoField().mesh().V()[cellI]*(1-particleCloud_.voidfractions()[index][subCell]);//mod by alice
-                	force += -g_.value()*forceSubM(0).rhoField()[cellI]*particleCloud_.mesh().V()[cellI]*(1-voidfractions_[cellI]);//mod by alice
-        	    }
+                    force += -g_.value()*forceSubM(0).rhoField()[cellI]*particleCloud_.mesh().V()[cellI]*(1-voidfractions_[cellI]);//mod by alice
+                }
             }
 
             //Set value fields and write the probe

src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C

@@ -59,7 +59,7 @@ BeetstraDrag::BeetstraDrag
     scaleDrag_(1.)
 {
     //Append the field names to be probed
-    particleCloud_.probeM().initialize(typeName, "BeetstraDrag.logDat");
+    particleCloud_.probeM().initialize(typeName, typeName+".logDat");
     particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must  be the force
     particleCloud_.probeM().vectorFields_.append("Urel");
     particleCloud_.probeM().scalarFields_.append("Rep");
@@ -70,11 +70,11 @@ BeetstraDrag::BeetstraDrag
     // init force sub model
     setForceSubModels(propsDict_);
     // define switches which can be read from dict
-    forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(2,true); // activate implDEM switch
-    forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch
-    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
-    forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity switch
+    forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
+    forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM,true); // activate implDEM switch
+    forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(SW_SCALAR_VISCOSITY,true); // activate scalarViscosity switch
     forceSubM(0).readSwitches();
 
     particleCloud_.checkCG(true);
@@ -97,8 +97,11 @@ BeetstraDrag::~BeetstraDrag()
 void BeetstraDrag::setForce() const
 {
     if (scaleDia_ > 1)
+    {
         Info << "Beetstra using scale = " << scaleDia_ << endl;
-    else if (particleCloud_.cg() > 1){
+    }
+    else if (particleCloud_.cg() > 1)
+    {
         scaleDia_=particleCloud_.cg();
         Info << "Beetstra using scale from liggghts cg = " << scaleDia_ << endl;
     }
@@ -125,13 +128,13 @@ void BeetstraDrag::setForce() const
 
     vector dragExplicit(0,0,0);
     scalar dragCoefficient(0);
-    
+
     interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
     interpolationCellPoint<vector> UInterpolator_(U_);
 
     #include "setupProbeModel.H"
 
-    for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
+    for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
             cellI = particleCloud_.cellIDs()[index][0];
             drag = vector(0,0,0);
@@ -163,7 +166,7 @@ void BeetstraDrag::setForce() const
                 Ur = Ufluid-Us;
                 magUr = mag(Ur);
                 ds = 2*particleCloud_.radius(index);
-		ds_scaled = ds/scaleDia_;
+                ds_scaled = ds/scaleDia_;
                 rho = rhoField[cellI];
                 nuf = nufField[cellI];
 
@@ -171,7 +174,7 @@ void BeetstraDrag::setForce() const
                 localPhiP = 1.0f-voidfraction+SMALL;
 
                 // calc particle's drag coefficient (i.e., Force per unit slip velocity and Stokes drag)
-               
+
                 Rep=ds_scaled*voidfraction*magUr/nuf+SMALL;
                 dragCoefficient = 10.0*localPhiP/(voidfraction*voidfraction) +
                                   voidfraction*voidfraction*(1.0+1.5*Foam::sqrt(localPhiP)) +

src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H

@@ -51,7 +51,7 @@ private:
 
     const volScalarField& voidfraction_;
 
-    word UsFieldName_;   
+    word UsFieldName_;
 
     const volVectorField& UsField_;
 

src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.C

@@ -34,8 +34,6 @@ Description
 #include "DiFeliceDrag.H"
 #include "addToRunTimeSelectionTable.H"
 
-//#include <mpi.h>
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -74,12 +72,12 @@ DiFeliceDrag::DiFeliceDrag
     scaleDrag_(1.)
 {
     //Append the field names to be probed
-    particleCloud_.probeM().initialize(typeName, "diFeliceDrag.logDat");
-    particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must the be the force
-    particleCloud_.probeM().vectorFields_.append("Urel");        //other are debug
-    particleCloud_.probeM().scalarFields_.append("Rep");          //other are debug
-    particleCloud_.probeM().scalarFields_.append("Cd");                 //other are debug
-    particleCloud_.probeM().scalarFields_.append("voidfraction");       //other are debug
+    particleCloud_.probeM().initialize(typeName, typeName+".logDat");
+    particleCloud_.probeM().vectorFields_.append("dragForce");    // first entry must the be the force
+    particleCloud_.probeM().vectorFields_.append("Urel");         // other are debug
+    particleCloud_.probeM().scalarFields_.append("Rep");          // other are debug
+    particleCloud_.probeM().scalarFields_.append("Cd");           // other are debug
+    particleCloud_.probeM().scalarFields_.append("voidfraction"); // other are debug
     particleCloud_.probeM().writeHeader();
 
     particleCloud_.checkCG(true);
@@ -92,11 +90,11 @@ DiFeliceDrag::DiFeliceDrag
     setForceSubModels(propsDict_);
 
     // define switches which can be read from dict
-    forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(2,true); // activate implDEM switch
-    forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch
-    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
-    forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity switch
+    forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
+    forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM,true); // activate implDEM switch
+    forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(SW_SCALAR_VISCOSITY,true); // activate scalarViscosity switch
 
     // read those switches defined above, if provided in dict
     forceSubM(0).readSwitches();
@@ -114,15 +112,18 @@ DiFeliceDrag::~DiFeliceDrag()
 void DiFeliceDrag::setForce() const
 {
     if (scaleDia_ > 1)
+    {
         Info << "DiFeliceDrag using scale = " << scaleDia_ << endl;
-    else if (particleCloud_.cg() > 1){
+    }
+    else if (particleCloud_.cg() > 1)
+    {
         scaleDia_=particleCloud_.cg();
         Info << "DiFeliceDrag using scale from liggghts cg = " << scaleDia_ << endl;
     }
 
     const volScalarField& nufField = forceSubM(0).nuField();
     const volScalarField& rhoField = forceSubM(0).rhoField();
-    
+
     vector position(0,0,0);
     scalar voidfraction(1);
     vector Ufluid(0,0,0);
@@ -144,7 +145,7 @@ void DiFeliceDrag::setForce() const
 
     #include "setupProbeModel.H"
 
-    for(int index = 0;index <  particleCloud_.numberOfParticles(); index++)
+    for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
             cellI = particleCloud_.cellIDs()[index][0];
             drag = vector(0,0,0);
@@ -159,7 +160,8 @@ void DiFeliceDrag::setForce() const
                     position = particleCloud_.position(index);
                     voidfraction = voidfractionInterpolator_.interpolate(position,cellI);
                     Ufluid = UInterpolator_.interpolate(position,cellI);
-                }else
+                }
+                else
                 {
                     voidfraction = voidfraction_[cellI];
                     Ufluid = U_[cellI];
@@ -189,8 +191,8 @@ void DiFeliceDrag::setForce() const
                     // calc particle's drag
                     dragCoefficient = 0.125*Cd*rho
                                      *M_PI
-                                     *ds*ds     
-                                     *scaleDia_ 
+                                     *ds*ds
+                                     *scaleDia_
                                      *pow(voidfraction,(2-Xi))*magUr
                                      *scaleDrag_;
                     if (modelType_=="B")

src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.H

@@ -68,7 +68,7 @@ private:
 
     const volScalarField& voidfraction_;
 
-    word UsFieldName_;   
+    word UsFieldName_;
 
     const volVectorField& UsField_;  // the average particle velocity field (for implicit/expliti force split)
 

src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.C

@@ -66,7 +66,7 @@ ErgunStatFines::ErgunStatFines
     switchingVoidfraction_(0.8)
 {
     //Append the field names to be probed
-    particleCloud_.probeM().initialize(typeName, "ErgunStatFines.logDat");
+    particleCloud_.probeM().initialize(typeName, typeName+".logDat");
     particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must  be the force
     particleCloud_.probeM().vectorFields_.append("Urel");
     particleCloud_.probeM().scalarFields_.append("Rep");
@@ -77,26 +77,26 @@ ErgunStatFines::ErgunStatFines
     // init force sub model
     setForceSubModels(propsDict_);
     // define switches which can be read from dict
-    forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(2,true); // activate implDEM switch
-    forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch
-    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
-    forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity switch
+    forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
+    forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM,true); // activate implDEM switch
+    forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(SW_SCALAR_VISCOSITY,true); // activate scalarViscosity switch
     forceSubM(0).readSwitches();
 
     particleCloud_.checkCG(true);
     if (propsDict_.found("scale"))
-        scaleDia_=scalar(readScalar(propsDict_.lookup("scale")));
+        scaleDia_ = scalar(readScalar(propsDict_.lookup("scale")));
     if (propsDict_.found("scaleDrag"))
-        scaleDrag_=scalar(readScalar(propsDict_.lookup("scaleDrag")));
+        scaleDrag_ = scalar(readScalar(propsDict_.lookup("scaleDrag")));
 
     if (propsDict_.found("switchingVoidfraction"))
         switchingVoidfraction_ = readScalar(propsDict_.lookup("switchingVoidfraction"));
 
     dictionary SauterDict(dict.subDict("dSauterProps"));
     if (SauterDict.found("scaleDist"))
-        scaleDist_=scalar(readScalar(SauterDict.lookup("scaleDist")));
-    
+        scaleDist_ = scalar(readScalar(SauterDict.lookup("scaleDist")));
+
 }
 
 
@@ -117,7 +117,9 @@ scalar ErgunStatFines::dSauter(label cellI) const
 void ErgunStatFines::setForce() const
 {
     if (scaleDia_ > 1)
+    {
         Info << "ErgunStatFines using scale = " << scaleDia_ << endl;
+    }
     else if (particleCloud_.cg() > 1)
     {
         scaleDia_=particleCloud_.cg();
@@ -131,7 +133,7 @@ void ErgunStatFines::setForce() const
     scalar voidfraction(1);
     vector Ufluid(0,0,0);
     vector drag(0,0,0);
-    label cellI=0;
+    label cellI = 0;
 
     vector Us(0,0,0);
     vector Ur(0,0,0);
@@ -144,13 +146,13 @@ void ErgunStatFines::setForce() const
     scalar alphaPartEff(0);
 
     scalar CdMagUrLag(0);       //Cd of the very particle
-    scalar betaP(0);             //momentum exchange of the very particle
+    scalar betaP(0);            //momentum exchange of the very particle
 
     vector dragExplicit(0,0,0);
     scalar dragCoefficient(0);
 
     scalar scaleDia3 = scaleDia_*scaleDia_*scaleDia_;
-    
+
     interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
     interpolationCellPoint<vector> UInterpolator_(U_);
 
@@ -158,15 +160,15 @@ void ErgunStatFines::setForce() const
 
     if(forceSubM(0).verbose())
                 Info << "Entering force loop of ErgunStatFines.\n" << endl;
-    
-    for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
-    {	    
+
+    for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
+    {
             cellI = particleCloud_.cellIDs()[index][0];
             drag = vector(0,0,0);
             dragExplicit = vector(0,0,0);
             betaP = 0;
-            Ufluid =vector(0,0,0);
-            voidfraction=0;
+            Ufluid = vector(0,0,0);
+            voidfraction = 0;
             dragCoefficient = 0;
 
             if (cellI > -1) // particle found
@@ -174,34 +176,34 @@ void ErgunStatFines::setForce() const
 
                 if( forceSubM(0).interpolation() )
                 {
-	            position     = particleCloud_.position(index);
+                    position     = particleCloud_.position(index);
                     voidfraction = voidfractionInterpolator_.interpolate(position,cellI);
-                    Ufluid       = UInterpolator_.interpolate(position,cellI);      
+                    Ufluid       = UInterpolator_.interpolate(position,cellI);
                 }
                 else
                 {
-		    voidfraction = voidfraction_[cellI];
+                    voidfraction = voidfraction_[cellI];
                     Ufluid = U_[cellI];
                 }
 
                 // ensure voidfraction to be meaningful
                 // problems could arise from interpolation or empty cells
 
-                if(voidfraction>0.999) 
+                if(voidfraction > 0.999)
                     voidfraction = 0.999;
-                else if(voidfraction<0.05)
+                else if(voidfraction < 0.05)
                     voidfraction = 0.05;
 
                 Us = particleCloud_.velocity(index);
                 Ur = Ufluid-Us;
                 magUr = mag(Ur);
-		dSauterMix = dSauterMix_[cellI];
-		ds = 2*particleCloud_.radius(index);
+                dSauterMix = dSauterMix_[cellI];
+                ds = 2*particleCloud_.radius(index);
                 rho = rhoField[cellI];
                 nuf = nufField[cellI];
 
-                Rep=0.0;
-		alphaPartEff = 1.0 - voidfraction + alphaSt_[cellI] + SMALL;
+                Rep = 0.0;
+                alphaPartEff = 1.0 - voidfraction + alphaSt_[cellI] + SMALL;
 
                 // calc particle's drag coefficient (i.e., Force per unit slip velocity and per m³ PARTICLE)
                 if(voidfraction > switchingVoidfraction_) //dilute, no static hold-up present
@@ -210,7 +212,7 @@ void ErgunStatFines::setForce() const
                     CdMagUrLag = (24.0*nuf/(dSauterMix*voidfraction)) //1/magUr missing here, but compensated in expression for betaP!
                                  *(scalar(1.0)+0.15*Foam::pow(Rep, 0.687));
 
-                    betaP = 0.75* alphaPartEff * ( 
+                    betaP = 0.75* alphaPartEff * (
                                             rho*voidfraction*CdMagUrLag
                                           /
                                             (dSauterMix*Foam::pow(voidfraction,2.65))
@@ -224,12 +226,12 @@ void ErgunStatFines::setForce() const
                               (1.75 * magUr * rho * alphaPartEff)
                              /((dSauterMix*phi_));
                 }
-       
+
                 // calc particle's drag
                 betaP /= (1-alphaPartEff);
                 dragCoefficient = M_PI/6 * ds/scaleDia_ * ds/scaleDia_ * dSauter(cellI) * voidfraction / (1 - voidfraction) * betaP * scaleDrag_;
                 dragCoefficient *= scaleDia3;
-                if (modelType_=="B")
+                if (modelType_ == "B")
                     dragCoefficient /= voidfraction;
 
                 drag = dragCoefficient * Ur;
@@ -237,7 +239,7 @@ void ErgunStatFines::setForce() const
                 // explicitCorr
                 forceSubM(0).explicitCorr(drag,dragExplicit,dragCoefficient,Ufluid,U_[cellI],Us,UsField_[cellI],forceSubM(0).verbose());
 
-                if(forceSubM(0).verbose() && index >=0 && index <2)
+                if(forceSubM(0).verbose() && index >= 0 && index < 2)
                 {
                     Pout << "cellI = " << cellI << endl;
                     Pout << "index = " << index << endl;
@@ -271,7 +273,7 @@ void ErgunStatFines::setForce() const
             forceSubM(0).partToArray(index,drag,dragExplicit,Ufluid,dragCoefficient);
 
     }// end loop particles
-    
+
     if(forceSubM(0).verbose())
         Pout << "Leaving force loop of ErgunStatFines.\n" << endl;
 }

src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.H

@@ -50,18 +50,18 @@ private:
     word voidfractionFieldName_;
 
     const volScalarField& voidfraction_;
-    
+
     const volScalarField& dSauter_;
-    
+
     const volScalarField& dSauterMix_;
-    
+
     const volScalarField& alphaP_;
-    
+
     const volScalarField& alphaSt_;
 
     const scalar phi_;
 
-    word UsFieldName_;   
+    word UsFieldName_;
 
     const volVectorField& UsField_;  // the average particle velocity field
 

src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.C

@@ -65,16 +65,16 @@ FanningDynFines::FanningDynFines
     // init force sub model
     setForceSubModels(propsDict_);
     // define switches which can be read from dict
-    forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch
+    forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
+    forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
     forceSubM(0).readSwitches();
-    forceSubM(0).setSwitches(0,true);
+    forceSubM(0).setSwitches(SW_TREAT_FORCE_EXPLICIT,true);
 
     particleCloud_.checkCG(true);
     if (propsDict_.found("scale"))
-        scaleDia_=scalar(readScalar(propsDict_.lookup("scale")));
+        scaleDia_ = scalar(readScalar(propsDict_.lookup("scale")));
     if (propsDict_.found("scaleDrag"))
-        scaleDrag_=scalar(readScalar(propsDict_.lookup("scaleDrag")));
+        scaleDrag_ = scalar(readScalar(propsDict_.lookup("scaleDrag")));
 
 }
 
@@ -91,18 +91,20 @@ void FanningDynFines::setForce() const
 {
     if(forceSubM(0).verbose())
         Info << "Entering force loop of FanningDynFines.\n" << endl;
- 
+
     if (scaleDia_ > 1)
+    {
         Info << "FanningDynFines using scale = " << scaleDia_ << endl;
+    }
     else if (particleCloud_.cg() > 1)
     {
         scaleDia_=particleCloud_.cg();
         Info << "FanningDynFines using scale from liggghts cg = " << scaleDia_ << endl;
     }
-   
+
     vector UDyn(0,0,0);
     vector drag(0,0,0);
-    label cellI=0;
+    label cellI = 0;
 
     vector Us(0,0,0);
     vector Ur(0,0,0);
@@ -111,11 +113,11 @@ void FanningDynFines::setForce() const
     scalar scaleDia3 = scaleDia_*scaleDia_*scaleDia_;
 
     scalar dragCoefficient(0);
-    
+
 
     #include "setupProbeModel.H"
 
-    for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
+    for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
             cellI = particleCloud_.cellIDs()[index][0];
             drag = vector(0,0,0);
@@ -127,14 +129,14 @@ void FanningDynFines::setForce() const
                 UDyn = UDyn_[cellI];
                 Us = UsField_[cellI];
                 Ur = UDyn-Us;
-                ds = 2*particleCloud_.radius(index);
-		ds_scaled = ds/scaleDia_;
-		
+                ds = 2 * particleCloud_.radius(index);
+                ds_scaled = ds/scaleDia_;
+
                 dragCoefficient = FanningCoeff_[cellI];
 
                 // calc particle's drag
                 dragCoefficient *= M_PI/6 * ds_scaled * ds_scaled / alphaP_[cellI] * dSauter_[cellI] * scaleDia3 * scaleDrag_;
-                if (modelType_=="B")
+                if (modelType_ == "B")
                     dragCoefficient /= voidfraction_[cellI];
 
                 drag = dragCoefficient * Ur;
@@ -143,8 +145,8 @@ void FanningDynFines::setForce() const
             // write particle based data to global array
             forceSubM(0).partToArray(index,drag,vector::zero);
     }
-    
-    if(forceSubM(0).verbose())
+
+    if (forceSubM(0).verbose())
         Info << "Leaving force loop of FanningDynFines.\n" << endl;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.H

@@ -41,7 +41,7 @@ class FanningDynFines
     public forceModel
 {
 private:
-  
+
     dictionary propsDict_;
 
     word velFieldName_;
@@ -52,22 +52,22 @@ private:
 
     const volScalarField& voidfraction_;
 
-    word UsFieldName_;   
+    word UsFieldName_;
 
     const volVectorField& UsField_;
-    
+
     const volVectorField& UDyn_;
-    
+
     const volScalarField& FanningCoeff_;
-    
+
     const volScalarField& alphaP_;
-    
+
     const volScalarField& dSauter_;
 
     mutable scalar scaleDia_;
 
     mutable scalar scaleDrag_;
-    
+
 
 public:
 

src/lagrangian/cfdemParticle/subModels/forceModel/Fines/Fines.H

@@ -41,7 +41,7 @@ class Fines
     public forceModel
 {
 private:
-    mutable FinesFields finesFields_;    
+    mutable FinesFields finesFields_;
 
 public:
 

src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.C

@@ -56,7 +56,7 @@ FinesFields::FinesFields
     p_(sm.mesh().lookupObject<volScalarField> (pFieldName_)),
     rhoGFieldName_(propsDict_.lookupOrDefault<word>("rhoGFieldName","rho")),
     rhoG_(sm.mesh().lookupObject<volScalarField> (rhoGFieldName_)),
-    dSauter_(sm.mesh().lookupObject<volScalarField> ("dSauter")),   
+    dSauter_(sm.mesh().lookupObject<volScalarField> ("dSauter")),
     alphaG_
     (   IOobject
         (
@@ -126,7 +126,7 @@ FinesFields::FinesFields
         ),
         sm.mesh(),
         dimensionedScalar("zero", dimensionSet(0,1,0,0,0), 0),
-	"zeroGradient"
+        "zeroGradient"
     ),
     DragCoeff_
     (   IOobject
@@ -151,7 +151,7 @@ FinesFields::FinesFields
         ),
         sm.mesh(),
         dimensionedScalar("zero", dimensionSet(0,1,0,0,0), 0),
-	"zeroGradient"
+        "zeroGradient"
     ),
     FanningCoeff_
     (   IOobject
@@ -199,8 +199,8 @@ FinesFields::FinesFields
             IOobject::NO_WRITE
         ),
         sm.mesh(),
-	dimensionedScalar("zero", dimensionSet(1,0,-1,0,0), 0)
-	//dimensionedVector("zero", dimensionSet(1,-2,-1,0,0), vector::zero)
+        dimensionedScalar("zero", dimensionSet(1,0,-1,0,0), 0)
+        //dimensionedVector("zero", dimensionSet(1,-2,-1,0,0), vector::zero)
     ),
     uDyn_
     (   IOobject
@@ -243,12 +243,12 @@ FinesFields::FinesFields
     if (propsDict_.found("rhoFine"))
         rhoFine_.value()=readScalar(propsDict_.lookup ("rhoFine"));
     else
-        FatalError <<"Please specify rhoFine.\n" << abort(FatalError);  
+        FatalError <<"Please specify rhoFine.\n" << abort(FatalError);
     if (propsDict_.found("nuAve"))
         nuAve_.value()=readScalar(propsDict_.lookup ("nuAve"));
     if (propsDict_.found("alphaDynMax"))
         alphaDynMax_=readScalar(propsDict_.lookup ("alphaDynMax"));
-    
+
     if(verbose_)
     {
       alphaG_.writeOpt() = IOobject::AUTO_WRITE;
@@ -305,143 +305,144 @@ void FinesFields::update()
 }
 
 
-void FinesFields::calcSource() 
+void FinesFields::calcSource()
 {
-    Sds_.primitiveFieldRef()=0;
+    Sds_.primitiveFieldRef() = 0;
     deltaAlpha_.primitiveFieldRef() = 0.0;
     scalar f(0.0);
     scalar critpore(0.0);
     scalar dmean(0.0);
     scalar d1(0.0);
     scalar d2(0.0);
+
     forAll(Sds_,cellI)
     {
         // calculate everything in units auf dSauter
         critpore = nCrit_*dFine_.value()/dSauter_[cellI];
-	// pore size from hydraulic radius
-	dmean = 2 * (1 - alphaP_[cellI]) / ( (1 + poresizeWidth_*poresizeWidth_/3) * 3 * alphaP_[cellI] );
-	// Sweeney and Martin, Acta Materialia 51 (2003): ratio of hydraulic to pore throat radius
-	dmean /= ratioHydraulicPore_;
-	d1 = dmean * (1 - poresizeWidth_);
-	d2 = dmean * (1 + poresizeWidth_);
-        
+        // pore size from hydraulic radius
+        dmean = 2 * (1 - alphaP_[cellI]) / ( (1 + poresizeWidth_*poresizeWidth_/3) * 3 * alphaP_[cellI] );
+        // Sweeney and Martin, Acta Materialia 51 (2003): ratio of hydraulic to pore throat radius
+        dmean /= ratioHydraulicPore_;
+        d1 = dmean * (1 - poresizeWidth_);
+        d2 = dmean * (1 + poresizeWidth_);
+
         f = (critpore*critpore*critpore - d1 * d1 * d1) / (d2 * d2 * d2 - d1 * d1 * d1);
-	if (f<0)
-	{
-	    f=0.0;    
-	}
-	else if (f>1.0)
-	{
-	    f=1.0;
+        if (f < 0)
+        {
+            f = 0.0;
         }
-	
-	// at this point, voidfraction is still calculated from the true particle sizes
-	deltaAlpha_[cellI] = f * (alphaMax_ - alphaP_[cellI]) - alphaSt_[cellI];
-	// too much volume occupied: release it (50% per time step)
-	if (deltaAlpha_[cellI] < 0.0)
-	{
-	    Sds_[cellI] = 0.5*deltaAlpha_[cellI];
-	}
-	// volume too occupy available: deposit at most 80% of dyn hold up
-	else if (depRate_ * deltaAlpha_[cellI] > 0.8 * alphaDyn_[cellI])
-	{
-	    Sds_[cellI] = 0.8 * alphaDyn_[cellI];
-	}
-	else
-	{
-	    Sds_[cellI] = depRate_ * deltaAlpha_[cellI];
-	}
-    }   
+        else if (f > 1.0)
+        {
+            f = 1.0;
+        }
+
+        // at this point, voidfraction is still calculated from the true particle sizes
+        deltaAlpha_[cellI] = f * (alphaMax_ - alphaP_[cellI]) - alphaSt_[cellI];
+        // too much volume occupied: release it (50% per time step)
+        if (deltaAlpha_[cellI] < 0.0)
+        {
+            Sds_[cellI] = 0.5*deltaAlpha_[cellI];
+        }
+        // volume too occupy available: deposit at most 80% of dyn hold up
+        else if (depRate_ * deltaAlpha_[cellI] > 0.8 * alphaDyn_[cellI])
+        {
+            Sds_[cellI] = 0.8 * alphaDyn_[cellI];
+        }
+        else
+        {
+            Sds_[cellI] = depRate_ * deltaAlpha_[cellI];
+        }
+    }
 }
 
 
-void FinesFields::integrateFields() 
+void FinesFields::integrateFields()
 {
-  
+
     surfaceScalarField phiSt(linearInterpolate(UsField_) & particleCloud_.mesh().Sf());
     surfaceScalarField phiDyn(linearInterpolate(uDyn_) & particleCloud_.mesh().Sf());
-  
+
     fvScalarMatrix alphaStEqn
     (
           fvm::ddt(alphaSt_)
-	+ fvm::div(phiSt,alphaSt_)
-	==
-	Sds_
+        + fvm::div(phiSt,alphaSt_)
+        ==
+        Sds_
     );
     fvScalarMatrix alphaDynEqn
     (
         fvm::ddt(alphaDyn_)
-	+ fvm::div(phiDyn,alphaDyn_)
-	- fvm::laplacian(diffCoeff_,alphaDyn_)
-	==
-	-Sds_
+        + fvm::div(phiDyn,alphaDyn_)
+        - fvm::laplacian(diffCoeff_,alphaDyn_)
+        ==
+        -Sds_
     );
     alphaStEqn.solve();
     alphaDynEqn.solve();
-    
+
     if(smoothing_)
         particleCloud_.smoothingM().smoothen(alphaDyn_);
-    
+
     // limit hold-ups, should be done more elegantly
-    
+
     scalar alphaStErr(0.0);
     scalar alphaDynErr1(0.0);
     scalar alphaDynErr2(0.0);
     forAll(alphaSt_, cellI)
     {
         if (alphaSt_[cellI] < 0.0)
-	{
-	    alphaStErr += alphaSt_[cellI] * particleCloud_.mesh().V()[cellI];
-	    alphaSt_[cellI] = 0.0;
-	}
-	
-	if (alphaDyn_[cellI] < 0.0)
-	{
-	    alphaDynErr1 += alphaDyn_[cellI] * particleCloud_.mesh().V()[cellI];
-	    alphaDyn_[cellI] = 0.0;
-	}
-	else if (alphaDyn_[cellI] > alphaDynMax_)
-	{
-	    alphaDynErr2 +=  (alphaDyn_[cellI] - alphaDynMax_) * particleCloud_.mesh().V()[cellI];
-	    alphaDyn_[cellI] = alphaDynMax_;
-	}
+        {
+            alphaStErr += alphaSt_[cellI] * particleCloud_.mesh().V()[cellI];
+            alphaSt_[cellI] = 0.0;
+        }
+
+        if (alphaDyn_[cellI] < 0.0)
+        {
+            alphaDynErr1 += alphaDyn_[cellI] * particleCloud_.mesh().V()[cellI];
+            alphaDyn_[cellI] = 0.0;
+        }
+        else if (alphaDyn_[cellI] > alphaDynMax_)
+        {
+            alphaDynErr2 +=  (alphaDyn_[cellI] - alphaDynMax_) * particleCloud_.mesh().V()[cellI];
+            alphaDyn_[cellI] = alphaDynMax_;
+        }
     }
-    
+
     if (verbose_)
     {
         Sout << "[" << Pstream::myProcNo() << "] " << "amount of alphaSt added because of positivity requirement: " << -alphaStErr << endl;
         Sout << "[" << Pstream::myProcNo() << "] " << "amount of alphaDyn added because of positivity requirement: " << -alphaDynErr1 << endl;
         Sout << "[" << Pstream::myProcNo() << "] " << "amount of alphaDyn removed because of max. value: " << -alphaDynErr2 << endl;
     }
-    
+
     alphaSt_.correctBoundaryConditions();
     alphaDyn_.correctBoundaryConditions();
-    
+
     massFluxDyn_ = rhoFine_ * fvc::interpolate(alphaDyn_) * phiDyn;
 }
 
 
-void FinesFields::updateAlphaG() 
+void FinesFields::updateAlphaG()
 {
-  alphaG_ = max(voidfraction_ - alphaSt_ - alphaDyn_, critVoidfraction_);
+    alphaG_ = max(voidfraction_ - alphaSt_ - alphaDyn_, critVoidfraction_);
 }
 
 
-void FinesFields::updateAlphaP() 
+void FinesFields::updateAlphaP()
 {
     alphaP_ = 1.0 - voidfraction_ + SMALL;
 }
 
 
-void FinesFields::updateDHydMix() 
+void FinesFields::updateDHydMix()
 {
     forAll(dHydMix_,cellI)
     {
         scalar aPSt =  alphaP_[cellI] + alphaSt_[cellI];
-	if(aPSt < SMALL || aPSt > 1 - SMALL)
-	  dHydMix_[cellI] = SMALL;
-	else 
-	  dHydMix_[cellI] = 2*(1 - aPSt) / (3*aPSt ) * dSauterMix_[cellI];
+        if(aPSt < SMALL || aPSt > 1 - SMALL)
+            dHydMix_[cellI] = SMALL;
+        else
+            dHydMix_[cellI] = 2*(1 - aPSt) / (3*aPSt ) * dSauterMix_[cellI];
     }
     dHydMix_.correctBoundaryConditions();
 }
@@ -453,54 +454,54 @@ void FinesFields::updateDragCoeff()
     volScalarField Ref = dFine_ * alphaG_ / nuAve_ * mag(U_ - uDyn_);
     scalar Cd(0.0);
     scalar Ref1(0.0);
-    
+
     // calculate drag coefficient for cells
     forAll(DragCoeff_,cellI)
     {
         Ref1 = Ref[cellI];
         if(Ref1 <= SMALL)
-	    Cd = 24.0 / SMALL;
+            Cd = 24.0 / SMALL;
         else if(Ref1 <= 1.0)
-	    Cd = 24.0 / Ref1;
-	else if(Ref1 <= 1000)
-	    Cd = 24 * (1.0 + 0.15 * Foam::pow(Ref1,0.687) ) / Ref1;
-	else
-	    Cd = 0.44;
-	DragCoeff_[cellI] = Cd * beta[cellI];
+            Cd = 24.0 / Ref1;
+        else if(Ref1 <= 1000)
+            Cd = 24 * (1.0 + 0.15 * Foam::pow(Ref1,0.687) ) / Ref1;
+        else
+            Cd = 0.44;
+        DragCoeff_[cellI] = Cd * beta[cellI];
     }
-    
+
     // calculate drag coefficient for faces
     forAll(DragCoeff_.boundaryField(), patchI)
         forAll(DragCoeff_.boundaryField()[patchI], faceI)
         {
             Ref1 = Ref.boundaryField()[patchI][faceI];
             if(Ref1 <= SMALL)
-	        Cd = 24.0 / SMALL;
+                Cd = 24.0 / SMALL;
             else if(Ref1 <= 1.0)
-	        Cd = 24.0 / Ref1;
-	    else if(Ref1 <= 1000)
-	        Cd = 24 * (1.0 + 0.15 * Foam::pow(Ref1,0.687) ) / Ref1;
-	    else
-	        Cd = 0.44;
-	    DragCoeff_.boundaryFieldRef()[patchI][faceI] = Cd * beta.boundaryFieldRef()[patchI][faceI];
+                Cd = 24.0 / Ref1;
+            else if(Ref1 <= 1000)
+                Cd = 24 * (1.0 + 0.15 * Foam::pow(Ref1,0.687) ) / Ref1;
+            else
+                Cd = 0.44;
+            DragCoeff_.boundaryFieldRef()[patchI][faceI] = Cd * beta.boundaryFieldRef()[patchI][faceI];
         }
-        
+
     DragCoeff_ = max( DragCoeff_, dimensionedScalar("SMALL", dimensionSet(1,-3,-1,0,0), SMALL) );
 }
 
 
-void FinesFields::updateDSauter() 
+void FinesFields::updateDSauter()
 {
     forAll(dSauterMix_,cellI)
     {
         scalar aP = alphaP_[cellI];
-	scalar aSt = alphaSt_[cellI];
-	if(aSt < SMALL)
-	    dSauterMix_[cellI] = dSauter_[cellI];
-	else if(aP < SMALL)
-	    dSauterMix_[cellI] = dFine_.value();
-	else
-	    dSauterMix_[cellI] = (aP + aSt) / (aP / dSauter_[cellI] + aSt / dFine_.value() );
+        scalar aSt = alphaSt_[cellI];
+        if(aSt < SMALL)
+            dSauterMix_[cellI] = dSauter_[cellI];
+        else if(aP < SMALL)
+            dSauterMix_[cellI] = dFine_.value();
+        else
+            dSauterMix_[cellI] = (aP + aSt) / (aP / dSauter_[cellI] + aSt / dFine_.value() );
     }
     dSauterMix_.correctBoundaryConditions();
 }
@@ -514,7 +515,7 @@ void FinesFields::updateFanningCoeff()
 }
 
 
-void FinesFields::updateFroude() 
+void FinesFields::updateFroude()
 {
     // seems like different authors use different conventions for the Froude number
     // Chen et al. (1994) define it in terms of a superficial velocity,
@@ -534,27 +535,27 @@ void FinesFields::updateUDyn()
     volScalarField denom = FanningCoeff_ + DragCoeff_;
 
     uDyn_ = num / denom;
-    
+
     // limit uDyn for stability reasons
     forAll(uDyn_,cellI)
     {
         scalar mU(mag(U_[cellI]));
-	scalar muDyn(mag(uDyn_[cellI]));
+        scalar muDyn(mag(uDyn_[cellI]));
         if(muDyn > mU && muDyn > SMALL)
-	{
-	    uDyn_[cellI] *= mU / muDyn;
-	}
+        {
+            uDyn_[cellI] *= mU / muDyn;
+        }
     }
-    
+
     forAll(uDyn_.boundaryField(), patchI)
         forAll(uDyn_.boundaryField()[patchI], faceI)
         {
-	    scalar mU(mag(U_.boundaryField()[patchI][faceI]));
-	    scalar muDyn(mag(uDyn_.boundaryField()[patchI][faceI]));
+            scalar mU(mag(U_.boundaryField()[patchI][faceI]));
+            scalar muDyn(mag(uDyn_.boundaryField()[patchI][faceI]));
             if(muDyn > mU && muDyn > SMALL)
-	    {
-	        uDyn_.boundaryFieldRef()[patchI][faceI] *= mU / muDyn;
-	    }
+            {
+                uDyn_.boundaryFieldRef()[patchI][faceI] *= mU / muDyn;
+            }
         }
 }
 

src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.H

@@ -39,11 +39,11 @@ class FinesFields
 {
 private:
     cfdemCloud& particleCloud_;
-  
+
     dictionary propsDict_;
-    
+
     bool smoothing_;
-    
+
     bool verbose_;
 
     word velFieldName_;
@@ -54,95 +54,95 @@ private:
 
     volScalarField& voidfraction_;
 
-    word UsFieldName_;   
+    word UsFieldName_;
 
     const volVectorField& UsField_;
-    
+
     word pFieldName_;
 
     const volScalarField& p_;
-    
+
     word rhoGFieldName_;
 
     const volScalarField& rhoG_;
-    
+
     const volScalarField& dSauter_;
-    
+
     volScalarField alphaG_;
-    
+
     volScalarField alphaDyn_;
-    
+
     volScalarField alphaP_;
 
     volScalarField alphaSt_;
 
     volScalarField deltaAlpha_;
-    
+
     volScalarField dHydMix_;
-    
+
     volScalarField DragCoeff_;
-    
+
     volScalarField dSauterMix_;
-    
+
     volScalarField FanningCoeff_;
-    
+
     volScalarField Froude_;
-    
+
     volScalarField Sds_;
-    
+
     //volVectorField massFluxDyn_;
     surfaceScalarField massFluxDyn_;
-    
+
     volVectorField uDyn_;
 
     dimensionedScalar dFine_;
-    
+
     dimensionedScalar diffCoeff_;
-    
+
     dimensionedScalar nuAve_;
-    
+
     dimensionedScalar rhoFine_;
-    
+
     const dimensionedVector g_;
-    
+
     scalar alphaDynMax_;
-    
+
     scalar alphaMax_;
-    
+
     scalar critVoidfraction_;
-    
+
     scalar depRate_;
-    
+
     scalar exponent_;
-    
+
     scalar nCrit_;
-    
+
     scalar poresizeWidth_;
-    
+
     scalar prefactor_;
-    
+
     scalar ratioHydraulicPore_;
-       
+
     void calcSource();
-    
+
     void integrateFields();
-    
+
     void updateAlphaG();
-       
+
     void updateAlphaP();
-    
+
     void updateDHydMix();
-    
+
     void updateDragCoeff();
-    
+
     void updateDSauter();
-    
+
     void updateFanningCoeff();
-    
+
     void updateFroude();
-      
+
     void updateUDyn();
-    
+
 public:
 
     //- Runtime type information
@@ -165,7 +165,7 @@ public:
 
     // Member Functions
 
-    
+
     void update();
 };
 

src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.C

@@ -75,8 +75,8 @@ GidaspowDrag::GidaspowDrag
     switchingVoidfraction_(0.8)
 {
     //Append the field names to be probed
-    particleCloud_.probeM().initialize(typeName, "gidaspowDrag.logDat");
-    particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must  be the force
+    particleCloud_.probeM().initialize(typeName, typeName+".logDat");
+    particleCloud_.probeM().vectorFields_.append("dragForce"); // first entry must  be the force
     particleCloud_.probeM().vectorFields_.append("Urel");
     particleCloud_.probeM().scalarFields_.append("Rep");
     particleCloud_.probeM().scalarFields_.append("betaP");
@@ -86,18 +86,18 @@ GidaspowDrag::GidaspowDrag
     // init force sub model
     setForceSubModels(propsDict_);
     // define switches which can be read from dict
-    forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(2,true); // activate implDEM switch
-    forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch
-    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
-    forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity switch
+    forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
+    forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM,true); // activate implDEM switch
+    forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(SW_SCALAR_VISCOSITY,true); // activate scalarViscosity switch
     forceSubM(0).readSwitches();
 
     particleCloud_.checkCG(true);
     if (propsDict_.found("scale"))
-        scaleDia_=scalar(readScalar(propsDict_.lookup("scale")));
+        scaleDia_ = scalar(readScalar(propsDict_.lookup("scale")));
     if (propsDict_.found("scaleDrag"))
-        scaleDrag_=scalar(readScalar(propsDict_.lookup("scaleDrag")));
+        scaleDrag_ = scalar(readScalar(propsDict_.lookup("scaleDrag")));
 
     if (propsDict_.found("switchingVoidfraction"))
         switchingVoidfraction_ = readScalar(propsDict_.lookup("switchingVoidfraction"));
@@ -114,10 +114,13 @@ GidaspowDrag::~GidaspowDrag()
 
 void GidaspowDrag::setForce() const
 {
-    if (scaleDia_ > 1)
+    if (scaleDia_ > 1.)
+    {
         Info << "Gidaspow using scale = " << scaleDia_ << endl;
-    else if (particleCloud_.cg() > 1){
-        scaleDia_=particleCloud_.cg();
+    }
+    else if (particleCloud_.cg() > 1.)
+    {
+        scaleDia_ = particleCloud_.cg();
         Info << "Gidaspow using scale from liggghts cg = " << scaleDia_ << endl;
     }
 
@@ -128,7 +131,7 @@ void GidaspowDrag::setForce() const
     scalar voidfraction(1);
     vector Ufluid(0,0,0);
     vector drag(0,0,0);
-    label cellI=0;
+    label cellI = 0;
 
     vector Us(0,0,0);
     vector Ur(0,0,0);
@@ -145,13 +148,13 @@ void GidaspowDrag::setForce() const
 
     vector dragExplicit(0,0,0);
     scalar dragCoefficient(0);
-    
+
     interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
     interpolationCellPoint<vector> UInterpolator_(U_);
 
     #include "setupProbeModel.H"
 
-    for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
+    for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
         //if(mask[index][0])
         //{
@@ -160,26 +163,26 @@ void GidaspowDrag::setForce() const
             dragExplicit = vector(0,0,0);
             betaP = 0;
             Vs = 0;
-            Ufluid =vector(0,0,0);
-            voidfraction=0;
+            Ufluid = vector(0,0,0);
+            voidfraction = 0;
             dragCoefficient = 0;
 
             if (cellI > -1) // particle Found
             {
 
-                if( forceSubM(0).interpolation() )
+                if ( forceSubM(0).interpolation() )
                 {
-	                position     = particleCloud_.position(index);
+                    position     = particleCloud_.position(index);
                     voidfraction = voidfractionInterpolator_.interpolate(position,cellI);
                     Ufluid       = UInterpolator_.interpolate(position,cellI);
                     //Ensure interpolated void fraction to be meaningful
                     // Info << " --> voidfraction: " << voidfraction << endl;
-                    if(voidfraction>1.00) voidfraction = 1.0;
-                    if(voidfraction<0.10) voidfraction = 0.10;
+                    if (voidfraction > 1.0) voidfraction = 1.0;
+                    else if (voidfraction < 0.1) voidfraction = 0.10;
                 }
                 else
                 {
-					voidfraction = voidfraction_[cellI];
+                    voidfraction = voidfraction_[cellI];
                     Ufluid = U_[cellI];
                 }
 
@@ -190,14 +193,14 @@ void GidaspowDrag::setForce() const
                 rho = rhoField[cellI];
                 nuf = nufField[cellI];
 
-                Rep=0.0;
+                Rep = 0.0;
                 localPhiP = 1.0f-voidfraction+SMALL;
                 Vs = ds*ds*ds*M_PI/6;
 
                 // calc particle's drag coefficient (i.e., Force per unit slip velocity and per m³ PARTICLE)
-                if(voidfraction > switchingVoidfraction_) //dilute
+                if (voidfraction > switchingVoidfraction_) //dilute
                 {
-                    Rep=ds/scaleDia_*voidfraction*magUr/nuf;
+                    Rep = ds/scaleDia_*voidfraction*magUr/nuf;
                     CdMagUrLag = (24.0*nuf/(ds/scaleDia_*voidfraction)) //1/magUr missing here, but compensated in expression for betaP!
                                  *(scalar(1.0)+0.15*Foam::pow(Rep, 0.687));
 
@@ -218,7 +221,7 @@ void GidaspowDrag::setForce() const
 
                 // calc particle's drag
                 dragCoefficient = Vs*betaP*scaleDrag_;
-                if (modelType_=="B")
+                if (modelType_ == "B")
                     dragCoefficient /= voidfraction;
 
                 drag = dragCoefficient * Ur;
@@ -226,7 +229,7 @@ void GidaspowDrag::setForce() const
                 // explicitCorr
                 forceSubM(0).explicitCorr(drag,dragExplicit,dragCoefficient,Ufluid,U_[cellI],Us,UsField_[cellI],forceSubM(0).verbose());
 
-                if(forceSubM(0).verbose() && index >=0 && index <2)
+                if (forceSubM(0).verbose() && index >= 0 && index < 2)
                 {
                     Pout << "cellI = " << cellI << endl;
                     Pout << "index = " << index << endl;
@@ -244,7 +247,7 @@ void GidaspowDrag::setForce() const
                 }
 
                 //Set value fields and write the probe
-                if(probeIt_)
+                if (probeIt_)
                 {
                     #include "setupProbeModelfields.H"
                     vValues.append(drag);   //first entry must the be the force

src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.H

@@ -74,7 +74,7 @@ private:
 
     const scalar phi_;
 
-    word UsFieldName_;   
+    word UsFieldName_;
 
     const volVectorField& UsField_;  // the average particle velocity field
 

src/lagrangian/cfdemParticle/subModels/forceModel/KochHillDrag/KochHillDrag.C

@@ -74,7 +74,7 @@ KochHillDrag::KochHillDrag
     scaleDrag_(1.)
 {
     //Append the field names to be probed
-    particleCloud_.probeM().initialize(typeName, "kochHillDrag.logDat");
+    particleCloud_.probeM().initialize(typeName, typeName+".logDat");
     particleCloud_.probeM().vectorFields_.append("dragForce");      //first entry must the be the force
     particleCloud_.probeM().vectorFields_.append("Urel");           //other are debug
     particleCloud_.probeM().scalarFields_.append("Rep");            //other are debug
@@ -86,12 +86,12 @@ KochHillDrag::KochHillDrag
     setForceSubModels(propsDict_);
 
     // define switches which can be read from dict
-    forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(2,true); // activate implDEM switch
-    forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch
-    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
-    forceSubM(0).setSwitchesList(7,true); // activate implForceDEMacc switch
-    forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity switch
+    forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
+    forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM,true); // activate implDEM switch
+    forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM_ACCUMULATED,true); // activate implForceDEMacc switch
+    forceSubM(0).setSwitchesList(SW_SCALAR_VISCOSITY,true); // activate scalarViscosity switch
 
     // read those switches defined above, if provided in dict
     forceSubM(0).readSwitches();
@@ -99,9 +99,9 @@ KochHillDrag::KochHillDrag
     particleCloud_.checkCG(true);
 
     if (propsDict_.found("scale"))
-        scaleDia_=scalar(readScalar(propsDict_.lookup("scale")));
+        scaleDia_ = scalar(readScalar(propsDict_.lookup("scale")));
     if (propsDict_.found("scaleDrag"))
-        scaleDrag_=scalar(readScalar(propsDict_.lookup("scaleDrag")));
+        scaleDrag_ = scalar(readScalar(propsDict_.lookup("scaleDrag")));
 }
 
 
@@ -134,7 +134,7 @@ void KochHillDrag::setForce() const
     vector drag(0,0,0);
     vector dragExplicit(0,0,0);
     scalar dragCoefficient(0);
-    label cellI=0;
+    label cellI = 0;
 
     vector Us(0,0,0);
     vector Ur(0,0,0);
@@ -154,16 +154,16 @@ void KochHillDrag::setForce() const
 
     #include "setupProbeModel.H"
 
-    for (int index=0; index<particleCloud_.numberOfParticles(); index++)
+    for (int index = 0; index<particleCloud_.numberOfParticles(); ++index)
     {
             cellI = particleCloud_.cellIDs()[index][0];
             drag = vector(0,0,0);
             dragExplicit = vector(0,0,0);
-            dragCoefficient=0;
+            dragCoefficient = 0;
             betaP = 0;
             Vs = 0;
-            Ufluid =vector(0,0,0);
-            voidfraction=0;
+            Ufluid = vector(0,0,0);
+            voidfraction = 0;
 
             if (cellI > -1) // particle Found
             {
@@ -174,8 +174,8 @@ void KochHillDrag::setForce() const
                     Ufluid = UInterpolator_.interpolate(position,cellI);
                     //Ensure interpolated void fraction to be meaningful
                     // Info << " --> voidfraction: " << voidfraction << endl;
-                    if (voidfraction > 1.00) voidfraction = 1.00;
-                    if (voidfraction < 0.40) voidfraction = 0.40;
+                    if (voidfraction > 1.0) voidfraction = 1.0;
+                    else if (voidfraction < 0.4) voidfraction = 0.4;
                 }
                 else
                 {
@@ -227,10 +227,10 @@ void KochHillDrag::setForce() const
 
                     // calc particle's drag
                     dragCoefficient = Vs*betaP*scaleDrag_;
-                    if (modelType_=="B")
+                    if (modelType_ == "B")
                         dragCoefficient /= voidfraction;
 
-                    if (forceSubM(0).switches()[7]) // implForceDEMaccumulated=true
+                    if (forceSubM(0).switches()[SW_IMPL_FORCE_DEM_ACCUMULATED]) // implForceDEMaccumulated=true
                     {
                         //get drag from the particle itself
                         for (int j=0; j<3; j++) drag[j] = particleCloud_.fAccs()[index][j]/couplingInterval;
@@ -270,7 +270,7 @@ void KochHillDrag::setForce() const
                     sValues.append(betaP);
                     sValues.append(voidfraction);
                     particleCloud_.probeM().writeProbe(index, sValues, vValues);
-                }    
+                }
             }
 
             // write particle based data to global array

src/lagrangian/cfdemParticle/subModels/forceModel/KochHillDrag/KochHillDrag.H

@@ -29,7 +29,7 @@ Description
     and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
 
     Koch, Hill drag law
-    based on Koch Hill 2001,"Inertial effects in suspensions and porous-media 
+    based on Koch Hill 2001,"Inertial effects in suspensions and porous-media
     flows", Annual Review of fluid mechanics.
     including interpolation of the velocity to the exact position
     including drag coefficient for implicit drag for DEM
@@ -72,9 +72,9 @@ private:
 
     const volScalarField& voidfraction_;
 
-    word UsFieldName_;   
+    word UsFieldName_;
 
-    const volVectorField& UsField_; 
+    const volVectorField& UsField_;
 
     mutable scalar scaleDia_;
 

src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.C

@@ -79,16 +79,16 @@ KochHillRWDrag::KochHillRWDrag
     RanGen_(label(0))
 {
 
-    if (propsDict_.found("verbose")) verbose_=true;
-    if (propsDict_.found("interpolation")) interpolation_=true;
-    if (propsDict_.found("randomTauE")) randomTauE_=true;
+    if (propsDict_.found("verbose")) verbose_ = true;
+    if (propsDict_.found("interpolation")) interpolation_ = true;
+    if (propsDict_.found("randomTauE")) randomTauE_ = true;
 
     // init force sub model
     setForceSubModels(propsDict_);
 
     // define switches which can be read from dict
-    forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(2,true); // activate implDEM switch
+    forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
+    forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM,true); // activate implDEM switch
 
     // read those switches defined above, if provided in dict
     forceSubM(0).readSwitches();
@@ -120,8 +120,8 @@ KochHillRWDrag::KochHillRWDrag
 
 KochHillRWDrag::~KochHillRWDrag()
 {
-    delete partTime_;
-    delete partUfluct_;
+    particleCloud_.dataExchangeM().destroy(partTime_, 1);
+    particleCloud_.dataExchangeM().destroy(partUfluct_, 3);
 }
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -152,7 +152,7 @@ void KochHillRWDrag::setForce() const
     vector drag(0,0,0);
     vector dragExplicit(0,0,0);
     scalar dragCoefficient(0);
-    label cellI=0;
+    label cellI = 0;
 
     vector Us(0,0,0);
     vector Ur(0,0,0);
@@ -195,17 +195,17 @@ void KochHillRWDrag::setForce() const
 
     //Info << "RW-TEST: We are in setForce() at t = " << t << endl; // TEST-Output
 
-    for (int index=0; index<particleCloud_.numberOfParticles(); index++)
+    for (int index = 0; index<particleCloud_.numberOfParticles(); ++index)
     {
         //if (mask[index][0])
         //{
             cellI = particleCloud_.cellIDs()[index][0];
             drag = vector(0,0,0);
             dragExplicit = vector(0,0,0);
-            dragCoefficient=0;
+            dragCoefficient = 0;
             betaP = 0;
             Vs = 0;
-            Ufluid =vector(0,0,0);
+            Ufluid = vector(0,0,0);
 
             // Pout << "RW-TEST: cellI = " << cellI << endl; // TEST-Output
             if (cellI > -1) // particle Found
@@ -217,8 +217,8 @@ void KochHillRWDrag::setForce() const
                     Ufluid = UInterpolator_.interpolate(position,cellI);
                     //Ensure interpolated void fraction to be meaningful
                     // Info << " --> voidfraction: " << voidfraction << endl;
-                    if (voidfraction > 1.00) voidfraction = 1.00;
-                    if (voidfraction < 0.40) voidfraction = 0.40;
+                    if (voidfraction > 1.0) voidfraction = 1.0;
+                    else if (voidfraction < 0.4) voidfraction = 0.4;
                 }
                 else
                 {

src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.H

@@ -29,7 +29,7 @@ Description
     and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
 
     Koch, Hill drag law
-    based on Koch Hill 2001,"Inertial effects in suspensions and porous-media 
+    based on Koch Hill 2001,"Inertial effects in suspensions and porous-media
     flows", Annual Review of fluid mechanics.
     including interpolation of the velocity to the exact position
     including drag coefficient for implicit drag for DEM
@@ -75,7 +75,7 @@ private:
 
     const volScalarField& voidfraction_;
 
-    word UsFieldName_;   
+    word UsFieldName_;
 
     const volVectorField& UsField_;  // the average particle velocity field
 

src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.C

@@ -89,9 +89,9 @@ LaEuScalarTemp::LaEuScalarTemp
     setForceSubModels(propsDict_);
 
     // define switches which can be read from dict
-    forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch
-    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
-    forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity switch
+    forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(SW_SCALAR_VISCOSITY,true); // activate scalarViscosity switch
 
     // read those switches defined above, if provided in dict
     forceSubM(0).readSwitches();
@@ -105,15 +105,15 @@ LaEuScalarTemp::LaEuScalarTemp
 
 LaEuScalarTemp::~LaEuScalarTemp()
 {
-    delete partTemp_;
-    delete partHeatFlux_;
+    particleCloud_.dataExchangeM().destroy(partTemp_,1);
+    particleCloud_.dataExchangeM().destroy(partHeatFlux_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
 void LaEuScalarTemp::allocateMyArrays() const
 {
     // get memory for 2d arrays
-    double initVal=0.0;
+    double initVal = 0.0;
     particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1);  // field/initVal/with/lenghtFromLigghts
     particleCloud_.dataExchangeM().allocateArray(partHeatFlux_,initVal,1);
 }
@@ -156,14 +156,14 @@ void LaEuScalarTemp::manipulateScalarField(volScalarField& EuField) const
     interpolationCellPoint<vector> UInterpolator_(U_);
     interpolationCellPoint<scalar> TInterpolator_(tempField_);
 
-    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
         //if(particleCloud_.regionM().inRegion()[index][0])
         //{
             cellI = particleCloud_.cellIDs()[index][0];
-            if(cellI >= 0)
+            if (cellI >= 0)
             {
-                if(forceSubM(0).interpolation())
+                if (forceSubM(0).interpolation())
                 {
                     vector position = particleCloud_.position(index);
                     voidfraction = voidfractionInterpolator_.interpolate(position,cellI);
@@ -206,7 +206,7 @@ void LaEuScalarTemp::manipulateScalarField(volScalarField& EuField) const
                 partHeatFlux_[index][0] = partHeatFlux;
 
 
-                if(forceSubM(0).verbose() && index >=0 && index <2)
+                if(forceSubM(0).verbose() && index >= 0 && index < 2)
                 {
                     Info << "partHeatFlux = " << partHeatFlux << endl;
                     Info << "magUr = " << magUr << endl;

src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.H

@@ -71,7 +71,7 @@ private:
 
     word voidfractionFieldName_;
 
-    const volScalarField& voidfraction_;    // ref to voidfraction field
+    const volScalarField& voidfraction_;         // ref to voidfraction field
 
     scalar maxSource_;                           // max (limited) value of src field
 
@@ -81,7 +81,7 @@ private:
 
     word partTempName_;
 
-    mutable double **partTemp_;                          // Lagrangian array
+    mutable double **partTemp_;                  // Lagrangian array
 
     word partHeatFluxName_;
 

src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.C

@@ -5,7 +5,8 @@
     www.cfdem.com
                                 Christoph Goniva, christoph.goniva@cfdem.com
                                 Copyright 2009-2012 JKU Linz
-                                Copyright 2012-     DCS Computing GmbH, Linz
+                                Copyright 2012-2015 DCS Computing GmbH, Linz
+                                Copyright 2015-     JKU Linz
 -------------------------------------------------------------------------------
 License
     This file is part of CFDEMcoupling.
@@ -27,6 +28,20 @@ License
 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).
+
+    This function is based on the derivation in R. Mei,
+    An approximate expression for shear lift force on a spherical  particle at a
+    finite Reynolds number,
+    Int. J. Multiph. Flow 18 (1992) 145–147
+
+    The data for this functions is based on J.B. Mclaughlin,
+    Inertial migration of a small sphere in linear shear flows,
+    Journal of Fluid Mechanics. 224 (1991) 261-274.
+
+    The second order terms are based on E. Loth and A. J. Dorgan,
+    An equation of motion for particles of finite Reynolds number and size,
+    Environ. Fluid Mech. 9 (2009) 187–206
+    and can be added to the lift coefficient if desired
 \*---------------------------------------------------------------------------*/
 
 #include "error.H"
@@ -73,22 +88,22 @@ MeiLift::MeiLift
     // init force sub model
     setForceSubModels(propsDict_);
     // define switches which can be read from dict
-    forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch
-    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
-    forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity switch
+    forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
+    forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(SW_SCALAR_VISCOSITY,true); // activate scalarViscosity switch
     forceSubM(0).readSwitches();
 
     particleCloud_.checkCG(false);
 
     //Append the field names to be probed
-    particleCloud_.probeM().initialize(typeName, "meiLift.logDat");
+    particleCloud_.probeM().initialize(typeName, typeName+".logDat");
     particleCloud_.probeM().vectorFields_.append("liftForce"); //first entry must the be the force
-    particleCloud_.probeM().vectorFields_.append("Urel");        //other are debug
-    particleCloud_.probeM().vectorFields_.append("vorticity");  //other are debug
-    particleCloud_.probeM().scalarFields_.append("Rep");          //other are debug
-    particleCloud_.probeM().scalarFields_.append("Rew");          //other are debug
-    particleCloud_.probeM().scalarFields_.append("J_star");       //other are debug
+    particleCloud_.probeM().vectorFields_.append("Urel");      //other are debug
+    particleCloud_.probeM().vectorFields_.append("vorticity"); //other are debug
+    particleCloud_.probeM().scalarFields_.append("Rep");       //other are debug
+    particleCloud_.probeM().scalarFields_.append("Rew");       //other are debug
+    particleCloud_.probeM().scalarFields_.append("J_star");    //other are debug
     particleCloud_.probeM().writeHeader();
 }
 
@@ -124,6 +139,7 @@ void MeiLift::setForce() const
     scalar J_star(0);
     scalar Omega_eq(0);
     scalar alphaStar(0);
+    scalar epsilonSqr(0.0);
     scalar epsilon(0);
     scalar omega_star(0);
     vector vorticity(0,0,0);
@@ -134,96 +150,102 @@ void MeiLift::setForce() const
 
     #include "setupProbeModel.H"
 
-    for(int index = 0;index <  particleCloud_.numberOfParticles(); index++)
+    for (int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
         //if(mask[index][0])
         //{
-            lift           = vector::zero;
+            lift = vector::zero;
             label cellI = particleCloud_.cellIDs()[index][0];
 
             if (cellI > -1) // particle Found
             {
                 Us = particleCloud_.velocity(index);
 
-                if( forceSubM(0).interpolation() )
+                if (forceSubM(0).interpolation())
                 {
-	                position       = particleCloud_.position(index);
-                    Ur               = UInterpolator_.interpolate(position,cellI) 
-                                        - Us;
-                    vorticity       = VorticityInterpolator_.interpolate(position,cellI);
+                    position  = particleCloud_.position(index);
+                    Ur        = UInterpolator_.interpolate(position,cellI) - Us;
+                    vorticity = VorticityInterpolator_.interpolate(position,cellI);
                 }
                 else
                 {
-                    Ur =  U_[cellI]
-                          - Us;
-                    vorticity=vorticityField[cellI];
+                    Ur        = U_[cellI] - Us;
+                    vorticity = vorticityField[cellI];
                 }
 
-                magUr           = mag(Ur);
+                magUr        = mag(Ur);
                 magVorticity = mag(vorticity);
 
                 if (magUr > 0 && magVorticity > 0)
                 {
-                    ds  = 2*particleCloud_.radius(index);
+                    ds  = 2. * particleCloud_.radius(index);
                     nuf = nufField[cellI];
                     rho = rhoField[cellI];
 
                     // calc dimensionless properties
                     Rep = ds*magUr/nuf;
-		            Rew = magVorticity*ds*ds/nuf;
+                    Rew = magVorticity*ds*ds/nuf;
 
-                    alphaStar   = magVorticity*ds/magUr/2.0;
-                    epsilon       =  sqrt(2.0*alphaStar /Rep );
-                    omega_star=2.0*alphaStar;
+                    omega_star = magVorticity * ds / magUr;
+                    alphaStar = 0.5 * omega_star;
+                    epsilonSqr = omega_star / Rep;
+                    epsilon = sqrt(epsilonSqr);
 
                     //Basic model for the correction to the Saffman lift
-                    //Based on McLaughlin (1991)
-                    if(epsilon < 0.1)
+                    //McLaughlin (1991), Mei (1992), Loth and Dorgan (2009)
+                    //J_star = 0.443 * J
+                    if (epsilon < 0.1) //epsilon << 1
                     {
-                        J_star = -140 *epsilon*epsilon*epsilon*epsilon*epsilon 
-                                             *log( 1./(epsilon*epsilon+SMALL) );
+                        //McLaughlin (1991), Eq (3.27): J = 32 * pi^2 * epsilon^5 * ln(1 / epsilon^2)
+                        J_star = -140.0 * epsilonSqr * epsilonSqr * epsilon * log(1. / (epsilonSqr+SMALL));
                     }
-                    else if(epsilon > 20)
+                    else if (epsilon > 20.0) //epsilon >> 1
                     {
-                      J_star = 1.0-0.287/(epsilon*epsilon+SMALL);
+                        //McLaughlin (1991), Eq (3.26): J = 2.255 - 0.6463 / epsilon^2
+                        J_star = 1.0 - 0.287 / epsilonSqr;
                     }
                     else
                     {
-                     J_star = 0.3
-                                *(     1.0
-                                      +tanh(  2.5 * log10(epsilon+0.191)  )
-                                 )
-                                *(    0.667
-                                     +tanh(  6.0 * (epsilon-0.32)  )
-                                  );
+                        //Mei (1992), Eq (10)
+                        //Loth and Dorgan (2009), Eq (32)
+                        J_star = 0.3
+                                 * (1.0   + tanh(2.5 * (log10(epsilon) + 0.191)))
+                                 * (0.667 + tanh(6.0 * (      epsilon  - 0.32 )));
                     }
-                    Cl=J_star*4.11*epsilon; //multiply McLaughlin's correction to the basic Saffman model
+                    //Loth and Dorgan (2009), Eq (31): Saffman lift-coefficient: ClSaff = 12.92 / pi * epsilon ~ 4.11 * epsilon
+                    //Loth and Dorgan (2009), Eq (32)
+                    Cl = J_star * 4.11 * epsilon; //multiply correction to the basic Saffman model
 
-                    //Second order terms given by Loth and Dorgan 2009 
-                    if(useSecondOrderTerms_)
-                    {   
-                        Omega_eq = omega_star/2.0*(1.0-0.0075*Rew)*(1.0-0.062*sqrt(Rep)-0.001*Rep);
-                        Cl_star=1.0-(0.675+0.15*(1.0+tanh(0.28*(omega_star/2.0-2.0))))*tanh(0.18*sqrt(Rep));
-                        Cl += Omega_eq*Cl_star;
+                    //Second order terms given by Loth and Dorgan (2009)
+                    if (useSecondOrderTerms_)
+                    {
+                        scalar sqrtRep = sqrt(Rep);
+                        //Loth and Dorgan (2009), Eq (34)
+                        Cl_star = 1.0 - (0.675 + 0.15 * (1.0 + tanh(0.28 * (alphaStar - 2.0)))) * tanh(0.18 * sqrtRep);
+                        //Loth and Dorgan (2009), Eq (38)
+                        Omega_eq = alphaStar * (1.0 - 0.0075 * Rew) * (1.0 - 0.062 * sqrtRep - 0.001 * Rep);
+                        //Loth and Dorgan (2009), Eq (39)
+                        Cl += Omega_eq * Cl_star;
                     }
 
-                    lift =  0.125*M_PI
-                           *rho
-                           *Cl  
-                           *magUr*Ur^vorticity/magVorticity
-                           *ds*ds;
+                    //Loth and Dorgan (2009), Eq (27)
+                    lift = 0.125 * constant::mathematical::pi
+                           * rho
+                           * Cl
+                           * magUr * Ur ^ vorticity / magVorticity
+                           * ds * ds;
 
-                    if (modelType_=="B")
+                    if (modelType_ == "B")
                     {
                         voidfraction = particleCloud_.voidfraction(index);
                         lift /= voidfraction;
                     }
                 }
 
-                //**********************************        
+                //**********************************
                 //SAMPLING AND VERBOSE OUTOUT
-                if( forceSubM(0).verbose() )
-                {   
+                if ( forceSubM(0).verbose() )
+                {
                     Pout << "index = " << index << endl;
                     Pout << "Us = " << Us << endl;
                     Pout << "Ur = " << Ur << endl;
@@ -240,7 +262,7 @@ void MeiLift::setForce() const
                 }
 
                 //Set value fields and write the probe
-                if(probeIt_)
+                if (probeIt_)
                 {
                     #include "setupProbeModelfields.H"
                     vValues.append(lift);   //first entry must the be the force
@@ -252,7 +274,7 @@ void MeiLift::setForce() const
                     particleCloud_.probeM().writeProbe(index, sValues, vValues);
                 }
                 // END OF SAMPLING AND VERBOSE OUTOUT
-                //**********************************        
+                //**********************************
 
             }
             // write particle based data to global array

src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.H

@@ -5,7 +5,8 @@
     www.cfdem.com
                                 Christoph Goniva, christoph.goniva@cfdem.com
                                 Copyright 2009-2012 JKU Linz
-                                Copyright 2012-     DCS Computing GmbH, Linz
+                                Copyright 2012-2015 DCS Computing GmbH, Linz
+                                Copyright 2015-     JKU Linz
 -------------------------------------------------------------------------------
 License
     This file is part of CFDEMcoupling.
@@ -28,18 +29,20 @@ 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).
 
-     This function is based on the derivation in R. Kurose, S. Komori, 
-     Drag and lift forces on a rotating sphere in a linear shear flow, 
-     Journal of Fluid Mechanics. 384 (1999) 183-206.
+    This function is based on the derivation in R. Mei,
+    An approximate expression for shear lift force on a spherical  particle at a
+    finite Reynolds number,
+    Int. J. Multiph. Flow 18 (1992) 145–147
 
-     The data for this functions is based on J.B. Mclaughlin, 
-     Inertial migration of a small sphere in linear shear flows, 
-     Journal of Fluid Mechanics. 224 (1991) 261-274.
+    The data for this functions is based on J.B. Mclaughlin,
+    Inertial migration of a small sphere in linear shear flows,
+    Journal of Fluid Mechanics. 224 (1991) 261-274.
 
-     The second order terms are based on: 
-     Mei Lift force following Loth and Dorgan 2009,
-     and can be added to the lift coefficient if desired
-     (contribution from RQ)
+    The second order terms are based on E. Loth and A. J. Dorgan,
+    An equation of motion for particles of finite Reynolds number and size,
+    Environ. Fluid Mech. 9 (2009) 187–206
+    and can be added to the lift coefficient if desired
+    (contribution from RQ)
 
     - including interpolation of the velocity to the particle position (optional)
     - including output to file for testing/data analysis (optional)

src/lagrangian/cfdemParticle/subModels/forceModel/SchillerNaumannDrag/SchillerNaumannDrag.C

@@ -69,7 +69,7 @@ SchillerNaumannDrag::SchillerNaumannDrag
     U_(sm.mesh().lookupObject<volVectorField> (velFieldName_))
 {
     //Append the field names to be probed
-    particleCloud_.probeM().initialize(typeName, "schillerNaumannDrag.logDat");
+    particleCloud_.probeM().initialize(typeName, typeName+".logDat");
     particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must the be the force
     particleCloud_.probeM().vectorFields_.append("Urel");        //other are debug
     particleCloud_.probeM().scalarFields_.append("Rep");          //other are debug
@@ -82,7 +82,7 @@ SchillerNaumannDrag::SchillerNaumannDrag
     setForceSubModels(propsDict_);
 
     // define switches which can be read from dict
-    forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch
+    forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
 
     // read those switches defined above, if provided in dict
     forceSubM(0).readSwitches();