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) 1991-2009 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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "Cloud.H"
  28 #include "processorPolyPatch.H"
  29 #include "globalMeshData.H"
  30 #include "PstreamCombineReduceOps.H"
  31 #include "mapPolyMesh.H"
  32 #include "Time.H"
  33 #include "OFstream.H"
  34
  35 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  36
  37 template<class ParticleType>
  38 Foam::Cloud<ParticleType>::Cloud
  39 (
  40     const polyMesh& pMesh,
  41     const IDLList<ParticleType>& particles
  42 )
  43 :
  44     cloud(pMesh),
  45     IDLList<ParticleType>(),
  46     polyMesh_(pMesh),
  47     particleCount_(0)
  48 {
  49     IDLList<ParticleType>::operator=(particles);
  50 }
  51
  52
  53 template<class ParticleType>
  54 Foam::Cloud<ParticleType>::Cloud
  55 (
  56     const polyMesh& pMesh,
  57     const word& cloudName,
  58     const IDLList<ParticleType>& particles
  59 )
  60 :
  61     cloud(pMesh, cloudName),
  62     IDLList<ParticleType>(),
  63     polyMesh_(pMesh),
  64     particleCount_(0)
  65 {
  66     IDLList<ParticleType>::operator=(particles);
  67 }
  68
  69
  70 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  71
  72 template<class ParticleType>
  73 Foam::label Foam::Cloud<ParticleType>::getNewParticleID() const
  74 {
  75     label id = particleCount_++;
  76
  77     if (id == labelMax)
  78     {
  79         WarningIn("Cloud<ParticleType>::getNewParticleID() const")
  80             << "Particle counter has overflowed. This might cause problems"
  81             << " when reconstructing particle tracks." << endl;
  82     }
  83     return id;
  84 }
  85
  86
  87 template<class ParticleType>
  88 void Foam::Cloud<ParticleType>::addParticle(ParticleType* pPtr)
  89 {
  90     append(pPtr);
  91 }
  92
  93
  94 template<class ParticleType>
  95 void Foam::Cloud<ParticleType>::deleteParticle(ParticleType& p)
  96 {
  97     delete(this->remove(&p));
  98 }
  99
 100
 101 namespace Foam
 102 {
 103
 104 class combineNsTransPs
 105 {
 106
 107 public:
 108
 109     void operator()(labelListList& x, const labelListList& y) const
 110     {
 111         forAll(y, i)
 112         {
 113             if (y[i].size())
 114             {
 115                 x[i] = y[i];
 116             }
 117         }
 118     }
 119 };
 120
 121 } // End namespace Foam
 122
 123
 124 template<class ParticleType>
 125 template<class TrackingData>
 126 void Foam::Cloud<ParticleType>::move(TrackingData& td)
 127 {
 128     const globalMeshData& pData = polyMesh_.globalData();
 129     const labelList& processorPatches = pData.processorPatches();
 130     const labelList& processorPatchIndices = pData.processorPatchIndices();
 131     const labelList& processorPatchNeighbours =
 132         pData.processorPatchNeighbours();
 133
 134     // Initialise the setpFraction moved for the particles
 135     forAllIter(typename Cloud<ParticleType>, *this, pIter)
 136     {
 137         pIter().stepFraction() = 0;
 138     }
 139
 140     // Assume there will be particles to transfer
 141     bool transfered = true;
 142
 143     // While there are particles to transfer
 144     while (transfered)
 145     {
 146         // List of lists of particles to be transfered for all the processor
 147         // patches
 148         List<IDLList<ParticleType> > transferList(processorPatches.size());
 149
 150         // Loop over all particles
 151         forAllIter(typename Cloud<ParticleType>, *this, pIter)
 152         {
 153             ParticleType& p = pIter();
 154
 155             // Move the particle
 156             bool keepParticle = p.move(td);
 157
 158             // If the particle is to be kept
 159             // (i.e. it hasn't passed through an inlet or outlet)
 160             if (keepParticle)
 161             {
 162                 // If we are running in parallel and the particle is on a
 163                 // boundary face
 164                 if (Pstream::parRun() && p.facei_ >= pMesh().nInternalFaces())
 165                 {
 166                     label patchi = pMesh().boundaryMesh().whichPatch(p.facei_);
 167                     label n = processorPatchIndices[patchi];
 168
 169                     // ... and the face is on a processor patch
 170                     // prepare it for transfer
 171                     if (n != -1)
 172                     {
 173                         p.prepareForParallelTransfer(patchi, td);
 174                         transferList[n].append(this->remove(&p));
 175                     }
 176                 }
 177             }
 178             else
 179             {
 180                 deleteParticle(p);
 181             }
 182         }
 183
 184         if (Pstream::parRun())
 185         {
 186             // List of the numbers of particles to be transfered across the
 187             // processor patches
 188             labelList nsTransPs(transferList.size());
 189
 190             forAll(transferList, i)
 191             {
 192                 nsTransPs[i] = transferList[i].size();
 193             }
 194
 195             // List of the numbers of particles to be transfered across the
 196             // processor patches for all the processors
 197             labelListList allNTrans(Pstream::nProcs());
 198             allNTrans[Pstream::myProcNo()] = nsTransPs;
 199             combineReduce(allNTrans, combineNsTransPs());
 200
 201             transfered = false;
 202
 203             forAll(allNTrans, i)
 204             {
 205                 forAll(allNTrans[i], j)
 206                 {
 207                     if (allNTrans[i][j])
 208                     {
 209                         transfered = true;
 210                         break;
 211                     }
 212                 }
 213             }
 214
 215             if (!transfered)
 216             {
 217                 break;
 218             }
 219
 220             forAll(transferList, i)
 221             {
 222                 if (transferList[i].size())
 223                 {
 224                     OPstream particleStream
 225                     (
 226                         Pstream::blocking,
 227                         refCast<const processorPolyPatch>
 228                         (
 229                             pMesh().boundaryMesh()[processorPatches[i]]
 230                         ).neighbProcNo()
 231                     );
 232
 233                     particleStream << transferList[i];
 234                 }
 235             }
 236
 237             forAll(processorPatches, i)
 238             {
 239                 label patchi = processorPatches[i];
 240
 241                 const processorPolyPatch& procPatch =
 242                     refCast<const processorPolyPatch>
 243                     (pMesh().boundaryMesh()[patchi]);
 244
 245                 label neighbProci =
 246                     procPatch.neighbProcNo() - Pstream::masterNo();
 247
 248                 label neighbProcPatchi = processorPatchNeighbours[patchi];
 249
 250                 label nRecPs = allNTrans[neighbProci][neighbProcPatchi];
 251
 252                 if (nRecPs)
 253                 {
 254                     IPstream particleStream
 255                     (
 256                         Pstream::blocking,
 257                         procPatch.neighbProcNo()
 258                     );
 259                     IDLList<ParticleType> newParticles
 260                     (
 261                         particleStream,
 262                         typename ParticleType::iNew(*this)
 263                     );
 264
 265                     forAllIter
 266                     (
 267                         typename Cloud<ParticleType>,
 268                         newParticles,
 269                         newpIter
 270                     )
 271                     {
 272                         ParticleType& newp = newpIter();
 273                         newp.correctAfterParallelTransfer(patchi, td);
 274                         addParticle(newParticles.remove(&newp));
 275                     }
 276                 }
 277             }
 278         }
 279         else
 280         {
 281             transfered = false;
 282         }
 283     }
 284 }
 285
 286
 287 template<class ParticleType>
 288 void Foam::Cloud<ParticleType>::autoMap(const mapPolyMesh& mapper)
 289 {
 290     if (cloud::debug)
 291     {
 292         Info<< "Cloud<ParticleType>::autoMap(const morphFieldMapper& map) "
 293                "for lagrangian cloud " << cloud::name() << endl;
 294     }
 295
 296     const labelList& reverseCellMap = mapper.reverseCellMap();
 297     const labelList& reverseFaceMap = mapper.reverseFaceMap();
 298
 299     forAllIter(typename Cloud<ParticleType>, *this, pIter)
 300     {
 301         if (reverseCellMap[pIter().celli_] >= 0)
 302         {
 303             pIter().celli_ = reverseCellMap[pIter().celli_];
 304
 305             if (pIter().facei_ >= 0 && reverseFaceMap[pIter().facei_] >= 0)
 306             {
 307                 pIter().facei_ = reverseFaceMap[pIter().facei_];
 308             }
 309             else
 310             {
 311                 pIter().facei_ = -1;
 312             }
 313         }
 314         else
 315         {
 316             label trackStartCell = mapper.mergedCell(pIter().celli_);
 317
 318             if (trackStartCell < 0)
 319             {
 320                 trackStartCell = 0;
 321             }
 322
 323             vector p = pIter().position();
 324             const_cast<vector&>(pIter().position()) =
 325                 polyMesh_.cellCentres()[trackStartCell];
 326             pIter().stepFraction() = 0;
 327             pIter().track(p);
 328         }
 329     }
 330 }
 331
 332
 333 template<class ParticleType>
 334 void Foam::Cloud<ParticleType>::writePositions() const
 335 {
 336     OFstream pObj
 337     (
 338         this->db().time().path()/this->name() + "_positions.obj"
 339     );
 340
 341     forAllConstIter(typename Cloud<ParticleType>, *this, pIter)
 342     {
 343         const ParticleType& p = pIter();
 344         pObj<< "v " << p.position().x() << " " << p.position().y() << " "
 345             << p.position().z() << nl;
 346     }
 347
 348     pObj.flush();
 349 }
 350
 351
 352 // * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
 353
 354 #include "CloudIO.C"
 355
 356 // ************************************************************************* //