wallBoundedParticle.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) 2011-2016 OpenFOAM Foundation
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 "wallBoundedParticle.H"
27 
28 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
29 
30 const std::size_t Foam::wallBoundedParticle::sizeofFields_
31 (
32  sizeof(wallBoundedParticle) - sizeof(particle)
33 );
34 
35 
36 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
37 
39 {
40  if ((meshEdgeStart_ != -1) == (diagEdge_ != -1))
41  {
43  << "Particle:"
44  << info()
45  << "cannot both be on a mesh edge and a face-diagonal edge."
46  << " meshEdgeStart_:" << meshEdgeStart_
47  << " diagEdge_:" << diagEdge_
48  << abort(FatalError);
49  }
50 
51  const Foam::face& f = mesh_.faces()[tetFace()];
52 
53  if (meshEdgeStart_ != -1)
54  {
55  return edge(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_));
56  }
57  else
58  {
59  label faceBasePtI = mesh_.tetBasePtIs()[tetFace()];
60  label diagPtI = (faceBasePtI+diagEdge_)%f.size();
61 
62  return edge(f[faceBasePtI], f[diagPtI]);
63  }
64 }
65 
66 
68 (
69  const edge& meshEdge
70 )
71 {
72  // Update tetFace, tetPt
74 
75  // Update face to be same as tracking one
76  face() = tetFace();
77 
78  // And adapt meshEdgeStart_.
79  const Foam::face& f = mesh_.faces()[tetFace()];
80  label fp = findIndex(f, meshEdge[0]);
81 
82  if (f.nextLabel(fp) == meshEdge[1])
83  {
84  meshEdgeStart_ = fp;
85  }
86  else
87  {
88  label fpMin1 = f.rcIndex(fp);
89 
90  if (f[fpMin1] == meshEdge[1])
91  {
92  meshEdgeStart_ = fpMin1;
93  }
94  else
95  {
97  << "Problem :"
98  << " particle:"
99  << info()
100  << "face:" << tetFace()
101  << " verts:" << f
102  << " meshEdge:" << meshEdge
103  << abort(FatalError);
104  }
105  }
106 
107  diagEdge_ = -1;
108 
109  // Check that still on same mesh edge
110  const edge eNew(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_));
111  if (eNew != meshEdge)
112  {
114  << "Problem" << abort(FatalError);
115  }
116 }
117 
118 
120 {
121  if (diagEdge_ == -1)
122  {
124  << "Particle:"
125  << info()
126  << "not on a diagonal edge" << abort(FatalError);
127  }
128  if (meshEdgeStart_ != -1)
129  {
131  << "Particle:"
132  << info()
133  << "meshEdgeStart_:" << meshEdgeStart_ << abort(FatalError);
134  }
135 
136  const Foam::face& f = mesh_.faces()[tetFace()];
137 
138  // tetPtI starts from 1, goes up to f.size()-2
139 
140  if (tetPt() == diagEdge_)
141  {
142  tetPt() = f.rcIndex(tetPt());
143  }
144  else
145  {
146  label nextTetPt = f.fcIndex(tetPt());
147  if (diagEdge_ == nextTetPt)
148  {
149  tetPt() = nextTetPt;
150  }
151  else
152  {
154  << "Particle:"
155  << info()
156  << "tetPt:" << tetPt()
157  << " diagEdge:" << diagEdge_ << abort(FatalError);
158  }
159  }
160 
161  meshEdgeStart_ = -1;
162 }
163 
164 
166 (
167  const vector& endPosition,
168  label& minEdgeI
169 )
170 {
171  // Track p from position to endPosition
172  const triFace tri(currentTetIndices().faceTriIs(mesh_));
173  vector n = tri.normal(mesh_.points());
174  n /= mag(n)+VSMALL;
175 
176  // Check which edge intersects the trajectory.
177  // Project trajectory onto triangle
178  minEdgeI = -1;
179  scalar minS = 1; // end position
180 
181  edge currentE(-1, -1);
182  if (meshEdgeStart_ != -1 || diagEdge_ != -1)
183  {
184  currentE = currentEdge();
185  }
186 
187  // Determine path along line position+s*d to see where intersections
188  // are.
189 
190  forAll(tri, i)
191  {
192  label j = tri.fcIndex(i);
193 
194  const point& pt0 = mesh_.points()[tri[i]];
195  const point& pt1 = mesh_.points()[tri[j]];
196 
197  if (edge(tri[i], tri[j]) == currentE)
198  {
199  // Do not check particle is on
200  continue;
201  }
202 
203  // Outwards pointing normal
204  vector edgeNormal = (pt1-pt0)^n;
205 
206  edgeNormal /= mag(edgeNormal)+VSMALL;
207 
208  // Determine whether position and end point on either side of edge.
209  scalar sEnd = (endPosition-pt0)&edgeNormal;
210  if (sEnd >= 0)
211  {
212  // endPos is outside triangle. position() should always be
213  // inside.
214  scalar sStart = (position()-pt0)&edgeNormal;
215  if (mag(sEnd-sStart) > VSMALL)
216  {
217  scalar s = sStart/(sStart-sEnd);
218 
219  if (s >= 0 && s < minS)
220  {
221  minS = s;
222  minEdgeI = i;
223  }
224  }
225  }
226  }
227 
228  if (minEdgeI != -1)
229  {
230  position() += minS*(endPosition-position());
231  }
232  else
233  {
234  // Did not hit any edge so tracked to the end position. Set position
235  // without any calculation to avoid truncation errors.
236  position() = endPosition;
237  minS = 1.0;
238  }
239 
240  // Project position() back onto plane of triangle
241  const point& triPt = mesh_.points()[tri[0]];
242  position() -= ((position()-triPt)&n)*n;
243 
244  return minS;
245 }
246 
247 
249 (
250  const point& endPosition
251 ) const
252 {
253  const triFace triVerts(currentTetIndices().faceTriIs(mesh_));
254  const edge currentE = currentEdge();
255 
256  //if (debug)
257  {
258  if
259  (
260  currentE[0] == currentE[1]
261  || findIndex(triVerts, currentE[0]) == -1
262  || findIndex(triVerts, currentE[1]) == -1
263  )
264  {
266  << "Edge " << currentE << " not on triangle " << triVerts
267  << info()
268  << abort(FatalError);
269  }
270  }
271 
272 
273  const vector dir = endPosition-position();
274 
275  // Get normal of currentE
276  vector n = triVerts.normal(mesh_.points());
277  n /= mag(n);
278 
279  forAll(triVerts, i)
280  {
281  label j = triVerts.fcIndex(i);
282  const point& pt0 = mesh_.points()[triVerts[i]];
283  const point& pt1 = mesh_.points()[triVerts[j]];
284 
285  if (edge(triVerts[i], triVerts[j]) == currentE)
286  {
287  vector edgeNormal = (pt1-pt0)^n;
288  return (dir&edgeNormal) < 0;
289  }
290  }
291 
293  << "Problem" << abort(FatalError);
294 
295  return false;
296 }
297 
298 
299 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
300 
302 (
303  const polyMesh& mesh,
304  const vector& position,
305  const label celli,
306  const label tetFacei,
307  const label tetPtI,
308  const label meshEdgeStart,
309  const label diagEdge
310 )
311 :
312  particle(mesh, position, celli, tetFacei, tetPtI),
313  meshEdgeStart_(meshEdgeStart),
314  diagEdge_(diagEdge)
315 {}
316 
317 
319 (
320  const polyMesh& mesh,
321  Istream& is,
322  bool readFields
323 )
324 :
325  particle(mesh, is, readFields)
326 {
327  if (readFields)
328  {
329  if (is.format() == IOstream::ASCII)
330  {
331  is >> meshEdgeStart_ >> diagEdge_;
332  }
333  else
334  {
335  is.read
336  (
337  reinterpret_cast<char*>(&meshEdgeStart_),
338  sizeof(meshEdgeStart_)
339  + sizeof(diagEdge_)
340  );
341  }
342  }
343 
344  // Check state of Istream
345  is.check
346  (
347  "wallBoundedParticle::wallBoundedParticle"
348  "(const Cloud<wallBoundedParticle>&, Istream&, bool)"
349  );
350 }
351 
352 
354 (
355  const wallBoundedParticle& p
356 )
357 :
358  particle(p),
361 {}
362 
363 
364 // * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
365 
366 Foam::Ostream& Foam::operator<<
367 (
368  Ostream& os,
369  const wallBoundedParticle& p
370 )
371 {
372  if (os.format() == IOstream::ASCII)
373  {
374  os << static_cast<const particle&>(p)
376  << token::SPACE << p.diagEdge_;
377  }
378  else
379  {
380  os << static_cast<const particle&>(p);
381  os.write
382  (
383  reinterpret_cast<const char*>(&p.meshEdgeStart_),
384  wallBoundedParticle::sizeofFields_
385  );
386  }
387 
388  return os;
389 }
390 
391 
392 Foam::Ostream& Foam::operator<<
393 (
394  Ostream& os,
396 )
397 {
398  const wallBoundedParticle& p = ip.t_;
399 
400  tetPointRef tpr(p.currentTet());
401 
402  os << " " << static_cast<const particle&>(p) << nl
403  << " tet:" << nl;
404  os << " ";
405  meshTools::writeOBJ(os, tpr.a());
406  os << " ";
407  meshTools::writeOBJ(os, tpr.b());
408  os << " ";
409  meshTools::writeOBJ(os, tpr.c());
410  os << " ";
411  meshTools::writeOBJ(os, tpr.d());
412  os << " l 1 2" << nl
413  << " l 1 3" << nl
414  << " l 1 4" << nl
415  << " l 2 3" << nl
416  << " l 2 4" << nl
417  << " l 3 4" << nl;
418  os << " ";
419  meshTools::writeOBJ(os, p.position());
420 
421  return os;
422 }
423 
424 
425 // ************************************************************************* //
label meshEdgeStart_
Particle is on mesh edge:
label & tetPt()
Return current tet face particle is in.
Definition: particleI.H:628
A tetrahedron primitive.
Definition: tetrahedron.H:62
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
streamFormat format() const
Return current stream format.
Definition: IOstream.H:377
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
An Istream is an abstract base class for all input systems (streams, files, token lists etc)...
Definition: Istream.H:57
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
virtual bool check(const char *operation) const
Check IOstream status for given operation.
Definition: IOstream.C:92
void crossEdgeConnectedFace(const label &celli, label &tetFacei, label &tetPtI, const edge &e)
Cross the from the given face across the given edge of the.
Definition: particleI.H:458
tetIndices currentTetIndices() const
Return the indices of the current tet that the.
Definition: particleI.H:634
edge currentEdge() const
Construct current edge.
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
particle(const polyMesh &mesh, const vector &position, const label celli, const label tetFacei, const label tetPtI)
Construct from components.
Definition: particle.C:48
Base particle class.
Definition: particle.H:78
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:979
vector normal(const pointField &) const
Vector normal; magnitude is equal to area of face.
Definition: triFaceI.H:181
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
tetPointRef currentTet() const
Return the geometry of the current tet that the.
Definition: particleI.H:640
const vector & position() const
Return current particle position.
Definition: particleI.H:586
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
virtual Istream & read(token &)=0
Return next token from stream.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
A triangular face using a FixedList of labels corresponding to mesh vertices.
Definition: triFace.H:68
const polyMesh & mesh_
Reference to the polyMesh database.
Definition: particle.H:137
label diagEdge_
Particle is on diagonal edge:
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: FixedListI.H:115
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Particle class that tracks on triangles of boundary faces. Use trackToEdge similar to trackToFace on ...
bool isTriAlongTrack(const point &endPosition) const
Is current triangle in the track direction.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
label & face()
Return current face particle is on otherwise -1.
Definition: particleI.H:664
static const char nl
Definition: Ostream.H:262
const labelList & tetBasePtIs() const
Return the tetBasePtIs.
Definition: polyMesh.C:818
void crossDiagonalEdge()
Cross diagonal edge into different triangle on same face,cell.
labelList f(nPoints)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
label rcIndex(const label i) const
Return the reverse circular index, i.e. the previous index.
Definition: UListI.H:65
void crossEdgeConnectedFace(const edge &meshEdge)
Check if inside current tet.
A helper class for outputting values to Ostream.
Definition: InfoProxy.H:45
wallBoundedParticle(const polyMesh &c, const vector &position, const label celli, const label tetFacei, const label tetPtI, const label meshEdgeStart, const label diagEdge)
Construct from components.
scalar trackFaceTri(const vector &endPosition, label &minEdgeI)
Track through single triangle.
label & cell()
Return current cell particle is in.
Definition: particleI.H:604
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
label nextLabel(const label i) const
Next vertex on face.
Definition: faceI.H:117
virtual Ostream & write(const token &)=0
Write next token to stream.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
label & tetFace()
Return current tet face particle is in.
Definition: particleI.H:616
volScalarField & p
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1004
InfoProxy< wallBoundedParticle > info() const
Return info proxy.