81 fVertices =
new vector<G4ThreeVector>(3);
99 for (
G4int i = 0; i < 3; ++i) fIndices[i] = -1;
108 ostringstream message;
109 message <<
"Length of sides of facet are too small." <<
G4endl
113 <<
"Side lengths = fVertices[0]->fVertices[1]" << eMag1 <<
G4endl
114 <<
"Side lengths = fVertices[0]->fVertices[2]" << eMag2 <<
G4endl
115 <<
"Side lengths = fVertices[1]->fVertices[2]" << eMag3;
116 G4Exception(
"G4TriangularFacet::G4TriangularFacet()",
119 fSurfaceNormal.
set(0,0,0);
122 fArea = fRadius = 0.0;
131 fDet = fabs(fA*fC - fB*fB);
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);
142 fCircumcentre = p0 + lambda0*fE1 + lambda1*fE2;
143 G4double radiusSqr = (fCircumcentre-p0).mag2();
144 fRadius = sqrt(radiusSqr);
153 fVertices =
new vector<G4ThreeVector>(3);
158 for (
G4int i = 0; i < 3; ++i) fIndices[i] = -1;
160 fSurfaceNormal.
set(0,0,0);
165 fArea = fRadius = 0.0;
179 char *
p = (
char *) &rhs;
180 copy(p, p +
sizeof(*
this), (
char *)
this);
182 if (fIndices[0] < 0 && fVertices)
184 fVertices =
new vector<G4ThreeVector>(3);
185 for (
G4int i = 0; i < 3; ++i) (*fVertices)[i] = (*rhs.fVertices)[i];
273 if (-d >= fA) {q = 1.0; fSqrDist = fA + 2.0*d +
f;}
274 else {q = -d/fA; fSqrDist = d*q +
f;}
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;}
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;}
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;}
312 fSqrDist = q*(fA*q + fB*t + 2.0*
d) + t*(fB*q + fC*t + 2.0*e) +
f;
328 if (numer >= denom) {q = 1.0; t = 0.0; fSqrDist = fA + 2.0*d +
f;}
333 fSqrDist = q*(fA*q + fB*t +2.0*
d) + t*(fB*q + fC*t + 2.0*e) +
f;
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;}
355 if (numer >= denom) {t = 1.0; q = 0.0; fSqrDist = fC + 2.0*e +
f;}
360 fSqrDist = q*(fA*q + fB*t +2.0*
d) + t*(fB*q + fC*t + 2.0*e) +
f;
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;}
381 fSqrDist = fC + 2.0*e +
f;
386 if (numer >= denom) {q = 1.0; t = 0.0; fSqrDist = fA + 2.0*d +
f;}
391 fSqrDist = q*(fA*q + fB*t + 2.0*
d) + t*(fB*q + fC*t + 2.0*e) +
f;
413 if (fSqrDist < 0.0) fSqrDist = 0.;
420 if (fSqrDist > u2) fSqrDist = u2;
444 if ((p-fCircumcentre).mag()-fRadius < minDist)
480 if ((p-fCircumcentre).mag()-fRadius < minDist)
489 G4bool wrongSide = (dir > 0.0 && !outgoing) || (dir < 0.0 && outgoing);
496 if (wrongSide) dist = 0.0;
499 else if (!wrongSide) dist = dist1;
515 if (sp > ss) ss =
sp;
517 if (sp > ss) ss =
sp;
563 distance = kInfinity;
564 distFromSurface = kInfinity;
576 distFromSurface = D.
dot(fSurfaceNormal);
581 distance = kInfinity;
582 distFromSurface = kInfinity;
587 wrongSide = (outgoing && distFromSurface < 0.0)
588 || (!outgoing && distFromSurface > 0.0);
604 normal = fSurfaceNormal;
614 distance = kInfinity;
615 distFromSurface = kInfinity;
644 vprime, P0prime, E0prime, E1prime, loc))
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;
657 if ((normDist0 < 0.0 && normDist1 < 0.0)
658 || (normDist0 > 0.0 && normDist1 > 0.0)
659 || (normDist0 == 0.0 && normDist1 == 0.0) )
661 distance = kInfinity;
662 distFromSurface = kInfinity;
668 G4double dnormDist = normDist1 - normDist0;
672 normal = fSurfaceNormal;
673 if (!outgoing) distFromSurface = -distFromSurface;
678 distance = s0 - normDist0*(s1-
s0)/dnormDist;
679 normal = fSurfaceNormal;
680 if (!outgoing) distFromSurface = -distFromSurface;
687 distance = kInfinity;
688 distFromSurface = kInfinity;
700 distance = distFromSurface / w;
713 if (ss < -sTolerance || t < -tTolerance || ( ss+t - fDet ) > detTolerance)
718 distance = distFromSurface = kInfinity;
727 normal = fSurfaceNormal;
728 if (!outgoing) distFromSurface = -distFromSurface;
746 return GetVertex(0) + lambda0*fE1 + lambda1*fE2;
764 return "G4TriangularFacet";
771 return fSurfaceNormal;
778 fSurfaceNormal = normal;
void set(double x, double y, double z)
ThreeVector shoot(const G4int Ap, const G4int Af)
G4TriangularFacet & operator=(const G4TriangularFacet &right)
G4bool Intersect(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool outgoing, G4double &distance, G4double &distFromSurface, G4ThreeVector &normal)
static const G4double kCarTolerance
double dot(const Hep3Vector &) const
void SetVertices(std::vector< G4ThreeVector > *v)
void SetVertex(G4int i, const G4ThreeVector &val)
void copy(std::vector< T > &main, const std::vector< T > &data)
G4GeometryType GetEntityType() const
static const G4double dirTolerance
void SetSurfaceNormal(G4ThreeVector normal)
G4double Extent(const G4ThreeVector axis)
static G4bool IntersectLineAndTriangle2D(const G4TwoVector &p, const G4TwoVector &v, const G4TwoVector &p0, const G4TwoVector &e0, const G4TwoVector &e1, G4TwoVector location[2])
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4TriangularFacet * GetFlippedFacet()
G4ThreeVector GetVertex(G4int i) const
Hep3Vector cross(const Hep3Vector &) const
G4ThreeVector GetPointOnFace() const
G4ThreeVector GetSurfaceNormal() const
G4ThreeVector Distance(const G4ThreeVector &p)