Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4TriangularFacet Class Reference

#include <G4TriangularFacet.hh>

Inheritance diagram for G4TriangularFacet:
Collaboration 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 () 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
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

- 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::G4TriangularFacet ( )

Definition at line 167 of file G4TriangularFacet.cc.

168  : fSqrDist(0.)
169 {
170  fVertices = new vector<G4ThreeVector>(3);
171  G4ThreeVector zero(0,0,0);
172  SetVertex(0, zero);
173  SetVertex(1, zero);
174  SetVertex(2, zero);
175  for (G4int i = 0; i < 3; ++i) fIndices[i] = -1;
176  fIsDefined = false;
177  fSurfaceNormal.set(0,0,0);
178  fA = fB = fC = 0;
179  fE1 = zero;
180  fE2 = zero;
181  fDet = 0.0;
182  fArea = fRadius = 0.0;
183 }
void set(double x, double y, double z)
void SetVertex(G4int i, const G4ThreeVector &val)
int G4int
Definition: G4Types.hh:78

Here is the call graph for this function:

Here is the caller graph for this function:

G4TriangularFacet::~G4TriangularFacet ( )

Definition at line 187 of file G4TriangularFacet.cc.

188 {
189  SetVertices(0);
190 }
void SetVertices(std::vector< G4ThreeVector > *v)

Here is the call graph for this function:

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

Definition at line 79 of file G4TriangularFacet.cc.

83  : fSqrDist(0.)
84 {
85  fVertices = new vector<G4ThreeVector>(3);
86 
87  SetVertex(0, vt0);
88  if (vertexType == ABSOLUTE)
89  {
90  SetVertex(1, vt1);
91  SetVertex(2, vt2);
92  fE1 = vt1 - vt0;
93  fE2 = vt2 - vt0;
94  }
95  else
96  {
97  SetVertex(1, vt0 + vt1);
98  SetVertex(2, vt0 + vt2);
99  fE1 = vt1;
100  fE2 = vt2;
101  }
102 
103  G4ThreeVector E1xE2 = fE1.cross(fE2);
104  fArea = 0.5 * E1xE2.mag();
105  for (G4int i = 0; i < 3; ++i) fIndices[i] = -1;
106 
107  fIsDefined = true;
108  G4double delta = kCarTolerance; // Set tolerance for checking
109 
110  // Check length of edges
111  //
112  G4double leng1 = fE1.mag();
113  G4double leng2 = (fE2-fE1).mag();
114  G4double leng3 = fE2.mag();
115  if (leng1 <= delta || leng2 <= delta || leng3 <= delta)
116  {
117  fIsDefined = false;
118  }
119 
120  // Check min height of triangle
121  //
122  if (fIsDefined)
123  {
124  if (2.*fArea/std::max(std::max(leng1,leng2),leng3) <= delta)
125  {
126  fIsDefined = false;
127  }
128  }
129 
130  // Define facet
131  //
132  if (!fIsDefined)
133  {
134  ostringstream message;
135  message << "Facet is too small or too narrow." << G4endl
136  << "Triangle area = " << fArea << G4endl
137  << "P0 = " << GetVertex(0) << G4endl
138  << "P1 = " << GetVertex(1) << G4endl
139  << "P2 = " << GetVertex(2) << G4endl
140  << "Side1 length (P0->P1) = " << leng1 << G4endl
141  << "Side2 length (P1->P2) = " << leng2 << G4endl
142  << "Side3 length (P2->P0) = " << leng3;
143  G4Exception("G4TriangularFacet::G4TriangularFacet()",
144  "GeomSolids1001", JustWarning, message);
145  fSurfaceNormal.set(0,0,0);
146  fA = fB = fC = 0.0;
147  fDet = 0.0;
148  fCircumcentre = vt0 + 0.5*fE1 + 0.5*fE2;
149  fArea = fRadius = 0.0;
150  }
151  else
152  {
153  fSurfaceNormal = E1xE2.unit();
154  fA = fE1.mag2();
155  fB = fE1.dot(fE2);
156  fC = fE2.mag2();
157  fDet = std::fabs(fA*fC - fB*fB);
158 
159  fCircumcentre =
160  vt0 + (E1xE2.cross(fE1)*fC + fE2.cross(E1xE2)*fA) / (2.*E1xE2.mag2());
161  fRadius = (fCircumcentre - vt0).mag();
162  }
163 }
void set(double x, double y, double z)
static const G4double kCarTolerance
Definition: G4VFacet.hh:102
double dot(const Hep3Vector &) const
void SetVertex(G4int i, const G4ThreeVector &val)
int G4int
Definition: G4Types.hh:78
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Hep3Vector unit() const
double mag2() const
G4ThreeVector GetVertex(G4int i) const
#define G4endl
Definition: G4ios.hh:61
Hep3Vector cross(const Hep3Vector &) const
double G4double
Definition: G4Types.hh:76
double mag() const

