release on 2012-12-03_09-41-15

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.
+\*---------------------------------------------------------------------------*/

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

@@ -61,11 +61,19 @@ twoWayM2M::twoWayM2M
 )
 :
     dataExchangeModel(dict,sm),
-    propsDict_(dict.subDict(typeName + "Props"))/*,
+    propsDict_(dict.subDict(typeName + "Props")),
-    lmp2foam_(*new Many2Many(MPI_COMM_WORLD)), // init of many2many &
+    pbm_(sm.mesh().boundaryMesh()), 
+    pData_(sm.mesh().globalData()),
+    procPatches_(pData_.processorPatches()),
+    procPatchIndices_(pData_.processorPatchIndices()),
+    neighbourProcs_(pData_[Pstream::myProcNo()]),
+    neighbourProcIndices_(Pstream::nProcs(), -1)
+/*  lmp2foam_(*new Many2Many(MPI_COMM_WORLD)), // init of many2many &
     lmp2foam_vec_(*new Many2Many(MPI_COMM_WORLD)),
     foam2lmp_vec_(*new Many2Many(MPI_COMM_WORLD))*/
 {
+    forAll(neighbourProcs_, i) neighbourProcIndices_[neighbourProcs_[i]] = i;
+
     Info<<"Starting up LIGGGHTS for first time execution"<<endl;
 
     MPI_Comm_rank(MPI_COMM_WORLD,&me);
@@ -361,28 +369,21 @@ int Foam::twoWayM2M::getNumberOfClumps() const
 
 void Foam::twoWayM2M::syncIDs() const
 {
-    // if I do not creat new communicators it fails at liggghts porcessor transfer!
+    particleCloud_.clockM().start(5,"recv_DEM_ids");
-    // this should probably not be done that often!!!  
+
+    // update communication 
     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 << couplingStep_ << "st == syncIDs  " << endl;
-//if(couplingStep_==30){
-//FatalError<<"stop!!!"<< abort(FatalError);
-//}
-
-    particleCloud_.clockM().start(5,"recv_DEM_ids");
     // get data from lammps
     nlocal_lammps_ = *((int *) lammps_extract_global(lmp,"nlocal"));
     int*  id_lammps_sync;
     double** pos_lammps_sync;
-    if(firstRun_)  // do not forget iterator !
+    if(firstRun_)
     {
         // IDs for vectors
         if(nlocal_lammps_>0){
@@ -393,9 +394,6 @@ void Foam::twoWayM2M::syncIDs() const
             id_lammps_[0]=0;
         }
 
-        Pout << couplingStep_ << "st id_lammps_[0]=" << id_lammps_[0]<< endl;
-
-
         //delete [] id_lammps_vec_;
         Foam::dataExchangeModel::allocateArray(id_lammps_vec_,0,nlocal_lammps_*3);
         for (int i = 0; i < nlocal_lammps_; i++)
@@ -470,6 +468,12 @@ void Foam::twoWayM2M::syncIDs() const
         //for (int i = 0; i < nlocal_lammps_; i++)
         //    Pout << "getData3: " << "v" <<"=" << pos_lammps_[i][0]<<","<<pos_lammps_[i][1]<<","<<pos_lammps_[i][2] <<endl;
 
+//MPI_Barrier(MPI_COMM_WORLD);
+//Pout << couplingStep_ << "st == syncIDs  " << endl;
+//if(couplingStep_==30){
+//FatalError<<"stop!!!"<< abort(FatalError);
+//}
+
         // output
         /*Info << "LAMMPS " << endl;
         for (int i = 0; i < nlocal_lammps_; i++)
@@ -494,18 +498,20 @@ void Foam::twoWayM2M::syncIDs() const
         {
             Pout << couplingStep_ << "st id_foam_vec_[" << i << "]=" << id_foam_vec_[i] << "  -  "<<endl;
         }*/
-        Pout << couplingStep_ << "st nlocal_lammps_=" << nlocal_lammps_ << endl;
+        //Pout << couplingStep_ << "st nlocal_lammps_=" << nlocal_lammps_ << endl;
-        Pout << couplingStep_ << "st nlocal_foam_=" << nlocal_foam_ << endl;
+        //Pout << couplingStep_ << "st nlocal_foam_=" << nlocal_foam_ << endl;
 
+    // correct mapping
+    particleCloud_.clockM().start(11,"setup_Comm");
+
+    // update communication
     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);
-
+ 
-    // correct mapping
-    particleCloud_.clockM().start(11,"setup_Comm");
     if(firstRun_)
     {
         //lmp2foam_->setup(nlocal_lammps_,id_lammpsComm_,nlocal_foam_,id_foam_);
@@ -557,7 +563,8 @@ void Foam::twoWayM2M::locateParticle() const
         }
     }
 
-    // loop all lmp particles
+    // stage 1 - look on proc or send or prepare for all-to-all
+    particleCloud_.clockM().start(7,"locate_Stage1");
     int iterate;
     if(firstRun_) iterate=nlocal_lammps_;
     else iterate=nlocal_foam_;
@@ -566,39 +573,13 @@ void Foam::twoWayM2M::locateParticle() const
     nlocal_foam_lost_ = 0;
     vector pos;
     label cellID = 0;
-    label searchCellID;
+    label searchCellID=-1;
-
+    List< DynamicList<int> > particleTransferID(neighbourProcs_.size());
-    particleCloud_.clockM().start(7,"locate_Stage1");
+    List< DynamicList<vector> > particleTransferPos(neighbourProcs_.size());
-
-/*//===============
-    // move this to top level
-            const polyBoundaryMesh& pbm = particleCloud_.mesh().boundaryMesh(); 
-            const globalMeshData& pData = particleCloud_.mesh().globalData(); // polyMesh???
-
-            // Which patches are processor patches
-            const labelList& procPatches = pData.processorPatches();
-
-            // Indexing of patches into the procPatches list
-            const labelList& procPatchIndices = pData.processorPatchIndices();
-
-            // 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;
-            }
-
-            List< DynamicList<int> > particleTransferLists(neighbourProcs.size());
-//===============*/
 
     for (int i = 0; i < iterate; i++)
     {
         pos = vector(pos_lammps_[i][0],pos_lammps_[i][1],pos_lammps_[i][2]);
-        searchCellID = -1;
         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]);
 
