Geant4  10.00.p01
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 #include "G4Polyhedron.hh"
35 #include "G4VPVParameterisation.hh"
36 
37 using CLHEP::twopi;
38 
40 //
41 // Constructor (GEANT3 style parameters)
42 //
43 // GEANT3 PGON radii are specified in the distance to the norm of each face.
44 //
46  G4double phiStart,
47  G4double phiTotal,
48  G4int numSide,
49  G4int numZPlanes,
50  const G4double zPlane[],
51  const G4double rInner[],
52  const G4double rOuter[] )
53  : G4USolid(name, new UPolyhedra(name,phiStart, phiTotal, numSide,
54  numZPlanes, zPlane, rInner, rOuter))
55 {
56 }
57 
58 
60 //
61 // Constructor (generic parameters)
62 //
64  G4double phiStart,
65  G4double phiTotal,
66  G4int numSide,
67  G4int numRZ,
68  const G4double r[],
69  const G4double z[] )
70  : G4USolid(name, new UPolyhedra(name, phiStart, phiTotal, numSide,
71  numRZ, r, z))
72 {
73 }
74 
75 
77 //
78 // Fake default constructor - sets only member data and allocates memory
79 // for usage restricted to object persistency.
80 //
82  : G4USolid(a)
83 {
84 }
85 
86 
88 //
89 // Destructor
90 //
92 {
93 }
94 
95 
97 //
98 // Copy constructor
99 //
101  : G4USolid( source )
102 {
103 }
104 
105 
107 //
108 // Assignment operator
109 //
111 {
112  if (this == &source) return *this;
113 
114  G4USolid::operator=( source );
115 
116  return *this;
117 }
118 
119 
121 //
122 // Dispatch to parameterisation for replication mechanism dimension
123 // computation & modification.
124 //
126  const G4int n,
127  const G4VPhysicalVolume* pRep)
128 {
129  p->ComputeDimensions(*(G4Polyhedra*)this,n,pRep);
130 }
131 
132 
134 //
135 // Make a clone of the object
136 
138 {
139  return new G4UPolyhedra(*this);
140 }
141 
142 
144 //
145 // CreatePolyhedron
146 //
148 {
149  if (!IsGeneric())
150  {
151  G4PolyhedraHistorical* original_parameters = GetOriginalParameters();
153  polyhedron = new G4PolyhedronPgon( GetOriginalParameters()->Start_angle,
154  GetOriginalParameters()->Opening_angle,
155  GetOriginalParameters()->numSide,
156  GetOriginalParameters()->Num_z_planes,
157  GetOriginalParameters()->Z_values,
158  GetOriginalParameters()->Rmin,
159  GetOriginalParameters()->Rmax);
160  delete original_parameters; // delete local copy
161  return polyhedron;
162  }
163  else
164  {
165  // The following code prepares for:
166  // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
167  // const double xyz[][3],
168  // const int faces_vec[][4])
169  // Here is an extract from the header file HepPolyhedron.h:
187  G4int nNodes;
188  G4int nFaces;
189  typedef G4double double3[3];
190  double3* xyz;
191  typedef G4int int4[4];
192  int4* faces_vec;
193  if (IsOpen())
194  {
195  // Triangulate open ends. Simple ear-chopping algorithm...
196  // I'm not sure how robust this algorithm is (J.Allison).
197  //
198  std::vector<G4bool> chopped(GetNumRZCorner(), false);
199  std::vector<G4int*> triQuads;
200  G4int remaining = GetNumRZCorner();
201  G4int iStarter = 0;
202  while (remaining >= 3)
203  {
204  // Find unchopped corners...
205  //
206  G4int A = -1, B = -1, C = -1;
207  G4int iStepper = iStarter;
208  do
209  {
210  if (A < 0) { A = iStepper; }
211  else if (B < 0) { B = iStepper; }
212  else if (C < 0) { C = iStepper; }
213  do
214  {
215  if (++iStepper >= GetNumRZCorner()) iStepper = 0;
216  }
217  while (chopped[iStepper]);
218  }
219  while (C < 0 && iStepper != iStarter);
220 
221  // Check triangle at B is pointing outward (an "ear").
222  // Sign of z cross product determines...
223 
224  G4double BAr = GetCorner(A).r - GetCorner(B).r;
225  G4double BAz = GetCorner(A).z - GetCorner(B).z;
226  G4double BCr = GetCorner(C).r - GetCorner(B).r;
227  G4double BCz = GetCorner(C).z - GetCorner(B).z;
228  if (BAr * BCz - BAz * BCr < kCarTolerance)
229  {
230  G4int* tq = new G4int[3];
231  tq[0] = A + 1;
232  tq[1] = B + 1;
233  tq[2] = C + 1;
234  triQuads.push_back(tq);
235  chopped[B] = true;
236  --remaining;
237  }
238  else
239  {
240  do
241  {
242  if (++iStarter >= GetNumRZCorner()) { iStarter = 0; }
243  }
244  while (chopped[iStarter]);
245  }
246  }
247 
248  // Transfer to faces...
249  G4int numSide=GetNumSide();
250  nNodes = (numSide + 1) * GetNumRZCorner();
251  nFaces = numSide * GetNumRZCorner() + 2 * triQuads.size();
252  faces_vec = new int4[nFaces];
253  G4int iface = 0;
254  G4int addition = GetNumRZCorner() * numSide;
255  G4int d = GetNumRZCorner() - 1;
256  for (G4int iEnd = 0; iEnd < 2; ++iEnd)
257  {
258  for (size_t i = 0; i < triQuads.size(); ++i)
259  {
260  // Negative for soft/auxiliary/normally invisible edges...
261  //
262  G4int a, b, c;
263  if (iEnd == 0)
264  {
265  a = triQuads[i][0];
266  b = triQuads[i][1];
267  c = triQuads[i][2];
268  }
269  else
270  {
271  a = triQuads[i][0] + addition;
272  b = triQuads[i][2] + addition;
273  c = triQuads[i][1] + addition;
274  }
275  G4int ab = std::abs(b - a);
276  G4int bc = std::abs(c - b);
277  G4int ca = std::abs(a - c);
278  faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
279  faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
280  faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
281  faces_vec[iface][3] = 0;
282  ++iface;
283  }
284  }
285 
286  // Continue with sides...
287 
288  xyz = new double3[nNodes];
289  const G4double dPhi = (GetEndPhi() - GetStartPhi()) / numSide;
290  G4double phi = GetStartPhi();
291  G4int ixyz = 0;
292  for (G4int iSide = 0; iSide < numSide; ++iSide)
293  {
294  for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
295  {
296  xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
297  xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
298  xyz[ixyz][2] = GetCorner(iCorner).z;
299  if (iCorner < GetNumRZCorner() - 1)
300  {
301  faces_vec[iface][0] = ixyz + 1;
302  faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
303  faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
304  faces_vec[iface][3] = ixyz + 2;
305  }
306  else
307  {
308  faces_vec[iface][0] = ixyz + 1;
309  faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
310  faces_vec[iface][2] = ixyz + 2;
311  faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
312  }
313  ++iface;
314  ++ixyz;
315  }
316  phi += dPhi;
317  }
318 
319  // Last GetCorner...
320 
321  for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
322  {
323  xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
324  xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
325  xyz[ixyz][2] = GetCorner(iCorner).z;
326  ++ixyz;
327  }
328  }
329  else // !phiIsOpen - i.e., a complete 360 degrees.
330  {
331  nNodes = GetNumSide() * GetNumRZCorner();
332  nFaces = GetNumSide() * GetNumRZCorner();;
333  xyz = new double3[nNodes];
334  faces_vec = new int4[nFaces];
335  // const G4double dPhi = (endPhi - startPhi) / numSide;
336  const G4double dPhi = twopi / GetNumSide(); // !phiIsOpen endPhi-startPhi = 360 degrees.
337  G4double phi = GetStartPhi();
338  G4int ixyz = 0, iface = 0;
339  for (G4int iSide = 0; iSide < GetNumSide(); ++iSide)
340  {
341  for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
342  {
343  xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
344  xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
345  xyz[ixyz][2] = GetCorner(iCorner).z;
346  if (iSide < GetNumSide() - 1)
347  {
348  if (iCorner < GetNumRZCorner() - 1)
349  {
350  faces_vec[iface][0] = ixyz + 1;
351  faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
352  faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
353  faces_vec[iface][3] = ixyz + 2;
354  }
355  else
356  {
357  faces_vec[iface][0] = ixyz + 1;
358  faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
359  faces_vec[iface][2] = ixyz + 2;
360  faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
361  }
362  }
363  else // Last side joins ends...
364  {
365  if (iCorner < GetNumRZCorner() - 1)
366  {
367  faces_vec[iface][0] = ixyz + 1;
368  faces_vec[iface][1] = ixyz + GetNumRZCorner() - nFaces + 1;
369  faces_vec[iface][2] = ixyz + GetNumRZCorner() - nFaces + 2;
370  faces_vec[iface][3] = ixyz + 2;
371  }
372  else
373  {
374  faces_vec[iface][0] = ixyz + 1;
375  faces_vec[iface][1] = ixyz - nFaces + GetNumRZCorner() + 1;
376  faces_vec[iface][2] = ixyz - nFaces + 2;
377  faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
378  }
379  }
380  ++ixyz;
381  ++iface;
382  }
383  phi += dPhi;
384  }
385  }
386  G4Polyhedron* polyhedron = new G4Polyhedron;
387  G4int problem = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
388  delete [] faces_vec;
389  delete [] xyz;
390  if (problem)
391  {
392  std::ostringstream message;
393  message << "Problem creating G4Polyhedron for: " << GetName();
394  G4Exception("G4Polyhedra::CreatePolyhedron()", "GeomSolids1002",
395  JustWarning, message);
396  delete polyhedron;
397  return 0;
398  }
399  else
400  {
401  return polyhedron;
402  }
403  }
404 }
G4String GetName() const
G4UPolyhedra & operator=(const G4UPolyhedra &source)
G4double z
Definition: TRTMaterials.hh:39
G4String name
Definition: TRTMaterials.hh:40
G4double a
Definition: TRTMaterials.hh:39
int G4int
Definition: G4Types.hh:78
G4int GetNumRZCorner() const
G4UPolyhedra(const G4String &name, G4double phiStart, G4double phiTotal, G4int numSide, G4int numZPlanes, const G4double zPlane[], const G4double rInner[], const G4double rOuter[])
Definition: G4UPolyhedra.cc:45
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
G4double GetEndPhi() const
G4double GetStartPhi() const
const G4int n
G4PolyhedraHistorical * GetOriginalParameters() const
static const G4double A[nN]
G4bool IsGeneric() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4PolyhedraSideRZ GetCorner(const G4int index) const
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
G4VSolid * Clone() const
static const G4double ab
G4Polyhedron * CreatePolyhedron() const
G4int GetNumSide() const
G4USolid & operator=(const G4USolid &rhs)
Definition: G4USolid.cc:370
G4double kCarTolerance
Definition: G4VSolid.hh:305
double G4double
Definition: G4Types.hh:76
G4bool IsOpen() const