release on 2012-11-28_19-59-46

README

@@ -0,0 +1,82 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Christoph Goniva, christoph.goniva@cfdem.com
+                                Copyright 2009-2012 JKU Linz
+                                Copyright 2012-     DCS Computing GmbH, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling.
+
+    CFDEMcoupling is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 3 of the License, or (at your
+    option) any later version.
+
+    CFDEMcoupling is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with CFDEMcoupling; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
+    and OpenFOAM. Note: this code is not part of OpenFOAM (see DISCLAIMER).
+\*---------------------------------------------------------------------------*/
+
+
+CFDEM coupling provides an open source parallel coupled CFD-DEM framework 
+combining the strengths of LIGGGHTS DEM code and the Open Source 
+CFD package OpenFOAM(R)(*). The CFDEMcoupling toolbox allows to expand 
+standard CFD solvers of OpenFOAM(R)(*) to include a coupling to the DEM 
+code LIGGGHTS. In this toolbox the particle representation within the 
+CFD solver is organized by "cloud" classes. Key functionalities are organised 
+in sub-models (e.g. force models, data exchange models, etc.) which can easily 
+be selected and combined by dictionary settings.
+
+The coupled solvers run fully parallel on distributed-memory clusters. 
+
+Features are:
+
+- its modular approach allows users to easily implement new models
+- its MPI parallelization enables to use it for large scale problems
+- the "forum"_lws on CFD-DEM gives the possibility to exchange with other 
+  users / developers
+- the use of GIT allows to easily update to the latest version
+- basic documentation is provided
+
+The file structure:
+
+- "src" directory including the source files of the coupling toolbox and models
+- "applications" directory including the solver files for coupled CFD-DEM simulations
+- "doc" directory including the documentation of CFDEMcoupling
+- "tutorials" directory including basic tutorial cases showing the functionality
+
+
+
+Details on installation are given on the "www.cfdem.com"
+
+The functionality of this CFD-DEM framwork is described via "tutorial cases" showing 
+how to use different solvers and models.
+
+CFDEMcoupling stands for Computational Fluid Dynamics (CFD) -
+Discrete Element Method (DEM) coupling.
+
+CFDEMcoupling is an open-source code, distributed freely under the terms of the 
+GNU Public License (GPL).
+
+Core development of CFDEMcoupling is done by 
+Christoph Goniva and Christoph Kloss, both at DCS Computing GmbH, 2012
+
+
+\*---------------------------------------------------------------------------*/
+(*) "OpenFOAM(R)"_of is a registered trade mark of Silicon Graphics 
+International Corp. This offering is not affiliated, approved or endorsed by 
+Silicon Graphics International Corp., the producer of the OpenFOAM(R) software 
+and owner of the OpenFOAM(R) trademark.
+\*---------------------------------------------------------------------------*/

doc/locateModel_engineSearch.html

@@ -28,7 +28,7 @@ engineProps
 <PRE>locateModel engine;
 engineProps
 {
-    treeSearch false;
+    treeSearch true;
 } 
 </PRE>
 <P><B>Description:</B>
@@ -36,9 +36,9 @@ engineProps
 <P>The locateModel "engine" locates the CFD cell and cellID corresponding to a given position.
 The engineSearch locate Model can be used with different settings to use different algorithms:
 </P>
-<UL><LI>treeSearch false;     will execute some geometric (linear) search using the last known cellID (recommended) 
+<UL><LI>treeSearch false;     will execute some geometric (linear) search using the last known cellID 
 
-<LI>treeSearch true;      will use a recursive tree structure to find the cell. 
+<LI>treeSearch true;      will use a recursive tree structure to find the cell (recommended). 
 
 
 </UL>

doc/locateModel_engineSearch.txt

@@ -25,7 +25,7 @@ engineProps
 locateModel engine;
 engineProps
 \{
-    treeSearch false;
+    treeSearch true;
 \} :pre
 
 [Description:]
@@ -33,8 +33,8 @@ engineProps
 The locateModel "engine" locates the CFD cell and cellID corresponding to a given position.
 The engineSearch locate Model can be used with different settings to use different algorithms:
 
