Geant4  10.01.p03
UTriangularFacet.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * This Software is part of the AIDA Unified Solids Library package *
4 // * See: https://aidasoft.web.cern.ch/USolids *
5 // ********************************************************************
6 //
7 // $Id:$
8 //
9 // --------------------------------------------------------------------
10 //
11 // UTriangularFacet
12 //
13 // 22.08.12 Marek Gayer
14 // Created from original implementation in Geant4
15 // --------------------------------------------------------------------
16 
17 #include <sstream>
18 #include <algorithm>
19 
20 #include "UVector2.hh"
21 #include "UTessellatedSolid.hh"
22 #include "VUSolid.hh"
23 #include "UUtils.hh"
24 #include "UTriangularFacet.hh"
25 
26 using namespace std;
27 
29 //
30 // Definition of triangular facet using absolute vectors to vertices.
31 // From this for first vector is retained to define the facet location and
32 // two relative vectors (E0 and E1) define the sides and orientation of
33 // the outward surface normal.
34 //
35 UTriangularFacet::UTriangularFacet(const UVector3& vt0, const UVector3& vt1, const UVector3& vt2, UFacetVertexType vertexType)
36  : fSqrDist(0.)
37 {
38  fVertices = new vector<UVector3>(3);
39 
40  SetVertex(0, vt0);
41  if (vertexType == UABSOLUTE)
42  {
43  SetVertex(1, vt1);
44  SetVertex(2, vt2);
45  fE1 = vt1 - vt0;
46  fE2 = vt2 - vt0;
47  }
48  else
49  {
50  SetVertex(1, vt0 + vt1);
51  SetVertex(2, vt0 + vt2);
52  fE1 = vt1;
53  fE2 = vt2;
54  }
55 
56  for (int i = 0; i < 3; ++i) fIndices[i] = -1;
57 
58  double eMag1 = fE1.Mag();
59  double eMag2 = fE2.Mag();
60  double eMag3 = (fE2 - fE1).Mag();
61 
62  if (eMag1 <= kCarTolerance || eMag2 <= kCarTolerance || eMag3 <= kCarTolerance)
63  {
64  ostringstream message;
65  message << "Length of sides of facet are too small." << endl
66  << "P[0] = " << GetVertex(0) << endl
67  << "P[1] = " << GetVertex(1) << endl
68  << "P[2] = " << GetVertex(2) << endl
69  << "Side lengths = P[0]->P[1]" << eMag1 << endl
70  << "Side lengths = P[0]->P[2]" << eMag2 << endl
71  << "Side lengths = P[1]->P[2]" << eMag3;
72  UUtils::Exception("UTriangularFacet::UTriangularFacet()", "GeomSolids1001",
73  UWarning, 1, message.str().c_str());
74  fIsDefined = false;
76  fA = fB = fC = 0.0;
77  fDet = 0.0;
78  fArea = fRadius = 0.0;
79  }
80  else
81  {
82  fIsDefined = true;
84  fA = fE1.Mag2();
85  fB = fE1.Dot(fE2);
86  fC = fE2.Mag2();
87  fDet = fabs(fA * fC - fB * fB);
88 
89  // sMin = -0.5*kCarTolerance/sqrt(fA);
90  // sMax = 1.0 - sMin;
91  // tMin = -0.5*kCarTolerance/sqrt(fC);
92  // UVector3 vtmp = 0.25 * (fE1 + fE2);
93 
94  fArea = 0.5 * (fE1.Cross(fE2)).Mag();
95  double lambda0 , lambda1 ;
96  if(std::fabs(fArea) < kCarTolerance*kCarTolerance)
97  {
98  ostringstream message;
99  message << "Area of Facet is too small, possible flat triangle!" << endl
100  << " fVertices[0] = " << GetVertex(0) << endl
101  << " fVertices[1] = " << GetVertex(1) << endl
102  << " fVertices[2] = " << GetVertex(2) << endl
103  << "Area = " << fArea;
104  UUtils::Exception("UTriangularFacet::UTriangularFacet()",
105  "GeomSolids1001", UWarning, 1, message.str().c_str());
106  lambda0 = 0.5;
107  lambda1 = 0.5;
108  }
109  else
110  {
111  lambda0 = (fA-fB) * fC / (8.0*fArea*fArea);
112  lambda1 = (fC-fB) * fA / (8.0*fArea*fArea);
113  }
114 
115  UVector3 p0 = GetVertex(0);
116  fCircumcentre = p0 + lambda0 * fE1 + lambda1 * fE2;
117  double radiusSqr = (fCircumcentre - p0).Mag2();
118  fRadius = sqrt(radiusSqr);
119  }
120 }
121 
122 
124  : fSqrDist(0.)
125 {
126  fVertices = new vector<UVector3>(3);
127  UVector3 zero(0, 0, 0);
128  SetVertex(0, zero);
129  SetVertex(1, zero);
130  SetVertex(2, zero);
131  for (int i = 0; i < 3; ++i) fIndices[i] = -1;
132  fIsDefined = false;
133  fSurfaceNormal.Set(0);
134  fA = fB = fC = 0;
135  fE1 = zero;
136  fE2 = zero;
137  fDet = 0.0;
138  fArea = fRadius = 0;
139 }
140 
142 {
143  SetVertices(NULL);
144 }
145 
146 
148 {
149  char* p = (char*) &rhs;
150  copy(p, p + sizeof(*this), (char*)this);
151 
152  if (fIndices[0] < 0 && fVertices)
153  {
154  fVertices = new vector<UVector3>(3);
155  for (int i = 0; i < 3; ++i)(*fVertices)[i] = (*rhs.fVertices)[i];
156  }
157  fIsDefined = rhs.fIsDefined;
159  fA = rhs.fA; fB = rhs.fB; fC = rhs.fC;
160  fE1 = rhs.fE1;
161  fE2 = rhs.fE2;
162  fDet = rhs.fDet;
163  fArea = rhs.fArea;
164  fRadius = rhs.fRadius;
165  fSqrDist = rhs.fSqrDist;
166 }
167 
169 {
170  CopyFrom(rhs);
171 }
172 
174 {
175  SetVertices(NULL);
176 
177  if (this != &rhs)
178  CopyFrom(rhs);
179 
180  return *this;
181 }
182 
184 //
185 // GetClone
186 //
187 // Simple member function to generate fA duplicate of the triangular facet.
188 //
190 {
192  return fc;
193 }
194 
196 //
197 // GetFlippedFacet
198 //
199 // Member function to generate an identical facet, but with the normal vector
200 // pointing at 180 degrees.
201 //
203 {
205  return flipped;
206 }
207 
209 //
210 // Distance (UVector3)
211 //
212 // Determines the vector between p and the closest point on the facet to p.
213 // This is based on the algorithm published in "Geometric Tools for Computer
214 // Graphics," Philip J Scheider and David H Eberly, Elsevier Science (USA),
215 // 2003. at the time of writing, the algorithm is also available in fA
216 // technical note "Distance between point and triangle in 3D," by David Eberly
217 // at http://www.geometrictools.com/Documentation/DistancePoint3Triangle3.pdf
218 //
219 // The by-product is the square-distance fSqrDist, which is retained
220 // in case needed by the other "Distance" member functions.
221 //
223 {
224  UVector3 D = GetVertex(0) - p;
225  double d = fE1.Dot(D);
226  double e = fE2.Dot(D);
227  double f = D.Mag2();
228  double q = fB * e - fC * d;
229  double t = fB * d - fA * e;
230  fSqrDist = 0.;
231 
232  if (q + t <= fDet)
233  {
234  if (q < 0.0)
235  {
236  if (t < 0.0)
237  {
238  //
239  // We are in region 4.
240  //
241  if (d < 0.0)
242  {
243  t = 0.0;
244  if (-d >= fA)
245  {
246  q = 1.0;
247  fSqrDist = fA + 2.0 * d + f;
248  }
249  else
250  {
251  q = -d / fA;
252  fSqrDist = d * q + f;
253  }
254  }
255  else
256  {
257  q = 0.0;
258  if (e >= 0.0)
259  {
260  t = 0.0;
261  fSqrDist = f;
262  }
263  else if (-e >= fC)
264  {
265  t = 1.0;
266  fSqrDist = fC + 2.0 * e + f;
267  }
268  else
269  {
270  t = -e / fC;
271  fSqrDist = e * t + f;
272  }
273  }
274  }
275  else
276  {
277  //
278  // We are in region 3.
279  //
280  q = 0.0;
281  if (e >= 0.0)
282  {
283  t = 0.0;
284  fSqrDist = f;
285  }
286  else if (-e >= fC)
287  {
288  t = 1.0;
289  fSqrDist = fC + 2.0 * e + f;
290  }
291  else
292  {
293  t = -e / fC;
294  fSqrDist = e * t + f;
295  }
296  }
297  }
298  else if (t < 0.0)
299  {
300  //
301  // We are in region 5.
302  //
303  t = 0.0;
304  if (d >= 0.0)
305  {
306  q = 0.0;
307  fSqrDist = f;
308  }
309  else if (-d >= fA)
310  {
311  q = 1.0;
312  fSqrDist = fA + 2.0 * d + f;
313  }
314  else
315  {
316  q = -d / fA;
317  fSqrDist = d * q + f;
318  }
319  }
320  else
321  {
322  //
323  // We are in region 0.
324  //
325  q = q / fDet;
326  t = t / fDet;
327  fSqrDist = q * (fA * q + fB * t + 2.0 * d) + t * (fB * q + fC * t + 2.0 * e) + f;
328  }
329  }
330  else
331  {
332  if (q < 0.0)
333  {
334  //
335  // We are in region 2.
336  //
337  double tmp0 = fB + d;
338  double tmp1 = fC + e;
339  if (tmp1 > tmp0)
340  {
341  double numer = tmp1 - tmp0;
342  double denom = fA - 2.0 * fB + fC;
343  if (numer >= denom)
344  {
345  q = 1.0;
346  t = 0.0;
347  fSqrDist = fA + 2.0 * d + f;
348  }
349  else
350  {
351  q = numer / denom;
352  t = 1.0 - q;
353  fSqrDist = q * (fA * q + fB * t + 2.0 * d) + t * (fB * q + fC * t + 2.0 * e) + f;
354  }
355  }
356  else
357  {
358  q = 0.0;
359  if (tmp1 <= 0.0)
360  {
361  t = 1.0;
362  fSqrDist = fC + 2.0 * e + f;
363  }
364  else if (e >= 0.0)
365  {
366  t = 0.0;
367  fSqrDist = f;
368  }
369  else
370  {
371  t = -e / fC;
372  fSqrDist = e * t + f;
373  }
374  }
375  }
376  else if (t < 0.0)
377  {
378  //
379  // We are in region 6.
380  //
381  double tmp0 = fB + e;
382  double tmp1 = fA + d;
383  if (tmp1 > tmp0)
384  {
385  double numer = tmp1 - tmp0;
386  double denom = fA - 2.0 * fB + fC;
387  if (numer >= denom)
388  {
389  t = 1.0;
390  q = 0.0;
391  fSqrDist = fC + 2.0 * e + f;
392  }
393  else
394  {
395  t = numer / denom;
396  q = 1.0 - t;
397  fSqrDist = q * (fA * q + fB * t + 2.0 * d) + t * (fB * q + fC * t + 2.0 * e) + f;
398  }
399  }
400  else
401  {
402  t = 0.0;
403  if (tmp1 <= 0.0)
404  {
405  q = 1.0;
406  fSqrDist = fA + 2.0 * d + f;
407  }
408  else if (d >= 0.0)
409  {
410  q = 0.0;
411  fSqrDist = f;
412  }
413  else
414  {
415  q = -d / fA;
416  fSqrDist = d * q + f;
417  }
418  }
419  }
420  else
421  //
422  // We are in region 1.
423  //
424  {
425  double numer = fC + e - fB - d;
426  if (numer <= 0.0)
427  {
428  q = 0.0;
429  t = 1.0;
430  fSqrDist = fC + 2.0 * e + f;
431  }
432  else
433  {
434  double denom = fA - 2.0 * fB + fC;
435  if (numer >= denom)
436  {
437  q = 1.0;
438  t = 0.0;
439  fSqrDist = fA + 2.0 * d + f;
440  }
441  else
442  {
443  q = numer / denom;
444  t = 1.0 - q;
445  fSqrDist = q * (fA * q + fB * t + 2.0 * d) + t * (fB * q + fC * t + 2.0 * e) + f;
446  }
447  }
448  }
449  }
450  //
451  //
452  // Do fA check for rounding errors in the distance-squared. It appears that
453  // the conventional methods for calculating fSqrDist breaks down when very near
454  // to or at the surface (as required by transport). We'll therefore also use
455  // the magnitude-squared of the vector displacement. (Note that I've also
456  // tried to get around this problem by using the existing equations for
457  //
458  // fSqrDist = function(fA,fB,fC,d,q,t)
459  //
460  // and use fA more accurate addition process which minimises errors and
461  // breakdown of cummutitivity [where (A+B)+C != A+(B+C)] but this still
462  // doesn't work.
463  // Calculation from u = D + q*fE1 + t*fE2 is less efficient, but appears
464  // more robust.
465  //
466  if (fSqrDist < 0.0) fSqrDist = 0.;
467  UVector3 u = D + q * fE1 + t * fE2;
468  double u2 = u.Mag2();
469  //
470  // The following (part of the roundoff error check) is from Oliver Merle'q
471  // updates.
472  //
473  if (fSqrDist > u2) fSqrDist = u2;
474 
475  return u;
476 }
477 
479 //
480 // Distance (UVector3, double)
481 //
482 // Determines the closest distance between point p and the facet. This makes
483 // use of UVector3 UTriangularFacet::Distance, which stores the
484 // square of the distance in variable fSqrDist. If approximate methods show
485 // the distance is to be greater than minDist, then forget about further
486 // computation and return fA very large number.
487 //
488 double UTriangularFacet::Distance(const UVector3& p, const double minDist)
489 {
490  //
491  // Start with quicky test to determine if the surface of the sphere enclosing
492  // the triangle is any closer to p than minDist. If not, then don't bother
493  // about more accurate test.
494  //
495  double dist = UUtils::kInfinity;
496  if ((p - fCircumcentre).Mag() - fRadius < minDist)
497  {
498  //
499  // It's possible that the triangle is closer than minDist, so do more accurate
500  // assessment.
501  //
502  dist = Distance(p).Mag();
503  // dist = sqrt(fSqrDist);
504  }
505  return dist;
506 }
507 
509 //
510 // Distance (UVector3, double, bool)
511 //
512 // Determine the distance to point p. UUtils::kInfinity is returned if either:
513 // (1) outgoing is TRUE and the dot product of the normal vector to the facet
514 // and the displacement vector from p to the triangle is negative.
515 // (2) outgoing is FALSE and the dot product of the normal vector to the facet
516 // and the displacement vector from p to the triangle is positive.
517 // If approximate methods show the distance is to be greater than minDist, then
518 // forget about further computation and return fA very large number.
519 //
520 // This method has been heavily modified thanks to the valuable comments and
521 // corrections of Rickard Holmberg.
522 //
523 double UTriangularFacet::Distance(const UVector3& p, const double minDist, const bool outgoing)
524 {
525  //
526  // Start with quicky test to determine if the surface of the sphere enclosing
527  // the triangle is any closer to p than minDist. If not, then don't bother
528  // about more accurate test.
529  //
530  double dist = UUtils::kInfinity;
531  if ((p - fCircumcentre).Mag() - fRadius < minDist)
532  {
533  //
534  // It's possible that the triangle is closer than minDist, so do more accurate
535  // assessment.
536  //
537  UVector3 v = Distance(p);
538  double dist1 = sqrt(fSqrDist);
539  double dir = v.Dot(fSurfaceNormal);
540  bool wrongSide = (dir > 0.0 && !outgoing) || (dir < 0.0 && outgoing);
541  if (dist1 <= kCarTolerance)
542  {
543  //
544  // Point p is very close to triangle. Check if it's on the wrong side, in
545  // which case return distance of 0.0 otherwise .
546  //
547  if (wrongSide) dist = 0.0;
548  else dist = dist1;
549  }
550  else if (!wrongSide) dist = dist1;
551  }
552  return dist;
553 }
554 
556 //
557 // Extent
558 //
559 // Calculates the furthest the triangle extends in fA particular direction
560 // defined by the vector axis.
561 //
563 {
564  double ss = GetVertex(0).Dot(axis);
565  double sp = GetVertex(1).Dot(axis);
566  if (sp > ss) ss = sp;
567  sp = GetVertex(2).Dot(axis);
568  if (sp > ss) ss = sp;
569  return ss;
570 }
571 
573 //
574 // Intersect
575 //
576 // Member function to find the next intersection when going from p in the
577 // direction of v. If:
578 // (1) "outgoing" is TRUE, only consider the face if we are going out through
579 // the face.
580 // (2) "outgoing" is FALSE, only consider the face if we are going in through
581 // the face.
582 // Member functions returns TRUE if there is an intersection, FALSE otherwise.
583 // Sets the distance (distance along w), distFromSurface (orthogonal distance)
584 // and normal.
585 //
586 // Also considers intersections that happen with negative distance for small
587 // distances of distFromSurface = 0.5*kCarTolerance in the wrong direction.
588 // This is to detect kSurface without doing fA full Inside(p) in
589 // UTessellatedSolid::Distance(p,v) calculation.
590 //
591 // This member function is thanks the valuable work of Rickard Holmberg. PT.
592 // However, "gotos" are the Work of the Devil have been exorcised with
593 // extreme prejudice!!
594 //
595 // IMPORTANT NOTE: These calculations are predicated on v being fA unit
596 // vector. If UTessellatedSolid or other classes call this member function
597 // with |v| != 1 then there will be errors.
598 //
599 bool UTriangularFacet::Intersect(const UVector3& p, const UVector3& v, bool outgoing, double& distance, double& distFromSurface, UVector3& normal)
600 {
601  //
602  // Check whether the direction of the facet is consistent with the vector v
603  // and the need to be outgoing or ingoing. If inconsistent, disregard and
604  // return false.
605  //
606  double w = v.Dot(fSurfaceNormal);
607  if ((outgoing && w < -dirTolerance) || (!outgoing && w > dirTolerance))
608  {
609  distance = UUtils::kInfinity;
610  distFromSurface = UUtils::kInfinity;
611  normal.Set(0);
612  return false;
613  }
614  //
615  // Calculate the orthogonal distance from p to the surface containing the
616  // triangle. Then determine if we're on the right or wrong side of the
617  // surface (at fA distance greater than kCarTolerance to be consistent with
618  // "outgoing".
619  //
620  const UVector3& p0 = GetVertex(0);
621  UVector3 D = p0 - p;
622  distFromSurface = D.Dot(fSurfaceNormal);
623  bool wrongSide = (outgoing && distFromSurface < -0.5 * kCarTolerance) ||
624  (!outgoing && distFromSurface > 0.5 * kCarTolerance);
625  if (wrongSide)
626  {
627  distance = UUtils::kInfinity;
628  distFromSurface = UUtils::kInfinity;
629  normal.Set(0);
630  return false;
631  }
632 
633  wrongSide = (outgoing && distFromSurface < 0.0) || (!outgoing && distFromSurface > 0.0);
634  if (wrongSide)
635  {
636  //
637  // We're slightly on the wrong side of the surface. Check if we're close
638  // enough using fA precise distance calculation.
639  //
640  //UVector3 u = Distance(p);
642  {
643  //
644  // We're very close. Therefore return fA small negative number to pretend
645  // we intersect.
646  //
647  // distance = -0.5*kCarTolerance
648  distance = 0.0;
649  normal = fSurfaceNormal;
650  return true;
651  }
652  else
653  {
654  //
655  // We're close to the surface containing the triangle, but sufficiently
656  // far from the triangle, and on the wrong side compared to the directions
657  // of the surface normal and v. There is no intersection.
658  //
659  distance = UUtils::kInfinity;
660  distFromSurface = UUtils::kInfinity;
661  normal.Set(0);
662  return false;
663  }
664  }
665  if (w < dirTolerance && w > -dirTolerance)
666  {
667  //
668  // The ray is within the plane of the triangle. Project the problem into 2D
669  // in the plane of the triangle. First try to create orthogonal unit vectors
670  // mu and nu, where mu is fE1/|fE1|. This is kinda like
671  // the original algorithm due to Rickard Holmberg, but with better mathematical
672  // justification than the original method ... however, beware Rickard's was less
673  // time-consuming.
674  //
675  // Note that vprime is not fA unit vector. We need to keep it unnormalised
676  // since the values of distance along vprime (s0 and s1) for intersection with
677  // the triangle will be used to determine if we cut the plane at the same
678  // time.
679  //
680  UVector3 mu = fE1.Unit();
681  UVector3 nu = fSurfaceNormal.Cross(mu);
682  UVector2 pprime(p.Dot(mu), p.Dot(nu));
683  UVector2 vprime(v.Dot(mu), v.Dot(nu));
684  UVector2 P0prime(p0.Dot(mu), p0.Dot(nu));
685  UVector2 E0prime(fE1.Mag(), 0.0);
686  UVector2 E1prime(fE2.Dot(mu), fE2.Dot(nu));
687  UVector2 loc[2];
688  if (UTessellatedGeometryAlgorithms::IntersectLineAndTriangle2D(pprime, vprime, P0prime, E0prime, E1prime, loc))
689  {
690  //
691  // There is an intersection between the line and triangle in 2D. Now check
692  // which part of the line intersects with the plane containing the triangle
693  // in 3D.
694  //
695  double vprimemag = vprime.mag();
696  double s0 = (loc[0] - pprime).mag() / vprimemag;
697  double s1 = (loc[1] - pprime).mag() / vprimemag;
698  double normDist0 = fSurfaceNormal.Dot(s0 * v) - distFromSurface;
699  double normDist1 = fSurfaceNormal.Dot(s1 * v) - distFromSurface;
700 
701  if ((normDist0 < 0.0 && normDist1 < 0.0) || (normDist0 > 0.0 && normDist1 > 0.0) ||
702  (normDist0 == 0.0 && normDist1 == 0.0))
703  {
704  distance = UUtils::kInfinity;
705  distFromSurface = UUtils::kInfinity;
706  normal.Set(0);
707  return false;
708  }
709  else
710  {
711  double dnormDist = normDist1 - normDist0;
712  if (fabs(dnormDist) < DBL_EPSILON)
713  {
714  distance = s0;
715  normal = fSurfaceNormal;
716  if (!outgoing) distFromSurface = -distFromSurface;
717  return true;
718  }
719  else
720  {
721  distance = s0 - normDist0 * (s1 - s0) / dnormDist;
722  normal = fSurfaceNormal;
723  if (!outgoing) distFromSurface = -distFromSurface;
724  return true;
725  }
726  }
727 
728  // UVector3 dloc = loc1 - loc0;
729  // UVector3 dlocXv = dloc.cross(v);
730  // double dlocXvmag = dlocXv.mag();
731  // if (dloc.mag() <= 0.5*kCarTolerance || dlocXvmag <= DBL_EPSILON)
732  // {
733  // distance = loc0.mag();
734  // normal = fSurfaceNormal;
735  // if (!outgoing) distFromSurface = -distFromSurface;
736  // return true;
737  // }
738 
739  // UVector3 loc0Xv = loc0.cross(v);
740  // UVector3 loc1Xv = loc1.cross(v);
741  // double sameDir = -loc0Xv.dot(loc1Xv);
742  // if (sameDir < 0.0)
743  // {
744  // distance = UUtils::kInfinity;
745  // distFromSurface = UUtils::kInfinity;
746  // normal = UVector3(0.0,0.0,0.0);
747  // return false;
748  // }
749  // else
750  // {
751  // distance = loc0.mag() + loc0Xv.mag() * dloc.mag()/dlocXvmag;
752  // normal = fSurfaceNormal;
753  // if (!outgoing) distFromSurface = -distFromSurface;
754  // return true;
755  // }
756  }
757  else
758  {
759  distance = UUtils::kInfinity;
760  distFromSurface = UUtils::kInfinity;
761  normal.Set(0);
762  return false;
763  }
764  }
765  //
766  //
767  // Use conventional algorithm to determine the whether there is an
768  // intersection. This involves determining the point of intersection of the
769  // line with the plane containing the triangle, and then calculating if the
770  // point is within the triangle.
771  //
772  distance = distFromSurface / w;
773  UVector3 pp = p + v * distance;
774  UVector3 DD = p0 - pp;
775  double d = fE1.Dot(DD);
776  double e = fE2.Dot(DD);
777  double ss = fB * e - fC * d;
778  double t = fB * d - fA * e;
779 
780  double sTolerance = (fabs(fB) + fabs(fC) + fabs(d) + fabs(e)) * kCarTolerance;
781  double tTolerance = (fabs(fA) + fabs(fB) + fabs(d) + fabs(e)) * kCarTolerance;
782  double detTolerance = (fabs(fA) + fabs(fC) + 2 * fabs(fB)) * kCarTolerance;
783 
784  //if (ss < 0.0 || t < 0.0 || ss+t > fDet)
785  if (ss < -sTolerance || t < -tTolerance || (ss + t - fDet) > detTolerance)
786  {
787  //
788  // The intersection is outside of the triangle.
789  //
790  distance = distFromSurface = UUtils::kInfinity;
791  normal.Set(0);
792  return false;
793  }
794  else
795  {
796  //
797  // There is an intersection. Now we only need to set the surface normal.
798  //
799  normal = fSurfaceNormal;
800  if (!outgoing) distFromSurface = -distFromSurface;
801  return true;
802  }
803 }
804 
806 //
807 // GetPointOnFace
808 //```
809 // Auxiliary method for get fA random point on surface
810 
812 {
813  double alpha = UUtils::Random(0., 1.);
814  double beta = UUtils::Random(0., 1.);
815  double lambda1 = alpha * beta;
816  double lambda0 = alpha - lambda1;
817 
818  return GetVertex(0) + lambda0 * fE1 + lambda1 * fE2;
819 }
820 
822 //
823 // GetArea
824 //
825 // Auxiliary method for returning the surface fArea
826 
828 {
829  return fArea;
830 }
831 
833 {
834  return "TriangularFacet";
835 }
836 
838 {
839  return fSurfaceNormal;
840 }
841 
843 {
845 }
double Mag2() const
Definition: UVector3.hh:287
UTriangularFacet & operator=(const UTriangularFacet &right)
static const double dirTolerance
Definition: VUFacet.hh:88
UVector3 GetPointOnFace() const
UVector3 GetSurfaceNormal() const
void SetSurfaceNormal(UVector3 normal)
UVector3 Cross(const UVector3 &) const
Definition: UVector3.hh:281
static bool IntersectLineAndTriangle2D(const UVector2 &p, const UVector2 &v, const UVector2 &p0, const UVector2 &e0, const UVector2 &e1, UVector2 location[2])
UFacetVertexType
Definition: VUFacet.hh:49
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
void copy(std::vector< T > &main, const std::vector< T > &data)
Definition: DicomRun.hh:91
static const double kInfinity
Definition: UUtils.hh:54
void CopyFrom(const UTriangularFacet &rhs)
void SetVertex(int i, const UVector3 &val)
#define DBL_EPSILON
Definition: templates.hh:87
double Dot(const UVector3 &) const
Definition: UVector3.hh:276
double Mag() const
Definition: UVector3.cc:48
double mag() const
Definition: UVector2.hh:310
UTriangularFacet * GetFlippedFacet()
void Exception(const char *originOfException, const char *exceptionCode, UExceptionSeverity severity, int level, const char *description)
Definition: UUtils.cc:122
UVector3 GetVertex(int i) const
const G4double p0
void Set(double xx, double yy, double zz)
Definition: UVector3.hh:264
std::vector< UVector3 > * fVertices
UVector3 Unit() const
Definition: UVector3.cc:80
std::string UGeometryType
Definition: UTypes.hh:39
double Extent(const UVector3 axis)
bool Intersect(const UVector3 &p, const UVector3 &v, const bool outgoing, double &distance, double &distFromSurface, UVector3 &normal)
UGeometryType GetEntityType() const
void SetVertices(std::vector< UVector3 > *v)
static const G4double alpha
static const double kCarTolerance
Definition: VUFacet.hh:89
double Random(double min=0.0, double max=1.0)
Definition: UUtils.cc:69
UVector3 Distance(const UVector3 &p)