Geant4 9.6.0
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)
 
G4TriangularFacetoperator= (const 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 ()
 
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
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 ()=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

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

Detailed Description

Definition at line 65 of file G4TriangularFacet.hh.

Constructor & Destructor Documentation

◆ G4TriangularFacet() [1/3]

G4TriangularFacet::G4TriangularFacet ( )

Definition at line 150 of file G4TriangularFacet.cc.

151 : fSqrDist(0.)
152{
153 fVertices = new vector<G4ThreeVector>(3);
154 G4ThreeVector zero(0,0,0);
155 SetVertex(0, zero);
156 SetVertex(1, zero);
157 SetVertex(2, zero);
158 for (G4int i = 0; i < 3; ++i) fIndices[i] = -1;
159 fIsDefined = false;
160 fSurfaceNormal.set(0,0,0);
161 fA = fB = fC = 0;
162 fE1 = zero;
163 fE2 = zero;
164 fDet = 0.0;
165 fArea = fRadius = 0.0;
166}
int G4int
Definition: G4Types.hh:66
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 170 of file G4TriangularFacet.cc.

171{
172 SetVertices(0);
173}
void SetVertices(std::vector< G4ThreeVector > *v)

◆ G4TriangularFacet() [2/3]

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

Definition at line 75 of file G4TriangularFacet.cc.

79 : fSqrDist(0.)
80{
81 fVertices = new vector<G4ThreeVector>(3);
82
83 SetVertex(0, vt0);
84 if (vertexType == ABSOLUTE)
85 {
86 SetVertex(1, vt1);
87 SetVertex(2, vt2);
88 fE1 = vt1 - vt0;
89 fE2 = vt2 - vt0;
90 }
91 else
92 {
93 SetVertex(1, vt0 + vt1);
94 SetVertex(2, vt0 + vt2);
95 fE1 = vt1;
96 fE2 = vt2;
97 }
98
99 for (G4int i = 0; i < 3; ++i) fIndices[i] = -1;
100
101 G4double eMag1 = fE1.mag();
102 G4double eMag2 = fE2.mag();
103 G4double eMag3 = (fE2-fE1).mag();
104
105 if (eMag1 <= kCarTolerance || eMag2 <= kCarTolerance
106 || eMag3 <= kCarTolerance)
107 {
108 ostringstream message;
109 message << "Length of sides of facet are too small." << endl
110 << "fVertices[0] = " << GetVertex(0) << endl
111 << "fVertices[1] = " << GetVertex(1) << endl
112 << "fVertices[2] = " << GetVertex(2) << endl
113 << "Side lengths = fVertices[0]->fVertices[1]" << eMag1 << endl
114 << "Side lengths = fVertices[0]->fVertices[2]" << eMag2 << endl
115 << "Side lengths = fVertices[1]->fVertices[2]" << eMag3;
116 G4Exception("G4TriangularFacet::G4TriangularFacet()",
117 "GeomSolids1001", JustWarning, message);
118 fIsDefined = false;
119 fSurfaceNormal.set(0,0,0);
120 fA = fB = fC = 0.0;
121 fDet = 0.0;
122 fArea = fRadius = 0.0;
123 }
124 else
125 {
126 fIsDefined = true;
127 fSurfaceNormal = fE1.cross(fE2).unit();
128 fA = fE1.mag2();
129 fB = fE1.dot(fE2);
130 fC = fE2.mag2();
131 fDet = fabs(fA*fC - fB*fB);
132
133 // sMin = -0.5*kCarTolerance/sqrt(fA);
134 // sMax = 1.0 - sMin;
135 // tMin = -0.5*kCarTolerance/sqrt(fC);
136 // G4ThreeVector vtmp = 0.25 * (fE1 + fE2);
137
138 fArea = 0.5 * (fE1.cross(fE2)).mag();
139 G4double lambda0 = (fA-fB) * fC / (8.0*fArea*fArea);
140 G4double lambda1 = (fC-fB) * fA / (8.0*fArea*fArea);
141 G4ThreeVector p0 = GetVertex(0);
142 fCircumcentre = p0 + lambda0*fE1 + lambda1*fE2;
143 G4double radiusSqr = (fCircumcentre-p0).mag2();
144 fRadius = sqrt(radiusSqr);
145 }
146}
@ JustWarning
double G4double
Definition: G4Types.hh:64
@ ABSOLUTE
Definition: G4VFacet.hh:56
Hep3Vector unit() const
double mag2() const
Hep3Vector cross(const Hep3Vector &) const
double dot(const Hep3Vector &) const
double mag() const
G4ThreeVector GetVertex(G4int i) const
static const G4double kCarTolerance
Definition: G4VFacet.hh:102
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41

◆ G4TriangularFacet() [3/3]

G4TriangularFacet::G4TriangularFacet ( const G4TriangularFacet right)

Definition at line 191 of file G4TriangularFacet.cc.

192 : G4VFacet(rhs)
193{
194 CopyFrom(rhs);
195}

Member Function Documentation

◆ AllocatedMemory()

G4int G4TriangularFacet::AllocatedMemory ( )
inlinevirtual

Implements G4VFacet.

Definition at line 165 of file G4TriangularFacet.hh.

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

◆ Distance() [1/3]

G4ThreeVector G4TriangularFacet::Distance ( const G4ThreeVector p)

Definition at line 251 of file G4TriangularFacet.cc.

252{
253 G4ThreeVector D = GetVertex(0) - p;
254 G4double d = fE1.dot(D);
255 G4double e = fE2.dot(D);
256 G4double f = D.mag2();
257 G4double q = fB*e - fC*d;
258 G4double t = fB*d - fA*e;
259 fSqrDist = 0.;
260
261 if (q+t <= fDet)
262 {
263 if (q < 0.0)
264 {
265 if (t < 0.0)
266 {
267 //
268 // We are in region 4.
269 //
270 if (d < 0.0)
271 {
272 t = 0.0;
273 if (-d >= fA) {q = 1.0; fSqrDist = fA + 2.0*d + f;}
274 else {q = -d/fA; fSqrDist = d*q + f;}
275 }
276 else
277 {
278 q = 0.0;
279 if (e >= 0.0) {t = 0.0; fSqrDist = f;}
280 else if (-e >= fC) {t = 1.0; fSqrDist = fC + 2.0*e + f;}
281 else {t = -e/fC; fSqrDist = e*t + f;}
282 }
283 }
284 else
285 {
286 //
287 // We are in region 3.
288 //
289 q = 0.0;
290 if (e >= 0.0) {t = 0.0; fSqrDist = f;}
291 else if (-e >= fC) {t = 1.0; fSqrDist = fC + 2.0*e + f;}
292 else {t = -e/fC; fSqrDist = e*t + f;}
293 }
294 }
295 else if (t < 0.0)
296 {
297 //
298 // We are in region 5.
299 //
300 t = 0.0;
301 if (d >= 0.0) {q = 0.0; fSqrDist = f;}
302 else if (-d >= fA) {q = 1.0; fSqrDist = fA + 2.0*d + f;}
303 else {q = -d/fA; fSqrDist = d*q + f;}
304 }
305 else
306 {
307 //
308 // We are in region 0.
309 //
310 q = q / fDet;
311 t = t / fDet;
312 fSqrDist = q*(fA*q + fB*t + 2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
313 }
314 }
315 else
316 {
317 if (q < 0.0)
318 {
319 //
320 // We are in region 2.
321 //
322 G4double tmp0 = fB + d;
323 G4double tmp1 = fC + e;
324 if (tmp1 > tmp0)
325 {
326 G4double numer = tmp1 - tmp0;
327 G4double denom = fA - 2.0*fB + fC;
328 if (numer >= denom) {q = 1.0; t = 0.0; fSqrDist = fA + 2.0*d + f;}
329 else
330 {
331 q = numer/denom;
332 t = 1.0 - q;
333 fSqrDist = q*(fA*q + fB*t +2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
334 }
335 }
336 else
337 {
338 q = 0.0;
339 if (tmp1 <= 0.0) {t = 1.0; fSqrDist = fC + 2.0*e + f;}
340 else if (e >= 0.0) {t = 0.0; fSqrDist = f;}
341 else {t = -e/fC; fSqrDist = e*t + f;}
342 }
343 }
344 else if (t < 0.0)
345 {
346 //
347 // We are in region 6.
348 //
349 G4double tmp0 = fB + e;
350 G4double tmp1 = fA + d;
351 if (tmp1 > tmp0)
352 {
353 G4double numer = tmp1 - tmp0;
354 G4double denom = fA - 2.0*fB + fC;
355 if (numer >= denom) {t = 1.0; q = 0.0; fSqrDist = fC + 2.0*e + f;}
356 else
357 {
358 t = numer/denom;
359 q = 1.0 - t;
360 fSqrDist = q*(fA*q + fB*t +2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
361 }
362 }
363 else
364 {
365 t = 0.0;
366 if (tmp1 <= 0.0) {q = 1.0; fSqrDist = fA + 2.0*d + f;}
367 else if (d >= 0.0) {q = 0.0; fSqrDist = f;}
368 else {q = -d/fA; fSqrDist = d*q + f;}
369 }
370 }
371 else
372 //
373 // We are in region 1.
374 //
375 {
376 G4double numer = fC + e - fB - d;
377 if (numer <= 0.0)
378 {
379 q = 0.0;
380 t = 1.0;
381 fSqrDist = fC + 2.0*e + f;
382 }
383 else
384 {
385 G4double denom = fA - 2.0*fB + fC;
386 if (numer >= denom) {q = 1.0; t = 0.0; fSqrDist = fA + 2.0*d + f;}
387 else
388 {
389 q = numer/denom;
390 t = 1.0 - q;
391 fSqrDist = q*(fA*q + fB*t + 2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
392 }
393 }
394 }
395 }
396 //
397 //
398 // Do fA check for rounding errors in the distance-squared. It appears that
399 // the conventional methods for calculating fSqrDist breaks down when very
400 // near to or at the surface (as required by transport).
401 // We'll therefore also use the magnitude-squared of the vector displacement.
402 // (Note that I've also tried to get around this problem by using the
403 // existing equations for
404 //
405 // fSqrDist = function(fA,fB,fC,d,q,t)
406 //
407 // and use fA more accurate addition process which minimises errors and
408 // breakdown of cummutitivity [where (A+B)+C != A+(B+C)] but this still
409 // doesn't work.
410 // Calculation from u = D + q*fE1 + t*fE2 is less efficient, but appears
411 // more robust.
412 //
413 if (fSqrDist < 0.0) fSqrDist = 0.;
414 G4ThreeVector u = D + q*fE1 + t*fE2;
415 G4double u2 = u.mag2();
416 //
417 // The following (part of the roundoff error check) is from Oliver Merle'q
418 // updates.
419 //
420 if (fSqrDist > u2) fSqrDist = u2;
421
422 return u;
423}

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

◆ Distance() [2/3]

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

Implements G4VFacet.

Definition at line 435 of file G4TriangularFacet.cc.

437{
438 //
439 // Start with quicky test to determine if the surface of the sphere enclosing
440 // the triangle is any closer to p than minDist. If not, then don't bother
441 // about more accurate test.
442 //
443 G4double dist = kInfinity;
444 if ((p-fCircumcentre).mag()-fRadius < minDist)
445 {
446 //
447 // It's possible that the triangle is closer than minDist,
448 // so do more accurate assessment.
449 //
450 dist = Distance(p).mag();
451 }
452 return dist;
453}
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 470 of file G4TriangularFacet.cc.

473{
474 //
475 // Start with quicky test to determine if the surface of the sphere enclosing
476 // the triangle is any closer to p than minDist. If not, then don't bother
477 // about more accurate test.
478 //
479 G4double dist = kInfinity;
480 if ((p-fCircumcentre).mag()-fRadius < minDist)
481 {
482 //
483 // It's possible that the triangle is closer than minDist,
484 // so do more accurate assessment.
485 //
486 G4ThreeVector v = Distance(p);
487 G4double dist1 = sqrt(fSqrDist);
488 G4double dir = v.dot(fSurfaceNormal);
489 G4bool wrongSide = (dir > 0.0 && !outgoing) || (dir < 0.0 && outgoing);
490 if (dist1 <= kCarTolerance)
491 {
492 //
493 // Point p is very close to triangle. Check if it's on the wrong side,
494 // in which case return distance of 0.0 otherwise .
495 //
496 if (wrongSide) dist = 0.0;
497 else dist = dist1;
498 }
499 else if (!wrongSide) dist = dist1;
500 }
501 return dist;
502}
bool G4bool
Definition: G4Types.hh:67

◆ Extent()

G4double G4TriangularFacet::Extent ( const G4ThreeVector  axis)
virtual

Implements G4VFacet.

Definition at line 511 of file G4TriangularFacet.cc.

512{
513 G4double ss = GetVertex(0).dot(axis);
514 G4double sp = GetVertex(1).dot(axis);
515 if (sp > ss) ss = sp;
516 sp = GetVertex(2).dot(axis);
517 if (sp > ss) ss = sp;
518 return ss;
519}

◆ GetArea()

G4double G4TriangularFacet::GetArea ( )
virtual

Implements G4VFacet.

Definition at line 755 of file G4TriangularFacet.cc.

756{
757 return fArea;
758}

Referenced by G4QuadrangularFacet::GetArea().

◆ GetCircumcentre()

G4ThreeVector G4TriangularFacet::GetCircumcentre ( ) const
inlinevirtual

Implements G4VFacet.

Definition at line 155 of file G4TriangularFacet.hh.

156{
157 return fCircumcentre;
158}

◆ GetClone()

G4VFacet * G4TriangularFacet::GetClone ( )
virtual

Implements G4VFacet.

Definition at line 216 of file G4TriangularFacet.cc.

◆ GetEntityType()

G4GeometryType G4TriangularFacet::GetEntityType ( ) const
virtual

Implements G4VFacet.

Definition at line 762 of file G4TriangularFacet.cc.

763{
764 return "G4TriangularFacet";
765}

◆ GetFlippedFacet()

G4TriangularFacet * G4TriangularFacet::GetFlippedFacet ( )

Definition at line 230 of file G4TriangularFacet.cc.

231{
232 G4TriangularFacet *flipped =
234 return flipped;
235}

◆ GetNumberOfVertices()

G4int G4TriangularFacet::GetNumberOfVertices ( ) const
inlinevirtual

Implements G4VFacet.

Definition at line 139 of file G4TriangularFacet.hh.

140{
141 return 3;
142}

Referenced by AllocatedMemory().

◆ GetPointOnFace()

G4ThreeVector G4TriangularFacet::GetPointOnFace ( ) const
virtual

Implements G4VFacet.

Definition at line 739 of file G4TriangularFacet.cc.

740{
741 G4double alpha = G4RandFlat::shoot(0., 1.);
742 G4double beta = G4RandFlat::shoot(0., 1.);
743 G4double lambda1 = alpha*beta;
744 G4double lambda0 = alpha-lambda1;
745
746 return GetVertex(0) + lambda0*fE1 + lambda1*fE2;
747}

Referenced by G4QuadrangularFacet::GetPointOnFace().

◆ GetRadius()

G4double G4TriangularFacet::GetRadius ( ) const
inlinevirtual

Implements G4VFacet.

Definition at line 160 of file G4TriangularFacet.hh.

161{
162 return fRadius;
163}

◆ GetSurfaceNormal()

G4ThreeVector G4TriangularFacet::GetSurfaceNormal ( ) const
virtual

Implements G4VFacet.

Definition at line 769 of file G4TriangularFacet.cc.

770{
771 return fSurfaceNormal;
772}

Referenced by G4QuadrangularFacet::G4QuadrangularFacet(), and G4QuadrangularFacet::GetSurfaceNormal().

◆ GetVertex()

G4ThreeVector G4TriangularFacet::GetVertex ( G4int  i) const
inlinevirtual

Implements G4VFacet.

Definition at line 144 of file G4TriangularFacet.hh.

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

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 172 of file G4TriangularFacet.hh.

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

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

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

Referenced by G4QuadrangularFacet::Intersect().

◆ IsDefined()

G4bool G4TriangularFacet::IsDefined ( ) const
inlinevirtual

Implements G4VFacet.

Definition at line 134 of file G4TriangularFacet.hh.

135{
136 return fIsDefined;
137}

Referenced by G4QuadrangularFacet::IsDefined().

◆ operator=()

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

Definition at line 200 of file G4TriangularFacet.cc.

201{
202 SetVertices(0);
203
204 if (this != &rhs)
205 CopyFrom(rhs);
206
207 return *this;
208}

◆ SetSurfaceNormal()

void G4TriangularFacet::SetSurfaceNormal ( G4ThreeVector  normal)

Definition at line 776 of file G4TriangularFacet.cc.

777{
778 fSurfaceNormal = normal;
779}

Referenced by G4QuadrangularFacet::G4QuadrangularFacet().

◆ SetVertex()

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

Implements G4VFacet.

Definition at line 150 of file G4TriangularFacet.hh.

151{
152 (*fVertices)[i] = val;
153}

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

◆ SetVertexIndex()

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

Implements G4VFacet.

Definition at line 178 of file G4TriangularFacet.hh.

179{
180 fIndices[i] = j;
181}

◆ SetVertices()

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

Implements G4VFacet.

Definition at line 183 of file G4TriangularFacet.hh.

184{
185 if (fIndices[0] < 0 && fVertices)
186 {
187 delete fVertices;
188 fVertices = 0;
189 }
190 fVertices = v;
191}

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


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