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