Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4TriangularFacet Class Reference

#include <G4TriangularFacet.hh>

+ Inheritance diagram for G4TriangularFacet:

Public Member Functions

 G4TriangularFacet ()
 
 ~G4TriangularFacet ()
 
 G4TriangularFacet (const G4ThreeVector &vt0, const G4ThreeVector &vt1, const G4ThreeVector &vt2, G4FacetVertexType)
 
 G4TriangularFacet (const G4TriangularFacet &right)
 
 G4TriangularFacet (G4TriangularFacet &&right)
 
G4TriangularFacetoperator= (const G4TriangularFacet &right)
 
G4TriangularFacetoperator= (G4TriangularFacet &&right)
 
G4VFacetGetClone ()
 
G4TriangularFacetGetFlippedFacet ()
 
G4ThreeVector Distance (const G4ThreeVector &p)
 
G4double Distance (const G4ThreeVector &p, G4double minDist)
 
G4double Distance (const G4ThreeVector &p, G4double minDist, const G4bool outgoing)
 
G4double Extent (const G4ThreeVector axis)
 
G4bool Intersect (const G4ThreeVector &p, const G4ThreeVector &v, const G4bool outgoing, G4double &distance, G4double &distFromSurface, G4ThreeVector &normal)
 
G4double GetArea () const
 
G4ThreeVector GetPointOnFace () const
 
G4ThreeVector GetSurfaceNormal () const
 
void SetSurfaceNormal (G4ThreeVector normal)
 
G4GeometryType GetEntityType () const
 
G4bool IsDefined () const
 
G4int GetNumberOfVertices () const
 
G4ThreeVector GetVertex (G4int i) const
 
void SetVertex (G4int i, const G4ThreeVector &val)
 
G4ThreeVector GetCircumcentre () const
 
G4double GetRadius () const
 
G4int AllocatedMemory ()
 
G4int GetVertexIndex (G4int i) const
 
void SetVertexIndex (G4int i, G4int j)
 
void SetVertices (std::vector< G4ThreeVector > *v)
 
- Public Member Functions inherited from G4VFacet
 G4VFacet ()
 
virtual ~G4VFacet ()
 
G4bool operator== (const G4VFacet &right) const
 
virtual G4int GetNumberOfVertices () const =0
 
virtual G4ThreeVector GetVertex (G4int i) const =0
 
virtual void SetVertex (G4int i, const G4ThreeVector &val)=0
 
virtual G4GeometryType GetEntityType () const =0
 
virtual G4ThreeVector GetSurfaceNormal () const =0
 
virtual G4bool IsDefined () const =0
 
virtual G4ThreeVector GetCircumcentre () const =0
 
virtual G4double GetRadius () const =0
 
virtual G4VFacetGetClone ()=0
 
virtual G4double Distance (const G4ThreeVector &, G4double)=0
 
virtual G4double Distance (const G4ThreeVector &, G4double, const G4bool)=0
 
virtual G4double Extent (const G4ThreeVector)=0
 
virtual G4bool Intersect (const G4ThreeVector &, const G4ThreeVector &, const G4bool, G4double &, G4double &, G4ThreeVector &)=0
 
virtual G4double GetArea () const =0
 
virtual G4ThreeVector GetPointOnFace () const =0
 
void ApplyTranslation (const G4ThreeVector v)
 
std::ostream & StreamInfo (std::ostream &os) const
 
G4bool IsInside (const G4ThreeVector &p) const
 
virtual G4int AllocatedMemory ()=0
 
virtual void SetVertexIndex (G4int i, G4int j)=0
 
virtual G4int GetVertexIndex (G4int i) const =0
 
virtual void SetVertices (std::vector< G4ThreeVector > *vertices)=0
 

Additional Inherited Members

- Protected Attributes inherited from G4VFacet
G4double kCarTolerance
 
- Static Protected Attributes inherited from G4VFacet
static const G4double dirTolerance = 1.0E-14
 

Detailed Description

Definition at line 60 of file G4TriangularFacet.hh.

Constructor & Destructor Documentation

