Merge pull request #27 from ParticulateFlow/feature/cherry_pick_PUBLIC

Feature/cherry pick public

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.

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/cfdTools/setupProbeModel.H

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

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C

@@ -126,6 +126,7 @@ cfdemCloud::cfdemCloud
         mesh,
         dimensionedScalar("zero", dimensionSet(0,0,-1,0,0), 0)  // 1/s
     ),
+    checkPeriodicCells_(false),
     turbulence_
     (
         mesh.lookupObject<turbulenceModel>
@@ -196,7 +197,7 @@ cfdemCloud::cfdemCloud
         clockModel::New
         (
             couplingProperties_,
-            *this
+            mesh.time()
         )
     ),
     smoothingModel_
@@ -225,6 +226,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;
@@ -289,7 +298,37 @@ 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;
+    }
+    else 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 +451,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 +469,32 @@ void cfdemCloud::setPos(double**& pos)
         }
     }
 }
+
 // * * * * * * * * * * * * * * * ACCESS  * * * * * * * * * * * * * //
 
-label cfdemCloud::particleCell(int index)
+label cfdemCloud::particleCell(int index) const
 {
-    label cellI = cellIDs()[index][0];
+    return cellIDs()[index][0];
-    return cellI;
 }
 
-vector cfdemCloud::position(int index)
+vector cfdemCloud::position(int index) const
 {
-    vector pos;
+    return vector(positions()[index][0],positions()[index][1],positions()[index][2]);
-    for(int i=0;i<3;i++) pos[i] = positions()[index][i];
-    return pos;
 }
 
-vector cfdemCloud::velocity(int index)
+vector cfdemCloud::velocity(int index) const
 {
-    vector vel;
+    return vector(velocities()[index][0],velocities()[index][1],velocities()[index][2]);
-    for(int i=0;i<3;i++) vel[i] = velocities()[index][i];
-    return vel;
 }
 
-vector cfdemCloud::expForce(int index)
+vector cfdemCloud::expForce(int index) const
 {
-    vector force;
+    return vector(DEMForces()[index][0],DEMForces()[index][1],DEMForces()[index][2]);
-    for(int i=0;i<3;i++) force[i] = DEMForces()[index][i];
-    return force;
 }
 
-vector cfdemCloud::fluidVel(int index)
+vector cfdemCloud::fluidVel(int index) const
 {
-    vector vel;
+    return vector(fluidVels()[index][0],fluidVels()[index][1],fluidVels()[index][2]);
-    for(int i=0;i<3;i++) vel[i] = fluidVels()[index][i];
-    return vel;
 }
 
 const forceModel& cfdemCloud::forceM(int i)
@@ -468,43 +502,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;
+            return i;
-            break;
         }
     }
-    return index;
+    return -1;
-}
-
-std::vector< std::vector<double*> >* cfdemCloud::getVprobe()
-{
- return probeModel_->getVprobe();
-}
-
-std::vector< std::vector<double> >* cfdemCloud::getSprobe()
-{
- return probeModel_->getSprobe();
 }
 
 // * * * * * * * * * * * * * * * WRITE  * * * * * * * * * * * * * //
@@ -675,29 +697,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"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -160,6 +160,8 @@ protected:
 
     mutable volScalarField ddtVoidfraction_;
 
+    mutable Switch checkPeriodicCells_;
+
     const turbulenceModel& turbulence_;
 
     autoPtr<forceModel>* forceModel_;
@@ -228,41 +230,41 @@ public:
 
        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 +302,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 +324,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 +341,7 @@ public:
 
        inline bool arraysReallocated() const;
 
-       inline const wordList& forceModels();
+       inline const wordList& forceModels() const;
 
        inline const voidFractionModel& voidFractionM() const;
 
@@ -347,11 +349,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 +378,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 +396,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;
+    return vector(angularVelocities_[index][0],angularVelocities_[index][1],angularVelocities_[index][2]);
-    for(int i=0;i<3;i++) vel[i] = angularVelocities_[index][i];
-    return vel;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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(positionsCM_,0.,3,"nbodies");
-        dataExchangeM().allocateArray(velocitiesCM_,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(clumpWeights_,1.,1,"nbodies");
-        dataExchangeM().allocateArray(impForcesCM_,0,3,"nbodies");
+        dataExchangeM().allocateArray(impForcesCM_,0.,3,"nbodies");
-        dataExchangeM().allocateArray(expForcesCM_,0,3,"nbodies");
+        dataExchangeM().allocateArray(expForcesCM_,0.,3,"nbodies");
-        dataExchangeM().allocateArray(DEMForcesCM_,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/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,22 @@ void Foam::averagingModel::undoWeightFields(double**const& mask) const
     }
 }
 
-tmp<volVectorField> Foam::averagingModel::UsInterp() const
+tmp<volVectorField> averagingModel::UsInterp() const
 {
-    tmp<volVectorField> tsource
-    (
-        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
-            )
-        )
-    );
-
     if (particleCloud_.dataExchangeM().couplingStep() > 1)
     {
-        tsource.ref() = (1 - particleCloud_.dataExchangeM().timeStepFraction()) * UsPrev_
+        return tmp<volVectorField>
-                    + particleCloud_.dataExchangeM().timeStepFraction() * UsNext_;
+        (
+            new volVectorField("Us_averagingModel", (1. - particleCloud_.dataExchangeM().timeStepFraction()) * UsPrev_ + particleCloud_.dataExchangeM().timeStepFraction() * UsNext_)
+        );
     }
     else
     {
-        tsource.ref() = UsNext_;
+        return tmp<volVectorField>
+        (
+            new volVectorField("Us_averagingModel", UsNext_)
+        );
     }
-
-    return tsource;
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

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

@@ -147,11 +147,11 @@ species::species
 
 species::~species()
 {
-    delete partTemp_;
+    particleCloud_.dataExchangeM().destroy(partTemp_,1);
-    delete partRho_;
+    particleCloud_.dataExchangeM().destroy(partRho_,1);
 
-    for (int i=0; i<speciesNames_.size();i++) delete [] concentrations_[i];
+    for (int i=0; i<speciesNames_.size();i++) particleCloud_.dataExchangeM().destroy(concentrations_[i],1);
-    for (int i=0; i<speciesNames_.size();i++) delete [] changeOfSpeciesMass_[i];
+    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];
+                Yfluid_[i] = Y_[i][cellI];
-            concentrations_[i][index][0]=Yfluid_[i];
+                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;
+    int myrank = -1;
-    MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
+    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 SwapMB = SwapMem/1024.0; //kB -> MB
-    double RssMB = static_cast<double>(RssMem)/1024.0;
+    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_(time.startTime().value()+time.deltaT().value()+SMALL),  // delay start of measurement by deltaT
-    //startTime_(0),                                //no delay
+    n_(32),
-    n_(30),
     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
+        void Hist() const;                                      //calc Histogram
-        virtual void normHist() const;                      //calc normalized Histogram
+        virtual void normHist() const;                          //calc normalized Histogram
-        void plotHist(double,const std::string&,int,int) const;    //plot histogramm to terminal
+        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,7 +206,7 @@ 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_;
@@ -220,25 +217,25 @@ scalar Foam::dataExchangeModel::timeStepFraction() const
     return frac;
 }
 
-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

@@ -245,16 +245,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"));
+            partRefTemp_.value()=readScalar(propsDict_.lookup ("partRefTemp"));
-	partTempField_.writeOpt() = IOobject::AUTO_WRITE;
+        partTempField_.writeOpt() = IOobject::AUTO_WRITE;
-	partRelTempField_.writeOpt() = IOobject::AUTO_WRITE;
+        partRelTempField_.writeOpt() = IOobject::AUTO_WRITE;
-	partTempField_.write();
+        partTempField_.write();
-	partRelTempField_.write();
+        partRelTempField_.write();
-	Info <<  "Particle temperature field activated." << endl;
+        Info <<  "Particle temperature field activated." << endl;
     }
     if (verbose_)
     {
         ReField_.writeOpt() = IOobject::AUTO_WRITE;
-	NuField_.writeOpt() = IOobject::AUTO_WRITE;
+        NuField_.writeOpt() = IOobject::AUTO_WRITE;
-	ReField_.write();
+        ReField_.write();
-	NuField_.write();
+        NuField_.write();
     }
 }
 
@@ -155,10 +155,10 @@ heatTransferGunn::heatTransferGunn
 
 heatTransferGunn::~heatTransferGunn()
 {
-    delete partTemp_;
+    particleCloud_.dataExchangeM().destroy(partTemp_,1);
-    delete partHeatFlux_;
+    particleCloud_.dataExchangeM().destroy(partHeatFlux_,1);
-    delete partRe_;
+    particleCloud_.dataExchangeM().destroy(partRe_,1);
-    delete partNu_;
+    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_);
+        volScalarField sumTp (particleCloud_.averagingM().UsWeightField() * partTempField_);
-	dimensionedScalar aveTemp("aveTemp",dimensionSet(0,0,0,1,0,0,0), gSum(sumTp) / particleCloud_.numberOfParticles());
+        dimensionedScalar aveTemp("aveTemp",dimensionSet(0,0,0,1,0,0,0), gSum(sumTp) / particleCloud_.numberOfParticles());
-	partRelTempField_ = (partTempField_ - aveTemp) / (aveTemp - partRefTemp_);
+        partRelTempField_ = (partTempField_ - aveTemp) / (aveTemp - partRefTemp_);
-	Info << "heatTransferGunn: average part. temp = " << aveTemp.value() << endl;
+        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_)
+                if(verbose_)
-		{
+                {
-		    partRe_[index][0] = Rep;
+                    partRe_[index][0] = Rep;
-		    partNu_[index][0] = Nup;
+                    partNu_[index][0] = Nup;
-		}
+                }
-		
+
-		if(particleCloud_.verbose() && index >=0 && index <2)
+                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;
+        NuField_.primitiveFieldRef() = 0.0;
-	particleCloud_.averagingM().resetWeightFields();
+        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);
+    return (7 - 10 * voidfraction + 5 * voidfraction * voidfraction) *
-    Nup = (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) *
+                        (1.33 - 2.4 * voidfraction + 1.2 * voidfraction * voidfraction) *
-			Foam::pow(Rep,0.7) * Foam::pow(Pr,0.33);
+                        Foam::pow(Rep,0.7) * Foam::pow(Pr,0.33);
-                        
-    return Nup;
 }
 
 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(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(1,true); // activate treatForceDEM switch
+    forceSubM(0).setSwitchesList(SW_TREAT_FORCE_DEM,true); // activate treatForceDEM switch
     forceSubM(0).readSwitches();
 
-    if (modelType_=="A" || modelType_=="Bfull"){
+    if (modelType_=="A" || modelType_=="Bfull")
-        if(!forceSubM(0).switches()[1]) // treatDEM != true
+    {
+        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"){
+    else if (modelType_=="B")
-        if(forceSubM(0).switches()[1]) // treatDEM = true
+    {
+        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(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(2,true); // activate implDEM switch
+    forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM,true); // activate implDEM switch
-    forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch
+    forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
-    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
-    forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity 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().initialize(typeName, typeName+".logDat");
-    particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must the be the force
+    particleCloud_.probeM().vectorFields_.append("dragForce");    // first entry must the be the force
-    particleCloud_.probeM().vectorFields_.append("Urel");        //other are debug
+    particleCloud_.probeM().vectorFields_.append("Urel");         // other are debug
-    particleCloud_.probeM().scalarFields_.append("Rep");          //other are debug
+    particleCloud_.probeM().scalarFields_.append("Rep");          // other are debug
-    particleCloud_.probeM().scalarFields_.append("Cd");                 //other are debug
+    particleCloud_.probeM().scalarFields_.append("Cd");           // other are debug
-    particleCloud_.probeM().scalarFields_.append("voidfraction");       //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(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(2,true); // activate implDEM switch
+    forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM,true); // activate implDEM switch
-    forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch
+    forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
-    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
-    forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity 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     
+                                     *ds*ds
-                                     *scaleDia_ 
+                                     *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(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(2,true); // activate implDEM switch
+    forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM,true); // activate implDEM switch
-    forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch
+    forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
-    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
-    forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity 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);
+            Ufluid = vector(0,0,0);
-            voidfraction=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];
+                dSauterMix = dSauterMix_[cellI];
-		ds = 2*particleCloud_.radius(index);
+                ds = 2*particleCloud_.radius(index);
                 rho = rhoField[cellI];
                 nuf = nufField[cellI];
 
-                Rep=0.0;
+                Rep = 0.0;
-		alphaPartEff = 1.0 - voidfraction + alphaSt_[cellI] + SMALL;
+                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(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(3,true); // activate search for verbose 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 = 2 * particleCloud_.radius(index);
-		ds_scaled = ds/scaleDia_;
+                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)
+        dimensionedScalar("zero", dimensionSet(1,0,-1,0,0), 0)
-	//dimensionedVector("zero", dimensionSet(1,-2,-1,0,0), vector::zero)
+        //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
+        // pore size from hydraulic radius
-	dmean = 2 * (1 - alphaP_[cellI]) / ( (1 + poresizeWidth_*poresizeWidth_/3) * 3 * alphaP_[cellI] );
+        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
+        // Sweeney and Martin, Acta Materialia 51 (2003): ratio of hydraulic to pore throat radius
-	dmean /= ratioHydraulicPore_;
+        dmean /= ratioHydraulicPore_;
-	d1 = dmean * (1 - poresizeWidth_);
+        d1 = dmean * (1 - poresizeWidth_);
-	d2 = dmean * (1 + poresizeWidth_);
+        d2 = dmean * (1 + poresizeWidth_);
-        
+
         f = (critpore*critpore*critpore - d1 * d1 * d1) / (d2 * d2 * d2 - d1 * d1 * d1);
-	if (f<0)
+        if (f < 0)
-	{
+        {
-	    f=0.0;    
+            f = 0.0;
-	}
-	else if (f>1.0)
-	{
-	    f=1.0;
         }
-	
+        else if (f > 1.0)
-	// at this point, voidfraction is still calculated from the true particle sizes
+        {
-	deltaAlpha_[cellI] = f * (alphaMax_ - alphaP_[cellI]) - alphaSt_[cellI];
+            f = 1.0;
-	// too much volume occupied: release it (50% per time step)
+        }
-	if (deltaAlpha_[cellI] < 0.0)
+
-	{
+        // at this point, voidfraction is still calculated from the true particle sizes
-	    Sds_[cellI] = 0.5*deltaAlpha_[cellI];
+        deltaAlpha_[cellI] = f * (alphaMax_ - alphaP_[cellI]) - alphaSt_[cellI];
-	}
+        // too much volume occupied: release it (50% per time step)
-	// volume too occupy available: deposit at most 80% of dyn hold up
+        if (deltaAlpha_[cellI] < 0.0)
-	else if (depRate_ * deltaAlpha_[cellI] > 0.8 * alphaDyn_[cellI])
+        {
-	{
+            Sds_[cellI] = 0.5*deltaAlpha_[cellI];
-	    Sds_[cellI] = 0.8 * alphaDyn_[cellI];
+        }
-	}
+        // volume too occupy available: deposit at most 80% of dyn hold up
-	else
+        else if (depRate_ * deltaAlpha_[cellI] > 0.8 * alphaDyn_[cellI])
-	{
+        {
-	    Sds_[cellI] = depRate_ * deltaAlpha_[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_)
+        + fvm::div(phiSt,alphaSt_)
-	==
+        ==
-	Sds_
+        Sds_
     );
     fvScalarMatrix alphaDynEqn
     (
         fvm::ddt(alphaDyn_)
-	+ fvm::div(phiDyn,alphaDyn_)
+        + fvm::div(phiDyn,alphaDyn_)
-	- fvm::laplacian(diffCoeff_,alphaDyn_)
+        - fvm::laplacian(diffCoeff_,alphaDyn_)
-	==
+        ==
-	-Sds_
+        -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];
+            alphaStErr += alphaSt_[cellI] * particleCloud_.mesh().V()[cellI];
-	    alphaSt_[cellI] = 0.0;
+            alphaSt_[cellI] = 0.0;
-	}
+        }
-	
+
-	if (alphaDyn_[cellI] < 0.0)
+        if (alphaDyn_[cellI] < 0.0)
-	{
+        {
-	    alphaDynErr1 += alphaDyn_[cellI] * particleCloud_.mesh().V()[cellI];
+            alphaDynErr1 += alphaDyn_[cellI] * particleCloud_.mesh().V()[cellI];
-	    alphaDyn_[cellI] = 0.0;
+            alphaDyn_[cellI] = 0.0;
-	}
+        }
-	else if (alphaDyn_[cellI] > alphaDynMax_)
+        else if (alphaDyn_[cellI] > alphaDynMax_)
-	{
+        {
-	    alphaDynErr2 +=  (alphaDyn_[cellI] - alphaDynMax_) * particleCloud_.mesh().V()[cellI];
+            alphaDynErr2 +=  (alphaDyn_[cellI] - alphaDynMax_) * particleCloud_.mesh().V()[cellI];
-	    alphaDyn_[cellI] = alphaDynMax_;
+            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)
+        if(aPSt < SMALL || aPSt > 1 - SMALL)
-	  dHydMix_[cellI] = SMALL;
+            dHydMix_[cellI] = SMALL;
-	else 
+        else
-	  dHydMix_[cellI] = 2*(1 - aPSt) / (3*aPSt ) * dSauterMix_[cellI];
+            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;
+            Cd = 24.0 / Ref1;
-	else if(Ref1 <= 1000)
+        else if(Ref1 <= 1000)
-	    Cd = 24 * (1.0 + 0.15 * Foam::pow(Ref1,0.687) ) / Ref1;
+            Cd = 24 * (1.0 + 0.15 * Foam::pow(Ref1,0.687) ) / Ref1;
-	else
+        else
-	    Cd = 0.44;
+            Cd = 0.44;
-	DragCoeff_[cellI] = Cd * beta[cellI];
+        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;
+                Cd = 24.0 / Ref1;
-	    else if(Ref1 <= 1000)
+            else if(Ref1 <= 1000)
-	        Cd = 24 * (1.0 + 0.15 * Foam::pow(Ref1,0.687) ) / Ref1;
+                Cd = 24 * (1.0 + 0.15 * Foam::pow(Ref1,0.687) ) / Ref1;
-	    else
+            else
-	        Cd = 0.44;
+                Cd = 0.44;
-	    DragCoeff_.boundaryFieldRef()[patchI][faceI] = Cd * beta.boundaryFieldRef()[patchI][faceI];
+            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];
+        scalar aSt = alphaSt_[cellI];
-	if(aSt < SMALL)
+        if(aSt < SMALL)
-	    dSauterMix_[cellI] = dSauter_[cellI];
+            dSauterMix_[cellI] = dSauter_[cellI];
-	else if(aP < SMALL)
+        else if(aP < SMALL)
-	    dSauterMix_[cellI] = dFine_.value();
+            dSauterMix_[cellI] = dFine_.value();
-	else
+        else
-	    dSauterMix_[cellI] = (aP + aSt) / (aP / dSauter_[cellI] + aSt / dFine_.value() );
+            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 mU(mag(U_.boundaryField()[patchI][faceI]));
-	    scalar muDyn(mag(uDyn_.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().initialize(typeName, typeName+".logDat");
-    particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must  be the force
+    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(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(2,true); // activate implDEM switch
+    forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM,true); // activate implDEM switch
-    forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch
+    forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
-    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
-    forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity 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);
+            Ufluid = vector(0,0,0);
-            voidfraction=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 > 1.0) voidfraction = 1.0;
-                    if(voidfraction<0.10) voidfraction = 0.10;
+                    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(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(2,true); // activate implDEM switch
+    forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM,true); // activate implDEM switch
-    forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch
+    forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
-    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
-    forceSubM(0).setSwitchesList(7,true); // activate implForceDEMacc switch
+    forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM_ACCUMULATED,true); // activate implForceDEMacc switch
-    forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity 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);
+            Ufluid = vector(0,0,0);
-            voidfraction=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 > 1.0) voidfraction = 1.0;
-                    if (voidfraction < 0.40) voidfraction = 0.40;
+                    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("verbose")) verbose_ = true;
-    if (propsDict_.found("interpolation")) interpolation_=true;
+    if (propsDict_.found("interpolation")) interpolation_ = true;
-    if (propsDict_.found("randomTauE")) randomTauE_=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(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(2,true); // activate implDEM 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_;
+    particleCloud_.dataExchangeM().destroy(partTime_, 1);
-    delete partUfluct_;
+    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 > 1.0) voidfraction = 1.0;
-                    if (voidfraction < 0.40) voidfraction = 0.40;
+                    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(SW_VERBOSE,true); // activate search for verbose switch
-    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
-    forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity 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_;
+    particleCloud_.dataExchangeM().destroy(partTemp_,1);
-    delete partHeatFlux_;
+    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

