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

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 143 of file G4TriangularFacet.cc.

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

Referenced by GetClone(), and GetFlippedFacet().

◆ ~G4TriangularFacet()

G4TriangularFacet::~G4TriangularFacet ( )
override

Definition at line 162 of file G4TriangularFacet.cc.

163{
164 SetVertices(nullptr);
165}
void SetVertices(std::vector< G4ThreeVector > *v) override

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

◆ G4TriangularFacet() [3/4]

G4TriangularFacet::G4TriangularFacet ( const G4TriangularFacet & right)

Definition at line 201 of file G4TriangularFacet.cc.

202 : G4VFacet(rhs)
203{
204 CopyFrom(rhs);
205}

◆ G4TriangularFacet() [4/4]

G4TriangularFacet::G4TriangularFacet ( G4TriangularFacet && right)
noexcept

Definition at line 209 of file G4TriangularFacet.cc.

210 : G4VFacet(rhs)
211{
212 MoveFrom(rhs);
213}

Member Function Documentation

◆ AllocatedMemory()

G4int G4TriangularFacet::AllocatedMemory ( )
inlineoverridevirtual

Implements G4VFacet.

Definition at line 163 of file G4TriangularFacet.hh.

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

◆ Distance() [1/3]

G4ThreeVector G4TriangularFacet::Distance ( const G4ThreeVector & p)

Definition at line 288 of file G4TriangularFacet.cc.

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

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

◆ Distance() [2/3]

G4double G4TriangularFacet::Distance ( const G4ThreeVector & p,
G4double minDist )
overridevirtual

Implements G4VFacet.

Definition at line 472 of file G4TriangularFacet.cc.

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

◆ Distance() [3/3]

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

Implements G4VFacet.

Definition at line 507 of file G4TriangularFacet.cc.

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

◆ Extent()

G4double G4TriangularFacet::Extent ( const G4ThreeVector axis)
overridevirtual

Implements G4VFacet.

Definition at line 548 of file G4TriangularFacet.cc.

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

◆ GetArea()

G4double G4TriangularFacet::GetArea ( ) const
overridevirtual

Implements G4VFacet.

Definition at line 791 of file G4TriangularFacet.cc.

792{
793 return fArea;
794}

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

◆ GetCircumcentre()

G4ThreeVector G4TriangularFacet::GetCircumcentre ( ) const
inlineoverridevirtual

Implements G4VFacet.

Definition at line 153 of file G4TriangularFacet.hh.

154{
155 return fCircumcentre;
156}

◆ GetClone()

G4VFacet * G4TriangularFacet::GetClone ( )
overridevirtual

Implements G4VFacet.

Definition at line 253 of file G4TriangularFacet.cc.

254{
255 auto fc = new G4TriangularFacet (GetVertex(0), GetVertex(1),
256 GetVertex(2), ABSOLUTE);
257 return fc;
258}

◆ GetEntityType()

G4GeometryType G4TriangularFacet::GetEntityType ( ) const
overridevirtual

Implements G4VFacet.

Definition at line 798 of file G4TriangularFacet.cc.

799{
800 return "G4TriangularFacet";
801}

◆ GetFlippedFacet()

G4TriangularFacet * G4TriangularFacet::GetFlippedFacet ( )

Definition at line 267 of file G4TriangularFacet.cc.

268{
269 auto flipped = new G4TriangularFacet (GetVertex(0), GetVertex(1),
270 GetVertex(2), ABSOLUTE);
271 return flipped;
272}

◆ GetNumberOfVertices()

G4int G4TriangularFacet::GetNumberOfVertices ( ) const
inlineoverridevirtual

Implements G4VFacet.

Definition at line 137 of file G4TriangularFacet.hh.

138{
139 return 3;
140}

Referenced by AllocatedMemory().

◆ GetPointOnFace()

G4ThreeVector G4TriangularFacet::GetPointOnFace ( ) const
overridevirtual

Implements G4VFacet.

Definition at line 777 of file G4TriangularFacet.cc.

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

Referenced by G4QuadrangularFacet::GetPointOnFace().

◆ GetRadius()

G4double G4TriangularFacet::GetRadius ( ) const
inlineoverridevirtual

Implements G4VFacet.

Definition at line 158 of file G4TriangularFacet.hh.

159{
160 return fRadius;
161}

◆ GetSurfaceNormal()

G4ThreeVector G4TriangularFacet::GetSurfaceNormal ( ) const
overridevirtual

Implements G4VFacet.

Definition at line 805 of file G4TriangularFacet.cc.

806{
807 return fSurfaceNormal;
808}

Referenced by G4QuadrangularFacet::GetSurfaceNormal().

◆ GetVertex()

G4ThreeVector G4TriangularFacet::GetVertex ( G4int i) const
inlineoverridevirtual

Implements G4VFacet.

Definition at line 142 of file G4TriangularFacet.hh.

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

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

◆ GetVertexIndex()

G4int G4TriangularFacet::GetVertexIndex ( G4int i) const
inlineoverridevirtual

Implements G4VFacet.

Definition at line 170 of file G4TriangularFacet.hh.

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

◆ Intersect()

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

Implements G4VFacet.

Definition at line 585 of file G4TriangularFacet.cc.

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

Implements G4VFacet.

Definition at line 132 of file G4TriangularFacet.hh.

133{
134 return fIsDefined;
135}

Referenced by G4QuadrangularFacet::IsDefined().

◆ operator=() [1/2]

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

Definition at line 218 of file G4TriangularFacet.cc.

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

◆ operator=() [2/2]

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

Definition at line 234 of file G4TriangularFacet.cc.

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

◆ SetSurfaceNormal()

void G4TriangularFacet::SetSurfaceNormal ( const G4ThreeVector & normal)

Definition at line 812 of file G4TriangularFacet.cc.

813{
814 fSurfaceNormal = normal;
815}

Referenced by G4QuadrangularFacet::G4QuadrangularFacet().

◆ SetVertex()

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

Implements G4VFacet.

Definition at line 148 of file G4TriangularFacet.hh.

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

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

◆ SetVertexIndex()

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

Implements G4VFacet.

Definition at line 176 of file G4TriangularFacet.hh.

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

◆ SetVertices()

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

Implements G4VFacet.

Definition at line 181 of file G4TriangularFacet.hh.

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

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


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