35#if !defined(G4GEOM_USE_UGENERICTRAP)
60const G4int G4GenericTrap::fgkNofVertices = 8;
61const G4double G4GenericTrap::fgkTolerance = 1E-3;
66 const std::vector<G4TwoVector>& vertices )
67 :
G4VSolid(name), fDz(halfZ), fVertices(),
74 G4String errorDescription =
"InvalidSetup in \" ";
75 errorDescription += name;
76 errorDescription +=
"\"";
82 if (
G4int(vertices.size()) != fgkNofVertices )
84 G4Exception(
"G4GenericTrap::G4GenericTrap()",
"GeomSolids0002",
92 G4Exception(
"G4GenericTrap::G4GenericTrap()",
"GeomSolids0002",
98 if(CheckOrder(vertices))
100 for (
G4int i=0; i<fgkNofVertices; ++i) {fVertices.push_back(vertices[i]);}
104 for (
auto i=0; i <4; ++i) {fVertices.push_back(vertices[3-i]);}
105 for (
auto i=0; i <4; ++i) {fVertices.push_back(vertices[7-i]);}
110 for (
auto j=0; j<2; ++j)
112 for (
auto i=1; i<4; ++i)
115 length = (fVertices[k]-fVertices[k-1]).mag();
118 std::ostringstream message;
119 message <<
"Length segment is too small." <<
G4endl
120 <<
"Distance between " << fVertices[k-1] <<
" and "
121 << fVertices[k] <<
" is only " << length <<
" mm !";
122 G4Exception(
"G4GenericTrap::G4GenericTrap()",
"GeomSolids1001",
123 JustWarning, message,
"Vertices will be collapsed.");
124 fVertices[k]=fVertices[k-1];
131 for(
auto i=0; i<4; ++i) { fTwist[i]=0.; }
132 fIsTwisted = ComputeIsTwisted();
142 if ( !fIsTwisted ) { fTessellatedSolid = CreateTessellatedSolid(); }
149 :
G4VSolid(a), halfCarTolerance(0.), fDz(0.), fVertices(),
160 delete fTessellatedSolid;
167 halfCarTolerance(rhs.halfCarTolerance),
168 fDz(rhs.fDz), fVertices(rhs.fVertices),
169 fIsTwisted(rhs.fIsTwisted),
170 fMinBBoxVector(rhs.fMinBBoxVector), fMaxBBoxVector(rhs.fMaxBBoxVector),
171 fVisSubdivisions(rhs.fVisSubdivisions),
172 fSurfaceArea(rhs.fSurfaceArea), fCubicVolume(rhs.fCubicVolume)
174 for (
auto i=0; i<4; ++i) { fTwist[i] = rhs.fTwist[i]; }
176 if (rhs.fTessellatedSolid && !fIsTwisted )
177 { fTessellatedSolid = CreateTessellatedSolid(); }
187 if (
this == &rhs) {
return *
this; }
195 halfCarTolerance = rhs.halfCarTolerance;
196 fDz = rhs.fDz; fVertices = rhs.fVertices;
197 fIsTwisted = rhs.fIsTwisted; fTessellatedSolid =
nullptr;
198 fMinBBoxVector = rhs.fMinBBoxVector; fMaxBBoxVector = rhs.fMaxBBoxVector;
199 fVisSubdivisions = rhs.fVisSubdivisions;
200 fSurfaceArea = rhs.fSurfaceArea; fCubicVolume = rhs.fCubicVolume;
202 for (
auto i=0; i<4; ++i) { fTwist[i] = rhs.fTwist[i]; }
204 if (rhs.fTessellatedSolid && !fIsTwisted )
205 {
delete fTessellatedSolid; fTessellatedSolid = CreateTessellatedSolid(); }
217 const std::vector<G4TwoVector>& poly)
const
223 for (
G4int i=0; i<4; ++i)
227 cross = (p.
x()-poly[i].x())*(poly[j].y()-poly[i].y())-
228 (p.
y()-poly[i].y())*(poly[j].x()-poly[i].x());
230 len2=(poly[i]-poly[j]).mag2();
233 if(cross*cross<=len2*halfCarTolerance*halfCarTolerance)
242 if (poly[j].x() > poly[i].x())
251 if ( p.
x() > poly[iMax].x()+halfCarTolerance
252 || p.
x() < poly[iMin].x()-halfCarTolerance )
257 if (poly[j].y() > poly[i].y())
267 if ( p.
y() > poly[iMax].y()+halfCarTolerance
268 || p.
y() < poly[iMin].y()-halfCarTolerance )
273 if ( poly[iMax].x() != poly[iMin].x() )
275 test = (p.
x()-poly[iMin].x())/(poly[iMax].x()-poly[iMin].x())
276 * (poly[iMax].y()-poly[iMin].y())+poly[iMin].y();
285 if( (test>=(poly[iMin].y()-halfCarTolerance))
286 && (test<=(poly[iMax].y()+halfCarTolerance)) )
295 else if (cross<0.) {
return kOutside; }
307 if ( (std::fabs(p.
x()-poly[0].x())
308 +std::fabs(p.
y()-poly[0].y())) > halfCarTolerance )
323 if ( fTessellatedSolid )
325 return fTessellatedSolid->
Inside(p);
330 std::vector<G4TwoVector> xy;
332 if (std::fabs(p.
z()) <= fDz+halfCarTolerance)
337 for (
auto i=0; i<4; ++i)
339 xy.push_back(fVertices[i+4]+cf*( fVertices[i]-fVertices[i+4]));
342 innew=InsidePolygone(p,xy);
346 if(std::fabs(p.
z()) > fDz-halfCarTolerance) { innew=
kSurface; }
360 if ( fTessellatedSolid )
367 p0, p1, p2, r1, r2, r3, r4;
368 G4int noSurfaces = 0;
372 distz = fDz-std::fabs(p.
z());
373 if (distz < halfCarTolerance)
389 std:: vector<G4TwoVector> vertices;
391 for (
auto i=0; i<4; ++i)
393 vertices.push_back(fVertices[i+4]+cf*(fVertices[i]-fVertices[i+4]));
398 for (
G4int q=0; q<4; ++q)
409 p2=
G4ThreeVector(vertices[(q+1)%4].x(),vertices[(q+1)%4].y(),p.
z());
417 p2=
G4ThreeVector(fVertices[(q+1)%4].x(),fVertices[(q+1)%4].y(),-fDz);
421 p2=
G4ThreeVector(fVertices[(q+1)%4+4].x(),fVertices[(q+1)%4+4].y(),fDz);
424 lnorm = (p1-p0).cross(p2-p0);
425 lnorm = lnorm.
unit();
426 if(zPlusSide) { lnorm=-lnorm; }
435 G4double proj=(p-p0).dot(p2-p0)/normP;
436 if(proj<0) { proj=0; }
437 if(proj>normP) { proj=normP; }
443 r1=r1+proj*(r2-r1)/normP;
444 r3=r3+proj*(r4-r3)/normP;
451 distxy=std::fabs((p0-p).dot(lnorm));
452 if ( distxy<halfCarTolerance )
458 sumnorm=sumnorm+lnorm;
472 if ( noSurfaces == 0 )
475 G4Exception(
"G4GenericTrap::SurfaceNormal(p)",
"GeomSolids1002",
481 else if ( noSurfaces == 1 ) { ; }
482 else { sumnorm = sumnorm.unit(); }
490 const G4int ipl )
const
495 if ( fTessellatedSolid )
511 u=fVertices[i+4]+cf*(fVertices[i]-fVertices[i+4]);
512 v=fVertices[j+4]+cf*(fVertices[j]-fVertices[j+4]);
518 if (std::fabs(distz)<halfCarTolerance)
533 if ( std::fabs(p.
z()+fDz) > halfCarTolerance )
542 lnorm=-(p1-p0).cross(p2-p0);
543 if (distz>-halfCarTolerance) { lnorm=-lnorm.
unit(); }
544 else { lnorm=lnorm.
unit(); }
553 G4double proj=(p-p0).dot(p2-p0)/normP;
554 if (proj<0) { proj=0; }
555 if (proj>normP) { proj=normP; }
561 r1=r1+proj*(r2-r1)/normP;
562 r3=r3+proj*(r4-r3)/normP;
576 const G4int ipl)
const
588 xa=fVertices[ipl].x();
589 ya=fVertices[ipl].y();
590 xb=fVertices[ipl+4].x();
591 yb=fVertices[ipl+4].y();
594 xd=fVertices[4+j].x();
595 yd=fVertices[4+j].y();
612 G4double a = (dtx*v.
y()-dty*v.
x()+(tx1*ty2-tx2*ty1)*v.
z())*v.
z();
613 G4double b = dxs*v.
y()-dys*v.
x()+(dtx*p.
y()-dty*p.
x()+ty2*xs1-ty1*xs2
614 + tx1*ys2-tx2*ys1)*v.
z();
615 G4double c=dxs*p.
y()-dys*p.
x()+xs1*ys2-xs2*ys1;
626 if (q>-halfCarTolerance)
628 if (q<halfCarTolerance)
630 if (NormalToPlane(p,ipl).
dot(v)<=0)
633 {
return kInfinity; }
639 if (std::fabs(zi)<fDz)
647 zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
648 if (zi<=halfCarTolerance) {
return q; }
656 if (a>0) { q=0.5*(-b-std::sqrt(d))/a; }
657 else { q=0.5*(-b+std::sqrt(d))/a; }
661 if (q>-halfCarTolerance)
663 if(q<halfCarTolerance)
665 if (NormalToPlane(p,ipl).
dot(v)<=0)
671 if (a>0) { q=0.5*(-b+std::sqrt(d))/a; }
672 else { q=0.5*(-b-std::sqrt(d))/a; }
673 if (q<=halfCarTolerance) {
return kInfinity; }
679 if (std::fabs(zi)<fDz)
687 zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
688 if (zi<=halfCarTolerance) {
return q; }
691 if (a>0) { q=0.5*(-b+std::sqrt(d))/a; }
692 else { q=0.5*(-b-std::sqrt(d))/a; }
696 if (q>-halfCarTolerance)
698 if(q<halfCarTolerance)
700 if (NormalToPlane(p,ipl).
dot(v)<=0)
706 if (a>0) { q=0.5*(-b-std::sqrt(d))/a; }
707 else { q=0.5*(-b+std::sqrt(d))/a; }
708 if (q<=halfCarTolerance) {
return kInfinity; }
714 if (std::fabs(zi)<fDz)
722 zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
723 if (zi<=halfCarTolerance) {
return q; }
736 if ( fTessellatedSolid )
750 dist[i]=DistToPlane(p, v, i);
756 if (std::fabs(p.
z())>fDz-halfCarTolerance)
763 dist[4] = (fDz-p.
z())/v.
z();
767 dist[4] = (-fDz-p.
z())/v.
z();
769 if (dist[4]<-halfCarTolerance)
775 if(dist[4]<halfCarTolerance)
779 if (n.dot(v)<0) { dist[4]=0.; }
780 else { dist[4]=kInfinity; }
790 if (dist[i] < distmin) { distmin = dist[i]; }
793 if (distmin<halfCarTolerance) { distmin=0.; }
805 if ( fTessellatedSolid )
812 if(safz<0) { safz=0; }
818 for (iseg=0; iseg<4; ++iseg)
820 safxy = SafetyToFace(p,iseg);
821 if (safxy>safe) { safe=safxy; }
838 p1=
G4ThreeVector(fVertices[iseg].x(),fVertices[iseg].y(),-fDz);
840 norm=NormalToPlane(p,iseg);
841 safe = (p-p1).dot(norm);
862 if ( (std::fabs(xa-xc)+std::fabs(ya-yc)) < halfCarTolerance )
864 xc=fVertices[j+4].x();
865 yc=fVertices[j+4].y();
871 if ( (std::fabs(xb-xc)+std::fabs(yb-yc)) < halfCarTolerance )
878 G4double c=(xb-xa)*(yc-ya)-(xc-xa)*(yb-ya);
884 t=-(a*p.
x()+b*p.
y()+c*p.
z()+d)/t;
886 if ( (t<halfCarTolerance) && (t>-halfCarTolerance) )
911 if ( fTessellatedSolid )
913 return fTessellatedSolid->
DistanceToOut(p, v, calcNorm, validNorm, n);
918 G4bool lateral_cross =
false;
921 if (calcNorm) { *validNorm =
true; }
925 distmin=(-fDz-p.
z())/v.
z();
932 distmin = (fDz-p.
z())/v.
z();
935 else { distmin = kInfinity; }
942 for (
G4int ipl=0; ipl<4; ++ipl)
945 xa=fVertices[ipl].x();
946 ya=fVertices[ipl].y();
947 xb=fVertices[ipl+4].x();
948 yb=fVertices[ipl+4].y();
951 xd=fVertices[4+j].x();
952 yd=fVertices[4+j].y();
954 if ( ((std::fabs(xb-xd)+std::fabs(yb-yd))<halfCarTolerance)
955 || ((std::fabs(xa-xc)+std::fabs(ya-yc))<halfCarTolerance) )
957 G4double q=DistToTriangle(p,v,ipl) ;
958 if ( (q>=0) && (q<distmin) )
979 G4double a = (dtx*v.
y()-dty*v.
x()+(tx1*ty2-tx2*ty1)*v.
z())*v.
z();
980 G4double b = dxs*v.
y()-dys*v.
x()+(dtx*p.
y()-dty*p.
x()+ty2*xs1-ty1*xs2
981 + tx1*ys2-tx2*ys1)*v.
z();
982 G4double c=dxs*p.
y()-dys*p.
x()+xs1*ys2-xs2*ys1;
992 if ((q > -halfCarTolerance) && (q < distmin))
994 if (q < halfCarTolerance)
996 if (NormalToPlane(p,ipl).
dot(v)<0.) {
continue; }
1007 if (a > 0) { q=0.5*(-b-std::sqrt(d))/a; }
1008 else { q=0.5*(-b+std::sqrt(d))/a; }
1012 if (q > -halfCarTolerance )
1016 if(q < halfCarTolerance)
1018 if (NormalToPlane(p,ipl).
dot(v)<0.)
1020 if (a > 0) { q=0.5*(-b+std::sqrt(d))/a; }
1021 else { q=0.5*(-b-std::sqrt(d))/a; }
1022 if (( q > halfCarTolerance) && (q < distmin))
1025 lateral_cross =
true;
1032 lateral_cross =
true;
1038 if (a > 0) { q=0.5*(-b+std::sqrt(d))/a; }
1039 else { q=0.5*(-b-std::sqrt(d))/a; }
1043 if ((q > -halfCarTolerance) && (q < distmin))
1045 if (q < halfCarTolerance)
1047 if (NormalToPlane(p,ipl).
dot(v)<0.)
1049 if (a > 0) { q=0.5*(-b-std::sqrt(d))/a; }
1050 else { q=0.5*(-b+std::sqrt(d))/a; }
1051 if ( ( q > halfCarTolerance) && (q < distmin) )
1054 lateral_cross =
true;
1061 lateral_cross =
true;
1075 if (v.z()>0.) { i=4; }
1076 std::vector<G4TwoVector> xy;
1077 for (
G4int j=0; j<4; j++) { xy.push_back(fVertices[i+j]); }
1081 if (InsidePolygone(pt,xy)==
kOutside)
1092 if(v.z()>0) {side=kPZ;}
1103 *n=NormalToPlane(pt,0);
1106 *n=NormalToPlane(pt,1);
1109 *n=NormalToPlane(pt,2);
1112 *n=NormalToPlane(pt,3);
1122 std::ostringstream message;
1123 G4int oldprc = message.precision(16);
1124 message <<
"Undefined side for valid surface normal to solid." <<
G4endl
1126 <<
" p.x() = " << p.
x()/mm <<
" mm" <<
G4endl
1127 <<
" p.y() = " << p.
y()/mm <<
" mm" <<
G4endl
1128 <<
" p.z() = " << p.
z()/mm <<
" mm" <<
G4endl
1129 <<
"Direction:" <<
G4endl
1130 <<
" v.x() = " << v.
x() <<
G4endl
1131 <<
" v.y() = " << v.
y() <<
G4endl
1132 <<
" v.z() = " << v.
z() <<
G4endl
1133 <<
"Proposed distance :" <<
G4endl
1134 <<
" distmin = " << distmin/mm <<
" mm";
1135 message.precision(oldprc);
1136 G4Exception(
"G4GenericTrap::DistanceToOut(p,v,..)",
1142 if (distmin<halfCarTolerance) { distmin=0.; }
1153 if ( fTessellatedSolid )
1160 if (safz<0) { safz = 0; }
1165 for (
G4int iseg=0; iseg<4; ++iseg)
1167 safxy = std::fabs(SafetyToFace(p,iseg));
1168 if (safxy < safe) { safe = safxy; }
1179 pMin = GetMinimumBBox();
1180 pMax = GetMaximumBBox();
1184 if (pMin.
x() >= pMax.
x() || pMin.
y() >= pMax.
y() || pMin.
z() >= pMax.
z())
1186 std::ostringstream message;
1187 message <<
"Bad bounding box (min >= max) for solid: "
1189 <<
"\npMin = " << pMin
1190 <<
"\npMax = " << pMax;
1191 G4Exception(
"G4GenericTrap::BoundingLimits()",
"GeomMgt0001",
1217 return exist = (pMin < pMax) ?
true :
false;
1229 for (
G4int i=0; i<4; ++i)
1233 baseA[2*i].set(va.
x(),va.
y(),-dz);
1234 baseB[2*i].set(vb.
x(),vb.
y(), dz);
1236 for (
G4int i=0; i<4; ++i)
1238 G4int k1=2*i, k2=(2*i+2)%8;
1239 G4double ax = (baseA[k2].x()-baseA[k1].x());
1240 G4double ay = (baseA[k2].y()-baseA[k1].y());
1241 G4double bx = (baseB[k2].x()-baseB[k1].x());
1242 G4double by = (baseB[k2].y()-baseB[k1].y());
1244 baseA[k1+1] = (znorm < 0.0) ? baseA[k2] : baseA[k1];
1245 baseB[k1+1] = (znorm < 0.0) ? baseB[k1] : baseB[k2];
1248 std::vector<const G4ThreeVectorList *> polygons(2);
1249 polygons[0] = &baseA;
1250 polygons[1] = &baseB;
1275 G4int oldprc = os.precision(16);
1276 os <<
"-----------------------------------------------------------\n"
1277 <<
" *** Dump for solid - " <<
GetName() <<
" *** \n"
1278 <<
" =================================================== \n"
1280 <<
" half length Z: " << fDz/mm <<
" mm \n"
1281 <<
" list of vertices:\n";
1283 for (
G4int i=0; i<fgkNofVertices; ++i )
1285 os << std::setw(5) <<
"#" << i
1286 <<
" vx = " << fVertices[i].x()/mm <<
" mm"
1287 <<
" vy = " << fVertices[i].y()/mm <<
" mm" <<
G4endl;
1289 os.precision(oldprc);
1300 if ( fTessellatedSolid )
1308 G4double rand,area,chose,cf,lambda0,lambda1,alfa,beta,zp;
1311 std::vector<G4ThreeVector> vertices;
1312 for (
auto i=0; i<4; ++i)
1314 vertices.push_back(
G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz));
1316 for (
auto i=4; i<8; ++i)
1318 vertices.push_back(
G4ThreeVector(fVertices[i].x(),fVertices[i].y(),fDz));
1323 G4double Surface0=GetFaceSurfaceArea(vertices[0],vertices[1],
1324 vertices[2],vertices[3]);
1325 G4double Surface1=GetFaceSurfaceArea(vertices[0],vertices[1],
1326 vertices[5],vertices[4]);
1327 G4double Surface2=GetFaceSurfaceArea(vertices[3],vertices[0],
1328 vertices[4],vertices[7]);
1329 G4double Surface3=GetFaceSurfaceArea(vertices[2],vertices[3],
1330 vertices[7],vertices[6]);
1331 G4double Surface4=GetFaceSurfaceArea(vertices[2],vertices[1],
1332 vertices[5],vertices[6]);
1333 G4double Surface5=GetFaceSurfaceArea(vertices[4],vertices[5],
1334 vertices[6],vertices[7]);
1336 area = Surface0+Surface1+Surface2+Surface3+Surface4+Surface5;
1339 if ( ( chose < Surface0)
1340 || ( chose > (Surface0+Surface1+Surface2+Surface3+Surface4)) )
1344 if(chose < Surface0)
1347 u = fVertices[ipl]; v = fVertices[j];
1348 w = fVertices[(ipl+3)%4];
1353 u = fVertices[ipl+4]; v = fVertices[j+4];
1354 w = fVertices[(ipl+3)%4+4];
1359 lambda0=alfa-lambda1;
1362 v = u+lambda0*v+lambda1*w;
1366 if (chose < Surface0+Surface1) { ipl=0; }
1367 else if (chose < Surface0+Surface1+Surface2) { ipl=1; }
1368 else if (chose < Surface0+Surface1+Surface2+Surface3) { ipl=2; }
1372 cf = 0.5*(fDz-zp)/fDz;
1373 u = fVertices[ipl+4]+cf*( fVertices[ipl]-fVertices[ipl+4]);
1374 v = fVertices[j+4]+cf*(fVertices[j]-fVertices[j+4]);
1397 if (fSurfaceArea == 0.0)
1401 fSurfaceArea = GetFaceSurfaceArea(v0,v1,v2,v3)
1402 + GetTwistedFaceSurfaceArea(v1,v0,v4,v5)
1403 + GetTwistedFaceSurfaceArea(v2,v1,v5,v6)
1404 + GetTwistedFaceSurfaceArea(v3,v2,v6,v7)
1405 + GetTwistedFaceSurfaceArea(v0,v3,v7,v4)
1406 + GetFaceSurfaceArea(v7,v6,v5,v4);
1410 fSurfaceArea = GetFaceSurfaceArea(v0,v1,v2,v3)
1411 + GetFaceSurfaceArea(v1,v0,v4,v5)
1412 + GetFaceSurfaceArea(v2,v1,v5,v6)
1413 + GetFaceSurfaceArea(v3,v2,v6,v7)
1414 + GetFaceSurfaceArea(v0,v3,v7,v4)
1415 + GetFaceSurfaceArea(v7,v6,v5,v4);
1418 return fSurfaceArea;
1425 if (fCubicVolume == 0.0)
1444 fCubicVolume = GetFaceCubicVolume(v0,v1,v2,v3)
1445 + GetFaceCubicVolume(v1,v0,v4,v5)
1446 + GetFaceCubicVolume(v2,v1,v5,v6)
1447 + GetFaceCubicVolume(v3,v2,v6,v7)
1448 + GetFaceCubicVolume(v0,v3,v7,v4)
1449 + GetFaceCubicVolume(v7,v6,v5,v4);
1452 return fCubicVolume;
1463 return 0.5*((p2-p0).cross(p3-p1)).mag();
1469G4GenericTrap::GetTwistedFaceSurfaceArea(
const G4ThreeVector& p0,
1478 for (
G4int is = 0; is < nstep; ++is)
1486 for (
G4int it = 0; it < nstep; ++it)
1492 area += 0.5*((t2-t0).cross(t3-t1)).mag();
1507 return (((p2-p0).cross(p3-p1)).dot(p0)) / 6.;
1512G4bool G4GenericTrap::ComputeIsTwisted()
1519 G4int nv = fgkNofVertices/2;
1521 for (
G4int i=0; i<4; ++i )
1523 dx1 = fVertices[(i+1)%nv].x()-fVertices[i].x();
1524 dy1 = fVertices[(i+1)%nv].y()-fVertices[i].y();
1525 if ( (dx1 == 0) && (dy1 == 0) ) {
continue; }
1527 dx2 = fVertices[nv+(i+1)%nv].x()-fVertices[nv+i].x();
1528 dy2 = fVertices[nv+(i+1)%nv].y()-fVertices[nv+i].y();
1530 if ( dx2 == 0 && dy2 == 0 ) {
continue; }
1531 G4double twist_angle = std::fabs(dy1*dx2 - dx1*dy2);
1532 if ( twist_angle < fgkTolerance ) {
continue; }
1534 SetTwistAngle(i,twist_angle);
1538 twist_angle = std::acos( (dx1*dx2 + dy1*dy2)
1539 / (std::sqrt(dx1*dx1+dy1*dy1)
1540 * std::sqrt(dx2*dx2+dy2*dy2)) );
1544 std::ostringstream message;
1545 message <<
"Twisted Angle is bigger than 90 degrees - " <<
GetName()
1547 <<
" Potential problem of malformed Solid !" <<
G4endl
1548 <<
" TwistANGLE = " << twist_angle
1549 <<
"*rad for lateral plane N= " << i;
1550 G4Exception(
"G4GenericTrap::ComputeIsTwisted()",
"GeomSolids1002",
1560G4bool G4GenericTrap::CheckOrder(
const std::vector<G4TwoVector>& vertices)
const
1565 G4bool clockwise_order=
true;
1570 for (
G4int i=0; i<4; ++i)
1573 sum1 += vertices[i].x()*vertices[j].y() - vertices[j].x()*vertices[i].y();
1574 sum2 += vertices[i+4].x()*vertices[j+4].y()
1575 - vertices[j+4].x()*vertices[i+4].y();
1577 if (sum1*sum2 < -fgkTolerance)
1579 std::ostringstream message;
1580 message <<
"Lower/upper faces defined with opposite clockwise - "
1582 G4Exception(
"G4GenericTrap::CheckOrder()",
"GeomSolids0002",
1586 if ((sum1 > 0.)||(sum2 > 0.))
1588 std::ostringstream message;
1589 message <<
"Vertices must be defined in clockwise XY planes - "
1591 G4Exception(
"G4GenericTrap::CheckOrder()",
"GeomSolids1001",
1593 clockwise_order =
false;
1598 G4bool illegal_cross =
false;
1599 illegal_cross = IsSegCrossingZ(vertices[0],vertices[4],
1600 vertices[1],vertices[5]);
1604 illegal_cross = IsSegCrossingZ(vertices[2],vertices[6],
1605 vertices[3],vertices[7]);
1610 illegal_cross = IsSegCrossing(vertices[0],vertices[1],
1611 vertices[2],vertices[3]);
1615 illegal_cross = IsSegCrossing(vertices[0],vertices[3],
1616 vertices[1],vertices[2]);
1620 illegal_cross = IsSegCrossing(vertices[4],vertices[5],
1621 vertices[6],vertices[7]);
1625 illegal_cross = IsSegCrossing(vertices[4],vertices[7],
1626 vertices[5],vertices[6]);
1631 std::ostringstream message;
1632 message <<
"Malformed polygone with opposite sides - " <<
GetName();
1633 G4Exception(
"G4GenericTrap::CheckOrderAndSetup()",
1636 return clockwise_order;
1641void G4GenericTrap::ReorderVertices(std::vector<G4ThreeVector>& vertices)
const
1645 std::vector<G4ThreeVector> oldVertices(vertices);
1647 for (
size_t i=0; i<oldVertices.size(); ++i )
1649 vertices[i] = oldVertices[oldVertices.size()-1-i];
1663 G4double dx1,dx2,xm=0.,ym=0.,a1=0.,a2=0.,b1=0.,b2=0.;
1667 if( std::fabs(dx1) < fgkTolerance ) { stand1 =
true; }
1668 if( std::fabs(dx2) < fgkTolerance ) { stand2 =
true; }
1671 a1 = (b.
x()*a.
y()-a.
x()*b.
y())/dx1;
1676 a2 = (d.
x()*c.
y()-c.
x()*d.
y())/dx2;
1679 if (stand1 && stand2)
1683 if (std::fabs(a.
x()-c.
x())<fgkTolerance)
1687 if ( ((c.
y()-a.
y())*(c.
y()-b.
y())<-fgkTolerance)
1688 || ((d.
y()-a.
y())*(d.
y()-b.
y())<-fgkTolerance)
1689 || ((a.
y()-c.
y())*(a.
y()-d.
y())<-fgkTolerance)
1690 || ((b.
y()-c.
y())*(b.
y()-d.
y())<-fgkTolerance) ) {
return true; }
1713 if (std::fabs(b1-b2) < fgkTolerance)
1717 if (std::fabs(c.
y()-(a1+b1*c.
x())) > fgkTolerance) {
return false; }
1721 if ( ((c.
x()-a.
x())*(c.
x()-b.
x())<-fgkTolerance)
1722 || ((d.
x()-a.
x())*(d.
x()-b.
x())<-fgkTolerance)
1723 || ((a.
x()-c.
x())*(a.
x()-d.
x())<-fgkTolerance)
1724 || ((b.
x()-c.
x())*(b.
x()-d.
x())<-fgkTolerance) ) {
return true; }
1728 xm = (a1-a2)/(b2-b1);
1729 ym = (a1*b2-a2*b1)/(b2-b1);
1735 G4double check = (xm-a.
x())*(xm-b.
x())+(ym-a.
y())*(ym-b.
y());
1736 if (check > -fgkTolerance) {
return false; }
1737 check = (xm-c.
x())*(xm-d.
x())+(ym-c.
y())*(ym-d.
y());
1738 if (check > -fgkTolerance) {
return false; }
1771 (std::fabs((p4-p3).y()) <
kCarTolerance ) ) {
return false; }
1775 det = dv.
x()*v1.
y()*v2.
z()+dv.
y()*v1.
z()*v2.
x()
1776 - dv.
x()*v1.
z()*v2.
y()-dv.
y()*v1.
x()*v2.
z();
1780 temp1 = v1.
cross(v2);
1781 temp2 = (p2-p1).cross(v2);
1782 if (temp1.
dot(temp2) < 0) {
return false; }
1786 q = ((dv).cross(v2)).mag()/q;
1796G4GenericTrap::MakeDownFacet(
const std::vector<G4ThreeVector>& fromVertices,
1803 if ( (fromVertices[ind1] == fromVertices[ind2]) ||
1804 (fromVertices[ind2] == fromVertices[ind3]) ||
1805 (fromVertices[ind1] == fromVertices[ind3]) ) {
return 0; }
1807 std::vector<G4ThreeVector> vertices;
1808 vertices.push_back(fromVertices[ind1]);
1809 vertices.push_back(fromVertices[ind2]);
1810 vertices.push_back(fromVertices[ind3]);
1814 G4ThreeVector cross=(vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
1816 if ( cross.
z() > 0.0 )
1820 std::ostringstream message;
1821 message <<
"Vertices in wrong order - " <<
GetName();
1822 G4Exception(
"G4GenericTrap::MakeDownFacet",
"GeomSolids0002",
1832G4GenericTrap::MakeUpFacet(
const std::vector<G4ThreeVector>& fromVertices,
1840 if ( (fromVertices[ind1] == fromVertices[ind2]) ||
1841 (fromVertices[ind2] == fromVertices[ind3]) ||
1842 (fromVertices[ind1] == fromVertices[ind3]) ) {
return nullptr; }
1844 std::vector<G4ThreeVector> vertices;
1845 vertices.push_back(fromVertices[ind1]);
1846 vertices.push_back(fromVertices[ind2]);
1847 vertices.push_back(fromVertices[ind3]);
1851 G4ThreeVector cross=(vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
1853 if ( cross.
z() < 0.0 )
1857 std::ostringstream message;
1858 message <<
"Vertices in wrong order - " <<
GetName();
1859 G4Exception(
"G4GenericTrap::MakeUpFacet",
"GeomSolids0002",
1869G4GenericTrap::MakeSideFacet(
const G4ThreeVector& downVertex0,
1877 if ( (downVertex0 == downVertex1) && (upVertex0 == upVertex1) )
1882 if ( downVertex0 == downVertex1 )
1887 if ( upVertex0 == upVertex1 )
1902 G4int nv = fgkNofVertices/2;
1903 std::vector<G4ThreeVector> downVertices;
1904 for (
G4int i=0; i<nv; ++i )
1907 fVertices[i].y(), -fDz));
1910 std::vector<G4ThreeVector> upVertices;
1911 for (
G4int i=nv; i<2*nv; ++i )
1914 fVertices[i].y(), fDz));
1920 = (downVertices[1]-downVertices[0]).cross(downVertices[2]-downVertices[1]);
1922 = (upVertices[1]-upVertices[0]).cross(upVertices[2]-upVertices[1]);
1923 if ( (cross.
z() > 0.0) || (cross1.
z() > 0.0) )
1925 ReorderVertices(downVertices);
1926 ReorderVertices(upVertices);
1932 facet = MakeDownFacet(downVertices, 0, 1, 2);
1933 if (facet) { tessellatedSolid->
AddFacet( facet ); }
1934 facet = MakeDownFacet(downVertices, 0, 2, 3);
1935 if (facet) { tessellatedSolid->
AddFacet( facet ); }
1936 facet = MakeUpFacet(upVertices, 0, 2, 1);
1937 if (facet) { tessellatedSolid->
AddFacet( facet ); }
1938 facet = MakeUpFacet(upVertices, 0, 3, 2);
1939 if (facet) { tessellatedSolid->
AddFacet( facet ); }
1943 for (
G4int i = 0; i < nv; ++i )
1945 G4int j = (i+1) % nv;
1946 facet = MakeSideFacet(downVertices[j], downVertices[i],
1947 upVertices[i], upVertices[j]);
1949 if ( facet ) { tessellatedSolid->
AddFacet( facet ); }
1954 return tessellatedSolid;
1959void G4GenericTrap::ComputeBBox()
1964 minX = maxX = fVertices[0].x();
1965 minY = maxY = fVertices[0].y();
1967 for (
G4int i=1; i< fgkNofVertices; i++)
1969 if (minX>fVertices[i].x()) { minX=fVertices[i].x(); }
1970 if (maxX<fVertices[i].x()) { maxX=fVertices[i].x(); }
1971 if (minY>fVertices[i].y()) { minY=fVertices[i].y(); }
1972 if (maxY<fVertices[i].y()) { maxY=fVertices[i].y(); }
1984 if ( fTessellatedSolid )
2010 if ( fTessellatedSolid )
2026 if ( fTessellatedSolid !=
nullptr )
2035 minVec.
y(), maxVec.
y(),
2036 minVec.
z(), maxVec.
z());
2045 if ( fTessellatedSolid !=
nullptr )
2055 size_t nVertices, nFacets;
2057 G4int subdivisions=0;
2080 Dx = 0.5*(maxVec.
x()- minVec.
y());
2081 Dy = 0.5*(maxVec.
y()- minVec.
y());
2082 if (Dy > Dx) { Dx=Dy; }
2084 subdivisions=8*
G4int(maxTwist/(Dx*Dx*Dx)*fDz);
2085 if (subdivisions<4) { subdivisions=4; }
2086 if (subdivisions>30) { subdivisions=30; }
2089 G4int sub4=4*subdivisions;
2090 nVertices = 8+subdivisions*4;
2091 nFacets = 6+subdivisions*4;
2100 fVertices[i].y(),-fDz));
2102 for( i=0; i<subdivisions; ++i)
2104 for(
G4int j=0;j<4;j++)
2106 G4TwoVector u=fVertices[j]+cf*(i+1)*( fVertices[j+4]-fVertices[j]);
2113 fVertices[i].y(),fDz));
2119 for (i=0; i<subdivisions+1; ++i)
2122 polyhedron->
AddFacet(5+is,8+is,4+is,1+is);
2123 polyhedron->
AddFacet(8+is,7+is,3+is,4+is);
2124 polyhedron->
AddFacet(7+is,6+is,2+is,3+is);
2125 polyhedron->
AddFacet(6+is,5+is,1+is,2+is);
2127 polyhedron->
AddFacet(5+sub4,6+sub4,7+sub4,8+sub4);
std::vector< G4ThreeVector > G4ThreeVectorList
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
#define G4MUTEX_INITIALIZER
CLHEP::Hep3Vector G4ThreeVector
Hep3Vector cross(const Hep3Vector &) const
double dot(const Hep3Vector &) const
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4Polyhedron * CreatePolyhedron() const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4GenericTrap & operator=(const G4GenericTrap &rhs)
G4double GetZHalfLength() const
G4double GetCubicVolume()
G4TwoVector GetVertex(G4int index) const
G4double GetSurfaceArea()
G4bool fRebuildPolyhedron
G4ThreeVector GetPointOnSurface() const
G4double GetTwistAngle(G4int index) const
G4Polyhedron * GetPolyhedron() const
std::ostream & StreamInfo(std::ostream &os) const
G4VisExtent GetExtent() const
EInside Inside(const G4ThreeVector &p) const
G4GeometryType GetEntityType() const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
G4int GetVisSubdivisions() const
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=nullptr, G4ThreeVector *n=nullptr) const
void DescribeYourselfTo(G4VGraphicsScene &scene) const
G4GenericTrap(const G4String &name, G4double halfZ, const std::vector< G4TwoVector > &vertices)
G4Polyhedron * fpPolyhedron
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
void AddFacet(const G4int iv1, const G4int iv2, const G4int iv3, const G4int iv4=0)
void AddVertex(const G4ThreeVector &v)
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
virtual G4Polyhedron * GetPolyhedron() const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4bool AddFacet(G4VFacet *aFacet)
virtual G4double DistanceToOut(const G4ThreeVector &p) const
virtual void DescribeYourselfTo(G4VGraphicsScene &scene) const
void SetSolidClosed(const G4bool t)
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
virtual G4VisExtent GetExtent() const
virtual G4Polyhedron * CreatePolyhedron() const
virtual EInside Inside(const G4ThreeVector &p) const
virtual G4ThreeVector GetPointOnSurface() const
virtual void AddSolid(const G4Box &)=0
G4VSolid & operator=(const G4VSolid &rhs)
virtual G4double GetCubicVolume()
static G4int GetNumberOfRotationSteps()