@@ -73,16 +73,16 @@ 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(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch
+    forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
-    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
-    forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity 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
@@ -134,11 +134,11 @@ 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
@@ -147,10 +147,10 @@ void MeiLift::setForce() const
 
                 if( forceSubM(0).interpolation() )
                 {
-	                position       = particleCloud_.position(index);
+                    position  = particleCloud_.position(index);
-                    Ur               = UInterpolator_.interpolate(position,cellI) 
+                    Ur        = UInterpolator_.interpolate(position,cellI)
-                                        - Us;
+                                - Us;
-                    vorticity       = VorticityInterpolator_.interpolate(position,cellI);
+                    vorticity = VorticityInterpolator_.interpolate(position,cellI);
                 }
                 else
                 {
@@ -159,7 +159,7 @@ void MeiLift::setForce() const
                     vorticity=vorticityField[cellI];
                 }
 
-                magUr           = mag(Ur);
+                magUr        = mag(Ur);
                 magVorticity = mag(vorticity);
 
                 if (magUr > 0 && magVorticity > 0)
@@ -170,26 +170,26 @@ void MeiLift::setForce() const
 
                     // calc dimensionless properties
                     Rep = ds*magUr/nuf;
-		            Rew = magVorticity*ds*ds/nuf;
+                    Rew = magVorticity*ds*ds/nuf;
 
-                    alphaStar   = magVorticity*ds/magUr/2.0;
+                    alphaStar  = magVorticity*ds/magUr/2.0;
-                    epsilon       =  sqrt(2.0*alphaStar /Rep );
+                    epsilon    = sqrt(2.0*alphaStar /Rep );
-                    omega_star=2.0*alphaStar;
+                    omega_star = 2.0*alphaStar;
 
                     //Basic model for the correction to the Saffman lift
                     //Based on McLaughlin (1991)
                     if(epsilon < 0.1)
                     {
-                        J_star = -140 *epsilon*epsilon*epsilon*epsilon*epsilon 
+                        J_star = -140 *epsilon*epsilon*epsilon*epsilon*epsilon
                                              *log( 1./(epsilon*epsilon+SMALL) );
                     }
                     else if(epsilon > 20)
                     {
-                      J_star = 1.0-0.287/(epsilon*epsilon+SMALL);
+                        J_star = 1.0-0.287/(epsilon*epsilon+SMALL);
                     }
                     else
                     {
-                     J_star = 0.3
+                        J_star = 0.3
                                 *(     1.0
                                       +tanh(  2.5 * log10(epsilon+0.191)  )
                                  )
@@ -197,11 +197,11 @@ void MeiLift::setForce() const
                                      +tanh(  6.0 * (epsilon-0.32)  )
                                   );
                     }
-                    Cl=J_star*4.11*epsilon; //multiply McLaughlin's correction to the basic Saffman model
+                    Cl = J_star * 4.11 * epsilon; //multiply McLaughlin's correction to the basic Saffman model
 
-                    //Second order terms given by Loth and Dorgan 2009 
+                    //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;
@@ -209,21 +209,21 @@ void MeiLift::setForce() const
 
                     lift =  0.125*M_PI
                            *rho
-                           *Cl  
+                           *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 +240,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 +252,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/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();

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

@@ -72,7 +72,7 @@ ShirgaonkarIB::ShirgaonkarIB
     p_(sm.mesh().lookupObject<volScalarField> (pressureFieldName_))
 {
     //Append the field names to be probed
-    particleCloud_.probeM().initialize(typeName, "shirgaonkarIB.logDat");
+    particleCloud_.probeM().initialize(typeName, typeName+".logDat");
     particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must the be the force
     particleCloud_.probeM().writeHeader();
 
@@ -89,7 +89,7 @@ ShirgaonkarIB::ShirgaonkarIB
     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();

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

@@ -92,20 +92,20 @@ dSauter::dSauter
         ),
         sm.mesh(),
         dimensionedScalar("zero", dimensionSet(0,1,0,0,0), 0),
-	"zeroGradient"
+        "zeroGradient"
     )
 {
     allocateMyArrays();
     dSauter_.write();
-    
+
 
     // init force sub model
     setForceSubModels(propsDict_);
 
     if (propsDict_.found("scaleCG"))
-        scaleDia_=scalar(readScalar(propsDict_.lookup("scaleCG")));
+        scaleDia_ = scalar(readScalar(propsDict_.lookup("scaleCG")));
     if (propsDict_.found("scaleDist"))
-        scaleDiaDist_=scalar(readScalar(propsDict_.lookup("scaleDist")));
+        scaleDiaDist_ = scalar(readScalar(propsDict_.lookup("scaleDist")));
 }
 
 
@@ -113,27 +113,31 @@ dSauter::dSauter
 
 dSauter::~dSauter()
 {
-    delete d2_;
+    particleCloud_.dataExchangeM().destroy(d2_,1);
-    delete d3_;
+    particleCloud_.dataExchangeM().destroy(d3_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+
 void dSauter::allocateMyArrays() const
 {
     // get memory for 2d arrays
-    double initVal=0.0;
+    double initVal = 0.0;
     particleCloud_.dataExchangeM().allocateArray(d2_,initVal,1);  // field/initVal/with/lenghtFromLigghts
     particleCloud_.dataExchangeM().allocateArray(d3_,initVal,1);
 }
+
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
 
 void dSauter::setForce() const
 {
     if (scaleDia_ > 1)
+    {
         Info << "dSauter using scaleCG = " << scaleDia_ << endl;
+    }
     else if (particleCloud_.cg() > 1)
     {
-        scaleDia_=particleCloud_.cg();
+        scaleDia_ = particleCloud_.cg();
         Info << "dSauter using scaleCG from liggghts cg = " << scaleDia_ << endl;
     }
 
@@ -142,21 +146,21 @@ void dSauter::setForce() const
     label cellI=0;
     scalar ds(0);
     scalar scale = scaleDiaDist_/scaleDia_;
-    
+
-    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
         cellI = particleCloud_.cellIDs()[index][0];
-        if(cellI >= 0)
+        if (cellI >= 0)
         {
-	    ds = particleCloud_.d(index);
+            ds = particleCloud_.d(index);
-	    d2_[index][0] = ds*ds;
+            d2_[index][0] = ds*ds;
-	    d3_[index][0] = ds*ds*ds;
+            d3_[index][0] = ds*ds*ds;
-	}
+        }
     }
-    
+
     d2Field_.primitiveFieldRef() = 0.0;
     d3Field_.primitiveFieldRef() = 0.0;
-    
+
     particleCloud_.averagingM().setScalarSum
     (
         d2Field_,
@@ -164,7 +168,7 @@ void dSauter::setForce() const
         particleCloud_.particleWeights(),
         NULL
     );
-    
+
     particleCloud_.averagingM().setScalarSum
     (
         d3Field_,
@@ -172,19 +176,19 @@ void dSauter::setForce() const
         particleCloud_.particleWeights(),
         NULL
     );
-    
+
     forAll(dSauter_,cellI)
     {
-        if(d2Field_[cellI] > ROOTVSMALL)
+        if (d2Field_[cellI] > ROOTVSMALL)
         {
-             dSauter_[cellI] = d3Field_[cellI] / d2Field_[cellI] * scale;
+            dSauter_[cellI] = d3Field_[cellI] / d2Field_[cellI] * scale;
         }
         else
-	{
+        {
-	  dSauter_[cellI] = SMALL;
+            dSauter_[cellI] = SMALL;
-	}
+        }
     }
-    
+
     dSauter_.correctBoundaryConditions();
 }
 

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

@@ -88,14 +88,16 @@ forceModel::forceModel
     probeIt_(sm.probeM().active()),
     requiresEx_(false),
     forceSubModels_(0),
-    forceSubModel_(new autoPtr<forceSubModel>[nrForceSubModels()])
+    forceSubModel_(NULL)
 {}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 forceModel::~forceModel()
-{}
+{
+  delete [] forceSubModel_;
+}
 
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
 /*tmp<volScalarField> forceModel::provideScalarField()

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

@@ -131,7 +131,7 @@ public:
         virtual void manipulateScalarField(volScalarField&) const;
 
     // Access
-        word modelType() { return modelType_; }
+        const word& modelType() const { return modelType_; }
 
         inline volVectorField& impParticleForces() const { return impParticleForces_; }
 
@@ -147,7 +147,7 @@ public:
 
         inline double ** fluidVel() const { return particleCloud_.fluidVel_; }
 
-        inline const bool& coupleForce() const { return coupleForce_; }
+        inline bool coupleForce() const { return coupleForce_; }
 
         virtual inline bool& requiresEx() { return requiresEx_; }
 
@@ -155,11 +155,11 @@ public:
 
         void treatVoidCells() const;
 
-        inline const wordList& forceSubModels(){ return forceSubModels_; }
+        inline const wordList& forceSubModels() { return forceSubModels_; }
 
-        inline const forceSubModel& forceSubM(int i) const { return forceSubModel_[i]; }
+        inline forceSubModel& forceSubM(int i) const { return forceSubModel_[i](); }
 
-        inline int nrForceSubModels(){ return forceSubModels_.size(); }
+        inline int nrForceSubModels() const { return forceSubModels_.size(); }
 
         void setForceSubModels(dictionary& dict);
 };

src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/ScaleForce/ScaleForce.C

@@ -62,7 +62,7 @@ ScaleForce::~ScaleForce()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 void ScaleForce::partToArray
 (
-    label& index,
+    label index,
     vector& dragTot,
     const vector& dragEx,
     const vector& Ufluid,

src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/ScaleForce/ScaleForce.H

@@ -41,9 +41,9 @@ class ScaleForce
     public forceSubModel
 {
 private:
-  
+
     dictionary propsDict_;
-  
+
     word scaleFieldName_;
 
     const volScalarField& scaleField_;
@@ -70,8 +70,8 @@ public:
 
 
     // Member Functions
-	
+
-	void partToArray(label&, vector&, const vector&, const vector& Ufluid=vector::zero, scalar Cd=scalar(0)) const;
+        void partToArray(label, vector&, const vector&, const vector& Ufluid=vector::zero, scalar Cd=scalar(0)) const;
 
         word myType() const{return typeName; };
 };

src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/forceSubModel/forceSubModel.C

@@ -59,7 +59,7 @@ forceSubModel::forceSubModel
     dict_(dict),
     particleCloud_(sm),
     forceModel_(fm),
-    nrDefaultSwitches_(9),                                          // !!!
+    nrDefaultSwitches_(SW_MAX),
     switchesNameList_(nrDefaultSwitches_),
     switchesList_(nrDefaultSwitches_),
     switches_(nrDefaultSwitches_),
@@ -106,16 +106,15 @@ forceSubModel::forceSubModel
     rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_))
 {
     // init standard switch list
-    int iCounter(0);
+    switchesNameList_[SW_TREAT_FORCE_EXPLICIT] = "treatForceExplicit";
-    switchesNameList_[iCounter]="treatForceExplicit"; iCounter++;   //0
+    switchesNameList_[SW_TREAT_FORCE_DEM] = "treatForceDEM";
-    switchesNameList_[iCounter]="treatForceDEM";iCounter++;         //1
+    switchesNameList_[SW_IMPL_FORCE_DEM] = "implForceDEM";
-    switchesNameList_[iCounter]="implForceDEM";iCounter++;          //2
+    switchesNameList_[SW_VERBOSE] = "verbose";
-    switchesNameList_[iCounter]="verbose";iCounter++;               //3
+    switchesNameList_[SW_INTERPOLATION] = "interpolation";
-    switchesNameList_[iCounter]="interpolation";iCounter++;         //4
+    switchesNameList_[SW_FILTERED_DRAG_MODEL] ="useFilteredDragModel";
-    switchesNameList_[iCounter]="useFilteredDragModel";iCounter++;  //5
+    switchesNameList_[SW_PARCEL_SIZE_DEPENDENT_FILTERED_DRAG] = "useParcelSizeDependentFilteredDrag";
-    switchesNameList_[iCounter]="useParcelSizeDependentFilteredDrag";iCounter++;  //6
+    switchesNameList_[SW_IMPL_FORCE_DEM_ACCUMULATED] = "implForceDEMaccumulated";
-    switchesNameList_[iCounter]="implForceDEMaccumulated";iCounter++;             //7
+    switchesNameList_[SW_SCALAR_VISCOSITY] = "scalarViscosity";
-    switchesNameList_[iCounter]="scalarViscosity";iCounter++;                     //8
 
     // sanity check of what is defined above
     if(switchesNameList_.size() != nrDefaultSwitches_)
@@ -131,7 +130,7 @@ forceSubModel::~forceSubModel()
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
 void forceSubModel::partToArray
 (
-    label& index,
+    label index,
     vector& dragTot,
     const vector& dragEx,
     const vector& Ufluid,
@@ -139,17 +138,17 @@ void forceSubModel::partToArray
 ) const
 {
     // forces for CFD
-    if(!switches_[1])// !treatDEM
+    if(!switches_[SW_TREAT_FORCE_DEM])// !treatDEM
     {
-        if(switches_[0]) // treatExplicit
+        if(switches_[SW_TREAT_FORCE_EXPLICIT]) // treatExplicit
         {
             for(int j=0;j<3;j++)
                 myForceM().expForces()[index][j] += dragTot[j];
-        }    
+        }
         else   //implicit treatment, taking explicit force contribution into account
         {
-            for(int j=0;j<3;j++) 
+            for(int j=0;j<3;j++)
-            { 
+            {
                 myForceM().impForces()[index][j] += dragTot[j] - dragEx[j]; //only consider implicit part!
                 myForceM().expForces()[index][j] += dragEx[j];
             }
@@ -157,16 +156,16 @@ void forceSubModel::partToArray
     }
 
     // forces for DEM
-    if(switches_[2]) // implForceDEM
+    if(switches_[SW_IMPL_FORCE_DEM]) // implForceDEM
     {
         for(int j=0;j<3;j++)
             myForceM().fluidVel()[index][j]=Ufluid[j];
 
-        myForceM().Cds()[index][0]=Cd;
+        myForceM().Cds()[index][0] = Cd;
     }
     else
     {
-        for(int j=0;j<3;j++) 
+        for(int j=0;j<3;j++)
             myForceM().DEMForces()[index][j] += dragTot[j];
     }
 }
@@ -183,7 +182,7 @@ void forceSubModel::explicitCorr
     vector& Us,
     const vector& UsCell,
     bool verbose,
-    label index    
+    label index
 ) const
 {
     dragExplicit=vector::zero;
@@ -191,7 +190,7 @@ void forceSubModel::explicitCorr
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-void forceSubModel::readSwitches() const
+void forceSubModel::readSwitches()
 {
     Info << "\nreading switches for forceSubModel:" << myType() << endl;
     forAll(switchesNameList_,i)
@@ -200,23 +199,23 @@ void forceSubModel::readSwitches() const
         {
             Info << "  looking for " << switchesNameList_[i] << " ..." << endl;
             if (dict_.found(switchesNameList_[i]))
-                switches_[i]=Switch(dict_.lookup(switchesNameList_[i]));
+                switches_[i] = Switch(dict_.lookup(switchesNameList_[i]));
-                
+
             Info << "\t" << switchesNameList_[i] << " = " << switches_[i] << endl;
-        }        
+        }
     }
     Info << endl;
 
-    if(switches_[2]) // implForceDEM=true
+    if(switches_[SW_IMPL_FORCE_DEM]) // implForceDEM=true
     {
         // communicate implForceDEM to particleCloud
-        particleCloud_.impDEMdrag_=true;
+        particleCloud_.impDEMdrag_ = true;
 
         // do sanity check
         // This can work if the accumulator is used, but is explicitely applied on the CFD side
         // Sanity check is therefore not necessary here
         /*
-        if(switches_[0]) // treatExplicit=true
+        if(switches_[SW_TREAT_FORCE_EXPLICIT]) // treatExplicit=true
         {
             FatalError << "Please check your settings, treatExplicit together with implForceDEM does not work!."
                        << abort(FatalError);
@@ -224,24 +223,24 @@ void forceSubModel::readSwitches() const
         */
     }
 
