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