1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
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
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 \*---------------------------------------------------------------------------*/
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
41 const vector& position
44 const polyMesh& mesh = cloud_.polyMesh_;
45 const labelList& faces = mesh.cells()[celli_];
46 const vector& C = mesh.cellCentres()[celli_];
48 labelList faceList(0);
51 label facei = faces[i];
52 scalar lam = lambda(C, position, facei);
54 if ((lam > 0) && (lam < 1.0))
56 label n = faceList.size();
57 faceList.setSize(n+1);
66 template<class ParticleType>
67 Foam::labelList Foam::Particle<ParticleType>::findFaces
69 const vector& position,
71 const scalar stepFraction
74 const polyMesh& mesh = cloud_.polyMesh_;
75 const labelList& faces = mesh.cells()[celli];
76 const vector& C = mesh.cellCentres()[celli];
78 labelList faceList(0);
81 label facei = faces[i];
82 scalar lam = lambda(C, position, facei, stepFraction);
84 if ((lam > 0) && (lam < 1.0))
86 label n = faceList.size();
87 faceList.setSize(n+1);
96 template<class ParticleType>
97 template<class TrackData>
98 void Foam::Particle<ParticleType>::prepareForParallelTransfer
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
117 const processorPolyPatch& ppp =
118 refCast<const processorPolyPatch>
119 (cloud_.pMesh().boundaryMesh()[patchi]);
121 celli_ = ppp.faceCells()[facei_];
125 if (ppp.forwardT().size() == 1)
127 const tensor& T = ppp.forwardT()[0];
128 transformPosition(T);
129 static_cast<ParticleType&>(*this).transformProperties(T);
133 const tensor& T = ppp.forwardT()[facei_];
134 transformPosition(T);
135 static_cast<ParticleType&>(*this).transformProperties(T);
138 else if (ppp.separated())
140 if (ppp.separation().size() == 1)
142 position_ -= ppp.separation()[0];
143 static_cast<ParticleType&>(*this).transformProperties
150 position_ -= ppp.separation()[facei_];
151 static_cast<ParticleType&>(*this).transformProperties
153 -ppp.separation()[facei_]
158 // Reset the face index for the next tracking operation
159 if (stepFraction_ > (1.0 - SMALL))
166 facei_ += ppp.start();
171 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
173 template<class ParticleType>
174 Foam::Particle<ParticleType>::Particle
176 const Cloud<ParticleType>& cloud,
177 const vector& position,
189 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
191 template<class ParticleType>
192 template<class TrackData>
193 Foam::label Foam::Particle<ParticleType>::track
195 const vector& endPosition,
201 // Tracks to endPosition or stop on boundary
202 while(!onBoundary() && stepFraction_ < 1.0 - SMALL)
204 stepFraction_ += trackToFace(endPosition, td)*(1.0 - stepFraction_);
212 template<class ParticleType>
213 Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
216 return track(endPosition, dummyTd);
219 template<class ParticleType>
220 template<class TrackData>
221 Foam::scalar Foam::Particle<ParticleType>::trackToFace
223 const vector& endPosition,
227 const polyMesh& mesh = cloud_.polyMesh_;
229 labelList faces = findFaces(endPosition);
232 scalar trackFraction = 0.0;
234 if (faces.size() == 0) // inside cell
237 position_ = endPosition;
241 scalar lambdaMin = GREAT;
243 if (faces.size() == 1)
245 lambdaMin = lambda(position_, endPosition, faces[0], stepFraction_);
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
258 lambda(position_, endPosition, faces[i], stepFraction_);
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.
277 if (lambdaMin <= 1.0)
279 trackFraction = lambdaMin;
280 position_ += trackFraction*(endPosition - position_);
285 position_ = endPosition;
288 else if (static_cast<ParticleType&>(*this).softImpact())
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
294 position_ = endPosition;
298 if (internalFace) // Internal face
300 if (celli_ == cloud_.owner_[facei_])
302 celli_ = cloud_.neighbour_[facei_];
304 else if (celli_ == cloud_.neighbour_[facei_])
306 celli_ = cloud_.owner_[facei_];
312 "Particle::trackToFace(const vector&, TrackData&)"
313 )<< "addressing failure" << nl
314 << abort(FatalError);
319 ParticleType& p = static_cast<ParticleType&>(*this);
321 // Soft-sphere algorithm ignores the boundary
325 position_ = endPosition;
328 label patchi = patch(facei_);
329 const polyPatch& patch = mesh.boundaryMesh()[patchi];
331 if (isA<wedgePolyPatch>(patch))
335 static_cast<const wedgePolyPatch&>(patch), td
338 else if (isA<symmetryPolyPatch>(patch))
342 static_cast<const symmetryPolyPatch&>(patch), td
345 else if (isA<cyclicPolyPatch>(patch))
349 static_cast<const cyclicPolyPatch&>(patch), td
352 else if (isA<processorPolyPatch>(patch))
356 static_cast<const processorPolyPatch&>(patch), td
359 else if (isA<wallPolyPatch>(patch))
363 static_cast<const wallPolyPatch&>(patch), td
366 else if (isA<polyPatch>(patch))
370 static_cast<const polyPatch&>(patch), td
377 "Particle::trackToFace"
378 "(const vector& endPosition, scalar& trackFraction)"
379 )<< "patch type " << patch.type() << " not suported" << nl
380 << abort(FatalError);
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)
393 position_ += 1.0e-6*(mesh.cellCentres()[celli_] - position_);
396 return trackFraction;
399 template<class ParticleType>
400 Foam::scalar Foam::Particle<ParticleType>::trackToFace
402 const vector& endPosition
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&)
421 template<class ParticleType>
422 void Foam::Particle<ParticleType>::transformProperties(const vector&)
426 template<class ParticleType>
427 template<class TrackData>
428 void Foam::Particle<ParticleType>::hitWedgePatch
430 const wedgePolyPatch& wpp,
434 vector nf = wpp.faceAreas()[wpp.whichFace(facei_)];
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,
449 vector nf = spp.faceAreas()[spp.whichFace(facei_)];
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,
464 label patchFacei_ = cpp.whichFace(facei_);
466 facei_ = cpp.transformGlobalFace(facei_);
468 celli_ = cloud_.polyMesh_.faceOwner()[facei_];
472 const tensor& T = cpp.transformT(patchFacei_);
474 transformPosition(T);
475 static_cast<ParticleType&>(*this).transformProperties(T);
477 else if (cpp.separated())
479 position_ += cpp.separation(patchFacei_);
480 static_cast<ParticleType&>(*this).transformProperties
482 cpp.separation(patchFacei_)
488 template<class ParticleType>
489 template<class TrackData>
490 void Foam::Particle<ParticleType>::hitProcessorPatch
492 const processorPolyPatch& spp,
498 template<class ParticleType>
499 template<class TrackData>
500 void Foam::Particle<ParticleType>::hitWallPatch
502 const wallPolyPatch& spp,
508 template<class ParticleType>
509 template<class TrackData>
510 void Foam::Particle<ParticleType>::hitPatch
512 const polyPatch& spp,
518 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
520 #include "ParticleIO.C"
522 // ************************************************************************* //