-    if(switches_[7]) // implForceDEMaccumulated=true
+    if(switches_[SW_IMPL_FORCE_DEM_ACCUMULATED]) // implForceDEMaccumulated=true
     {
         // sanity check for implForceDEMaccumulated
-        if(!switches_[2]) //implForceDEM=false
+        if(!switches_[SW_IMPL_FORCE_DEM]) //implForceDEM=false
         {
             Warning<< "please check your settings, implForceDEMaccumulated without implForceDEM does not work! (using implForceDEMaccumulated=false)" << endl;
-            switches_[3]=false;
+            switches_[SW_VERBOSE] = false;
         }else
         {
-            particleCloud_.impDEMdragAcc_=true;
+            particleCloud_.impDEMdragAcc_ = true;
         }
     }
 
-    if(switches_[8]) // scalarViscosity=true
+    if(switches_[SW_SCALAR_VISCOSITY]) // scalarViscosity=true
     {
         Info << "Using a constant viscosity for this force model." << endl;
         dimensionedScalar  nu0_("nu", dimensionSet(0, 2, -1, 0, 0), dict_.lookup("nu"));
-        nu_=volScalarField
+        nu_ = volScalarField
         (
             IOobject
             (
@@ -259,7 +258,7 @@ void forceSubModel::readSwitches() const
     // look for old nomenclature
     if (dict_.found("treatExplicit") || dict_.found("treatDEM") || dict_.found("implDEM"))
         FatalError<< "You are using an old nomenclature for force model settings, please have a look at the forceSubModel doc." << abort(FatalError);
-        
+
     // look for old nomenclature
     if (dict_.found("verbose"))
         Warning<< "Please make sure you use the new nomenclature for verbose force model settings, please have a look at the forceSubModel doc." << endl;
@@ -274,7 +273,7 @@ const volScalarField& forceSubModel::nuField() const
         nu_=particleCloud_.turbulence().mu() / rho_;
         return nu_;
     #else
-        if(switches_[8]) // scalarViscosity=true
+        if(switches_[SW_SCALAR_VISCOSITY]) // scalarViscosity=true
             return nu_;
         else
             return particleCloud_.turbulence().nu();
@@ -287,9 +286,9 @@ const volScalarField& forceSubModel::muField() const
         return particleCloud_.turbulence().mu();
     #else
         // passing the ref to nu*rho will not work->generate a mu_ field like nu_
-        FatalError<< "implementation not complete!" << abort(FatalError);
+        FatalError << "implementation not complete!" << abort(FatalError);
 
-        if(switches_[8]) // scalarViscosity=true
+        if(switches_[SW_SCALAR_VISCOSITY]) // scalarViscosity=true
             return nu_*rho_;
         else
             return particleCloud_.turbulence().nu()*rho_;

src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/forceSubModel/forceSubModel.H

@@ -47,6 +47,19 @@ SourceFiles
 namespace Foam
 {
 
+enum {
+    SW_TREAT_FORCE_EXPLICIT = 0,
+    SW_TREAT_FORCE_DEM,
+    SW_IMPL_FORCE_DEM,
+    SW_VERBOSE,
+    SW_INTERPOLATION,
+    SW_FILTERED_DRAG_MODEL,
+    SW_PARCEL_SIZE_DEPENDENT_FILTERED_DRAG,
+    SW_IMPL_FORCE_DEM_ACCUMULATED,
+    SW_SCALAR_VISCOSITY,
+    SW_MAX
+};
+
 /*---------------------------------------------------------------------------*\
                            Class forceSubModel Declaration
 \*---------------------------------------------------------------------------*/
@@ -67,9 +80,9 @@ protected:
 
         wordList switchesNameList_;                 // names of switches available
 
-        mutable List<Switch> switchesList_;         // switches which are requested in dict
+        List<Switch> switchesList_;         // switches which are requested from dict
 
-        mutable List<Switch> switches_;
+        List<Switch> switches_;             // switch status
 
         mutable volScalarField nu_;
 
@@ -130,35 +143,33 @@ public:
 
 
     // Member Functions
-        virtual void partToArray(label&, vector&, const vector&, const vector& Ufluid=vector::zero, scalar Cd=scalar(0)) const;
+        virtual void partToArray(label, vector&, const vector&, const vector& Ufluid=vector::zero, scalar Cd=scalar(0)) const;
 
         virtual void explicitCorr(vector&, vector&, scalar&, vector&, const vector&, vector&, const vector&, bool,label index=100) const;
 
     // Access
 
-        inline bool verbose() const      { return switches_[3]; }
+        inline bool verbose() const      { return switches_[SW_VERBOSE]; }
 
-        inline bool interpolation() const { return switches_[4]; }
+        inline bool interpolation() const { return switches_[SW_INTERPOLATION]; }
 
-        inline bool useFilteredDragModel() const { return switches_[5]; }
+        inline bool useFilteredDragModel() const { return switches_[SW_FILTERED_DRAG_MODEL]; }
 
-        inline bool useParcelSizeDependentFilteredDrag() const { return switches_[6]; }
+        inline bool useParcelSizeDependentFilteredDrag() const { return switches_[SW_PARCEL_SIZE_DEPENDENT_FILTERED_DRAG]; }
 
-        virtual word myType() const=0;
+        virtual word myType() const = 0;
 
-        inline forceModel& myForceM() const { return forceModel_;}
+        inline forceModel& myForceM() const { return forceModel_; }
 
-        inline const List<Switch>& switches() const { return switches_;}
+        inline const List<Switch>& switches() const { return switches_; }
 
-        inline const wordList& switchesNameList() const { return switchesNameList_;}
+        inline const wordList& switchesNameList() const { return switchesNameList_; }
 
-        void setSwitchesList(label i, bool v) const { switchesList_[i] = v;}
+        void setSwitchesList(label i, Switch v) { switchesList_[i] = v; }
 
-        void setSwitches(label i, Switch v) const { switches_[i] = v;}
+        void setSwitches(label i, Switch v) { switches_[i] = v; }
 
-        virtual void readSwitches() const;
+        virtual void readSwitches();
-
-        const label& nrDefaultSwitches() const {return nrDefaultSwitches_;}
 
         const volScalarField& nuField() const;
 

src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/scaleForceBoundary/scaleForceBoundary.C

@@ -56,35 +56,35 @@ scaleForceBoundary::scaleForceBoundary
     if (propsDict_.found("x1") && propsDict_.found("x2"))
     {
         coordinateInner_=readScalar(propsDict_.lookup ("x1"));
-	coordinateOuter_=readScalar(propsDict_.lookup ("x2"));
+        coordinateOuter_=readScalar(propsDict_.lookup ("x2"));
-	dim_ = 0;
+        dim_ = 0;
         Info << "scaleForceBoundary: Limiting force in x direction."  << endl;
     }
     else if (propsDict_.found("y1") && propsDict_.found("y2"))
     {
         coordinateInner_=readScalar(propsDict_.lookup ("y1"));
-	coordinateOuter_=readScalar(propsDict_.lookup ("y2"));
+        coordinateOuter_=readScalar(propsDict_.lookup ("y2"));
-	dim_ = 1;
+        dim_ = 1;
         Info << "scaleForceBoundary: Limiting force in y direction."  << endl;
     }
     else if (propsDict_.found("z1") && propsDict_.found("z2"))
     {
         coordinateInner_=readScalar(propsDict_.lookup ("z1"));
-	coordinateOuter_=readScalar(propsDict_.lookup ("z2"));
+        coordinateOuter_=readScalar(propsDict_.lookup ("z2"));
-	dim_ = 2;
+        dim_ = 2;
         Info << "scaleForceBoundary: Limiting force in z direction."  << endl;
     }
-    
+
     if (propsDict_.found("outerValue"))
         outerVal_=readScalar(propsDict_.lookup ("outerValue"));
-    
+
     if(coordinateOuter_ > coordinateInner_)
         orientation_ = 1;
     else
         orientation_ = -1;
-  
+
     dist_ = fabs(coordinateOuter_ - coordinateInner_);
-  
+
 }
 
 
@@ -97,7 +97,7 @@ scaleForceBoundary::~scaleForceBoundary()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 void scaleForceBoundary::partToArray
 (
-    label& index,
+    label index,
     vector& dragTot,
     const vector& dragEx,
     const vector& Ufluid,

src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/scaleForceBoundary/scaleForceBoundary.H

@@ -45,19 +45,19 @@ class scaleForceBoundary
     public forceSubModel
 {
 private:
-  
+
     dictionary propsDict_;
-  
+
     scalar coordinateInner_;
-    
+
     scalar coordinateOuter_;
-    
+
     scalar dist_;
-    
+
     scalar outerVal_;
-    
+
     int dim_;
-    
+
     int orientation_;
 
 public:
@@ -82,10 +82,10 @@ public:
 
 
     // Member Functions
-	
-	void partToArray(label&, vector&, const vector&, const vector& Ufluid=vector::zero, scalar Cd=scalar(0)) const;
 
-        word myType() const{return typeName; };
+        void partToArray(label, vector&, const vector&, const vector& Ufluid=vector::zero, scalar Cd=scalar(0)) const;
+
+        word myType() const {return typeName; }
 };
 
 

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

@@ -74,9 +74,9 @@ gradPForce::gradPForce
     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
-    forceSubM(0).setSwitchesList(1,true); // activate treatForceDEM switch
+    forceSubM(0).setSwitchesList(SW_TREAT_FORCE_DEM,true); // activate treatForceDEM switch
-    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
 
     // read those switches defined above, if provided in dict
     forceSubM(0).readSwitches();
@@ -84,26 +84,30 @@ gradPForce::gradPForce
     if (modelType_ == "B")
     {
         FatalError <<"using  model gradPForce with model type B is not valid\n" << abort(FatalError);
-    }else if (modelType_ == "Bfull")
+    }
+    else if (modelType_ == "Bfull")
     {
-        if(forceSubM(0).switches()[1])
+        if(forceSubM(0).switches()[SW_TREAT_FORCE_DEM])
         {
             Info << "Using treatForceDEM false!" << endl;
-            forceSubM(0).setSwitches(1,false); // treatForceDEM = false
+            forceSubM(0).setSwitches(SW_TREAT_FORCE_DEM,false); // treatForceDEM = false
         }
-    }else // modelType_=="A"
+    }
+    else // modelType_=="A"
     {
-        if(!forceSubM(0).switches()[1])
+        if(!forceSubM(0).switches()[SW_TREAT_FORCE_DEM])
         {
             Info << "Using treatForceDEM true!" << endl;
-            forceSubM(0).setSwitches(1,true); // treatForceDEM = true
+            forceSubM(0).setSwitches(SW_TREAT_FORCE_DEM,true); // treatForceDEM = true
         }
     }
 
-    if (propsDict_.found("useU")) useU_=true;
+    if (propsDict_.found("useU"))
-    if (propsDict_.found("useAddedMass")) 
+        useU_ = true;
+
+    if (propsDict_.found("useAddedMass"))
     {
-        addedMassCoeff_ =  readScalar(propsDict_.lookup("useAddedMass"));
+        addedMassCoeff_ = readScalar(propsDict_.lookup("useAddedMass"));
         Info << "gradP will also include added mass with coefficient: " << addedMassCoeff_ << endl;
         Info << "WARNING: use fix nve/sphere/addedMass in LIGGGHTS input script to correctly account for added mass effects!" << endl;
     }
@@ -113,7 +117,7 @@ gradPForce::gradPForce
 
     particleCloud_.checkCG(true);
 
-    particleCloud_.probeM().initialize(typeName, "gradP.logDat");
+    particleCloud_.probeM().initialize(typeName, typeName+".logDat");
     particleCloud_.probeM().vectorFields_.append("gradPForce"); //first entry must the be the force
     particleCloud_.probeM().scalarFields_.append("Vs");
     particleCloud_.probeM().scalarFields_.append("rho");
@@ -152,21 +156,22 @@ void gradPForce::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])
         //{
-            force=vector(0,0,0);
+            force = vector(0,0,0);
             cellI = particleCloud_.cellIDs()[index][0];
 
             if (cellI > -1) // particle Found
             {
                 position = particleCloud_.position(index);
 
-                if(forceSubM(0).interpolation()) // use intepolated values for alpha (normally off!!!)
+                if (forceSubM(0).interpolation()) // use intepolated values for alpha (normally off!!!)
                 {
                     gradP = gradPInterpolator_.interpolate(position,cellI);
-                }else
+                }
+                else
                 {
                     gradP = gradPField[cellI];
                 }
@@ -180,7 +185,7 @@ void gradPForce::setForce() const
                 else
                     force = -Vs*gradP*(1.0+addedMassCoeff_);
 
-                if(forceSubM(0).verbose() && index >=0 && index <2)
+                if (forceSubM(0).verbose() && index >= 0 && index < 2)
                 {
                     Info << "index = " << index << endl;
                     Info << "gradP = " << gradP << endl;
@@ -188,7 +193,7 @@ void gradPForce::setForce() const
                 }
 
                 //Set value fields and write the probe
-                if(probeIt_)
+                if (probeIt_)
                 {
                     #include "setupProbeModelfields.H"
                     vValues.append(force);           //first entry must the be the force

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

@@ -67,12 +67,12 @@ granKineticEnergy::granKineticEnergy
         ),
         sm.mesh(),
         dimensionedScalar("zero", dimensionSet(0,2,-2,0,0), 0),
-	"zeroGradient"
+        "zeroGradient"
     )
 {
     allocateMyArrays();
     granKineticEnergy_.write();
-    
+
 
     // init force sub model
     setForceSubModels(propsDict_);
@@ -82,14 +82,14 @@ granKineticEnergy::granKineticEnergy
 
 granKineticEnergy::~granKineticEnergy()
 {
-    delete vfluc_;
+    particleCloud_.dataExchangeM().destroy(vfluc_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
 void granKineticEnergy::allocateMyArrays() const
 {
     // get memory for 2d arrays
-    double initVal=0.0;
+    double initVal = 0.0;
     particleCloud_.dataExchangeM().allocateArray(vfluc_,initVal,1);
 }
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
@@ -98,20 +98,20 @@ void granKineticEnergy::setForce() const
 {
     allocateMyArrays();
 
-    label cellI=0;
+    label cellI = 0;
     vector velfluc(0,0,0);
 
-    
+
-    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
         cellI = particleCloud_.cellIDs()[index][0];
         if(cellI >= 0)
         {
-	    velfluc = particleCloud_.velocity(index) - UsField_[cellI];
+            velfluc = particleCloud_.velocity(index) - UsField_[cellI];
-	    vfluc_[index][0] = magSqr(velfluc);
+            vfluc_[index][0] = magSqr(velfluc);
-	}
+        }
     }
-    
+
     granKineticEnergy_.primitiveFieldRef() = 0.0;
 
     particleCloud_.averagingM().resetWeightFields();
@@ -120,10 +120,10 @@ void granKineticEnergy::setForce() const
         granKineticEnergy_,
         vfluc_,
         particleCloud_.particleWeights(),
-	particleCloud_.averagingM().UsWeightField(),
+        particleCloud_.averagingM().UsWeightField(),
         NULL
     );
-    
+
     granKineticEnergy_ *= 0.5;
 
 }

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

@@ -83,7 +83,7 @@ interface::interface
     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();

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

@@ -72,7 +72,7 @@ noDrag::noDrag
     setForceSubModels(propsDict_);
 
     // define switches which can be read from dict
-    forceSubM(0).setSwitchesList(0,true); // activate treatExplicit
+    forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit
 
     // read those switches defined above, if provided in dict
     forceSubM(0).readSwitches();
@@ -98,7 +98,7 @@ void noDrag::setForce() const
     // Do nothing
     Info << "noDrag::setForce" << endl;
     label cellI=0;
-    bool treatExplicit=forceSubM(0).switches()[0];
+    bool treatExplicit = forceSubM(0).switches()[SW_TREAT_FORCE_EXPLICIT];
     for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
     {
         cellI = particleCloud_.cellIDs()[index][0];
@@ -125,7 +125,7 @@ void noDrag::setForce() const
                 }
             }
             //==========================
-        }        
+        }
     }
 }
 

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

@@ -83,28 +83,30 @@ virtualMassForce::virtualMassForce
     // 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(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
     forceSubM(0).readSwitches();
 
     //Extra switches/settings
     if(propsDict_.found("splitUrelCalculation"))
     {
-        splitUrelCalculation_ = readBool(propsDict_.lookup("splitUrelCalculation")); 
+        splitUrelCalculation_ = readBool(propsDict_.lookup("splitUrelCalculation"));
         if(splitUrelCalculation_)
+        {
             Info << "Virtual mass model: will split the Urel calculation\n";
             Info << "WARNING: be sure that LIGGGHTS integration takes ddtv_p implicitly into account! \n";
+        }
     }
     if(propsDict_.found("Cadd"))
     {
-        Cadd_ = readScalar(propsDict_.lookup("Cadd")); 
+        Cadd_ = readScalar(propsDict_.lookup("Cadd"));
         Info << "Virtual mass model: using non-standard Cadd = " << Cadd_ << endl;
     }
 
     particleCloud_.checkCG(true);
 
     //Append the field names to be probed
-    particleCloud_.probeM().initialize(typeName, "virtualMass.logDat");
+    particleCloud_.probeM().initialize(typeName, typeName+".logDat");
     particleCloud_.probeM().vectorFields_.append("virtualMassForce"); //first entry must the be the force
     particleCloud_.probeM().vectorFields_.append("Urel");
     particleCloud_.probeM().vectorFields_.append("UrelOld");
@@ -119,7 +121,7 @@ virtualMassForce::virtualMassForce
 
 virtualMassForce::~virtualMassForce()
 {
-    delete UrelOld_;
+    particleCloud_.dataExchangeM().destroy(UrelOld_,3);
 }
 
 
@@ -146,7 +148,7 @@ void virtualMassForce::setForce() const
 
     #include "setupProbeModel.H"
 