Here is the call graph for this function:

G4TriangularFacet::G4TriangularFacet ( const G4TriangularFacet right)

Definition at line 208 of file G4TriangularFacet.cc.

209  : G4VFacet(rhs)
210 {
211  CopyFrom(rhs);
212 }

Member Function Documentation

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
int G4int
Definition: G4Types.hh:78
G4int GetNumberOfVertices() const

Here is the call graph for this function:

G4ThreeVector G4TriangularFacet::Distance ( const G4ThreeVector p)

Definition at line 268 of file G4TriangularFacet.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Implements G4VFacet.

Definition at line 452 of file G4TriangularFacet.cc.

454 {
455  //
456  // Start with quicky test to determine if the surface of the sphere enclosing
457  // the triangle is any closer to p than minDist. If not, then don't bother
458  // about more accurate test.
459  //
460  G4double dist = kInfinity;
461  if ((p-fCircumcentre).mag()-fRadius < minDist)
462  {
463  //
464  // It's possible that the triangle is closer than minDist,
465  // so do more accurate assessment.
466  //
467  dist = Distance(p).mag();
468  }
469  return dist;
470 }
static const G4double kInfinity
Definition: geomdefs.hh:42
double G4double
Definition: G4Types.hh:76
double mag() const
G4ThreeVector Distance(const G4ThreeVector &p)

Here is the call graph for this function:

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

Implements G4VFacet.

Definition at line 487 of file G4TriangularFacet.cc.

490 {
491  //
492  // Start with quicky test to determine if the surface of the sphere enclosing
493  // the triangle is any closer to p than minDist. If not, then don't bother
494  // about more accurate test.
495  //
496  G4double dist = kInfinity;
497  if ((p-fCircumcentre).mag()-fRadius < minDist)
498  {
499  //
500  // It's possible that the triangle is closer than minDist,
501  // so do more accurate assessment.
502  //
503  G4ThreeVector v = Distance(p);
504  G4double dist1 = sqrt(fSqrDist);
505  G4double dir = v.dot(fSurfaceNormal);
506  G4bool wrongSide = (dir > 0.0 && !outgoing) || (dir < 0.0 && outgoing);
507  if (dist1 <= kCarTolerance)
508  {
509  //
510  // Point p is very close to triangle. Check if it's on the wrong side,
511  // in which case return distance of 0.0 otherwise .
512  //
513  if (wrongSide) dist = 0.0;
514  else dist = dist1;
515  }
516  else if (!wrongSide) dist = dist1;
517  }
518  return dist;
519 }
static const G4double kInfinity
Definition: geomdefs.hh:42
static const G4double kCarTolerance
Definition: G4VFacet.hh:102
double dot(const Hep3Vector &) const
bool G4bool
Definition: G4Types.hh:79
tuple v
Definition: test.py:18
double G4double
Definition: G4Types.hh:76
G4ThreeVector Distance(const G4ThreeVector &p)

Here is the call graph for this function:

G4double G4TriangularFacet::Extent ( const G4ThreeVector  axis)
virtual

Implements G4VFacet.

Definition at line 528 of file G4TriangularFacet.cc.

529 {
530  G4double ss = GetVertex(0).dot(axis);
531  G4double sp = GetVertex(1).dot(axis);
532  if (sp > ss) ss = sp;
533  sp = GetVertex(2).dot(axis);
534  if (sp > ss) ss = sp;
535  return ss;
536 }
double dot(const Hep3Vector &) const
G4ThreeVector GetVertex(G4int i) const
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4TriangularFacet::GetArea ( ) const
virtual

