triangleFuncs.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2022 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 "triangleFuncs.H"
27 #include "pointField.H"
28 #include "treeBoundBox.H"
29 #include "SortableList.H"
30 #include "boolList.H"
31 
32 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
33 
34 void Foam::triangleFuncs::setIntersection
35 (
36  const point& oppositeSidePt,
37  const scalar oppositeSign,
38 
39  const point& thisSidePt,
40  const scalar thisSign,
41 
42  const scalar tol,
43 
44  point& pt
45 )
46 {
47  scalar denom = oppositeSign - thisSign;
48 
49  if (mag(denom) < tol)
50  {
51  // If almost does not cut choose one which certainly cuts.
52  pt = oppositeSidePt;
53  }
54  else
55  {
56  pt = oppositeSidePt + oppositeSign/denom*(thisSidePt - oppositeSidePt);
57  }
58 }
59 
60 
61 void Foam::triangleFuncs::selectPt
62 (
63  const bool select0,
64  const point& p0,
65  const point& p1,
66  point& min
67 )
68 {
69  if (select0)
70  {
71  min = p0;
72  }
73  else
74  {
75  min = p1;
76  }
77 }
78 
79 
80 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
81 
82 // Intersect triangle with parallel edges aligned with axis i0.
83 // Returns true (and intersection in pInter) if any of them intersects triangle.
85 (
86  const point& V0,
87  const point& V10,
88  const point& V20,
89  const label i0,
90  const pointField& origin,
91  const scalar maxLength,
92  point& pInter
93 )
94 {
95  // Based on Graphics Gems - Fast Ray Triangle intersection.
96  // Since direction is coordinate axis there is no need to do projection,
97  // we can directly check u,v components for inclusion in triangle.
98 
99  // Get other components
100  const label i1 = (i0 + 1) % 3;
101  const label i2 = (i1 + 1) % 3;
102 
103  const scalar u1 = V10[i1];
104  const scalar v1 = V10[i2];
105 
106  const scalar u2 = V20[i1];
107  const scalar v2 = V20[i2];
108 
109  const scalar localScale = mag(u1) + mag(v1) + mag(u2) + mag(v2);
110 
111  const scalar det = v2*u1 - u2*v1;
112 
113  // Fix for V0:(-31.71428 0 -15.10714)
114  // V10:(-1.285715 8.99165e-16 -1.142858)
115  // V20:(0 0 -1.678573)
116  // i0:0
117  if (localScale < vSmall || Foam::mag(det)/localScale < small)
118  {
119  // Triangle parallel to dir
120  return false;
121  }
122 
123  forAll(origin, originI)
124  {
125  const point& P = origin[originI];
126 
127  const scalar u0 = P[i1] - V0[i1];
128  const scalar v0 = P[i2] - V0[i2];
129 
130  scalar alpha = 0;
131  scalar beta = 0;
132  bool inter = false;
133 
134  if (mag(u1) < rootVSmall)
135  {
136  beta = u0/u2;
137  if ((beta >= 0) && (beta <= 1))
138  {
139  alpha = (v0 - beta*v2)/v1;
140  inter = ((alpha >= 0) && ((alpha + beta) <= 1));
141  }
142  }
143  else
144  {
145  beta = (v0*u1 - u0*v1)/det;
146  if ((beta >= 0) && (beta <= 1))
147  {
148  alpha = (u0 - beta*u2)/u1;
149  inter = ((alpha >= 0) && ((alpha + beta) <= 1));
150  }
151  }
152 
153  if (inter)
154  {
155  pInter = V0 + alpha*V10 + beta*V20;
156  const scalar s = pInter[i0] - P[i0];
157 
158  if ((s >= 0) && (s <= maxLength))
159  {
160  return true;
161  }
162  }
163  }
164 
165  return false;
166 }
167 
168 
169 // Intersect triangle with bounding box. Return true if
170 // any of the faces of bb intersect triangle.
171 // Note: so returns false if triangle inside bb.
173 (
174  const point& p0,
175  const point& p1,
176  const point& p2,
177  const treeBoundBox& cubeBb
178 )
179 {
180  const vector p10 = p1 - p0;
181  const vector p20 = p2 - p0;
182 
183  // cubeBb points; counted as if cell with vertex0 at cubeBb.min().
184  const point& min = cubeBb.min();
185  const point& max = cubeBb.max();
186 
187  const point& cube0 = min;
188  const point cube1(min.x(), min.y(), max.z());
189  const point cube2(max.x(), min.y(), max.z());
190  const point cube3(max.x(), min.y(), min.z());
191 
192  const point cube4(min.x(), max.y(), min.z());
193  const point cube5(min.x(), max.y(), max.z());
194  const point cube7(max.x(), max.y(), min.z());
195 
196  //
197  // Intersect all 12 edges of cube with triangle
198  //
199 
200  point pInter;
201  pointField origin(4);
202  // edges in x direction
203  origin[0] = cube0;
204  origin[1] = cube1;
205  origin[2] = cube5;
206  origin[3] = cube4;
207 
208  scalar maxSx = max.x() - min.x();
209 
210  if (intersectAxesBundle(p0, p10, p20, 0, origin, maxSx, pInter))
211  {
212  return true;
213  }
214 
215  // edges in y direction
216  origin[0] = cube0;
217  origin[1] = cube1;
218  origin[2] = cube2;
219  origin[3] = cube3;
220 
221  scalar maxSy = max.y() - min.y();
222 
223  if (intersectAxesBundle(p0, p10, p20, 1, origin, maxSy, pInter))
224  {
225  return true;
226  }
227 
228  // edges in z direction
229  origin[0] = cube0;
230  origin[1] = cube3;
231  origin[2] = cube7;
232  origin[3] = cube4;
233 
234  scalar maxSz = max.z() - min.z();
235 
236  if (intersectAxesBundle(p0, p10, p20, 2, origin, maxSz, pInter))
237  {
238  return true;
239  }
240 
241 
242  // Intersect triangle edges with bounding box
243  if (cubeBb.intersects(p0, p1, pInter))
244  {
245  return true;
246  }
247  if (cubeBb.intersects(p1, p2, pInter))
248  {
249  return true;
250  }
251  if (cubeBb.intersects(p2, p0, pInter))
252  {
253  return true;
254  }
255 
256  return false;
257 }
258 
259 
263 //bool Foam::triangleFuncs::intersectBbExact
264 //(
265 // const point& p0,
266 // const point& p1,
267 // const point& p2,
268 // const treeBoundBox& cubeBb
269 //)
270 //{
271 // const point& min = cubeBb.min();
272 // const point& max = cubeBb.max();
273 //
274 // const point& cube0 = min;
275 // const point cube1(min.x(), min.y(), max.z());
276 // const point cube2(max.x(), min.y(), max.z());
277 // const point cube3(max.x(), min.y(), min.z());
278 //
279 // const point cube4(min.x(), max.y(), min.z());
280 // const point cube5(min.x(), max.y(), max.z());
281 // const point& cube6 = max;
282 // const point cube7(max.x(), max.y(), min.z());
283 //
284 // // Test intersection of triangle with twelve edges of box.
285 // {
286 // triPointRef tri(p0, p1, p2);
287 // if (tri.intersectionExact(cube0, cube1).hit())
288 // {
289 // return true;
290 // }
291 // if (tri.intersectionExact(cube1, cube2).hit())
292 // {
293 // return true;
294 // }
295 // if (tri.intersectionExact(cube2, cube3).hit())
296 // {
297 // return true;
298 // }
299 // if (tri.intersectionExact(cube3, cube0).hit())
300 // {
301 // return true;
302 // }
303 //
304 // if (tri.intersectionExact(cube4, cube5).hit())
305 // {
306 // return true;
307 // }
308 // if (tri.intersectionExact(cube5, cube6).hit())
309 // {
310 // return true;
311 // }
312 // if (tri.intersectionExact(cube6, cube7).hit())
313 // {
314 // return true;
315 // }
316 // if (tri.intersectionExact(cube7, cube4).hit())
317 // {
318 // return true;
319 // }
320 //
321 // if (tri.intersectionExact(cube0, cube4).hit())
322 // {
323 // return true;
324 // }
325 // if (tri.intersectionExact(cube1, cube5).hit())
326 // {
327 // return true;
328 // }
329 // if (tri.intersectionExact(cube2, cube6).hit())
330 // {
331 // return true;
332 // }
333 // if (tri.intersectionExact(cube3, cube7).hit())
334 // {
335 // return true;
336 // }
337 // }
338 // // Test intersection of triangle edges with bounding box
339 // {
340 // triPointRef tri(cube0, cube1, cube2);
341 // if (tri.intersectionExact(p0, p1).hit())
342 // {
343 // return true;
344 // }
345 // if (tri.intersectionExact(p1, p2).hit())
346 // {
347 // return true;
348 // }
349 // if (tri.intersectionExact(p2, p0).hit())
350 // {
351 // return true;
352 // }
353 // }
354 // {
355 // triPointRef tri(cube2, cube3, cube0);
356 // if (tri.intersectionExact(p0, p1).hit())
357 // {
358 // return true;
359 // }
360 // if (tri.intersectionExact(p1, p2).hit())
361 // {
362 // return true;
363 // }
364 // if (tri.intersectionExact(p2, p0).hit())
365 // {
366 // return true;
367 // }
368 // }
369 // {
370 // triPointRef tri(cube4, cube5, cube6);
371 // if (tri.intersectionExact(p0, p1).hit())
372 // {
373 // return true;
374 // }
375 // if (tri.intersectionExact(p1, p2).hit())
376 // {
377 // return true;
378 // }
379 // if (tri.intersectionExact(p2, p0).hit())
380 // {
381 // return true;
382 // }
383 // }
384 // {
385 // triPointRef tri(cube6, cube7, cube4);
386 // if (tri.intersectionExact(p0, p1).hit())
387 // {
388 // return true;
389 // }
390 // if (tri.intersectionExact(p1, p2).hit())
391 // {
392 // return true;
393 // }
394 // if (tri.intersectionExact(p2, p0).hit())
395 // {
396 // return true;
397 // }
398 // }
399 //
400 //
401 // {
402 // triPointRef tri(cube4, cube5, cube1);
403 // if (tri.intersectionExact(p0, p1).hit())
404 // {
405 // return true;
406 // }
407 // if (tri.intersectionExact(p1, p2).hit())
408 // {
409 // return true;
410 // }
411 // if (tri.intersectionExact(p2, p0).hit())
412 // {
413 // return true;
414 // }
415 // }
416 // {
417 // triPointRef tri(cube1, cube0, cube4);
418 // if (tri.intersectionExact(p0, p1).hit())
419 // {
420 // return true;
421 // }
422 // if (tri.intersectionExact(p1, p2).hit())
423 // {
424 // return true;
425 // }
426 // if (tri.intersectionExact(p2, p0).hit())
427 // {
428 // return true;
429 // }
430 // }
431 // {
432 // triPointRef tri(cube7, cube6, cube2);
433 // if (tri.intersectionExact(p0, p1).hit())
434 // {
435 // return true;
436 // }
437 // if (tri.intersectionExact(p1, p2).hit())
438 // {
439 // return true;
440 // }
441 // if (tri.intersectionExact(p2, p0).hit())
442 // {
443 // return true;
444 // }
445 // }
446 // {
447 // triPointRef tri(cube2, cube3, cube7);
448 // if (tri.intersectionExact(p0, p1).hit())
449 // {
450 // return true;
451 // }
452 // if (tri.intersectionExact(p1, p2).hit())
453 // {
454 // return true;
455 // }
456 // if (tri.intersectionExact(p2, p0).hit())
457 // {
458 // return true;
459 // }
460 // }
461 //
462 // {
463 // triPointRef tri(cube0, cube4, cube7);
464 // if (tri.intersectionExact(p0, p1).hit())
465 // {
466 // return true;
467 // }
468 // if (tri.intersectionExact(p1, p2).hit())
469 // {
470 // return true;
471 // }
472 // if (tri.intersectionExact(p2, p0).hit())
473 // {
474 // return true;
475 // }
476 // }
477 // {
478 // triPointRef tri(cube7, cube3, cube0);
479 // if (tri.intersectionExact(p0, p1).hit())
480 // {
481 // return true;
482 // }
483 // if (tri.intersectionExact(p1, p2).hit())
484 // {
485 // return true;
486 // }
487 // if (tri.intersectionExact(p2, p0).hit())
488 // {
489 // return true;
490 // }
491 // }
492 // {
493 // triPointRef tri(cube1, cube5, cube6);
494 // if (tri.intersectionExact(p0, p1).hit())
495 // {
496 // return true;
497 // }
498 // if (tri.intersectionExact(p1, p2).hit())
499 // {
500 // return true;
501 // }
502 // if (tri.intersectionExact(p2, p0).hit())
503 // {
504 // return true;
505 // }
506 // }
507 // {
508 // triPointRef tri(cube6, cube2, cube1);
509 // if (tri.intersectionExact(p0, p1).hit())
510 // {
511 // return true;
512 // }
513 // if (tri.intersectionExact(p1, p2).hit())
514 // {
515 // return true;
516 // }
517 // if (tri.intersectionExact(p2, p0).hit())
518 // {
519 // return true;
520 // }
521 // }
522 // return false;
523 //}
524 
525 
527 (
528  const point& va0,
529  const point& va10,
530  const point& va20,
531 
532  const point& base,
533  const point& normal,
534 
535  point& pInter0,
536  point& pInter1
537 )
538 {
539  // Get triangle normal
540  vector na = va10 ^ va20;
541  scalar magArea = mag(na);
542  na/magArea;
543 
544  if (mag(na & normal) > (1 - small))
545  {
546  // Parallel
547  return false;
548  }
549 
550  const point va1 = va0 + va10;
551  const point va2 = va0 + va20;
552 
553  // Find the triangle point on the other side.
554  scalar sign0 = (va0 - base) & normal;
555  scalar sign1 = (va1 - base) & normal;
556  scalar sign2 = (va2 - base) & normal;
557 
558  label oppositeVertex = -1;
559 
560  if (sign0 < 0)
561  {
562  if (sign1 < 0)
563  {
564  if (sign2 < 0)
565  {
566  // All on same side of plane
567  return false;
568  }
569  else // sign2 >= 0
570  {
571  // 2 on opposite side.
572  oppositeVertex = 2;
573  }
574  }
575  else // sign1 >= 0
576  {
577  if (sign2 < 0)
578  {
579  // 1 on opposite side.
580  oppositeVertex = 1;
581  }
582  else
583  {
584  // 0 on opposite side.
585  oppositeVertex = 0;
586  }
587  }
588  }
589  else // sign0 >= 0
590  {
591  if (sign1 < 0)
592  {
593  if (sign2 < 0)
594  {
595  // 0 on opposite side.
596  oppositeVertex = 0;
597  }
598  else // sign2 >= 0
599  {
600  // 1 on opposite side.
601  oppositeVertex = 1;
602  }
603  }
604  else // sign1 >= 0
605  {
606  if (sign2 < 0)
607  {
608  // 2 on opposite side.
609  oppositeVertex = 2;
610  }
611  else // sign2 >= 0
612  {
613  // All on same side of plane
614  return false;
615  }
616  }
617  }
618 
619  scalar tol = small*Foam::sqrt(magArea);
620 
621  if (oppositeVertex == 0)
622  {
623  // 0 on opposite side. Cut edges 01 and 02
624  setIntersection(va0, sign0, va1, sign1, tol, pInter0);
625  setIntersection(va0, sign0, va2, sign2, tol, pInter1);
626  }
627  else if (oppositeVertex == 1)
628  {
629  // 1 on opposite side. Cut edges 10 and 12
630  setIntersection(va1, sign1, va0, sign0, tol, pInter0);
631  setIntersection(va1, sign1, va2, sign2, tol, pInter1);
632  }
633  else // oppositeVertex == 2
634  {
635  // 2 on opposite side. Cut edges 20 and 21
636  setIntersection(va2, sign2, va0, sign0, tol, pInter0);
637  setIntersection(va2, sign2, va1, sign1, tol, pInter1);
638  }
639 
640  return true;
641 }
642 
643 
645 (
646  const point& va0,
647  const point& va10,
648  const point& va20,
649 
650  const point& vb0,
651  const point& vb10,
652  const point& vb20,
653 
654  point& pInter0,
655  point& pInter1
656 )
657 {
658  // Get triangle normals
659  vector na = va10 ^ va20;
660  na/mag(na);
661 
662  vector nb = vb10 ^ vb20;
663  nb/mag(nb);
664 
665  // Calculate intersection of triangle a with plane of b
666  point planeB0;
667  point planeB1;
668  if (!intersect(va0, va10, va20, vb0, nb, planeB0, planeB1))
669  {
670  return false;
671  }
672 
673  // ,, triangle b with plane of a
674  point planeA0;
675  point planeA1;
676  if (!intersect(vb0, vb10, vb20, va0, na, planeA0, planeA1))
677  {
678  return false;
679  }
680 
681  // Now check if intersections overlap (w.r.t. intersection of the two
682  // planes)
683 
684  vector intersection(na ^ nb);
685 
686  scalar coordB0 = planeB0 & intersection;
687  scalar coordB1 = planeB1 & intersection;
688 
689  scalar coordA0 = planeA0 & intersection;
690  scalar coordA1 = planeA1 & intersection;
691 
692  // Put information in indexable form for sorting.
693  List<point*> pts(4);
694  boolList isFromB(4);
695  SortableList<scalar> sortCoords(4);
696 
697  pts[0] = &planeB0;
698  isFromB[0] = true;
699  sortCoords[0] = coordB0;
700 
701  pts[1] = &planeB1;
702  isFromB[1] = true;
703  sortCoords[1] = coordB1;
704 
705  pts[2] = &planeA0;
706  isFromB[2] = false;
707  sortCoords[2] = coordA0;
708 
709  pts[3] = &planeA1;
710  isFromB[3] = false;
711  sortCoords[3] = coordA1;
712 
713  sortCoords.sort();
714 
715  const labelList& indices = sortCoords.indices();
716 
717  if (isFromB[indices[0]] == isFromB[indices[1]])
718  {
719  // Entry 0 and 1 are of same region (both a or both b). Hence that
720  // region does not overlap.
721  return false;
722  }
723  else
724  {
725  // Different type. Start of overlap at indices[1], end at indices[2]
726  pInter0 = *pts[indices[1]];
727  pInter1 = *pts[indices[2]];
728 
729  return true;
730  }
731 }
732 
733 
734 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:91
A list that is sorted upon construction or when explicitly requested with the sort() method.
Definition: SortableList.H:55
const labelList & indices() const
Return the list of sorted indices. Updated every sort.
Definition: SortableList.H:96
void sort()
(stable) sort the list (if changed after construction time)
Definition: SortableList.C:112
const point & min() const
Minimum point defining the bounding box.
Definition: boundBoxI.H:60
const point & max() const
Maximum point defining the bounding box.
Definition: boundBoxI.H:66
Foam::intersection.
Definition: intersection.H:50
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:90
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:229
static bool intersectBb(const point &p0, const point &p1, const point &p2, const treeBoundBox &cubeBb)
Does triangle intersect bounding box.
static bool intersect(const point &va0, const point &va10, const point &va20, const point &basePoint, const vector &normal, point &pInter0, point &pInter1)
Does triangle intersect plane. Return bool and set intersection segment.
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:85
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.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
volScalarField alpha(IOobject("alpha", runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
dimensionedScalar det(const dimensionedSphericalTensor &dt)
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
vector point
Point is a vector.
Definition: point.H:41
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< scalar > mag(const dimensioned< Type > &)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
scalar u0