-treeSearch false;     will execute some geometric (linear) search using the last known cellID (recommended) :ulb,l
+treeSearch false;     will execute some geometric (linear) search using the last known cellID :ulb,l
-treeSearch true;      will use a recursive tree structure to find the cell. :l
+treeSearch true;      will use a recursive tree structure to find the cell (recommended). :l
 :ule
 
 [Restrictions:] none.

doc/locateModel_turboEngineSearch.html

@@ -28,7 +28,7 @@ turboEngineProps
 <PRE>locateModel turboEngine;
 turboEngineProps
 {
-    treeSearch false;
+    treeSearch true;
 } 
 </PRE>
 <P><B>Description:</B>
@@ -37,10 +37,9 @@ turboEngineProps
 </P>
 <P>The turboEngineSearch locate Model can be used with different settings to use different algorithms:
 </P>
-<UL><LI>faceDecomp false; treeSearch false;     will execute some geometric (linear) search using the last known cellID (recommended) 
+<UL><LI>faceDecomp false; treeSearch false;     will execute some geometric (linear) search using the last known cellID 
-
-<LI>faceDecomp false; treeSearch true;      will use a recursive tree structure to find the cell. 
 
+faceDecomp false; treeSearch true;      will use a recursive tree structure to find the cell. (recommended):l
 
 </UL>
 <P><B>Restrictions:</B> none.

doc/locateModel_turboEngineSearch.txt

@@ -25,7 +25,7 @@ turboEngineProps
 locateModel turboEngine;
 turboEngineProps
 \{
-    treeSearch false;
+    treeSearch true;
 \} :pre
 
 [Description:]
@@ -34,8 +34,8 @@ The locateModel "turboEngine" locates the CFD cell and cellID corresponding to a
 
 The turboEngineSearch locate Model can be used with different settings to use different algorithms:
 
-faceDecomp false; treeSearch false;     will execute some geometric (linear) search using the last known cellID (recommended) :ulb,l
+faceDecomp false; treeSearch false;     will execute some geometric (linear) search using the last known cellID :ulb,l
-faceDecomp false; treeSearch true;      will use a recursive tree structure to find the cell. :l
+faceDecomp false; treeSearch true;      will use a recursive tree structure to find the cell. (recommended):l
 :ule
 
 [Restrictions:] none.

src/lagrangian/cfdemParticle/Make/files

@@ -20,18 +20,33 @@ $(forceModels)/forceModel/forceModel.C
 $(forceModels)/forceModel/newForceModel.C
 $(forceModels)/noDrag/noDrag.C
 $(forceModels)/DiFeliceDrag/DiFeliceDrag.C
+$(forceModels)/DiFeliceDragNLift/DiFeliceDragNLift.C
 $(forceModels)/GidaspowDrag/GidaspowDrag.C
 $(forceModels)/SchillerNaumannDrag/SchillerNaumannDrag.C
 $(forceModels)/Archimedes/Archimedes.C
 $(forceModels)/ArchimedesIB/ArchimedesIB.C
 $(forceModels)/interface/interface.C
 $(forceModels)/ShirgaonkarIB/ShirgaonkarIB.C
+$(forceModels)/fieldTimeAverage/fieldTimeAverage.C
+$(forceModels)/fieldBound/fieldBound.C
+$(forceModels)/volWeightedAverage/volWeightedAverage.C
+$(forceModels)/totalMomentumExchange/totalMomentumExchange.C
 $(forceModels)/KochHillDrag/KochHillDrag.C
+$(forceModels)/BeetstraDrag/multiphaseFlowBasic/multiphaseFlowBasic.C
+$(forceModels)/BeetstraDrag/BeetstraDrag.C
+$(forceModels)/LaEuScalarLiquid/LaEuScalarLiquid.C
 $(forceModels)/LaEuScalarTemp/LaEuScalarTemp.C
+$(forceModels)/LaEuScalarDust/LaEuScalarDust.C
 $(forceModels)/virtualMassForce/virtualMassForce.C
 $(forceModels)/gradPForce/gradPForce.C
