G4TriangularFacet Class Reference

#include <G4TriangularFacet.hh>

Inheritance diagram for G4TriangularFacet:

G4VFacet

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)

Detailed Description

Definition at line 65 of file G4TriangularFacet.hh.


Constructor & Destructor Documentation

G4TriangularFacet::G4TriangularFacet (  ) 

Definition at line 150 of file G4TriangularFacet.cc.

References SetVertex().

Referenced by GetClone(), and GetFlippedFacet().

00151   : fSqrDist(0.)
00152 {
00153   fVertices = new vector<G4ThreeVector>(3);
00154   G4ThreeVector zero(0,0,0);
00155   SetVertex(0, zero);
00156   SetVertex(1, zero);
00157   SetVertex(2, zero);
00158   for (G4int i = 0; i < 3; ++i) fIndices[i] = -1;
00159   fIsDefined = false;
00160   fSurfaceNormal.set(0,0,0);
00161   fA = fB = fC = 0;
00162   fE1 = zero;
00163   fE2 = zero;
00164   fDet = 0.0;
00165   fArea = fRadius = 0.0;
00166 }

G4TriangularFacet::~G4TriangularFacet (  ) 

Definition at line 170 of file G4TriangularFacet.cc.

References SetVertices().

00171 {
00172   SetVertices(0);
00173 }

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

Definition at line 75 of file G4TriangularFacet.cc.

References ABSOLUTE, G4endl, G4Exception(), GetVertex(), JustWarning, G4VFacet::kCarTolerance, and SetVertex().

00079   : fSqrDist(0.)
00080 {
00081   fVertices = new vector<G4ThreeVector>(3);
00082 
00083   SetVertex(0, vt0);
00084   if (vertexType == ABSOLUTE)
00085   {
00086     SetVertex(1, vt1);
00087     SetVertex(2, vt2);
00088     fE1 = vt1 - vt0;
00089     fE2 = vt2 - vt0;
00090   }
00091   else
00092   {
00093     SetVertex(1, vt0 + vt1);
00094     SetVertex(2, vt0 + vt2);
00095     fE1 = vt1;
00096     fE2 = vt2;
00097   }
00098 
00099   for (G4int i = 0; i < 3; ++i) fIndices[i] = -1;
00100 
00101   G4double eMag1 = fE1.mag();
00102   G4double eMag2 = fE2.mag();
00103   G4double eMag3 = (fE2-fE1).mag();
00104 
00105   if (eMag1 <= kCarTolerance || eMag2 <= kCarTolerance
00106                              || eMag3 <= kCarTolerance)
00107   {
00108     ostringstream message;
00109     message << "Length of sides of facet are too small." << G4endl
00110       << "fVertices[0] = " << GetVertex(0) << G4endl
00111       << "fVertices[1] = " << GetVertex(1) << G4endl
00112       << "fVertices[2] = " << GetVertex(2) << G4endl
00113       << "Side lengths = fVertices[0]->fVertices[1]" << eMag1 << G4endl
00114       << "Side lengths = fVertices[0]->fVertices[2]" << eMag2 << G4endl
00115       << "Side lengths = fVertices[1]->fVertices[2]" << eMag3;
00116     G4Exception("G4TriangularFacet::G4TriangularFacet()",
00117                 "GeomSolids1001", JustWarning, message);
00118     fIsDefined     = false;
00119     fSurfaceNormal.set(0,0,0);
00120     fA = fB = fC = 0.0;
00121     fDet = 0.0;
00122     fArea = fRadius = 0.0;
00123   }
00124   else
00125   { 
00126     fIsDefined     = true;
00127     fSurfaceNormal = fE1.cross(fE2).unit();
00128     fA   = fE1.mag2();
00129     fB   = fE1.dot(fE2);
00130     fC   = fE2.mag2();
00131     fDet = fabs(fA*fC - fB*fB);
00132 
00133     //    sMin = -0.5*kCarTolerance/sqrt(fA);
00134     //    sMax = 1.0 - sMin;
00135     //    tMin = -0.5*kCarTolerance/sqrt(fC);
00136     //    G4ThreeVector vtmp = 0.25 * (fE1 + fE2);
00137 
00138     fArea = 0.5 * (fE1.cross(fE2)).mag();
00139     G4double lambda0 = (fA-fB) * fC / (8.0*fArea*fArea);
00140     G4double lambda1 = (fC-fB) * fA / (8.0*fArea*fArea);
00141     G4ThreeVector p0 = GetVertex(0);
00142     fCircumcentre = p0 + lambda0*fE1 + lambda1*fE2;
00143     G4double radiusSqr = (fCircumcentre-p0).mag2();
00144     fRadius = sqrt(radiusSqr);
00145   }
00146 }