-    bool haveUrelOld_(false); 
+    bool haveUrelOld_(false);
 
     for(int index = 0;index <  particleCloud_.numberOfParticles(); index++)
     {
@@ -156,19 +158,19 @@ void virtualMassForce::setForce() const
             if (cellI > -1) // particle Found
             {
 
-                if(forceSubM(0).interpolation()) 
+                if(forceSubM(0).interpolation())
                 {
-	                position     = particleCloud_.position(index);
+                    position    = particleCloud_.position(index);
-                    Ufluid       = UInterpolator_.interpolate(position,cellI);
+                    Ufluid      = UInterpolator_.interpolate(position,cellI);
                 }
                 else
                 {
                     Ufluid = U_[cellI];
                 }
-           
+
 
                 if(splitUrelCalculation_)  //if split, just use total derivative of fluid velocity
-                if(forceSubM(0).interpolation()) 
+                if(forceSubM(0).interpolation())
                 {
                     DDtU = DDtUInterpolator_.interpolate(position,cellI);
                 }
@@ -182,7 +184,7 @@ void virtualMassForce::setForce() const
                     Ur = Ufluid - Us;
                 }
 
-                
+
                 //Check of particle was on this CPU the last step
                 if(UrelOld_[index][0]==NOTONCPU) //use 1. element to indicate that particle was on this CPU the last time step
                     haveUrelOld_ = false;
@@ -229,7 +231,7 @@ void virtualMassForce::setForce() const
             }
             else    //particle not on this CPU
                 UrelOld_[index][0]=NOTONCPU;
-    
+
             // write particle based data to global array
             forceSubM(0).partToArray(index,virtualMassForce,vector::zero);
     }

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

@@ -71,10 +71,10 @@ viscForce::viscForce
     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
-    forceSubM(0).setSwitchesList(1,true); // activate treatForceDEM switch
+    forceSubM(0).setSwitchesList(SW_TREAT_FORCE_DEM,true); // activate treatForceDEM switch
-    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
-    forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity switch
+    forceSubM(0).setSwitchesList(SW_SCALAR_VISCOSITY,true); // activate scalarViscosity switch
 
     // read those switches defined above, if provided in dict
     forceSubM(0).readSwitches();
@@ -84,22 +84,22 @@ viscForce::viscForce
         FatalError <<"using  model viscForce with model type B is not valid\n" << abort(FatalError);
     }else if (modelType_ == "Bfull")
     {
-        if(forceSubM(0).switches()[1])
+        if(forceSubM(0).switches()[SW_TREAT_FORCE_DEM])
         {
             Info << "Using treatForceDEM false!" << endl;
-            forceSubM(0).setSwitches(1,false); // treatForceDEM = false
+            forceSubM(0).setSwitches(SW_TREAT_FORCE_DEM,false); // treatForceDEM = false
         }
 
     }else // modelType_=="A"
     {
-        if(!forceSubM(0).switches()[1])
+        if(!forceSubM(0).switches()[SW_TREAT_FORCE_DEM])
         {
             Info << "Using treatForceDEM true!" << endl;
-            forceSubM(0).setSwitches(1,true); // treatForceDEM = true
+            forceSubM(0).setSwitches(SW_TREAT_FORCE_DEM,true); // treatForceDEM = true
         }
     }
 
-    if (propsDict_.found("useAddedMass")) 
+    if (propsDict_.found("useAddedMass"))
     {
         addedMassCoeff_ =  readScalar(propsDict_.lookup("useAddedMass"));
         Info << "viscForce will also include added mass with coefficient: " << addedMassCoeff_ << endl;
@@ -109,7 +109,7 @@ viscForce::viscForce
     particleCloud_.checkCG(true);
 
     //Append the field names to be probed
-    particleCloud_.probeM().initialize(typeName, "visc.logDat");
+    particleCloud_.probeM().initialize(typeName, typeName+".logDat");
     particleCloud_.probeM().vectorFields_.append("viscForce"); //first entry must the be the force
     particleCloud_.probeM().scalarFields_.append("Vs");
     particleCloud_.probeM().writeHeader();
@@ -160,7 +160,7 @@ void viscForce::setForce() const
 
                 Vs = particleCloud_.particleVolume(index);
 
-                // calc the contribution of the deviatoric stress 
+                // calc the contribution of the deviatoric stress
                 // to the generalized buoyancy force
                 force = -Vs*divTau*(1.0+addedMassCoeff_);
 

src/lagrangian/cfdemParticle/subModels/forceModelMS/DiFeliceDragMS/DiFeliceDragMS.C

@@ -75,21 +75,21 @@ DiFeliceDragMS::DiFeliceDragMS
     //dH_(readScalar(propsDict_.lookup("hydraulicDiameter")))
 {
     //Append the field names to be probed
-    particleCloud_.probeM().initialize(typeName, "diFeliceDrag.logDat");
+    particleCloud_.probeM().initialize(typeName, typeName+".logDat");
-    particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must the be the force
+    particleCloud_.probeM().vectorFields_.append("dragForce");    // first entry must the be the force
-    particleCloud_.probeM().vectorFields_.append("Urel");        //other are debug
+    particleCloud_.probeM().vectorFields_.append("Urel");         // other are debug
-    particleCloud_.probeM().scalarFields_.append("Rep");          //other are debug
+    particleCloud_.probeM().scalarFields_.append("Rep");          // other are debug
-    particleCloud_.probeM().scalarFields_.append("Cd");                 //other are debug
+    particleCloud_.probeM().scalarFields_.append("Cd");           // other are debug
-    particleCloud_.probeM().scalarFields_.append("voidfraction");       //other are debug
+    particleCloud_.probeM().scalarFields_.append("voidfraction"); // other are debug
     particleCloud_.probeM().writeHeader();
 
     // 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(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
-    forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch
+    forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
-    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
 
     // read those switches defined above, if provided in dict
     forceSubM(0).readSwitches();
@@ -98,11 +98,12 @@ DiFeliceDragMS::DiFeliceDragMS
     {
         Warning << " interpolation is commented for this force model - it seems to be unstable with AMI!" << endl;
     }
+
     if (propsDict_.found("splitImplicitExplicit"))
     {
         Info << "will split implicit / explicit force contributions." << endl;
         splitImplicitExplicit_ = true;
-        if(!forceSubM(0).interpolation()) 
+        if(!forceSubM(0).interpolation())
             Info << "WARNING: will only consider fluctuating particle velocity in implicit / explicit force split!" << endl;
     }
     particleCloud_.checkCG(false);
@@ -124,11 +125,11 @@ void DiFeliceDragMS::setForce() const
     const volScalarField& nufField = forceSubM(0).nuField();
     const volScalarField& rhoField = forceSubM(0).rhoField();
 
-    vector position(0,0,0);
+    //vector position(0,0,0);
     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);
     scalar ds(0);
@@ -138,17 +139,17 @@ void DiFeliceDragMS::setForce() const
     scalar Rep(0);
     scalar Cd(0);
 
-	vector UfluidFluct(0,0,0);
+    vector UfluidFluct(0,0,0);
     vector UsFluct(0,0,0);
     vector dragExplicit(0,0,0);
-  	scalar dragCoefficient(0);
+    scalar dragCoefficient(0);
 
     //interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
     //interpolationCellPoint<vector> UInterpolator_(U_);
 
     #include "setupProbeModel.H"
 
-    for(int index = 0;index <  cloudRefMS().numberOfClumps(); index++)
+    for(int index = 0; index < cloudRefMS().numberOfClumps(); ++index)
     {
 
         //if(mask[index][0])  // would have to be transformed from body ID to particle ID
@@ -214,7 +215,7 @@ void DiFeliceDragMS::setForce() const
                     }
                 }
 
-                if(forceSubM(0).verbose() && index >=0 && index <10)
+                if(forceSubM(0).verbose() && index >= 0 && index < 10)
                 {
                     Pout << "index = " << index << endl;
                     Pout << "Us = " << Us << endl;
@@ -246,8 +247,13 @@ void DiFeliceDragMS::setForce() const
                     particleCloud_.probeM().writeProbe(index, sValues, vValues);
                 }
             }
+
             // set force on bodies
-            if(forceSubM(0).switches()[0]) for(int j=0;j<3;j++) cloudRefMS().expForcesCM()[index][j] += drag[j];
+            if (forceSubM(0).switches()[SW_TREAT_FORCE_EXPLICIT])
+            {
+                for(int j=0;j<3;j++)
+                    cloudRefMS().expForcesCM()[index][j] += drag[j];
+            }
             else   //implicit treatment, taking explicit force contribution into account
             {
                 for(int j=0;j<3;j++)
@@ -256,7 +262,9 @@ void DiFeliceDragMS::setForce() const
                     cloudRefMS().expForcesCM()[index][j] += dragExplicit[j];
                 }
             }
-            for(int j=0;j<3;j++) cloudRefMS().DEMForcesCM()[index][j] += drag[j];
+
+            for(int j=0;j<3;j++)
+                cloudRefMS().DEMForcesCM()[index][j] += drag[j];
         //}
     }
 

src/lagrangian/cfdemParticle/subModels/forceModelMS/DiFeliceDragMS/DiFeliceDragMS.H

@@ -71,7 +71,7 @@ private:
 
     bool splitImplicitExplicit_; // use splitting of implicit and explict force contribution
 
-    word UsFieldName_;   
+    word UsFieldName_;
 
     const volVectorField& UsField_;  // the average particle velocity field (for implicit/expliti force split)
 

src/lagrangian/cfdemParticle/subModels/momCoupleModel/explicitCouple/explicitCouple.C

@@ -108,38 +108,15 @@ explicitCouple::~explicitCouple()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 tmp<volVectorField> explicitCouple::expMomSource() const
 {
-    tmp<volVectorField> tsource
-    (
-        new volVectorField
-        (
-            IOobject
-            (
-                "f_explicitCouple",
-                particleCloud_.mesh().time().timeName(),
-                particleCloud_.mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            particleCloud_.mesh(),
-            dimensionedVector
-            (
-                "zero",
-                dimensionSet(1, -2, -2, 0, 0), // N/m3
-                vector::zero
-            ),
-            "zeroGradient"
-        )
-    );
-
     scalar tsf = particleCloud_.dataExchangeM().timeStepFraction();
 
-    if(1-tsf < 1e-4) //tsf==1
+    if (1. - tsf < 1e-4) //tsf==1
     {
         // calc fNext
         forAll(fNext_,cellI)
         {
             fNext_[cellI] = arrayToField(cellI);
-    
+
             // limiter
             for (int i=0;i<3;i++)
             {
@@ -147,21 +124,27 @@ tmp<volVectorField> explicitCouple::expMomSource() const
                 if (magF > fLimit_[i]) fNext_[cellI][i] *= fLimit_[i]/magF;
             }
         }
-        tsource.ref() = fPrev_;
+        return tmp<volVectorField>
-    }else
+        (
-    {
+            new volVectorField("f_explicitCouple", fPrev_)
-        tsource.ref() = (1 - tsf) * fPrev_ + tsf * fNext_;
+        );
+    }
+    else
+    {
+        return tmp<volVectorField>
+        (
+            new volVectorField("f_explicitCouple", (1. - tsf) * fPrev_ + tsf * fNext_)
+        );
     }
-    return tsource;
 }
 
-void Foam::explicitCouple::resetMomSourceField() const
+void explicitCouple::resetMomSourceField() const
 {
     fPrev_.ref() = fNext_.ref();
     fNext_.primitiveFieldRef() = vector::zero;
 }
 
-inline vector Foam::explicitCouple::arrayToField(label cellI) const
+inline vector explicitCouple::arrayToField(label cellI) const
 {
     return particleCloud_.forceM(0).expParticleForces()[cellI] / particleCloud_.mesh().V()[cellI];
 }

src/lagrangian/cfdemParticle/subModels/momCoupleModel/implicitCouple/implicitCouple.C

@@ -122,40 +122,19 @@ implicitCouple::~implicitCouple()
 
 tmp<volScalarField> implicitCouple::impMomSource() const
 {
-    tmp<volScalarField> tsource
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "Ksl_implicitCouple",
-                particleCloud_.mesh().time().timeName(),
-                particleCloud_.mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            particleCloud_.mesh(),
-            dimensionedScalar
-            (
-                "zero",
-                dimensionSet(1, -3, -1, 0, 0), // N/m3 / m/s
-                0
-            )
-        )
-    );
-
     scalar tsf = particleCloud_.dataExchangeM().timeStepFraction();
 
     // calc Ksl
-    scalar Ur;
 
-    if(1-tsf < 1e-4) //tsf==1
+    if (1. - tsf < 1e-4) //tsf==1
     {
+        scalar Ur;
+
         forAll(KslNext_,cellI)
         {
             Ur = mag(U_[cellI] - Us_[cellI]);
 
-            if(Ur > SMALL && alpha_[cellI] < maxAlpha_) //momentum exchange switched off if alpha too big
+            if (Ur > SMALL && alpha_[cellI] < maxAlpha_) //momentum exchange switched off if alpha too big
             {
                 KslNext_[cellI] = mag(particleCloud_.forceM(0).impParticleForces()[cellI])
                             / Ur
@@ -166,16 +145,21 @@ tmp<volScalarField> implicitCouple::impMomSource() const
             // limiter
             if (KslNext_[cellI] > KslLimit_) KslNext_[cellI] = KslLimit_;
         }
-        tsource.ref() = KslPrev_;
+        return tmp<volScalarField>
-    }else
+        (
-    {
+            new volScalarField("Ksl_implicitCouple", KslPrev_)
-        tsource.ref() = (1 - tsf) * KslPrev_ + tsf * KslNext_;
+        );
+    }
+    else
+    {
+        return tmp<volScalarField>
+        (
+            new volScalarField("Ksl_implicitCouple", (1. - tsf) * KslPrev_ + tsf * KslNext_)
+        );
     }
-
-    return tsource;
 }
 