+$(forceModels)/gradULiftForce/gradULiftForce.C
 $(forceModels)/viscForce/viscForce.C
 $(forceModels)/MeiLift/MeiLift.C
+$(forceModels)/KochHillDragNLift/KochHillDragNLift.C
+$(forceModels)/solidsPressureForce/solidsPressureForce.C
+$(forceModels)/periodicPressure/periodicPressure.C
+$(forceModels)/periodicPressureControl/periodicPressureControl.C
+$(forceModels)/averageSlipVel/averageSlipVel.C
 
 $(forceModelsMS)/forceModelMS/forceModelMS.C
 $(forceModelsMS)/forceModelMS/newForceModelMS.C
@@ -42,6 +57,7 @@ $(IOModels)/IOModel/newIOModel.C
 $(IOModels)/noIO/noIO.C
 $(IOModels)/basicIO/basicIO.C
 $(IOModels)/trackIO/trackIO.C
+$(IOModels)/sophIO/sophIO.C
 
 $(voidFractionModels)/voidFractionModel/voidFractionModel.C
 $(voidFractionModels)/voidFractionModel/newVoidFractionModel.C
@@ -60,20 +76,22 @@ $(locateModels)/turboEngineSearch/turboEngineSearch.C
 $(locateModels)/turboEngineSearchM2M/turboEngineSearchM2M.C
 $(locateModels)/engineSearchIB/engineSearchIB.C
 
-
 $(meshMotionModels)/meshMotionModel/meshMotionModel.C
 $(meshMotionModels)/meshMotionModel/newMeshMotionModel.C
 $(meshMotionModels)/noMeshMotion/noMeshMotion.C
+$(meshMotionModels)/DEMdrivenMeshMotion/DEMdrivenMeshMotion.C
 
 $(momCoupleModels)/momCoupleModel/momCoupleModel.C
 $(momCoupleModels)/momCoupleModel/newMomCoupleModel.C
 $(momCoupleModels)/explicitCouple/explicitCouple.C
+$(momCoupleModels)/explicitCoupleSource/explicitCoupleSource.C
 $(momCoupleModels)/implicitCouple/implicitCouple.C
 $(momCoupleModels)/noCouple/noCouple.C
 
 $(regionModels)/regionModel/regionModel.C
 $(regionModels)/regionModel/newRegionModel.C
 $(regionModels)/allRegion/allRegion.C
+$(regionModels)/differentialRegion/differentialRegion.C
 
 $(dataExchangeModels)/dataExchangeModel/dataExchangeModel.C
 $(dataExchangeModels)/dataExchangeModel/newDataExchangeModel.C

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

@@ -148,7 +148,9 @@ twoWayM2M::~twoWayM2M()
     delete[] lost_pos_;
     delete[] cellID_foam_;
     delete[] pos_foam_;
-    //delete[] lmp2foam_;
+    delete lmp2foam_;
+    delete lmp2foam_vec_;
+    delete foam2lmp_vec_;
     //delete lmp;
 }
 
@@ -247,7 +249,7 @@ void Foam::twoWayM2M::allocateArray
 ) const
 {
     //if(length==-1) then LIGGGHTS uses own length data
-    allocate_external_double(array, width,length,initVal,lmp);
+    allocate_external_double(array, width,max(length,1),initVal,lmp);
 }
 
 void Foam::twoWayM2M::allocateArray
@@ -272,7 +274,7 @@ void Foam::twoWayM2M::allocateArray
 ) const
 {
     //if(length==-1) then LIGGGHTS uses own length data
-    allocate_external_int(array, width,length,initVal,lmp);
+    allocate_external_int(array, width,max(length,1),initVal,lmp);
 }
 
 void Foam::twoWayM2M::allocateArray