@@ -621,98 +602,183 @@ void Foam::twoWayM2M::locateParticle() const
         }
         else
         {
-            //-----------------
+            // find out where particle has migrated (must have passed a CFD proc border)
-            id_foam_lost_[nlocal_foam_lost_] = id_lammps_[i];
+            bool commPart=false;
-          
-            for (int j=0; j<3; j++)
-                lost_pos_[nlocal_foam_lost_*3+j] = pos[j];
-
-            nlocal_foam_lost_ += 1;
-            //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;
+            //Pout << couplingStep_ << "st nearestFace="<< nearestFace << " oldPos="<< oldPos <<" at pos=" << pos << endl;
 
-
-            // If we hit a boundary face
             if (nearestFace >= particleCloud_.mesh().nInternalFaces())
             {
-                Pout << " face=" << nearestFace << " , is a boundary face!" << endl;
+                label patchI = pbm_.whichPatch(nearestFace);
-                label patchI = pbm.whichPatch(nearestFace);
 
-                // ... and the face is on a processor patch
+                if (procPatchIndices_[patchI] != -1)
-                // prepare it for transfer
-                if (procPatchIndices[patchI] != -1)
                 {
-                    Pout << " face=" << nearestFace << " , is a proc face!" << endl;
+                    label n = neighbourProcIndices_
-                    label n = neighbourProcIndices
                     [
                         refCast<const processorPolyPatch>
                         (
-                            pbm[patchI]
+                            pbm_[patchI]
                         ).neighbProcNo()
                     ];
-                    Pout << " communicate to n=" << n << endl;
+                    particleTransferID[n].append(id_lammps_[i]);
-                    particleTransferLists[n].append(id_lammps_[i]);
+                    particleTransferPos[n].append(pos);
+                    commPart=true;
                 }
-            }*/
+            }
-            //-----------------
+            if (!commPart)
+            {
+                // prepare for all to all comm
+                id_foam_lost_[nlocal_foam_lost_] = id_lammps_[i];
+          
+                for (int j=0; j<3; j++)
+                    lost_pos_[nlocal_foam_lost_*3+j] = pos[j];
 
+                nlocal_foam_lost_ += 1;
+                //Pout << couplingStep_ << "st cellID="<< cellID << " lost particle id="<< id_lammps_[i] <<" at pos=" << pos << endl;
+            }
+        }
+    }
+    particleCloud_.clockM().stop("locate_Stage1");
+
+    // stage 2 - recv particle, locate or all-to-all
+    particleCloud_.clockM().start(9,"locate_Stage2");
+ 
+    // Allocate transfer buffers
+    PstreamBuffers pBufs(Pstream::nonBlocking);
+
+    // Stream into send buffers
+    forAll(particleTransferID, i)
+    {
+        if (particleTransferID[i].size())
+        {
+            UOPstream particleStream
+            (
+                neighbourProcs_[i],
+                pBufs
+            );
+            particleStream << particleTransferID[i]<<particleTransferPos[i];
         }
     }
 
