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