1270fbac96117884cd50d2c8946a91a0928bb885
[OpenFOAM-1.5.x.git] / src / lagrangian / basic / Particle / Particle.C
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.
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 "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"
35
36 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
37
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_];
47
48     labelList faceList(0);
49     forAll(faces, i)
50     {
51         label facei = faces[i];
52         scalar lam = lambda(C, position, facei);
53
54         if ((lam > 0) && (lam < 1.0))
55         {
56             label n = faceList.size();
57             faceList.setSize(n+1);
58             faceList[n] = facei;
59         }
60     }
61
62     return faceList;
63 }
64
65
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];
77
78     labelList faceList(0);
79     forAll(faces, i)
80     {
81         label facei = faces[i];
82         scalar lam = lambda(C, position, facei, stepFraction);
83
84         if ((lam > 0) && (lam < 1.0))
85         {
86             label n = faceList.size();
87             faceList.setSize(n+1);
88             faceList[n] = facei;
89         }
90     }
91
92     return faceList;
93 }
94
95
96 template<class ParticleType>
97 template<class TrackData>
98 void Foam::Particle<ParticleType>::prepareForParallelTransfer
99 (
100     const label patchi,
101     TrackData& td
102 )
103 {
104     // Convert the face index to be local to the processor patch
105     facei_ = patchFace(patchi, facei_);
106 }
107
108
109 template<class ParticleType>
110 template<class TrackData>
111 void Foam::Particle<ParticleType>::correctAfterParallelTransfer
112 (
113     const label patchi,
114     TrackData& td
115 )
116 {
117     const processorPolyPatch& ppp =
118         refCast<const processorPolyPatch>
119         (cloud_.pMesh().boundaryMesh()[patchi]);
120
121     celli_ = ppp.faceCells()[facei_];
122
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     }
157
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     }
168 }
169
170
171 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
172
173 template<class ParticleType>
174 Foam::Particle<ParticleType>::Particle
175 (
176     const Cloud<ParticleType>& cloud,
177     const vector& position,
178     const label celli
179 )
180 :
181     cloud_(cloud),
182     position_(position),
183     celli_(celli),
184     facei_(-1),
185     stepFraction_(0.0)
186 {}
187
188
189 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
190
191 template<class ParticleType>
192 template<class TrackData>
193 Foam::label Foam::Particle<ParticleType>::track
194 (
195     const vector& endPosition,
196     TrackData& td
197 )
198 {
199     facei_ = -1;
200
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     }
206
207     return facei_;
208 }
209
210
211
212 template<class ParticleType>
213 Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
214 {
215     int dummyTd;
216     return track(endPosition, dummyTd);
217 }
218
219 template<class ParticleType>
220 template<class TrackData>
221 Foam::scalar Foam::Particle<ParticleType>::trackToFace
222 (
223     const vector& endPosition,
224     TrackData& td
225 )
226 {
227     const polyMesh& mesh = cloud_.polyMesh_;
228
229     labelList faces = findFaces(endPosition);
230
231     facei_ = -1;
232     scalar trackFraction = 0.0;
233
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;
242
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_);
259
260                 if (lam < lambdaMin)
261                 {
262                     lambdaMin = lam;
263                     facei_ = faces[i];
264                 }
265             }
266         }
267
268         bool internalFace = cloud_.internalFace(facei_);
269
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         }
296
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);
320
321             // Soft-sphere algorithm ignores the boundary
322             if (p.softImpact())
323             {
324                 trackFraction = 1.0;
325                 position_ = endPosition;
326             }
327
328             label patchi = patch(facei_);
329             const polyPatch& patch = mesh.boundaryMesh()[patchi];
330
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     }
384
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     }
395
396     return trackFraction;
397 }
398
399 template<class ParticleType>
400 Foam::scalar Foam::Particle<ParticleType>::trackToFace
401 (
402     const vector& endPosition
403 )
404 {
405     int dummyTd;
406     return trackToFace(endPosition, dummyTd);
407 }
408
409 template<class ParticleType>
410 void Foam::Particle<ParticleType>::transformPosition(const tensor& T)
411 {
412     position_ = transform(T, position_);
413 }
414
415
416 template<class ParticleType>
417 void Foam::Particle<ParticleType>::transformProperties(const tensor&)
418 {}
419
420
421 template<class ParticleType>
422 void Foam::Particle<ParticleType>::transformProperties(const vector&)
423 {}
424
425
426 template<class ParticleType>
427 template<class TrackData>
428 void Foam::Particle<ParticleType>::hitWedgePatch
429 (
430     const wedgePolyPatch& wpp,
431     TrackData&
432 )
433 {
434     vector nf = wpp.faceAreas()[wpp.whichFace(facei_)];
435     nf /= mag(nf);
436
437     static_cast<ParticleType&>(*this).transformProperties(I - 2.0*nf*nf);
438 }
439
440
441 template<class ParticleType>
442 template<class TrackData>
443 void Foam::Particle<ParticleType>::hitSymmetryPatch
444 (
445     const symmetryPolyPatch& spp,
446     TrackData&
447 )
448 {
449     vector nf = spp.faceAreas()[spp.whichFace(facei_)];
450     nf /= mag(nf);
451
452     static_cast<ParticleType&>(*this).transformProperties(I - 2.0*nf*nf);
453 }
454
455
456 template<class ParticleType>
457 template<class TrackData>
458 void Foam::Particle<ParticleType>::hitCyclicPatch
459 (
460     const cyclicPolyPatch& cpp,
461     TrackData&
462 )
463 {
464     label patchFacei_ = cpp.whichFace(facei_);
465
466     facei_ = cpp.transformGlobalFace(facei_);
467
468     celli_ = cloud_.polyMesh_.faceOwner()[facei_];
469
470     if (!cpp.parallel())
471     {
472         const tensor& T = cpp.transformT(patchFacei_);
473
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     }
485 }
486
487
488 template<class ParticleType>
489 template<class TrackData>
490 void Foam::Particle<ParticleType>::hitProcessorPatch
491 (
492     const processorPolyPatch& spp,
493     TrackData& td
494 )
495 {}
496
497
498 template<class ParticleType>
499 template<class TrackData>
500 void Foam::Particle<ParticleType>::hitWallPatch
501 (
502     const wallPolyPatch& spp,
503     TrackData&
504 )
505 {}
506
507
508 template<class ParticleType>
509 template<class TrackData>
510 void Foam::Particle<ParticleType>::hitPatch
511 (
512     const polyPatch& spp,
513     TrackData&
514 )
515 {}
516
517
518 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
519
520 #include "ParticleIO.C"
521
522 // ************************************************************************* //