-void Foam::implicitCouple::resetMomSourceField() const
+void implicitCouple::resetMomSourceField() const
 {
     KslPrev_.ref() = KslNext_.ref();
     KslNext_.primitiveFieldRef() = 0;

src/lagrangian/cfdemParticle/subModels/momCoupleModel/noCouple/noCouple.C

@@ -72,7 +72,7 @@ noCouple::~noCouple()
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-void Foam::noCouple::resetMomSourceField() const
+void noCouple::resetMomSourceField() const
 {
     FatalError<<"the solver calls for resetMomSourceField() although you use IB method where this is not needed!\n"
               <<", check your solver! - PANIC -\n";

src/lagrangian/cfdemParticle/subModels/probeModel/particleProbe/particleProbe.C

@@ -87,12 +87,12 @@ particleProbe::particleProbe
     if (propsDict_.found("verbose")) verbose_=true;
     if (propsDict_.found("verboseToFile")) verboseToFile_=true;
 
-    if (propsDict_.found("printEvery")) printEvery_= readScalar(propsDict_.lookup("printEvery"));
+    if (propsDict_.found("printEvery")) printEvery_ = readScalar(propsDict_.lookup("printEvery"));
-    if (propsDict_.found("sampleAll")) sampleAll_=true;
+    if (propsDict_.found("sampleAll")) sampleAll_ = true;
-    if (propsDict_.found("probeDebug")) probeDebug_=true;
+    if (propsDict_.found("probeDebug")) probeDebug_ = true;
-    if (propsDict_.found("includePosition")) includePosition_=true;
+    if (propsDict_.found("includePosition")) includePosition_ = true;
 
-    if (propsDict_.found("writePrecision")) writePrecision_= readScalar(propsDict_.lookup("writePrecision"));
+    if (propsDict_.found("writePrecision")) writePrecision_ = readScalar(propsDict_.lookup("writePrecision"));
 
 }
 
@@ -100,25 +100,36 @@ particleProbe::particleProbe
 
 particleProbe::~particleProbe()
 {
-    clearProbes();
+    forAll(sPtrList_, i)
+        delete sPtrList_[i];
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void particleProbe::setOutputFile() const
+void particleProbe::setOutputFile(const word& logFileName)
 {
-    //set the current item ID
+    if (itemCounter_ > 0 && verboseToFile_)
-    if(currItemId_== itemCounter_)
+    {
-        currItemId_=1;
+        bool foundFile = false;
-    else
+        forAll(itemsToSample_, i)
-        currItemId_+=1;
+        {
-    sPtr = sPtrList_[currItemId_-1]; //set the pointer to the output file from list
+            if (itemsToSample_[i] == logFileName)
-    probeIndex_=currItemId_-1;
+            {
+                probeIndex_ = i;
+                foundFile = true;
+            }
+        }
+
+        if(!foundFile)
+             FatalError << "particleProbe::setOutputFile for logFileName " << logFileName << " : " << "File not found" << abort(FatalError);
+        currItemId_ = probeIndex_ + 1;
+        setCounter();
+    }
 }
 
 
-void particleProbe::initialize(word typeName, word logFileName) const
+void particleProbe::initialize(const word& modelName, const word& logFileName)
 {
     //update the list of items to be sampled
     ++itemCounter_;
@@ -126,7 +137,7 @@ void particleProbe::initialize(word typeName, word logFileName) const
 
     // init environment
     //propsDict_ = particleCloud_.couplingProperties().subDict(typeName + "Props");
-    name_ = typeName;
+    name_ = modelName;
 
     if (verboseToFile_)
     {
@@ -136,8 +147,8 @@ void particleProbe::initialize(word typeName, word logFileName) const
         MPI_Comm_rank(MPI_COMM_WORLD, &rank_);
 
         //open a separate file for each processor
-        char* filecurrent_ = new char[strlen(logFileName.c_str()) + 4]; //reserve 4 chars for processor name
+        char* filecurrent_ = new char[logFileName.length() + 1 + 4 + 1]; //reserve 4 chars for processor name
-        sprintf(filecurrent_,"%s%s%d", logFileName.c_str(), ".", rank_);
+        sprintf(filecurrent_,"%s.%d", logFileName.c_str(), rank_);
 
         Info << "particleProbe for model " << name_ << " will write to file " << filecurrent_ << endl;
 
@@ -171,120 +182,44 @@ void particleProbe::initialize(word typeName, word logFileName) const
         scalarFields_.clear();
         vectorFields_.clear();
     }
-
-    return;
 }
 
+
 void particleProbe::writeHeader() const
 {
+    if (verboseToFile_)
+    {
+        *sPtr << "#processor: " << rank_ << endl;
+        *sPtr << "#index   time   " << "   ";
+        *sPtr << "||  vectorData:  " << "   ";
 
-     if(verboseToFile_ )
+        forAll(vectorFields_, iter)
-     {
+        {
-         *sPtr<<"#processor: " << rank_ << endl;
+            if (!probeDebug_ && iter > 0) break;
-         *sPtr<<"#index   time   "  << "   ";  
+            *sPtr << vectorFields_(iter) << "   ";
+        }
 
-
+        if (probeDebug_)
-        *sPtr<<"||  vectorData:  "  << "   "; 
+        {
-    
+            *sPtr << "||   scalarData:  "  << "   ";
-         forAll(vectorFields_, iter)
-         {
-             if(!probeDebug_ && iter>0) break;
-             *sPtr << vectorFields_(iter) << "   ";
-         }
-
-         if(probeDebug_)
-         {
-            *sPtr<<"||   scalarData:  "  << "   ";
             forAll(scalarFields_, iter)
             {
-                 *sPtr << scalarFields_(iter)  << "   ";
+                *sPtr << scalarFields_(iter)  << "   ";
             }
-         }
-
-         if(includePosition_) *sPtr<<" ||  position" << endl;
-         else *sPtr << endl;
-     }
-
-}
-
-void particleProbe::clearProbes() const
-{
-    for (unsigned int i=0; i<vProbes_.size(); i++)
-        vProbes_[i].clear();
-  
-    for (unsigned int j=0; j<sProbes_.size(); j++)
-        sProbes_[j].clear();
-
-    sProbes_.clear();
-    vProbes_.clear();
-}
-
-void particleProbe::updateProbes(int index, Field<scalar> sValues, Field<vector> vValues) const
-{
-    int vSize_=vProbes_.size();
-    int sSize_=sProbes_.size();
-
-    //check if the particle already has an allocated vector. If not, create it. It should be only called at the beginning.
-    while(index >= vSize_)
-    {
-        std::vector<double*> particleVector_;
-        vProbes_.push_back(particleVector_);
-        vSize_=vProbes_.size();
-    }
-
-    while(index >= sSize_)
-    {
-        std::vector<double> particleScalar_;
-        sProbes_.push_back(particleScalar_);
-        sSize_=sProbes_.size();
-    }
-
-    //register vector probes on the corresponding vector
-    forAll(vValues, iter)
-    {
-        int ProbeSize_=vProbes_[index].size();
-
-        if(probeIndex_<ProbeSize_) //The corresponding probe for this particle already exists, values are overwritten.
-        {
-            vProbes_[index][probeIndex_][0]=vValues[iter][0];
-            vProbes_[index][probeIndex_][1]=vValues[iter][1];
-            vProbes_[index][probeIndex_][2]=vValues[iter][2];
         }
-        else //The corresponding probe for this particle has to be created
-        {
-            double * probe_= new double[3];
 
-            probe_[0]=vValues[iter][0];
+        if (includePosition_) *sPtr << " ||  position" << endl;
-            probe_[1]=vValues[iter][1];
+        else *sPtr << endl;
-            probe_[2]=vValues[iter][2];
-
-            vProbes_[index].push_back(probe_);
-        }
-    }
-
-    //register scalar probes on the corresponding vector
-    forAll(sValues, iter)
-    {
-        int ProbeSize_=sProbes_[index].size();
-
-        if(probeIndex_<ProbeSize_) //The corresponding probe for this particle already exists, values are overwritten.
-        {
-            sProbes_[index][probeIndex_]=sValues[iter];
-        }
-        else //The corresponding probe for this particle has to be created
-        {
-            sProbes_[index].push_back(sValues[iter]);
-        }
     }
 }
 
-void particleProbe::writeProbe(int index, Field<scalar> sValues, Field<vector> vValues) const
+
+void particleProbe::writeProbe(int index, Field<scalar> sValues, Field<vector> vValues)
 {
-    updateProbes(index,sValues,vValues); //update probe vectors
+    if (printNow_ && verboseToFile_ && checkIDForPrint(index))
-
-
-    if(printNow_ && checkIDForPrint(index) &&  verboseToFile_)
     {
+        sPtr = sPtrList_[probeIndex_]; //set the pointer to the output file from list
+
         //index and time
         *sPtr << setprecision(IOstream::defaultPrecision()+7);
         *sPtr << index  << tab << particleCloud_.mesh().time().value()  << "   ";
@@ -298,7 +233,6 @@ void particleProbe::writeProbe(int index, Field<scalar> sValues, Field<vector> v
             *sPtr << vValues[iter][0] << "   ";
             *sPtr << vValues[iter][1] << "   ";
             *sPtr << vValues[iter][2] << "   ";
-
         }
 
         //scalarFields
@@ -324,42 +258,40 @@ void particleProbe::writeProbe(int index, Field<scalar> sValues, Field<vector> v
             *sPtr << endl;
         }
     }
-
-    return;
 }
 
+
 bool particleProbe::checkIDForPrint(int index) const
 {
-      bool sampleThisId_ = false;
+    if(sampleAll_)
-      if(sampleAll_) sampleThisId_ = true;
+    {
-      else
+        return true;
-      {
+    }
-            forAll(particleIDsToSample_, iSample)
+    else
-            {
+    {
-                  if(index==particleIDsToSample_[iSample]) sampleThisId_ = true;
+        forAll(particleIDsToSample_, iSample)
-            }
+        {
-      }
+            if (index == particleIDsToSample_[iSample]) return true;
-      return sampleThisId_;
+        }
+    }
+    return false;
 }
 
-void particleProbe::setCounter() const
+void particleProbe::setCounter()
 {
-
+    //reset or increment counter for printing to file
-    //reset or increment counter for printing to file 
     //Do only if called by first item in the list of items!
-    if(currItemId_==1)
+    if (currItemId_ == 1)
     {
-        printCounter_++;
+        ++printCounter_;
-        if(  printCounter_ >= printEvery_ )
+
+        if (printCounter_ >= printEvery_)
         {
-            printCounter_=0;
+            printCounter_ = 0;
             printNow_ = true;
         }
         else printNow_ = false;
-
     }
-    return;    
-
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/cfdemParticle/subModels/probeModel/particleProbe/particleProbe.H

@@ -63,7 +63,7 @@ private:
 
     dictionary propsDict_;
 
-    mutable word name_;
+    word name_;
 
     cfdemCloud& particleCloud_;
 
@@ -75,7 +75,7 @@ private:
 
     word dirName_;
 
-    mutable int rank_;
+    int rank_;
 
     mutable OFstream* sPtr;
 
@@ -89,25 +89,23 @@ private:
 
     const labelList particleIDsToSample_;
 
-    mutable wordList itemsToSample_;
+    wordList itemsToSample_;
 
-    mutable List<OFstream*> sPtrList_;
+    List<OFstream*> sPtrList_;
 
-    mutable int itemCounter_;
+    int itemCounter_;
 
-    mutable int currItemId_;
+    int currItemId_;
 
-    mutable int printCounter_;
+    int printCounter_;
 
-    mutable bool printNow_;
+    bool printNow_;
 
-    mutable std::vector< std::vector<double> > sProbes_;
+    std::vector<std::string> probeName_;
 
-    mutable std::vector< std::vector<double*> > vProbes_;
+    int probeIndex_;
 
-    mutable  std::vector<std::string> probeName_;
+    void setCounter();
-
-    mutable int probeIndex_;
 
 public:
 
@@ -131,17 +129,11 @@ public:
         ~particleProbe();
 
     // Member Functions
-    void updateProbes(int index, Field<scalar> sValues, Field<vector> vValues) const;
+    void initialize(const word& modelName, const word& logFileName);
-    void initialize(word typeName, word logFileName) const;
+    void setOutputFile(const word& logFileName);
-    void setOutputFile() const;
     void writeHeader() const;
-    void writeProbe(int index, Field<scalar> sValues, Field<vector> vValues) const;
+    void writeProbe(int index, Field<scalar> sValues, Field<vector> vValues);
     bool checkIDForPrint(int) const;
-    void setCounter() const;
-    void clearProbes() const;
-    std::vector< std::vector<double*> >* getVprobe() { return &vProbes_; }
-    std::vector< std::vector<double> >* getSprobe() { return &sProbes_; }
-
 };
 
 

src/lagrangian/cfdemParticle/subModels/probeModel/probeModel/probeModel.H

@@ -133,16 +133,12 @@ public:
 
     // Member Functions
 
-    virtual void initialize(word typeName, word  logFileName) const {}
+    virtual void initialize(const word& modelName, const word& logFileName) {}
-    virtual void setOutputFile() const {}
+    virtual void setOutputFile(const word& logFileName) {}
     virtual void writeHeader() const {}
-    virtual void writeProbe(int index, Field<scalar> sValues, Field<vector> vValues) const {}
+    virtual void writeProbe(int index, Field<scalar> sValues, Field<vector> vValues) {}
     virtual bool checkIDForPrint(int) const { return false; }
-    virtual void setCounter() const {}
     virtual bool active() const { return true; }
-    virtual std::vector< std::vector<double*> >* getVprobe() { return NULL; }
-    virtual std::vector< std::vector<double> >* getSprobe() { return NULL; }
-    virtual void clearProbes() const {}
 
     // Access
 

src/lagrangian/cfdemParticle/subModels/regionModel/regionModel/regionModel.C

@@ -46,13 +46,13 @@ defineRunTimeSelectionTable(regionModel, dictionary);
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::regionModel::reAllocArrays() const
+void regionModel::reAllocArrays() const
 {
     if(particleCloud_.numberOfParticlesChanged())
     {
         // get arrays of new length
-        particleCloud_.dataExchangeM().allocateArray(inRegion_,1,1);
+        particleCloud_.dataExchangeM().allocateArray(inRegion_,1.,1);
-        particleCloud_.dataExchangeM().allocateArray(outRegion_,1,1);
+        particleCloud_.dataExchangeM().allocateArray(outRegion_,1.,1);
     }
 }
 
@@ -83,8 +83,8 @@ regionModel::regionModel
 
 regionModel::~regionModel()
 {
-    free(inRegion_);
+    particleCloud_.dataExchangeM().destroy(inRegion_,1);
-    free(outRegion_);
+    particleCloud_.dataExchangeM().destroy(outRegion_,1);
 }
 
 

src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.C

@@ -6,7 +6,7 @@
                                 Christoph Goniva, christoph.goniva@cfdem.com
                                 Copyright 2009-2012 JKU Linz
                                 Copyright 2012-     DCS Computing GmbH, Linz
-                                Copyright (C) 2013-     Graz University of  
+                                Copyright (C) 2013-     Graz University of
                                                         Technology, IPPT
 -------------------------------------------------------------------------------
 License
@@ -71,10 +71,10 @@ constDiffSmoothing::constDiffSmoothing
     verbose_(false)
 {
 
-    if(propsDict_.found("verbose"))  
+    if(propsDict_.found("verbose"))
         verbose_ = true;
 
-    if(propsDict_.found("smoothingLengthReferenceField"))  
+    if(propsDict_.found("smoothingLengthReferenceField"))
        smoothingLengthReferenceField_.value() = double(readScalar(propsDict_.lookup("smoothingLengthReferenceField")));
 
     checkFields(sSmoothField_);
@@ -94,7 +94,7 @@ bool constDiffSmoothing::doSmoothing() const
 }
 
 
-void Foam::constDiffSmoothing::smoothen(volScalarField& fieldSrc) const
+void constDiffSmoothing::smoothen(volScalarField& fieldSrc) const
 {
     // Create scalar smooth field from virgin scalar smooth field template
     volScalarField sSmoothField = sSmoothField_;
@@ -120,11 +120,11 @@ void Foam::constDiffSmoothing::smoothen(volScalarField& fieldSrc) const
     forAll(sSmoothField,cellI)
     {
         sSmoothField[cellI]=max(lowerLimit_,min(upperLimit_,sSmoothField[cellI]));
-    }  
+    }
 
     // get data from working sSmoothField - will copy only values at new time
     fieldSrc=sSmoothField;
-    fieldSrc.correctBoundaryConditions(); 
+    fieldSrc.correctBoundaryConditions();
 
     if(verbose_)
     {
@@ -135,7 +135,7 @@ void Foam::constDiffSmoothing::smoothen(volScalarField& fieldSrc) const
 
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-void Foam::constDiffSmoothing::smoothen(volVectorField& fieldSrc) const
+void constDiffSmoothing::smoothen(volVectorField& fieldSrc) const
 {
     // Create scalar smooth field from virgin scalar smooth field template
     volVectorField vSmoothField = vSmoothField_;
@@ -159,7 +159,7 @@ void Foam::constDiffSmoothing::smoothen(volVectorField& fieldSrc) const
 
     // get data from working vSmoothField
     fieldSrc=vSmoothField;
-    fieldSrc.correctBoundaryConditions(); 
+    fieldSrc.correctBoundaryConditions();
 
     if(verbose_)
     {
@@ -170,7 +170,7 @@ void Foam::constDiffSmoothing::smoothen(volVectorField& fieldSrc) const
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-void Foam::constDiffSmoothing::smoothenReferenceField(volVectorField& fieldSrc) const
+void constDiffSmoothing::smoothenReferenceField(volVectorField& fieldSrc) const
 {
     // Create scalar smooth field from virgin scalar smooth field template
     volVectorField vSmoothField = vSmoothField_;
@@ -185,7 +185,7 @@ void Foam::constDiffSmoothing::smoothenReferenceField(volVectorField& fieldSrc)
     double sourceStrength = 1e5; //large number to keep reference values constant
 
     dimensionedScalar deltaT = vSmoothField.mesh().time().deltaT();
-    DT_.value() = smoothingLengthReferenceField_.value() 
+    DT_.value() = smoothingLengthReferenceField_.value()
                          * smoothingLengthReferenceField_.value() / deltaT.value();
 
     tmp<volScalarField> NLarge
@@ -194,7 +194,7 @@ void Foam::constDiffSmoothing::smoothenReferenceField(volVectorField& fieldSrc)
         (
             IOobject
             (
-                "xxx",
+                "NLarge",
                 particleCloud_.mesh().time().timeName(),
                 particleCloud_.mesh(),
                 IOobject::NO_READ,
@@ -218,14 +218,14 @@ void Foam::constDiffSmoothing::smoothenReferenceField(volVectorField& fieldSrc)
     (
         fvm::ddt(vSmoothField)
        -fvm::laplacian( DT_, vSmoothField)
-       == 
+       ==
         NLarge() / deltaT * vSmoothField.oldTime()  //add source to keep cell values constant
        -fvm::Sp( NLarge() / deltaT, vSmoothField)   //add sink to keep cell values constant
     );
 
     // get data from working vSmoothField
     fieldSrc=vSmoothField;
-    fieldSrc.correctBoundaryConditions(); 
+    fieldSrc.correctBoundaryConditions();
 
     if(verbose_)
     {

src/lagrangian/cfdemParticle/subModels/smoothingModel/noSmoothing/noSmoothing.C

@@ -70,13 +70,13 @@ noSmoothing::~noSmoothing()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::noSmoothing::smoothen(volScalarField& field) const
+void noSmoothing::smoothen(volScalarField& field) const
 {}
 
-void Foam::noSmoothing::smoothen(volVectorField& field) const
+void noSmoothing::smoothen(volVectorField& field) const
 {}
 
-void Foam::noSmoothing::smoothenReferenceField(volVectorField& field) const
+void noSmoothing::smoothenReferenceField(volVectorField& field) const
 {}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/cfdemParticle/subModels/smoothingModel/smoothingModel/smoothingModel.C

@@ -50,7 +50,7 @@ defineRunTimeSelectionTable(smoothingModel, dictionary);
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
-Foam::smoothingModel::smoothingModel
+smoothingModel::smoothingModel
 (
     const dictionary& dict,
     cfdemCloud& sm
@@ -59,7 +59,7 @@ Foam::smoothingModel::smoothingModel
     dict_(dict),
     particleCloud_(sm),
     vSmoothField_
-    (   
+    (
         IOobject
         (
             "vSmoothField",
@@ -72,7 +72,7 @@ Foam::smoothingModel::smoothingModel
         dimensionedVector("zero", dimensionSet(0,0,0,0,0), vector::zero)
     ),
     sSmoothField_
-    (   
+    (
         IOobject
         (
             "sSmoothField",
@@ -89,7 +89,7 @@ Foam::smoothingModel::smoothingModel
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-Foam::smoothingModel::~smoothingModel()
+smoothingModel::~smoothingModel()
 {}
 
 
@@ -110,7 +110,7 @@ void smoothingModel::checkFields(volVectorField& vSmoothField_) const
 {
     // currently it is detected if field was auto generated or defined
     // improvement would be changing the type here automatically
-    forAll(vSmoothField_.boundaryField(),patchI)      
+    forAll(vSmoothField_.boundaryField(),patchI)
         if(vSmoothField_.boundaryField()[patchI].type()=="calculated")
             FatalError <<"Vector field:"<< vSmoothField_.name() << " must be defined.\n" << abort(FatalError);
 

src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C

@@ -137,7 +137,7 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract
                 if (hashSetLength > maxCellsPerParticle_)
                 {
                     FatalError<< "big particle algo found more cells ("<< hashSetLength
-                              <<") than storage is prepered ("<<maxCellsPerParticle_<<")" << abort(FatalError);
+                              <<") than storage is prepared ("<<maxCellsPerParticle_<<")" << abort(FatalError);
                 }
                 else if (hashSetLength > 0)
                 {

src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.C

@@ -36,8 +36,6 @@ Description
 #include "locateModel.H"
 #include "dataExchangeModel.H"
 
-#include <mpi.h>
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -68,17 +66,15 @@ IBVoidFraction::IBVoidFraction
     propsDict_(dict.subDict(typeName + "Props")),
     alphaMin_(readScalar(propsDict_.lookup("alphaMin"))),
     alphaLimited_(0),
-    scaleUpVol_(readScalar(propsDict_.lookup("scaleUpVol"))),
+    scaleUpVol_(readScalar(propsDict_.lookup("scaleUpVol")))
-    checkPeriodicCells_(false)
 {
-    Info << "\n\n W A R N I N G - do not use in combination with differentialRegion model! \n\n" << endl;
+    Info << "\n\n W A R N I N G - do not use in combination with differentialRegion model!\n\n" << endl;
-    //Info << "\n\n W A R N I N G - this model does not yet work properly! \n\n" << endl;
+    maxCellsPerParticle_ = readLabel(propsDict_.lookup("maxCellsPerParticle"));
-    maxCellsPerParticle_=readLabel(propsDict_.lookup("maxCellsPerParticle"));
 
-    if(scaleUpVol_ < 1){ FatalError<< "scaleUpVol shloud be > 1."<< abort(FatalError); }
+    if (scaleUpVol_ < 1.0)
-    if(alphaMin_ > 1 || alphaMin_ < 0.01){ FatalError<< "alphaMin shloud be > 1 and < 0.01." << abort(FatalError); }
+        FatalError << "scaleUpVol shloud be > 1." << abort(FatalError);
-    
+    if (alphaMin_ > 1.0 || alphaMin_ < 0.01)
-    if(propsDict_.found("checkPeriodicCells")) checkPeriodicCells_=true;
+        FatalError << "alphaMin shloud be > 1 and < 0.01." << abort(FatalError);
 }
 
 
@@ -96,25 +92,26 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
 
     reAllocArrays();
 
-    voidfractionNext_.ref()=1;
+    voidfractionNext_.ref() = 1.0;
 
-    for(int index=0; index< particleCloud_.numberOfParticles(); index++)
+    for (int index=0; index < particleCloud_.numberOfParticles(); index++)
     {
         //if(mask[index][0])
         //{
             //reset
-            for(int subcell=0;subcell<cellsPerParticle_[index][0];subcell++)
+            for (int subcell=0; subcell < cellsPerParticle_[index][0]; subcell++)
             {
-                particleWeights[index][subcell]=0;
+                particleWeights[index][subcell] = 0.0;
-                particleVolumes[index][subcell]=0;
+                particleVolumes[index][subcell] = 0.0;
             }
+
             cellsPerParticle_[index][0]=1;
             particleV[index][0]=0;
 
             //collecting data
             label particleCenterCellID=particleCloud_.cellIDs()[index][0];
             scalar radius =  particleCloud_.radius(index);
-            vector positionCenter=particleCloud_.position(index);
+            vector positionCenter = particleCloud_.position(index);
 
             if (particleCenterCellID >= 0)
             {
@@ -122,92 +119,66 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
 
                 //compute the voidfraction for the cell "particleCentreCellID
                 vector cellCentrePosition = particleCloud_.mesh().C()[particleCenterCellID];
-                scalar centreDist=mag(cellCentrePosition-positionCenter);
+                scalar fc = pointInParticle(index, positionCenter, cellCentrePosition);
+                vector minPeriodicParticlePos = positionCenter;
 
-                vector minPeriodicParticlePos=positionCenter;
+                if (particleCloud_.checkPeriodicCells()) //consider minimal distance to all periodic images of this particle
-                if(checkPeriodicCells_) //consider minimal distance to all periodic images of this particle
                 {
-                    centreDist = minPeriodicDistance(cellCentrePosition, positionCenter, globalBb,
+                    fc = minPeriodicDistance(index,cellCentrePosition, positionCenter, globalBb, minPeriodicParticlePos);
-                                                        minPeriodicParticlePos);
                 }
 
-                if(centreDist + 0.5*sqrt(3.0)*pow(particleCloud_.mesh().V()[particleCenterCellID],0.33333) < radius)
+                scalar centreDist = mag(cellCentrePosition - minPeriodicParticlePos);
+                scalar corona = 0.5 * sqrt(3.0) * pow(particleCloud_.mesh().V()[particleCenterCellID], 1./3.);
+                vector coronaPoint = cellCentrePosition + (cellCentrePosition - minPeriodicParticlePos) * (corona / centreDist);
+
+                if (pointInParticle(index, minPeriodicParticlePos, coronaPoint) < 0.0)
                 {
-                    voidfractionNext_[particleCenterCellID] = 0;
+                    voidfractionNext_[particleCenterCellID] = 0.0;
                 }
                 else
                 {
                     const labelList& vertices = particleCloud_.mesh().cellPoints()[particleCenterCellID];
+                    const double ratio = 0.125;
+
                     forAll(vertices, i)
                     {
                         vector vertexPosition = particleCloud_.mesh().points()[vertices[i]];
-                        scalar centreVertexDist = mag(vertexPosition-positionCenter);
+                        scalar fv = pointInParticle(index, positionCenter, vertexPosition);
-                        if(checkPeriodicCells_) //consider minimal distance to all periodic images of this particle
+
+                        if (particleCloud_.checkPeriodicCells()) //consider minimal distance to all periodic images of this particle
                         {
-                            centreVertexDist = minPeriodicDistance(vertexPosition, positionCenter, globalBb,
+                            fv = minPeriodicDistance(index, vertexPosition, positionCenter, globalBb, minPeriodicParticlePos);
-                                                                        minPeriodicParticlePos);
                         }
 
-                        if(centreDist<radius &&  centreVertexDist<radius)
+                        if (fc < 0.0 && fv < 0.0)
                         {
-                            voidfractionNext_[particleCenterCellID]-=0.125;
+                            voidfractionNext_[particleCenterCellID] -= ratio;
                         }
-                        else if(centreDist<radius && centreVertexDist>radius)
+                        else if (fc < 0.0 && fv > 0.0)
                         {
                             //compute lambda
-                            scalar a = (vertexPosition - cellCentrePosition)
+                            scalar lambda = segmentParticleIntersection(index, minPeriodicParticlePos, cellCentrePosition, vertexPosition);
-                                     & (vertexPosition - cellCentrePosition);
+                            voidfractionNext_[particleCenterCellID] -= ratio * lambda;
-                            scalar b = 2. * (vertexPosition - cellCentrePosition)
-                                     &      (cellCentrePosition-minPeriodicParticlePos);
-                            scalar c = ( (cellCentrePosition-minPeriodicParticlePos)
-                                        &(cellCentrePosition-minPeriodicParticlePos)
-                                       )
-                                        - radius*radius;
-
-                            scalar lambda = 0.;
-
-                            if (b*b-4.*a*c>=0.)  lambda =  (-b+sqrt(b*b-4.*a*c))/(2.*a);
-                            if (lambda > 0 && lambda <=1) voidfractionNext_[particleCenterCellID]-=lambda*0.125;
-                            else
-                            {
-                                lambda =  (-b-sqrt(b*b-4.*a*c))/(2.*a);
-                                if (lambda > 0. && lambda <=1.) voidfractionNext_[particleCenterCellID]-=lambda*0.125;
-                            }
                         }
-                        else if(centreDist>radius && centreVertexDist<radius)
+                        else if (fc > 0.0 && fv < 0.0)
                         {
-                            //compute another lambda too
+                            //compute lambda
-                            scalar a = (vertexPosition - cellCentrePosition)
+                            scalar lambda = segmentParticleIntersection(index, minPeriodicParticlePos, vertexPosition, cellCentrePosition);
-                                     & (vertexPosition - cellCentrePosition);
+                            voidfractionNext_[particleCenterCellID] -= ratio * lambda;
-                            scalar b = 2.* (vertexPosition - cellCentrePosition)
-                                     &     (cellCentrePosition-minPeriodicParticlePos);
-                            scalar c = ( (cellCentrePosition-minPeriodicParticlePos)
-                                        &(cellCentrePosition-minPeriodicParticlePos)
-                                       )
-                                     - radius*radius;
-                            scalar lambda = 0.;
-
-                            if(b*b-4.*a*c>=0.)  lambda =  (-b+sqrt(b*b-4.*a*c))/(2.*a);
-                            if(lambda > 0. && lambda <=1.) voidfractionNext_[particleCenterCellID]-=(1.-lambda)*0.125;
-                            else
-                            {
-                                lambda =  (-b-sqrt(b*b-4.*a*c))/(2.*a);
-                                if (lambda > 0. && lambda <=1.) voidfractionNext_[particleCenterCellID]-=(1.-lambda)*0.125;
-                            }
                         }
                     }
                 } //end particle partially overlapping with cell
 
                 //generating list with cell and subcells
-                buildLabelHashSet(radius, minPeriodicParticlePos, particleCenterCellID, hashSett, true);
+                buildLabelHashSet(index, minPeriodicParticlePos, particleCenterCellID, hashSett, true);
 
                 //Add cells of periodic particle images on same processor
-                if(checkPeriodicCells_)
+                if (particleCloud_.checkPeriodicCells())
                 {
                     int doPeriodicImage[3];
-                    for (int iDir=0; iDir<3; iDir++)
+                    for (int iDir=0; iDir < 3; iDir++)
                     {
-                        doPeriodicImage[iDir]= 0;
+                        doPeriodicImage[iDir] = 0;
                         if ((minPeriodicParticlePos[iDir]+radius) > globalBb.max()[iDir])
                         {
                             doPeriodicImage[iDir] = -1;
@@ -238,10 +209,10 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
                         copyCounter++;
                     }
                     //y-direction
-                    int currCopyCounter=copyCounter;
+                    int currCopyCounter = copyCounter;
                     if (doPeriodicImage[1] != 0)
                     {
-                        for(int yDirCop=0; yDirCop<=currCopyCounter; yDirCop++)
+                        for (int yDirCop=0; yDirCop <= currCopyCounter; yDirCop++)
                         {
                             particlePosList.append( particlePosList[yDirCop]
                                                     + vector(
@@ -254,10 +225,10 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
                         }
                     }
                     //z-direction
-                    currCopyCounter=copyCounter;
+                    currCopyCounter = copyCounter;
                     if (doPeriodicImage[2] != 0)
                     {
-                        for(int zDirCop=0; zDirCop<=currCopyCounter; zDirCop++)
+                        for (int zDirCop=0; zDirCop <= currCopyCounter; zDirCop++)
                         {
                             particlePosList.append( particlePosList[zDirCop]
                                                     + vector(
@@ -273,7 +244,7 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
 
                     //add the nearest cell labels
                     particleLabelList.append(particleCenterCellID);
-                    for(int iPeriodicImage=1; iPeriodicImage<=copyCounter; iPeriodicImage++)
+                    for (int iPeriodicImage=1; iPeriodicImage <= copyCounter; iPeriodicImage++)
                     {
                         label partCellId = particleCloud_.mesh().findNearestCell(particlePosList[iPeriodicImage]);
                         particleLabelList.append(partCellId);
@@ -287,17 +258,17 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
                 label hashSetLength = hashSett.size();
                 if (hashSetLength > maxCellsPerParticle_)
                 {
-                    FatalError<< "big particle algo found more cells ("<< hashSetLength 
+                    FatalError << "big particle algo found more cells (" << hashSetLength
-                              <<") than storage is prepared ("<<maxCellsPerParticle_<<")" << abort(FatalError);
+                               << ") than storage is prepared (" << maxCellsPerParticle_ << ")" << abort(FatalError);
                 }
                 else if (hashSetLength > 0)
                 {
                     cellsPerParticle_[index][0]=hashSetLength;
                     hashSett.erase(particleCenterCellID);
 
-                    for(label i=0;i<hashSetLength-1;i++)
+                    for (label i=0; i < hashSetLength-1; i++)
                     {
-                        label cellI=hashSett.toc()[i];
+                        label cellI = hashSett.toc()[i];
                         particleCloud_.cellIDs()[index][i+1]=cellI; //adding subcell represenation
                     }
                 }//end cells found on this proc
@@ -305,14 +276,16 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
         //}// end if masked
     }// end loop all particles
 
-    for(label index=0; index< particleCloud_.numberOfParticles(); index++)
+    for (label index=0; index < particleCloud_.numberOfParticles(); index++)
     {
-        for(label subcell=0;subcell<cellsPerParticle_[index][0];subcell++)
+        for (label subcell=0; subcell < cellsPerParticle_[index][0]; subcell++)
         {
             label cellID = particleCloud_.cellIDs()[index][subcell];
 
             if(cellID >= 0)
             {
+                if (voidfractionNext_[cellID] < 0.0)
+                    voidfractionNext_[cellID] = 0.0;
                  voidfractions[index][subcell] = voidfractionNext_[cellID];
             }
             else
@@ -325,14 +298,15 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
 
 void IBVoidFraction::buildLabelHashSet
 (
-    const scalar radius,
+    int index,
     const vector position,
     const label cellID,
     labelHashSet& hashSett,
     bool initialInsert //initial insertion of own cell
 )const
 {
-    if(initialInsert)  hashSett.insert(cellID);
+    if(initialInsert)
+        hashSett.insert(cellID);
 
     const labelList& nc = particleCloud_.mesh().cellCells()[cellID];
     forAll(nc,i)
@@ -341,115 +315,99 @@ void IBVoidFraction::buildLabelHashSet
         vector cellCentrePosition = particleCloud_.mesh().C()[neighbor];
         scalar centreDist = mag(cellCentrePosition-position);
 
-        if(!hashSett.found(neighbor) && centreDist + 0.5*sqrt(3.0)*pow(particleCloud_.mesh().V()[neighbor],0.33333) < radius)
+        scalar fc = pointInParticle(index, position, cellCentrePosition);
+        scalar corona = 0.5 * sqrt(3.0) * pow(particleCloud_.mesh().V()[neighbor], 1./3.);
+        vector coronaPoint = cellCentrePosition + (cellCentrePosition - position) * (corona / centreDist);
+
+        if (!hashSett.found(neighbor) && pointInParticle(index, position, coronaPoint) < 0.0)
         {
             voidfractionNext_[neighbor] = 0.;
-            buildLabelHashSet(radius,position,neighbor,hashSett,true);
+            buildLabelHashSet(index,position,neighbor,hashSett,true);
         }
-        else if(!hashSett.found(neighbor) && centreDist < radius + sqrt(3.0)*pow(particleCloud_.mesh().V()[neighbor],0.33333))
+        else if(!hashSett.found(neighbor))
         {
             scalar scale = 1.;
             const labelList& vertexPoints = particleCloud_.mesh().cellPoints()[neighbor];
+            const scalar ratio = 0.125;
 
             forAll(vertexPoints, j)
             {
                 vector vertexPosition = particleCloud_.mesh().points()[vertexPoints[j]];
-                scalar vertexDist = mag(vertexPosition - position);
+                scalar fv = pointInParticle(index, position, vertexPosition);
 
-                if (centreDist < radius)
+                if (fc < 0.0)
                 {
-                    if (vertexDist < radius)
+                    if (fv < 0.0)
                     {
-                        scale -= 0.125;
+                        scale -= ratio;
                     }
                     else
                     {
-                        scalar a = (vertexPosition - cellCentrePosition)&(vertexPosition - cellCentrePosition);
+                        scalar lambda = segmentParticleIntersection(index, position, cellCentrePosition, vertexPosition);
-                        scalar b = 2.* (vertexPosition - cellCentrePosition)&(cellCentrePosition-position);
+                        scale -= lambda * ratio;
-                        scalar c = ((cellCentrePosition-position)&(cellCentrePosition-position))-radius*radius;
-                        scalar lambda = 0.;
-
-                        if(b*b-4.*a*c >= 0.)  lambda =  (-b+sqrt(b*b-4.*a*c))/(2.*a);
-                        if (lambda > 0. && lambda <=1.)
-                        {
-                            scale -=lambda * 0.125;
-                        }
-                        else
-                        {
-                            lambda =  (-b-sqrt(b*b-4.*a*c))/(2.*a);
-                            if (lambda > 0. && lambda <= 1.) scale -=lambda * 0.125;
-                        }
                     }
                 }
-                else if (vertexDist < radius)
+                else if (fv < 0.0)
                 {
-                    scalar a = (vertexPosition - cellCentrePosition)&(vertexPosition - cellCentrePosition);
+                    scalar lambda = segmentParticleIntersection(index, position, vertexPosition, cellCentrePosition);
-                    scalar b = 2.* (vertexPosition - cellCentrePosition)&(cellCentrePosition-position);
+                    scale -= lambda * ratio;
-                    scalar c = ((cellCentrePosition-position)&(cellCentrePosition-position))-radius*radius;
-                    scalar lambda = 0.;
-
-                    if (b*b-4.*a*c >= 0.)  lambda =  (-b+sqrt(b*b-4.*a*c))/(2.*a);
-                    if (lambda > 0. && lambda <= 1.)
-                    {
-                        scale -=(1.-lambda) * 0.125;
-                    }
-                    else
-                    {
-                        lambda =  (-b-sqrt(b*b-4.*a*c))/(2.*a);
-                        if (lambda > 0. && lambda <= 1.) scale -= (1.-lambda) * 0.125;
-                    }
                 }
             }
 
-            if(voidfractionNext_[neighbor] == 1.0)
+            if (scale < 0.0)
+                scale = 0.0;
+
+            if (voidfractionNext_[neighbor] == 1.0)
             {
                 voidfractionNext_[neighbor] = scale;
             }
             else
             {
                 voidfractionNext_[neighbor] -= (1.0-scale);
-                if(voidfractionNext_[neighbor] < 0.) voidfractionNext_[neighbor] = 0.0;
+                if (voidfractionNext_[neighbor] < 0.)
+                    voidfractionNext_[neighbor] = 0.0;
             }
-            if(!(scale == 1.0))  buildLabelHashSet(radius,position,neighbor,hashSett, true);
+
+            if (!(scale == 1.0))
+                buildLabelHashSet(index,position,neighbor,hashSett, true);
         }
     }
 }
 
-//Function to determine minimal distance of point
+double IBVoidFraction::segmentParticleIntersection(int index, vector positionCenter, vector pointInside, vector pointOutside) const
-//to one of the periodic images of a particle
-double IBVoidFraction::minPeriodicDistance(vector    cellCentrePosition, 
-                                           vector    positionCenter, 
-                                           boundBox  globalBb,
-                                           vector&   minPeriodicPos)const
 {
-    double centreDist = 999e32;
+    const scalar radius = particleCloud_.radius(index);
-    vector positionCenterPeriodic;
+    const scalar a = (pointOutside - pointInside) & (pointOutside - pointInside);
+    const scalar b = 2.*(pointOutside - pointInside) & (pointInside - positionCenter);
+    const scalar c = ((pointInside - positionCenter) & (pointInside - positionCenter)) - radius*radius;
+    const scalar D = b*b - 4.0*a*c;
+    const scalar eps = 1e-12;
+    scalar lambda_ = 0.0;
+    scalar lambda = 0.0;
 
-    for (int xDir=-1; xDir<=1; xDir++)
+    if (D >= 0.0)
     {
-        positionCenterPeriodic[0] =  positionCenter[0]
+        scalar sqrtD = sqrt(D);
-                                  + static_cast<double>(xDir)
+        lambda_ = (-b + sqrtD)/(2.0*a);
-                                  * (globalBb.max()[0] - globalBb.min()[0]);
+
-        for (int yDir=-1; yDir<=1; yDir++)
+        if (lambda_ >= -eps && lambda_ <= 1.0+eps)
         {
-            positionCenterPeriodic[1] =  positionCenter[1]
+            lambda = lambda_;
-                                      + static_cast<double>(yDir)
+        }
-                                      * (globalBb.max()[1] - globalBb.min()[1]);
+        else
-            for (int zDir=-1; zDir<=1; zDir++)
+        {
-            {
+            lambda_ = (-b - sqrtD)/(2.0*a);
-                positionCenterPeriodic[2] =  positionCenter[2]
+            if (lambda_ >= -eps && lambda_ <= 1.0+eps)
-                                          + static_cast<double>(zDir)
+                lambda = lambda_;
-                                          * (globalBb.max()[2] - globalBb.min()[2]);
-                if (mag(cellCentrePosition-positionCenterPeriodic) < centreDist)
-                {
-                    centreDist     = mag(cellCentrePosition-positionCenterPeriodic);
-                    minPeriodicPos = positionCenterPeriodic;
-                }
-            }
         }
     }
 
-    return centreDist;
+    if (lambda < 0.0)
+        return 0.0;
+    else if (lambda > 1.0)
+        return 1.0;
+
+    return lambda;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.H

@@ -28,7 +28,7 @@ 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).
 
-    void fraction model for the smooth representation of spheres with 
+    void fraction model for the smooth representation of spheres with
     radius > cell length.
     contribution from Alice Hager
 
@@ -64,11 +64,9 @@ private:
 
         const scalar alphaMin_;          //NP min value of voidFraction
 
-	mutable bool alphaLimited_;
+        mutable bool alphaLimited_;
 
         const scalar scaleUpVol_;        //NP scaling radius, keeping volume of particle
-        
-    	mutable bool checkPeriodicCells_;
 
 public:
 
@@ -95,22 +93,14 @@ public:
 
         void buildLabelHashSet
         (
-            const scalar radius,
+            int index,
             const vector position,
             const label cellID,
             labelHashSet& hashSett,
             bool          initialInsert
         ) const;
-        
+
-        double minPeriodicDistance
+        virtual double segmentParticleIntersection(int index, vector positionCenter, vector pointInside, vector pointOutside) const;
-        (
-            vector    cellCentrePosition, 
-            vector    positionCenter, 
-            boundBox  globalBb,
-            vector&   minPeriodicPos
-        ) const;
-        
-        
 };
 
 

src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C

@@ -132,7 +132,7 @@ void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voi
                 if (hashSetLength > maxCellsPerParticle_)
                 {
                     FatalError<< "big particle algo found more cells ("<< hashSetLength
-                              <<") than storage is prepered ("<<maxCellsPerParticle_<<")" << abort(FatalError);
+                              <<") than storage is prepared ("<<maxCellsPerParticle_<<")" << abort(FatalError);
                 }
                 else if (hashSetLength > 0)
                 {

src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C

@@ -140,6 +140,7 @@ void centreVoidFraction::setvoidFraction(double** const& mask,double**& voidfrac
             }
         //}
     }
+    voidfractionNext_.correctBoundaryConditions();
 
     // bring voidfraction from Eulerian Field to particle array
     for(int index=0; index< particleCloud_.numberOfParticles(); index++)

src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.C

@@ -35,9 +35,7 @@ Description
 #include "addToRunTimeSelectionTable.H"
 #include "locateModel.H"
 #include "dataExchangeModel.H"
-#include <math.h>
 
-//#include <mpi.h>
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -67,29 +65,95 @@ dividedVoidFraction::dividedVoidFraction
     voidFractionModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
     verbose_(false),
+    procBoundaryCorrection_(propsDict_.lookupOrDefault<Switch>("procBoundaryCorrection", false)),
     alphaMin_(readScalar(propsDict_.lookup("alphaMin"))),
     alphaLimited_(0),
     tooMuch_(0.0),
     interpolation_(false),
     cfdemUseOnly_(false)
 {
-    maxCellsPerParticle_ = 29;
+    maxCellsPerParticle_ = numberOfMarkerPoints;
 
-    if(alphaMin_ > 1 || alphaMin_ < 0.01){ Warning << "alphaMin should be < 1 and > 0.01 !!!" << endl; }
+    if (alphaMin_ > 1.0 || alphaMin_ < 0.01)
-    if (propsDict_.found("interpolation")){
+        Warning << "alphaMin should be < 1 and > 0.01 !!!" << endl;
-        interpolation_=true;
+
+    if (propsDict_.found("interpolation"))
+    {
+        interpolation_ = true;
         Warning << "interpolation for dividedVoidFraction does not yet work correctly!" << endl;
-        Info << "Using interpolated voidfraction field - do not use this in combination with interpolation in drag model!"<< endl;
+        Info << "Using interpolated voidfraction field - do not use this in combination with interpolation in drag model!" << endl;
     }
 
     checkWeightNporosity(propsDict_);
 
-    if (propsDict_.found("verbose")) verbose_=true;
+    if (propsDict_.found("verbose"))
+        verbose_ = true;
 
     if (propsDict_.found("cfdemUseOnly"))
     {
         cfdemUseOnly_ = readBool(propsDict_.lookup("cfdemUseOnly"));
     }
+
+    if (procBoundaryCorrection_)
+    {
+        if (!(particleCloud_.locateM().type() == "engineIB"))
+        {
+            FatalError << typeName << ": You are requesting procBoundaryCorrection, this requires the use of engineIB!\n"
+                       << abort(FatalError);
+        }
+    }
+    else
+    {
+        if (particleCloud_.locateM().type() == "engineIB")
+        {
+            FatalError << typeName << ": You are using engineIB, this requires using procBoundaryCorrection=true!\n"
+                       << abort(FatalError);
+        }
+    }
+
+    // generate marker points on unit sphere
+    label m = 0;
+    offsets[m][0] = offsets[m][1] = offsets[m][2] = 0.0;
+    ++m;
+
+    // for 2 different radii
+    scalar r1 = cbrt(1.0/numberOfMarkerPoints);
+    scalar r2 = cbrt(15.0/numberOfMarkerPoints);
+    scalar r[] = { 0.75 * (r2*r2*r2*r2 - r1*r1*r1*r1)/(r2*r2*r2 - r1*r1*r1),
+                   0.75 * (1.0 - r2*r2*r2*r2)/(1.0 - r2*r2*r2) };
+
+    for (label ir = 0; ir < 2; ++ir)
+    {
+        // try 8 subpoints derived from spherical coordinates
+        for (scalar zeta = M_PI_4; zeta < constant::mathematical::twoPi; zeta += constant::mathematical::piByTwo)
+        {
+            for (scalar theta = M_PI_4; theta < constant::mathematical::pi; theta += constant::mathematical::piByTwo)
+            {
+                offsets[m][0] = r[ir] * Foam::sin(theta) * Foam::cos(zeta);
+                offsets[m][1] = r[ir] * Foam::sin(theta) * Foam::sin(zeta);
+                offsets[m][2] = r[ir] * Foam::cos(theta);
+                ++m;
+            }
+        }
+        // try 2 more subpoints for each coordinate direction (6 total)
+        for (label j = -1; j <= 1; j += 2)
+        {
+            offsets[m][0] = r[ir] * j;
+            offsets[m][1] = 0.;
+            offsets[m][2] = 0.;
+            ++m;
+
+            offsets[m][0] = 0.;
+            offsets[m][1] = r[ir] * j;
+            offsets[m][2] = 0.;
+            ++m;
+
+            offsets[m][0] = 0.;
+            offsets[m][1] = 0.;
+            offsets[m][2] = r[ir] * j;
+            ++m;
+        }
+    }
 }
 
 
@@ -103,7 +167,7 @@ dividedVoidFraction::~dividedVoidFraction()
 
 void dividedVoidFraction::setvoidFraction(double** const& mask,double**& voidfractions,double**& particleWeights,double**& particleVolumes, double**& particleV) const
 {
-    if(cfdemUseOnly_)
+    if (cfdemUseOnly_)
         reAllocArrays(particleCloud_.numberOfParticles());
     else
         reAllocArrays();
@@ -115,19 +179,20 @@ void dividedVoidFraction::setvoidFraction(double** const& mask,double**& voidfra
     scalar cellVol(0.);
     scalar scaleVol = weight();
     scalar scaleRadius = pow(porosity(),1./3.);
+    const boundBox& globalBb = particleCloud_.mesh().bounds();
 
-    for(int index=0; index< particleCloud_.numberOfParticles(); index++)
+    for (int index=0; index < particleCloud_.numberOfParticles(); index++)
     {
-        if(!checkParticleType(index)) continue; //skip this particle if not correct type
+        if (!checkParticleType(index)) continue; //skip this particle if not correct type
 
         //if(mask[index][0])
         //{
             // reset
 
-            for(int subcell=0;subcell<cellsPerParticle_[index][0];subcell++)
+            for (int subcell=0; subcell < cellsPerParticle_[index][0]; subcell++)
             {
-                particleWeights[index][subcell]=0;
+                particleWeights[index][subcell] = 0.;
-                particleVolumes[index][subcell]=0;
+                particleVolumes[index][subcell] = 0.;
             }
             particleV[index][0] = 0.;
 
@@ -137,95 +202,71 @@ void dividedVoidFraction::setvoidFraction(double** const& mask,double**& voidfra
             radius = particleCloud_.radius(index);
             volume = Vp(index,radius,scaleVol);
             radius *= scaleRadius;
-            cellVol=0;
+            cellVol = 0;
 
             //--variables for sub-search
-            int nPoints = 29;
+            int nPoints = numberOfMarkerPoints;
-            int nNotFound=0,nUnEqual=0,nTotal=0;
+            int nNotFound = 0, nUnEqual = 0, nTotal = 0;
             vector offset(0.,0.,0.);
             int cellsSet = 0;
 
+            if (procBoundaryCorrection_)
+            {
+                label cellWithCenter(-1);
+                // switch off cellIDs for force calc if steming from parallel search success
+                cellWithCenter = particleCloud_.locateM().findSingleCell(position,cellID);
+                particleCloud_.cellIDs()[index][0] = cellWithCenter;
+            }
+
             if (cellID >= 0)  // particel centre is in domain
             {
                 cellVol = particleCloud_.mesh().V()[cellID];
 
-                // for 2 different radii
+                if (procBoundaryCorrection_)
-                const scalar delta_r = 0.293976*radius;
-                for (scalar r = 0.623926*radius; r < radius; r+=delta_r)
                 {
-                    // try 8 subpoint derived from spherical coordinates
+                    offset = radius * offsets[0];
-                    for (scalar zeta=M_PI_4; zeta<2.*M_PI; zeta+=M_PI_2)
+                    #include "setWeightedSource.H"   // set source terms at position+offset
-                    {
+                }
-                        for (scalar theta=M_PI_4; theta<M_PI; theta+=M_PI_2)
-                        {
-                            offset[0] = r * Foam::sin(theta) * Foam::cos(zeta);
-                            offset[1] = r * Foam::sin(theta) * Foam::sin(zeta);
-                            offset[2] = r * Foam::cos(theta);
-                            #include "setWeightedSource.H"   // set source terms at position+offset
-                        }
-                    }
-                    // try 2 more subpoints for each coordinate direction (6 total)
-                    for (int j=-1; j<=1; j+=2)
-                    {
-                        offset[0] = r * static_cast<double>(j);
-                        offset[1] = 0.;
-                        offset[2] = 0.;
-                        #include "setWeightedSource.H"   //NP set source terms at position+offset
 
-                        offset[0] = 0.;
+                for (label i = 1; i < numberOfMarkerPoints; ++i)
-                        offset[1] = r * static_cast<double>(j);
+                {
-                        offset[2] = 0.;
+                    offset = radius * offsets[i];
-                        #include "setWeightedSource.H"   //NP set source terms at position+offset
+                    #include "setWeightedSource.H"   // set source terms at position+offset
+                }
 
-                        offset[0] = 0.;
+                if (cellsSet > maxCellsPerParticle_ || cellsSet < 0)
-                        offset[1] = 0.;
-                        offset[2] = r * static_cast<double>(j);
-                        #include "setWeightedSource.H"   //NP set source terms at position+offset
-                    }
-                }// end loop radiivoidfractions
-
-                if(cellsSet>29 || cellsSet<0)
                 {
                     Info << "ERROR  cellsSet =" << cellsSet << endl;
                 }
 
-                // set source for particle center; source 1/nPts+weight of all subpoints that have not been found
+                if (!procBoundaryCorrection_)
-                scalar centreWeight = 1./nPoints*(nPoints-cellsSet);
-
-                // update voidfraction for each particle read
-                scalar newAlpha = voidfractionNext_[cellID]- volume*centreWeight/cellVol;
-                if(newAlpha > alphaMin_) voidfractionNext_[cellID] = newAlpha;
-                else
                 {
-                    voidfractionNext_[cellID] = alphaMin_;
+                    // set source for particle center; source 1/nPts+weight of all subpoints that have not been found
-                    tooMuch_ += (alphaMin_-newAlpha) * cellVol;
+                    scalar centreWeight = 1./nPoints*(nPoints-cellsSet);
-                }
 
-                // store cellweight for each particle --- this should be done for subpoints as well!!
+                    // update voidfraction for each particle read
-                particleWeights[index][0] += centreWeight;
+                    scalar newAlpha = voidfractionNext_[cellID]- volume*centreWeight/cellVol;
-
+                    if (newAlpha > alphaMin_)
-                // store particleVolume for each particle
-                particleVolumes[index][0] += volume*centreWeight;
-                particleV[index][0] += volume*centreWeight;
-
-                /*//OUTPUT
-                if (index==0 && verbose_)
-                {
-                    Info << "centre cellID = " << cellID << endl;
-                    Info << "cellsPerParticle_=" << cellsPerParticle_[index][0] << endl;
-
-                    for(int i=0;i<cellsPerParticle_[index][0];i++)
                     {
-                       if(i==0)Info << "cellids, voidfractions, particleWeights, : \n";
+                        voidfractionNext_[cellID] = newAlpha;
-                       Info << particleCloud_.cellIDs()[index][i] << " ," << endl;
+                    }
-                       Info << voidfractions[index][i] << " ," << endl;
+                    else
-                       Info << particleWeights[index][i] << " ," << endl;
+                    {
-                     }
+                        voidfractionNext_[cellID] = alphaMin_;
-                }*/
+                        tooMuch_ += (alphaMin_-newAlpha) * cellVol;
+                    }
 
+                    // store cellweight for each particle --- this should be done for subpoints as well!!
+                    particleWeights[index][0] += centreWeight;
+
+                    // store particleVolume for each particle
+                    particleVolumes[index][0] += volume*centreWeight;
+                    particleV[index][0] += volume*centreWeight;
+                }
             }// end if in cell
         //}// end if in mask
     }// end loop all particles
+    voidfractionNext_.correctBoundaryConditions();
 
     // reset counter of lost volume
     if (verbose_) Pout << "Total particle volume neglected: " << tooMuch_ << endl;
@@ -234,46 +275,17 @@ void dividedVoidFraction::setvoidFraction(double** const& mask,double**& voidfra
     // bring voidfraction from Eulerian Field to particle array
     //interpolationCellPoint<scalar> voidfractionInterpolator_(voidfractionNext_);
     //scalar voidfractionAtPos(0);
-    for(int index=0; index< particleCloud_.numberOfParticles(); index++)
+    for(int index=0; index < particleCloud_.numberOfParticles(); index++)
     {
-        /*if(interpolation_)
         {
-            label cellI = particleCloud_.cellIDs()[index][0];
+            for (int subcell=0; subcell < cellsPerParticle_[index][0]; subcell++)
-            if(cellI >= 0)
-            {
-                position = particleCloud_.position(index);
-                voidfractionAtPos=voidfractionInterpolator_.interpolate(position,cellI);
-            }else{
-                voidfractionAtPos=-1;
-            }
-    
-            for(int subcell=0;subcell<cellsPerParticle_[index][0];subcell++)
             {
                 label cellID = particleCloud_.cellIDs()[index][subcell];
 
-                if(cellID >= 0)
+                if (cellID >= 0)
-                {
-                    if(voidfractionAtPos > 0)
-                        voidfractions[index][subcell] = voidfractionAtPos;
-                    else
-                        voidfractions[index][subcell] = voidfractionNext_[cellID];
-                } 
-                else
-                {
-                    voidfractions[index][subcell] = -1.;
-                }
-            }
-        }
-        else*/
-        {
-            for(int subcell=0;subcell<cellsPerParticle_[index][0];subcell++)
-            {
-                label cellID = particleCloud_.cellIDs()[index][subcell];
-
-                if(cellID >= 0)
                 {
                     voidfractions[index][subcell] = voidfractionNext_[cellID];
-                } 
+                }
                 else
                 {
                     voidfractions[index][subcell] = -1.;

src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.H

@@ -57,13 +57,17 @@ class dividedVoidFraction
 {
 
 private:
+        static const int numberOfMarkerPoints = 29;
+
         dictionary propsDict_;
 
         bool verbose_;
 
+        Switch procBoundaryCorrection_;
+
         const scalar alphaMin_;          // min value of voidFraction
 
-	    mutable bool alphaLimited_;
+        mutable bool alphaLimited_;
 
         mutable scalar tooMuch_;         // particle volume which is lost due to voidFraction limitation
 
@@ -71,6 +75,8 @@ private:
 
         bool cfdemUseOnly_;
 
+        vector offsets[numberOfMarkerPoints];
+
         virtual inline scalar Vp(int index, scalar radius, scalar scaleVol) const
         {
             return 4.188790205*radius*radius*radius*scaleVol; //4/3*pi=4.188790205

src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/setWeightedSource.H

@@ -3,6 +3,18 @@
 {
     // locate subPoint
     vector subPosition = position + offset;
+
+    if (particleCloud_.checkPeriodicCells())
+    {
+        for (int iDir=0; iDir<3; iDir++)
+        {
+            if (subPosition[iDir] > globalBb.max()[iDir])
+                subPosition[iDir] -= globalBb.max()[iDir]-globalBb.min()[iDir];
+            else if (subPosition[iDir] < globalBb.min()[iDir])
+                subPosition[iDir] += globalBb.max()[iDir]-globalBb.min()[iDir];
+      }
+    }
+
     label partCellId = particleCloud_.locateM().findSingleCell(subPosition,cellID);
                
     //NP fprintf(lmp->screen,"cellID=%d, partCellId=%d\n",static_cast<int>(cellID),static_cast<int>(partCellId));

src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C

@@ -0,0 +1,392 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Christoph Goniva, christoph.goniva@cfdem.com
+                                Copyright 2012-     DCS Computing GmbH, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling.
+
+    CFDEMcoupling is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 3 of the License, or (at your
+    option) any later version.
+
+    CFDEMcoupling is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with CFDEMcoupling; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
+    and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "trilinearVoidFraction.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(trilinearVoidFraction, 0);
+
+addToRunTimeSelectionTable
+(
+    voidFractionModel,
+    trilinearVoidFraction,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+trilinearVoidFraction::trilinearVoidFraction
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    voidFractionModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    alphaMin_(readScalar(propsDict_.lookup("alphaMin"))),
+    bb_(particleCloud_.mesh().points(),false),
+    cellVol_(particleCloud_.mesh().V()[0]),
+    cellLength_(pow(cellVol_,1./3.)),
+    nCellXYZ_
+    (
+        round((bb_.max()[0] - bb_.min()[0]) / cellLength_),
+        round((bb_.max()[1] - bb_.min()[1]) / cellLength_),
+        round((bb_.max()[2] - bb_.min()[2]) / cellLength_)
+    )
+{
+    maxCellsPerParticle_ = 8;
+    checkWeightNporosity(propsDict_);
+    if (porosity() != 1.) FatalError << "porosity not used in trilinearVoidFraction" << abort(FatalError);
+
+    Warning << "trilinearVoidFraction model is not yet complete and does not work near boundaries" << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+trilinearVoidFraction::~trilinearVoidFraction()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void trilinearVoidFraction::setvoidFraction(double** const& mask,double**& voidfractions,double**& particleWeights,double**& particleVolumes,double**& particleV) const
+{
+    reAllocArrays();
+
+    scalar radius(-1.);
+    scalar volume(0.);
+    const scalar scaleVol = weight();
+    const scalar fourPiByThree = 4.0*constant::mathematical::pi/3.0;
+
+    vector partPos(0.,0.,0.);
+    vector pt(0.,0.,0.);
+    vector posShift(0.,0.,0.);
+    vector offsetCell(0.,0.,0.);
+    vector offsetOrigin(0.,0.,0.);
+
+    label i000(0);
+    label i100(0);
+    label i110(0);
+    label i101(0);
+    label i111(0);
+    label i010(0);
+    label i011(0);
+    label i001(0);
+
+    scalar C000(0.);
+    scalar C100(0.);
+    scalar C110(0.);
+    scalar C101(0.);
+    scalar C111(0.);
+    scalar C010(0.);
+    scalar C011(0.);
+    scalar C001(0.);
+
+    scalar x(0.);
+    scalar y(0.);
+    scalar z(0.);
+
+    scalar a(0.);
+    scalar b(0.);
+    scalar c(0.);
+
+    for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
+    {
+        // reset
+        cellsPerParticle_[index][0] = 8;
+        //TODO do we need to set particleVolumes, particleV?
+        // ===
+
+        label cellI = particleCloud_.cellIDs()[index][0];
+
+        if (cellI >= 0)  // particel centre is in domain
+        {
+            radius = particleCloud_.radius(index);
+            volume = fourPiByThree * radius * radius * radius * scaleVol;
+
+            // store volume for each particle
+            particleVolumes[index][0] = volume;
+            particleV[index][0] = volume;
+
+            // find a,b,c
+            partPos = particleCloud_.position(index);
+            offsetCell = partPos - particleCloud_.mesh().C()[cellI];
+            a = offsetCell[0];
+            b = offsetCell[1];
+            c = offsetCell[2];
+
+            // find "origin" index for mapping
+            if ( a > 0.)
+            {
+                if (b > 0.)
+                {
+                    if (c > 0.) //FNE
+                        i000 = cellI;
+                    else       //BNE
+                        i000 = cellI - nCellXYZ_[0] * nCellXYZ_[1];
+                }
+                else
+                {
+                    if (c > 0.) //FSE
+                        i000 = cellI - nCellXYZ_[0];
+                    else       //BSE
+                        i000 = cellI - nCellXYZ_[0] - nCellXYZ_[0] * nCellXYZ_[1];
+                }
+            }
+            else
+            {
+                if (b > 0.)
+                {
+                    if (c > 0.) //FNW
+                        i000 = cellI - 1;
+                    else       //BNW
+                        i000 = cellI - 1 - nCellXYZ_[0] * nCellXYZ_[1];
+                }
+                else
+                {
+                    if (c > 0.) //FSW
+                        i000 = cellI - 1 - nCellXYZ_[0];
+                    else       //BSW
+                        i000 = cellI - 1 - nCellXYZ_[0] - nCellXYZ_[0] * nCellXYZ_[1];
+                }
+            }
+
+            // check boundaries
+            // TODO different handling for periodic and processor boundaries
+            pt = particleCloud_.mesh().C()[cellI];
+            posShift = vector::zero;
+            if (a > 0.)
+            {
+                pt += vector(cellLength_,0.,0.);
+                if (pt[0] > bb_.max()[0])
+                {
+                    --i000;
+                    posShift[0] = -a;
+                }
+            }
+            else
+            {
+                pt -= vector(cellLength_,0.,0.);
+                if (pt[0] < bb_.min()[0])
+                {
+                    ++i000;
+                    posShift[0] = -a;
+                }
+            }
+            if (b > 0.)
+            {
+                pt += vector(0.,cellLength_,0.);
+                if (pt[1] > bb_.max()[1])
+                {
+                    i000 -= nCellXYZ_[0];
+                    posShift[1] = -b;
+                }
+            }
+            else
+            {
+                pt -= vector(0.,cellLength_,0.);
+                if (pt[1] < bb_.min()[1])
+                {
+                    i000 += nCellXYZ_[0];
+                    posShift[1] = -b;
+                }
+            }
+            if (c > 0.)
+            {
+                pt += vector(0.,0.,cellLength_);
+                if (pt[2] > bb_.max()[2])
+                {
+                    i000 -= nCellXYZ_[0] * nCellXYZ_[1];
+                    posShift[2] = -c;
+                }
+            }
+            else
+            {
+                pt -= vector(0.,0.,cellLength_);
+                if (pt[2] < bb_.min()[2])
+                {
+                    i000 += nCellXYZ_[0] * nCellXYZ_[1];
+                    posShift[2] = -c;
+                }
+            }
+
+            // define other 7 indices
+            i100 = i000 + 1;
+            i110 = i100 + nCellXYZ_[0];
+            i010 = i000 + nCellXYZ_[0];
+            i001 = i000 + nCellXYZ_[0] * nCellXYZ_[1];
+            i101 = i100 + nCellXYZ_[0] * nCellXYZ_[1];
+            i111 = i110 + nCellXYZ_[0] * nCellXYZ_[1];
+            i011 = i010 + nCellXYZ_[0] * nCellXYZ_[1];
+
+            // find x,y,z
+            // TODO changes needed here when generalized for quader cells
+            offsetOrigin = particleCloud_.mesh().C()[i000] - (partPos + posShift);
+            x = mag(offsetOrigin[0]) / cellLength_;
+            y = mag(offsetOrigin[1]) / cellLength_;
+            z = mag(offsetOrigin[2]) / cellLength_;
+
+            // calculate the mapping coeffs
+            C000 = (1 - x) * (1 - y) * (1 - z);
+            C100 = x       * (1 - y) * (1 - z);
+            C110 = x       * y       * (1 - z);
+            C010 = (1 - x) * y       * (1 - z);
+            C001 = (1 - x) * (1 - y) * z;
+            C101 = x       * (1 - y) * z;
+            C111 = x       * y       * z;
+            C011 = (1 - x) * y       * z;
+
+            // set weights
+            particleWeights[index][0] = C000;
+            particleWeights[index][1] = C100;
+            particleWeights[index][2] = C110;
+            particleWeights[index][3] = C010;
+            particleWeights[index][4] = C001;
+            particleWeights[index][5] = C101;
+            particleWeights[index][6] = C111;
+            particleWeights[index][7] = C011;
+
+            // set cellIDs
+            particleCloud_.cellIDs()[index][0] = i000;
+            particleCloud_.cellIDs()[index][1] = i100;
+            particleCloud_.cellIDs()[index][2] = i110;
+            particleCloud_.cellIDs()[index][3] = i010;
+            particleCloud_.cellIDs()[index][4] = i001;
+            particleCloud_.cellIDs()[index][5] = i101;
+            particleCloud_.cellIDs()[index][6] = i111;
+            particleCloud_.cellIDs()[index][7] = i011;
+
+            //distribute volume
+            // TODO use different cell volume when generalized for quader cells
+            voidfractionNext_[i000] -= volume*C000 / cellVol_;
+            voidfractionNext_[i100] -= volume*C100 / cellVol_;
+            voidfractionNext_[i010] -= volume*C010 / cellVol_;
+            voidfractionNext_[i001] -= volume*C001 / cellVol_;
+            voidfractionNext_[i101] -= volume*C101 / cellVol_;
+            voidfractionNext_[i011] -= volume*C011 / cellVol_;
+            voidfractionNext_[i110] -= volume*C110 / cellVol_;
+            voidfractionNext_[i111] -= volume*C111 / cellVol_;
+
+            // debugging
+            /*Pout << "cellI=" << cellI << endl;
+            Pout << "a=" << a << endl;
+            Pout << "b=" << b << endl;
+            Pout << "c=" << c << endl;
+            Pout << "x=" << x << endl;
+            Pout << "y=" << y << endl;
+            Pout << "z=" << z << endl;
+            Pout << "i000=" << i000 << endl;
+            Pout << "i100=" << i100 << endl;
+            Pout << "i010=" << i010 << endl;
+            Pout << "i001=" << i001 << endl;
+            Pout << "i101=" << i101 << endl;
+            Pout << "i011=" << i011 << endl;
+            Pout << "i110=" << i110 << endl;
+            Pout << "i111=" << i111 << endl;
+
+            Pout << "C000=" << C000 << endl;
+            Pout << "C100=" << C100 << endl;
+            Pout << "C010=" << C010 << endl;
+            Pout << "C001=" << C001 << endl;
+            Pout << "C101=" << C101 << endl;
+            Pout << "C011=" << C011 << endl;
+            Pout << "C110=" << C110 << endl;
+            Pout << "C111=" << C111 << endl;
+            Pout << "sum(Cijk)=" << C000+C100+C010+C001+C101+C011+C110+C111 << endl;*/
+
+            /*voidfractionNext_[i000]=0.999;
+            voidfractionNext_[i100]=0.999;
+            voidfractionNext_[i010]=0.999;
+            voidfractionNext_[i001]=0.999;
+            voidfractionNext_[i101]=0.999;
+            voidfractionNext_[i011]=0.999;
+            voidfractionNext_[i110]=0.999;
+            voidfractionNext_[i111]=0.999;*/
+
+            // limit volumefraction
+            // TODO implement limiter for all 8 indices
+            /*if(voidfractionNext_[cellI] < alphaMin_ )
+            {
+                voidfractionNext_[cellI] = alphaMin_;
+                alphaLimited = true;
+            }
+            if(index==0 && alphaLimited) Info<<"alpha limited to" <<alphaMin_<<endl;*/
+
+            // store voidFraction for each particle
+            voidfractions[index][0] = voidfractionNext_[cellI];
+
+            // store cellweight for each particle  - this should not live here
+            particleWeights[index][0] = 1.;
+        }
+    }
+    voidfractionNext_.correctBoundaryConditions();
+
+    // bring voidfraction from Eulerian Field to particle array
+    for (int index = 0; index < particleCloud_.numberOfParticles(); ++index)
+    {
+        label cellID = particleCloud_.cellIDs()[index][0];
+
+        if (cellID >= 0)
+        {
+            voidfractions[index][0] = voidfractionNext_[i000];
+            voidfractions[index][1] = voidfractionNext_[i100];
+            voidfractions[index][2] = voidfractionNext_[i110];
+            voidfractions[index][3] = voidfractionNext_[i010];
+            voidfractions[index][4] = voidfractionNext_[i001];
+            voidfractions[index][5] = voidfractionNext_[i101];
+            voidfractions[index][6] = voidfractionNext_[i111];
+            voidfractions[index][7] = voidfractionNext_[i011];
+        }
+        else
+        {
+            for (int i = 0; i < 8; ++i)
+                voidfractions[index][i] = -1.;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.H

@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Christoph Goniva, christoph.goniva@cfdem.com
+                                Copyright 2012-     DCS Computing GmbH, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling.
+
+    CFDEMcoupling is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 3 of the License, or (at your
+    option) any later version.
+
+    CFDEMcoupling is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with CFDEMcoupling; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
+    and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
+
+Class
+    trilinearVoidFraction
+
+SourceFiles
+    trilinearVoidFraction.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef trilinearVoidFraction_H
+#define trilinearVoidFraction_H
+
+#include "voidFractionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class trilinearVoidFraction Declaration
+\*---------------------------------------------------------------------------*/
+
+class trilinearVoidFraction
+:
+    public voidFractionModel
+{
+
+private:
+        dictionary propsDict_;
+
+        const scalar alphaMin_;
+
+        boundBox bb_;
+
+        scalar cellVol_;
+
+        scalar cellLength_;
+
+        vector nCellXYZ_;
+
+public:
+
+    //- Runtime type information
+    TypeName("trilinear");
+
+
+    // Constructors
+
+        //- Construct from components
+        trilinearVoidFraction
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+    // Destructor
+
+        ~trilinearVoidFraction();
+
+
+    // Member Functions
+        void setvoidFraction(double** const& ,double**&, double**&, double**&, double**&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/newVoidFractionModel.C

@@ -32,7 +32,6 @@ Description
 #include "error.H"
 
 #include "voidFractionModel.H"
-#include "centreVoidFraction.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C

@@ -100,70 +100,38 @@ voidFractionModel::~voidFractionModel()
 }
 
 // * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
-tmp<volScalarField> Foam::voidFractionModel::voidFractionInterp() const
+tmp<volScalarField> voidFractionModel::voidFractionInterp() const
 {
-    tmp<volScalarField> tsource
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "alpha_voidFractionModel",
-                particleCloud_.mesh().time().timeName(),
-                particleCloud_.mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            particleCloud_.mesh(),
-            dimensionedScalar
-            (
-                "zero",
-                dimensionSet(0, 0, 0, 0, 0),
-                0
-            )
-        )
-    );
-
     scalar tsf = particleCloud_.dataExchangeM().timeStepFraction();
-    if(1-tsf < 1e-4 && particleCloud_.dataExchangeM().couplingStep() > 1) //tsf==1
+
+    if (1. - tsf < 1e-4 && particleCloud_.dataExchangeM().couplingStep() > 1) //tsf==1
     {
-        tsource.ref() = voidfractionPrev_;
+        return tmp<volScalarField>
+        (
+            new volScalarField("alpha_voidFractionModel", voidfractionPrev_)
+        );
     }
     else
     {
-        tsource.ref() = (1 - tsf) * voidfractionPrev_ + tsf * voidfractionNext_;
+        return tmp<volScalarField>
+        (
+            new volScalarField("alpha_voidFractionModel", (1. - tsf) * voidfractionPrev_ + tsf * voidfractionNext_)
+        );
     }
-    return tsource;
 }
 
-void Foam::voidFractionModel::resetVoidFractions() const
+void voidFractionModel::resetVoidFractions() const
 {
     voidfractionPrev_.ref() = voidfractionNext_.ref();
-    voidfractionNext_.ref() = 1;
+    voidfractionNext_.ref() = 1.;
 }
 
-/*void Foam::voidFractionModel::undoVoidFractions(double**const& mask) const
-{
-    voidfractionPrev_.internalField() = voidfractionNext_.internalField();
-
-    for(int index=0; index< particleCloud_.numberOfParticles(); index++)
-    {
-        if(mask[index][0])
-        {
-            // undo voidfraction cause by particle
-            label cellI = particleCloud_.cellIDs()[index][0];
-            scalar cellVolume=voidfractionNext_.mesh().V()[cellI];
-            voidfractionNext_[cellI] += particleCloud_.particleVolumes()[index][0]/cellVolume;
-        }
-    }
-}*/
-
 int** const& voidFractionModel::cellsPerParticle() const
 {
     return cellsPerParticle_;
 }
 
-int Foam::voidFractionModel::maxCellsPerParticle() const
+int voidFractionModel::maxCellsPerParticle() const
 {
     return maxCellsPerParticle_;
 }
@@ -186,6 +154,60 @@ void voidFractionModel::reAllocArrays(int nP) const
     }
 }
 
+scalar voidFractionModel::pointInParticle(int index, const vector& positionCenter, const vector& point, double scale) const
+{
+    const scalar radius = particleCloud_.radius(index);
+
+    if(radius > SMALL)
+    {
+        scalar pointDistSq = magSqr(point - positionCenter);
+        return pointDistSq / (scale*scale * radius*radius) - 1.0;
+    }
+    else
+    {
+        return 0.;
+    }
+}
+
+//Function to determine minimal distance of point
+//to one of the periodic images of a particle
+scalar voidFractionModel::minPeriodicDistance(int index,
+                                           const vector& cellCentrePosition,
+                                           const vector& positionCenter,
+                                           const boundBox& globalBb,
+                                           vector& minPeriodicPos) const
+{
+    scalar f = VGREAT;
+    vector positionCenterPeriodic;
+
+    for(label xDir=-1; xDir<=1; ++xDir)
+    {
+        positionCenterPeriodic[0] =  positionCenter[0]
+                                  + static_cast<scalar>(xDir)
+                                  * (globalBb.max()[0]-globalBb.min()[0]);
+        for(label yDir=-1; yDir<=1; ++yDir)
+        {
+            positionCenterPeriodic[1] =  positionCenter[1]
+                                      + static_cast<scalar>(yDir)
+                                      * (globalBb.max()[1]-globalBb.min()[1]);
+            for(label zDir=-1; zDir<=1; ++zDir)
+            {
+                positionCenterPeriodic[2] =  positionCenter[2]
+                                          + static_cast<scalar>(zDir)
+                                          * (globalBb.max()[2]-globalBb.min()[2]);
+
+                if(pointInParticle(index, positionCenterPeriodic, cellCentrePosition) < f)
+                {
+                    f = pointInParticle(index, positionCenterPeriodic, cellCentrePosition);
+                    minPeriodicPos = positionCenterPeriodic;
+                }
+            }
+        }
+    }
+
+    return f;
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.H

@@ -155,6 +155,17 @@ public:
         virtual void setParticleType(label type) const {}
 
         virtual bool checkParticleType(label) const { return true; } //consider all particles by default
+
+        virtual scalar minPeriodicDistance
+        (
+            int index,
+            const vector& cellCentrePosition,
+            const vector& positionCenter,
+            const boundBox& globalBb,
+            vector& minPeriodicPos
+        ) const;
+
+        virtual scalar pointInParticle(int index, const vector& positionCenter, const vector& point, double scale=1.0) const;
 };