53 std::vector<G4TwoVector> polygon,
54 std::vector<ZSection> zsections)
57 fNz(zsections.size()),
62 fGeometryType(
"G4ExtrudedSolid")
71 std::ostringstream message;
72 message <<
"Number of polygon vertices < 3 - " << pName;
73 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0002",
79 std::ostringstream message;
80 message <<
"Number of z-sides < 2 - " << pName;
81 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0002",
85 for (
G4int i=0; i<fNz-1; ++i )
87 if ( zsections[i].fZ > zsections[i+1].fZ )
89 std::ostringstream message;
90 message <<
"Z-sections have to be ordered by z value (z0 < z1 < z2...) - "
92 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0002",
95 if ( std::fabs( zsections[i+1].fZ - zsections[i].fZ ) <
kCarTolerance * 0.5 )
97 std::ostringstream message;
98 message <<
"Z-sections with the same z position are not supported - "
100 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0001",
109 for (
G4int i=0; i<fNv; ++i ) {
111 if ( j == fNv ) j = 0;
112 area += 0.5 * ( polygon[i].x()*polygon[j].y() - polygon[j].x()*polygon[i].y());
119 for (
G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[i]); }
126 for (
G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[fNv-i-1]); }
132 for (
G4int i=0; i<fNz; ++i ) { fZSections.push_back(zsections[i]); }
138 std::ostringstream message;
139 message <<
"Making facets failed - " << pName;
140 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0003",
143 fIsConvex = IsConvex();
146 ComputeProjectionParameters();
152 std::vector<G4TwoVector> polygon,
163 fGeometryType(
"G4ExtrudedSolid")
172 std::ostringstream message;
173 message <<
"Number of polygon vertices < 3 - " << pName;
174 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0002",
182 for (
G4int i=0; i<fNv; ++i )
185 if ( j == fNv ) { j = 0; }
186 area += 0.5 * ( polygon[i].x()*polygon[j].y()
187 - polygon[j].x()*polygon[i].y());
195 for (
G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[i]); }
203 for (
G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[fNv-i-1]); }
208 fZSections.push_back(
ZSection(-dz, off1, scale1));
209 fZSections.push_back(
ZSection( dz, off2, scale2));
214 std::ostringstream message;
215 message <<
"Making facets failed - " << pName;
216 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0003",
219 fIsConvex = IsConvex();
221 ComputeProjectionParameters();
228 fTriangles(), fIsConvex(false), fGeometryType(
"G4ExtrudedSolid")
238 fPolygon(rhs.fPolygon), fZSections(rhs.fZSections),
239 fTriangles(rhs.fTriangles), fIsConvex(rhs.fIsConvex),
240 fGeometryType(rhs.fGeometryType), fKScales(rhs.fKScales),
241 fScale0s(rhs.fScale0s), fKOffsets(rhs.fKOffsets), fOffset0s(rhs.fOffset0s)
252 if (
this == &rhs) {
return *
this; }
260 fNv = rhs.fNv; fNz = rhs.fNz;
261 fPolygon = rhs.fPolygon; fZSections = rhs.fZSections;
262 fTriangles = rhs.fTriangles; fIsConvex = rhs.fIsConvex;
263 fGeometryType = rhs.fGeometryType; fKScales = rhs.fKScales;
264 fScale0s = rhs.fScale0s; fKOffsets = rhs.fKOffsets;
265 fOffset0s = rhs.fOffset0s;
279 void G4ExtrudedSolid::ComputeProjectionParameters()
298 G4double kscale = (scale2 - scale1)/(z2 - z1);
299 G4double scale0 = scale2 - kscale*(z2 - z1)/2.0;
303 fKScales.push_back(kscale);
304 fScale0s.push_back(scale0);
305 fKOffsets.push_back(koff);
306 fOffset0s.push_back(off0);
318 + fZSections[iz].fOffset.x(),
319 fPolygon[ind].y() * fZSections[
iz].fScale
320 + fZSections[
iz].fOffset.y(), fZSections[
iz].fZ);
337 while ( point.
z() > fZSections[iz+1].fZ && iz < fNz-2 ) { ++
iz; }
339 G4double z0 = ( fZSections[iz+1].fZ + fZSections[
iz].fZ )/2.0;
341 G4double pscale = fKScales[
iz]*(point.
z()-
z0) + fScale0s[iz];
351 return (p2 - poffset)/pscale;
361 if ( l1.
x() == l2.
x() )
377 G4bool squareComp= (dy*dy < (1+slope*slope) * tol * tol);
399 return IsSameLine(p, l1, l2);
409 return ( (p1.
x() - l1.
x()) * (l2.
y() - l1.
y())
410 - (l2.
x() - l1.
x()) * (p1.
y() - l1.
y()) )
411 * ( (p2.
x() - l1.
x()) * (l2.
y() - l1.
y())
412 - (l2.
x() - l1.
x()) * (p2.
y() - l1.
y()) ) > 0;
425 if ( ( p.
x() < a.
x() && p.
x() < b.
x() && p.
x() < c.
x() ) ||
426 ( p.
x() > a.
x() && p.
x() > b.
x() && p.
x() > c.
x() ) ||
427 ( p.
y() < a.
y() && p.
y() < b.
y() && p.
y() < c.
y() ) ||
428 ( p.
y() > a.
y() && p.
y() > b.
y() && p.
y() > c.
y() ) )
return false;
431 = IsSameSide(p, a, b, c)
432 && IsSameSide(p, b, a, c)
433 && IsSameSide(p, c, a, b);
436 = IsSameLineSegment(p, a, b)
437 || IsSameLineSegment(p, b, c)
438 || IsSameLineSegment(p, c, a);
440 return inside || onEdge;
455 if ( result < 0 ) result += 2*
pi;
463 G4ExtrudedSolid::MakeDownFacet(
G4int ind1,
G4int ind2,
G4int ind3)
const
468 std::vector<G4ThreeVector> vertices;
476 = (vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
478 if ( cross.
z() > 0.0 )
486 vertices[1] = vertices[2];
502 std::vector<G4ThreeVector> vertices;
503 vertices.push_back(
GetVertex(fNz-1, ind1));
504 vertices.push_back(
GetVertex(fNz-1, ind2));
505 vertices.push_back(
GetVertex(fNz-1, ind3));
510 = (vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
512 if ( cross.
z() < 0.0 )
520 vertices[1] = vertices[2];
530 G4bool G4ExtrudedSolid::AddGeneralPolygonFacets()
534 typedef std::pair < G4TwoVector, G4int > Vertex;
538 std::vector< Vertex > verticesToBeDone;
539 for (
G4int i=0; i<fNv; ++i )
541 verticesToBeDone.push_back(Vertex(fPolygon[i], i));
543 std::vector< Vertex > ears;
545 std::vector< Vertex >::iterator
c1 = verticesToBeDone.begin();
546 std::vector< Vertex >::iterator c2 = c1+1;
547 std::vector< Vertex >::iterator c3 = c1+2;
548 while ( verticesToBeDone.size()>2 )
557 G4double angle = GetAngle(c2->first, c3->first, c1->first);
570 if ( c3 == verticesToBeDone.end() ) { c3 = verticesToBeDone.begin(); }
576 angle = GetAngle(c2->first, c3->first, c1->first);
581 if ( counter > fNv) {
582 G4Exception(
"G4ExtrudedSolid::AddGeneralPolygonFacets",
584 "Triangularisation has failed.");
590 std::vector< Vertex >::iterator it;
591 for ( it=verticesToBeDone.begin(); it != verticesToBeDone.end(); ++it )
595 if ( it == c1 || it == c2 || it == c3 ) {
continue; }
597 if ( IsPointInside(c1->first, c2->first, c3->first, it->first) )
607 if ( c3 == verticesToBeDone.end() ) { c3 = verticesToBeDone.begin(); }
622 result =
AddFacet( MakeDownFacet(c1->second, c2->second, c3->second) );
623 if ( ! result ) {
return false; }
625 result =
AddFacet( MakeUpFacet(c1->second, c2->second, c3->second) );
626 if ( ! result ) {
return false; }
628 std::vector<G4int> triangle(3);
629 triangle[0] = c1->second;
630 triangle[1] = c2->second;
631 triangle[2] = c3->second;
632 fTriangles.push_back(triangle);
636 verticesToBeDone.erase(c2);
637 c1 = verticesToBeDone.begin();
647 G4bool G4ExtrudedSolid::MakeFacets()
659 if ( ! good ) {
return false; }
663 if ( ! good ) {
return false; }
665 std::vector<G4int> triangle(3);
669 fTriangles.push_back(triangle);
677 if ( ! good ) {
return false; }
682 if ( ! good ) {
return false; }
684 std::vector<G4int> triangle1(3);
688 fTriangles.push_back(triangle1);
690 std::vector<G4int> triangle2(3);
694 fTriangles.push_back(triangle2);
698 good = AddGeneralPolygonFacets();
699 if ( ! good ) {
return false; }
704 for (
G4int iz = 0; iz < fNz-1; ++
iz )
706 for (
G4int i = 0; i < fNv; ++i )
708 G4int j = (i+1) % fNv;
712 if ( ! good ) {
return false; }
723 G4bool G4ExtrudedSolid::IsConvex()
const
727 for (
G4int i=0; i< fNv; ++i )
729 G4int j = ( i + 1 ) % fNv;
730 G4int k = ( i + 2 ) % fNv;
734 if ( dphi < 0. ) { dphi += 2.*
pi; }
736 if ( dphi >=
pi ) {
return false; }
748 return fGeometryType;
785 for (
G4int i=0; i<fNv; ++i )
787 G4int j = (i+1) % fNv;
788 if ( IsSameLineSegment(pscaled, fPolygon[i], fPolygon[j]) )
799 std::vector< std::vector<G4int> >::const_iterator it = fTriangles.begin();
803 if ( IsPointInside(fPolygon[(*it)[0]], fPolygon[(*it)[1]],
804 fPolygon[(*it)[2]], pscaled) ) { inside =
true; }
806 }
while ( (inside ==
false) && (it != fTriangles.end()) );
812 if ( std::fabs( p.
z() - fZSections[0].fZ ) <
kCarTolerance * 0.5 ||
844 if (validNorm) { *validNorm = fIsConvex; }
863 G4int oldprc = os.precision(16);
864 os <<
"-----------------------------------------------------------\n"
865 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
866 <<
" ===================================================\n"
867 <<
" Solid geometry type: " << fGeometryType <<
G4endl;
870 { os <<
" Convex polygon; list of vertices:" <<
G4endl; }
872 { os <<
" Concave polygon; list of vertices:" <<
G4endl; }
874 for (
G4int i=0; i<fNv; ++i )
876 os << std::setw(5) <<
"#" << i
877 <<
" vx = " << fPolygon[i].x()/
mm <<
" mm"
878 <<
" vy = " << fPolygon[i].y()/
mm <<
" mm" <<
G4endl;
881 os <<
" Sections:" <<
G4endl;
882 for (
G4int iz=0; iz<fNz; ++
iz )
884 os <<
" z = " << fZSections[
iz].fZ/
mm <<
" mm "
885 <<
" x0= " << fZSections[
iz].fOffset.x()/
mm <<
" mm "
886 <<
" y0= " << fZSections[
iz].fOffset.y()/
mm <<
" mm "
887 <<
" scale= " << fZSections[
iz].fScale <<
G4endl;
905 os.precision(oldprc);
G4GeometryType GetEntityType() const
void SetSolidClosed(const G4bool t)
CLHEP::Hep3Vector G4ThreeVector
G4ExtrudedSolid(const G4String &pName, std::vector< G4TwoVector > polygon, std::vector< ZSection > zsections)
G4double GetMaxXExtent() const
G4double G4NeutronHPJENDLHEData::G4double result
G4double GetMinXExtent() const
virtual ~G4ExtrudedSolid()
EInside Inside(const G4ThreeVector &p) const
G4double GetMaxZExtent() const
G4TwoVector GetVertex(G4int index) const
G4double GetMinZExtent() const
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::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
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double GetMaxYExtent() const
G4ExtrudedSolid & operator=(const G4ExtrudedSolid &rhs)
G4double GetMinYExtent() const
G4TessellatedSolid & operator=(const G4TessellatedSolid &right)