Merge branch 'master' of ssh://git.code.sf.net/p/foam-extend/foam-extend-3.2
[foam-extend-3.2.git] / src / lagrangian / basic / CloudTemplate / CloudTemplate.C
blobaaac8a3b6f747cab9a5c84b0a99825b79c37f660
1 /*---------------------------------------------------------------------------*\
2   =========                 |
3   \\      /  F ield         | foam-extend: Open Source CFD
4    \\    /   O peration     | Version:     3.2
5     \\  /    A nd           | Web:         http://www.foam-extend.org
6      \\/     M anipulation  | For copyright notice see file Copyright
7 -------------------------------------------------------------------------------
8 License
9     This file is part of foam-extend.
11     foam-extend 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 3 of the License, or (at your
14     option) any later version.
16     foam-extend is distributed in the hope that it will be useful, but
17     WITHOUT ANY WARRANTY; without even the implied warranty of
18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
19     General Public License for more details.
21     You should have received a copy of the GNU General Public License
22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
24 \*---------------------------------------------------------------------------*/
26 #include "CloudTemplate.H"
27 #include "processorPolyPatch.H"
28 #include "globalMeshData.H"
29 #include "PstreamCombineReduceOps.H"
30 #include "mapPolyMesh.H"
31 #include "foamTime.H"
32 #include "OFstream.H"
34 #include "profiling.H"
36 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
38 template<class ParticleType>
39 Foam::Cloud<ParticleType>::Cloud
41     const polyMesh& pMesh,
42     const IDLList<ParticleType>& particles
45     cloud(pMesh),
46     IDLList<ParticleType>(),
47     polyMesh_(pMesh),
48     particleCount_(0)
50     IDLList<ParticleType>::operator=(particles);
54 template<class ParticleType>
55 Foam::Cloud<ParticleType>::Cloud
57     const polyMesh& pMesh,
58     const word& cloudName,
59     const IDLList<ParticleType>& particles
62     cloud(pMesh, cloudName),
63     IDLList<ParticleType>(),
64     polyMesh_(pMesh),
65     particleCount_(0)
67     IDLList<ParticleType>::operator=(particles);
71 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
73 template<class ParticleType>
74 Foam::label Foam::Cloud<ParticleType>::getNewParticleID() const
76     label id = particleCount_++;
78     if (id == labelMax)
79     {
80         WarningIn("Cloud<ParticleType>::getNewParticleID() const")
81             << "Particle counter has overflowed. This might cause problems"
82             << " when reconstructing particle tracks." << endl;
83     }
84     return id;
88 template<class ParticleType>
89 void Foam::Cloud<ParticleType>::addParticle(ParticleType* pPtr)
91     this->append(pPtr);
95 template<class ParticleType>
96 void Foam::Cloud<ParticleType>::deleteParticle(ParticleType& p)
98     delete(this->remove(&p));
102 namespace Foam
105 class combineNsTransPs
108 public:
110     void operator()(labelListList& x, const labelListList& y) const
111     {
112         forAll(y, i)
113         {
114             if (y[i].size())
115             {
116                 x[i] = y[i];
117             }
118         }
119     }
122 } // End namespace Foam
125 template<class ParticleType>
126 template<class TrackingData>
127 void Foam::Cloud<ParticleType>::move(TrackingData& td)
129     addProfile2(moveProfile,"Cloud<ParticleType>::move_"+this->name());
131     const globalMeshData& pData = polyMesh_.globalData();
132     const labelList& processorPatches = pData.processorPatches();
133     const labelList& processorPatchIndices = pData.processorPatchIndices();
134     const labelList& processorPatchNeighbours =
135         pData.processorPatchNeighbours();
137     // Initialise the setpFraction moved for the particles
138     forAllIter(typename Cloud<ParticleType>, *this, pIter)
139     {
140         pIter().stepFraction() = 0;
141     }
143     // Assume there will be particles to transfer
144     bool transfered = true;
146     // While there are particles to transfer
147     while (transfered)
148     {
149         // List of lists of particles to be transfered for all the processor
150         // patches
151         List<IDLList<ParticleType> > transferList(processorPatches.size());
153         // Loop over all particles
154         forAllIter(typename Cloud<ParticleType>, *this, pIter)
155         {
156             ParticleType& p = pIter();
158             // Move the particle
159             bool keepParticle = p.move(td);
161             // If the particle is to be kept
162             // (i.e. it hasn't passed through an inlet or outlet)
163             if (keepParticle)
164             {
165                 // If we are running in parallel and the particle is on a
166                 // boundary face
167                 if (Pstream::parRun() && p.facei_ >= pMesh().nInternalFaces())
168                 {
169                     label patchi = pMesh().boundaryMesh().whichPatch(p.facei_);
170                     label n = processorPatchIndices[patchi];
172                     // ... and the face is on a processor patch
173                     // prepare it for transfer
174                     if (n != -1)
175                     {
176                         p.prepareForParallelTransfer(patchi, td);
177                         transferList[n].append(this->remove(&p));
178                     }
179                 }
180             }
181             else
182             {
183                 deleteParticle(p);
184             }
185         }
187         if (Pstream::parRun())
188         {
189             // List of the numbers of particles to be transfered across the
190             // processor patches
191             labelList nsTransPs(transferList.size());
193             forAll(transferList, i)
194             {
195                 nsTransPs[i] = transferList[i].size();
196             }
198             // List of the numbers of particles to be transfered across the
199             // processor patches for all the processors
200             labelListList allNTrans(Pstream::nProcs());
201             allNTrans[Pstream::myProcNo()] = nsTransPs;
202             combineReduce(allNTrans, combineNsTransPs());
204             transfered = false;
206             forAll(allNTrans, i)
207             {
208                 forAll(allNTrans[i], j)
209                 {
210                     if (allNTrans[i][j])
211                     {
212                         transfered = true;
213                         break;
214                     }
215                 }
216             }
218             if (!transfered)
219             {
220                 break;
221             }
223             forAll(transferList, i)
224             {
225                 if (transferList[i].size())
226                 {
227                     OPstream particleStream
228                     (
229                         Pstream::blocking,
230                         refCast<const processorPolyPatch>
231                         (
232                             pMesh().boundaryMesh()[processorPatches[i]]
233                         ).neighbProcNo()
234                     );
236                     particleStream << transferList[i];
237                 }
238             }
240             forAll(processorPatches, i)
241             {
242                 label patchi = processorPatches[i];
244                 const processorPolyPatch& procPatch =
245                     refCast<const processorPolyPatch>
246                     (pMesh().boundaryMesh()[patchi]);
248                 label neighbProci =
249                     procPatch.neighbProcNo() - Pstream::masterNo();
251                 label neighbProcPatchi = processorPatchNeighbours[patchi];
253                 label nRecPs = allNTrans[neighbProci][neighbProcPatchi];
255                 if (nRecPs)
256                 {
257                     IPstream particleStream
258                     (
259                         Pstream::blocking,
260                         procPatch.neighbProcNo()
261                     );
262                     IDLList<ParticleType> newParticles
263                     (
264                         particleStream,
265                         typename ParticleType::iNew(*this)
266                     );
268                     forAllIter
269                     (
270                         typename Cloud<ParticleType>,
271                         newParticles,
272                         newpIter
273                     )
274                     {
275                         ParticleType& newp = newpIter();
276                         newp.correctAfterParallelTransfer(patchi, td);
277                         addParticle(newParticles.remove(&newp));
278                     }
279                 }
280             }
281         }
282         else
283         {
284             transfered = false;
285         }
286     }
290 template<class ParticleType>
291 void Foam::Cloud<ParticleType>::autoMap(const mapPolyMesh& mapper)
293     if (cloud::debug)
294     {
295         Info<< "Cloud<ParticleType>::autoMap(const morphFieldMapper& map) "
296                "for lagrangian cloud " << cloud::name() << endl;
297     }
299     const labelList& reverseCellMap = mapper.reverseCellMap();
300     const labelList& reverseFaceMap = mapper.reverseFaceMap();
302     forAllIter(typename Cloud<ParticleType>, *this, pIter)
303     {
304         if (reverseCellMap[pIter().celli_] >= 0)
305         {
306             pIter().celli_ = reverseCellMap[pIter().celli_];
308             if (pIter().facei_ >= 0 && reverseFaceMap[pIter().facei_] >= 0)
309             {
310                 pIter().facei_ = reverseFaceMap[pIter().facei_];
311             }
312             else
313             {
314                 pIter().facei_ = -1;
315             }
316         }
317         else
318         {
319             label trackStartCell = mapper.mergedCell(pIter().celli_);
321 //          Tommaso, bug fix 27/02/2008
322 //          when the topology of the mesh is changed
323 //          it is necessary to update the cell in which the parcel is,
324 //          in particular when cells are removed
326             // HJ, merge 1.5.x  20/Oct/2008
327             if (trackStartCell < 0)
328             {
329                 pIter().celli_ = polyMesh_.findCell(pIter().position());
330             }
332             vector p = pIter().position();
333             const_cast<vector&>(pIter().position()) =
334                 polyMesh_.cellCentres()[trackStartCell];
335             pIter().stepFraction() = 0;
336             pIter().track(p);
337         }
338     }
342 template<class ParticleType>
343 void Foam::Cloud<ParticleType>::writePositions() const
345     OFstream pObj
346     (
347         this->db().time().path()/this->name() + "_positions.obj"
348     );
350     forAllConstIter(typename Cloud<ParticleType>, *this, pIter)
351     {
352         const ParticleType& p = pIter();
353         pObj<< "v " << p.position().x() << " " << p.position().y() << " "
354             << p.position().z() << nl;
355     }
357     pObj.flush();
361 // * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
363 #include "CloudTemplateIO.C"
365 // ************************************************************************* //