36 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G4GEOM_USE_PARTIAL_USOLIDS) )
45 using namespace CLHEP;
51 G4UGenericPolycone::G4UGenericPolycone(
const G4String&
name,
57 : G4USolid(name, new UGenericPolycone(name, phiStart, phiTotal, numRZ, r, z))
67 G4UGenericPolycone::G4UGenericPolycone(__void__& a)
77 G4UGenericPolycone::~G4UGenericPolycone()
86 G4UGenericPolycone::G4UGenericPolycone(
const G4UGenericPolycone &source)
97 G4UGenericPolycone::operator=(
const G4UGenericPolycone &source)
99 if (
this == &source)
return *
this;
101 G4USolid::operator=( source );
106 G4double G4UGenericPolycone::GetStartPhi()
const
108 return GetShape()->GetStartPhi();
110 G4double G4UGenericPolycone::GetEndPhi()
const
112 return GetShape()->GetEndPhi();
114 G4double G4UGenericPolycone::GetSinStartPhi()
const
116 if (!GetShape()->IsOpen())
return 0;
117 G4double phi = GetShape()->GetStartPhi();
118 return std::sin(phi);
120 G4double G4UGenericPolycone::GetCosStartPhi()
const
122 if (!GetShape()->IsOpen())
return 1;
123 G4double phi = GetShape()->GetStartPhi();
124 return std::cos(phi);
126 G4double G4UGenericPolycone::GetSinEndPhi()
const
128 if (!GetShape()->IsOpen())
return 0;
129 G4double phi = GetShape()->GetEndPhi();
130 return std::sin(phi);
132 G4double G4UGenericPolycone::GetCosEndPhi()
const
134 if (!GetShape()->IsOpen())
return 1;
135 G4double phi = GetShape()->GetEndPhi();
136 return std::cos(phi);
138 G4bool G4UGenericPolycone::IsOpen()
const
140 return GetShape()->IsOpen();
142 G4int G4UGenericPolycone::GetNumRZCorner()
const
144 return GetShape()->GetNumRZCorner();
148 UPolyconeSideRZ pside = GetShape()->GetCorner(index);
164 for (
G4int i=0; i<GetNumRZCorner(); ++i)
167 if (corner.
r < rmin) rmin = corner.
r;
168 if (corner.
r > rmax) rmax = corner.
r;
169 if (corner.
z < zmin) zmin = corner.
z;
170 if (corner.
z > zmax) zmax = corner.
z;
177 GetSinStartPhi(),GetCosStartPhi(),
178 GetSinEndPhi(),GetCosEndPhi(),
180 pMin.
set(vmin.
x(),vmin.
y(),zmin);
181 pMax.
set(vmax.
x(),vmax.
y(),zmax);
185 pMin.
set(-rmax,-rmax, zmin);
186 pMax.
set( rmax, rmax, zmax);
191 if (pMin.
x() >= pMax.
x() || pMin.
y() >= pMax.
y() || pMin.
z() >= pMax.
z())
193 std::ostringstream message;
194 message <<
"Bad bounding box (min >= max) for solid: "
196 <<
"\npMin = " << pMin
197 <<
"\npMax = " << pMax;
198 G4Exception(
"G4UGenericPolycone::Extent()",
"GeomMgt0001",
209 G4UGenericPolycone::CalculateExtent(
const EAxis pAxis,
222 if (
true)
return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
224 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
226 return exist = (pMin < pMax) ?
true :
false;
239 for (
G4int i=0; i<GetNumRZCorner(); ++i)
245 if (area < 0.) std::reverse(contourRZ.begin(),contourRZ.end());
250 std::ostringstream message;
251 message <<
"Triangulation of RZ contour has failed for solid: "
253 <<
"\nExtent has been calculated using boundary box";
254 G4Exception(
"G4UGenericPolycone::CalculateExtent()",
256 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
260 const G4int NSTEPS = 24;
266 G4int ksteps = (dphi <= astep) ? 1 : (
G4int)((dphi-
deg)/astep) + 1;
269 G4double sinHalf = std::sin(0.5*ang);
270 G4double cosHalf = std::cos(0.5*ang);
271 G4double sinStep = 2.*sinHalf*cosHalf;
272 G4double cosStep = 1. - 2.*sinHalf*sinHalf;
274 G4double sinStart = GetSinStartPhi();
275 G4double cosStart = GetCosStartPhi();
280 std::vector<const G4ThreeVectorList *> polygons;
281 polygons.resize(ksteps+2);
283 for (
G4int k=0; k<ksteps+2; ++k) pols[k].resize(6);
284 for (
G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k];
291 G4int ntria = triangles.size()/3;
292 for (
G4int i=0; i<ntria; ++i)
295 for (
G4int k=0; k<3; ++k)
297 G4int e0 = i3+k, e1 = (k<2) ? e0+1 : i3;
300 r0[k2+0] = triangles[e0].x(); z0[k2+0] = triangles[e0].y();
301 r0[k2+1] = triangles[e1].x(); z0[k2+1] = triangles[e1].y();
305 if (z0[k2+1] - z0[k2+0] <= 0)
continue;
311 G4double sinCur = sinStart*cosHalf + cosStart*sinHalf;
312 G4double cosCur = cosStart*cosHalf - sinStart*sinHalf;
313 for (
G4int j=0; j<6; ++j)
315 pols[0][j].set(r0[j]*cosStart,r0[j]*sinStart,z0[j]);
317 for (
G4int k=1; k<ksteps+1; ++k)
319 for (
G4int j=0; j<6; ++j)
321 pols[k][j].set(r1[j]*cosCur,r1[j]*sinCur,z0[j]);
324 sinCur = sinCur*cosStep + cosCur*sinStep;
325 cosCur = cosCur*cosStep - sinTmp*sinStep;
327 for (
G4int j=0; j<6; ++j)
329 pols[ksteps+1][j].set(r0[j]*cosEnd,r0[j]*sinEnd,z0[j]);
335 if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax))
continue;
336 if (emin < pMin) pMin = emin;
337 if (emax > pMax) pMax =
emax;
338 if (eminlim > pMin && emaxlim < pMax)
return true;
340 return (pMin < pMax);
347 G4Polyhedron* G4UGenericPolycone::CreatePolyhedron()
const
373 const G4int numSide =
375 * (GetEndPhi() - GetStartPhi()) /
twopi) + 1;
380 typedef G4int int4[4];
387 std::vector<G4bool> chopped(GetNumRZCorner(),
false);
388 std::vector<G4int*> triQuads;
389 G4int remaining = GetNumRZCorner();
391 while (remaining >= 3)
396 G4int iStepper = iStarter;
399 if (A < 0) { A = iStepper; }
400 else if (
B < 0) {
B = iStepper; }
401 else if (
C < 0) {
C = iStepper; }
404 if (++iStepper >= GetNumRZCorner()) { iStepper = 0; }
406 while (chopped[iStepper]);
408 while (
C < 0 && iStepper != iStarter);
413 G4double BAr = GetCorner(A).r - GetCorner(
B).r;
414 G4double BAz = GetCorner(A).z - GetCorner(
B).z;
415 G4double BCr = GetCorner(
C).r - GetCorner(
B).r;
416 G4double BCz = GetCorner(
C).z - GetCorner(
B).z;
423 triQuads.push_back(tq);
431 if (++iStarter >= GetNumRZCorner()) { iStarter = 0; }
433 while (chopped[iStarter]);
438 nNodes = (numSide + 1) * GetNumRZCorner();
439 nFaces = numSide * GetNumRZCorner() + 2 * triQuads.size();
440 faces_vec =
new int4[nFaces];
442 G4int addition = GetNumRZCorner() * numSide;
443 G4int d = GetNumRZCorner() - 1;
444 for (
G4int iEnd = 0; iEnd < 2; ++iEnd)
446 for (
size_t i = 0; i < triQuads.size(); ++i)
459 a = triQuads[i][0] + addition;
460 b = triQuads[i][2] + addition;
461 c = triQuads[i][1] + addition;
464 G4int bc = std::abs(c - b);
465 G4int ca = std::abs(a - c);
466 faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
467 faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
468 faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
469 faces_vec[iface][3] = 0;
476 xyz =
new double3[nNodes];
477 const G4double dPhi = (GetEndPhi() - GetStartPhi()) / numSide;
480 for (
G4int iSide = 0; iSide < numSide; ++iSide)
482 for (
G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
484 xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
485 xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
486 xyz[ixyz][2] = GetCorner(iCorner).z;
489 if (iCorner < GetNumRZCorner() - 1)
491 faces_vec[iface][0] = ixyz + 1;
492 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
493 faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
494 faces_vec[iface][3] = ixyz + 2;
498 faces_vec[iface][0] = ixyz + 1;
499 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
500 faces_vec[iface][2] = ixyz + 2;
501 faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
504 else if (iSide == numSide - 1)
506 if (iCorner < GetNumRZCorner() - 1)
508 faces_vec[iface][0] = ixyz + 1;
509 faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
510 faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
511 faces_vec[iface][3] = -(ixyz + 2);
515 faces_vec[iface][0] = ixyz + 1;
516 faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
517 faces_vec[iface][2] = ixyz + 2;
518 faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
523 if (iCorner < GetNumRZCorner() - 1)
525 faces_vec[iface][0] = ixyz + 1;
526 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
527 faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
528 faces_vec[iface][3] = -(ixyz + 2);
532 faces_vec[iface][0] = ixyz + 1;
533 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
534 faces_vec[iface][2] = ixyz + 2;
535 faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
546 for (
G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
548 xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
549 xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
550 xyz[ixyz][2] = GetCorner(iCorner).z;
556 nNodes = numSide * GetNumRZCorner();
557 nFaces = numSide * GetNumRZCorner();;
558 xyz =
new double3[nNodes];
559 faces_vec =
new int4[nFaces];
560 const G4double dPhi = (GetEndPhi() - GetStartPhi()) / numSide;
562 G4int ixyz = 0, iface = 0;
563 for (
G4int iSide = 0; iSide < numSide; ++iSide)
565 for (
G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
567 xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
568 xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
569 xyz[ixyz][2] = GetCorner(iCorner).z;
571 if (iSide < numSide - 1)
573 if (iCorner < GetNumRZCorner() - 1)
575 faces_vec[iface][0] = ixyz + 1;
576 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
577 faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
578 faces_vec[iface][3] = -(ixyz + 2);
582 faces_vec[iface][0] = ixyz + 1;
583 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
584 faces_vec[iface][2] = ixyz + 2;
585 faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
590 if (iCorner < GetNumRZCorner() - 1)
592 faces_vec[iface][0] = ixyz + 1;
593 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() - nFaces + 1);
594 faces_vec[iface][2] = ixyz + GetNumRZCorner() - nFaces + 2;
595 faces_vec[iface][3] = -(ixyz + 2);
599 faces_vec[iface][0] = ixyz + 1;
600 faces_vec[iface][1] = -(ixyz - nFaces + GetNumRZCorner() + 1);
601 faces_vec[iface][2] = ixyz - nFaces + 2;
602 faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
617 std::ostringstream message;
618 message <<
"Problem creating G4Polyhedron for: " << GetName();
619 G4Exception(
"G4GenericPolycone::CreatePolyhedron()",
"GeomSolids1002",
630 #endif // G4GEOM_USE_USOLIDS
void set(double x, double y, double z)
static const G4double kInfinity
G4int createPolyhedron(G4int Nnodes, G4int Nfaces, const G4double xyz[][3], const G4int faces[][4])
double B(double temperature)
static constexpr double twopi
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
const G4double kCarTolerance
std::vector< G4ThreeVector > G4ThreeVectorList
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
static const G4double emax
static G4int GetNumberOfRotationSteps()
CLHEP::Hep2Vector G4TwoVector
static constexpr double deg
G4double GetMaxExtent(const EAxis pAxis) const
G4double GetMinExtent(const EAxis pAxis) const