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