FreeFOAM The Cross-Platform CFD Toolkit
ExactParticle.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 <autoMesh/ExactParticle.H>
27 
28 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
29 
30 template<class ParticleType>
31 template<class TrackingData>
33 (
34  const vector& endPosition,
35  TrackingData& td
36 )
37 {
38  this->facei_ = -1;
39 
40  // Tracks to endPosition or stop on boundary
41  while (!this->onBoundary() && this->stepFraction_ < 1.0 - SMALL)
42  {
43  this->stepFraction_ +=
44  trackToFace(endPosition, td)
45  *(1.0 - this->stepFraction_);
46  }
47 
48  return this->facei_;
49 }
50 
51 
52 template<class ParticleType>
54 (
55  const vector& endPosition
56 )
57 {
58  int dummyTd;
59  return track(endPosition, dummyTd);
60 }
61 
62 
63 template<class ParticleType>
64 template<class TrackingData>
66 (
67  const vector& endPosition,
68  TrackingData& td
69 )
70 {
71  const polyMesh& mesh = this->cloud().pMesh();
72  const labelList& cFaces = mesh.cells()[this->celli_];
73 
74  point intersection(vector::zero);
75  scalar trackFraction = VGREAT;
76  label hitFacei = -1;
77 
78  const vector vec = endPosition-this->position_;
79 
80  forAll(cFaces, i)
81  {
82  label facei = cFaces[i];
83 
84  if (facei != this->face())
85  {
86  pointHit inter = mesh.faces()[facei].intersection
87  (
88  this->position_,
89  vec,
90  mesh.faceCentres()[facei],
91  mesh.points(),
92  intersection::HALF_RAY
93  );
94 
95  if (inter.hit() && inter.distance() < trackFraction)
96  {
97  trackFraction = inter.distance();
98  hitFacei = facei;
99  intersection = inter.hitPoint();
100  }
101  }
102  }
103 
104  if (hitFacei == -1)
105  {
106  // Did not find any intersection. Fall back to original approximate
107  // algorithm
109  (
110  endPosition,
111  td
112  );
113  }
114 
115  if (trackFraction >= (1.0-SMALL))
116  {
117  // Nearest intersection beyond endPosition so we hit endPosition.
118  trackFraction = 1.0;
119  this->position_ = endPosition;
120  this->facei_ = -1;
121  return 1.0;
122  }
123  else
124  {
125  this->position_ = intersection;
126  this->facei_ = hitFacei;
127  }
128 
129 
130  // Normal situation (trackFraction 0..1). Straight copy
131  // of Particle::trackToFace.
132 
133  bool internalFace = this->cloud().internalFace(this->facei_);
134 
135  // change cell
136  if (internalFace) // Internal face
137  {
138  if (this->celli_ == mesh.faceOwner()[this->facei_])
139  {
140  this->celli_ = mesh.faceNeighbour()[this->facei_];
141  }
142  else if (this->celli_ == mesh.faceNeighbour()[this->facei_])
143  {
144  this->celli_ = mesh.faceOwner()[this->facei_];
145  }
146  else
147  {
149  (
150  "ExactParticle::trackToFace"
151  "(const vector&, TrackingData&)"
152  )<< "addressing failure" << nl
153  << abort(FatalError);
154  }
155  }
156  else
157  {
158  ParticleType& p = static_cast<ParticleType&>(*this);
159 
160  // Soft-sphere algorithm ignores the boundary
161  if (p.softImpact())
162  {
163  trackFraction = 1.0;
164  this->position_ = endPosition;
165  }
166 
167  label patchi = patch(this->facei_);
168  const polyPatch& patch = mesh.boundaryMesh()[patchi];
169 
170  if (isA<wedgePolyPatch>(patch))
171  {
172  p.hitWedgePatch
173  (
174  static_cast<const wedgePolyPatch&>(patch), td
175  );
176  }
177  else if (isA<symmetryPolyPatch>(patch))
178  {
179  p.hitSymmetryPatch
180  (
181  static_cast<const symmetryPolyPatch&>(patch), td
182  );
183  }
184  else if (isA<cyclicPolyPatch>(patch))
185  {
186  p.hitCyclicPatch
187  (
188  static_cast<const cyclicPolyPatch&>(patch), td
189  );
190  }
191  else if (isA<processorPolyPatch>(patch))
192  {
193  p.hitProcessorPatch
194  (
195  static_cast<const processorPolyPatch&>(patch), td
196  );
197  }
198  else if (isA<wallPolyPatch>(patch))
199  {
200  p.hitWallPatch
201  (
202  static_cast<const wallPolyPatch&>(patch), td
203  );
204  }
205  else if (isA<polyPatch>(patch))
206  {
207  p.hitPatch
208  (
209  static_cast<const polyPatch&>(patch), td
210  );
211  }
212  else
213  {
215  (
216  "ExactParticle::trackToFace"
217  "(const vector& endPosition, scalar& trackFraction)"
218  )<< "patch type " << patch.type() << " not suported" << nl
219  << abort(FatalError);
220  }
221  }
222 
223  // If the trackFraction = 0 something went wrong.
224  // Either the particle is flipping back and forth across a face perhaps
225  // due to velocity interpolation errors or it is in a "hole" in the mesh
226  // caused by face warpage.
227  // In both cases resolve the positional ambiguity by moving the particle
228  // slightly towards the cell-centre.
229  if (trackFraction < SMALL)
230  {
231  this->position_ +=
232  1.0e-6*(mesh.cellCentres()[this->celli_] - this->position_);
233  }
234 
235  return trackFraction;
236 }
237 
238 
239 template<class ParticleType>
241 (
242  const vector& endPosition
243 )
244 {
245  int dummyTd;
246  return trackToFace(endPosition, dummyTd);
247 }
248 
249 
250 template<class ParticleType>
251 Foam::Ostream& Foam::operator<<
252 (
253  Ostream& os,
255 )
256 {
257  return operator<<(os, static_cast<const Particle<ParticleType>&>(p));
258 }
259 
260 
261 // ************************ vim: set sw=4 sts=4 et: ************************ //