-    /*forAll(particleTransferLists, i)
+    // Set up transfers when in non-blocking mode. Returns sizes (in bytes)
+    // to be sent/received.
+    labelListList allNTrans(Pstream::nProcs());
+
+    pBufs.finishedSends(allNTrans);
+
+    forAll(allNTrans, i)
     {
-        forAll(particleTransferLists[i],j)
+        forAll(allNTrans[i], j)
         {
-            Pout << "communicate particle   i="<< particleTransferLists[i][j]<<" to proc "<< i << endl;
+            if (allNTrans[i][j])
-        }
-    }*/
-
-    particleCloud_.clockM().stop("locate_Stage1");
-
-    // using allgather to allreduce lost particles
-    particleCloud_.clockM().start(8,"locate_Stage2");
-
-    int nlocal_foam_lost_all = LAMMPS_NS::MPI_Allgather_Vector(lost_pos_, nlocal_foam_lost_*3, lost_pos_all, MPI_COMM_WORLD)/3;
-    LAMMPS_NS::MPI_Allgather_Vector(id_foam_lost_, nlocal_foam_lost_, id_foam_lost_all, MPI_COMM_WORLD);
-    Info << couplingStep_ << "st nlocal_foam_lost_all=" << nlocal_foam_lost_all << endl;
-
-    // locate lost particles
-    for (int i = 0; i < nlocal_foam_lost_all; i++)
-    {
-        pos = vector(lost_pos_all[i*3+0],lost_pos_all[i*3+1],lost_pos_all[i*3+2]);
-        //Pout << "stage2 looking for particle pos="<< pos << endl;
-        searchCellID = -1;
-        particleCloud_.clockM().start(9,"findSingleCell");   
-        cellID = particleCloud_.locateM().findSingleCell(pos,searchCellID);
-        particleCloud_.clockM().stop("findSingleCell");
-        
-        // found particle on cfd proc
-        if (cellID >= 0)
-        {
-            // IDs for scalars
-            id_foam_[nlocal_foam_] = id_foam_lost_all[i];
-
-            // IDs for vectors
-            for (int j=0;j<3;j++)
             {
-                id_foam_vec_[nlocal_foam_*3+j] = id_foam_lost_all[i]*3+j;
+                break;
-                pos_foam_[nlocal_foam_*3+j] = pos[j];
             }
-            cellID_foam_[nlocal_foam_] = cellID;
+        }
+    }
 