◆ G4TriangularFacet() [1/4]

G4TriangularFacet::G4TriangularFacet ( )

Definition at line 144 of file G4TriangularFacet.cc.

145 : fSqrDist(0.)
146{
147 fVertices = new vector<G4ThreeVector>(3);
148 G4ThreeVector zero(0,0,0);
149 SetVertex(0, zero);
150 SetVertex(1, zero);
151 SetVertex(2, zero);
152 for (G4int i = 0; i < 3; ++i) fIndices[i] = -1;
153 fIsDefined = false;
154 fSurfaceNormal.set(0,0,0);
155 fA = fB = fC = 0;
156 fE1 = zero;
157 fE2 = zero;
158 fDet = 0.0;
159 fArea = fRadius = 0.0;
160}
int G4int
Definition: G4Types.hh:85
void set(double x, double y, double z)
void SetVertex(G4int i, const G4ThreeVector &val)

Referenced by GetClone(), and GetFlippedFacet().

◆ ~G4TriangularFacet()

G4TriangularFacet::~G4TriangularFacet ( )

Definition at line 164 of file G4TriangularFacet.cc.

165{
166 SetVertices(nullptr);
167}
void SetVertices(std::vector< G4ThreeVector > *v)

◆ G4TriangularFacet() [2/4]

G4TriangularFacet::G4TriangularFacet ( const G4ThreeVector vt0,
const G4ThreeVector vt1,
const G4ThreeVector vt2,
G4FacetVertexType  vertexType 
)

Definition at line 56 of file G4TriangularFacet.cc.

60 : G4VFacet()
61{
62 fVertices = new vector<G4ThreeVector>(3);
63
64 SetVertex(0, vt0);
65 if (vertexType == ABSOLUTE)
66 {
67 SetVertex(1, vt1);
68 SetVertex(2, vt2);
69 fE1 = vt1 - vt0;
70 fE2 = vt2 - vt0;
71 }
72 else
73 {
74 SetVertex(1, vt0 + vt1);
75 SetVertex(2, vt0 + vt2);
76 fE1 = vt1;
77 fE2 = vt2;
78 }
79
80 G4ThreeVector E1xE2 = fE1.cross(fE2);
81 fArea = 0.5 * E1xE2.mag();
82 for (G4int i = 0; i < 3; ++i) fIndices[i] = -1;
83
84 fIsDefined = true;
85 G4double delta = kCarTolerance; // Set tolerance for checking
86
87 // Check length of edges
88 //
89 G4double leng1 = fE1.mag();
90 G4double leng2 = (fE2-fE1).mag();
91 G4double leng3 = fE2.mag();
92 if (leng1 <= delta || leng2 <= delta || leng3 <= delta)
93 {
94 fIsDefined = false;
95 }
96
97 // Check min height of triangle
98 //
99 if (fIsDefined)
100 {
101 if (2.*fArea/std::max(std::max(leng1,leng2),leng3) <= delta)
102 {
103 fIsDefined = false;
104 }
105 }
106
107 // Define facet
108 //
109 if (!fIsDefined)
110 {
111 ostringstream message;
112 message << "Facet is too small or too narrow." << G4endl
113 << "Triangle area = " << fArea << G4endl
114 << "P0 = " << GetVertex(0) << G4endl
115 << "P1 = " << GetVertex(1) << G4endl
116 << "P2 = " << GetVertex(2) << G4endl
117 << "Side1 length (P0->P1) = " << leng1 << G4endl
118 << "Side2 length (P1->P2) = " << leng2 << G4endl
119 << "Side3 length (P2->P0) = " << leng3;
120 G4Exception("G4TriangularFacet::G4TriangularFacet()",
121 "GeomSolids1001", JustWarning, message);
122 fSurfaceNormal.set(0,0,0);
123 fA = fB = fC = 0.0;
124 fDet = 0.0;
125 fCircumcentre = vt0 + 0.5*fE1 + 0.5*fE2;
126 fArea = fRadius = 0.0;
127 }
128 else
129 {
130 fSurfaceNormal = E1xE2.unit();
131 fA = fE1.mag2();
132 fB = fE1.dot(fE2);
133 fC = fE2.mag2();
134 fDet = std::fabs(fA*fC - fB*fB);
135
136 fCircumcentre =
137 vt0 + (E1xE2.cross(fE1)*fC + fE2.cross(E1xE2)*fA) / (2.*E1xE2.mag2());
138 fRadius = (fCircumcentre - vt0).mag();
139 }
140}
@ JustWarning
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:59
double G4double
Definition: G4Types.hh:83
@ ABSOLUTE
Definition: G4VFacet.hh:48
#define G4endl
Definition: G4ios.hh:57
Hep3Vector unit() const
double mag2() const
Hep3Vector cross(const Hep3Vector &) const
double dot(const Hep3Vector &) const
double mag() const
G4ThreeVector GetVertex(G4int i) const
G4VFacet()
Definition: G4VFacet.cc:44
G4double kCarTolerance
Definition: G4VFacet.hh:93

