| 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 // G4GenericPolycone implementation
27 //
28 // Authors: T.Nikitina, G.Cosmo - CERN
29 // --------------------------------------------------------------------
30
31 #include "G4GenericPolycone.hh"
32
33 #if !defined(G4GEOM_USE_UGENERICPOLYCONE)
34
35 #include "G4PolyconeSide.hh"
36 #include "G4PolyPhiFace.hh"
37
38 #include "G4GeomTools.hh"
39 #include "G4VoxelLimits.hh"
40 #include "G4AffineTransform.hh"
41 #include "G4BoundingEnvelope.hh"
42
43 #include "G4QuickRand.hh"
44
45 #include "G4Polyhedron.hh"
46 #include "G4EnclosingCylinder.hh"
47 #include "G4ReduciblePolygon.hh"
48 #include "G4VPVParameterisation.hh"
49
50 namespace
51 {
52 G4Mutex surface_elementsMutex = G4MUTEX_INITIALIZER;
53 }
54
55 using namespace CLHEP;
56
57 // Constructor (generic parameters)
58 //
59 G4GenericPolycone::G4GenericPolycone( const G4String& name,
60 G4double phiStart,
61 G4double phiTotal,
62 G4int numRZ,
63 const G4double r[],
64 const G4double z[] )
65 : G4VCSGfaceted( name )
66 {
67
68 auto rz = new G4ReduciblePolygon( r, z, numRZ );
69
70 Create( phiStart, phiTotal, rz );
71
72 // Set original_parameters struct for consistency
73 //
74 //SetOriginalParameters(rz);
75
76 delete rz;
77 }
78
79 // Create
80 //
81 // Generic create routine, called by each constructor after
82 // conversion of arguments
83 //
84 void G4GenericPolycone::Create( G4double phiStart,
85 G4double phiTotal,
86 G4ReduciblePolygon *rz )
87 {
88 //
89 // Perform checks of rz values
90 //
91 if (rz->Amin() < 0.0)
92 {
93 std::ostringstream message;
94 message << "Illegal input parameters - " << GetName() << G4endl
95 << " All R values must be >= 0 !";
96 G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
97 FatalErrorInArgument, message);
98 }
99
100 G4double rzArea = rz->Area();
101 if (rzArea < -kCarTolerance)
102 {
103 rz->ReverseOrder();
104 }
105 else if (rzArea < kCarTolerance)
106 {
107 std::ostringstream message;
108 message << "Illegal input parameters - " << GetName() << G4endl
109 << " R/Z cross section is zero or near zero: " << rzArea;
110 G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
111 FatalErrorInArgument, message);
112 }
113
114 if ( (!rz->RemoveDuplicateVertices( kCarTolerance ))
115 || (!rz->RemoveRedundantVertices( kCarTolerance )) )
116 {
117 std::ostringstream message;
118 message << "Illegal input parameters - " << GetName() << G4endl
119 << " Too few unique R/Z values !";
120 G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
121 FatalErrorInArgument, message);
122 }
123
124 if (rz->CrossesItself(1/kInfinity))
125 {
126 std::ostringstream message;
127 message << "Illegal input parameters - " << GetName() << G4endl
128 << " R/Z segments cross !";
129 G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
130 FatalErrorInArgument, message);
131 }
132
133 numCorner = rz->NumVertices();
134
135 //
136 // Phi opening? Account for some possible roundoff, and interpret
137 // nonsense value as representing no phi opening
138 //
139 if (phiTotal <= 0 || phiTotal > twopi-1E-10)
140 {
141 phiIsOpen = false;
142 startPhi = 0;
143 endPhi = twopi;
144 }
145 else
146 {
147 phiIsOpen = true;
148
149 //
150 // Convert phi into our convention
151 //
152 startPhi = phiStart;
153 while( startPhi < 0 ) // Loop checking, 13.08.2015, G.Cosmo
154 startPhi += twopi;
155
156 endPhi = phiStart+phiTotal;
157 while( endPhi < startPhi ) // Loop checking, 13.08.2015, G.Cosmo
158 endPhi += twopi;
159 }
160
161 //
162 // Allocate corner array.
163 //
164 corners = new G4PolyconeSideRZ[numCorner];
165
166 //
167 // Copy corners
168 //
169 G4ReduciblePolygonIterator iterRZ(rz);
170
171 G4PolyconeSideRZ* next = corners;
172 iterRZ.Begin();
173 do // Loop checking, 13.08.2015, G.Cosmo
174 {
175 next->r = iterRZ.GetA();
176 next->z = iterRZ.GetB();
177 } while( ++next, iterRZ.Next() );
178
179 //
180 // Allocate face pointer array
181 //
182 numFace = phiIsOpen ? numCorner+2 : numCorner;
183 faces = new G4VCSGface*[numFace];
184
185 //
186 // Construct conical faces
187 //
188 // But! Don't construct a face if both points are at zero radius!
189 //
190 G4PolyconeSideRZ *corner = corners,
191 *prev = corners + numCorner-1,
192 *nextNext;
193 G4VCSGface** face = faces;
194 do // Loop checking, 13.08.2015, G.Cosmo
195 {
196 next = corner+1;
197 if (next >= corners+numCorner) next = corners;
198 nextNext = next+1;
199 if (nextNext >= corners+numCorner) nextNext = corners;
200
201 if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue;
202
203 //
204 // We must decide here if we can dare declare one of our faces
205 // as having a "valid" normal (i.e. allBehind = true). This
206 // is never possible if the face faces "inward" in r.
207 //
208 G4bool allBehind;
209 if (corner->z > next->z)
210 {
211 allBehind = false;
212 }
213 else
214 {
215 //
216 // Otherwise, it is only true if the line passing
217 // through the two points of the segment do not
218 // split the r/z cross section
219 //
220 allBehind = !rz->BisectedBy( corner->r, corner->z,
221 next->r, next->z, kCarTolerance );
222 }
223
224 *face++ = new G4PolyconeSide( prev, corner, next, nextNext,
225 startPhi, endPhi-startPhi, phiIsOpen, allBehind );
226 } while( prev=corner, corner=next, corner > corners );
227
228 if (phiIsOpen)
229 {
230 //
231 // Construct phi open edges
232 //
233 *face++ = new G4PolyPhiFace( rz, startPhi, 0, endPhi );
234 *face++ = new G4PolyPhiFace( rz, endPhi, 0, startPhi );
235 }
236
237 //
238 // We might have dropped a face or two: recalculate numFace
239 //
240 numFace = (G4int)(face-faces);
241
242 //
243 // Make enclosingCylinder
244 //
245 enclosingCylinder =
246 new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal );
247 }
248
249 // Fake default constructor - sets only member data and allocates memory
250 // for usage restricted to object persistency.
251 //
252 G4GenericPolycone::G4GenericPolycone( __void__& a )
253 : G4VCSGfaceted(a), startPhi(0.), endPhi(0.), numCorner(0)
254 {
255 }
256
257 // Destructor
258 //
259 G4GenericPolycone::~G4GenericPolycone()
260 {
261 delete [] corners;
262 delete enclosingCylinder;
263 delete fElements;
264 delete fpPolyhedron;
265 corners = nullptr;
266 enclosingCylinder = nullptr;
267 fElements = nullptr;
268 fpPolyhedron = nullptr;
269 }
270
271 // Copy constructor
272 //
273 G4GenericPolycone::G4GenericPolycone( const G4GenericPolycone& source )
274 : G4VCSGfaceted( source )
275 {
276 CopyStuff( source );
277 }
278
279 // Assignment operator
280 //
281 G4GenericPolycone&
282 G4GenericPolycone::operator=( const G4GenericPolycone& source )
283 {
284 if (this == &source) return *this;
285
286 G4VCSGfaceted::operator=( source );
287
288 delete [] corners;
289 // if (original_parameters) delete original_parameters;
290
291 delete enclosingCylinder;
292
293 CopyStuff( source );
294
295 return *this;
296 }
297
298 // CopyStuff
299 //
300 void G4GenericPolycone::CopyStuff( const G4GenericPolycone& source )
301 {
302 //
303 // Simple stuff
304 //
305 startPhi = source.startPhi;
306 endPhi = source.endPhi;
307 phiIsOpen = source.phiIsOpen;
308 numCorner = source.numCorner;
309
310 //
311 // The corner array
312 //
313 corners = new G4PolyconeSideRZ[numCorner];
314
315 G4PolyconeSideRZ *corn = corners,
316 *sourceCorn = source.corners;
317 do // Loop checking, 13.08.2015, G.Cosmo
318 {
319 *corn = *sourceCorn;
320 } while( ++sourceCorn, ++corn < corners+numCorner );
321
322 //
323 // Enclosing cylinder
324 //
325 enclosingCylinder = new G4EnclosingCylinder( *source.enclosingCylinder );
326
327 //
328 // Surface elements
329 //
330 delete fElements;
331 fElements = nullptr;
332
333 // Polyhedron
334 //
335 fRebuildPolyhedron = false;
336 delete fpPolyhedron;
337 fpPolyhedron = nullptr;
338 }
339
340 // Reset
341 //
342 G4bool G4GenericPolycone::Reset()
343 {
344 std::ostringstream message;
345 message << "Solid " << GetName() << " built using generic construct."
346 << G4endl << "Not applicable to the generic construct !";
347 G4Exception("G4GenericPolycone::Reset()", "GeomSolids1001",
348 JustWarning, message, "Parameters NOT resetted.");
349 return true;
350 }
351
352 // Inside
353 //
354 // This is an override of G4VCSGfaceted::Inside, created in order
355 // to speed things up by first checking with G4EnclosingCylinder.
356 //
357 EInside G4GenericPolycone::Inside( const G4ThreeVector& p ) const
358 {
359 //
360 // Quick test
361 //
362 if (enclosingCylinder->MustBeOutside(p)) return kOutside;
363
364 //
365 // Long answer
366 //
367 return G4VCSGfaceted::Inside(p);
368 }
369
370 // DistanceToIn
371 //
372 // This is an override of G4VCSGfaceted::Inside, created in order
373 // to speed things up by first checking with G4EnclosingCylinder.
374 //
375 G4double G4GenericPolycone::DistanceToIn( const G4ThreeVector& p,
376 const G4ThreeVector& v ) const
377 {
378 //
379 // Quick test
380 //
381 if (enclosingCylinder->ShouldMiss(p,v))
382 return kInfinity;
383
384 //
385 // Long answer
386 //
387 return G4VCSGfaceted::DistanceToIn( p, v );
388 }
389
390 // DistanceToIn
391 //
392 G4double G4GenericPolycone::DistanceToIn( const G4ThreeVector& p ) const
393 {
394 return G4VCSGfaceted::DistanceToIn(p);
395 }
396
397 // BoundingLimits
398 //
399 // Get bounding box
400 //
401 void
402 G4GenericPolycone::BoundingLimits(G4ThreeVector& pMin,
403 G4ThreeVector& pMax) const
404 {
405 G4double rmin = kInfinity, rmax = -kInfinity;
406 G4double zmin = kInfinity, zmax = -kInfinity;
407
408 for (G4int i=0; i<GetNumRZCorner(); ++i)
409 {
410 G4PolyconeSideRZ corner = GetCorner(i);
411 if (corner.r < rmin) rmin = corner.r;
412 if (corner.r > rmax) rmax = corner.r;
413 if (corner.z < zmin) zmin = corner.z;
414 if (corner.z > zmax) zmax = corner.z;
415 }
416
417 if (IsOpen())
418 {
419 G4TwoVector vmin,vmax;
420 G4GeomTools::DiskExtent(rmin,rmax,
421 GetSinStartPhi(),GetCosStartPhi(),
422 GetSinEndPhi(),GetCosEndPhi(),
423 vmin,vmax);
424 pMin.set(vmin.x(),vmin.y(),zmin);
425 pMax.set(vmax.x(),vmax.y(),zmax);
426 }
427 else
428 {
429 pMin.set(-rmax,-rmax, zmin);
430 pMax.set( rmax, rmax, zmax);
431 }
432
433 // Check correctness of the bounding box
434 //
435 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
436 {
437 std::ostringstream message;
438 message << "Bad bounding box (min >= max) for solid: "
439 << GetName() << " !"
440 << "\npMin = " << pMin
441 << "\npMax = " << pMax;
442 G4Exception("GenericG4Polycone::BoundingLimits()", "GeomMgt0001",
443 JustWarning, message);
444 DumpInfo();
445 }
446 }
447
448 // CalculateExtent
449 //
450 // Calculate extent under transform and specified limit
451 //
452 G4bool
453 G4GenericPolycone::CalculateExtent(const EAxis pAxis,
454 const G4VoxelLimits& pVoxelLimit,
455 const G4AffineTransform& pTransform,
456 G4double& pMin, G4double& pMax) const
457 {
458 G4ThreeVector bmin, bmax;
459 G4bool exist;
460
461 // Check bounding box (bbox)
462 //
463 BoundingLimits(bmin,bmax);
464 G4BoundingEnvelope bbox(bmin,bmax);
465 #ifdef G4BBOX_EXTENT
466 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
467 #endif
468 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
469 {
470 return exist = pMin < pMax;
471 }
472
473 // To find the extent, RZ contour of the polycone is subdivided
474 // in triangles. The extent is calculated as cumulative extent of
475 // all sub-polycones formed by rotation of triangles around Z
476 //
477 G4TwoVectorList contourRZ;
478 G4TwoVectorList triangles;
479 G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
480 G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
481
482 // get RZ contour, ensure anticlockwise order of corners
483 for (G4int i=0; i<GetNumRZCorner(); ++i)
484 {
485 G4PolyconeSideRZ corner = GetCorner(i);
486 contourRZ.emplace_back(corner.r,corner.z);
487 }
488 G4double area = G4GeomTools::PolygonArea(contourRZ);
489 if (area < 0.) std::reverse(contourRZ.begin(),contourRZ.end());
490
491 // triangulate RZ countour
492 if (!G4GeomTools::TriangulatePolygon(contourRZ,triangles))
493 {
494 std::ostringstream message;
495 message << "Triangulation of RZ contour has failed for solid: "
496 << GetName() << " !"
497 << "\nExtent has been calculated using boundary box";
498 G4Exception("G4GenericPolycone::CalculateExtent()",
499 "GeomMgt1002", JustWarning, message);
500 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
501 }
502
503 // set trigonometric values
504 const G4int NSTEPS = 24; // number of steps for whole circle
505 G4double astep = twopi/NSTEPS; // max angle for one step
506
507 G4double sphi = GetStartPhi();
508 G4double ephi = GetEndPhi();
509 G4double dphi = IsOpen() ? ephi-sphi : twopi;
510 G4int ksteps = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1;
511 G4double ang = dphi/ksteps;
512
513 G4double sinHalf = std::sin(0.5*ang);
514 G4double cosHalf = std::cos(0.5*ang);
515 G4double sinStep = 2.*sinHalf*cosHalf;
516 G4double cosStep = 1. - 2.*sinHalf*sinHalf;
517
518 G4double sinStart = GetSinStartPhi();
519 G4double cosStart = GetCosStartPhi();
520 G4double sinEnd = GetSinEndPhi();
521 G4double cosEnd = GetCosEndPhi();
522
523 // define vectors and arrays
524 std::vector<const G4ThreeVectorList *> polygons;
525 polygons.resize(ksteps+2);
526 G4ThreeVectorList pols[NSTEPS+2];
527 for (G4int k=0; k<ksteps+2; ++k) pols[k].resize(6);
528 for (G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k];
529 G4double r0[6],z0[6]; // contour with original edges of triangle
530 G4double r1[6]; // shifted radii of external edges of triangle
531
532 // main loop along triangles
533 pMin = kInfinity;
534 pMax =-kInfinity;
535 G4int ntria = (G4int)triangles.size()/3;
536 for (G4int i=0; i<ntria; ++i)
537 {
538 G4int i3 = i*3;
539 for (G4int k=0; k<3; ++k)
540 {
541 G4int e0 = i3+k, e1 = (k<2) ? e0+1 : i3;
542 G4int k2 = k*2;
543 // set contour with original edges of triangle
544 r0[k2+0] = triangles[e0].x(); z0[k2+0] = triangles[e0].y();
545 r0[k2+1] = triangles[e1].x(); z0[k2+1] = triangles[e1].y();
546 // set shifted radii
547 r1[k2+0] = r0[k2+0];
548 r1[k2+1] = r0[k2+1];
549 if (z0[k2+1] - z0[k2+0] <= 0) continue;
550 r1[k2+0] /= cosHalf;
551 r1[k2+1] /= cosHalf;
552 }
553
554 // rotate countour, set sequence of 6-sided polygons
555 G4double sinCur = sinStart*cosHalf + cosStart*sinHalf;
556 G4double cosCur = cosStart*cosHalf - sinStart*sinHalf;
557 for (G4int j=0; j<6; ++j)
558 {
559 pols[0][j].set(r0[j]*cosStart,r0[j]*sinStart,z0[j]);
560 }
561 for (G4int k=1; k<ksteps+1; ++k)
562 {
563 for (G4int j=0; j<6; ++j)
564 {
565 pols[k][j].set(r1[j]*cosCur,r1[j]*sinCur,z0[j]);
566 }
567 G4double sinTmp = sinCur;
568 sinCur = sinCur*cosStep + cosCur*sinStep;
569 cosCur = cosCur*cosStep - sinTmp*sinStep;
570 }
571 for (G4int j=0; j<6; ++j)
572 {
573 pols[ksteps+1][j].set(r0[j]*cosEnd,r0[j]*sinEnd,z0[j]);
574 }
575
576 // set sub-envelope and adjust extent
577 G4double emin,emax;
578 G4BoundingEnvelope benv(polygons);
579 if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
580 if (emin < pMin) pMin = emin;
581 if (emax > pMax) pMax = emax;
582 if (eminlim > pMin && emaxlim < pMax) return true; // max possible extent
583 }
584 return (pMin < pMax);
585 }
586
587 // GetEntityType
588 //
589 G4GeometryType G4GenericPolycone::GetEntityType() const
590 {
591 return {"G4GenericPolycone"};
592 }
593
594 // Clone
595 //
596 // Make a clone of the object
597 //
598 G4VSolid* G4GenericPolycone::Clone() const
599 {
600 return new G4GenericPolycone(*this);
601 }
602
603 // StreamInfo
604 //
605 // Stream object contents to an output stream
606 //
607 std::ostream& G4GenericPolycone::StreamInfo( std::ostream& os ) const
608 {
609 G4long oldprc = os.precision(16);
610 os << "-----------------------------------------------------------\n"
611 << " *** Dump for solid - " << GetName() << " ***\n"
612 << " ===================================================\n"
613 << " Solid type: G4GenericPolycone\n"
614 << " Parameters: \n"
615 << " starting phi angle : " << startPhi/degree << " degrees \n"
616 << " ending phi angle : " << endPhi/degree << " degrees \n";
617 G4int i=0;
618
619 os << " number of RZ points: " << numCorner << "\n"
620 << " RZ values (corners): \n";
621 for (i=0; i<numCorner; i++)
622 {
623 os << " "
624 << corners[i].r << ", " << corners[i].z << "\n";
625 }
626 os << "-----------------------------------------------------------\n";
627 os.precision(oldprc);
628
629 return os;
630 }
631
632 //////////////////////////////////////////////////////////////////////////
633 //
634 // Return volume
635
636 G4double G4GenericPolycone::GetCubicVolume()
637 {
638 if (fCubicVolume == 0.)
639 {
640 G4double total = 0.;
641 G4int nrz = GetNumRZCorner();
642 G4PolyconeSideRZ a = GetCorner(nrz - 1);
643 for (G4int i=0; i<nrz; ++i)
644 {
645 G4PolyconeSideRZ b = GetCorner(i);
646 total += (b.r*b.r + b.r*a.r + a.r*a.r)*(b.z - a.z);
647 a = b;
648 }
649 fCubicVolume = std::abs(total)*(GetEndPhi() - GetStartPhi())/6.;
650 }
651 return fCubicVolume;
652 }
653
654 //////////////////////////////////////////////////////////////////////////
655 //
656 // Return surface area
657
658 G4double G4GenericPolycone::GetSurfaceArea()
659 {
660 if (fSurfaceArea == 0.)
661 {
662 // phi cut area
663 G4int nrz = GetNumRZCorner();
664 G4double scut = 0.;
665 if (IsOpen())
666 {
667 G4PolyconeSideRZ a = GetCorner(nrz - 1);
668 for (G4int i=0; i<nrz; ++i)
669 {
670 G4PolyconeSideRZ b = GetCorner(i);
671 scut += a.r*b.z - a.z*b.r;
672 a = b;
673 }
674 scut = std::abs(scut);
675 }
676 // lateral surface area
677 G4double slat = 0;
678 G4PolyconeSideRZ a = GetCorner(nrz - 1);
679 for (G4int i=0; i<nrz; ++i)
680 {
681 G4PolyconeSideRZ b = GetCorner(i);
682 G4double h = std::sqrt((b.r - a.r)*(b.r - a.r) + (b.z - a.z)*(b.z - a.z));
683 slat += (b.r + a.r)*h;
684 a = b;
685 }
686 slat *= (GetEndPhi() - GetStartPhi())/2.;
687 fSurfaceArea = scut + slat;
688 }
689 return fSurfaceArea;
690 }
691
692 //////////////////////////////////////////////////////////////////////////
693 //
694 // Set vector of surface elements, auxiliary method for sampling
695 // random points on surface
696
697 void G4GenericPolycone::SetSurfaceElements() const
698 {
699 fElements = new std::vector<G4GenericPolycone::surface_element>;
700 G4double sarea = 0.;
701 G4int nrz = GetNumRZCorner();
702
703 // set lateral surface elements
704 G4double dphi = GetEndPhi() - GetStartPhi();
705 G4int ia = nrz - 1;
706 for (G4int ib=0; ib<nrz; ++ib)
707 {
708 G4PolyconeSideRZ a = GetCorner(ia);
709 G4PolyconeSideRZ b = GetCorner(ib);
710 G4GenericPolycone::surface_element selem;
711 selem.i0 = ia;
712 selem.i1 = ib;
713 selem.i2 = -1;
714 ia = ib;
715 if (a.r == 0. && b.r == 0.) continue;
716 G4double h = std::sqrt((b.r - a.r)*(b.r - a.r) + (b.z - a.z)*(b.z - a.z));
717 sarea += 0.5*dphi*(b.r + a.r)*h;
718 selem.area = sarea;
719 fElements->push_back(selem);
720 }
721
722 // set elements for phi cuts
723 if (IsOpen())
724 {
725 G4TwoVectorList contourRZ;
726 std::vector<G4int> triangles;
727 for (G4int i=0; i<nrz; ++i)
728 {
729 G4PolyconeSideRZ corner = GetCorner(i);
730 contourRZ.emplace_back(corner.r, corner.z);
731 }
732 G4GeomTools::TriangulatePolygon(contourRZ, triangles);
733 auto ntria = (G4int)triangles.size();
734 for (G4int i=0; i<ntria; i+=3)
735 {
736 G4GenericPolycone::surface_element selem;
737 selem.i0 = triangles[i];
738 selem.i1 = triangles[i+1];
739 selem.i2 = triangles[i+2];
740 G4PolyconeSideRZ a = GetCorner(selem.i0);
741 G4PolyconeSideRZ b = GetCorner(selem.i1);
742 G4PolyconeSideRZ c = GetCorner(selem.i2);
743 G4double stria =
744 std::abs(G4GeomTools::TriangleArea(a.r, a.z, b.r, b.z, c.r, c.z));
745 sarea += stria;
746 selem.area = sarea;
747 fElements->push_back(selem); // start phi
748 sarea += stria;
749 selem.area = sarea;
750 selem.i0 += nrz;
751 fElements->push_back(selem); // end phi
752 }
753 }
754 }
755
756 //////////////////////////////////////////////////////////////////////////
757 //
758 // Generate random point on surface
759
760 G4ThreeVector G4GenericPolycone::GetPointOnSurface() const
761 {
762 // Set surface elements
763 if (fElements == nullptr)
764 {
765 G4AutoLock l(&surface_elementsMutex);
766 SetSurfaceElements();
767 l.unlock();
768 }
769
770 // Select surface element
771 G4GenericPolycone::surface_element selem;
772 selem = fElements->back();
773 G4double select = selem.area*G4QuickRand();
774 auto it = std::lower_bound(fElements->begin(), fElements->end(), select,
775 [](const G4GenericPolycone::surface_element& x, G4double val)
776 -> G4bool { return x.area < val; });
777
778 // Generate random point
779 G4double r = 0, z = 0, phi = 0;
780 G4double u = G4QuickRand();
781 G4double v = G4QuickRand();
782 G4int i0 = (*it).i0;
783 G4int i1 = (*it).i1;
784 G4int i2 = (*it).i2;
785 if (i2 < 0) // lateral surface
786 {
787 G4PolyconeSideRZ p0 = GetCorner(i0);
788 G4PolyconeSideRZ p1 = GetCorner(i1);
789 if (p1.r < p0.r)
790 {
791 p0 = GetCorner(i1);
792 p1 = GetCorner(i0);
793 }
794 if (p1.r - p0.r < kCarTolerance) // cylindrical surface
795 {
796 r = (p1.r - p0.r)*u + p0.r;
797 z = (p1.z - p0.z)*u + p0.z;
798 }
799 else // conical surface
800 {
801 r = std::sqrt(p1.r*p1.r*u + p0.r*p0.r*(1. - u));
802 z = p0.z + (p1.z - p0.z)*(r - p0.r)/(p1.r - p0.r);
803 }
804 phi = (GetEndPhi() - GetStartPhi())*v + GetStartPhi();
805 }
806 else // phi cut
807 {
808 G4int nrz = GetNumRZCorner();
809 phi = (i0 < nrz) ? GetStartPhi() : GetEndPhi();
810 if (i0 >= nrz) { i0 -= nrz; }
811 G4PolyconeSideRZ p0 = GetCorner(i0);
812 G4PolyconeSideRZ p1 = GetCorner(i1);
813 G4PolyconeSideRZ p2 = GetCorner(i2);
814 if (u + v > 1.) { u = 1. - u; v = 1. - v; }
815 r = (p1.r - p0.r)*u + (p2.r - p0.r)*v + p0.r;
816 z = (p1.z - p0.z)*u + (p2.z - p0.z)*v + p0.z;
817 }
818 return {r*std::cos(phi), r*std::sin(phi), z};
819 }
820
821 //////////////////////////////////////////////////////////////////////////
822 //
823 // CreatePolyhedron
824
825 G4Polyhedron* G4GenericPolycone::CreatePolyhedron() const
826 {
827 std::vector<G4TwoVector> rz(numCorner);
828 for (G4int i = 0; i < numCorner; ++i)
829 rz[i].set(corners[i].r, corners[i].z);
830 return new G4PolyhedronPcon(startPhi, endPhi - startPhi, rz);
831 }
832
833 #endif
834