Geant4  10.02.p02
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 "G4VPVParameterisation.hh"
38 
39 using CLHEP::twopi;
40 
42 //
43 // Constructor (GEANT3 style parameters)
44 //
45 // GEANT3 PGON radii are specified in the distance to the norm of each face.
46 //
47 G4UPolyhedra::G4UPolyhedra(const G4String& name,
48  G4double phiStart,
49  G4double phiTotal,
50  G4int numSide,
51  G4int numZPlanes,
52  const G4double zPlane[],
53  const G4double rInner[],
54  const G4double rOuter[] )
55  : G4USolid(name, new UPolyhedra(name,phiStart, phiTotal, numSide,
56  numZPlanes, zPlane, rInner, rOuter))
57 {
58 }
59 
60 
62 //
63 // Constructor (generic parameters)
64 //
65 G4UPolyhedra::G4UPolyhedra(const G4String& name,
66  G4double phiStart,
67  G4double phiTotal,
68  G4int numSide,
69  G4int numRZ,
70  const G4double r[],
71  const G4double z[] )
72  : G4USolid(name, new UPolyhedra(name, phiStart, phiTotal, numSide,
73  numRZ, r, z))
74 {
75 }
76 
77 
79 //
80 // Fake default constructor - sets only member data and allocates memory
81 // for usage restricted to object persistency.
82 //
83 G4UPolyhedra::G4UPolyhedra( __void__& a )
84  : G4USolid(a)
85 {
86 }
87 
88 
90 //
91 // Destructor
92 //
93 G4UPolyhedra::~G4UPolyhedra()
94 {
95 }
96 
97 
99 //
100 // Copy constructor
101 //
102 G4UPolyhedra::G4UPolyhedra( const G4UPolyhedra &source )
103  : G4USolid( source )
104 {
105 }
106 
107 
109 //
110 // Assignment operator
111 //
112 G4UPolyhedra& G4UPolyhedra::operator=( const G4UPolyhedra &source )
113 {
114  if (this == &source) return *this;
115 
116  G4USolid::operator=( source );
117 
118  return *this;
119 }
120 
121 
123 //
124 // Accessors & modifiers
125 //
126 G4int G4UPolyhedra::GetNumSide() const
127 {
128  return GetShape()->GetNumSide();
129 }
130 G4double G4UPolyhedra::GetStartPhi() const
131 {
132  return GetShape()->GetStartPhi();
133 }
134 G4double G4UPolyhedra::GetEndPhi() const
135 {
136  return GetShape()->GetEndPhi();
137 }
138 G4bool G4UPolyhedra::IsOpen() const
139 {
140  return GetShape()->IsOpen();
141 }
142 G4bool G4UPolyhedra::IsGeneric() const
143 {
144  return GetShape()->IsGeneric();
145 }
146 G4int G4UPolyhedra::GetNumRZCorner() const
147 {
148  return GetShape()->GetNumRZCorner();
149 }
150 G4PolyhedraSideRZ G4UPolyhedra::GetCorner(G4int index) const
151 {
152  UPolyhedraSideRZ pside = GetShape()->GetCorner(index);
153  G4PolyhedraSideRZ psiderz = { pside.r, pside.z };
154 
155  return psiderz;
156 }
157 G4PolyhedraHistorical* G4UPolyhedra::GetOriginalParameters() const
158 {
159  UPolyhedraHistorical* pars = GetShape()->GetOriginalParameters();
160  G4PolyhedraHistorical* pdata = new G4PolyhedraHistorical(pars->fNumZPlanes);
161  pdata->Start_angle = pars->fStartAngle;
162  pdata->Opening_angle = pars->fOpeningAngle;
163  pdata->numSide = pars->fNumSide;
164  for (G4int i=0; i<pars->fNumZPlanes; ++i)
165  {
166  pdata->Z_values[i] = pars->fZValues[i];
167  pdata->Rmin[i] = pars->Rmin[i];
168  pdata->Rmax[i] = pars->Rmax[i];
169  }
170  return pdata;
171 }
172 void G4UPolyhedra::SetOriginalParameters(G4PolyhedraHistorical* pars)
173 {
174  UPolyhedraHistorical* pdata = GetShape()->GetOriginalParameters();
175  pdata->fStartAngle = pars->Start_angle;
176  pdata->fOpeningAngle = pars->Opening_angle;
177  pdata->fNumSide = pars->numSide;
178  pdata->fNumZPlanes = pars->Num_z_planes;
179  for (G4int i=0; i<pdata->fNumZPlanes; ++i)
180  {
181  pdata->fZValues[i] = pars->Z_values[i];
182  pdata->Rmin[i] = pars->Rmin[i];
183  pdata->Rmax[i] = pars->Rmax[i];
184  }
185  fRebuildPolyhedron = true;
186 }
187 G4bool G4UPolyhedra::Reset()
188 {
189  return GetShape()->Reset();
190 }
191 
192 
194 //
195 // Dispatch to parameterisation for replication mechanism dimension
196 // computation & modification.
197 //
198 void G4UPolyhedra::ComputeDimensions(G4VPVParameterisation* p,
199  const G4int n,
200  const G4VPhysicalVolume* pRep)
201 {
202  p->ComputeDimensions(*(G4Polyhedra*)this,n,pRep);
203 }
204 
205 
207 //
208 // Make a clone of the object
209 
210 G4VSolid* G4UPolyhedra::Clone() const
211 {
212  return new G4UPolyhedra(*this);
213 }
214 
215 
217 //
218 // CreatePolyhedron
219 //
220 G4Polyhedron* G4UPolyhedra::CreatePolyhedron() const
221 {
222  if (!IsGeneric())
223  {
224  G4PolyhedraHistorical* original_parameters = GetOriginalParameters();
226  polyhedron = new G4PolyhedronPgon( GetOriginalParameters()->Start_angle,
227  GetOriginalParameters()->Opening_angle,
228  GetOriginalParameters()->numSide,
229  GetOriginalParameters()->Num_z_planes,
230  GetOriginalParameters()->Z_values,
231  GetOriginalParameters()->Rmin,
232  GetOriginalParameters()->Rmax);
233  delete original_parameters; // delete local copy
234  return polyhedron;
235  }
236  else
237  {
238  // The following code prepares for:
239  // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
240  // const double xyz[][3],
241  // const int faces_vec[][4])
242  // Here is an extract from the header file HepPolyhedron.h:
260  G4int nNodes;
261  G4int nFaces;
262  typedef G4double double3[3];
263  double3* xyz;
264  typedef G4int int4[4];
265  int4* faces_vec;
266  if (IsOpen())
267  {
268  // Triangulate open ends. Simple ear-chopping algorithm...
269  // I'm not sure how robust this algorithm is (J.Allison).
270  //
271  std::vector<G4bool> chopped(GetNumRZCorner(), false);
272  std::vector<G4int*> triQuads;
273  G4int remaining = GetNumRZCorner();
274  G4int iStarter = 0;
275  while (remaining >= 3) // Loop checking, 13.08.2015, G.Cosmo
276  {
277  // Find unchopped corners...
278  //
279  G4int A = -1, B = -1, C = -1;
280  G4int iStepper = iStarter;
281  do // Loop checking, 13.08.2015, G.Cosmo
282  {
283  if (A < 0) { A = iStepper; }
284  else if (B < 0) { B = iStepper; }
285  else if (C < 0) { C = iStepper; }
286  do // Loop checking, 13.08.2015, G.Cosmo
287  {
288  if (++iStepper >= GetNumRZCorner()) iStepper = 0;
289  }
290  while (chopped[iStepper]);
291  }
292  while (C < 0 && iStepper != iStarter);
293 
294  // Check triangle at B is pointing outward (an "ear").
295  // Sign of z cross product determines...
296 
297  G4double BAr = GetCorner(A).r - GetCorner(B).r;
298  G4double BAz = GetCorner(A).z - GetCorner(B).z;
299  G4double BCr = GetCorner(C).r - GetCorner(B).r;
300  G4double BCz = GetCorner(C).z - GetCorner(B).z;
301  if (BAr * BCz - BAz * BCr < kCarTolerance)
302  {
303  G4int* tq = new G4int[3];
304  tq[0] = A + 1;
305  tq[1] = B + 1;
306  tq[2] = C + 1;
307  triQuads.push_back(tq);
308  chopped[B] = true;
309  --remaining;
310  }
311  else
312  {
313  do // Loop checking, 13.08.2015, G.Cosmo
314  {
315  if (++iStarter >= GetNumRZCorner()) { iStarter = 0; }
316  }
317  while (chopped[iStarter]);
318  }
319  }
320 
321  // Transfer to faces...
322  G4int numSide=GetNumSide();
323  nNodes = (numSide + 1) * GetNumRZCorner();
324  nFaces = numSide * GetNumRZCorner() + 2 * triQuads.size();
325  faces_vec = new int4[nFaces];
326  G4int iface = 0;
327  G4int addition = GetNumRZCorner() * numSide;
328  G4int d = GetNumRZCorner() - 1;
329  for (G4int iEnd = 0; iEnd < 2; ++iEnd)
330  {
331  for (size_t i = 0; i < triQuads.size(); ++i)
332  {
333  // Negative for soft/auxiliary/normally invisible edges...
334  //
335  G4int a, b, c;
336  if (iEnd == 0)
337  {
338  a = triQuads[i][0];
339  b = triQuads[i][1];
340  c = triQuads[i][2];
341  }
342  else
343  {
344  a = triQuads[i][0] + addition;
345  b = triQuads[i][2] + addition;
346  c = triQuads[i][1] + addition;
347  }
348  G4int ab = std::abs(b - a);
349  G4int bc = std::abs(c - b);
350  G4int ca = std::abs(a - c);
351  faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
352  faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
353  faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
354  faces_vec[iface][3] = 0;
355  ++iface;
356  }
357  }
358 
359  // Continue with sides...
360 
361  xyz = new double3[nNodes];
362  const G4double dPhi = (GetEndPhi() - GetStartPhi()) / numSide;
363  G4double phi = GetStartPhi();
364  G4int ixyz = 0;
365  for (G4int iSide = 0; iSide < numSide; ++iSide)
366  {
367  for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
368  {
369  xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
370  xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
371  xyz[ixyz][2] = GetCorner(iCorner).z;
372  if (iCorner < GetNumRZCorner() - 1)
373  {
374  faces_vec[iface][0] = ixyz + 1;
375  faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
376  faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
377  faces_vec[iface][3] = ixyz + 2;
378  }
379  else
380  {
381  faces_vec[iface][0] = ixyz + 1;
382  faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
383  faces_vec[iface][2] = ixyz + 2;
384  faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
385  }
386  ++iface;
387  ++ixyz;
388  }
389  phi += dPhi;
390  }
391 
392  // Last GetCorner...
393 
394  for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
395  {
396  xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
397  xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
398  xyz[ixyz][2] = GetCorner(iCorner).z;
399  ++ixyz;
400  }
401  }
402  else // !phiIsOpen - i.e., a complete 360 degrees.
403  {
404  nNodes = GetNumSide() * GetNumRZCorner();
405  nFaces = GetNumSide() * GetNumRZCorner();;
406  xyz = new double3[nNodes];
407  faces_vec = new int4[nFaces];
408  // const G4double dPhi = (endPhi - startPhi) / numSide;
409  const G4double dPhi = twopi / GetNumSide(); // !phiIsOpen endPhi-startPhi = 360 degrees.
410  G4double phi = GetStartPhi();
411  G4int ixyz = 0, iface = 0;
412  for (G4int iSide = 0; iSide < GetNumSide(); ++iSide)
413  {
414  for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
415  {
416  xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
417  xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
418  xyz[ixyz][2] = GetCorner(iCorner).z;
419  if (iSide < GetNumSide() - 1)
420  {
421  if (iCorner < GetNumRZCorner() - 1)
422  {
423  faces_vec[iface][0] = ixyz + 1;
424  faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
425  faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
426  faces_vec[iface][3] = ixyz + 2;
427  }
428  else
429  {
430  faces_vec[iface][0] = ixyz + 1;
431  faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
432  faces_vec[iface][2] = ixyz + 2;
433  faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
434  }
435  }
436  else // Last side joins ends...
437  {
438  if (iCorner < GetNumRZCorner() - 1)
439  {
440  faces_vec[iface][0] = ixyz + 1;
441  faces_vec[iface][1] = ixyz + GetNumRZCorner() - nFaces + 1;
442  faces_vec[iface][2] = ixyz + GetNumRZCorner() - nFaces + 2;
443  faces_vec[iface][3] = ixyz + 2;
444  }
445  else
446  {
447  faces_vec[iface][0] = ixyz + 1;
448  faces_vec[iface][1] = ixyz - nFaces + GetNumRZCorner() + 1;
449  faces_vec[iface][2] = ixyz - nFaces + 2;
450  faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
451  }
452  }
453  ++ixyz;
454  ++iface;
455  }
456  phi += dPhi;
457  }
458  }
459  G4Polyhedron* polyhedron = new G4Polyhedron;
460  G4int problem = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
461  delete [] faces_vec;
462  delete [] xyz;
463  if (problem)
464  {
465  std::ostringstream message;
466  message << "Problem creating G4Polyhedron for: " << GetName();
467  G4Exception("G4Polyhedra::CreatePolyhedron()", "GeomSolids1002",
468  JustWarning, message);
469  delete polyhedron;
470  return 0;
471  }
472  else
473  {
474  return polyhedron;
475  }
476  }
477 }
478 
479 #endif // G4GEOM_USE_USOLIDS
G4double z
Definition: TRTMaterials.hh:39
G4String name
Definition: TRTMaterials.hh:40
G4double a
Definition: TRTMaterials.hh:39
double C(double temp)
double B(double temperature)
int G4int
Definition: G4Types.hh:78
double A(double temperature)
bool G4bool
Definition: G4Types.hh:79
static const double twopi
Definition: G4SIunits.hh:75
const G4int n
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
static const G4double ab
double G4double
Definition: G4Types.hh:76