FreeFOAM The Cross-Platform CFD Toolkit
Particle.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 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 "Particle.H"
27 #include <lagrangian/Cloud.H>
32 #include <OpenFOAM/wallPolyPatch.H>
33 #include <OpenFOAM/transform.H>
34 
35 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
36 
37 template<class ParticleType>
39 (
40  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  faceList.clear();
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  faceList.append(facei);
57  }
58  }
59 }
60 
61 
62 template<class ParticleType>
64 (
65  const vector& position,
66  const label celli,
67  const scalar stepFraction,
68  DynamicList<label>& faceList
69 ) const
70 {
71  const polyMesh& mesh = cloud_.pMesh();
72  const labelList& faces = mesh.cells()[celli];
73  const vector& C = mesh.cellCentres()[celli];
74 
75  faceList.clear();
76  forAll(faces, i)
77  {
78  label facei = faces[i];
79  scalar lam = lambda(C, position, facei, stepFraction);
80 
81  if ((lam > 0) && (lam < 1.0))
82  {
83  faceList.append(facei);
84  }
85  }
86 }
87 
88 
89 template<class ParticleType>
90 template<class TrackData>
92 (
93  const label patchi,
94  TrackData& td
95 )
96 {
97  // Convert the face index to be local to the processor patch
98  facei_ = patchFace(patchi, facei_);
99 }
100 
101 
102 template<class ParticleType>
103 template<class TrackData>
105 (
106  const label patchi,
107  TrackData& td
108 )
109 {
110  const processorPolyPatch& ppp =
111  refCast<const processorPolyPatch>
112  (cloud_.pMesh().boundaryMesh()[patchi]);
113 
114  celli_ = ppp.faceCells()[facei_];
115 
116  if (!ppp.parallel())
117  {
118  if (ppp.forwardT().size() == 1)
119  {
120  const tensor& T = ppp.forwardT()[0];
121  transformPosition(T);
122  static_cast<ParticleType&>(*this).transformProperties(T);
123  }
124  else
125  {
126  const tensor& T = ppp.forwardT()[facei_];
127  transformPosition(T);
128  static_cast<ParticleType&>(*this).transformProperties(T);
129  }
130  }
131  else if (ppp.separated())
132  {
133  if (ppp.separation().size() == 1)
134  {
135  position_ -= ppp.separation()[0];
136  static_cast<ParticleType&>(*this).transformProperties
137  (
138  -ppp.separation()[0]
139  );
140  }
141  else
142  {
143  position_ -= ppp.separation()[facei_];
144  static_cast<ParticleType&>(*this).transformProperties
145  (
146  -ppp.separation()[facei_]
147  );
148  }
149  }
150 
151  // Reset the face index for the next tracking operation
152  if (stepFraction_ > (1.0 - SMALL))
153  {
154  stepFraction_ = 1.0;
155  facei_ = -1;
156  }
157  else
158  {
159  facei_ += ppp.start();
160  }
161 }
162 
163 
164 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
165 
166 template<class ParticleType>
168 (
169  const Cloud<ParticleType>& cloud,
170  const vector& position,
171  const label celli
172 )
173 :
174  cloud_(cloud),
175  position_(position),
176  celli_(celli),
177  facei_(-1),
178  stepFraction_(0.0),
179  origProc_(Pstream::myProcNo()),
180  origId_(cloud_.getNewParticleID())
181 {}
182 
183 
184 template<class ParticleType>
186 :
187  cloud_(p.cloud_),
188  position_(p.position_),
189  celli_(p.celli_),
190  facei_(p.facei_),
191  stepFraction_(p.stepFraction_),
192  origProc_(p.origProc_),
193  origId_(p.origId_)
194 {}
195 
196 
197 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
198 
199 template<class ParticleType>
200 template<class TrackData>
202 (
203  const vector& endPosition,
204  TrackData& td
205 )
206 {
207  facei_ = -1;
208 
209  // Tracks to endPosition or stop on boundary
210  while(!onBoundary() && stepFraction_ < 1.0 - SMALL)
211  {
212  stepFraction_ += trackToFace(endPosition, td)*(1.0 - stepFraction_);
213  }
214 
215  return facei_;
216 }
217 
218 
219 
220 template<class ParticleType>
221 Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
222 {
223  int dummyTd;
224  return track(endPosition, dummyTd);
225 }
226 
227 template<class ParticleType>
228 template<class TrackData>
230 (
231  const vector& endPosition,
232  TrackData& td
233 )
234 {
235  const polyMesh& mesh = cloud_.polyMesh_;
236 
237  DynamicList<label>& faces = cloud_.labels_;
238  findFaces(endPosition, faces);
239 
240  facei_ = -1;
241  scalar trackFraction = 0.0;
242 
243  if (faces.empty()) // inside cell
244  {
245  trackFraction = 1.0;
246  position_ = endPosition;
247  }
248  else // hit face
249  {
250  scalar lambdaMin = GREAT;
251 
252  if (faces.size() == 1)
253  {
254  lambdaMin = lambda(position_, endPosition, faces[0], stepFraction_);
255  facei_ = faces[0];
256  }
257  else
258  {
259  // If the particle has to cross more than one cell to reach the
260  // endPosition, we check which way to go.
261  // If one of the faces is a boundary face and the particle is
262  // outside, we choose the boundary face.
263  // The particle is outside if one of the lambda's is > 1 or < 0
264  forAll(faces, i)
265  {
266  scalar lam =
267  lambda(position_, endPosition, faces[i], stepFraction_);
268 
269  if (lam < lambdaMin)
270  {
271  lambdaMin = lam;
272  facei_ = faces[i];
273  }
274  }
275  }
276 
277  bool internalFace = cloud_.internalFace(facei_);
278 
279  // For warped faces the particle can be 'outside' the cell.
280  // This will yield a lambda larger than 1, or smaller than 0
281  // For values < 0, the particle travels away from the cell
282  // and we don't move the particle, only change cell.
283  // For values larger than 1, we move the particle to endPosition only.
284  if (lambdaMin > 0.0)
285  {
286  if (lambdaMin <= 1.0)
287  {
288  trackFraction = lambdaMin;
289  position_ += trackFraction*(endPosition - position_);
290  }
291  else
292  {
293  trackFraction = 1.0;
294  position_ = endPosition;
295  }
296  }
297  else if (static_cast<ParticleType&>(*this).softImpact())
298  {
299  // Soft-sphere particles can travel outside the domain
300  // but we don't use lambda since this the particle
301  // is going away from face
302  trackFraction = 1.0;
303  position_ = endPosition;
304  }
305 
306  // change cell
307  if (internalFace) // Internal face
308  {
309  if (celli_ == mesh.faceOwner()[facei_])
310  {
311  celli_ = mesh.faceNeighbour()[facei_];
312  }
313  else if (celli_ == mesh.faceNeighbour()[facei_])
314  {
315  celli_ = mesh.faceOwner()[facei_];
316  }
317  else
318  {
320  (
321  "Particle::trackToFace(const vector&, TrackData&)"
322  )<< "addressing failure" << nl
323  << abort(FatalError);
324  }
325  }
326  else
327  {
328  ParticleType& p = static_cast<ParticleType&>(*this);
329 
330  // Soft-sphere algorithm ignores the boundary
331  if (p.softImpact())
332  {
333  trackFraction = 1.0;
334  position_ = endPosition;
335  }
336 
337  label patchi = patch(facei_);
338  const polyPatch& patch = mesh.boundaryMesh()[patchi];
339 
340  if (!p.hitPatch(patch, td, patchi))
341  {
342  if (isA<wedgePolyPatch>(patch))
343  {
344  p.hitWedgePatch
345  (
346  static_cast<const wedgePolyPatch&>(patch), td
347  );
348  }
349  else if (isA<symmetryPolyPatch>(patch))
350  {
351  p.hitSymmetryPatch
352  (
353  static_cast<const symmetryPolyPatch&>(patch), td
354  );
355  }
356  else if (isA<cyclicPolyPatch>(patch))
357  {
358  p.hitCyclicPatch
359  (
360  static_cast<const cyclicPolyPatch&>(patch), td
361  );
362  }
363  else if (isA<processorPolyPatch>(patch))
364  {
365  p.hitProcessorPatch
366  (
367  static_cast<const processorPolyPatch&>(patch), td
368  );
369  }
370  else if (isA<wallPolyPatch>(patch))
371  {
372  p.hitWallPatch
373  (
374  static_cast<const wallPolyPatch&>(patch), td
375  );
376  }
377  else
378  {
379  p.hitPatch(patch, td);
380  }
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-3*(mesh.cellCentres()[celli_] - position_);
394  }
395 
396  return trackFraction;
397 }
398 
399 template<class ParticleType>
401 (
402  const vector& endPosition
403 )
404 {
405  int dummyTd;
406  return trackToFace(endPosition, dummyTd);
407 }
408 
409 template<class ParticleType>
411 {
412  position_ = transform(T, position_);
413 }
414 
415 
416 template<class ParticleType>
418 {}
419 
420 
421 template<class ParticleType>
423 {}
424 
425 
426 template<class ParticleType>
427 template<class TrackData>
429 (
430  const polyPatch&,
431  TrackData&,
432  const label
433 )
434 {
435  return false;
436 }
437 
438 
439 template<class ParticleType>
440 template<class TrackData>
442 (
443  const wedgePolyPatch& wpp,
444  TrackData&
445 )
446 {
447  vector nf = wpp.faceAreas()[wpp.whichFace(facei_)];
448  nf /= mag(nf);
449 
450  static_cast<ParticleType&>(*this).transformProperties(I - 2.0*nf*nf);
451 }
452 
453 
454 template<class ParticleType>
455 template<class TrackData>
457 (
458  const symmetryPolyPatch& spp,
459  TrackData&
460 )
461 {
462  vector nf = spp.faceAreas()[spp.whichFace(facei_)];
463  nf /= mag(nf);
464 
465  static_cast<ParticleType&>(*this).transformProperties(I - 2.0*nf*nf);
466 }
467 
468 
469 template<class ParticleType>
470 template<class TrackData>
472 (
473  const cyclicPolyPatch& cpp,
474  TrackData&
475 )
476 {
477  label patchFacei_ = cpp.whichFace(facei_);
478 
479  facei_ = cpp.transformGlobalFace(facei_);
480 
481  celli_ = cloud_.polyMesh_.faceOwner()[facei_];
482 
483  if (!cpp.parallel())
484  {
485  const tensor& T = cpp.transformT(patchFacei_);
486 
487  transformPosition(T);
488  static_cast<ParticleType&>(*this).transformProperties(T);
489  }
490  else if (cpp.separated())
491  {
492  position_ += cpp.separation(patchFacei_);
493  static_cast<ParticleType&>(*this).transformProperties
494  (
495  cpp.separation(patchFacei_)
496  );
497  }
498 }
499 
500 
501 template<class ParticleType>
502 template<class TrackData>
504 (
505  const processorPolyPatch& spp,
506  TrackData& td
507 )
508 {}
509 
510 
511 template<class ParticleType>
512 template<class TrackData>
514 (
515  const wallPolyPatch& spp,
516  TrackData&
517 )
518 {}
519 
520 
521 template<class ParticleType>
522 template<class TrackData>
524 (
525  const polyPatch& spp,
526  TrackData&
527 )
528 {}
529 
530 
531 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
532 
533 #include "ParticleIO.C"
534 
535 // ************************ vim: set sw=4 sts=4 et: ************************ //