1270fbac96117884cd50d2c8946a91a0928bb885
[OpenFOAM-1.5.x.git] / src / lagrangian / basic / Particle / Particle.C
blob1270fbac96117884cd50d2c8946a91a0928bb885
1 /*---------------------------------------------------------------------------*\
2   =========                 |
3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
4    \\    /   O peration     |
5     \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
6      \\/     M anipulation  |
7 -------------------------------------------------------------------------------
8 License
9     This file is part of OpenFOAM.
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.
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.
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
25 \*---------------------------------------------------------------------------*/
27 #include "Particle.H"
28 #include "Cloud.H"
29 #include "wedgePolyPatch.H"
30 #include "symmetryPolyPatch.H"
31 #include "cyclicPolyPatch.H"
32 #include "processorPolyPatch.H"
33 #include "wallPolyPatch.H"
34 #include "transform.H"
36 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
38 template<class ParticleType>
39 Foam::labelList Foam::Particle<ParticleType>::findFaces
40 (
41     const vector& position
42 ) const
43 {
44     const polyMesh& mesh = cloud_.polyMesh_;
45     const labelList& faces = mesh.cells()[celli_];
46     const vector& C = mesh.cellCentres()[celli_];
48     labelList faceList(0);
49     forAll(faces, i)
50     {
51         label facei = faces[i];
52         scalar lam = lambda(C, position, facei);
54         if ((lam > 0) && (lam < 1.0))
55         {
56             label n = faceList.size();
57             faceList.setSize(n+1);
58             faceList[n] = facei;
59         }
60     }
62     return faceList;
63 }
66 template<class ParticleType>
67 Foam::labelList Foam::Particle<ParticleType>::findFaces
68 (
69     const vector& position,
70     const label celli,
71     const scalar stepFraction
72 ) const
73 {
74     const polyMesh& mesh = cloud_.polyMesh_;
75     const labelList& faces = mesh.cells()[celli];
76     const vector& C = mesh.cellCentres()[celli];
78     labelList faceList(0);
79     forAll(faces, i)
80     {
81         label facei = faces[i];
82         scalar lam = lambda(C, position, facei, stepFraction);
84         if ((lam > 0) && (lam < 1.0))
85         {
86             label n = faceList.size();
87             faceList.setSize(n+1);
88             faceList[n] = facei;
89         }
90     }
92     return faceList;
93 }
96 template<class ParticleType>
97 template<class TrackData>
98 void Foam::Particle<ParticleType>::prepareForParallelTransfer
99 (
100     const label patchi,
101     TrackData& td
104     // Convert the face index to be local to the processor patch
105     facei_ = patchFace(patchi, facei_);
109 template<class ParticleType>
110 template<class TrackData>
111 void Foam::Particle<ParticleType>::correctAfterParallelTransfer
113     const label patchi,
114     TrackData& td
117     const processorPolyPatch& ppp =
118         refCast<const processorPolyPatch>
119         (cloud_.pMesh().boundaryMesh()[patchi]);
121     celli_ = ppp.faceCells()[facei_];
123     if (!ppp.parallel())
124     {
125         if (ppp.forwardT().size() == 1)
126         {
127             const tensor& T = ppp.forwardT()[0];
128             transformPosition(T);
129             static_cast<ParticleType&>(*this).transformProperties(T);
130         }
131         else
132         {
133             const tensor& T = ppp.forwardT()[facei_];
134             transformPosition(T);
135             static_cast<ParticleType&>(*this).transformProperties(T);
136         }
137     }
138     else if (ppp.separated())
139     {
140         if (ppp.separation().size() == 1)
141         {
142             position_ -= ppp.separation()[0];
143             static_cast<ParticleType&>(*this).transformProperties
144             (
145                 -ppp.separation()[0]
146             );
147         }
148         else
149         {
150             position_ -= ppp.separation()[facei_];
151             static_cast<ParticleType&>(*this).transformProperties
152             (
153                 -ppp.separation()[facei_]
154             );
155         }
156     }
158     // Reset the face index for the next tracking operation
159     if (stepFraction_ > (1.0 - SMALL))
160     {
161         stepFraction_ = 1.0;
162         facei_ = -1;
163     }
164     else
165     {
166         facei_ += ppp.start();
167     }
171 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
173 template<class ParticleType>
174 Foam::Particle<ParticleType>::Particle
176     const Cloud<ParticleType>& cloud,
177     const vector& position,
178     const label celli
181     cloud_(cloud),
182     position_(position),
183     celli_(celli),
184     facei_(-1),
185     stepFraction_(0.0)
186 {}
189 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
191 template<class ParticleType>
192 template<class TrackData>
193 Foam::label Foam::Particle<ParticleType>::track
195     const vector& endPosition,
196     TrackData& td
199     facei_ = -1;
201     // Tracks to endPosition or stop on boundary
202     while(!onBoundary() && stepFraction_ < 1.0 - SMALL)
203     {
204         stepFraction_ += trackToFace(endPosition, td)*(1.0 - stepFraction_);
205     }
207     return facei_;
212 template<class ParticleType>
213 Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
215     int dummyTd;
216     return track(endPosition, dummyTd);
219 template<class ParticleType>
220 template<class TrackData>
221 Foam::scalar Foam::Particle<ParticleType>::trackToFace
223     const vector& endPosition,
224     TrackData& td
227     const polyMesh& mesh = cloud_.polyMesh_;
229     labelList faces = findFaces(endPosition);
231     facei_ = -1;
232     scalar trackFraction = 0.0;
234     if (faces.size() == 0) // inside cell
235     {
236         trackFraction = 1.0;
237         position_ = endPosition;
238     }
239     else // hit face
240     {
241         scalar lambdaMin = GREAT;
243         if (faces.size() == 1)
244         {
245             lambdaMin = lambda(position_, endPosition, faces[0], stepFraction_);
246             facei_ = faces[0];
247         }
248         else
249         {
250             // If the particle has to cross more than one cell to reach the
251             // endPosition, we check which way to go.
252             // If one of the faces is a boundary face and the particle is
253             // outside, we choose the boundary face.
254             // The particle is outside if one of the lambda's is > 1 or < 0
255             forAll(faces, i)
256             {
257                 scalar lam =
258                     lambda(position_, endPosition, faces[i], stepFraction_);
260                 if (lam < lambdaMin)
261                 {
262                     lambdaMin = lam;
263                     facei_ = faces[i];
264                 }
265             }
266         }
268         bool internalFace = cloud_.internalFace(facei_);
270         // For warped faces the particle can be 'outside' the cell.
271         // This will yield a lambda larger than 1, or smaller than 0
272         // For values < 0, the particle travels away from the cell
273         // and we don't move the particle, only change cell.
274         // For values larger than 1, we move the particle to endPosition only.
275         if (lambdaMin > 0.0)
276         {
277             if (lambdaMin <= 1.0)
278             {
279                 trackFraction = lambdaMin;
280                 position_ += trackFraction*(endPosition - position_);
281             }
282             else
283             {
284                 trackFraction = 1.0;
285                 position_ = endPosition;
286             }
287         }
288         else if (static_cast<ParticleType&>(*this).softImpact())
289         {
290             // Soft-sphere particles can travel outside the domain
291             // but we don't use lambda since this the particle
292             // is going away from face
293             trackFraction = 1.0;
294             position_ = endPosition;
295         }
297         // change cell
298         if (internalFace) // Internal face
299         {
300             if (celli_ == cloud_.owner_[facei_])
301             {
302                 celli_ = cloud_.neighbour_[facei_];
303             }
304             else if (celli_ == cloud_.neighbour_[facei_])
305             {
306                 celli_ = cloud_.owner_[facei_];
307             }
308             else
309             {
310                 FatalErrorIn
311                 (
312                     "Particle::trackToFace(const vector&, TrackData&)"
313                 )<< "addressing failure" << nl
314                  << abort(FatalError);
315             }
316         }
317         else
318         {
319             ParticleType& p = static_cast<ParticleType&>(*this);
321             // Soft-sphere algorithm ignores the boundary
322             if (p.softImpact())
323             {
324                 trackFraction = 1.0;
325                 position_ = endPosition;
326             }
328             label patchi = patch(facei_);
329             const polyPatch& patch = mesh.boundaryMesh()[patchi];
331             if (isA<wedgePolyPatch>(patch))
332             {
333                 p.hitWedgePatch
334                 (
335                     static_cast<const wedgePolyPatch&>(patch), td
336                 );
337             }
338             else if (isA<symmetryPolyPatch>(patch))
339             {
340                 p.hitSymmetryPatch
341                 (
342                     static_cast<const symmetryPolyPatch&>(patch), td
343                 );
344             }
345             else if (isA<cyclicPolyPatch>(patch))
346             {
347                 p.hitCyclicPatch
348                 (
349                     static_cast<const cyclicPolyPatch&>(patch), td
350                 );
351             }
352             else if (isA<processorPolyPatch>(patch))
353             {
354                 p.hitProcessorPatch
355                 (
356                     static_cast<const processorPolyPatch&>(patch), td
357                 );
358             }
359             else if (isA<wallPolyPatch>(patch))
360             {
361                 p.hitWallPatch
362                 (
363                     static_cast<const wallPolyPatch&>(patch), td
364                 );
365             }
366             else if (isA<polyPatch>(patch))
367             {
368                 p.hitPatch
369                 (
370                     static_cast<const polyPatch&>(patch), td
371                 );
372             }
373             else
374             {
375                 FatalErrorIn
376                 (
377                     "Particle::trackToFace"
378                     "(const vector& endPosition, scalar& trackFraction)"
379                 )<< "patch type " << patch.type() << " not suported" << nl
380                  << abort(FatalError);
381             }
382         }
383     }
385     // If the trackFraction = 0 something went wrong.
386     // Either the particle is flipping back and forth across a face perhaps
387     // due to velocity interpolation errors or it is in a "hole" in the mesh
388     // caused by face warpage.
389     // In both cases resolve the positional ambiguity by moving the particle
390     // slightly towards the cell-centre.
391     if (trackFraction < SMALL)
392     {
393         position_ += 1.0e-6*(mesh.cellCentres()[celli_] - position_);
394     }
396     return trackFraction;
399 template<class ParticleType>
400 Foam::scalar Foam::Particle<ParticleType>::trackToFace
402     const vector& endPosition
405     int dummyTd;
406     return trackToFace(endPosition, dummyTd);
409 template<class ParticleType>
410 void Foam::Particle<ParticleType>::transformPosition(const tensor& T)
412     position_ = transform(T, position_);
416 template<class ParticleType>
417 void Foam::Particle<ParticleType>::transformProperties(const tensor&)
418 {}
421 template<class ParticleType>
422 void Foam::Particle<ParticleType>::transformProperties(const vector&)
423 {}
426 template<class ParticleType>
427 template<class TrackData>
428 void Foam::Particle<ParticleType>::hitWedgePatch
430     const wedgePolyPatch& wpp,
431     TrackData&
434     vector nf = wpp.faceAreas()[wpp.whichFace(facei_)];
435     nf /= mag(nf);
437     static_cast<ParticleType&>(*this).transformProperties(I - 2.0*nf*nf);
441 template<class ParticleType>
442 template<class TrackData>
443 void Foam::Particle<ParticleType>::hitSymmetryPatch
445     const symmetryPolyPatch& spp,
446     TrackData&
449     vector nf = spp.faceAreas()[spp.whichFace(facei_)];
450     nf /= mag(nf);
452     static_cast<ParticleType&>(*this).transformProperties(I - 2.0*nf*nf);
456 template<class ParticleType>
457 template<class TrackData>
458 void Foam::Particle<ParticleType>::hitCyclicPatch
460     const cyclicPolyPatch& cpp,
461     TrackData&
464     label patchFacei_ = cpp.whichFace(facei_);
466     facei_ = cpp.transformGlobalFace(facei_);
468     celli_ = cloud_.polyMesh_.faceOwner()[facei_];
470     if (!cpp.parallel())
471     {
472         const tensor& T = cpp.transformT(patchFacei_);
474         transformPosition(T);
475         static_cast<ParticleType&>(*this).transformProperties(T);
476     }
477     else if (cpp.separated())
478     {
479         position_ += cpp.separation(patchFacei_);
480         static_cast<ParticleType&>(*this).transformProperties
481         (
482             cpp.separation(patchFacei_)
483         );
484     }
488 template<class ParticleType>
489 template<class TrackData>
490 void Foam::Particle<ParticleType>::hitProcessorPatch
492     const processorPolyPatch& spp,
493     TrackData& td
495 {}
498 template<class ParticleType>
499 template<class TrackData>
500 void Foam::Particle<ParticleType>::hitWallPatch
502     const wallPolyPatch& spp,
503     TrackData&
505 {}
508 template<class ParticleType>
509 template<class TrackData>
510 void Foam::Particle<ParticleType>::hitPatch
512     const polyPatch& spp,
513     TrackData&
515 {}
518 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
520 #include "ParticleIO.C"
522 // ************************************************************************* //