@@ -352,19 +354,22 @@ int Foam::twoWayM2M::getNumberOfClumps() const
 
 void Foam::twoWayM2M::syncIDs() const
 {
+    // if I do not creat new communicators it fails at liggghts porcessor transfer!
+    // this should probably not be done that often!!!  
+    delete lmp2foam_;
+    delete lmp2foam_vec_;
+    delete foam2lmp_vec_;
+    lmp2foam_ = new Many2Many(MPI_COMM_WORLD);
+    lmp2foam_vec_ = new Many2Many(MPI_COMM_WORLD);
+    foam2lmp_vec_ = new Many2Many(MPI_COMM_WORLD);
+    //?=======
+
     //MPI_Barrier(MPI_COMM_WORLD);
-    //Pout << "== syncIDs  " << endl;
+    //Pout << couplingStep_ << "st == syncIDs  " << endl;
 
     particleCloud_.clockM().start(5,"recv_DEM_ids");
     // get data from lammps
-    int nlocal_lammps_old=nlocal_lammps_;
     nlocal_lammps_ = *((int *) lammps_extract_global(lmp,"nlocal"));
-
-    // weak form to detect lmp proc change - this will cause problems!!!
-    //bool comRun=false;
-    //if(nlocal_lammps_old != nlocal_lammps_) comRun=true;
-    //MPI_Allreduce(&comRun, &firstRun_, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
-
     int*  id_lammps_sync;
     double** pos_lammps_sync;
     if(firstRun_)  // do not forget iterator !
@@ -380,17 +385,11 @@ void Foam::twoWayM2M::syncIDs() const
         
         Foam::dataExchangeModel::allocateArray(pos_lammps_,-1.,3,nlocal_lammps_);  // do I need this???
         pos_lammps_ = (double **) lammps_extract_atom(lmp,"x");                
-
+    }
-        //for (int i = 0; i < nlocal_lammps_; i++){
+    else
-        //    vector pp;
-        //    for (int j = 0; j < 3; j++) pp[j]=pos_lammps_[i][j];
-        //    Pout << "pp=" << pp << endl;
-        //}
-                
-    }else
     {
        // re-arrange data using map
-        Foam::dataExchangeModel::allocateArray(id_lammps_sync,-1.,nlocal_lammps_);
+        Foam::dataExchangeModel::allocateArray(id_lammps_sync,-1.,nlocal_lammps_); // probably not necessary
         id_lammps_sync = (int *) lammps_extract_atom(lmp,"id");
 
         delete [] id_lammps_vec_;
@@ -404,10 +403,10 @@ void Foam::twoWayM2M::syncIDs() const
         lmp2foam_vec_->setup(nlocal_lammps_*3,id_lammps_vec_,nlocal_foam_*3,id_foam_vec_);
         foam2lmp_vec_->setup(nlocal_foam_*3,id_foam_vec_,nlocal_lammps_*3,id_lammps_vec_);
 
+        // map data according to last TS
         Foam::dataExchangeModel::allocateArray(id_lammps_,-1,nlocal_foam_);  // probably not necessary
         lmp2foam_->exchange(id_lammps_sync, id_lammps_);
 
-
         Foam::dataExchangeModel::allocateArray(pos_lammps_sync,-1.,3,nlocal_lammps_);  // probably not necessary
         pos_lammps_sync = (double **) lammps_extract_atom(lmp,"x");
 
@@ -417,66 +416,14 @@ void Foam::twoWayM2M::syncIDs() const
         Foam::dataExchangeModel::allocateArray(pos_lammps_,-1.,3,nlocal_foam_);
         lmp2foam_vec_->exchange(&(pos_lammps_sync[0][0]), &(gugu[0][0]));
 
-        // try the same with getData
-        //double** gugu;
-        //Foam::dataExchangeModel::allocateArray(gugu,-1.,3,nlocal_foam_);
-        //Foam::dataExchangeModel::allocateArray(pos_lammps_,-1.,3,nlocal_foam_);
-        //getData("x","vector-atom",gugu,1);
-
         // conversion of array (should not be necessary if above problem is solved)
         for (int i = 0; i < nlocal_foam_; i++)
             for (int j = 0; j < 3; j++)
                 pos_lammps_[i][j]=gugu[0][i*3+j];
 
-        /*// test communication - works!!!
-        int* llv;
-        int* ovv;
-        Many2Many* l2f = new Many2Many(MPI_COMM_WORLD);
-        double** gaga;
-        if(Pstream::myProcNo()==0)
-        {
-            Foam::dataExchangeModel::allocateArray(llv,0,6);
-            Foam::dataExchangeModel::allocateArray(ovv,0,1);
-            for(int i=0;i<6;i++)
-                llv[i]=i+1;
-            ovv[0]=-1;
-            l2f->setup(6,llv,0,ovv);
-            l2f->exchange(&(pos_lammps_sync[0][0]), &(gaga[0][0]));
-        }else{
-            Foam::dataExchangeModel::allocateArray(llv,0,1);
-            Foam::dataExchangeModel::allocateArray(ovv,0,6);
-            for(int i=0;i<6;i++)
-                ovv[i]=i+1;
-            llv[0]=-1;
-            l2f->setup(0,llv,6,ovv);
-
-            double** gugu;
-            Foam::dataExchangeModel::allocateArray(gaga,-1.,3,2);
-            Foam::dataExchangeModel::allocateArray(gugu,-1.,3,2);
-            l2f->exchange(&(pos_lammps_sync[0][0]), &(gugu[0][0]));
-
-            // this does not work! -conversion needed!
-            //for (int i = 0; i < 2; i++)
-            //    Pout << "getDataXX: " << "v" <<"=" << gugu[i][0]<<","<< gugu[i][1]<<","<< gugu[i][2]<<endl;
-
-            // conversion of array
-            for (int i = 0; i < 2; i++)
-                for (int j = 0; j < 3; j++)
-                    gaga[i][j]=gugu[0][i*3+j];
-
-            // output
-            for (int i = 0; i < 2; i++)
-                Pout << "getData3: " << "v" <<"=" << gaga[i][0]<<","<< gaga[i][1]<<","<< gaga[i][2]<<endl;
-        }*/
-
     }
     particleCloud_.clockM().stop("recv_DEM_ids");
 