Implements G4VFacet.

Definition at line 771 of file G4TriangularFacet.cc.

772 {
773  return fArea;
774 }

Here is the caller graph for this function:

G4ThreeVector G4TriangularFacet::GetCircumcentre ( ) const
inlinevirtual

Implements G4VFacet.

Definition at line 155 of file G4TriangularFacet.hh.

156 {
157  return fCircumcentre;
158 }
G4VFacet * G4TriangularFacet::GetClone ( )
virtual

Implements G4VFacet.

Definition at line 233 of file G4TriangularFacet.cc.

234 {
235  G4TriangularFacet *fc =
237  return fc;
238 }
G4ThreeVector GetVertex(G4int i) const

Here is the call graph for this function:

G4GeometryType G4TriangularFacet::GetEntityType ( ) const
virtual

Implements G4VFacet.

Definition at line 778 of file G4TriangularFacet.cc.

779 {
780  return "G4TriangularFacet";
781 }
G4TriangularFacet * G4TriangularFacet::GetFlippedFacet ( )

Definition at line 247 of file G4TriangularFacet.cc.

248 {
249  G4TriangularFacet *flipped =
251  return flipped;
252 }
G4ThreeVector GetVertex(G4int i) const

Here is the call graph for this function:

G4int G4TriangularFacet::GetNumberOfVertices ( ) const
inlinevirtual

Implements G4VFacet.

Definition at line 139 of file G4TriangularFacet.hh.

140 {
141  return 3;
142 }

Here is the caller graph for this function:

G4ThreeVector G4TriangularFacet::GetPointOnFace ( ) const
virtual

Implements G4VFacet.

Definition at line 757 of file G4TriangularFacet.cc.

758 {
759  G4double u = G4UniformRand();
761  if (u+v > 1.) { u = 1. - u; v = 1. - v; }
762  return GetVertex(0) + u*fE1 + v*fE2;
763 }
#define G4UniformRand()
Definition: Randomize.hh:97
tuple v
Definition: test.py:18
G4ThreeVector GetVertex(G4int i) const
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4TriangularFacet::GetRadius ( ) const
inlinevirtual

Implements G4VFacet.

Definition at line 160 of file G4TriangularFacet.hh.

161 {
162  return fRadius;
163 }
G4ThreeVector G4TriangularFacet::GetSurfaceNormal ( ) const
virtual

Implements G4VFacet.

Definition at line 785 of file G4TriangularFacet.cc.

786 {
787  return fSurfaceNormal;
788 }

Here is the caller graph for this function:

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 }
int G4int
Definition: G4Types.hh:78

Here is the caller graph for this function:

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 }
G4bool G4TriangularFacet::Intersect ( const G4ThreeVector p,
const G4ThreeVector v,
const G4bool  outgoing,
G4double distance,
G4double distFromSurface,
G4ThreeVector normal 
)
virtual

Implements G4VFacet.

Definition at line 565 of file G4TriangularFacet.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

G4bool G4TriangularFacet::IsDefined ( ) const
inlinevirtual

Implements G4VFacet.

Definition at line 134 of file G4TriangularFacet.hh.

135 {
136  return fIsDefined;
137 }

Here is the caller graph for this function:

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

Definition at line 217 of file G4TriangularFacet.cc.

218 {
219  SetVertices(0);
220 
221  if (this != &rhs)
222  CopyFrom(rhs);
223 
224  return *this;
225 }
void SetVertices(std::vector< G4ThreeVector > *v)

Here is the call graph for this function:

void G4TriangularFacet::SetSurfaceNormal ( G4ThreeVector  normal)

Definition at line 792 of file G4TriangularFacet.cc.

793 {
794  fSurfaceNormal = normal;
795 }
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77

Here is the call graph for this function:

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 }

Here is the caller graph for this function:

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

Implements G4VFacet.

Definition at line 178 of file G4TriangularFacet.hh.

179 {
180  fIndices[i] = j;
181 }
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 }
tuple v
Definition: test.py:18

Here is the caller graph for this function:


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