FreeFOAM The Cross-Platform CFD Toolkit
ParticleI.H
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 <OpenFOAM/polyMesh.H>
27 
28 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
29 
30 template<class ParticleType>
31 inline Foam::scalar Foam::Particle<ParticleType>::lambda
32 (
33  const vector& from,
34  const vector& to,
35  const label facei,
36  const scalar stepFraction
37 ) const
38 {
39  const polyMesh& mesh = cloud_.polyMesh_;
40  bool movingMesh = mesh.moving();
41 
42  if (movingMesh)
43  {
44  vector Sf = mesh.faceAreas()[facei];
45  Sf /= mag(Sf);
46  vector Cf = mesh.faceCentres()[facei];
47 
48  // move reference point for wall
49  if (!cloud_.internalFace(facei))
50  {
51  label patchi = patch(facei);
52  const polyPatch& patch = mesh.boundaryMesh()[patchi];
53 
54  // move reference point for wall
55  if (isA<wallPolyPatch>(patch))
56  {
57  const vector& C = mesh.cellCentres()[celli_];
58  scalar CCf = mag((C - Cf) & Sf);
59  // check if distance between cell centre and face centre
60  // is larger than wallImpactDistance
61  const ParticleType& p =
62  static_cast<const ParticleType&>(*this);
63  if (CCf > p.wallImpactDistance(Sf))
64  {
65  Cf -=p.wallImpactDistance(Sf)*Sf;
66  }
67  }
68  }
69 
70  // for a moving mesh we need to reconstruct the old
71  // Sf and Cf from oldPoints (they aren't stored)
72 
73  const vectorField& oldPoints = mesh.oldPoints();
74 
75  vector Cf00 = mesh.faces()[facei].centre(oldPoints);
76  vector Cf0 = Cf00 + stepFraction*(Cf - Cf00);
77 
78  vector Sf00 = mesh.faces()[facei].normal(oldPoints);
79 
80  // for layer addition Sf00 = vector::zero and we use Sf
81  if (mag(Sf00) > SMALL)
82  {
83  Sf00 /= mag(Sf00);
84  }
85  else
86  {
87  Sf00 = Sf;
88  }
89 
90  scalar magSfDiff = mag(Sf - Sf00);
91 
92  // check if the face is rotating
93  if (magSfDiff > SMALL)
94  {
95  vector Sf0 = Sf00 + stepFraction*(Sf - Sf00);
96 
97  // find center of rotation
98  vector omega = Sf0 ^ Sf;
99  scalar magOmega = mag(omega);
100  omega /= magOmega + SMALL;
101  vector n0 = omega ^ Sf0;
102  scalar lam = ((Cf - Cf0) & Sf)/(n0 & Sf);
103  vector r0 = Cf0 + lam*n0;
104 
105  // solve '(p - r0) & Sfp = 0', where
106  // p = from + lambda*(to - from)
107  // Sfp = Sf0 + lambda*(Sf - Sf0)
108  // which results in the quadratic eq.
109  // a*lambda^2 + b*lambda + c = 0
110  vector alpha = from - r0;
111  vector beta = to - from;
112  scalar a = beta & (Sf - Sf0);
113  scalar b = (alpha & (Sf - Sf0)) + (beta & Sf0);
114  scalar c = alpha & Sf0;
115 
116  if (mag(a) > SMALL)
117  {
118  // solve the second order polynomial
119  scalar ap = b/a;
120  scalar bp = c/a;
121  scalar cp = ap*ap - 4.0*bp;
122 
123  if (cp < 0.0)
124  {
125  // imaginary roots only
126  return GREAT;
127  }
128  else
129  {
130  scalar l1 = -0.5*(ap - ::sqrt(cp));
131  scalar l2 = -0.5*(ap + ::sqrt(cp));
132 
133  // one root is around 0-1, while
134  // the other is very large in mag
135  if (mag(l1) < mag(l2))
136  {
137  return l1;
138  }
139  else
140  {
141  return l2;
142  }
143  }
144  }
145  else
146  {
147  // when a==0, solve the first order polynomial
148  return (-c/b);
149  }
150  }
151  else // no rotation
152  {
153  vector alpha = from - Cf0;
154  vector beta = to - from - (Cf - Cf0);
155  scalar lambdaNominator = alpha & Sf;
156  scalar lambdaDenominator = beta & Sf;
157 
158  // check if trajectory is parallel to face
159  if (mag(lambdaDenominator) < SMALL)
160  {
161  if (lambdaDenominator < 0.0)
162  {
163  lambdaDenominator = -SMALL;
164  }
165  else
166  {
167  lambdaDenominator = SMALL;
168  }
169  }
170 
171  return (-lambdaNominator/lambdaDenominator);
172  }
173  }
174  else
175  {
176  // mesh is static and stepFraction is not needed
177  return lambda(from, to, facei);
178  }
179 }
180 
181 
182 template<class ParticleType>
183 inline Foam::scalar Foam::Particle<ParticleType>::lambda
184 (
185  const vector& from,
186  const vector& to,
187  const label facei
188 ) const
189 {
190  const polyMesh& mesh = cloud_.polyMesh_;
191 
192  vector Sf = mesh.faceAreas()[facei];
193  Sf /= mag(Sf);
194  vector Cf = mesh.faceCentres()[facei];
195 
196  // move reference point for wall
197  if (!cloud_.internalFace(facei))
198  {
199  label patchi = patch(facei);
200  const polyPatch& patch = mesh.boundaryMesh()[patchi];
201 
202  // move reference point for wall
203  if (isA<wallPolyPatch>(patch))
204  {
205  const vector& C = mesh.cellCentres()[celli_];
206  scalar CCf = mag((C - Cf) & Sf);
207  // check if distance between cell centre and face centre
208  // is larger than wallImpactDistance
209  const ParticleType& p = static_cast<const ParticleType&>(*this);
210  if (CCf > p.wallImpactDistance(Sf))
211  {
212  Cf -=p.wallImpactDistance(Sf)*Sf;
213  }
214  }
215  }
216 
217  scalar lambdaNominator = (Cf - from) & Sf;
218  scalar lambdaDenominator = (to - from) & Sf;
219 
220  // check if trajectory is parallel to face
221  if (mag(lambdaDenominator) < SMALL)
222  {
223  if (lambdaDenominator < 0.0)
224  {
225  lambdaDenominator = -SMALL;
226  }
227  else
228  {
229  lambdaDenominator = SMALL;
230  }
231  }
232 
233  return lambdaNominator/lambdaDenominator;
234 }
235 
236 
237 template<class ParticleType>
239 {
240  DynamicList<label>& faces = cloud_.labels_;
241  findFaces(position_, faces);
242 
243  return (!faces.size());
244 }
245 
246 
247 template<class ParticleType>
249 (
250  const vector& position,
251  const label celli,
252  const scalar stepFraction
253 ) const
254 {
255  DynamicList<label>& faces = cloud_.labels_;
256  findFaces(position, celli, stepFraction, faces);
257 
258  return (!faces.size());
259 }
260 
261 
262 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
263 
264 template<class ParticleType>
266 (
268 )
269 :
270  cloud_(cloud)
271 {}
272 
273 
274 template<class ParticleType>
277 {
278  return cloud_;
279 }
280 
281 
282 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
283 
284 template<class ParticleType>
285 inline const Foam::Cloud<ParticleType>&
287 {
288  return cloud_;
289 }
290 
291 
292 template<class ParticleType>
294 {
295  return position_;
296 }
297 
298 
299 template<class ParticleType>
301 {
302  return position_;
303 }
304 
305 
306 template<class ParticleType>
307 inline Foam::label Foam::Particle<ParticleType>::cell() const
308 {
309  return celli_;
310 }
311 
312 
313 template<class ParticleType>
315 {
316  return celli_;
317 }
318 
319 
320 template<class ParticleType>
321 inline Foam::label Foam::Particle<ParticleType>::face() const
322 {
323  return facei_;
324 }
325 
326 
327 template<class ParticleType>
329 {
330  return facei_ != -1 && facei_ >= cloud_.pMesh().nInternalFaces();
331 }
332 
333 
334 template<class ParticleType>
336 {
337  return stepFraction_;
338 }
339 
340 
341 template<class ParticleType>
343 {
344  return stepFraction_;
345 }
346 
347 
348 template<class ParticleType>
349 inline Foam::label Foam::Particle<ParticleType>::origProc() const
350 {
351  return origProc_;
352 }
353 
354 
355 template<class ParticleType>
356 inline Foam::label Foam::Particle<ParticleType>::origId() const
357 {
358  return origId_;
359 }
360 
361 
362 template<class ParticleType>
364 {
365  return false;
366 }
367 
368 
369 template<class ParticleType>
370 inline Foam::scalar Foam::Particle<ParticleType>::currentTime() const
371 {
372  return
373  cloud_.pMesh().time().value()
374  + stepFraction_*cloud_.pMesh().time().deltaT().value();
375 }
376 
377 
378 template<class ParticleType>
379 inline Foam::label Foam::Particle<ParticleType>::patch(const label facei) const
380 {
381  return cloud_.facePatch(facei);
382 }
383 
384 
385 template<class ParticleType>
387 (
388  const label patchi,
389  const label facei
390 ) const
391 {
392  return cloud_.patchFace(patchi, facei);
393 }
394 
395 
396 template<class ParticleType>
397 inline Foam::scalar
399 {
400  return 0.0;
401 }
402 
403 
404 template<class ParticleType>
406 {
407  return facei_;
408 }
409 
410 
411 // ************************ vim: set sw=4 sts=4 et: ************************ //