-    // make a copy of liggghts ID array
-    /*delete [] id_lammpsComm_;
-    Foam::dataExchangeModel::allocateArray(id_lammpsComm_,0,nlocal_lammps_);
-    for(int i=0;i<nlocal_lammps_;i++) id_lammpsComm_[i]=id_lammps_[i];*/
-
     particleCloud_.clockM().start(6,"locateParticle()");
     locateParticle();
     particleCloud_.clockM().stop("locateParticle()");
@@ -490,28 +437,28 @@ void Foam::twoWayM2M::syncIDs() const
         {
             if(firstRun_)
             {
-                Pout << "id_lammps_[" << i << "]=" << id_lammps_[i] << "  -  "<<endl;
+                Pout << couplingStep_ << "st id_lammps_[" << i << "]=" << id_lammps_[i] << "  -  "<<endl;
             }else{
-                Pout << "id_lammps_sync[" << i << "]=" << id_lammps_sync[i] << "  -  "<<endl;
+                Pout << couplingStep_ << "st id_lammps_sync[" << i << "]=" << id_lammps_sync[i] << "  -  "<<endl;
             }
         }*/
         /*for (int i = 0; i < nlocal_lammps_*3; i++)
         {
-            Pout << "id_lammps_vec_[" << i << "]=" << id_lammps_vec_[i] << "  -  "<<endl;
+            Pout << couplingStep_ << "st id_lammps_vec_[" << i << "]=" << id_lammps_vec_[i] << "  -  "<<endl;
         }*/
         /*Info << "FOAM "<< endl;
         for (int i = 0; i < nlocal_foam_; i++)
         {
-            Pout << "id_foam_[" << i << "]=" << id_foam_[i] << "  -  "<<endl;
+            Pout << couplingStep_ << "st id_foam_[" << i << "]=" << id_foam_[i] << "  -  "<<endl;
         }*/
         /*for (int i = 0; i < nlocal_foam_*3; i++)
         {
-            Pout << "id_foam_vec_[" << i << "]=" << id_foam_vec_[i] << "  -  "<<endl;
+            Pout << couplingStep_ << "st id_foam_vec_[" << i << "]=" << id_foam_vec_[i] << "  -  "<<endl;
         }*/
         Pout << couplingStep_ << "st nlocal_lammps_=" << nlocal_lammps_ << endl;
         Pout << couplingStep_ << "st nlocal_foam_=" << nlocal_foam_ << endl;
 