-            // mark that ID was finally found
+    // Retrieve from receive buffers
-            id_foam_lost_all[i]=-1;
+    label neighbProci;
+    label nRec;
+    forAll(neighbourProcs_, i)
+    {
+        neighbProci = neighbourProcs_[i];
+        nRec = allNTrans[neighbProci][Pstream::myProcNo()];
 
-            nlocal_foam_ += 1;
+        if (nRec)
-            //Pout << "stage2 found particle at pos=" << pos << endl;
+        {
+            UIPstream particleStream(neighbProci, pBufs);
+
+            labelList recvParticleTransferID(particleStream);
+            List<vector> recvParticleTransferPos(particleStream);
+
+            forAll(recvParticleTransferID,i)
+            {
+                //Pout << "received id="<<recvParticleTransferID[i] << "received pos="<<recvParticleTransferPos[i] << endl;
+                pos = recvParticleTransferPos[i];
+                cellID = particleCloud_.locateM().findSingleCell(pos,searchCellID);
+        
+                // found particle on cfd proc
+                if (cellID >= 0)
+                {
+                    // IDs for scalars
+                    id_foam_[nlocal_foam_] = recvParticleTransferID[i];
+
+                    // IDs for vectors
+                    for (int j=0;j<3;j++)
+                    {
+                        id_foam_vec_[nlocal_foam_*3+j] = recvParticleTransferID[i]*3+j;
+                        pos_foam_[nlocal_foam_*3+j] = pos[j];
+                    }
+                    cellID_foam_[nlocal_foam_] = cellID;
+
+                    // mark that ID was finally found
+                    //id_foam_lost_all[i]=-1;
+
+                    nlocal_foam_ += 1;
+                }
+                else // might have been comm to wrong proc
+                {
+                    // prepare for all to all comm
+                    id_foam_lost_[nlocal_foam_lost_] = recvParticleTransferID[i];
+          
+                    for (int j=0; j<3; j++)
+                        lost_pos_[nlocal_foam_lost_*3+j] = pos[j];
+
+                    nlocal_foam_lost_ += 1;
+                }
+
+            }
         }
     }
     particleCloud_.clockM().stop("locate_Stage2");
 
+    // stage 3 - all-to-all
+    particleCloud_.clockM().start(9,"locate_Stage3");
+
+    // check if all-to-all is necessary
+    int nlocal_foam_lost_all;
+    MPI_Allreduce(&nlocal_foam_lost_, &nlocal_foam_lost_all, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+
+    if (nlocal_foam_lost_all > 0)
+    {
+        //Info << "all-to-all necessary: nlocal_foam_lost_all=" << nlocal_foam_lost_all << endl;
+        int nlocal_foam_lost_all = LAMMPS_NS::MPI_Allgather_Vector(lost_pos_, nlocal_foam_lost_*3, lost_pos_all, MPI_COMM_WORLD)/3;
+        LAMMPS_NS::MPI_Allgather_Vector(id_foam_lost_, nlocal_foam_lost_, id_foam_lost_all, MPI_COMM_WORLD);
+        //Info << couplingStep_ << "st nlocal_foam_lost_all=" << nlocal_foam_lost_all << endl;
+
+        // locate lost particles
+        for (int i = 0; i < nlocal_foam_lost_all; i++)
+        {
+            pos = vector(lost_pos_all[i*3+0],lost_pos_all[i*3+1],lost_pos_all[i*3+2]);
+            //Pout << "stage3 look for particle at pos=" << pos << endl;
+            cellID = particleCloud_.locateM().findSingleCell(pos,searchCellID);
+        
+            // found particle on cfd proc
+            if (cellID >= 0)
+            {
+                // IDs for scalars
+                id_foam_[nlocal_foam_] = id_foam_lost_all[i];
+
+                // IDs for vectors
+                for (int j=0;j<3;j++)
+                {
+                    id_foam_vec_[nlocal_foam_*3+j] = id_foam_lost_all[i]*3+j;
+                    pos_foam_[nlocal_foam_*3+j] = pos[j];
+                }
+                cellID_foam_[nlocal_foam_] = cellID;
+
+                // mark that ID was finally found
+                id_foam_lost_all[i]=-1;
+
+                nlocal_foam_ += 1;
+                //Pout << "stage3 found particle at pos=" << pos << " ,id="<< id_foam_lost_all[i] << endl;
+            }
+        }
+    }
+    particleCloud_.clockM().stop("locate_Stage3");
+
   /*  // check if really all particles were found
     particleCloud_.clockM().start(10,"locate_Stage3");
     Foam::dataExchangeModel::allocateArray(id_foam_nowhere_all,1,nlocal_foam_lost_all);
@@ -759,20 +825,6 @@ Info << "nlocal_lammps_ALL=" << gugu << endl;*/
     //delete[] lost_pos_all;
 }
 
