FreeFOAM The Cross-Platform CFD Toolkit
triangleI.H
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) 1991-2010 OpenCFD Ltd.
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 <OpenFOAM/IOstreams.H>
27 #include <OpenFOAM/pointHit.H>
29 
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 
35 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
36 
37 template<class Point, class PointRef>
38 pointHit triangle<Point, PointRef>::nearestPoint
39 (
40  const Point& baseVertex,
41  const vector& E0,
42  const vector& E1,
43  const point& P
44 )
45 {
46  // Distance vector
47  const vector D(baseVertex - P);
48 
49  // Some geometrical factors
50  const scalar a = E0 & E0;
51  const scalar b = E0 & E1;
52  const scalar c = E1 & E1;
53 
54  // Precalculate distance factors
55  const scalar d = E0 & D;
56  const scalar e = E1 & D;
57  const scalar f = D & D;
58 
59  // Do classification
60  const scalar det = a*c - b*b;
61  scalar s = b*e - c*d;
62  scalar t = b*d - a*e;
63 
64  bool inside = false;
65 
66  if (s+t < det)
67  {
68  if (s < 0)
69  {
70  if (t < 0)
71  {
72  // Region 4
73  if (e > 0)
74  {
75  // min on edge t = 0
76  t = 0;
77  s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
78  }
79  else
80  {
81  // min on edge s=0
82  s = 0;
83  t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
84  }
85  }
86  else
87  {
88  // Region 3. Min on edge s = 0
89  s = 0;
90  t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
91  }
92  }
93  else if (t < 0)
94  {
95  // Region 5
96  t = 0;
97  s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
98  }
99  else
100  {
101  // Region 0
102  const scalar invDet = 1/det;
103  s *= invDet;
104  t *= invDet;
105 
106  inside = true;
107  }
108  }
109  else
110  {
111  if (s < 0)
112  {
113  // Region 2
114  const scalar tmp0 = b + d;
115  const scalar tmp1 = c + e;
116  if (tmp1 > tmp0)
117  {
118  // min on edge s+t=1
119  const scalar numer = tmp1 - tmp0;
120  const scalar denom = a-2*b+c;
121  s = (numer >= denom ? 1 : numer/denom);
122  t = 1 - s;
123  }
124  else
125  {
126  // min on edge s=0
127  s = 0;
128  t = (tmp1 <= 0 ? 1 : (e >= 0 ? 0 : - e/c));
129  }
130  }
131  else if (t < 0)
132  {
133  // Region 6
134  const scalar tmp0 = b + d;
135  const scalar tmp1 = c + e;
136  if (tmp1 > tmp0)
137  {
138  // min on edge s+t=1
139  const scalar numer = tmp1 - tmp0;
140  const scalar denom = a-2*b+c;
141  s = (numer >= denom ? 1 : numer/denom);
142  t = 1 - s;
143  }
144  else
145  {
146  // min on edge t=0
147  t = 0;
148  s = (tmp1 <= 0 ? 1 : (d >= 0 ? 0 : - d/a));
149  }
150  }
151  else
152  {
153  // Region 1
154  const scalar numer = c+e-(b+d);
155  if (numer <= 0)
156  {
157  s = 0;
158  }
159  else
160  {
161  const scalar denom = a-2*b+c;
162  s = (numer >= denom ? 1 : numer/denom);
163  }
164  }
165 
166  t = 1 - s;
167  }
168 
169  // Calculate distance.
170  // Note: Foam::mag used since truncation error causes negative distances
171  // with points very close to one of the triangle vertices.
172  // (Up to -2.77556e-14 seen). Could use +SMALL but that not large enough.
173 
174  return pointHit
175  (
176  inside,
177  baseVertex + s*E0 + t*E1,
178  Foam::sqrt
179  (
180  Foam::mag(a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f)
181  ),
182  !inside
183  );
184 }
185 
186 
187 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
188 
189 template<class Point, class PointRef>
191 (
192  const Point& a,
193  const Point& b,
194  const Point& c
195 )
196 :
197  a_(a),
198  b_(b),
199  c_(c)
200 {}
201 
202 
203 template<class Point, class PointRef>
205 {
206  // Read beginning of triangle point pair
207  is.readBegin("triangle");
208 
209  is >> a_ >> b_ >> c_;
210 
211  // Read end of triangle point pair
212  is.readEnd("triangle");
213 
214  // Check state of Istream
215  is.check("triangle::triangle(Istream& is)");
216 }
217 
218 
219 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
220 
221 template<class Point, class PointRef>
222 inline const Point& triangle<Point, PointRef>::a() const
223 {
224  return a_;
225 }
226 
227 template<class Point, class PointRef>
228 inline const Point& triangle<Point, PointRef>::b() const
229 {
230  return b_;
231 }
232 
233 template<class Point, class PointRef>
234 inline const Point& triangle<Point, PointRef>::c() const
235 {
236  return c_;
237 }
238 
239 
240 template<class Point, class PointRef>
242 {
243  return (1.0/3.0)*(a_ + b_ + c_);
244 }
245 
246 
247 template<class Point, class PointRef>
248 inline scalar triangle<Point, PointRef>::mag() const
249 {
251 }
252 
253 
254 template<class Point, class PointRef>
256 {
257  return 0.5*((b_ - a_)^(c_ - a_));
258 }
259 
260 
261 template<class Point, class PointRef>
263 {
264  scalar d1 = (c_ - a_)&(b_ - a_);
265  scalar d2 = -(c_ - b_)&(b_ - a_);
266  scalar d3 = (c_ - a_)&(c_ - b_);
267 
268  scalar c1 = d2*d3;
269  scalar c2 = d3*d1;
270  scalar c3 = d1*d2;
271 
272  scalar c = c1 + c2 + c3;
273 
274  return
275  (
276  ((c2 + c3)*a_ + (c3 + c1)*b_ + (c1 + c2)*c_)/(2*c)
277  );
278 }
279 
280 
281 template<class Point, class PointRef>
283 {
284  scalar d1 = (c_ - a_) & (b_ - a_);
285  scalar d2 = - (c_ - b_) & (b_ - a_);
286  scalar d3 = (c_ - a_) & (c_ - b_);
287 
288  scalar denom = d2*d3 + d3*d1 + d1*d2;
289 
290  if (Foam::mag(denom) < VSMALL)
291  {
292  return GREAT;
293  }
294  else
295  {
296  scalar a = (d1 + d2)*(d2 + d3)*(d3 + d1) / denom;
297 
298  return 0.5*Foam::sqrt(min(GREAT, max(0, a)));
299  }
300 }
301 
302 
303 template<class Point, class PointRef>
305 {
306  return
307  mag()
308  / (
310  * Foam::sqr(circumRadius())
311  * 0.413497
312  + VSMALL
313  );
314 }
315 
316 
317 template<class Point, class PointRef>
318 inline scalar triangle<Point, PointRef>::sweptVol(const triangle& t) const
319 {
320  return (1.0/12.0)*
321  (
322  ((t.a_ - a_) & ((b_ - a_)^(c_ - a_)))
323  + ((t.b_ - b_) & ((c_ - b_)^(t.a_ - b_)))
324  + ((c_ - t.c_) & ((t.b_ - t.c_)^(t.a_ - t.c_)))
325 
326  + ((t.a_ - a_) & ((b_ - a_)^(c_ - a_)))
327  + ((b_ - t.b_) & ((t.a_ - t.b_)^(t.c_ - t.b_)))
328  + ((c_ - t.c_) & ((b_ - t.c_)^(t.a_ - t.c_)))
329  );
330 }
331 
332 
333 template<class Point, class PointRef>
335 (
336  const point& p,
337  const vector& q,
338  const intersection::algorithm alg,
339  const intersection::direction dir
340 ) const
341 {
342  // Express triangle in terms of baseVertex (point a_) and
343  // two edge vectors
344  vector E0 = b_ - a_;
345  vector E1 = c_ - a_;
346 
347  // Initialize intersection to miss.
348  pointHit inter(p);
349 
350  vector n(0.5*(E0 ^ E1));
351  scalar area = Foam::mag(n);
352 
353  if (area < VSMALL)
354  {
355  // Ineligible miss.
356  inter.setMiss(false);
357 
358  // The miss point is the nearest point on the triangle. Take any one.
359  inter.setPoint(a_);
360 
361  // The distance to the miss is the distance between the
362  // original point and plane of intersection. No normal so use
363  // distance from p to a_
364  inter.setDistance(Foam::mag(a_ - p));
365 
366  return inter;
367  }
368 
369  vector q1 = q/Foam::mag(q);
370 
371  if (dir == intersection::CONTACT_SPHERE)
372  {
373  n /= area;
374 
375  return ray(p, q1 - n, alg, intersection::VECTOR);
376  }
377 
378  // Intersection point with triangle plane
379  point pInter;
380  // Is intersection point inside triangle
381  bool hit;
382  {
383  // Reuse the fast ray intersection routine below in FULL_RAY
384  // mode since the original intersection routine has rounding problems.
385  pointHit fastInter = intersection(p, q1, intersection::FULL_RAY);
386  hit = fastInter.hit();
387 
388  if (hit)
389  {
390  pInter = fastInter.rawPoint();
391  }
392  else
393  {
394  // Calculate intersection of ray with triangle plane
395  vector v = a_ - p;
396  pInter = p + (q1&v)*q1;
397  }
398  }
399 
400  // Distance to intersection point
401  scalar dist = q1 & (pInter - p);
402 
403  const scalar planarPointTol =
404  Foam::min
405  (
406  Foam::min
407  (
408  Foam::mag(E0),
409  Foam::mag(E1)
410  ),
411  Foam::mag(c_ - b_)
413 
414  bool eligible =
416  || (alg == intersection::HALF_RAY && dist > -planarPointTol)
417  || (
418  alg == intersection::VISIBLE
419  && ((q1 & normal()) < -VSMALL)
420  );
421 
422  if (hit && eligible)
423  {
424  // Hit. Set distance to intersection.
425  inter.setHit();
426  inter.setPoint(pInter);
427  inter.setDistance(dist);
428  }
429  else
430  {
431  // Miss or ineligible hit.
432  inter.setMiss(eligible);
433 
434  // The miss point is the nearest point on the triangle
435  inter.setPoint(nearestPoint(a_, E0, E1, p).rawPoint());
436 
437  // The distance to the miss is the distance between the
438  // original point and plane of intersection
439  inter.setDistance(Foam::mag(pInter - p));
440  }
441 
442 
443  return inter;
444 }
445 
446 
447 // From "Fast, Minimum Storage Ray/Triangle Intersection"
448 // Moeller/Trumbore.
449 template<class Point, class PointRef>
451 (
452  const point& orig,
453  const vector& dir,
454  const intersection::algorithm alg,
455  const scalar tol
456 ) const
457 {
458  const vector edge1 = b_ - a_;
459  const vector edge2 = c_ - a_;
460 
461  // begin calculating determinant - also used to calculate U parameter
462  const vector pVec = dir ^ edge2;
463 
464  // if determinant is near zero, ray lies in plane of triangle
465  const scalar det = edge1 & pVec;
466 
467  // Initialise to miss
468  pointHit intersection(false, vector::zero, GREAT, false);
469 
470  if (alg == intersection::VISIBLE)
471  {
472  // Culling branch
473  if (det < ROOTVSMALL)
474  {
475  // Ray on wrong side of triangle. Return miss
476  return intersection;
477  }
478  }
479  else if (alg == intersection::HALF_RAY || alg == intersection::FULL_RAY)
480  {
481  // Non-culling branch
482  if (det > -ROOTVSMALL && det < ROOTVSMALL)
483  {
484  // Ray parallel to triangle. Return miss
485  return intersection;
486  }
487  }
488 
489  const scalar inv_det = 1.0 / det;
490 
491  /* calculate distance from a_ to ray origin */
492  const vector tVec = orig-a_;
493 
494  /* calculate U parameter and test bounds */
495  const scalar u = (tVec & pVec)*inv_det;
496 
497  if (u < -tol || u > 1.0+tol)
498  {
499  // return miss
500  return intersection;
501  }
502 
503  /* prepare to test V parameter */
504  const vector qVec = tVec ^ edge1;
505 
506  /* calculate V parameter and test bounds */
507  const scalar v = (dir & qVec) * inv_det;
508 
509  if (v < -tol || u + v > 1.0+tol)
510  {
511  // return miss
512  return intersection;
513  }
514 
515  /* calculate t, scale parameters, ray intersects triangle */
516  const scalar t = (edge2 & qVec) * inv_det;
517 
518  if (alg == intersection::HALF_RAY && t < -tol)
519  {
520  // Wrong side of orig. Return miss
521  return intersection;
522  }
523 
524  intersection.setHit();
525  intersection.setPoint(a_ + u*edge1 + v*edge2);
526  intersection.setDistance(t);
527 
528  return intersection;
529 }
530 
531 
532 template<class Point, class PointRef>
534 (
535  const point& p
536 ) const
537 {
538  // Express triangle in terms of baseVertex (point a_) and
539  // two edge vectors
540  vector E0 = b_ - a_;
541  vector E1 = c_ - a_;
542 
543  return nearestPoint(a_, E0, E1, p);
544 }
545 
546 
547 template<class Point, class PointRef>
549 (
550  const point& p,
551  const scalar tol,
552  label& nearType,
553  label& nearLabel
554 ) const
555 {
556  const vector E0 = b_ - a_;
557  const vector E1 = c_ - a_;
558  const vector n = 0.5*(E0 ^ E1);
559 
560  // Get largest component of normal
561  scalar magX = Foam::mag(n.x());
562  scalar magY = Foam::mag(n.y());
563  scalar magZ = Foam::mag(n.z());
564 
565  label i0 = -1;
566  if ((magX >= magY) && (magX >= magZ))
567  {
568  i0 = 0;
569  }
570  else if ((magY >= magX) && (magY >= magZ))
571  {
572  i0 = 1;
573  }
574  else
575  {
576  i0 = 2;
577  }
578 
579  // Get other components
580  label i1 = (i0 + 1) % 3;
581  label i2 = (i1 + 1) % 3;
582 
583 
584  scalar u1 = E0[i1];
585  scalar v1 = E0[i2];
586 
587  scalar u2 = E1[i1];
588  scalar v2 = E1[i2];
589 
590  scalar det = v2*u1 - u2*v1;
591 
592  scalar u0 = p[i1] - a_[i1];
593  scalar v0 = p[i2] - a_[i2];
594 
595  scalar alpha = 0;
596  scalar beta = 0;
597 
598  bool hit = false;
599 
600  if (Foam::mag(u1) < ROOTVSMALL)
601  {
602  beta = u0/u2;
603 
604  alpha = (v0 - beta*v2)/v1;
605 
606  hit = ((alpha >= 0) && ((alpha + beta) <= 1));
607  }
608  else
609  {
610  beta = (v0*u1 - u0*v1)/det;
611 
612  alpha = (u0 - beta*u2)/u1;
613 
614  hit = ((alpha >= 0) && ((alpha + beta) <= 1));
615  }
616 
617  //
618  // Now alpha, beta are the coordinates in the triangle local coordinate
619  // system E0, E1
620  //
621 
622  //Pout<< "alpha:" << alpha << endl;
623  //Pout<< "beta:" << beta << endl;
624  //Pout<< "hit:" << hit << endl;
625  //Pout<< "tol:" << tol << endl;
626 
627  if (hit)
628  {
629  // alpha,beta might get negative due to precision errors
630  alpha = max(0.0, min(1.0, alpha));
631  beta = max(0.0, min(1.0, beta));
632  }
633 
634  nearType = NONE;
635  nearLabel = -1;
636 
637  if (Foam::mag(alpha+beta-1) <= tol)
638  {
639  // On edge between vert 1-2 (=edge 1)
640 
641  if (Foam::mag(alpha) <= tol)
642  {
643  nearType = POINT;
644  nearLabel = 2;
645  }
646  else if (Foam::mag(beta) <= tol)
647  {
648  nearType = POINT;
649  nearLabel = 1;
650  }
651  else if ((alpha >= 0) && (alpha <= 1) && (beta >= 0) && (beta <= 1))
652  {
653  nearType = EDGE;
654  nearLabel = 1;
655  }
656  }
657  else if (Foam::mag(alpha) <= tol)
658  {
659  // On edge between vert 2-0 (=edge 2)
660 
661  if (Foam::mag(beta) <= tol)
662  {
663  nearType = POINT;
664  nearLabel = 0;
665  }
666  else if (Foam::mag(beta-1) <= tol)
667  {
668  nearType = POINT;
669  nearLabel = 2;
670  }
671  else if ((beta >= 0) && (beta <= 1))
672  {
673  nearType = EDGE;
674  nearLabel = 2;
675  }
676  }
677  else if (Foam::mag(beta) <= tol)
678  {
679  // On edge between vert 0-1 (= edge 0)
680 
681  if (Foam::mag(alpha) <= tol)
682  {
683  nearType = POINT;
684  nearLabel = 0;
685  }
686  else if (Foam::mag(alpha-1) <= tol)
687  {
688  nearType = POINT;
689  nearLabel = 1;
690  }
691  else if ((alpha >= 0) && (alpha <= 1))
692  {
693  nearType = EDGE;
694  nearLabel = 0;
695  }
696  }
697 
698  return hit;
699 }
700 
701 
702 
703 
704 
705 // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
706 
707 template<class point, class pointRef>
709 {
710  // Read beginning of triangle point pair
711  is.readBegin("triangle");
712 
713  is >> t.a_ >> t.b_ >> t.c_;
714 
715  // Read end of triangle point pair
716  is.readEnd("triangle");
717 
718  // Check state of Ostream
719  is.check("Istream& operator>>(Istream&, triangle&)");
720 
721  return is;
722 }
723 
724 
725 template<class Point, class PointRef>
726 inline Ostream& operator<<(Ostream& os, const triangle<Point, PointRef>& t)
727 {
728  os << nl
730  << t.a_ << token::SPACE << t.b_ << token::SPACE << t.c_
731  << token::END_LIST;
732 
733  return os;
734 }
735 
736 
737 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
738 
739 } // End namespace Foam
740 
741 // ************************ vim: set sw=4 sts=4 et: ************************ //
742