| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Implementation of G4Polyhedra, a CSG polyhedra,
27 // as an inherited class of G4VCSGfaceted.
28 //
29 // Utility classes:
30 // * G4EnclosingCylinder: decided a quick check of geometry would be a
31 // good idea (for CPU speed). If the quick check fails, the regular
32 // full-blown G4VCSGfaceted version is invoked.
33 // * G4ReduciblePolygon: Really meant as a check of input parameters,
34 // this utility class also "converts" the GEANT3-like PGON/PCON
35 // arguments into the newer ones.
36 // Both these classes are implemented outside this file because they are
37 // shared with G4Polycone.
38 //
39 // Author: David C. Williams (davidw@scipp.ucsc.edu)
40 // --------------------------------------------------------------------
41
42 #include "G4Polyhedra.hh"
43
44 #if !defined(G4GEOM_USE_UPOLYHEDRA)
45
46 #include "G4PolyhedraSide.hh"
47 #include "G4PolyPhiFace.hh"
48
49 #include "G4GeomTools.hh"
50 #include "G4VoxelLimits.hh"
51 #include "G4AffineTransform.hh"
52 #include "G4BoundingEnvelope.hh"
53
54 #include "G4QuickRand.hh"
55
56 #include "G4EnclosingCylinder.hh"
57 #include "G4ReduciblePolygon.hh"
58 #include "G4VPVParameterisation.hh"
59
60 namespace
61 {
62 G4Mutex surface_elementsMutex = G4MUTEX_INITIALIZER;
63 }
64
65 using namespace CLHEP;
66
67 // Constructor (GEANT3 style parameters)
68 //
69 // GEANT3 PGON radii are specified in the distance to the norm of each face.
70 //
71 G4Polyhedra::G4Polyhedra( const G4String& name,
72 G4double phiStart,
73 G4double thePhiTotal,
74 G4int theNumSide,
75 G4int numZPlanes,
76 const G4double zPlane[],
77 const G4double rInner[],
78 const G4double rOuter[] )
79 : G4VCSGfaceted( name )
80 {
81 if (theNumSide <= 0)
82 {
83 std::ostringstream message;
84 message << "Solid must have at least one side - " << GetName() << G4endl
85 << " No sides specified !";
86 G4Exception("G4Polyhedra::G4Polyhedra()", "GeomSolids0002",
87 FatalErrorInArgument, message);
88 }
89
90 //
91 // Calculate conversion factor from G3 radius to G4 radius
92 //
93 G4double phiTotal = thePhiTotal;
94 if ( (phiTotal <=0) || (phiTotal >= twopi*(1-DBL_EPSILON)) )
95 { phiTotal = twopi; }
96 G4double convertRad = std::cos(0.5*phiTotal/theNumSide);
97
98 //
99 // Some historical stuff
100 //
101 original_parameters = new G4PolyhedraHistorical;
102
103 original_parameters->numSide = theNumSide;
104 original_parameters->Start_angle = phiStart;
105 original_parameters->Opening_angle = phiTotal;
106 original_parameters->Num_z_planes = numZPlanes;
107 original_parameters->Z_values = new G4double[numZPlanes];
108 original_parameters->Rmin = new G4double[numZPlanes];
109 original_parameters->Rmax = new G4double[numZPlanes];
110
111 for (G4int i=0; i<numZPlanes; ++i)
112 {
113 if (( i < numZPlanes-1) && ( zPlane[i] == zPlane[i+1] ))
114 {
115 if( (rInner[i] > rOuter[i+1])
116 ||(rInner[i+1] > rOuter[i]) )
117 {
118 DumpInfo();
119 std::ostringstream message;
120 message << "Cannot create a Polyhedra with no contiguous segments."
121 << G4endl
122 << " Segments are not contiguous !" << G4endl
123 << " rMin[" << i << "] = " << rInner[i]
124 << " -- rMax[" << i+1 << "] = " << rOuter[i+1] << G4endl
125 << " rMin[" << i+1 << "] = " << rInner[i+1]
126 << " -- rMax[" << i << "] = " << rOuter[i];
127 G4Exception("G4Polyhedra::G4Polyhedra()", "GeomSolids0002",
128 FatalErrorInArgument, message);
129 }
130 }
131 original_parameters->Z_values[i] = zPlane[i];
132 original_parameters->Rmin[i] = rInner[i]/convertRad;
133 original_parameters->Rmax[i] = rOuter[i]/convertRad;
134 }
135
136
137 //
138 // Build RZ polygon using special PCON/PGON GEANT3 constructor
139 //
140 auto rz = new G4ReduciblePolygon( rInner, rOuter, zPlane, numZPlanes );
141 rz->ScaleA( 1/convertRad );
142
143 //
144 // Do the real work
145 //
146 Create( phiStart, phiTotal, theNumSide, rz );
147
148 delete rz;
149 }
150
151 // Constructor (generic parameters)
152 //
153 G4Polyhedra::G4Polyhedra( const G4String& name,
154 G4double phiStart,
155 G4double phiTotal,
156 G4int theNumSide,
157 G4int numRZ,
158 const G4double r[],
159 const G4double z[] )
160 : G4VCSGfaceted( name ), genericPgon(true)
161 {
162 if (theNumSide <= 0)
163 {
164 std::ostringstream message;
165 message << "Solid must have at least one side - " << GetName() << G4endl
166 << " No sides specified !";
167 G4Exception("G4Polyhedra::G4Polyhedra()", "GeomSolids0002",
168 FatalErrorInArgument, message);
169 }
170
171 auto rz = new G4ReduciblePolygon( r, z, numRZ );
172
173 Create( phiStart, phiTotal, theNumSide, rz );
174
175 // Set original_parameters struct for consistency
176 //
177 SetOriginalParameters(rz);
178
179 delete rz;
180 }
181
182 // Create
183 //
184 // Generic create routine, called by each constructor
185 // after conversion of arguments
186 //
187 void G4Polyhedra::Create( G4double phiStart,
188 G4double phiTotal,
189 G4int theNumSide,
190 G4ReduciblePolygon* rz )
191 {
192 //
193 // Perform checks of rz values
194 //
195 if (rz->Amin() < 0.0)
196 {
197 std::ostringstream message;
198 message << "Illegal input parameters - " << GetName() << G4endl
199 << " All R values must be >= 0 !";
200 G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
201 FatalErrorInArgument, message);
202 }
203
204 G4double rzArea = rz->Area();
205 if (rzArea < -kCarTolerance)
206 {
207 rz->ReverseOrder();
208 }
209 else if (rzArea < kCarTolerance)
210 {
211 std::ostringstream message;
212 message << "Illegal input parameters - " << GetName() << G4endl
213 << " R/Z cross section is zero or near zero: " << rzArea;
214 G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
215 FatalErrorInArgument, message);
216 }
217
218 if ( (!rz->RemoveDuplicateVertices( kCarTolerance ))
219 || (!rz->RemoveRedundantVertices( kCarTolerance )) )
220 {
221 std::ostringstream message;
222 message << "Illegal input parameters - " << GetName() << G4endl
223 << " Too few unique R/Z values !";
224 G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
225 FatalErrorInArgument, message);
226 }
227
228 if (rz->CrossesItself( 1/kInfinity ))
229 {
230 std::ostringstream message;
231 message << "Illegal input parameters - " << GetName() << G4endl
232 << " R/Z segments cross !";
233 G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
234 FatalErrorInArgument, message);
235 }
236
237 numCorner = rz->NumVertices();
238
239
240 startPhi = phiStart;
241 while( startPhi < 0 ) // Loop checking, 13.08.2015, G.Cosmo
242 startPhi += twopi;
243 //
244 // Phi opening? Account for some possible roundoff, and interpret
245 // nonsense value as representing no phi opening
246 //
247 if ( (phiTotal <= 0) || (phiTotal > twopi*(1-DBL_EPSILON)) )
248 {
249 phiIsOpen = false;
250 endPhi = startPhi + twopi;
251 }
252 else
253 {
254 phiIsOpen = true;
255 endPhi = startPhi + phiTotal;
256 }
257
258 //
259 // Save number sides
260 //
261 numSide = theNumSide;
262
263 //
264 // Allocate corner array.
265 //
266 corners = new G4PolyhedraSideRZ[numCorner];
267
268 //
269 // Copy corners
270 //
271 G4ReduciblePolygonIterator iterRZ(rz);
272
273 G4PolyhedraSideRZ *next = corners;
274 iterRZ.Begin();
275 do // Loop checking, 13.08.2015, G.Cosmo
276 {
277 next->r = iterRZ.GetA();
278 next->z = iterRZ.GetB();
279 } while( ++next, iterRZ.Next() );
280
281 //
282 // Allocate face pointer array
283 //
284 numFace = phiIsOpen ? numCorner+2 : numCorner;
285 faces = new G4VCSGface*[numFace];
286
287 //
288 // Construct side faces
289 //
290 // To do so properly, we need to keep track of four successive RZ
291 // corners.
292 //
293 // But! Don't construct a face if both points are at zero radius!
294 //
295 G4PolyhedraSideRZ* corner = corners,
296 * prev = corners + numCorner-1,
297 * nextNext;
298 G4VCSGface** face = faces;
299 do // Loop checking, 13.08.2015, G.Cosmo
300 {
301 next = corner+1;
302 if (next >= corners+numCorner) next = corners;
303 nextNext = next+1;
304 if (nextNext >= corners+numCorner) nextNext = corners;
305
306 if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue;
307 /*
308 // We must decide here if we can dare declare one of our faces
309 // as having a "valid" normal (i.e. allBehind = true). This
310 // is never possible if the face faces "inward" in r *unless*
311 // we have only one side
312 //
313 G4bool allBehind;
314 if ((corner->z > next->z) && (numSide > 1))
315 {
316 allBehind = false;
317 }
318 else
319 {
320 //
321 // Otherwise, it is only true if the line passing
322 // through the two points of the segment do not
323 // split the r/z cross section
324 //
325 allBehind = !rz->BisectedBy( corner->r, corner->z,
326 next->r, next->z, kCarTolerance );
327 }
328 */
329 *face++ = new G4PolyhedraSide( prev, corner, next, nextNext,
330 numSide, startPhi, endPhi-startPhi, phiIsOpen );
331 } while( prev=corner, corner=next, corner > corners );
332
333 if (phiIsOpen)
334 {
335 //
336 // Construct phi open edges
337 //
338 *face++ = new G4PolyPhiFace( rz, startPhi, phiTotal/numSide, endPhi );
339 *face++ = new G4PolyPhiFace( rz, endPhi, phiTotal/numSide, startPhi );
340 }
341
342 //
343 // We might have dropped a face or two: recalculate numFace
344 //
345 numFace = (G4int)(face-faces);
346
347 //
348 // Make enclosingCylinder
349 //
350 enclosingCylinder =
351 new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal );
352 }
353
354 // Fake default constructor - sets only member data and allocates memory
355 // for usage restricted to object persistency.
356 //
357 G4Polyhedra::G4Polyhedra( __void__& a )
358 : G4VCSGfaceted(a), startPhi(0.), endPhi(0.)
359 {
360 }
361
362 // Destructor
363 //
364 G4Polyhedra::~G4Polyhedra()
365 {
366 delete [] corners;
367 delete original_parameters;
368 delete enclosingCylinder;
369 delete fElements;
370 delete fpPolyhedron;
371 corners = nullptr;
372 original_parameters = nullptr;
373 enclosingCylinder = nullptr;
374 fElements = nullptr;
375 fpPolyhedron = nullptr;
376 }
377
378 // Copy constructor
379 //
380 G4Polyhedra::G4Polyhedra( const G4Polyhedra& source )
381 : G4VCSGfaceted( source )
382 {
383 CopyStuff( source );
384 }
385
386 // Assignment operator
387 //
388 G4Polyhedra &G4Polyhedra::operator=( const G4Polyhedra& source )
389 {
390 if (this == &source) return *this;
391
392 G4VCSGfaceted::operator=( source );
393
394 delete [] corners;
395 delete original_parameters;
396 delete enclosingCylinder;
397
398 CopyStuff( source );
399
400 return *this;
401 }
402
403 // CopyStuff
404 //
405 void G4Polyhedra::CopyStuff( const G4Polyhedra& source )
406 {
407 //
408 // Simple stuff
409 //
410 numSide = source.numSide;
411 startPhi = source.startPhi;
412 endPhi = source.endPhi;
413 phiIsOpen = source.phiIsOpen;
414 numCorner = source.numCorner;
415 genericPgon= source.genericPgon;
416
417 //
418 // The corner array
419 //
420 corners = new G4PolyhedraSideRZ[numCorner];
421
422 G4PolyhedraSideRZ* corn = corners,
423 * sourceCorn = source.corners;
424 do // Loop checking, 13.08.2015, G.Cosmo
425 {
426 *corn = *sourceCorn;
427 } while( ++sourceCorn, ++corn < corners+numCorner );
428
429 //
430 // Original parameters
431 //
432 if (source.original_parameters != nullptr)
433 {
434 original_parameters =
435 new G4PolyhedraHistorical( *source.original_parameters );
436 }
437
438 //
439 // Enclosing cylinder
440 //
441 enclosingCylinder = new G4EnclosingCylinder( *source.enclosingCylinder );
442
443 //
444 // Surface elements
445 //
446 delete fElements;
447 fElements = nullptr;
448
449 //
450 // Polyhedron
451 //
452 fRebuildPolyhedron = false;
453 delete fpPolyhedron;
454 fpPolyhedron = nullptr;
455 }
456
457 // Reset
458 //
459 // Recalculates and reshapes the solid, given pre-assigned scaled
460 // original_parameters.
461 //
462 G4bool G4Polyhedra::Reset()
463 {
464 if (genericPgon)
465 {
466 std::ostringstream message;
467 message << "Solid " << GetName() << " built using generic construct."
468 << G4endl << "Not applicable to the generic construct !";
469 G4Exception("G4Polyhedra::Reset()", "GeomSolids1001",
470 JustWarning, message, "Parameters NOT resetted.");
471 return true;
472 }
473
474 //
475 // Clear old setup
476 //
477 G4VCSGfaceted::DeleteStuff();
478 delete [] corners;
479 delete enclosingCylinder;
480 delete fElements;
481 corners = nullptr;
482 fElements = nullptr;
483 enclosingCylinder = nullptr;
484
485 //
486 // Rebuild polyhedra
487 //
488 auto rz = new G4ReduciblePolygon( original_parameters->Rmin,
489 original_parameters->Rmax,
490 original_parameters->Z_values,
491 original_parameters->Num_z_planes );
492 Create( original_parameters->Start_angle,
493 original_parameters->Opening_angle,
494 original_parameters->numSide, rz );
495 delete rz;
496
497 return false;
498 }
499
500 // Inside
501 //
502 // This is an override of G4VCSGfaceted::Inside, created in order
503 // to speed things up by first checking with G4EnclosingCylinder.
504 //
505 EInside G4Polyhedra::Inside( const G4ThreeVector& p ) const
506 {
507 //
508 // Quick test
509 //
510 if (enclosingCylinder->MustBeOutside(p)) return kOutside;
511
512 //
513 // Long answer
514 //
515 return G4VCSGfaceted::Inside(p);
516 }
517
518 // DistanceToIn
519 //
520 // This is an override of G4VCSGfaceted::Inside, created in order
521 // to speed things up by first checking with G4EnclosingCylinder.
522 //
523 G4double G4Polyhedra::DistanceToIn( const G4ThreeVector& p,
524 const G4ThreeVector& v ) const
525 {
526 //
527 // Quick test
528 //
529 if (enclosingCylinder->ShouldMiss(p,v))
530 return kInfinity;
531
532 //
533 // Long answer
534 //
535 return G4VCSGfaceted::DistanceToIn( p, v );
536 }
537
538 // DistanceToIn
539 //
540 G4double G4Polyhedra::DistanceToIn( const G4ThreeVector& p ) const
541 {
542 return G4VCSGfaceted::DistanceToIn(p);
543 }
544
545 // Get bounding box
546 //
547 void G4Polyhedra::BoundingLimits(G4ThreeVector& pMin,
548 G4ThreeVector& pMax) const
549 {
550 G4double rmin = kInfinity, rmax = -kInfinity;
551 G4double zmin = kInfinity, zmax = -kInfinity;
552 for (G4int i=0; i<GetNumRZCorner(); ++i)
553 {
554 G4PolyhedraSideRZ corner = GetCorner(i);
555 if (corner.r < rmin) rmin = corner.r;
556 if (corner.r > rmax) rmax = corner.r;
557 if (corner.z < zmin) zmin = corner.z;
558 if (corner.z > zmax) zmax = corner.z;
559 }
560
561 G4double sphi = GetStartPhi();
562 G4double ephi = GetEndPhi();
563 G4double dphi = IsOpen() ? ephi-sphi : twopi;
564 G4int ksteps = GetNumSide();
565 G4double astep = dphi/ksteps;
566 G4double sinStep = std::sin(astep);
567 G4double cosStep = std::cos(astep);
568
569 G4double sinCur = GetSinStartPhi();
570 G4double cosCur = GetCosStartPhi();
571 if (!IsOpen()) rmin = 0.;
572 G4double xmin = rmin*cosCur, xmax = xmin;
573 G4double ymin = rmin*sinCur, ymax = ymin;
574 for (G4int k=0; k<ksteps+1; ++k)
575 {
576 G4double x = rmax*cosCur;
577 if (x < xmin) xmin = x;
578 if (x > xmax) xmax = x;
579 G4double y = rmax*sinCur;
580 if (y < ymin) ymin = y;
581 if (y > ymax) ymax = y;
582 if (rmin > 0)
583 {
584 G4double xx = rmin*cosCur;
585 if (xx < xmin) xmin = xx;
586 if (xx > xmax) xmax = xx;
587 G4double yy = rmin*sinCur;
588 if (yy < ymin) ymin = yy;
589 if (yy > ymax) ymax = yy;
590 }
591 G4double sinTmp = sinCur;
592 sinCur = sinCur*cosStep + cosCur*sinStep;
593 cosCur = cosCur*cosStep - sinTmp*sinStep;
594 }
595 pMin.set(xmin,ymin,zmin);
596 pMax.set(xmax,ymax,zmax);
597
598 // Check correctness of the bounding box
599 //
600 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
601 {
602 std::ostringstream message;
603 message << "Bad bounding box (min >= max) for solid: "
604 << GetName() << " !"
605 << "\npMin = " << pMin
606 << "\npMax = " << pMax;
607 G4Exception("G4Polyhedra::BoundingLimits()", "GeomMgt0001",
608 JustWarning, message);
609 DumpInfo();
610 }
611 }
612
613 // Calculate extent under transform and specified limit
614 //
615 G4bool G4Polyhedra::CalculateExtent(const EAxis pAxis,
616 const G4VoxelLimits& pVoxelLimit,
617 const G4AffineTransform& pTransform,
618 G4double& pMin, G4double& pMax) const
619 {
620 G4ThreeVector bmin, bmax;
621 G4bool exist;
622
623 // Check bounding box (bbox)
624 //
625 BoundingLimits(bmin,bmax);
626 G4BoundingEnvelope bbox(bmin,bmax);
627 #ifdef G4BBOX_EXTENT
628 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
629 #endif
630 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
631 {
632 return exist = pMin < pMax;
633 }
634
635 // To find the extent, RZ contour of the polycone is subdivided
636 // in triangles. The extent is calculated as cumulative extent of
637 // all sub-polycones formed by rotation of triangles around Z
638 //
639 G4TwoVectorList contourRZ;
640 G4TwoVectorList triangles;
641 std::vector<G4int> iout;
642 G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
643 G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
644
645 // get RZ contour, ensure anticlockwise order of corners
646 for (G4int i=0; i<GetNumRZCorner(); ++i)
647 {
648 G4PolyhedraSideRZ corner = GetCorner(i);
649 contourRZ.emplace_back(corner.r,corner.z);
650 }
651 G4GeomTools::RemoveRedundantVertices(contourRZ,iout,2*kCarTolerance);
652 G4double area = G4GeomTools::PolygonArea(contourRZ);
653 if (area < 0.) std::reverse(contourRZ.begin(),contourRZ.end());
654
655 // triangulate RZ countour
656 if (!G4GeomTools::TriangulatePolygon(contourRZ,triangles))
657 {
658 std::ostringstream message;
659 message << "Triangulation of RZ contour has failed for solid: "
660 << GetName() << " !"
661 << "\nExtent has been calculated using boundary box";
662 G4Exception("G4Polyhedra::CalculateExtent()",
663 "GeomMgt1002",JustWarning,message);
664 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
665 }
666
667 // set trigonometric values
668 G4double sphi = GetStartPhi();
669 G4double ephi = GetEndPhi();
670 G4double dphi = IsOpen() ? ephi-sphi : twopi;
671 G4int ksteps = GetNumSide();
672 G4double astep = dphi/ksteps;
673 G4double sinStep = std::sin(astep);
674 G4double cosStep = std::cos(astep);
675 G4double sinStart = GetSinStartPhi();
676 G4double cosStart = GetCosStartPhi();
677
678 // allocate vector lists
679 std::vector<const G4ThreeVectorList *> polygons;
680 polygons.resize(ksteps+1);
681 for (G4int k=0; k<ksteps+1; ++k)
682 {
683 polygons[k] = new G4ThreeVectorList(3);
684 }
685
686 // main loop along triangles
687 pMin = kInfinity;
688 pMax = -kInfinity;
689 G4int ntria = (G4int)triangles.size()/3;
690 for (G4int i=0; i<ntria; ++i)
691 {
692 G4double sinCur = sinStart;
693 G4double cosCur = cosStart;
694 G4int i3 = i*3;
695 for (G4int k=0; k<ksteps+1; ++k) // rotate triangle
696 {
697 auto ptr = const_cast<G4ThreeVectorList*>(polygons[k]);
698 auto iter = ptr->begin();
699 iter->set(triangles[i3+0].x()*cosCur,
700 triangles[i3+0].x()*sinCur,
701 triangles[i3+0].y());
702 iter++;
703 iter->set(triangles[i3+1].x()*cosCur,
704 triangles[i3+1].x()*sinCur,
705 triangles[i3+1].y());
706 iter++;
707 iter->set(triangles[i3+2].x()*cosCur,
708 triangles[i3+2].x()*sinCur,
709 triangles[i3+2].y());
710
711 G4double sinTmp = sinCur;
712 sinCur = sinCur*cosStep + cosCur*sinStep;
713 cosCur = cosCur*cosStep - sinTmp*sinStep;
714 }
715
716 // set sub-envelope and adjust extent
717 G4double emin,emax;
718 G4BoundingEnvelope benv(polygons);
719 if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
720 if (emin < pMin) pMin = emin;
721 if (emax > pMax) pMax = emax;
722 if (eminlim > pMin && emaxlim < pMax) break; // max possible extent
723 }
724 // free memory
725 for (G4int k=0; k<ksteps+1; ++k) { delete polygons[k]; polygons[k]=nullptr;}
726 return (pMin < pMax);
727 }
728
729 // ComputeDimensions
730 //
731 void G4Polyhedra::ComputeDimensions( G4VPVParameterisation* p,
732 const G4int n,
733 const G4VPhysicalVolume* pRep )
734 {
735 p->ComputeDimensions(*this,n,pRep);
736 }
737
738 // GetEntityType
739 //
740 G4GeometryType G4Polyhedra::GetEntityType() const
741 {
742 return {"G4Polyhedra"};
743 }
744
745 // IsFaceted
746 //
747 G4bool G4Polyhedra::IsFaceted() const
748 {
749 return true;
750 }
751
752 // Make a clone of the object
753 //
754 G4VSolid* G4Polyhedra::Clone() const
755 {
756 return new G4Polyhedra(*this);
757 }
758
759 // Stream object contents to an output stream
760 //
761 std::ostream& G4Polyhedra::StreamInfo( std::ostream& os ) const
762 {
763 G4long oldprc = os.precision(16);
764 os << "-----------------------------------------------------------\n"
765 << " *** Dump for solid - " << GetName() << " ***\n"
766 << " ===================================================\n"
767 << " Solid type: G4Polyhedra\n"
768 << " Parameters: \n"
769 << " starting phi angle : " << startPhi/degree << " degrees \n"
770 << " ending phi angle : " << endPhi/degree << " degrees \n"
771 << " number of sides : " << numSide << " \n";
772 G4int i=0;
773 if (!genericPgon)
774 {
775 G4int numPlanes = original_parameters->Num_z_planes;
776 os << " number of Z planes: " << numPlanes << "\n"
777 << " Z values: \n";
778 for (i=0; i<numPlanes; ++i)
779 {
780 os << " Z plane " << i << ": "
781 << original_parameters->Z_values[i] << "\n";
782 }
783 os << " Tangent distances to inner surface (Rmin): \n";
784 for (i=0; i<numPlanes; ++i)
785 {
786 os << " Z plane " << i << ": "
787 << original_parameters->Rmin[i] << "\n";
788 }
789 os << " Tangent distances to outer surface (Rmax): \n";
790 for (i=0; i<numPlanes; ++i)
791 {
792 os << " Z plane " << i << ": "
793 << original_parameters->Rmax[i] << "\n";
794 }
795 }
796 os << " number of RZ points: " << numCorner << "\n"
797 << " RZ values (corners): \n";
798 for (i=0; i<numCorner; ++i)
799 {
800 os << " "
801 << corners[i].r << ", " << corners[i].z << "\n";
802 }
803 os << "-----------------------------------------------------------\n";
804 os.precision(oldprc);
805
806 return os;
807 }
808
809 //////////////////////////////////////////////////////////////////////////
810 //
811 // Return volume
812
813 G4double G4Polyhedra::GetCubicVolume()
814 {
815 if (fCubicVolume == 0.)
816 {
817 G4double total = 0.;
818 G4int nrz = GetNumRZCorner();
819 G4PolyhedraSideRZ a = GetCorner(nrz - 1);
820 for (G4int i=0; i<nrz; ++i)
821 {
822 G4PolyhedraSideRZ b = GetCorner(i);
823 total += (b.r*b.r + b.r*a.r + a.r*a.r)*(b.z - a.z);
824 a = b;
825 }
826 fCubicVolume = std::abs(total)*
827 std::sin((GetEndPhi() - GetStartPhi())/GetNumSide())*GetNumSide()/6.;
828 }
829 return fCubicVolume;
830 }
831
832 //////////////////////////////////////////////////////////////////////////
833 //
834 // Return surface area
835
836 G4double G4Polyhedra::GetSurfaceArea()
837 {
838 if (fSurfaceArea == 0.)
839 {
840 G4double total = 0.;
841 G4int nrz = GetNumRZCorner();
842 if (IsOpen())
843 {
844 G4PolyhedraSideRZ a = GetCorner(nrz - 1);
845 for (G4int i=0; i<nrz; ++i)
846 {
847 G4PolyhedraSideRZ b = GetCorner(i);
848 total += a.r*b.z - a.z*b.r;
849 a = b;
850 }
851 total = std::abs(total);
852 }
853 G4double alp = (GetEndPhi() - GetStartPhi())/GetNumSide();
854 G4double cosa = std::cos(alp);
855 G4double sina = std::sin(alp);
856 G4PolyhedraSideRZ a = GetCorner(nrz - 1);
857 for (G4int i=0; i<nrz; ++i)
858 {
859 G4PolyhedraSideRZ b = GetCorner(i);
860 G4ThreeVector p1(a.r, 0, a.z);
861 G4ThreeVector p2(a.r*cosa, a.r*sina, a.z);
862 G4ThreeVector p3(b.r*cosa, b.r*sina, b.z);
863 G4ThreeVector p4(b.r, 0, b.z);
864 total += GetNumSide()*(G4GeomTools::QuadAreaNormal(p1, p2, p3, p4)).mag();
865 a = b;
866 }
867 fSurfaceArea = total;
868 }
869 return fSurfaceArea;
870 }
871
872 //////////////////////////////////////////////////////////////////////////
873 //
874 // Set vector of surface elements, auxiliary method for sampling
875 // random points on surface
876
877 void G4Polyhedra::SetSurfaceElements() const
878 {
879 fElements = new std::vector<G4Polyhedra::surface_element>;
880 G4double total = 0.;
881 G4int nrz = GetNumRZCorner();
882
883 // set lateral surface elements
884 G4double dphi = (GetEndPhi() - GetStartPhi())/GetNumSide();
885 G4double cosa = std::cos(dphi);
886 G4double sina = std::sin(dphi);
887 G4int ia = nrz - 1;
888 for (G4int ib=0; ib<nrz; ++ib)
889 {
890 G4PolyhedraSideRZ a = GetCorner(ia);
891 G4PolyhedraSideRZ b = GetCorner(ib);
892 G4Polyhedra::surface_element selem;
893 selem.i0 = ia;
894 selem.i1 = ib;
895 ia = ib;
896 if (a.r == 0. && b.r == 0.) continue;
897 G4ThreeVector p1(a.r, 0, a.z);
898 G4ThreeVector p2(a.r*cosa, a.r*sina, a.z);
899 G4ThreeVector p3(b.r*cosa, b.r*sina, b.z);
900 G4ThreeVector p4(b.r, 0, b.z);
901 if (a.r > 0.)
902 {
903 selem.i2 = -1;
904 total += GetNumSide()*(G4GeomTools::TriangleAreaNormal(p1, p2, p3)).mag();
905 selem.area = total;
906 fElements->push_back(selem);
907 }
908 if (b.r > 0.)
909 {
910 selem.i2 = -2;
911 total += GetNumSide()*(G4GeomTools::TriangleAreaNormal(p1, p3, p4)).mag();
912 selem.area = total;
913 fElements->push_back(selem);
914 }
915 }
916
917 // set elements for phi cuts
918 if (IsOpen())
919 {
920 G4TwoVectorList contourRZ;
921 std::vector<G4int> triangles;
922 for (G4int i=0; i<nrz; ++i)
923 {
924 G4PolyhedraSideRZ corner = GetCorner(i);
925 contourRZ.emplace_back(corner.r, corner.z);
926 }
927 G4GeomTools::TriangulatePolygon(contourRZ, triangles);
928 auto ntria = (G4int)triangles.size();
929 for (G4int i=0; i<ntria; i+=3)
930 {
931 G4Polyhedra::surface_element selem;
932 selem.i0 = triangles[i];
933 selem.i1 = triangles[i+1];
934 selem.i2 = triangles[i+2];
935 G4PolyhedraSideRZ a = GetCorner(selem.i0);
936 G4PolyhedraSideRZ b = GetCorner(selem.i1);
937 G4PolyhedraSideRZ c = GetCorner(selem.i2);
938 G4double stria =
939 std::abs(G4GeomTools::TriangleArea(a.r, a.z, b.r, b.z, c.r, c.z));
940 total += stria;
941 selem.area = total;
942 fElements->push_back(selem); // start phi
943 total += stria;
944 selem.area = total;
945 selem.i0 += nrz;
946 fElements->push_back(selem); // end phi
947 }
948 }
949 }
950
951 //////////////////////////////////////////////////////////////////////////
952 //
953 // Generate random point on surface
954
955 G4ThreeVector G4Polyhedra::GetPointOnSurface() const
956 {
957 // Set surface elements
958 if (fElements == nullptr)
959 {
960 G4AutoLock l(&surface_elementsMutex);
961 SetSurfaceElements();
962 l.unlock();
963 }
964
965 // Select surface element
966 G4Polyhedra::surface_element selem;
967 selem = fElements->back();
968 G4double select = selem.area*G4QuickRand();
969 auto it = std::lower_bound(fElements->begin(), fElements->end(), select,
970 [](const G4Polyhedra::surface_element& x, G4double val)
971 -> G4bool { return x.area < val; });
972
973 // Generate random point
974 G4double x = 0, y = 0, z = 0;
975 G4double u = G4QuickRand();
976 G4double v = G4QuickRand();
977 if (u + v > 1.) { u = 1. - u; v = 1. - v; }
978 G4int i0 = (*it).i0;
979 G4int i1 = (*it).i1;
980 G4int i2 = (*it).i2;
981 if (i2 < 0) // lateral surface
982 {
983 // sample point
984 G4int nside = GetNumSide();
985 G4double dphi = (GetEndPhi() - GetStartPhi())/nside;
986 G4double cosa = std::cos(dphi);
987 G4double sina = std::sin(dphi);
988 G4PolyhedraSideRZ a = GetCorner(i0);
989 G4PolyhedraSideRZ b = GetCorner(i1);
990 G4ThreeVector p0(a.r, 0, a.z);
991 G4ThreeVector p1(b.r, 0, b.z);
992 G4ThreeVector p2(b.r*cosa, b.r*sina, b.z);
993 if (i2 == -1) p1.set(a.r*cosa, a.r*sina, a.z);
994 p0 += (p1 - p0)*u + (p2 - p0)*v;
995 // find selected side and rotate point
996 G4double scurr = (*it).area;
997 G4double sprev = (it == fElements->begin()) ? 0. : (*(--it)).area;
998 G4int iside = nside*(select - sprev)/(scurr - sprev);
999 if (iside == 0 && GetStartPhi() == 0.)
1000 {
1001 x = p0.x();
1002 y = p0.y();
1003 z = p0.z();
1004 }
1005 else
1006 {
1007 if (iside == nside) --iside; // iside must be less then nside
1008 G4double phi = iside*dphi + GetStartPhi();
1009 G4double cosphi = std::cos(phi);
1010 G4double sinphi = std::sin(phi);
1011 x = p0.x()*cosphi - p0.y()*sinphi;
1012 y = p0.x()*sinphi + p0.y()*cosphi;
1013 z = p0.z();
1014 }
1015 }
1016 else // phi cut
1017 {
1018 G4int nrz = GetNumRZCorner();
1019 G4double phi = (i0 < nrz) ? GetStartPhi() : GetEndPhi();
1020 if (i0 >= nrz) { i0 -= nrz; }
1021 G4PolyhedraSideRZ p0 = GetCorner(i0);
1022 G4PolyhedraSideRZ p1 = GetCorner(i1);
1023 G4PolyhedraSideRZ p2 = GetCorner(i2);
1024 G4double r = (p1.r - p0.r)*u + (p2.r - p0.r)*v + p0.r;
1025 x = r*std::cos(phi);
1026 y = r*std::sin(phi);
1027 z = (p1.z - p0.z)*u + (p2.z - p0.z)*v + p0.z;
1028 }
1029 return {x, y, z};
1030 }
1031
1032 //////////////////////////////////////////////////////////////////////////
1033 //
1034 // CreatePolyhedron
1035
1036 G4Polyhedron* G4Polyhedra::CreatePolyhedron() const
1037 {
1038 std::vector<G4TwoVector> rz(numCorner);
1039 for (G4int i = 0; i < numCorner; ++i)
1040 rz[i].set(corners[i].r, corners[i].z);
1041 return new G4PolyhedronPgon(startPhi, endPhi - startPhi, numSide, rz);
1042 }
1043
1044 // SetOriginalParameters
1045 //
1046 void G4Polyhedra::SetOriginalParameters(G4ReduciblePolygon* rz)
1047 {
1048 G4int numPlanes = numCorner;
1049 G4bool isConvertible = true;
1050 G4double Zmax=rz->Bmax();
1051 rz->StartWithZMin();
1052
1053 // Prepare vectors for storage
1054 //
1055 std::vector<G4double> Z;
1056 std::vector<G4double> Rmin;
1057 std::vector<G4double> Rmax;
1058
1059 G4int countPlanes=1;
1060 G4int icurr=0;
1061 G4int icurl=0;
1062
1063 // first plane Z=Z[0]
1064 //
1065 Z.push_back(corners[0].z);
1066 G4double Zprev=Z[0];
1067 if (Zprev == corners[1].z)
1068 {
1069 Rmin.push_back(corners[0].r);
1070 Rmax.push_back (corners[1].r);icurr=1;
1071 }
1072 else if (Zprev == corners[numPlanes-1].z)
1073 {
1074 Rmin.push_back(corners[numPlanes-1].r);
1075 Rmax.push_back (corners[0].r);
1076 icurl=numPlanes-1;
1077 }
1078 else
1079 {
1080 Rmin.push_back(corners[0].r);
1081 Rmax.push_back (corners[0].r);
1082 }
1083
1084 // next planes until last
1085 //
1086 G4int inextr=0, inextl=0;
1087 for (G4int i=0; i < numPlanes-2; ++i)
1088 {
1089 inextr=1+icurr;
1090 inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1091
1092 if((corners[inextr].z >= Zmax) & (corners[inextl].z >= Zmax)) { break; }
1093
1094 G4double Zleft = corners[inextl].z;
1095 G4double Zright = corners[inextr].z;
1096 if(Zright>Zleft)
1097 {
1098 Z.push_back(Zleft);
1099 countPlanes++;
1100 G4double difZr=corners[inextr].z - corners[icurr].z;
1101 G4double difZl=corners[inextl].z - corners[icurl].z;
1102
1103 if(std::fabs(difZl) < kCarTolerance)
1104 {
1105 if(std::fabs(difZr) < kCarTolerance)
1106 {
1107 Rmin.push_back(corners[inextl].r);
1108 Rmax.push_back(corners[icurr].r);
1109 }
1110 else
1111 {
1112 Rmin.push_back(corners[inextl].r);
1113 Rmax.push_back(corners[icurr].r + (Zleft-corners[icurr].z)/difZr
1114 *(corners[inextr].r - corners[icurr].r));
1115 }
1116 }
1117 else if (difZl >= kCarTolerance)
1118 {
1119 if(std::fabs(difZr) < kCarTolerance)
1120 {
1121 Rmin.push_back(corners[icurl].r);
1122 Rmax.push_back(corners[icurr].r);
1123 }
1124 else
1125 {
1126 Rmin.push_back(corners[icurl].r);
1127 Rmax.push_back(corners[icurr].r + (Zleft-corners[icurr].z)/difZr
1128 *(corners[inextr].r - corners[icurr].r));
1129 }
1130 }
1131 else
1132 {
1133 isConvertible=false; break;
1134 }
1135 icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1136 }
1137 else if(std::fabs(Zright-Zleft)<kCarTolerance) // Zright=Zleft
1138 {
1139 Z.push_back(Zleft);
1140 ++countPlanes;
1141 ++icurr;
1142
1143 icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1144
1145 Rmin.push_back(corners[inextl].r);
1146 Rmax.push_back (corners[inextr].r);
1147 }
1148 else // Zright<Zleft
1149 {
1150 Z.push_back(Zright);
1151 ++countPlanes;
1152
1153 G4double difZr=corners[inextr].z - corners[icurr].z;
1154 G4double difZl=corners[inextl].z - corners[icurl].z;
1155 if(std::fabs(difZr) < kCarTolerance)
1156 {
1157 if(std::fabs(difZl) < kCarTolerance)
1158 {
1159 Rmax.push_back(corners[inextr].r);
1160 Rmin.push_back(corners[icurr].r);
1161 }
1162 else
1163 {
1164 Rmin.push_back(corners[icurl].r + (Zright-corners[icurl].z)/difZl
1165 * (corners[inextl].r - corners[icurl].r));
1166 Rmax.push_back(corners[inextr].r);
1167 }
1168 ++icurr;
1169 } // plate
1170 else if (difZr >= kCarTolerance)
1171 {
1172 if(std::fabs(difZl) < kCarTolerance)
1173 {
1174 Rmax.push_back(corners[inextr].r);
1175 Rmin.push_back (corners[icurr].r);
1176 }
1177 else
1178 {
1179 Rmax.push_back(corners[inextr].r);
1180 Rmin.push_back (corners[icurl].r+(Zright-corners[icurl].z)/difZl
1181 * (corners[inextl].r - corners[icurl].r));
1182 }
1183 ++icurr;
1184 }
1185 else
1186 {
1187 isConvertible=false; break;
1188 }
1189 }
1190 } // end for loop
1191
1192 // last plane Z=Zmax
1193 //
1194 Z.push_back(Zmax);
1195 ++countPlanes;
1196 inextr=1+icurr;
1197 inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1198
1199 Rmax.push_back(corners[inextr].r);
1200 Rmin.push_back(corners[inextl].r);
1201
1202 // Set original parameters Rmin,Rmax,Z
1203 //
1204 if(isConvertible)
1205 {
1206 original_parameters = new G4PolyhedraHistorical;
1207 original_parameters->numSide = numSide;
1208 original_parameters->Z_values = new G4double[countPlanes];
1209 original_parameters->Rmin = new G4double[countPlanes];
1210 original_parameters->Rmax = new G4double[countPlanes];
1211
1212 for(G4int j=0; j < countPlanes; ++j)
1213 {
1214 original_parameters->Z_values[j] = Z[j];
1215 original_parameters->Rmax[j] = Rmax[j];
1216 original_parameters->Rmin[j] = Rmin[j];
1217 }
1218 original_parameters->Start_angle = startPhi;
1219 original_parameters->Opening_angle = endPhi-startPhi;
1220 original_parameters->Num_z_planes = countPlanes;
1221
1222 }
1223 else // Set parameters(r,z) with Rmin==0 as convention
1224 {
1225 #ifdef G4SPECSDEBUG
1226 std::ostringstream message;
1227 message << "Polyhedra " << GetName() << G4endl
1228 << "cannot be converted to Polyhedra with (Rmin,Rmaz,Z) parameters!";
1229 G4Exception("G4Polyhedra::SetOriginalParameters()",
1230 "GeomSolids0002", JustWarning, message);
1231 #endif
1232 original_parameters = new G4PolyhedraHistorical;
1233 original_parameters->numSide = numSide;
1234 original_parameters->Z_values = new G4double[numPlanes];
1235 original_parameters->Rmin = new G4double[numPlanes];
1236 original_parameters->Rmax = new G4double[numPlanes];
1237
1238 for(G4int j=0; j < numPlanes; ++j)
1239 {
1240 original_parameters->Z_values[j] = corners[j].z;
1241 original_parameters->Rmax[j] = corners[j].r;
1242 original_parameters->Rmin[j] = 0.0;
1243 }
1244 original_parameters->Start_angle = startPhi;
1245 original_parameters->Opening_angle = endPhi-startPhi;
1246 original_parameters->Num_z_planes = numPlanes;
1247 }
1248 }
1249
1250 #endif
1251