-    // communicate lmp->foam
+    // correct mapping
     particleCloud_.clockM().start(11,"setup_Comm");
     if(firstRun_)
     {
@@ -531,14 +478,6 @@ void Foam::twoWayM2M::syncIDs() const
 
 void Foam::twoWayM2M::locateParticle() const
 {
-    //?=======
-    // if I do not creat new communicators it fails at liggghts porcessor transfer!
-    // this should probably not be done that often!!!  
-    lmp2foam_ = new Many2Many(MPI_COMM_WORLD);
-    lmp2foam_vec_ = new Many2Many(MPI_COMM_WORLD);
-    foam2lmp_vec_ = new Many2Many(MPI_COMM_WORLD);
-    //?=======
-
     int nop = particleCloud_.numberOfParticles();
 
     // realloc array of lost particles
@@ -590,9 +529,10 @@ void Foam::twoWayM2M::locateParticle() const
         pos = vector(pos_lammps_[i][0],pos_lammps_[i][1],pos_lammps_[i][2]);
         searchCellID = -1;
         //Pout << "stage1 looking for particle pos="<< pos << endl;
-
         cellID = particleCloud_.locateM().findSingleCell(pos,searchCellID);
     
+        point oldPos(pos_foam_[nlocal_foam_*3+0],pos_foam_[nlocal_foam_*3+1],pos_foam_[nlocal_foam_*3+2]);
+
         // found particle on cfd proc
         if (cellID >= 0)
         {
@@ -608,7 +548,7 @@ void Foam::twoWayM2M::locateParticle() const
             cellID_foam_[nlocal_foam_] = cellID;
 
             nlocal_foam_ += 1;
-            //Pout << "stage1 found particle at pos=" << pos << endl;
+            //Pout << couplingStep_ << "st stage1 found particle at pos=" << pos << " ,id_lammps_[i]=" << id_lammps_[i] << endl;
         }
         else
         {
@@ -618,7 +558,77 @@ void Foam::twoWayM2M::locateParticle() const
                 lost_pos_[nlocal_foam_lost_*3+j] = pos[j];
 
             nlocal_foam_lost_ += 1;
-            //Pout << "cellID="<< cellID << " lost particle id="<< id_lammps_[i] <<" at pos=" << pos << endl;
+            //Pout << couplingStep_ << "st cellID="<< cellID << " lost particle id="<< id_lammps_[i] <<" at pos=" << pos << endl;
+
+            // find out where particle has migrated (must have passed a CFD proc border)
+            /*point newPos=pos;
+            label nearestFace = particleCloud_.locateM().intersection(oldPos,newPos);
+            Pout << "nearest face=" << nearestFace << endl;
+
+            // If we hit a boundary face
+            if (nearestFace >= particleCloud_.mesh().nInternalFaces())
+            {
+                Pout << " face=" << nearestFace << " , is a boundary face!" << endl;
+            }*/
+
+/*  // probably useful proc stuff from Cloud.C
+
+    const polyBoundaryMesh& pbm = pMesh().boundaryMesh();
+    const globalMeshData& pData = polyMesh_.globalData();
+
+    // Which patches are processor patches
+    const labelList& procPatches = pData.processorPatches();
+
+    // Indexing of patches into the procPatches list
+    const labelList& procPatchIndices = pData.processorPatchIndices();
+
+    // Indexing of equivalent patch on neighbour processor into the
+    // procPatches list on the neighbour
+    const labelList& procPatchNeighbours = pData.processorPatchNeighbours();
+
+    // Which processors this processor is connected to
+    const labelList& neighbourProcs = pData[Pstream::myProcNo()];
+
+    // Indexing from the processor number into the neighbourProcs list
+    labelList neighbourProcIndices(Pstream::nProcs(), -1);
+
+    forAll(neighbourProcs, i)
+    {
+        neighbourProcIndices[neighbourProcs[i]] = i;
+    }
+
+
+
+                // If we are running in parallel and the particle is on a
+                // boundary face
+                if (Pstream::parRun() && p.face() >= pMesh().nInternalFaces())
+                {
+                    label patchI = pbm.whichPatch(p.face());
+
+                    // ... and the face is on a processor patch
+                    // prepare it for transfer
+                    if (procPatchIndices[patchI] != -1)
+                    {
+                        label n = neighbourProcIndices
+                        [
+                            refCast<const processorPolyPatch>
+                            (
+                                pbm[patchI]
+                            ).neighbProcNo()
+                        ];
+
+                        p.prepareForParallelTransfer(patchI, td);
+
+                        particleTransferLists[n].append(this->remove(&p));
+
+                        patchIndexTransferLists[n].append
+                        (
+                            procPatchNeighbours[patchI]
+                        );
+                    }
+                }
+*/
+
         }
     }
     particleCloud_.clockM().stop("locate_Stage1");

src/lagrangian/cfdemParticle/subModels/locateModel/locateModel/locateModel.H

@@ -118,6 +118,12 @@ public:
             vector& position,
             label& oldCellID
         ) const = 0;
