Geant4  10.00.p02
G4UGenericPolycone.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 //
30 // Implementation of G4UGenericPolycone wrapper class
31 // --------------------------------------------------------------------
32 
33 #include "G4GenericPolycone.hh"
34 #include "G4UGenericPolycone.hh"
35 #include "G4Polyhedron.hh"
36 
37 using CLHEP::twopi;
38 
40 //
41 // Constructor (generic parameters)
42 //
44  G4double phiStart,
45  G4double phiTotal,
46  G4int numRZ,
47  const G4double r[],
48  const G4double z[] )
49  : G4USolid(name, new UGenericPolycone(name, phiStart, phiTotal, numRZ, r, z))
50 {
51 }
52 
53 
55 //
56 // Fake default constructor - sets only member data and allocates memory
57 // for usage restricted to object persistency.
58 //
60  : G4USolid(a)
61 {
62 }
63 
64 
66 //
67 // Destructor
68 //
70 {
71 }
72 
73 
75 //
76 // Copy constructor
77 //
79  : G4USolid(source)
80 {
81 }
82 
83 
85 //
86 // Assignment operator
87 //
90 {
91  if (this == &source) return *this;
92 
93  G4USolid::operator=( source );
94 
95  return *this;
96 }
97 
99 {
100 
101 
102  // The following code prepares for:
103  // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
104  // const double xyz[][3],
105  // const int faces_vec[][4])
106  // Here is an extract from the header file HepPolyhedron.h:
124  const G4int numSide =
125  G4int(G4Polyhedron::GetNumberOfRotationSteps()
126  * (GetEndPhi() - GetStartPhi()) / twopi) + 1;
127  G4int nNodes;
128  G4int nFaces;
129  typedef G4double double3[3];
130  double3* xyz;
131  typedef G4int int4[4];
132  int4* faces_vec;
133  if (IsOpen())
134  {
135  // Triangulate open ends. Simple ear-chopping algorithm...
136  // I'm not sure how robust this algorithm is (J.Allison).
137  //
138  std::vector<G4bool> chopped(GetNumRZCorner(), false);
139  std::vector<G4int*> triQuads;
140  G4int remaining = GetNumRZCorner();
141  G4int iStarter = 0;
142  while (remaining >= 3)
143  {
144  // Find unchopped corners...
145  //
146  G4int A = -1, B = -1, C = -1;
147  G4int iStepper = iStarter;
148  do
149  {
150  if (A < 0) { A = iStepper; }
151  else if (B < 0) { B = iStepper; }
152  else if (C < 0) { C = iStepper; }
153  do
154  {
155  if (++iStepper >= GetNumRZCorner()) { iStepper = 0; }
156  }
157  while (chopped[iStepper]);
158  }
159  while (C < 0 && iStepper != iStarter);
160 
161  // Check triangle at B is pointing outward (an "ear").
162  // Sign of z cross product determines...
163  //
164  G4double BAr = GetCorner(A).r - GetCorner(B).r;
165  G4double BAz = GetCorner(A).z - GetCorner(B).z;
166  G4double BCr = GetCorner(C).r - GetCorner(B).r;
167  G4double BCz = GetCorner(C).z - GetCorner(B).z;
168  if (BAr * BCz - BAz * BCr < kCarTolerance)
169  {
170  G4int* tq = new G4int[3];
171  tq[0] = A + 1;
172  tq[1] = B + 1;
173  tq[2] = C + 1;
174  triQuads.push_back(tq);
175  chopped[B] = true;
176  --remaining;
177  }
178  else
179  {
180  do
181  {
182  if (++iStarter >= GetNumRZCorner()) { iStarter = 0; }
183  }
184  while (chopped[iStarter]);
185  }
186  }
187  // Transfer to faces...
188  //
189  nNodes = (numSide + 1) * GetNumRZCorner();
190  nFaces = numSide * GetNumRZCorner() + 2 * triQuads.size();
191  faces_vec = new int4[nFaces];
192  G4int iface = 0;
193  G4int addition = GetNumRZCorner() * numSide;
194  G4int d = GetNumRZCorner() - 1;
195  for (G4int iEnd = 0; iEnd < 2; ++iEnd)
196  {
197  for (size_t i = 0; i < triQuads.size(); ++i)
198  {
199  // Negative for soft/auxiliary/normally invisible edges...
200  //
201  G4int a, b, c;
202  if (iEnd == 0)
203  {
204  a = triQuads[i][0];
205  b = triQuads[i][1];
206  c = triQuads[i][2];
207  }
208  else
209  {
210  a = triQuads[i][0] + addition;
211  b = triQuads[i][2] + addition;
212  c = triQuads[i][1] + addition;
213  }
214  G4int ab = std::abs(b - a);
215  G4int bc = std::abs(c - b);
216  G4int ca = std::abs(a - c);
217  faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
218  faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
219  faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
220  faces_vec[iface][3] = 0;
221  ++iface;
222  }
223  }
224 
225  // Continue with sides...
226 
227  xyz = new double3[nNodes];
228  const G4double dPhi = (GetEndPhi() - GetStartPhi()) / numSide;
229  G4double phi = GetStartPhi();
230  G4int ixyz = 0;
231  for (G4int iSide = 0; iSide < numSide; ++iSide)
232  {
233  for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
234  {
235  xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
236  xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
237  xyz[ixyz][2] = GetCorner(iCorner).z;
238  if (iSide == 0) // startPhi
239  {
240  if (iCorner < GetNumRZCorner() - 1)
241  {
242  faces_vec[iface][0] = ixyz + 1;
243  faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
244  faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
245  faces_vec[iface][3] = ixyz + 2;
246  }
247  else
248  {
249  faces_vec[iface][0] = ixyz + 1;
250  faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
251  faces_vec[iface][2] = ixyz + 2;
252  faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
253  }
254  }
255  else if (iSide == numSide - 1) // endPhi
256  {
257  if (iCorner < GetNumRZCorner() - 1)
258  {
259  faces_vec[iface][0] = ixyz + 1;
260  faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
261  faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
262  faces_vec[iface][3] = -(ixyz + 2);
263  }
264  else
265  {
266  faces_vec[iface][0] = ixyz + 1;
267  faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
268  faces_vec[iface][2] = ixyz + 2;
269  faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
270  }
271  }
272  else
273  {
274  if (iCorner < GetNumRZCorner() - 1)
275  {
276  faces_vec[iface][0] = ixyz + 1;
277  faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
278  faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
279  faces_vec[iface][3] = -(ixyz + 2);
280  }
281  else
282  {
283  faces_vec[iface][0] = ixyz + 1;
284  faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
285  faces_vec[iface][2] = ixyz + 2;
286  faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
287  }
288  }
289  ++iface;
290  ++ixyz;
291  }
292  phi += dPhi;
293  }
294 
295  // Last corners...
296 
297  for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
298  {
299  xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
300  xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
301  xyz[ixyz][2] = GetCorner(iCorner).z;
302  ++ixyz;
303  }
304  }
305  else // !phiIsOpen - i.e., a complete 360 degrees.
306  {
307  nNodes = numSide * GetNumRZCorner();
308  nFaces = numSide * GetNumRZCorner();;
309  xyz = new double3[nNodes];
310  faces_vec = new int4[nFaces];
311  const G4double dPhi = (GetEndPhi() - GetStartPhi()) / numSide;
312  G4double phi = GetStartPhi();
313  G4int ixyz = 0, iface = 0;
314  for (G4int iSide = 0; iSide < numSide; ++iSide)
315  {
316  for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
317  {
318  xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
319  xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
320  xyz[ixyz][2] = GetCorner(iCorner).z;
321 
322  if (iSide < numSide - 1)
323  {
324  if (iCorner < GetNumRZCorner() - 1)
325  {
326  faces_vec[iface][0] = ixyz + 1;
327  faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
328  faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
329  faces_vec[iface][3] = -(ixyz + 2);
330  }
331  else
332  {
333  faces_vec[iface][0] = ixyz + 1;
334  faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
335  faces_vec[iface][2] = ixyz + 2;
336  faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
337  }
338  }
339  else // Last side joins ends...
340  {
341  if (iCorner < GetNumRZCorner() - 1)
342  {
343  faces_vec[iface][0] = ixyz + 1;
344  faces_vec[iface][1] = -(ixyz + GetNumRZCorner() - nFaces + 1);
345  faces_vec[iface][2] = ixyz + GetNumRZCorner() - nFaces + 2;
346  faces_vec[iface][3] = -(ixyz + 2);
347  }
348  else
349  {
350  faces_vec[iface][0] = ixyz + 1;
351  faces_vec[iface][1] = -(ixyz - nFaces + GetNumRZCorner() + 1);
352  faces_vec[iface][2] = ixyz - nFaces + 2;
353  faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
354  }
355  }
356  ++ixyz;
357  ++iface;
358  }
359  phi += dPhi;
360  }
361  }
362  G4Polyhedron* polyhedron = new G4Polyhedron;
363  G4int problem = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
364  delete [] faces_vec;
365  delete [] xyz;
366  if (problem)
367  {
368  std::ostringstream message;
369  message << "Problem creating G4Polyhedron for: " << GetName();
370  G4Exception("G4GenericPolycone::CreatePolyhedron()", "GeomSolids1002",
371  JustWarning, message);
372  delete polyhedron;
373  return 0;
374  }
375  else
376  {
377  return polyhedron;
378  }
379 }
380 
G4double GetStartPhi() const
G4String GetName() const
G4double z
Definition: TRTMaterials.hh:39
G4bool IsOpen() const
G4String name
Definition: TRTMaterials.hh:40
G4double a
Definition: TRTMaterials.hh:39
G4Polyhedron * CreatePolyhedron() const
int G4int
Definition: G4Types.hh:78
G4PolyconeSideRZ GetCorner(G4int index) const
G4double GetEndPhi() const
G4int GetNumRZCorner() const
G4UGenericPolycone(const G4String &name, G4double phiStart, G4double phiTotal, G4int numRZ, const G4double r[], const G4double z[])
static const G4double A[nN]
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4UGenericPolycone & operator=(const G4UGenericPolycone &source)
static const G4double ab
G4USolid & operator=(const G4USolid &rhs)
Definition: G4USolid.cc:372
G4double kCarTolerance
Definition: G4VSolid.hh:305
double G4double
Definition: G4Types.hh:76