G4TriangularFacet::G4TriangularFacet ( const G4TriangularFacet right  ) 

Definition at line 191 of file G4TriangularFacet.cc.

00192   : G4VFacet(rhs)
00193 {
00194   CopyFrom(rhs);
00195 }


Member Function Documentation

G4int G4TriangularFacet::AllocatedMemory (  )  [inline, virtual]

Implements G4VFacet.

Definition at line 165 of file G4TriangularFacet.hh.

References GetNumberOfVertices().

00166 {
00167   G4int size = sizeof(*this);
00168   size += GetNumberOfVertices() * sizeof(G4ThreeVector);
00169   return size;
00170 }

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

Implements G4VFacet.

Definition at line 470 of file G4TriangularFacet.cc.

References Distance(), and G4VFacet::kCarTolerance.

00473 {
00474   //
00475   // Start with quicky test to determine if the surface of the sphere enclosing
00476   // the triangle is any closer to p than minDist.  If not, then don't bother
00477   // about more accurate test.
00478   //
00479   G4double dist = kInfinity;
00480   if ((p-fCircumcentre).mag()-fRadius < minDist)
00481   {
00482     //
00483     // It's possible that the triangle is closer than minDist,
00484     // so do more accurate assessment.
00485     //
00486     G4ThreeVector v  = Distance(p);
00487     G4double dist1 = sqrt(fSqrDist);
00488     G4double dir = v.dot(fSurfaceNormal);
00489     G4bool wrongSide = (dir > 0.0 && !outgoing) || (dir < 0.0 && outgoing);
00490     if (dist1 <= kCarTolerance)
00491     {
00492       //
00493       // Point p is very close to triangle.  Check if it's on the wrong side,
00494       // in which case return distance of 0.0 otherwise .
00495       //
00496       if (wrongSide) dist = 0.0;
00497       else dist = dist1;
00498     }
00499     else if (!wrongSide) dist = dist1;
00500   }
00501   return dist;
00502 }

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

Implements G4VFacet.

Definition at line 435 of file G4TriangularFacet.cc.

References Distance().

00437 {
00438   //
00439   // Start with quicky test to determine if the surface of the sphere enclosing
00440   // the triangle is any closer to p than minDist.  If not, then don't bother
00441   // about more accurate test.
00442   //
00443   G4double dist = kInfinity;
00444   if ((p-fCircumcentre).mag()-fRadius < minDist)
00445   {
00446     //
00447     // It's possible that the triangle is closer than minDist,
00448     // so do more accurate assessment.
00449     //
00450     dist = Distance(p).mag();
00451   }
00452   return dist;
00453 }

G4ThreeVector G4TriangularFacet::Distance ( const G4ThreeVector p  ) 

Definition at line 251 of file G4TriangularFacet.cc.

References GetVertex().

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

