faceIntersection.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 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 "face.H"
27 #include "pointHit.H"
28 #include "triPointRef.H"
29 #include "line.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 // Return potential intersection with face with a ray starting
34 // at p, direction n (does not need to be normalized)
35 // Does face-center decomposition and returns triangle intersection
36 // point closest to p.
37 
38 // In case of miss the point is the nearest point intersection of the
39 // face plane and the ray and the distance is the distance between the
40 // intersection point and the nearest point on the face
41 
43 (
44  const point& p,
45  const vector& n,
46  const pointField& meshPoints,
47  const intersection::algorithm alg,
48  const intersection::direction dir
49 ) const
50 {
51  // If the face is a triangle, do a direct calculation
52  if (size() == 3)
53  {
54  return triPointRef
55  (
56  meshPoints[operator[](0)],
57  meshPoints[operator[](1)],
58  meshPoints[operator[](2)]
59  ).ray(p, n, alg, dir);
60  }
61 
62  point ctr = Foam::average(points(meshPoints));
63 
64  scalar nearestHitDist = GREAT;
65 
66  scalar nearestMissDist = GREAT;
67  bool eligible = false;
68 
69  // Initialize to miss, distance = GREAT
70  pointHit nearest(p);
71 
72  const labelList& f = *this;
73 
74  label nPoints = size();
75 
76  point nextPoint = ctr;
77 
78  for (label pI = 0; pI < nPoints; pI++)
79  {
80  nextPoint = meshPoints[f[fcIndex(pI)]];
81 
82  // Note: for best accuracy, centre point always comes last
83  //
84  pointHit curHit = triPointRef
85  (
86  meshPoints[f[pI]],
87  nextPoint,
88  ctr
89  ).ray(p, n, alg, dir);
90 
91  if (curHit.hit())
92  {
93  if (Foam::mag(curHit.distance()) < Foam::mag(nearestHitDist))
94  {
95  nearestHitDist = curHit.distance();
96  nearest.setHit();
97  nearest.setPoint(curHit.hitPoint());
98  }
99  }
100  else if (!nearest.hit())
101  {
102  // Miss and no hit yet. Update miss statistics.
103  if (curHit.eligibleMiss())
104  {
105  eligible = true;
106 
107  // Miss distance is the distance between the plane intersection
108  // point and the nearest point of the triangle
109  scalar missDist =
110  Foam::mag
111  (
112  p + curHit.distance()*n
113  - curHit.missPoint()
114  );
115 
116  if (missDist < nearestMissDist)
117  {
118  nearestMissDist = missDist;
119  nearest.setDistance(curHit.distance());
120  nearest.setPoint(curHit.missPoint());
121  }
122  }
123  }
124  }
125 
126  if (nearest.hit())
127  {
128  nearest.setDistance(nearestHitDist);
129  }
130  else
131  {
132  // Haven't hit a single face triangle
133  nearest.setMiss(eligible);
134  }
135 
136  return nearest;
137 }
138 
139 
141 (
142  const point& p,
143  const vector& q,
144  const point& ctr,
145  const pointField& meshPoints,
146  const intersection::algorithm alg,
147  const scalar tol
148 ) const
149 {
150  // If the face is a triangle, do a direct calculation
151  if (size() == 3)
152  {
153  return triPointRef
154  (
155  meshPoints[operator[](0)],
156  meshPoints[operator[](1)],
157  meshPoints[operator[](2)]
158  ).intersection(p, q, alg, tol);
159  }
160 
161  scalar nearestHitDist = VGREAT;
162 
163  // Initialize to miss, distance = GREAT
164  pointHit nearest(p);
165 
166  const labelList& f = *this;
167 
168  forAll(f, pI)
169  {
170  // Note: for best accuracy, centre point always comes last
171  pointHit curHit = triPointRef
172  (
173  meshPoints[f[pI]],
174  meshPoints[f[fcIndex(pI)]],
175  ctr
176  ).intersection(p, q, alg, tol);
177 
178  if (curHit.hit())
179  {
180  if (Foam::mag(curHit.distance()) < Foam::mag(nearestHitDist))
181  {
182  nearestHitDist = curHit.distance();
183  nearest.setHit();
184  nearest.setPoint(curHit.hitPoint());
185  }
186  }
187  }
188 
189  if (nearest.hit())
190  {
191  nearest.setDistance(nearestHitDist);
192  }
193 
194  return nearest;
195 }
196 
197 
199 (
200  const point& p,
201  const pointField& meshPoints
202 ) const
203 {
204  // Dummy labels
205  label nearType = -1;
206  label nearLabel = -1;
207 
208  return nearestPointClassify(p, meshPoints, nearType, nearLabel);
209 }
210 
211 
213 (
214  const point& p,
215  const pointField& meshPoints,
216  label& nearType,
217  label& nearLabel
218 ) const
219 {
220  // If the face is a triangle, do a direct calculation
221  if (size() == 3)
222  {
223  return triPointRef
224  (
225  meshPoints[operator[](0)],
226  meshPoints[operator[](1)],
227  meshPoints[operator[](2)]
228  ).nearestPointClassify(p, nearType, nearLabel);
229  }
230 
231  const face& f = *this;
232  point ctr = centre(meshPoints);
233 
234  // Initialize to miss, distance=GREAT
235  pointHit nearest(p);
236 
237  nearType = -1;
238  nearLabel = -1;
239 
240  label nPoints = f.size();
241 
242  point nextPoint = ctr;
243 
244  for (label pI = 0; pI < nPoints; pI++)
245  {
246  nextPoint = meshPoints[f[fcIndex(pI)]];
247 
248  label tmpNearType = -1;
249  label tmpNearLabel = -1;
250 
251  // Note: for best accuracy, centre point always comes last
252  triPointRef tri
253  (
254  meshPoints[f[pI]],
255  nextPoint,
256  ctr
257  );
258 
259  pointHit curHit = tri.nearestPointClassify
260  (
261  p,
262  tmpNearType,
263  tmpNearLabel
264  );
265 
266  if (Foam::mag(curHit.distance()) < Foam::mag(nearest.distance()))
267  {
268  nearest.setDistance(curHit.distance());
269 
270  // Assume at first that the near type is NONE on the
271  // triangle (i.e. on the face of the triangle) then it is
272  // therefore also for the face.
273 
274  nearType = NONE;
275 
276  if (tmpNearType == triPointRef::EDGE && tmpNearLabel == 0)
277  {
278  // If the triangle edge label is 0, then this is also
279  // an edge of the face, if not, it is on the face
280 
281  nearType = EDGE;
282 
283  nearLabel = pI;
284  }
285  else if (tmpNearType == triPointRef::POINT && tmpNearLabel < 2)
286  {
287  // If the triangle point label is 0 or 1, then this is
288  // also a point of the face, if not, it is on the face
289 
290  nearType = POINT;
291 
292  nearLabel = pI + tmpNearLabel;
293  }
294 
295  if (curHit.hit())
296  {
297  nearest.setHit();
298  nearest.setPoint(curHit.hitPoint());
299  }
300  else
301  {
302  // In nearest point, miss is always eligible
303  nearest.setMiss(true);
304  nearest.setPoint(curHit.missPoint());
305  }
306  }
307  }
308 
309  return nearest;
310 }
311 
312 
313 // ************************************************************************* //
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
const pointField & points
pointHit intersection(const point &p, const vector &q, const point &ctr, const pointField &, const intersection::algorithm alg, const scalar tol=0.0) const
Fast intersection with a ray.
pointHit ray(const point &p, const vector &n, const pointField &, const intersection::algorithm alg=intersection::FULL_RAY, const intersection::direction dir=intersection::VECTOR) const
Return potential intersection with face with a ray starting.
dimensioned< scalar > mag(const dimensioned< Type > &)
labelList f(nPoints)
const Point & missPoint() const
Return miss point.
Definition: PointHit.H:145
const Point & hitPoint() const
Return hit point.
Definition: PointHit.H:126
void setPoint(const Point &p)
Definition: PointHit.H:181
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
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:57
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
void setDistance(const scalar d)
Definition: PointHit.H:186
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
#define forAll(list, i)
Definition: UList.H:421
void setHit()
Definition: PointHit.H:169
scalar distance() const
Return distance to hit.
Definition: PointHit.H:139
void setMiss(const bool eligible)
Definition: PointHit.H:175
pointHit nearestPoint(const point &p, const pointField &) const
Return nearest point to face.
bool hit() const
Is there a hit.
Definition: PointHit.H:120
pointHit nearestPointClassify(const point &p, const pointField &, label &nearType, label &nearLabel) const
Return nearest point to face and classify it:
pointHit nearestPointClassify(const point &p, label &nearType, label &nearLabel) const
Find the nearest point to p on the triangle and classify it:
Definition: triangleI.H:508
label nPoints
bool eligibleMiss() const
Is this an eligible miss.
Definition: PointHit.H:164