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-------------------------------------------------------------------------------
11License
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
38void 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// ************************************************************************* //
Describes the interaction of a face and a point. It carries the info of a successful hit and (if succ...
Definition: PointHit.H:54
scalar distance() const noexcept
Return distance to hit.
Definition: PointHit.H:139
bool hit() const noexcept
Is there a hit.
Definition: PointHit.H:121
A list that is sorted upon construction or when explicitly requested with the sort() method.
Definition: SortableList.H:63
const labelList & indices() const noexcept
Return the list of sorted indices. Updated every sort.
Definition: SortableList.H:108
void sort()
Forward (stable) sort the list (if changed after construction).
Definition: SortableList.C:124
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
Foam::intersection.
Definition: intersection.H:53
Standard boundBox with extra functionality for use in octree.
Definition: treeBoundBox.H:89
static const edgeList edges
Edge to point addressing.
Definition: treeBoundBox.H:154
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
tmp< pointField > points() const
Vertex coordinates. In octant coding.
Definition: treeBoundBox.C:92
static bool intersectBb(const point &p0, const point &p1, const point &p2, const treeBoundBox &cubeBb)
Intersect triangle with bounding box.
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.
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
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:80
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
const volScalarField & p0
Definition: EEqn.H:36
const pointField & points
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))
dimensionedScalar det(const dimensionedSphericalTensor &dt)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
scalar u0
volScalarField & alpha
volScalarField & e
Definition: createFields.H:11
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333