+
+        virtual label intersection
+        (
+                const point& pStart,
+                const point& pEnd
+        ) const {return -1;};
 };
 
 

src/lagrangian/cfdemParticle/subModels/locateModel/turboEngineSearch/turboEngineSearch.C

@@ -63,12 +63,12 @@ turboEngineSearch::turboEngineSearch
     locateModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
     treeSearch_(propsDict_.lookup("treeSearch")),
+    bb_(particleCloud_.mesh().points(),false),
     #ifdef version16ext
-        searchEngine_(particleCloud_.mesh(),false), //(particleCloud_.mesh(),faceDecomp_)
+        searchEngine_(particleCloud_.mesh(),false) //(particleCloud_.mesh(),faceDecomp_)
     #elif defined(version21)
-        searchEngine_(particleCloud_.mesh(),polyMesh::FACEPLANES), // FACEPLANES or FACECENTRETETS; FACEDIAGTETS not stable
+        searchEngine_(particleCloud_.mesh(),polyMesh::FACEPLANES) // FACEPLANES or FACECENTRETETS; FACEDIAGTETS not stable
     #endif
-    bb_(particleCloud_.mesh().points(),false)
 {}
 
 

src/lagrangian/cfdemParticle/subModels/locateModel/turboEngineSearch/turboEngineSearch.H

@@ -65,10 +65,12 @@ private:
 
     Switch treeSearch_;
 
-    meshSearch searchEngine_;
-
     boundBox bb_;
 
+protected:
+
+    meshSearch searchEngine_;
+
 public:
 
     //- Runtime type information

src/lagrangian/cfdemParticle/subModels/locateModel/turboEngineSearchM2M/turboEngineSearchM2M.C

@@ -85,6 +85,20 @@ label turboEngineSearchM2M::findCell
     return 1;
 }
 
+label turboEngineSearchM2M::intersection
+(
+    const point& pStart,
+    const point& pEnd
+) const
+{
+    // find intersection with boundary
+    //pointIndexHit hit=searchEngine_.intersection(pStart,pEnd);
+    //Info << "hit.index()=" << hit.index()<< endl;
+    
+    //return searchEngine_.findNearestBoundaryFace(pStart);
+    return searchEngine_.intersection(pStart,pEnd).index();
+}
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/lagrangian/cfdemParticle/subModels/locateModel/turboEngineSearchM2M/turboEngineSearchM2M.H

@@ -89,6 +89,12 @@ public:
             int size
         ) const;
 
+        label intersection
+        (
+                const point& pStart,
+                const point& pEnd
+        ) const;
+
 };