-/*void Foam::twoWayM2M::exchange(double* demDat, double* cfdDat) const
-{
-    lmp2foam_->exchange(demDat,cfdDat);//(pos_lammps,pos_foam);
-}*/
-
-/*template <typename T>
-T*& Foam::twoWayM2M::extract_save(T *& a,char name)
-{
-    a = (T *) lammps_extract_atom(lmp,name);
-    if(!a)
-        a=new T[1];
-    return a;//static_const<T>(a); 
-}*/
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

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

@@ -111,6 +111,13 @@ private:
     mutable int *cellID_foam_;
     mutable double *pos_foam_;
 
+    const polyBoundaryMesh& pbm_; 
+    const globalMeshData& pData_;
+    const labelList& procPatches_;
+    const labelList& procPatchIndices_;
+    const labelList& neighbourProcs_;
+    labelList neighbourProcIndices_;
+
     // variables
     int me;
 
@@ -191,9 +198,6 @@ public:
 
         void syncIDs() const;
         void locateParticle() const;
-        //void exchange(double*,double*) const;
-        //template <typename T>
-        //T*& extract_save(T *& a,char name);
 };
 
 

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

@@ -92,11 +92,25 @@ label turboEngineSearchM2M::intersection
 ) const
 {
     // find intersection with boundary
-    //pointIndexHit hit=searchEngine_.intersection(pStart,pEnd);
+    label face = searchEngine_.intersection(pStart,pEnd).index();
-    //Info << "hit.index()=" << hit.index()<< endl;
+
-    
+    // try alternative
-    //return searchEngine_.findNearestBoundaryFace(pStart);
+    if (face==-1)
-    return searchEngine_.intersection(pStart,pEnd).index();
+    {
+        // try alternative
+        face = searchEngine_.findNearestBoundaryFace(pEnd);
+        //Pout << "found face=" << face << " with findNearestBoundaryFace" << endl;
+
+        // might have been first search
+        if (face==-1 && mag(pStart-point(0,0,0))<SMALL)
+        {
+            point pStart2 = pEnd+0.0001*(pStart-pEnd)/mag(pStart-pEnd);
+            face = searchEngine_.intersection(pStart2,pEnd).index();
+            //Pout << "found face=" << face << " with pStart2="<< pStart2 << endl;
+        }
+
+    }
+    return face;
 }
 
 

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/constant/couplingProperties

@@ -32,7 +32,7 @@ couplingInterval 100;
 
 voidFractionModel divided;
 
-locateModel engine;//turboEngineM2M;//
+locateModel turboEngineM2M;//engine;//
 
 meshMotionModel noMeshMotion;
 
@@ -40,7 +40,7 @@ regionModel allRegion;
 
 IOModel "off";
 
-dataExchangeModel twoWayMPI;//twoWayM2M;//twoWayFiles;//oneWayVTK;//
+dataExchangeModel twoWayM2M;//twoWayMPI;//twoWayFiles;//oneWayVTK;//
 
 averagingModel dense;//dilute;//