triangleFuncs.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011 OpenFOAM Foundation
9  Copyright (C) 2017 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "triangleFuncs.H"
30 #include "pointField.H"
31 #include "treeBoundBox.H"
32 #include "SortableList.H"
33 #include "boolList.H"
34 #include "triPointRef.H"
35 
36 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
37 
38 void Foam::triangleFuncs::setIntersection
39 (
40  const point& oppositeSidePt,
41  const scalar oppositeSign,
42 
43  const point& thisSidePt,
44  const scalar thisSign,
45 
46  const scalar tol,
47 
48  point& pt
49 )
50 {
51  const scalar denom = oppositeSign - thisSign;
52 
53  if (mag(denom) < tol)
54  {
55  // If almost does not cut choose one which certainly cuts.
56  pt = oppositeSidePt;
57  }
58  else
59  {
60  pt = oppositeSidePt + oppositeSign/denom*(thisSidePt - oppositeSidePt);
61  }
62 }
63 
64 
65 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
66 
68 (
69  const point& V0,
70  const point& V10,
71  const point& V20,
72  const label i0,
73  const pointField& origin,
74  const scalar maxLength,
75  point& pInter
76 )
77 {
78  // Based on Graphics Gems - Fast Ray Triangle intersection.
79  // Since direction is coordinate axis there is no need to do projection,
80  // we can directly check u,v components for inclusion in triangle.
81 
82  // Get other components
83  label i1 = (i0 + 1) % 3;
84  label i2 = (i1 + 1) % 3;
85 
86  scalar u1 = V10[i1];
87  scalar v1 = V10[i2];
88 
89  scalar u2 = V20[i1];
90  scalar v2 = V20[i2];
91 
92  scalar localScale = mag(u1)+mag(v1)+mag(u2)+mag(v2);
93 
94  scalar det = v2*u1 - u2*v1;
95 
96  // Fix for V0:(-31.71428 0 -15.10714)
97  // V10:(-1.285715 8.99165e-16 -1.142858)
98  // V20:(0 0 -1.678573)
99  // i0:0
100  if (localScale < VSMALL || Foam::mag(det)/localScale < SMALL)
101  {
102  // Triangle parallel to dir
103  return false;
104  }
105 
106  forAll(origin, originI)
107  {
108  const point& P = origin[originI];
109 
110  scalar u0 = P[i1] - V0[i1];
111  scalar v0 = P[i2] - V0[i2];
112 
113  scalar alpha = 0;
114  scalar beta = 0;
115  bool inter = false;
116 
117  if (Foam::mag(u1) < ROOTVSMALL)
118  {
119  beta = u0/u2;
120  if ((beta >= 0) && (beta <= 1))
121  {
122  alpha = (v0 - beta*v2)/v1;
123  inter = ((alpha >= 0) && ((alpha + beta) <= 1));
124  }
125  }
126  else
127  {
128  beta = (v0*u1 - u0*v1)/det;
129  if ((beta >= 0) && (beta <= 1))
130  {
131  alpha = (u0 - beta*u2)/u1;
132  inter = ((alpha >= 0) && ((alpha + beta) <= 1));
133  }
134  }
135 
136  if (inter)
137  {
138  pInter = V0 + alpha*V10 + beta*V20;
139  scalar s = (pInter - origin[originI])[i0];
140 
141  if ((s >= 0) && (s <= maxLength))
142  {
143  return true;
144  }
145  }
146  }
147  return false;
148 }
149 
150 
152 (
153  const point& p0,
154  const point& p1,
155  const point& p2,
156  const treeBoundBox& cubeBb
157 )
158 {
159  // Slow (edge by edge) bounding box intersection. TBD: replace with call
160  // to above intersectAxesBundle. However this function is not fully
161  // correct and misses intersection between some triangles.
162  {
163  const triPointRef tri(p0, p1, p2);
164 
165  const edgeList& es = treeBoundBox::edges;
166  const pointField points(cubeBb.points());
167 
168  forAll(es, i)
169  {
170  const edge& e = es[i];
171  const point& start = points[e[0]];
172  const point& end = points[e[1]];
173 
174  pointHit inter = tri.intersection
175  (
176  start,
177  end-start,
179  );
180 
181  if (inter.hit() && inter.distance() <= 1)
182  {
183  return true;
184  }
185  }
186  }
187 
188 
189  // Intersect triangle edges with bounding box
190  point pInter;
191  if (cubeBb.intersects(p0, p1, pInter))
192  {
193  return true;
194  }
195  if (cubeBb.intersects(p1, p2, pInter))
196  {
197  return true;
198  }
199  if (cubeBb.intersects(p2, p0, pInter))
200  {
201  return true;
202  }
203 
204  return false;
205 }
206 
207 
209 (
210  const point& va0,
211  const point& va10,
212  const point& va20,
213 
214  const point& base,
215  const point& normal,
216 
217  point& pInter0,
218  point& pInter1
219 )
220 {
221  // Get triangle normal
222  vector na = va10 ^ va20;
223  scalar magArea = mag(na);
224  na/magArea;
225 
226  if (mag(na & normal) > (1 - SMALL))
227  {
228  // Parallel
229  return false;
230  }
231 
232  const point va1 = va0 + va10;
233  const point va2 = va0 + va20;
234 
235  // Find the triangle point on the other side.
236  scalar sign0 = (va0 - base) & normal;
237  scalar sign1 = (va1 - base) & normal;
238  scalar sign2 = (va2 - base) & normal;
239 
240  label oppositeVertex = -1;
241 
242  if (sign0 < 0)
243  {
244  if (sign1 < 0)
245  {
246  if (sign2 < 0)
247  {
248  // All on same side of plane
249  return false;
250  }
251  else // sign2 >= 0
252  {
253  // 2 on opposite side.
254  oppositeVertex = 2;
255  }
256  }
257  else // sign1 >= 0
258  {
259  if (sign2 < 0)
260  {
261  // 1 on opposite side.
262  oppositeVertex = 1;
263  }
264  else
265  {
266  // 0 on opposite side.
267  oppositeVertex = 0;
268  }
269  }
270  }
271  else // sign0 >= 0
272  {
273  if (sign1 < 0)
274  {
275  if (sign2 < 0)
276  {
277  // 0 on opposite side.
278  oppositeVertex = 0;
279  }
280  else // sign2 >= 0
281  {
282  // 1 on opposite side.
283  oppositeVertex = 1;
284  }
285  }
286  else // sign1 >= 0
287  {
288  if (sign2 < 0)
289  {
290  // 2 on opposite side.
291  oppositeVertex = 2;
292  }
293  else // sign2 >= 0
294  {
295  // All on same side of plane
296  return false;
297  }
298  }
299  }
300 
301  scalar tol = SMALL*Foam::sqrt(magArea);
302 
303  if (oppositeVertex == 0)
304  {
305  // 0 on opposite side. Cut edges 01 and 02
306  setIntersection(va0, sign0, va1, sign1, tol, pInter0);
307  setIntersection(va0, sign0, va2, sign2, tol, pInter1);
308  }
309  else if (oppositeVertex == 1)
310  {
311  // 1 on opposite side. Cut edges 10 and 12
312  setIntersection(va1, sign1, va0, sign0, tol, pInter0);
313  setIntersection(va1, sign1, va2, sign2, tol, pInter1);
314  }
315  else // oppositeVertex == 2
316  {
317  // 2 on opposite side. Cut edges 20 and 21
318  setIntersection(va2, sign2, va0, sign0, tol, pInter0);
319  setIntersection(va2, sign2, va1, sign1, tol, pInter1);
320  }
321 
322  return true;
323 }
324 
325 
327 (
328  const point& va0,
329  const point& va10,
330  const point& va20,
331 
332  const point& vb0,
333  const point& vb10,
334  const point& vb20,
335 
336  point& pInter0,
337  point& pInter1
338 )
339 {
340  // Get triangle normals
341  vector na = va10 ^ va20;
342  na/mag(na);
343 
344  vector nb = vb10 ^ vb20;
345  nb/mag(nb);
346 
347  // Calculate intersection of triangle a with plane of b
348  point planeB0;
349  point planeB1;
350  if (!intersect(va0, va10, va20, vb0, nb, planeB0, planeB1))
351  {
352  return false;
353  }
354 
355  // ,, triangle b with plane of a
356  point planeA0;
357  point planeA1;
358  if (!intersect(vb0, vb10, vb20, va0, na, planeA0, planeA1))
359  {
360  return false;
361  }
362 
363  // Now check if intersections overlap (w.r.t. intersection of the two
364  // planes)
365 
366  vector intersection(na ^ nb);
367 
368  scalar coordB0 = planeB0 & intersection;
369  scalar coordB1 = planeB1 & intersection;
370 
371  scalar coordA0 = planeA0 & intersection;
372  scalar coordA1 = planeA1 & intersection;
373 
374  // Put information in indexable form for sorting.
375  List<point*> pts(4);
376  boolList isFromB(4);
377  SortableList<scalar> sortCoords(4);
378 
379  pts[0] = &planeB0;
380  isFromB[0] = true;
381  sortCoords[0] = coordB0;
382 
383  pts[1] = &planeB1;
384  isFromB[1] = true;
385  sortCoords[1] = coordB1;
386 
387  pts[2] = &planeA0;
388  isFromB[2] = false;
389  sortCoords[2] = coordA0;
390 
391  pts[3] = &planeA1;
392  isFromB[3] = false;
393  sortCoords[3] = coordA1;
394 
395  sortCoords.sort();
396 
397  const labelList& indices = sortCoords.indices();
398 
399  if (isFromB[indices[0]] == isFromB[indices[1]])
400  {
401  // Entry 0 and 1 are of same region (both a or both b). Hence that
402  // region does not overlap.
403  return false;
404  }
405  else
406  {
407  // Different type. Start of overlap at indices[1], end at indices[2]
408  pInter0 = *pts[indices[1]];
409  pInter1 = *pts[indices[2]];
410 
411  return true;
412  }
413 }
414 
415 
416 // ************************************************************************* //
boolList.H
Foam::intersection::HALF_RAY
Definition: intersection.H:75
s
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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))
Definition: gmvOutputSpray.H:25
Foam::SortableList::sort
void sort()
Forward (stable) sort the list (if changed after construction).
Definition: SortableList.C:124
Foam::treeBoundBox::points
tmp< pointField > points() const
Vertex coordinates. In octant coding.
Definition: treeBoundBox.C:92
Foam::PointHit
Describes the interaction of a face and a point. It carries the info of a successful hit and (if succ...
Definition: PointHit.H:53
Foam::triangleFuncs::intersectAxesBundle
static bool intersectAxesBundle(const point &V0, const point &V10, const point &V20, const label i0, const pointField &origin, const scalar maxLength, point &pInter)
Intersect triangle with parallel edges aligned with axis i0.
Definition: triangleFuncs.C:68
Foam::constant::atomic::alpha
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
Definition: readThermalProperties.H:212
Foam::treeBoundBox
Standard boundBox with extra functionality for use in octree.
Definition: treeBoundBox.H:86
Foam::edge
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:63
Foam::PointHit::distance
scalar distance() const noexcept
Return distance to hit.
Definition: PointHit.H:139
Foam::treeBoundBox::edges
static const edgeList edges
Edge to point addressing.
Definition: treeBoundBox.H:154
triPointRef.H
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::triangleFuncs::intersectBb
static bool intersectBb(const point &p0, const point &p1, const point &p2, const treeBoundBox &cubeBb)
Intersect triangle with bounding box.
Definition: triangleFuncs.C:152
Foam::triangle
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:59
SortableList.H
Foam::Field< vector >
treeBoundBox.H
Foam::treeBoundBox::intersects
bool intersects(const point &overallStart, const vector &overallVec, const point &start, const point &end, point &pt, direction &ptBits) const
Intersects segment; set point to intersection position and face,.
Definition: treeBoundBox.C:162
beta
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)
Foam::SortableList
A list that is sorted upon construction or when explicitly requested with the sort() method.
Definition: List.H:60
stdFoam::end
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:121
triangleFuncs.H
pointField.H
u0
scalar u0
Definition: readInitialConditions.H:90
Foam::Vector< scalar >
Foam::List< edge >
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::intersection
Foam::intersection.
Definition: intersection.H:52
Foam::triangle::intersection
pointHit intersection(const point &p, const vector &q, const intersection::algorithm alg, const scalar tol=0.0) const
Fast intersection with a ray.
Definition: triangleI.H:439
Foam::SortableList::indices
const labelList & indices() const
Return the list of sorted indices. Updated every sort.
Definition: SortableList.H:108
Foam::det
dimensionedScalar det(const dimensionedSphericalTensor &dt)
Definition: dimensionedSphericalTensor.C:62
Foam::triangleFuncs::intersect
static bool intersect(const point &va0, const point &va10, const point &va20, const point &basePoint, const vector &normal, point &pInter0, point &pInter1)
Intersect triangle with plane.
Definition: triangleFuncs.C:209
p0
const volScalarField & p0
Definition: EEqn.H:36
Foam::point
vector point
Point is a vector.
Definition: point.H:43
Foam::PointHit::hit
bool hit() const noexcept
Is there a hit.
Definition: PointHit.H:121