ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4UPolyhedra.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4UPolyhedra.cc
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 G4UPolycone wrapper class
27 //
28 // 31.10.13 G.Cosmo, CERN
29 // --------------------------------------------------------------------
30 
31 #include "G4Polyhedra.hh"
32 #include "G4UPolyhedra.hh"
33 
34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G4GEOM_USE_PARTIAL_USOLIDS) )
35 
36 #include "G4GeomTools.hh"
37 #include "G4GeometryTolerance.hh"
38 #include "G4AffineTransform.hh"
39 #include "G4VPVParameterisation.hh"
40 #include "G4BoundingEnvelope.hh"
41 
42 using namespace CLHEP;
43 
45 //
46 // Constructor (GEANT3 style parameters)
47 //
48 // GEANT3 PGON radii are specified in the distance to the norm of each face.
49 //
50 G4UPolyhedra::G4UPolyhedra(const G4String& name,
51  G4double phiStart,
52  G4double phiTotal,
53  G4int numSide,
54  G4int numZPlanes,
55  const G4double zPlane[],
56  const G4double rInner[],
57  const G4double rOuter[] )
58  : Base_t(name, phiStart, phiTotal, numSide,
59  numZPlanes, zPlane, rInner, rOuter)
60 {
61  fGenericPgon = false;
62  SetOriginalParameters();
63  wrStart = phiStart;
64  while (wrStart < 0)
65  {
66  wrStart += twopi;
67  }
68  wrDelta = phiTotal;
69  if (wrDelta <= 0. || wrDelta >= twopi*(1-DBL_EPSILON))
70  {
71  wrDelta = twopi;
72  }
73  wrNumSide = numSide;
74  G4double convertRad = 1./std::cos(0.5*wrDelta/wrNumSide);
75  rzcorners.resize(0);
76  for (G4int i=0; i<numZPlanes; ++i)
77  {
78  G4double z = zPlane[i];
79  G4double r = rOuter[i]*convertRad;
80  rzcorners.push_back(G4TwoVector(r,z));
81  }
82  for (G4int i=numZPlanes-1; i>=0; --i)
83  {
84  G4double z = zPlane[i];
85  G4double r = rInner[i]*convertRad;
86  rzcorners.push_back(G4TwoVector(r,z));
87  }
88  std::vector<G4int> iout;
89  G4GeomTools::RemoveRedundantVertices(rzcorners,iout,2*kCarTolerance);
90 }
91 
92 
94 //
95 // Constructor (generic parameters)
96 //
97 G4UPolyhedra::G4UPolyhedra(const G4String& name,
98  G4double phiStart,
99  G4double phiTotal,
100  G4int numSide,
101  G4int numRZ,
102  const G4double r[],
103  const G4double z[] )
104  : Base_t(name, phiStart, phiTotal, numSide, numRZ, r, z)
105 {
106  fGenericPgon = true;
107  SetOriginalParameters();
108  wrStart = phiStart;
109  while (wrStart < 0.)
110  {
111  wrStart += twopi;
112  }
113  wrDelta = phiTotal;
114  if (wrDelta <= 0. || wrDelta >= twopi*(1-DBL_EPSILON))
115  {
116  wrDelta = twopi;
117  }
118  wrNumSide = numSide;
119  G4double convertRad = 1./std::cos(0.5*wrDelta/wrNumSide);
120  rzcorners.resize(0);
121  for (G4int i=0; i<numRZ; ++i)
122  {
123  rzcorners.push_back(G4TwoVector(r[i]*convertRad,z[i]));
124  }
125  std::vector<G4int> iout;
126  G4GeomTools::RemoveRedundantVertices(rzcorners,iout,2*kCarTolerance);
127 }
128 
129 
131 //
132 // Fake default constructor - sets only member data and allocates memory
133 // for usage restricted to object persistency.
134 //
135 G4UPolyhedra::G4UPolyhedra( __void__& a )
136  : Base_t(a)
137 {
138 }
139 
140 
142 //
143 // Destructor
144 //
145 G4UPolyhedra::~G4UPolyhedra()
146 {
147 }
148 
149 
151 //
152 // Copy constructor
153 //
154 G4UPolyhedra::G4UPolyhedra( const G4UPolyhedra& source )
155  : Base_t( source )
156 {
157  fGenericPgon = source.fGenericPgon;
158  fOriginalParameters = source.fOriginalParameters;
159  wrStart = source.wrStart;
160  wrDelta = source.wrDelta;
161  wrNumSide = source.wrNumSide;
162  rzcorners = source.rzcorners;
163 }
164 
165 
167 //
168 // Assignment operator
169 //
170 G4UPolyhedra& G4UPolyhedra::operator=( const G4UPolyhedra& source )
171 {
172  if (this == &source) return *this;
173 
174  Base_t::operator=( source );
175  fGenericPgon = source.fGenericPgon;
176  fOriginalParameters = source.fOriginalParameters;
177  wrStart = source.wrStart;
178  wrDelta = source.wrDelta;
179  wrNumSide = source.wrNumSide;
180  rzcorners = source.rzcorners;
181 
182  return *this;
183 }
184 
185 
187 //
188 // Accessors & modifiers
189 //
190 G4int G4UPolyhedra::GetNumSide() const
191 {
192  return wrNumSide;
193 }
194 G4double G4UPolyhedra::GetStartPhi() const
195 {
196  return wrStart;
197 }
198 G4double G4UPolyhedra::GetEndPhi() const
199 {
200  return (wrStart + wrDelta);
201 }
202 G4double G4UPolyhedra::GetSinStartPhi() const
203 {
204  G4double phi = GetStartPhi();
205  return std::sin(phi);
206 }
207 G4double G4UPolyhedra::GetCosStartPhi() const
208 {
209  G4double phi = GetStartPhi();
210  return std::cos(phi);
211 }
212 G4double G4UPolyhedra::GetSinEndPhi() const
213 {
214  G4double phi = GetEndPhi();
215  return std::sin(phi);
216 }
217 G4double G4UPolyhedra::GetCosEndPhi() const
218 {
219  G4double phi = GetEndPhi();
220  return std::cos(phi);
221 }
222 G4bool G4UPolyhedra::IsOpen() const
223 {
224  return (wrDelta < twopi);
225 }
226 G4bool G4UPolyhedra::IsGeneric() const
227 {
228  return fGenericPgon;
229 }
230 G4int G4UPolyhedra::GetNumRZCorner() const
231 {
232  return rzcorners.size();
233 }
234 G4PolyhedraSideRZ G4UPolyhedra::GetCorner(G4int index) const
235 {
236  G4TwoVector rz = rzcorners.at(index);
237  G4PolyhedraSideRZ psiderz = { rz.x(), rz.y() };
238 
239  return psiderz;
240 }
241 G4PolyhedraHistorical* G4UPolyhedra::GetOriginalParameters() const
242 {
243  return new G4PolyhedraHistorical(fOriginalParameters);
244 }
245 void G4UPolyhedra::SetOriginalParameters()
246 {
247  G4double startPhi = GetPhiStart();
248  G4double deltaPhi = GetPhiDelta();
249  G4int numPlanes = GetZSegmentCount() + 1;
250  G4int numSides = GetSideCount();
251 
252  fOriginalParameters.Start_angle = startPhi;
253  fOriginalParameters.Opening_angle = deltaPhi;
254  fOriginalParameters.Num_z_planes = numPlanes;
255  fOriginalParameters.numSide = numSides;
256 
257  delete [] fOriginalParameters.Z_values;
258  delete [] fOriginalParameters.Rmin;
259  delete [] fOriginalParameters.Rmax;
260  fOriginalParameters.Z_values = new G4double[numPlanes];
261  fOriginalParameters.Rmin = new G4double[numPlanes];
262  fOriginalParameters.Rmax = new G4double[numPlanes];
263 
264  G4double convertRad = std::cos(0.5*deltaPhi/numSides);
265  for (G4int i=0; i<numPlanes; ++i)
266  {
267  fOriginalParameters.Z_values[i] = GetZPlanes()[i];
268  fOriginalParameters.Rmax[i] = GetRMax()[i]/convertRad;
269  fOriginalParameters.Rmin[i] = GetRMin()[i]/convertRad;
270  }
271 }
272 void G4UPolyhedra::SetOriginalParameters(G4PolyhedraHistorical* pars)
273 {
274  fOriginalParameters = *pars;
275  fRebuildPolyhedron = true;
276  Reset();
277 }
278 
279 G4bool G4UPolyhedra::Reset()
280 {
281  if (fGenericPgon)
282  {
283  std::ostringstream message;
284  message << "Solid " << GetName() << " built using generic construct."
285  << G4endl << "Not applicable to the generic construct !";
286  G4Exception("G4UPolyhedra::Reset()", "GeomSolids1001",
287  JustWarning, message, "Parameters NOT resetted.");
288  return true; // error code set
289  }
290 
291  //
292  // Rebuild polyhedra based on original parameters
293  //
294  wrStart = fOriginalParameters.Start_angle;
295  while (wrStart < 0.)
296  {
297  wrStart += twopi;
298  }
299  wrDelta = fOriginalParameters.Opening_angle;
300  if (wrDelta <= 0. || wrDelta >= twopi*(1-DBL_EPSILON))
301  {
302  wrDelta = twopi;
303  }
304  wrNumSide = fOriginalParameters.numSide;
305  rzcorners.resize(0);
306  for (G4int i=0; i<fOriginalParameters.Num_z_planes; ++i)
307  {
308  G4double z = fOriginalParameters.Z_values[i];
309  G4double r = fOriginalParameters.Rmax[i];
310  rzcorners.push_back(G4TwoVector(r,z));
311  }
312  for (G4int i=fOriginalParameters.Num_z_planes-1; i>=0; --i)
313  {
314  G4double z = fOriginalParameters.Z_values[i];
315  G4double r = fOriginalParameters.Rmin[i];
316  rzcorners.push_back(G4TwoVector(r,z));
317  }
318  std::vector<G4int> iout;
319  G4GeomTools::RemoveRedundantVertices(rzcorners,iout,2*kCarTolerance);
320 
321  return false; // error code unset
322 }
323 
324 
326 //
327 // Dispatch to parameterisation for replication mechanism dimension
328 // computation & modification.
329 //
330 void G4UPolyhedra::ComputeDimensions(G4VPVParameterisation* p,
331  const G4int n,
332  const G4VPhysicalVolume* pRep)
333 {
334  p->ComputeDimensions(*(G4Polyhedra*)this,n,pRep);
335 }
336 
337 
339 //
340 // Make a clone of the object
341 
342 G4VSolid* G4UPolyhedra::Clone() const
343 {
344  return new G4UPolyhedra(*this);
345 }
346 
347 
349 //
350 // Get bounding box
351 
352 void G4UPolyhedra::BoundingLimits(G4ThreeVector& pMin,
353  G4ThreeVector& pMax) const
354 {
355  static G4bool checkBBox = true;
356  static G4bool checkPhi = true;
357 
358  G4double rmin = kInfinity, rmax = -kInfinity;
359  G4double zmin = kInfinity, zmax = -kInfinity;
360  for (G4int i=0; i<GetNumRZCorner(); ++i)
361  {
362  G4PolyhedraSideRZ corner = GetCorner(i);
363  if (corner.r < rmin) rmin = corner.r;
364  if (corner.r > rmax) rmax = corner.r;
365  if (corner.z < zmin) zmin = corner.z;
366  if (corner.z > zmax) zmax = corner.z;
367  }
368 
369  G4double sphi = GetStartPhi();
370  G4double ephi = GetEndPhi();
371  G4double dphi = IsOpen() ? ephi-sphi : twopi;
372  G4int ksteps = GetNumSide();
373  G4double astep = dphi/ksteps;
374  G4double sinStep = std::sin(astep);
375  G4double cosStep = std::cos(astep);
376 
377  G4double sinCur = GetSinStartPhi();
378  G4double cosCur = GetCosStartPhi();
379  if (!IsOpen()) rmin = 0.;
380  G4double xmin = rmin*cosCur, xmax = xmin;
381  G4double ymin = rmin*sinCur, ymax = ymin;
382  for (G4int k=0; k<ksteps+1; ++k)
383  {
384  G4double x = rmax*cosCur;
385  if (x < xmin) xmin = x;
386  if (x > xmax) xmax = x;
387  G4double y = rmax*sinCur;
388  if (y < ymin) ymin = y;
389  if (y > ymax) ymax = y;
390  if (rmin > 0.)
391  {
392  G4double xx = rmin*cosCur;
393  if (xx < xmin) xmin = xx;
394  if (xx > xmax) xmax = xx;
395  G4double yy = rmin*sinCur;
396  if (yy < ymin) ymin = yy;
397  if (yy > ymax) ymax = yy;
398  }
399  G4double sinTmp = sinCur;
400  sinCur = sinCur*cosStep + cosCur*sinStep;
401  cosCur = cosCur*cosStep - sinTmp*sinStep;
402  }
403  pMin.set(xmin,ymin,zmin);
404  pMax.set(xmax,ymax,zmax);
405 
406  // Check correctness of the bounding box
407  //
408  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
409  {
410  std::ostringstream message;
411  message << "Bad bounding box (min >= max) for solid: "
412  << GetName() << " !"
413  << "\npMin = " << pMin
414  << "\npMax = " << pMax;
415  G4Exception("G4UPolyhedra::BoundingLimits()", "GeomMgt0001",
416  JustWarning, message);
417  StreamInfo(G4cout);
418  }
419 
420  // Check consistency of bounding boxes
421  //
422  if (checkBBox)
423  {
424  U3Vector vmin, vmax;
425  Extent(vmin,vmax);
426  if (std::abs(pMin.x()-vmin.x()) > kCarTolerance ||
427  std::abs(pMin.y()-vmin.y()) > kCarTolerance ||
428  std::abs(pMin.z()-vmin.z()) > kCarTolerance ||
429  std::abs(pMax.x()-vmax.x()) > kCarTolerance ||
430  std::abs(pMax.y()-vmax.y()) > kCarTolerance ||
431  std::abs(pMax.z()-vmax.z()) > kCarTolerance)
432  {
433  std::ostringstream message;
434  message << "Inconsistency in bounding boxes for solid: "
435  << GetName() << " !"
436  << "\nBBox min: wrapper = " << pMin << " solid = " << vmin
437  << "\nBBox max: wrapper = " << pMax << " solid = " << vmax;
438  G4Exception("G4UPolyhedra::BoundingLimits()", "GeomMgt0001",
439  JustWarning, message);
440  checkBBox = false;
441  }
442  }
443 
444  // Check consistency of angles
445  //
446  if (checkPhi)
447  {
448  if (GetStartPhi() != GetPhiStart() ||
449  GetEndPhi() != GetPhiEnd() ||
450  GetNumSide() != GetSideCount() ||
451  IsOpen() != (Base_t::GetPhiDelta() < twopi))
452  {
453  std::ostringstream message;
454  message << "Inconsistency in Phi angles or # of sides for solid: "
455  << GetName() << " !"
456  << "\nPhi start : wrapper = " << GetStartPhi()
457  << " solid = " << GetPhiStart()
458  << "\nPhi end : wrapper = " << GetEndPhi()
459  << " solid = " << GetPhiEnd()
460  << "\nPhi # sides: wrapper = " << GetNumSide()
461  << " solid = " << GetSideCount()
462  << "\nPhi is open: wrapper = " << (IsOpen() ? "true" : "false")
463  << " solid = "
464  << ((Base_t::GetPhiDelta() < twopi) ? "true" : "false");
465  G4Exception("G4UPolyhedra::BoundingLimits()", "GeomMgt0001",
466  JustWarning, message);
467  checkPhi = false;
468  }
469  }
470 }
471 
473 //
474 // Calculate extent under transform and specified limit
475 
476 G4bool
477 G4UPolyhedra::CalculateExtent(const EAxis pAxis,
478  const G4VoxelLimits& pVoxelLimit,
479  const G4AffineTransform& pTransform,
480  G4double& pMin, G4double& pMax) const
481 {
482  G4ThreeVector bmin, bmax;
483  G4bool exist;
484 
485  // Check bounding box (bbox)
486  //
487  BoundingLimits(bmin,bmax);
488  G4BoundingEnvelope bbox(bmin,bmax);
489 #ifdef G4BBOX_EXTENT
490  if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
491 #endif
492  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
493  {
494  return exist = (pMin < pMax) ? true : false;
495  }
496 
497  // To find the extent, RZ contour of the polycone is subdivided
498  // in triangles. The extent is calculated as cumulative extent of
499  // all sub-polycones formed by rotation of triangles around Z
500  //
501  G4TwoVectorList contourRZ;
502  G4TwoVectorList triangles;
503  std::vector<G4int> iout;
504  G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
505  G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
506 
507  // get RZ contour, ensure anticlockwise order of corners
508  for (G4int i=0; i<GetNumRZCorner(); ++i)
509  {
510  G4PolyhedraSideRZ corner = GetCorner(i);
511  contourRZ.push_back(G4TwoVector(corner.r,corner.z));
512  }
513  G4GeomTools::RemoveRedundantVertices(contourRZ,iout,2*kCarTolerance);
514  G4double area = G4GeomTools::PolygonArea(contourRZ);
515  if (area < 0.) std::reverse(contourRZ.begin(),contourRZ.end());
516 
517  // triangulate RZ countour
518  if (!G4GeomTools::TriangulatePolygon(contourRZ,triangles))
519  {
520  std::ostringstream message;
521  message << "Triangulation of RZ contour has failed for solid: "
522  << GetName() << " !"
523  << "\nExtent has been calculated using boundary box";
524  G4Exception("G4UPolyhedra::CalculateExtent()",
525  "GeomMgt1002",JustWarning,message);
526  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
527  }
528 
529  // set trigonometric values
530  G4double sphi = GetStartPhi();
531  G4double ephi = GetEndPhi();
532  G4double dphi = IsOpen() ? ephi-sphi : twopi;
533  G4int ksteps = GetNumSide();
534  G4double astep = dphi/ksteps;
535  G4double sinStep = std::sin(astep);
536  G4double cosStep = std::cos(astep);
537  G4double sinStart = GetSinStartPhi();
538  G4double cosStart = GetCosStartPhi();
539 
540  // allocate vector lists
541  std::vector<const G4ThreeVectorList *> polygons;
542  polygons.resize(ksteps+1);
543  for (G4int k=0; k<ksteps+1; ++k)
544  {
545  polygons[k] = new G4ThreeVectorList(3);
546  }
547 
548  // main loop along triangles
549  pMin = kInfinity;
550  pMax = -kInfinity;
551  G4int ntria = triangles.size()/3;
552  for (G4int i=0; i<ntria; ++i)
553  {
554  G4double sinCur = sinStart;
555  G4double cosCur = cosStart;
556  G4int i3 = i*3;
557  for (G4int k=0; k<ksteps+1; ++k) // rotate triangle
558  {
559  G4ThreeVectorList* ptr = const_cast<G4ThreeVectorList*>(polygons[k]);
560  G4ThreeVectorList::iterator iter = ptr->begin();
561  iter->set(triangles[i3+0].x()*cosCur,
562  triangles[i3+0].x()*sinCur,
563  triangles[i3+0].y());
564  iter++;
565  iter->set(triangles[i3+1].x()*cosCur,
566  triangles[i3+1].x()*sinCur,
567  triangles[i3+1].y());
568  iter++;
569  iter->set(triangles[i3+2].x()*cosCur,
570  triangles[i3+2].x()*sinCur,
571  triangles[i3+2].y());
572 
573  G4double sinTmp = sinCur;
574  sinCur = sinCur*cosStep + cosCur*sinStep;
575  cosCur = cosCur*cosStep - sinTmp*sinStep;
576  }
577 
578  // set sub-envelope and adjust extent
579  G4double emin,emax;
580  G4BoundingEnvelope benv(polygons);
581  if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
582  if (emin < pMin) pMin = emin;
583  if (emax > pMax) pMax = emax;
584  if (eminlim > pMin && emaxlim < pMax) break; // max possible extent
585  }
586  // free memory
587  for (G4int k=0; k<ksteps+1; ++k) { delete polygons[k]; polygons[k]=0;}
588  return (pMin < pMax);
589 }
590 
591 
593 //
594 // CreatePolyhedron
595 //
596 G4Polyhedron* G4UPolyhedra::CreatePolyhedron() const
597 {
598  if (!IsGeneric())
599  {
600  return new G4PolyhedronPgon( fOriginalParameters.Start_angle,
601  fOriginalParameters.Opening_angle,
602  fOriginalParameters.numSide,
603  fOriginalParameters.Num_z_planes,
604  fOriginalParameters.Z_values,
605  fOriginalParameters.Rmin,
606  fOriginalParameters.Rmax);
607  }
608  else
609  {
610  // The following code prepares for:
611  // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
612  // const double xyz[][3],
613  // const int faces_vec[][4])
614  // Here is an extract from the header file HepPolyhedron.h:
632  G4int nNodes;
633  G4int nFaces;
634  typedef G4double double3[3];
635  double3* xyz;
636  typedef G4int int4[4];
637  int4* faces_vec;
638  if (IsOpen())
639  {
640  // Triangulate open ends. Simple ear-chopping algorithm...
641  // I'm not sure how robust this algorithm is (J.Allison).
642  //
643  std::vector<G4bool> chopped(GetNumRZCorner(), false);
644  std::vector<G4int*> triQuads;
645  G4int remaining = GetNumRZCorner();
646  G4int iStarter = 0;
647  while (remaining >= 3) // Loop checking, 13.08.2015, G.Cosmo
648  {
649  // Find unchopped corners...
650  //
651  G4int A = -1, B = -1, C = -1;
652  G4int iStepper = iStarter;
653  do // Loop checking, 13.08.2015, G.Cosmo
654  {
655  if (A < 0) { A = iStepper; }
656  else if (B < 0) { B = iStepper; }
657  else if (C < 0) { C = iStepper; }
658  do // Loop checking, 13.08.2015, G.Cosmo
659  {
660  if (++iStepper >= GetNumRZCorner()) iStepper = 0;
661  }
662  while (chopped[iStepper]);
663  }
664  while (C < 0 && iStepper != iStarter);
665 
666  // Check triangle at B is pointing outward (an "ear").
667  // Sign of z cross product determines...
668 
669  G4double BAr = GetCorner(A).r - GetCorner(B).r;
670  G4double BAz = GetCorner(A).z - GetCorner(B).z;
671  G4double BCr = GetCorner(C).r - GetCorner(B).r;
672  G4double BCz = GetCorner(C).z - GetCorner(B).z;
673  if (BAr * BCz - BAz * BCr < kCarTolerance)
674  {
675  G4int* tq = new G4int[3];
676  tq[0] = A + 1;
677  tq[1] = B + 1;
678  tq[2] = C + 1;
679  triQuads.push_back(tq);
680  chopped[B] = true;
681  --remaining;
682  }
683  else
684  {
685  do // Loop checking, 13.08.2015, G.Cosmo
686  {
687  if (++iStarter >= GetNumRZCorner()) { iStarter = 0; }
688  }
689  while (chopped[iStarter]);
690  }
691  }
692 
693  // Transfer to faces...
694  G4int numSide=GetNumSide();
695  nNodes = (numSide + 1) * GetNumRZCorner();
696  nFaces = numSide * GetNumRZCorner() + 2 * triQuads.size();
697  faces_vec = new int4[nFaces];
698  G4int iface = 0;
699  G4int addition = GetNumRZCorner() * numSide;
700  G4int d = GetNumRZCorner() - 1;
701  for (G4int iEnd = 0; iEnd < 2; ++iEnd)
702  {
703  for (size_t i = 0; i < triQuads.size(); ++i)
704  {
705  // Negative for soft/auxiliary/normally invisible edges...
706  //
707  G4int a, b, c;
708  if (iEnd == 0)
709  {
710  a = triQuads[i][0];
711  b = triQuads[i][1];
712  c = triQuads[i][2];
713  }
714  else
715  {
716  a = triQuads[i][0] + addition;
717  b = triQuads[i][2] + addition;
718  c = triQuads[i][1] + addition;
719  }
720  G4int ab = std::abs(b - a);
721  G4int bc = std::abs(c - b);
722  G4int ca = std::abs(a - c);
723  faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
724  faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
725  faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
726  faces_vec[iface][3] = 0;
727  ++iface;
728  }
729  }
730 
731  // Continue with sides...
732 
733  xyz = new double3[nNodes];
734  const G4double dPhi = (GetEndPhi() - GetStartPhi()) / numSide;
735  G4double phi = GetStartPhi();
736  G4int ixyz = 0;
737  for (G4int iSide = 0; iSide < numSide; ++iSide)
738  {
739  for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
740  {
741  xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
742  xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
743  xyz[ixyz][2] = GetCorner(iCorner).z;
744  if (iCorner < GetNumRZCorner() - 1)
745  {
746  faces_vec[iface][0] = ixyz + 1;
747  faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
748  faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
749  faces_vec[iface][3] = ixyz + 2;
750  }
751  else
752  {
753  faces_vec[iface][0] = ixyz + 1;
754  faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
755  faces_vec[iface][2] = ixyz + 2;
756  faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
757  }
758  ++iface;
759  ++ixyz;
760  }
761  phi += dPhi;
762  }
763 
764  // Last GetCorner...
765 
766  for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
767  {
768  xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
769  xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
770  xyz[ixyz][2] = GetCorner(iCorner).z;
771  ++ixyz;
772  }
773  }
774  else // !phiIsOpen - i.e., a complete 360 degrees.
775  {
776  nNodes = GetNumSide() * GetNumRZCorner();
777  nFaces = GetNumSide() * GetNumRZCorner();;
778  xyz = new double3[nNodes];
779  faces_vec = new int4[nFaces];
780  // const G4double dPhi = (endPhi - startPhi) / numSide;
781  const G4double dPhi = twopi / GetNumSide();
782  // !phiIsOpen endPhi-startPhi = 360 degrees.
783  G4double phi = GetStartPhi();
784  G4int ixyz = 0, iface = 0;
785  for (G4int iSide = 0; iSide < GetNumSide(); ++iSide)
786  {
787  for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
788  {
789  xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
790  xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
791  xyz[ixyz][2] = GetCorner(iCorner).z;
792  if (iSide < GetNumSide() - 1)
793  {
794  if (iCorner < GetNumRZCorner() - 1)
795  {
796  faces_vec[iface][0] = ixyz + 1;
797  faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
798  faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
799  faces_vec[iface][3] = ixyz + 2;
800  }
801  else
802  {
803  faces_vec[iface][0] = ixyz + 1;
804  faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
805  faces_vec[iface][2] = ixyz + 2;
806  faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
807  }
808  }
809  else // Last side joins ends...
810  {
811  if (iCorner < GetNumRZCorner() - 1)
812  {
813  faces_vec[iface][0] = ixyz + 1;
814  faces_vec[iface][1] = ixyz + GetNumRZCorner() - nFaces + 1;
815  faces_vec[iface][2] = ixyz + GetNumRZCorner() - nFaces + 2;
816  faces_vec[iface][3] = ixyz + 2;
817  }
818  else
819  {
820  faces_vec[iface][0] = ixyz + 1;
821  faces_vec[iface][1] = ixyz - nFaces + GetNumRZCorner() + 1;
822  faces_vec[iface][2] = ixyz - nFaces + 2;
823  faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
824  }
825  }
826  ++ixyz;
827  ++iface;
828  }
829  phi += dPhi;
830  }
831  }
832  G4Polyhedron* polyhedron = new G4Polyhedron;
833  G4int prob = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
834  delete [] faces_vec;
835  delete [] xyz;
836  if (prob)
837  {
838  std::ostringstream message;
839  message << "Problem creating G4Polyhedron for: " << GetName();
840  G4Exception("G4Polyhedra::CreatePolyhedron()", "GeomSolids1002",
841  JustWarning, message);
842  delete polyhedron;
843  return nullptr;
844  }
845  else
846  {
847  return polyhedron;
848  }
849  }
850 }
851 
852 #endif // G4GEOM_USE_USOLIDS