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