◆ G4TriangularFacet() [3/4]

G4TriangularFacet::G4TriangularFacet ( const G4TriangularFacet right)

Definition at line 203 of file G4TriangularFacet.cc.

204 : G4VFacet(rhs)
205{
206 CopyFrom(rhs);
207}

◆ G4TriangularFacet() [4/4]

G4TriangularFacet::G4TriangularFacet ( G4TriangularFacet &&  right)

Definition at line 211 of file G4TriangularFacet.cc.

212 : G4VFacet(rhs)
213{
214 MoveFrom(rhs);
215}

Member Function Documentation

◆ AllocatedMemory()

G4int G4TriangularFacet::AllocatedMemory ( )
inlinevirtual

Implements G4VFacet.

Definition at line 162 of file G4TriangularFacet.hh.

163{
164 G4int size = sizeof(*this);
165 size += GetNumberOfVertices() * sizeof(G4ThreeVector);
166 return size;
167}
CLHEP::Hep3Vector G4ThreeVector
G4int GetNumberOfVertices() const

◆ Distance() [1/3]

G4ThreeVector G4TriangularFacet::Distance ( const G4ThreeVector p)

Definition at line 290 of file G4TriangularFacet.cc.

291{
292 G4ThreeVector D = GetVertex(0) - p;
293 G4double d = fE1.dot(D);
294 G4double e = fE2.dot(D);
295 G4double f = D.mag2();
296 G4double q = fB*e - fC*d;
297 G4double t = fB*d - fA*e;
298 fSqrDist = 0.;
299
300 if (q+t <= fDet)
301 {
302 if (q < 0.0)
303 {
304 if (t < 0.0)
305 {
306 //
307 // We are in region 4.
308 //
309 if (d < 0.0)
310 {
311 t = 0.0;
312 if (-d >= fA) {q = 1.0; fSqrDist = fA + 2.0*d + f;}
313 else {q = -d/fA; fSqrDist = d*q + f;}
314 }
315 else
316 {
317 q = 0.0;
318 if (e >= 0.0) {t = 0.0; fSqrDist = f;}
319 else if (-e >= fC) {t = 1.0; fSqrDist = fC + 2.0*e + f;}
320 else {t = -e/fC; fSqrDist = e*t + f;}
321 }
322 }
323 else
324 {
325 //
326 // We are in region 3.
327 //
328 q = 0.0;
329 if (e >= 0.0) {t = 0.0; fSqrDist = f;}
330 else if (-e >= fC) {t = 1.0; fSqrDist = fC + 2.0*e + f;}
331 else {t = -e/fC; fSqrDist = e*t + f;}
332 }
333 }
334 else if (t < 0.0)
335 {
336 //
337 // We are in region 5.
338 //
339 t = 0.0;
340 if (d >= 0.0) {q = 0.0; fSqrDist = f;}
341 else if (-d >= fA) {q = 1.0; fSqrDist = fA + 2.0*d + f;}
342 else {q = -d/fA; fSqrDist = d*q + f;}
343 }
344 else
345 {
346 //
347 // We are in region 0.
348 //
349 G4double dist = fSurfaceNormal.dot(D);
350 fSqrDist = dist*dist;
351 return fSurfaceNormal*dist;
352 }
353 }
354 else
355 {
356 if (q < 0.0)
357 {
358 //
359 // We are in region 2.
360 //
361 G4double tmp0 = fB + d;
362 G4double tmp1 = fC + e;
363 if (tmp1 > tmp0)
364 {
365 G4double numer = tmp1 - tmp0;
366 G4double denom = fA - 2.0*fB + fC;
367 if (numer >= denom) {q = 1.0; t = 0.0; fSqrDist = fA + 2.0*d + f;}
368 else
369 {
370 q = numer/denom;
371 t = 1.0 - q;
372 fSqrDist = q*(fA*q + fB*t +2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
373 }
374 }
375 else
376 {
377 q = 0.0;
378 if (tmp1 <= 0.0) {t = 1.0; fSqrDist = fC + 2.0*e + f;}
379 else if (e >= 0.0) {t = 0.0; fSqrDist = f;}
380 else {t = -e/fC; fSqrDist = e*t + f;}
381 }
382 }
383 else if (t < 0.0)
384 {
385 //
386 // We are in region 6.
387 //
388 G4double tmp0 = fB + e;
389 G4double tmp1 = fA + d;
390 if (tmp1 > tmp0)
391 {
392 G4double numer = tmp1 - tmp0;
393 G4double denom = fA - 2.0*fB + fC;
394 if (numer >= denom) {t = 1.0; q = 0.0; fSqrDist = fC + 2.0*e + f;}
395 else
396 {
397 t = numer/denom;
398 q = 1.0 - t;
399 fSqrDist = q*(fA*q + fB*t +2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
400 }
401 }
402 else
403 {
404 t = 0.0;
405 if (tmp1 <= 0.0) {q = 1.0; fSqrDist = fA + 2.0*d + f;}
406 else if (d >= 0.0) {q = 0.0; fSqrDist = f;}
407 else {q = -d/fA; fSqrDist = d*q + f;}
408 }
409 }
410 else
411 //
412 // We are in region 1.
413 //
414 {
415 G4double numer = fC + e - fB - d;
416 if (numer <= 0.0)
417 {
418 q = 0.0;
419 t = 1.0;
420 fSqrDist = fC + 2.0*e + f;
421 }
422 else
423 {
424 G4double denom = fA - 2.0*fB + fC;
425 if (numer >= denom) {q = 1.0; t = 0.0; fSqrDist = fA + 2.0*d + f;}
426 else
427 {
428 q = numer/denom;
429 t = 1.0 - q;
430 fSqrDist = q*(fA*q + fB*t + 2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
431 }
432 }
433 }
434 }
435 //
436 //
437 // Do fA check for rounding errors in the distance-squared. It appears that
438 // the conventional methods for calculating fSqrDist breaks down when very
439 // near to or at the surface (as required by transport).
440 // We'll therefore also use the magnitude-squared of the vector displacement.
441 // (Note that I've also tried to get around this problem by using the
442 // existing equations for
443 //
444 // fSqrDist = function(fA,fB,fC,d,q,t)
445 //
446 // and use fA more accurate addition process which minimises errors and
447 // breakdown of cummutitivity [where (A+B)+C != A+(B+C)] but this still
448 // doesn't work.
449 // Calculation from u = D + q*fE1 + t*fE2 is less efficient, but appears
450 // more robust.
451 //
452 if (fSqrDist < 0.0) fSqrDist = 0.;
453 G4ThreeVector u = D + q*fE1 + t*fE2;
454 G4double u2 = u.mag2();
455 //
456 // The following (part of the roundoff error check) is from Oliver Merle'q
457 // updates.
458 //
459 if (fSqrDist > u2) fSqrDist = u2;
460
461 return u;
462}
G4double D(G4double temp)

Referenced by G4QuadrangularFacet::Distance(), Distance(), and Intersect().

◆ Distance() [2/3]

G4double G4TriangularFacet::Distance ( const G4ThreeVector p,
G4double  minDist 
)
virtual

Implements G4VFacet.

Definition at line 474 of file G4TriangularFacet.cc.

476{
477 //
478 // Start with quicky test to determine if the surface of the sphere enclosing
479 // the triangle is any closer to p than minDist. If not, then don't bother
480 // about more accurate test.
481 //
482 G4double dist = kInfinity;
483 if ((p-fCircumcentre).mag()-fRadius < minDist)
484 {
485 //
486 // It's possible that the triangle is closer than minDist,
487 // so do more accurate assessment.
488 //
489 dist = Distance(p).mag();
490 }
491 return dist;
492}
G4ThreeVector Distance(const G4ThreeVector &p)

◆ Distance() [3/3]

G4double G4TriangularFacet::Distance ( const G4ThreeVector p,
G4double  minDist,
const G4bool  outgoing 
)
virtual

Implements G4VFacet.

Definition at line 509 of file G4TriangularFacet.cc.

512{
513 //
514 // Start with quicky test to determine if the surface of the sphere enclosing
515 // the triangle is any closer to p than minDist. If not, then don't bother
516 // about more accurate test.
517 //
518 G4double dist = kInfinity;
519 if ((p-fCircumcentre).mag()-fRadius < minDist)
520 {
521 //
522 // It's possible that the triangle is closer than minDist,
523 // so do more accurate assessment.
524 //
525 G4ThreeVector v = Distance(p);
526 G4double dist1 = sqrt(fSqrDist);
527 G4double dir = v.dot(fSurfaceNormal);
528 G4bool wrongSide = (dir > 0.0 && !outgoing) || (dir < 0.0 && outgoing);
529 if (dist1 <= kCarTolerance)
530 {
531 //
532 // Point p is very close to triangle. Check if it's on the wrong side,
533 // in which case return distance of 0.0 otherwise .
534 //
535 if (wrongSide) dist = 0.0;
536 else dist = dist1;
537 }
538 else if (!wrongSide) dist = dist1;
539 }
540 return dist;
541}
bool G4bool
Definition: G4Types.hh:86

◆ Extent()

G4double G4TriangularFacet::Extent ( const G4ThreeVector  axis)
virtual

Implements G4VFacet.

Definition at line 550 of file G4TriangularFacet.cc.

551{
552 G4double ss = GetVertex(0).dot(axis);
553 G4double sp = GetVertex(1).dot(axis);
554 if (sp > ss) ss = sp;
555 sp = GetVertex(2).dot(axis);
556 if (sp > ss) ss = sp;
557 return ss;
558}

◆ GetArea()

G4double G4TriangularFacet::GetArea ( ) const
virtual

Implements G4VFacet.

Definition at line 793 of file G4TriangularFacet.cc.

794{
795 return fArea;
796}

Referenced by G4QuadrangularFacet::GetArea(), and G4QuadrangularFacet::GetPointOnFace().

◆ GetCircumcentre()

G4ThreeVector G4TriangularFacet::GetCircumcentre ( ) const
inlinevirtual

Implements G4VFacet.

Definition at line 152 of file G4TriangularFacet.hh.

153{
154 return fCircumcentre;
155}

◆ GetClone()

G4VFacet * G4TriangularFacet::GetClone ( )
virtual

Implements G4VFacet.

Definition at line 255 of file G4TriangularFacet.cc.

◆ GetEntityType()

G4GeometryType G4TriangularFacet::GetEntityType ( ) const
virtual

Implements G4VFacet.

Definition at line 800 of file G4TriangularFacet.cc.

801{
802 return "G4TriangularFacet";
803}

◆ GetFlippedFacet()

G4TriangularFacet * G4TriangularFacet::GetFlippedFacet ( )

Definition at line 269 of file G4TriangularFacet.cc.

270{
271 G4TriangularFacet* flipped =
273 return flipped;
274}

◆ GetNumberOfVertices()

G4int G4TriangularFacet::GetNumberOfVertices ( ) const
inlinevirtual

Implements G4VFacet.

Definition at line 136 of file G4TriangularFacet.hh.

137{
138 return 3;
139}

Referenced by AllocatedMemory().

◆ GetPointOnFace()

G4ThreeVector G4TriangularFacet::GetPointOnFace ( ) const
virtual

Implements G4VFacet.

Definition at line 779 of file G4TriangularFacet.cc.

780{
783 if (u+v > 1.) { u = 1. - u; v = 1. - v; }
784 return GetVertex(0) + u*fE1 + v*fE2;
785}
#define G4UniformRand()
Definition: Randomize.hh:52

Referenced by G4QuadrangularFacet::GetPointOnFace().

◆ GetRadius()

G4double G4TriangularFacet::GetRadius ( ) const
inlinevirtual

Implements G4VFacet.

Definition at line 157 of file G4TriangularFacet.hh.

158{
159 return fRadius;
160}

◆ GetSurfaceNormal()

G4ThreeVector G4TriangularFacet::GetSurfaceNormal ( ) const
virtual

Implements G4VFacet.

Definition at line 807 of file G4TriangularFacet.cc.

808{
809 return fSurfaceNormal;
810}

Referenced by G4QuadrangularFacet::GetSurfaceNormal().

◆ GetVertex()

G4ThreeVector G4TriangularFacet::GetVertex ( G4int  i) const
inlinevirtual

Implements G4VFacet.

Definition at line 141 of file G4TriangularFacet.hh.

142{
143 G4int indice = fIndices[i];
144 return indice < 0 ? (*fVertices)[i] : (*fVertices)[indice];
145}

Referenced by Distance(), Extent(), G4TriangularFacet(), GetClone(), GetFlippedFacet(), GetPointOnFace(), G4QuadrangularFacet::GetVertex(), and Intersect().

◆ GetVertexIndex()

G4int G4TriangularFacet::GetVertexIndex ( G4int  i) const
inlinevirtual

Implements G4VFacet.

Definition at line 169 of file G4TriangularFacet.hh.

170{
171 if (i < 3) return fIndices[i];
172 else return 999999999;
173}

◆ Intersect()

G4bool G4TriangularFacet::Intersect ( const G4ThreeVector p,
const G4ThreeVector v,
const G4bool  outgoing,
G4double distance,
G4double distFromSurface,
G4ThreeVector normal 
)
virtual

Implements G4VFacet.

Definition at line 587 of file G4TriangularFacet.cc.

593{
594 //
595 // Check whether the direction of the facet is consistent with the vector v
596 // and the need to be outgoing or ingoing. If inconsistent, disregard and
597 // return false.
598 //
599 G4double w = v.dot(fSurfaceNormal);
600 if ((outgoing && w < -dirTolerance) || (!outgoing && w > dirTolerance))
601 {
602 distance = kInfinity;
603 distFromSurface = kInfinity;
604 normal.set(0,0,0);
605 return false;
606 }
607 //
608 // Calculate the orthogonal distance from p to the surface containing the
609 // triangle. Then determine if we're on the right or wrong side of the
610 // surface (at fA distance greater than kCarTolerance to be consistent with
611 // "outgoing".
612 //
613 const G4ThreeVector& p0 = GetVertex(0);
614 G4ThreeVector D = p0 - p;
615 distFromSurface = D.dot(fSurfaceNormal);
616 G4bool wrongSide = (outgoing && distFromSurface < -0.5*kCarTolerance) ||
617 (!outgoing && distFromSurface > 0.5*kCarTolerance);
618
619 if (wrongSide)
620 {
621 distance = kInfinity;
622 distFromSurface = kInfinity;
623 normal.set(0,0,0);
624 return false;
625 }
626
627 wrongSide = (outgoing && distFromSurface < 0.0)
628 || (!outgoing && distFromSurface > 0.0);
629 if (wrongSide)
630 {
631 //
632 // We're slightly on the wrong side of the surface. Check if we're close
633 // enough using fA precise distance calculation.
634 //
635 G4ThreeVector u = Distance(p);
636 if (fSqrDist <= kCarTolerance*kCarTolerance)
637 {
638 //
639 // We're very close. Therefore return fA small negative number
640 // to pretend we intersect.
641 //
642 // distance = -0.5*kCarTolerance
643 distance = 0.0;
644 normal = fSurfaceNormal;
645 return true;
646 }
647 else
648 {
649 //
650 // We're close to the surface containing the triangle, but sufficiently
651 // far from the triangle, and on the wrong side compared to the directions
652 // of the surface normal and v. There is no intersection.
653 //
654 distance = kInfinity;
655 distFromSurface = kInfinity;
656 normal.set(0,0,0);
657 return false;
658 }
659 }
660 if (w < dirTolerance && w > -dirTolerance)
661 {
662 //
663 // The ray is within the plane of the triangle. Project the problem into 2D
664 // in the plane of the triangle. First try to create orthogonal unit vectors
665 // mu and nu, where mu is fE1/|fE1|. This is kinda like
666 // the original algorithm due to Rickard Holmberg, but with better
667 // mathematical justification than the original method ... however,
668 // beware Rickard's was less time-consuming.
669 //
670 // Note that vprime is not fA unit vector. We need to keep it unnormalised
671 // since the values of distance along vprime (s0 and s1) for intersection
672 // with the triangle will be used to determine if we cut the plane at the
673 // same time.
674 //
675 G4ThreeVector mu = fE1.unit();
676 G4ThreeVector nu = fSurfaceNormal.cross(mu);
677 G4TwoVector pprime(p.dot(mu), p.dot(nu));
678 G4TwoVector vprime(v.dot(mu), v.dot(nu));
679 G4TwoVector P0prime(p0.dot(mu), p0.dot(nu));
680 G4TwoVector E0prime(fE1.mag(), 0.0);
681 G4TwoVector E1prime(fE2.dot(mu), fE2.dot(nu));
682 G4TwoVector loc[2];
684 vprime, P0prime, E0prime, E1prime, loc))
685 {
686 //
687 // There is an intersection between the line and triangle in 2D.
688 // Now check which part of the line intersects with the plane
689 // containing the triangle in 3D.
690 //
691 G4double vprimemag = vprime.mag();
692 G4double s0 = (loc[0] - pprime).mag()/vprimemag;
693 G4double s1 = (loc[1] - pprime).mag()/vprimemag;
694 G4double normDist0 = fSurfaceNormal.dot(s0*v) - distFromSurface;
695 G4double normDist1 = fSurfaceNormal.dot(s1*v) - distFromSurface;
696
697 if ((normDist0 < 0.0 && normDist1 < 0.0)
698 || (normDist0 > 0.0 && normDist1 > 0.0)
699 || (normDist0 == 0.0 && normDist1 == 0.0) )
700 {
701 distance = kInfinity;
702 distFromSurface = kInfinity;
703 normal.set(0,0,0);
704 return false;
705 }
706 else
707 {
708 G4double dnormDist = normDist1 - normDist0;
709 if (fabs(dnormDist) < DBL_EPSILON)
710 {
711 distance = s0;
712 normal = fSurfaceNormal;
713 if (!outgoing) distFromSurface = -distFromSurface;
714 return true;
715 }
716 else
717 {
718 distance = s0 - normDist0*(s1-s0)/dnormDist;
719 normal = fSurfaceNormal;
720 if (!outgoing) distFromSurface = -distFromSurface;
721 return true;
722 }
723 }
724 }
725 else
726 {
727 distance = kInfinity;
728 distFromSurface = kInfinity;
729 normal.set(0,0,0);
730 return false;
731 }
732 }
733 //
734 //
735 // Use conventional algorithm to determine the whether there is an
736 // intersection. This involves determining the point of intersection of the
737 // line with the plane containing the triangle, and then calculating if the
738 // point is within the triangle.
739 //
740 distance = distFromSurface / w;
741 G4ThreeVector pp = p + v*distance;
742 G4ThreeVector DD = p0 - pp;
743 G4double d = fE1.dot(DD);
744 G4double e = fE2.dot(DD);
745 G4double ss = fB*e - fC*d;
746 G4double t = fB*d - fA*e;
747
748 G4double sTolerance = (fabs(fB)+ fabs(fC) + fabs(d) + fabs(e))*kCarTolerance;
749 G4double tTolerance = (fabs(fA)+ fabs(fB) + fabs(d) + fabs(e))*kCarTolerance;
750 G4double detTolerance = (fabs(fA)+ fabs(fC) + 2*fabs(fB) )*kCarTolerance;
751
752 //if (ss < 0.0 || t < 0.0 || ss+t > fDet)
753 if (ss < -sTolerance || t < -tTolerance || ( ss+t - fDet ) > detTolerance)
754 {
755 //
756 // The intersection is outside of the triangle.
757 //
758 distance = distFromSurface = kInfinity;
759 normal.set(0,0,0);
760 return false;
761 }
762 else
763 {
764 //
765 // There is an intersection. Now we only need to set the surface normal.
766 //
767 normal = fSurfaceNormal;
768 if (!outgoing) distFromSurface = -distFromSurface;
769 return true;
770 }
771}
double mag() const
static G4bool IntersectLineAndTriangle2D(const G4TwoVector &p, const G4TwoVector &v, const G4TwoVector &p0, const G4TwoVector &e0, const G4TwoVector &e1, G4TwoVector location[2])
static const G4double dirTolerance
Definition: G4VFacet.hh:92
#define DBL_EPSILON
Definition: templates.hh:66

Referenced by G4QuadrangularFacet::Intersect().

◆ IsDefined()

G4bool G4TriangularFacet::IsDefined ( ) const
inlinevirtual

Implements G4VFacet.

Definition at line 131 of file G4TriangularFacet.hh.

132{
133 return fIsDefined;
134}

Referenced by G4QuadrangularFacet::IsDefined().

◆ operator=() [1/2]

G4TriangularFacet & G4TriangularFacet::operator= ( const G4TriangularFacet right)

Definition at line 220 of file G4TriangularFacet.cc.

221{
222 SetVertices(nullptr);
223
224 if (this != &rhs)
225 {
226 delete fVertices;
227 CopyFrom(rhs);
228 }
229
230 return *this;
231}

◆ operator=() [2/2]

G4TriangularFacet & G4TriangularFacet::operator= ( G4TriangularFacet &&  right)

Definition at line 236 of file G4TriangularFacet.cc.

237{
238 SetVertices(nullptr);
239
240 if (this != &rhs)
241 {
242 delete fVertices;
243 MoveFrom(rhs);
244 }
245
246 return *this;
247}

◆ SetSurfaceNormal()

void G4TriangularFacet::SetSurfaceNormal ( G4ThreeVector  normal)

Definition at line 814 of file G4TriangularFacet.cc.

815{
816 fSurfaceNormal = normal;
817}

Referenced by G4QuadrangularFacet::G4QuadrangularFacet().

◆ SetVertex()

void G4TriangularFacet::SetVertex ( G4int  i,
const G4ThreeVector val 
)
inlinevirtual

Implements G4VFacet.

Definition at line 147 of file G4TriangularFacet.hh.

148{
149 (*fVertices)[i] = val;
150}

Referenced by G4TriangularFacet(), and G4QuadrangularFacet::SetVertex().

◆ SetVertexIndex()

void G4TriangularFacet::SetVertexIndex ( G4int  i,
G4int  j 
)
inlinevirtual

Implements G4VFacet.

Definition at line 175 of file G4TriangularFacet.hh.

176{
177 fIndices[i] = j;
178}

◆ SetVertices()

void G4TriangularFacet::SetVertices ( std::vector< G4ThreeVector > *  v)
inlinevirtual

Implements G4VFacet.

Definition at line 180 of file G4TriangularFacet.hh.

181{
182 if (fIndices[0] < 0 && fVertices)
183 {
184 delete fVertices;
185 fVertices = nullptr;
186 }
187 fVertices = v;
188}

Referenced by operator=(), G4QuadrangularFacet::SetVertices(), and ~G4TriangularFacet().


The documentation for this class was generated from the following files: