40 #if !defined(G4GEOM_USE_UEXTRUDEDSOLID)
56 std::vector<G4TwoVector> polygon,
57 std::vector<ZSection> zsections)
60 fNz(zsections.size()),
65 fGeometryType(
"G4ExtrudedSolid")
74 std::ostringstream message;
75 message <<
"Number of polygon vertices < 3 - " << pName;
76 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0002",
82 std::ostringstream message;
83 message <<
"Number of z-sides < 2 - " << pName;
84 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0002",
90 if ( zsections[i].fZ > zsections[i+1].fZ )
92 std::ostringstream message;
93 message <<
"Z-sections have to be ordered by z value (z0 < z1 < z2...) - "
95 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0002",
98 if ( std::fabs( zsections[i+1].fZ - zsections[i].fZ ) <
kCarTolerance * 0.5 )
100 std::ostringstream message;
101 message <<
"Z-sections with the same z position are not supported - "
103 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0001",
114 if ( j == fNv ) j = 0;
115 area += 0.5 * ( polygon[i].x()*polygon[j].y() - polygon[j].x()*polygon[i].y());
141 std::ostringstream message;
142 message <<
"Making facets failed - " << pName;
143 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0003",
155 std::vector<G4TwoVector> polygon,
166 fGeometryType(
"G4ExtrudedSolid")
175 std::ostringstream message;
176 message <<
"Number of polygon vertices < 3 - " << pName;
177 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0002",
188 if ( j == fNv ) { j = 0; }
189 area += 0.5 * ( polygon[i].x()*polygon[j].y()
190 - polygon[j].x()*polygon[i].y());
217 std::ostringstream message;
218 message <<
"Making facets failed - " << pName;
219 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0003",
231 fTriangles(), fIsConvex(false), fGeometryType(
"G4ExtrudedSolid")
241 fPolygon(rhs.fPolygon), fZSections(rhs.fZSections),
242 fTriangles(rhs.fTriangles), fIsConvex(rhs.fIsConvex),
243 fGeometryType(rhs.fGeometryType), fKScales(rhs.fKScales),
244 fScale0s(rhs.fScale0s), fKOffsets(rhs.fKOffsets), fOffset0s(rhs.fOffset0s)
255 if (
this == &rhs) {
return *
this; }
301 G4double kscale = (scale2 - scale1)/(z2 - z1);
302 G4double scale0 = scale2 - kscale*(z2 - z1)/2.0;
355 return (
p2 - poffset)/pscale;
365 if ( l1.x() == l2.x() )
369 G4double slope= ((l2.y() - l1.y())/(l2.x() - l1.x()));
370 G4double predy= l1.y() + slope *(p.x() - l1.x());
381 G4bool squareComp= (dy*dy < (1+slope*slope) * tol * tol);
413 return ( (p1.x() - l1.x()) * (l2.y() - l1.y())
414 - (l2.x() - l1.x()) * (p1.y() - l1.y()) )
415 * ( (p2.x() - l1.x()) * (l2.y() - l1.y())
416 - (l2.x() - l1.x()) * (p2.y() - l1.y()) ) > 0;
429 if ( ( p.x() < a.x() && p.x() < b.x() && p.x() < c.x() ) ||
430 ( p.x() > a.x() && p.x() > b.x() && p.x() > c.x() ) ||
431 ( p.y() < a.y() && p.y() < b.y() && p.y() < c.y() ) ||
432 ( p.y() > a.y() && p.y() > b.y() && p.y() > c.y() ) )
return false;
444 return inside || onEdge;
457 G4double result = (std::atan2(t1.y(), t1.x()) - std::atan2(t2.y(), t2.x()));
459 if ( result < 0 ) result += 2*
pi;
472 std::vector<G4ThreeVector> vertices;
480 = (vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
482 if ( cross.z() > 0.0 )
490 vertices[1] = vertices[2];
506 std::vector<G4ThreeVector> vertices;
514 = (vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
516 if ( cross.z() < 0.0 )
524 vertices[1] = vertices[2];
538 typedef std::pair < G4TwoVector, G4int > Vertex;
542 std::vector< Vertex > verticesToBeDone;
545 verticesToBeDone.push_back(Vertex(
fPolygon[i], i));
547 std::vector< Vertex > ears;
549 std::vector< Vertex >::iterator
c1 = verticesToBeDone.begin();
550 std::vector< Vertex >::iterator
c2 = c1+1;
551 std::vector< Vertex >::iterator
c3 = c1+2;
552 while ( verticesToBeDone.size()>2 )
569 while ( angle >=
pi )
578 if ( c3 == verticesToBeDone.end() ) { c3 = verticesToBeDone.begin(); }
584 angle =
GetAngle(c2->first, c3->first, c1->first);
589 if ( counter > fNv) {
590 G4Exception(
"G4ExtrudedSolid::AddGeneralPolygonFacets",
592 "Triangularisation has failed.");
598 std::vector< Vertex >::iterator it;
599 for ( it=verticesToBeDone.begin(); it != verticesToBeDone.end(); ++it )
603 if ( it == c1 || it == c2 || it == c3 ) {
continue; }
605 if (
IsPointInside(c1->first, c2->first, c3->first, it->first) )
615 if ( c3 == verticesToBeDone.end() ) { c3 = verticesToBeDone.begin(); }
631 if ( ! result ) {
return false; }
634 if ( ! result ) {
return false; }
636 std::vector<G4int> triangle(3);
637 triangle[0] = c1->second;
638 triangle[1] = c2->second;
639 triangle[2] = c3->second;
644 verticesToBeDone.erase(c2);
645 c1 = verticesToBeDone.begin();
667 if ( ! good ) {
return false; }
671 if ( ! good ) {
return false; }
673 std::vector<G4int> triangle(3);
685 if ( ! good ) {
return false; }
690 if ( ! good ) {
return false; }
692 std::vector<G4int> triangle1(3);
698 std::vector<G4int> triangle2(3);
707 if ( ! good ) {
return false; }
716 G4int j = (i+1) % fNv;
720 if ( ! good ) {
return false; }
737 G4int j = ( i + 1 ) % fNv;
738 G4int k = ( i + 2 ) % fNv;
741 G4double dphi = v2.phi() - v1.phi();
742 if ( dphi < 0. ) { dphi += 2.*
pi; }
744 if ( dphi >=
pi ) {
return false; }
795 G4int j = (i+1) % fNv;
807 std::vector< std::vector<G4int> >::const_iterator it =
fTriangles.begin();
812 fPolygon[(*it)[2]], pscaled) ) { inside =
true; }
814 }
while ( (inside ==
false) && (it !=
fTriangles.end()) );
852 if (validNorm) { *validNorm =
fIsConvex; }
871 G4int oldprc = os.precision(16);
872 os <<
"-----------------------------------------------------------\n"
873 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
874 <<
" ===================================================\n"
878 { os <<
" Convex polygon; list of vertices:" <<
G4endl; }
880 { os <<
" Concave polygon; list of vertices:" <<
G4endl; }
884 os << std::setw(5) <<
"#" << i
889 os <<
" Sections:" <<
G4endl;
913 os.precision(oldprc);
static c2_factory< G4double > c2
G4bool IsPointInside(G4TwoVector a, G4TwoVector b, G4TwoVector c, G4TwoVector p) const
G4GeometryType GetEntityType() const
void SetSolidClosed(const G4bool t)
CLHEP::Hep3Vector G4ThreeVector
std::vector< G4TwoVector > fOffset0s
std::vector< G4double > fScale0s
G4ExtrudedSolid(const G4String &pName, std::vector< G4TwoVector > polygon, std::vector< ZSection > zsections)
G4double GetMaxXExtent() const
std::vector< ZSection > fZSections
G4double GetMinXExtent() const
virtual ~G4ExtrudedSolid()
EInside Inside(const G4ThreeVector &p) const
G4double GetMaxZExtent() const
G4TwoVector ProjectPoint(const G4ThreeVector &point) const
G4bool AddGeneralPolygonFacets()
G4TwoVector GetVertex(G4int index) const
G4bool IsSameLine(G4TwoVector p, G4TwoVector l1, G4TwoVector l2) const
G4double GetAngle(G4TwoVector p0, G4TwoVector pa, G4TwoVector pb) const
G4double GetMinZExtent() const
void ComputeProjectionParameters()
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
virtual G4double DistanceToOut(const G4ThreeVector &p) const
G4bool AddFacet(G4VFacet *aFacet)
std::vector< G4double > fKScales
std::ostream & StreamInfo(std::ostream &os) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4bool IsSameLineSegment(G4TwoVector p, G4TwoVector l1, G4TwoVector l2) const
std::vector< G4TwoVector > fKOffsets
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4bool IsSameSide(G4TwoVector p1, G4TwoVector p2, G4TwoVector l1, G4TwoVector l2) const
CLHEP::Hep2Vector G4TwoVector
G4VFacet * MakeUpFacet(G4int ind1, G4int ind2, G4int ind3) const
G4GeometryType fGeometryType
G4double GetMaxYExtent() const
G4ExtrudedSolid & operator=(const G4ExtrudedSolid &rhs)
G4VFacet * MakeDownFacet(G4int ind1, G4int ind2, G4int ind3) const
G4double GetMinYExtent() const
std::vector< std::vector< G4int > > fTriangles
G4TessellatedSolid & operator=(const G4TessellatedSolid &right)
std::vector< G4TwoVector > fPolygon