00252 {
00253   G4ThreeVector D  = GetVertex(0) - p;
00254   G4double d = fE1.dot(D);
00255   G4double e = fE2.dot(D);
00256   G4double f = D.mag2();
00257   G4double q = fB*e - fC*d;
00258   G4double t = fB*d - fA*e;
00259   fSqrDist = 0.;
00260 
00261   if (q+t <= fDet)
00262   {
00263     if (q < 0.0)
00264     {
00265       if (t < 0.0)
00266       {
00267         //
00268         // We are in region 4.
00269         //
00270         if (d < 0.0)
00271         {
00272           t = 0.0;
00273           if (-d >= fA) {q = 1.0; fSqrDist = fA + 2.0*d + f;}
00274           else         {q = -d/fA; fSqrDist = d*q + f;}
00275         }
00276         else
00277         {
00278           q = 0.0;
00279           if       (e >= 0.0) {t = 0.0; fSqrDist = f;}
00280           else if (-e >= fC)   {t = 1.0; fSqrDist = fC + 2.0*e + f;}
00281           else                {t = -e/fC; fSqrDist = e*t + f;}
00282         }
00283       }
00284       else
00285       {
00286         //
00287         // We are in region 3.
00288         //
00289         q = 0.0;
00290         if      (e >= 0.0) {t = 0.0; fSqrDist = f;}
00291         else if (-e >= fC)  {t = 1.0; fSqrDist = fC + 2.0*e + f;}
00292         else               {t = -e/fC; fSqrDist = e*t + f;}
00293       }
00294     }
00295     else if (t < 0.0)
00296     {
00297       //
00298       // We are in region 5.
00299       //
00300       t = 0.0;
00301       if      (d >= 0.0) {q = 0.0; fSqrDist = f;}
00302       else if (-d >= fA)  {q = 1.0; fSqrDist = fA + 2.0*d + f;}
00303       else               {q = -d/fA; fSqrDist = d*q + f;}
00304     }
00305     else
00306     {
00307       //
00308       // We are in region 0.
00309       //
00310       q       = q / fDet;
00311       t       = t / fDet;
00312       fSqrDist = q*(fA*q + fB*t + 2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
00313     }
00314   }
00315   else
00316   {
00317     if (q < 0.0)
00318     {
00319       //
00320       // We are in region 2.
00321       //
00322       G4double tmp0 = fB + d;
00323       G4double tmp1 = fC + e;
00324       if (tmp1 > tmp0)
00325       {
00326         G4double numer = tmp1 - tmp0;
00327         G4double denom = fA - 2.0*fB + fC;
00328         if (numer >= denom) {q = 1.0; t = 0.0; fSqrDist = fA + 2.0*d + f;}
00329         else
00330         {
00331           q       = numer/denom;
00332           t       = 1.0 - q;
00333           fSqrDist = q*(fA*q + fB*t +2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
00334         }
00335       }
00336       else
00337       {
00338         q = 0.0;
00339         if      (tmp1 <= 0.0) {t = 1.0; fSqrDist = fC + 2.0*e + f;}
00340         else if (e >= 0.0)    {t = 0.0; fSqrDist = f;}
00341         else                  {t = -e/fC; fSqrDist = e*t + f;}
00342       }
00343     }
00344     else if (t < 0.0)
00345     {
00346       //
00347       // We are in region 6.
00348       //
00349       G4double tmp0 = fB + e;
00350       G4double tmp1 = fA + d;
00351       if (tmp1 > tmp0)
00352       {
00353         G4double numer = tmp1 - tmp0;
00354         G4double denom = fA - 2.0*fB + fC;
00355         if (numer >= denom) {t = 1.0; q = 0.0; fSqrDist = fC + 2.0*e + f;}
00356         else
00357         {
00358           t       = numer/denom;
00359           q       = 1.0 - t;
00360           fSqrDist = q*(fA*q + fB*t +2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
00361         }
00362       }
00363       else
00364       {
00365         t = 0.0;
00366         if      (tmp1 <= 0.0) {q = 1.0; fSqrDist = fA + 2.0*d + f;}
00367         else if (d >= 0.0)    {q = 0.0; fSqrDist = f;}
00368         else                  {q = -d/fA; fSqrDist = d*q + f;}
00369       }
00370     }
00371     else
00372       //
00373       // We are in region 1.
00374       //
00375     {
00376       G4double numer = fC + e - fB - d;
00377       if (numer <= 0.0)
00378       {
00379         q       = 0.0;
00380         t       = 1.0;
00381         fSqrDist = fC + 2.0*e + f;
00382       }
00383       else
00384       {
00385         G4double denom = fA - 2.0*fB + fC;
00386         if (numer >= denom) {q = 1.0; t = 0.0; fSqrDist = fA + 2.0*d + f;}
00387         else
00388         {
00389           q       = numer/denom;
00390           t       = 1.0 - q;
00391           fSqrDist = q*(fA*q + fB*t + 2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
00392         }
00393       }
00394     }
00395   } 
00396   //
00397   //
00398   // Do fA check for rounding errors in the distance-squared.  It appears that
00399   // the conventional methods for calculating fSqrDist breaks down when very
00400   // near to or at the surface (as required by transport).
00401   // We'll therefore also use the magnitude-squared of the vector displacement.
00402   // (Note that I've also tried to get around this problem by using the
00403   // existing equations for
00404   //
00405   //    fSqrDist = function(fA,fB,fC,d,q,t)
00406   //
00407   // and use fA more accurate addition process which minimises errors and
00408   // breakdown of cummutitivity [where (A+B)+C != A+(B+C)] but this still
00409   // doesn't work.
00410   // Calculation from u = D + q*fE1 + t*fE2 is less efficient, but appears
00411   // more robust.
00412   //
00413   if (fSqrDist < 0.0) fSqrDist = 0.;
00414   G4ThreeVector u = D + q*fE1 + t*fE2;
00415   G4double u2 = u.mag2();
00416   //
00417   // The following (part of the roundoff error check) is from Oliver Merle'q
00418   // updates.
00419   //
00420   if (fSqrDist > u2) fSqrDist = u2;
00421 
00422   return u;
00423 }

G4double G4TriangularFacet::Extent ( const G4ThreeVector  axis  )  [virtual]

Implements G4VFacet.

Definition at line 511 of file G4TriangularFacet.cc.

References GetVertex(), and G4InuclParticleNames::sp.

00512 {
00513   G4double ss = GetVertex(0).dot(axis);
00514   G4double sp = GetVertex(1).dot(axis);
00515   if (sp > ss) ss = sp;
00516   sp = GetVertex(2).dot(axis);
00517   if (sp > ss) ss = sp;
00518   return ss;
00519 }

G4double G4TriangularFacet::GetArea (  )  [virtual]

Implements G4VFacet.

Definition at line 755 of file G4TriangularFacet.cc.

Referenced by G4QuadrangularFacet::GetArea().

00756 {
00757   return fArea;
00758 }

G4ThreeVector G4TriangularFacet::GetCircumcentre (  )  const [inline, virtual]

Implements G4VFacet.

Definition at line 155 of file G4TriangularFacet.hh.

00156 {
00157   return fCircumcentre;
00158 }

G4VFacet * G4TriangularFacet::GetClone (  )  [virtual]

Implements G4VFacet.

Definition at line 216 of file G4TriangularFacet.cc.

References ABSOLUTE, G4TriangularFacet(), and GetVertex().

00217 {
00218   G4TriangularFacet *fc =
00219     new G4TriangularFacet (GetVertex(0), GetVertex(1), GetVertex(2), ABSOLUTE);
00220   return fc;
00221 }

G4GeometryType G4TriangularFacet::GetEntityType (  )  const [virtual]

Implements G4VFacet.

Definition at line 762 of file G4TriangularFacet.cc.

00763 {
00764   return "G4TriangularFacet";
00765 }

G4TriangularFacet * G4TriangularFacet::GetFlippedFacet (  ) 

Definition at line 230 of file G4TriangularFacet.cc.

References ABSOLUTE, G4TriangularFacet(), and GetVertex().

00231 {
00232   G4TriangularFacet *flipped =
00233     new G4TriangularFacet (GetVertex(0), GetVertex(1), GetVertex(2), ABSOLUTE);
00234   return flipped;
00235 }

G4int G4TriangularFacet::GetNumberOfVertices (  )  const [inline, virtual]

Implements G4VFacet.

Definition at line 139 of file G4TriangularFacet.hh.

Referenced by AllocatedMemory().

00140 {
00141   return 3;
00142 }

G4ThreeVector G4TriangularFacet::GetPointOnFace (  )  const [virtual]

Implements G4VFacet.

Definition at line 739 of file G4TriangularFacet.cc.

References G4InuclParticleNames::alpha, and GetVertex().

Referenced by G4QuadrangularFacet::GetPointOnFace().

00740 {
00741   G4double alpha = G4RandFlat::shoot(0., 1.);
00742   G4double beta = G4RandFlat::shoot(0., 1.);
00743   G4double lambda1 = alpha*beta;
00744   G4double lambda0 = alpha-lambda1;
00745 
00746   return GetVertex(0) + lambda0*fE1 + lambda1*fE2;
00747 }

G4double G4TriangularFacet::GetRadius (  )  const [inline, virtual]

Implements G4VFacet.

Definition at line 160 of file G4TriangularFacet.hh.

00161 {
00162   return fRadius;
00163 }

G4ThreeVector G4TriangularFacet::GetSurfaceNormal (  )  const [virtual]

Implements G4VFacet.

Definition at line 769 of file G4TriangularFacet.cc.

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

00770 {
00771   return fSurfaceNormal;
00772 }

G4ThreeVector G4TriangularFacet::GetVertex ( G4int  i  )  const [inline, virtual]

Implements G4VFacet.

Definition at line 144 of file G4TriangularFacet.hh.

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

00145 {      
00146   G4int indice = fIndices[i];
00147   return indice < 0 ? (*fVertices)[i] : (*fVertices)[indice];
00148 }

G4int G4TriangularFacet::GetVertexIndex ( G4int  i  )  const [inline, virtual]

Implements G4VFacet.

Definition at line 172 of file G4TriangularFacet.hh.

00173 {
00174   if (i < 3) return fIndices[i];
00175   else       return 999999999;
00176 }

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.

References DBL_EPSILON, G4VFacet::dirTolerance, Distance(), GetVertex(), G4TessellatedGeometryAlgorithms::IntersectLineAndTriangle2D(), G4VFacet::kCarTolerance, G4InuclParticleNames::pp, and G4InuclParticleNames::s0.

Referenced by G4QuadrangularFacet::Intersect().

00554 {
00555   //
00556   // Check whether the direction of the facet is consistent with the vector v
00557   // and the need to be outgoing or ingoing.  If inconsistent, disregard and
00558   // return false.
00559   //
00560   G4double w = v.dot(fSurfaceNormal);
00561   if ((outgoing && w < -dirTolerance) || (!outgoing && w > dirTolerance))
00562   {
00563     distance = kInfinity;
00564     distFromSurface = kInfinity;
00565     normal.set(0,0,0);
00566     return false;
00567   } 
00568   //
00569   // Calculate the orthogonal distance from p to the surface containing the
00570   // triangle.  Then determine if we're on the right or wrong side of the
00571   // surface (at fA distance greater than kCarTolerance to be consistent with
00572   // "outgoing".
00573   //
00574   const G4ThreeVector &p0 = GetVertex(0);
00575   G4ThreeVector D  = p0 - p;
00576   distFromSurface  = D.dot(fSurfaceNormal);
00577   G4bool wrongSide = (outgoing && distFromSurface < -0.5*kCarTolerance) ||
00578     (!outgoing && distFromSurface >  0.5*kCarTolerance);
00579   if (wrongSide)
00580   {
00581     distance = kInfinity;
00582     distFromSurface = kInfinity;
00583     normal.set(0,0,0);
00584     return false;
00585   }
00586 
00587   wrongSide = (outgoing && distFromSurface < 0.0)
00588            || (!outgoing && distFromSurface > 0.0);
00589   if (wrongSide)
00590   {
00591     //
00592     // We're slightly on the wrong side of the surface.  Check if we're close
00593     // enough using fA precise distance calculation.
00594     //
00595     G4ThreeVector u = Distance(p);
00596     if (fSqrDist <= kCarTolerance*kCarTolerance)
00597     {
00598       //
00599       // We're very close.  Therefore return fA small negative number
00600       // to pretend we intersect.
00601       //
00602       // distance = -0.5*kCarTolerance
00603       distance = 0.0;
00604       normal = fSurfaceNormal;
00605       return true;
00606     }
00607     else
00608     {
00609       //
00610       // We're close to the surface containing the triangle, but sufficiently
00611       // far from the triangle, and on the wrong side compared to the directions
00612       // of the surface normal and v.  There is no intersection.
00613       //
00614       distance = kInfinity;
00615       distFromSurface = kInfinity;
00616       normal.set(0,0,0);
00617       return false;
00618     }
00619   }
00620   if (w < dirTolerance && w > -dirTolerance)
00621   {
00622     //
00623     // The ray is within the plane of the triangle. Project the problem into 2D
00624     // in the plane of the triangle. First try to create orthogonal unit vectors
00625     // mu and nu, where mu is fE1/|fE1|.  This is kinda like
00626     // the original algorithm due to Rickard Holmberg, but with better
00627     // mathematical justification than the original method ... however,
00628     // beware Rickard's was less time-consuming.
00629     //
00630     // Note that vprime is not fA unit vector.  We need to keep it unnormalised
00631     // since the values of distance along vprime (s0 and s1) for intersection
00632     // with the triangle will be used to determine if we cut the plane at the
00633     // same time.
00634     //
00635     G4ThreeVector mu = fE1.unit();
00636     G4ThreeVector nu = fSurfaceNormal.cross(mu);
00637     G4TwoVector pprime(p.dot(mu), p.dot(nu));
00638     G4TwoVector vprime(v.dot(mu), v.dot(nu));
00639     G4TwoVector P0prime(p0.dot(mu), p0.dot(nu));
00640     G4TwoVector E0prime(fE1.mag(), 0.0);
00641     G4TwoVector E1prime(fE2.dot(mu), fE2.dot(nu));
00642     G4TwoVector loc[2];
00643     if (G4TessellatedGeometryAlgorithms::IntersectLineAndTriangle2D(pprime,
00644                                     vprime, P0prime, E0prime, E1prime, loc))
00645     {
00646       //
00647       // There is an intersection between the line and triangle in 2D.
00648       // Now check which part of the line intersects with the plane
00649       // containing the triangle in 3D.
00650       //
00651       G4double vprimemag = vprime.mag();
00652       G4double s0        = (loc[0] - pprime).mag()/vprimemag;
00653       G4double s1        = (loc[1] - pprime).mag()/vprimemag;
00654       G4double normDist0 = fSurfaceNormal.dot(s0*v) - distFromSurface;
00655       G4double normDist1 = fSurfaceNormal.dot(s1*v) - distFromSurface;
00656 
00657       if ((normDist0 < 0.0 && normDist1 < 0.0)
00658        || (normDist0 > 0.0 && normDist1 > 0.0)
00659        || (normDist0 == 0.0 && normDist1 == 0.0) ) 
00660       {
00661         distance        = kInfinity;
00662         distFromSurface = kInfinity;
00663         normal.set(0,0,0);
00664         return false;
00665       }
00666       else
00667       {
00668         G4double dnormDist = normDist1 - normDist0;
00669         if (fabs(dnormDist) < DBL_EPSILON)
00670         {
00671           distance = s0;
00672           normal   = fSurfaceNormal;
00673           if (!outgoing) distFromSurface = -distFromSurface;
00674           return true;
00675         }
00676         else
00677         {
00678           distance = s0 - normDist0*(s1-s0)/dnormDist;
00679           normal   = fSurfaceNormal;
00680           if (!outgoing) distFromSurface = -distFromSurface;
00681           return true;
00682         }
00683       }
00684     }
00685     else
00686     {
00687       distance = kInfinity;
00688       distFromSurface = kInfinity;
00689       normal.set(0,0,0);
00690       return false;
00691     }
00692   }
00693   //
00694   //
00695   // Use conventional algorithm to determine the whether there is an
00696   // intersection.  This involves determining the point of intersection of the
00697   // line with the plane containing the triangle, and then calculating if the
00698   // point is within the triangle.
00699   //
00700   distance = distFromSurface / w;
00701   G4ThreeVector pp = p + v*distance;
00702   G4ThreeVector DD = p0 - pp;
00703   G4double d = fE1.dot(DD);
00704   G4double e = fE2.dot(DD);
00705   G4double ss = fB*e - fC*d;
00706   G4double t = fB*d - fA*e;
00707 
00708   G4double sTolerance = (fabs(fB)+ fabs(fC) + fabs(d) + fabs(e))*kCarTolerance;
00709   G4double tTolerance = (fabs(fA)+ fabs(fB) + fabs(d) + fabs(e))*kCarTolerance;
00710   G4double detTolerance = (fabs(fA)+ fabs(fC) + 2*fabs(fB) )*kCarTolerance;
00711 
00712   //if (ss < 0.0 || t < 0.0 || ss+t > fDet)
00713   if (ss < -sTolerance || t < -tTolerance || ( ss+t - fDet ) > detTolerance)
00714   {
00715     //
00716     // The intersection is outside of the triangle.
00717     //
00718     distance = distFromSurface = kInfinity;
00719     normal.set(0,0,0);
00720     return false;
00721   }
00722   else
00723   {
00724     //
00725     // There is an intersection.  Now we only need to set the surface normal.
00726     //
00727     normal = fSurfaceNormal;
00728     if (!outgoing) distFromSurface = -distFromSurface;
00729     return true;
00730   }
00731 }

G4bool G4TriangularFacet::IsDefined (  )  const [inline, virtual]

Implements G4VFacet.

Definition at line 134 of file G4TriangularFacet.hh.

Referenced by G4QuadrangularFacet::IsDefined().

00135 {
00136   return fIsDefined;
00137 }

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

Definition at line 200 of file G4TriangularFacet.cc.

References SetVertices().

00201 {
00202   SetVertices(0);
00203 
00204   if (this != &rhs)
00205     CopyFrom(rhs);
00206 
00207   return *this;
00208 }

void G4TriangularFacet::SetSurfaceNormal ( G4ThreeVector  normal  ) 

Definition at line 776 of file G4TriangularFacet.cc.

Referenced by G4QuadrangularFacet::G4QuadrangularFacet().

00777 {
00778   fSurfaceNormal = normal;
00779 }

void G4TriangularFacet::SetVertex ( G4int  i,
const G4ThreeVector val 
) [inline, virtual]

Implements G4VFacet.

Definition at line 150 of file G4TriangularFacet.hh.

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

00151 {
00152   (*fVertices)[i] = val;
00153 }

void G4TriangularFacet::SetVertexIndex ( G4int  i,
G4int  j 
) [inline, virtual]

Implements G4VFacet.

Definition at line 178 of file G4TriangularFacet.hh.

00179 {
00180   fIndices[i] = j;
00181 }

void G4TriangularFacet::SetVertices ( std::vector< G4ThreeVector > *  v  )  [inline, virtual]

Implements G4VFacet.

Definition at line 183 of file G4TriangularFacet.hh.

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

00184 {
00185   if (fIndices[0] < 0 && fVertices)
00186   {
00187     delete fVertices;
00188     fVertices = 0;
00189   }
00190   fVertices = v;
00191 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:53:34 2013 for Geant4 by  doxygen 1.4.7