src/lagrangian/basic/Cloud/Cloud.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "Cloud.H"
  27 #include "processorPolyPatch.H"
  28 #include "globalMeshData.H"
  29 #include "PstreamCombineReduceOps.H"
  30 #include "mapPolyMesh.H"
  31 #include "Time.H"
  32 #include "OFstream.H"
  33 #include "wallPolyPatch.H"
  34
  35 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
  36
  37 template<class ParticleType>
  38 void Foam::Cloud<ParticleType>::calcCellWallFaces() const
  39 {
  40     cellWallFacesPtr_.reset(new PackedBoolList(pMesh().nCells(), false));
  41
  42     PackedBoolList& cellWallFaces = cellWallFacesPtr_();
  43
  44     const polyBoundaryMesh& patches = polyMesh_.boundaryMesh();
  45
  46     forAll(patches, patchI)
  47     {
  48         if (isA<wallPolyPatch>(patches[patchI]))
  49         {
  50             const polyPatch& patch = patches[patchI];
  51
  52             const labelList& pFaceCells = patch.faceCells();
  53
  54             forAll(pFaceCells, pFCI)
  55             {
  56                 cellWallFaces[pFaceCells[pFCI]] = true;
  57             }
  58         }
  59     }
  60 }
  61
  62
  63 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  64
  65 template<class ParticleType>
  66 Foam::Cloud<ParticleType>::Cloud
  67 (
  68     const polyMesh& pMesh,
  69     const IDLList<ParticleType>& particles
  70 )
  71 :
  72     cloud(pMesh),
  73     IDLList<ParticleType>(),
  74     polyMesh_(pMesh),
  75     labels_(),
  76     nTrackingRescues_(),
  77     cellWallFacesPtr_()
  78 {
  79     // Ask for the tetBasePtIs to trigger all processors to build
  80     // them, otherwise, if some processors have no particles then
  81     // there is a comms mismatch.
  82     polyMesh_.tetBasePtIs();
  83
  84     IDLList<ParticleType>::operator=(particles);
  85 }
  86
  87
  88 template<class ParticleType>
  89 Foam::Cloud<ParticleType>::Cloud
  90 (
  91     const polyMesh& pMesh,
  92     const word& cloudName,
  93     const IDLList<ParticleType>& particles
  94 )
  95 :
  96     cloud(pMesh, cloudName),
  97     IDLList<ParticleType>(),
  98     polyMesh_(pMesh),
  99     labels_(),
 100     nTrackingRescues_(),
 101     cellWallFacesPtr_()
 102 {
 103     // Ask for the tetBasePtIs to trigger all processors to build
 104     // them, otherwise, if some processors have no particles then
 105     // there is a comms mismatch.
 106     polyMesh_.tetBasePtIs();
 107
 108     IDLList<ParticleType>::operator=(particles);
 109 }
 110
 111
 112 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 113
 114 template<class ParticleType>
 115 const Foam::PackedBoolList& Foam::Cloud<ParticleType>::cellHasWallFaces()
 116 const
 117 {
 118     if (!cellWallFacesPtr_.valid())
 119     {
 120         calcCellWallFaces();
 121     }
 122
 123     return cellWallFacesPtr_();
 124 }
 125
 126
 127 template<class ParticleType>
 128 void Foam::Cloud<ParticleType>::addParticle(ParticleType* pPtr)
 129 {
 130     this->append(pPtr);
 131 }
 132
 133
 134 template<class ParticleType>
 135 void Foam::Cloud<ParticleType>::deleteParticle(ParticleType& p)
 136 {
 137     delete(this->remove(&p));
 138 }
 139
 140
 141 template<class ParticleType>
 142 void Foam::Cloud<ParticleType>::cloudReset(const Cloud<ParticleType>& c)
 143 {
 144     // Reset particle cound and particles only
 145     // - not changing the cloud object registry or reference to the polyMesh
 146     ParticleType::particleCount_ = 0;
 147     IDLList<ParticleType>::operator=(c);
 148 }
 149
 150
 151 template<class ParticleType>
 152 template<class TrackData>
 153 void Foam::Cloud<ParticleType>::move(TrackData& td, const scalar trackTime)
 154 {
 155     const polyBoundaryMesh& pbm = pMesh().boundaryMesh();
 156     const globalMeshData& pData = polyMesh_.globalData();
 157
 158     // Which patches are processor patches
 159     const labelList& procPatches = pData.processorPatches();
 160
 161     // Indexing of patches into the procPatches list
 162     const labelList& procPatchIndices = pData.processorPatchIndices();
 163
 164     // Indexing of equivalent patch on neighbour processor into the
 165     // procPatches list on the neighbour
 166     const labelList& procPatchNeighbours = pData.processorPatchNeighbours();
 167
 168     // Which processors this processor is connected to
 169     const labelList& neighbourProcs = pData[Pstream::myProcNo()];
 170
 171     // Indexing from the processor number into the neighbourProcs list
 172     labelList neighbourProcIndices(Pstream::nProcs(), -1);
 173
 174     forAll(neighbourProcs, i)
 175     {
 176         neighbourProcIndices[neighbourProcs[i]] = i;
 177     }
 178
 179     // Initialise the stepFraction moved for the particles
 180     forAllIter(typename Cloud<ParticleType>, *this, pIter)
 181     {
 182         pIter().stepFraction() = 0;
 183     }
 184
 185     // Reset nTrackingRescues
 186     nTrackingRescues_ = 0;
 187
 188     // While there are particles to transfer
 189     while (true)
 190     {
 191         // List of lists of particles to be transfered for all of the
 192         // neighbour processors
 193         List<IDLList<ParticleType> > particleTransferLists
 194         (
 195             neighbourProcs.size()
 196         );
 197
 198         // List of destination processorPatches indices for all of the
 199         // neighbour processors
 200         List<DynamicList<label> > patchIndexTransferLists
 201         (
 202             neighbourProcs.size()
 203         );
 204
 205         // Loop over all particles
 206         forAllIter(typename Cloud<ParticleType>, *this, pIter)
 207         {
 208             ParticleType& p = pIter();
 209
 210             // Move the particle
 211             bool keepParticle = p.move(td, trackTime);
 212
 213             // If the particle is to be kept
 214             // (i.e. it hasn't passed through an inlet or outlet)
 215             if (keepParticle)
 216             {
 217                 // If we are running in parallel and the particle is on a
 218                 // boundary face
 219                 if (Pstream::parRun() && p.face() >= pMesh().nInternalFaces())
 220                 {
 221                     label patchI = pbm.whichPatch(p.face());
 222
 223                     // ... and the face is on a processor patch
 224                     // prepare it for transfer
 225                     if (procPatchIndices[patchI] != -1)
 226                     {
 227                         label n = neighbourProcIndices
 228                         [
 229                             refCast<const processorPolyPatch>
 230                             (
 231                                 pbm[patchI]
 232                             ).neighbProcNo()
 233                         ];
 234
 235                         p.prepareForParallelTransfer(patchI, td);
 236
 237                         particleTransferLists[n].append(this->remove(&p));
 238
 239                         patchIndexTransferLists[n].append
 240                         (
 241                             procPatchNeighbours[patchI]
 242                         );
 243                     }
 244                 }
 245             }
 246             else
 247             {
 248                 deleteParticle(p);
 249             }
 250         }
 251
 252         if (!Pstream::parRun())
 253         {
 254             break;
 255         }
 256
 257         // Allocate transfer buffers
 258         PstreamBuffers pBufs(Pstream::nonBlocking);
 259
 260         // Stream into send buffers
 261         forAll(particleTransferLists, i)
 262         {
 263             if (particleTransferLists[i].size())
 264             {
 265                 UOPstream particleStream
 266                 (
 267                     neighbourProcs[i],
 268                     pBufs
 269                 );
 270
 271                 particleStream
 272                     << patchIndexTransferLists[i]
 273                     << particleTransferLists[i];
 274             }
 275         }
 276
 277         // Set up transfers when in non-blocking mode. Returns sizes (in bytes)
 278         // to be sent/received.
 279         labelListList allNTrans(Pstream::nProcs());
 280
 281         pBufs.finishedSends(allNTrans);
 282
 283         bool transfered = false;
 284
 285         forAll(allNTrans, i)
 286         {
 287             forAll(allNTrans[i], j)
 288             {
 289                 if (allNTrans[i][j])
 290                 {
 291                     transfered = true;
 292                     break;
 293                 }
 294             }
 295         }
 296
 297         if (!transfered)
 298         {
 299             break;
 300         }
 301
 302         // Retrieve from receive buffers
 303         forAll(neighbourProcs, i)
 304         {
 305             label neighbProci = neighbourProcs[i];
 306
 307             label nRec = allNTrans[neighbProci][Pstream::myProcNo()];
 308
 309             if (nRec)
 310             {
 311                 UIPstream particleStream(neighbProci, pBufs);
 312
 313                 labelList receivePatchIndex(particleStream);
 314
 315                 IDLList<ParticleType> newParticles
 316                 (
 317                     particleStream,
 318                     typename ParticleType::iNew(polyMesh_)
 319                 );
 320
 321                 label pI = 0;
 322
 323                 forAllIter(typename Cloud<ParticleType>, newParticles, newpIter)
 324                 {
 325                     ParticleType& newp = newpIter();
 326
 327                     label patchI = procPatches[receivePatchIndex[pI++]];
 328
 329                     newp.correctAfterParallelTransfer(patchI, td);
 330
 331                     addParticle(newParticles.remove(&newp));
 332                 }
 333             }
 334         }
 335     }
 336
 337     if (cloud::debug)
 338     {
 339         reduce(nTrackingRescues_, sumOp<label>());
 340
 341         if (nTrackingRescues_ > 0)
 342         {
 343             Info<< nTrackingRescues_ << " tracking rescue corrections" << endl;
 344         }
 345     }
 346 }
 347
 348
 349 template<class ParticleType>
 350 template<class TrackData>
 351 void Foam::Cloud<ParticleType>::autoMap
 352 (
 353     TrackData& td,
 354     const mapPolyMesh& mapper
 355 )
 356 {
 357     if (cloud::debug)
 358     {
 359         Info<< "Cloud<ParticleType>::autoMap(TrackData&, const mapPolyMesh&) "
 360             << "for lagrangian cloud " << cloud::name() << endl;
 361     }
 362
 363     const labelList& reverseCellMap = mapper.reverseCellMap();
 364     const labelList& reverseFaceMap = mapper.reverseFaceMap();
 365
 366     // Reset stored data that relies on the mesh
 367 //    polyMesh_.clearCellTree();
 368     cellWallFacesPtr_.clear();
 369
 370     forAllIter(typename Cloud<ParticleType>, *this, pIter)
 371     {
 372         ParticleType& p = pIter();
 373
 374         if (reverseCellMap[p.cell()] >= 0)
 375         {
 376             p.cell() = reverseCellMap[p.cell()];
 377
 378             if (p.face() >= 0 && reverseFaceMap[p.face()] >= 0)
 379             {
 380                 p.face() = reverseFaceMap[p.face()];
 381             }
 382             else
 383             {
 384                 p.face() = -1;
 385             }
 386
 387             p.initCellFacePt();
 388         }
 389         else
 390         {
 391             label trackStartCell = mapper.mergedCell(p.cell());
 392
 393             if (trackStartCell < 0)
 394             {
 395                 trackStartCell = 0;
 396             }
 397
 398             vector pos = p.position();
 399
 400             const_cast<vector&>(p.position()) =
 401                 polyMesh_.cellCentres()[trackStartCell];
 402
 403             p.stepFraction() = 0;
 404
 405             p.initCellFacePt();
 406
 407             p.track(pos, td);
 408         }
 409     }
 410 }
 411
 412
 413 template<class ParticleType>
 414 void Foam::Cloud<ParticleType>::writePositions() const
 415 {
 416     OFstream pObj
 417     (
 418         this->db().time().path()/this->name() + "_positions.obj"
 419     );
 420
 421     forAllConstIter(typename Cloud<ParticleType>, *this, pIter)
 422     {
 423         const ParticleType& p = pIter();
 424         pObj<< "v " << p.position().x() << " " << p.position().y() << " "
 425             << p.position().z() << nl;
 426     }
 427
 428     pObj.flush();
 429 }
 430
 431
 432 // * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
 433
 434 #include "CloudIO.C"
 435
 436 // ************************************************************************* //