| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // Implementation of G4UGenericPolycone wrappe <<
27 // 26 //
28 // 30.10.13 G.Cosmo, CERN << 27 // $Id:$
>> 28 //
>> 29 //
>> 30 // Implementation of G4UGenericPolycone wrapper class
29 // ------------------------------------------- 31 // --------------------------------------------------------------------
30 32
31 #include "G4GenericPolycone.hh" 33 #include "G4GenericPolycone.hh"
32 #include "G4UGenericPolycone.hh" 34 #include "G4UGenericPolycone.hh"
33 35
34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G 36 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G4GEOM_USE_PARTIAL_USOLIDS) )
35 37
36 #include "G4GeomTools.hh" <<
37 #include "G4AffineTransform.hh" <<
38 #include "G4VPVParameterisation.hh" <<
39 #include "G4BoundingEnvelope.hh" <<
40 <<
41 #include "G4Polyhedron.hh" 38 #include "G4Polyhedron.hh"
42 39
43 using namespace CLHEP; << 40 using CLHEP::twopi;
44 41
45 ////////////////////////////////////////////// 42 ////////////////////////////////////////////////////////////////////////
46 // 43 //
47 // Constructor (generic parameters) 44 // Constructor (generic parameters)
48 // 45 //
49 G4UGenericPolycone::G4UGenericPolycone(const G << 46 G4UGenericPolycone::G4UGenericPolycone(const G4String& name,
50 G 47 G4double phiStart,
51 G 48 G4double phiTotal,
52 G 49 G4int numRZ,
53 const G 50 const G4double r[],
54 const G 51 const G4double z[] )
55 : Base_t(name, phiStart, phiTotal, numRZ, r, << 52 : G4USolid(name, new UGenericPolycone(name, phiStart, phiTotal, numRZ, r, z))
56 { << 53 {
57 wrStart = phiStart; while (wrStart < 0) wrSt <<
58 wrDelta = phiTotal; <<
59 if (wrDelta <= 0 || wrDelta >= twopi*(1-DBL_ <<
60 { <<
61 wrStart = 0; <<
62 wrDelta = twopi; <<
63 } <<
64 rzcorners.resize(0); <<
65 for (G4int i=0; i<numRZ; ++i) <<
66 { <<
67 rzcorners.emplace_back(r[i],z[i]); <<
68 } <<
69 std::vector<G4int> iout; <<
70 G4GeomTools::RemoveRedundantVertices(rzcorne <<
71 } 54 }
72 55
73 56
74 ////////////////////////////////////////////// 57 ////////////////////////////////////////////////////////////////////////
75 // 58 //
76 // Fake default constructor - sets only member 59 // Fake default constructor - sets only member data and allocates memory
77 // for usage restri 60 // for usage restricted to object persistency.
78 // 61 //
79 G4UGenericPolycone::G4UGenericPolycone(__void_ 62 G4UGenericPolycone::G4UGenericPolycone(__void__& a)
80 : Base_t(a) << 63 : G4USolid(a)
81 { 64 {
82 } 65 }
83 66
84 67
85 ////////////////////////////////////////////// 68 //////////////////////////////////////////////////////////////////////////
86 // 69 //
87 // Destructor 70 // Destructor
88 // 71 //
89 G4UGenericPolycone::~G4UGenericPolycone() = de << 72 G4UGenericPolycone::~G4UGenericPolycone()
>> 73 {
>> 74 }
90 75
91 76
92 ////////////////////////////////////////////// 77 //////////////////////////////////////////////////////////////////////////
93 // 78 //
94 // Copy constructor 79 // Copy constructor
95 // 80 //
96 G4UGenericPolycone::G4UGenericPolycone(const G << 81 G4UGenericPolycone::G4UGenericPolycone(const G4UGenericPolycone &source)
97 : Base_t(source) << 82 : G4USolid(source)
98 { 83 {
99 wrStart = source.wrStart; <<
100 wrDelta = source.wrDelta; <<
101 rzcorners = source.rzcorners; <<
102 } 84 }
103 85
104 86
105 ////////////////////////////////////////////// 87 //////////////////////////////////////////////////////////////////////////
106 // 88 //
107 // Assignment operator 89 // Assignment operator
108 // 90 //
109 G4UGenericPolycone& 91 G4UGenericPolycone&
110 G4UGenericPolycone::operator=(const G4UGeneric << 92 G4UGenericPolycone::operator=(const G4UGenericPolycone &source)
111 { 93 {
112 if (this == &source) return *this; 94 if (this == &source) return *this;
113 << 95
114 Base_t::operator=( source ); << 96 G4USolid::operator=( source );
115 wrStart = source.wrStart; << 97
116 wrDelta = source.wrDelta; <<
117 rzcorners = source.rzcorners; <<
118 <<
119 return *this; 98 return *this;
120 } 99 }
121 100
122 G4double G4UGenericPolycone::GetStartPhi() con 101 G4double G4UGenericPolycone::GetStartPhi() const
123 { 102 {
124 return wrStart; << 103 return GetShape()->GetStartPhi();
125 } 104 }
126 G4double G4UGenericPolycone::GetEndPhi() const 105 G4double G4UGenericPolycone::GetEndPhi() const
127 { 106 {
128 return (wrStart + wrDelta); << 107 return GetShape()->GetEndPhi();
129 } <<
130 G4double G4UGenericPolycone::GetSinStartPhi() <<
131 { <<
132 if (IsOpen()) return 0.; <<
133 G4double phi = GetStartPhi(); <<
134 return std::sin(phi); <<
135 } <<
136 G4double G4UGenericPolycone::GetCosStartPhi() <<
137 { <<
138 if (IsOpen()) return 1.; <<
139 G4double phi = GetStartPhi(); <<
140 return std::cos(phi); <<
141 } <<
142 G4double G4UGenericPolycone::GetSinEndPhi() co <<
143 { <<
144 if (IsOpen()) return 0.; <<
145 G4double phi = GetEndPhi(); <<
146 return std::sin(phi); <<
147 } <<
148 G4double G4UGenericPolycone::GetCosEndPhi() co <<
149 { <<
150 if (IsOpen()) return 1.; <<
151 G4double phi = GetEndPhi(); <<
152 return std::cos(phi); <<
153 } 108 }
154 G4bool G4UGenericPolycone::IsOpen() const 109 G4bool G4UGenericPolycone::IsOpen() const
155 { 110 {
156 return (wrDelta < twopi); << 111 return GetShape()->IsOpen();
157 } 112 }
158 G4int G4UGenericPolycone::GetNumRZCorner() con 113 G4int G4UGenericPolycone::GetNumRZCorner() const
159 { 114 {
160 return rzcorners.size(); << 115 return GetShape()->GetNumRZCorner();
161 } 116 }
162 G4PolyconeSideRZ G4UGenericPolycone::GetCorner 117 G4PolyconeSideRZ G4UGenericPolycone::GetCorner(G4int index) const
163 { 118 {
164 G4TwoVector rz = rzcorners.at(index); << 119 UPolyconeSideRZ pside = GetShape()->GetCorner(index);
165 G4PolyconeSideRZ psiderz = { rz.x(), rz.y() << 120 G4PolyconeSideRZ psiderz = { pside.r, pside.z };
166 121
167 return psiderz; 122 return psiderz;
168 } 123 }
169 124
170 ////////////////////////////////////////////// << 125 G4Polyhedron* G4UGenericPolycone::CreatePolyhedron() const
171 // <<
172 // Make a clone of the object <<
173 <<
174 G4VSolid* G4UGenericPolycone::Clone() const <<
175 { 126 {
176 return new G4UGenericPolycone(*this); <<
177 } <<
178 <<
179 ////////////////////////////////////////////// <<
180 // <<
181 // Get bounding box <<
182 127
183 void <<
184 G4UGenericPolycone::BoundingLimits(G4ThreeVect <<
185 G4ThreeVect <<
186 { <<
187 G4double rmin = kInfinity, rmax = -kInfinity <<
188 G4double zmin = kInfinity, zmax = -kInfinity <<
189 128
190 for (G4int i=0; i<GetNumRZCorner(); ++i) << 129 // The following code prepares for:
191 { << 130 // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
192 G4PolyconeSideRZ corner = GetCorner(i); << 131 // const double xyz[][3],
193 if (corner.r < rmin) rmin = corner.r; << 132 // const int faces_vec[][4])
194 if (corner.r > rmax) rmax = corner.r; << 133 // Here is an extract from the header file HepPolyhedron.h:
195 if (corner.z < zmin) zmin = corner.z; << 134 /**
196 if (corner.z > zmax) zmax = corner.z; << 135 * Creates user defined polyhedron.
197 } << 136 * This function allows to the user to define arbitrary polyhedron.
198 << 137 * The faces of the polyhedron should be either triangles or planar
199 if (IsOpen()) << 138 * quadrilateral. Nodes of a face are defined by indexes pointing to
200 { << 139 * the elements in the xyz array. Numeration of the elements in the
201 G4TwoVector vmin,vmax; << 140 * array starts from 1 (like in fortran). The indexes can be positive
202 G4GeomTools::DiskExtent(rmin,rmax, << 141 * or negative. Negative sign means that the corresponding edge is
203 GetSinStartPhi(),G << 142 * invisible. The normal of the face should be directed to exterior
204 GetSinEndPhi(),Get << 143 * of the polyhedron.
205 vmin,vmax); << 144 *
206 pMin.set(vmin.x(),vmin.y(),zmin); << 145 * @param Nnodes number of nodes
207 pMax.set(vmax.x(),vmax.y(),zmax); << 146 * @param Nfaces number of faces
208 } << 147 * @param xyz nodes
209 else << 148 * @param faces_vec faces (quadrilaterals or triangles)
210 { << 149 * @return status of the operation - is non-zero in case of problem
211 pMin.set(-rmax,-rmax, zmin); << 150 */
212 pMax.set( rmax, rmax, zmax); << 151 const G4int numSide =
213 } << 152 G4int(G4Polyhedron::GetNumberOfRotationSteps()
214 << 153 * (GetEndPhi() - GetStartPhi()) / twopi) + 1;
215 // Check correctness of the bounding box << 154 G4int nNodes;
216 // << 155 G4int nFaces;
217 if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 156 typedef G4double double3[3];
218 { << 157 double3* xyz;
219 std::ostringstream message; << 158 typedef G4int int4[4];
220 message << "Bad bounding box (min >= max) << 159 int4* faces_vec;
221 << GetName() << " !" << 160 if (IsOpen())
222 << "\npMin = " << pMin << 161 {
223 << "\npMax = " << pMax; << 162 // Triangulate open ends. Simple ear-chopping algorithm...
224 G4Exception("G4UGenericPolycone::BoundingL << 163 // I'm not sure how robust this algorithm is (J.Allison).
225 JustWarning, message); << 164 //
226 StreamInfo(G4cout); << 165 std::vector<G4bool> chopped(GetNumRZCorner(), false);
227 } << 166 std::vector<G4int*> triQuads;
228 } << 167 G4int remaining = GetNumRZCorner();
>> 168 G4int iStarter = 0;
>> 169 while (remaining >= 3) // Loop checking, 13.08.2015, G.Cosmo
>> 170 {
>> 171 // Find unchopped corners...
>> 172 //
>> 173 G4int A = -1, B = -1, C = -1;
>> 174 G4int iStepper = iStarter;
>> 175 do // Loop checking, 13.08.2015, G.Cosmo
>> 176 {
>> 177 if (A < 0) { A = iStepper; }
>> 178 else if (B < 0) { B = iStepper; }
>> 179 else if (C < 0) { C = iStepper; }
>> 180 do // Loop checking, 13.08.2015, G.Cosmo
>> 181 {
>> 182 if (++iStepper >= GetNumRZCorner()) { iStepper = 0; }
>> 183 }
>> 184 while (chopped[iStepper]);
>> 185 }
>> 186 while (C < 0 && iStepper != iStarter);
>> 187
>> 188 // Check triangle at B is pointing outward (an "ear").
>> 189 // Sign of z cross product determines...
>> 190 //
>> 191 G4double BAr = GetCorner(A).r - GetCorner(B).r;
>> 192 G4double BAz = GetCorner(A).z - GetCorner(B).z;
>> 193 G4double BCr = GetCorner(C).r - GetCorner(B).r;
>> 194 G4double BCz = GetCorner(C).z - GetCorner(B).z;
>> 195 if (BAr * BCz - BAz * BCr < kCarTolerance)
>> 196 {
>> 197 G4int* tq = new G4int[3];
>> 198 tq[0] = A + 1;
>> 199 tq[1] = B + 1;
>> 200 tq[2] = C + 1;
>> 201 triQuads.push_back(tq);
>> 202 chopped[B] = true;
>> 203 --remaining;
>> 204 }
>> 205 else
>> 206 {
>> 207 do // Loop checking, 13.08.2015, G.Cosmo
>> 208 {
>> 209 if (++iStarter >= GetNumRZCorner()) { iStarter = 0; }
>> 210 }
>> 211 while (chopped[iStarter]);
>> 212 }
>> 213 }
>> 214 // Transfer to faces...
>> 215 //
>> 216 nNodes = (numSide + 1) * GetNumRZCorner();
>> 217 nFaces = numSide * GetNumRZCorner() + 2 * triQuads.size();
>> 218 faces_vec = new int4[nFaces];
>> 219 G4int iface = 0;
>> 220 G4int addition = GetNumRZCorner() * numSide;
>> 221 G4int d = GetNumRZCorner() - 1;
>> 222 for (G4int iEnd = 0; iEnd < 2; ++iEnd)
>> 223 {
>> 224 for (size_t i = 0; i < triQuads.size(); ++i)
>> 225 {
>> 226 // Negative for soft/auxiliary/normally invisible edges...
>> 227 //
>> 228 G4int a, b, c;
>> 229 if (iEnd == 0)
>> 230 {
>> 231 a = triQuads[i][0];
>> 232 b = triQuads[i][1];
>> 233 c = triQuads[i][2];
>> 234 }
>> 235 else
>> 236 {
>> 237 a = triQuads[i][0] + addition;
>> 238 b = triQuads[i][2] + addition;
>> 239 c = triQuads[i][1] + addition;
>> 240 }
>> 241 G4int ab = std::abs(b - a);
>> 242 G4int bc = std::abs(c - b);
>> 243 G4int ca = std::abs(a - c);
>> 244 faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
>> 245 faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
>> 246 faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
>> 247 faces_vec[iface][3] = 0;
>> 248 ++iface;
>> 249 }
>> 250 }
229 251
230 ////////////////////////////////////////////// << 252 // Continue with sides...
231 // <<
232 // Calculate extent under transform and specif <<
233 253
234 G4bool << 254 xyz = new double3[nNodes];
235 G4UGenericPolycone::CalculateExtent(const EAxi << 255 const G4double dPhi = (GetEndPhi() - GetStartPhi()) / numSide;
236 const G4Vo << 256 G4double phi = GetStartPhi();
237 const G4Af << 257 G4int ixyz = 0;
238 G4do << 258 for (G4int iSide = 0; iSide < numSide; ++iSide)
239 { << 259 {
240 G4ThreeVector bmin, bmax; << 260 for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
241 G4bool exist; << 261 {
>> 262 xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
>> 263 xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
>> 264 xyz[ixyz][2] = GetCorner(iCorner).z;
>> 265 if (iSide == 0) // startPhi
>> 266 {
>> 267 if (iCorner < GetNumRZCorner() - 1)
>> 268 {
>> 269 faces_vec[iface][0] = ixyz + 1;
>> 270 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
>> 271 faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
>> 272 faces_vec[iface][3] = ixyz + 2;
>> 273 }
>> 274 else
>> 275 {
>> 276 faces_vec[iface][0] = ixyz + 1;
>> 277 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
>> 278 faces_vec[iface][2] = ixyz + 2;
>> 279 faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
>> 280 }
>> 281 }
>> 282 else if (iSide == numSide - 1) // endPhi
>> 283 {
>> 284 if (iCorner < GetNumRZCorner() - 1)
>> 285 {
>> 286 faces_vec[iface][0] = ixyz + 1;
>> 287 faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
>> 288 faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
>> 289 faces_vec[iface][3] = -(ixyz + 2);
>> 290 }
>> 291 else
>> 292 {
>> 293 faces_vec[iface][0] = ixyz + 1;
>> 294 faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
>> 295 faces_vec[iface][2] = ixyz + 2;
>> 296 faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
>> 297 }
>> 298 }
>> 299 else
>> 300 {
>> 301 if (iCorner < GetNumRZCorner() - 1)
>> 302 {
>> 303 faces_vec[iface][0] = ixyz + 1;
>> 304 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
>> 305 faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
>> 306 faces_vec[iface][3] = -(ixyz + 2);
>> 307 }
>> 308 else
>> 309 {
>> 310 faces_vec[iface][0] = ixyz + 1;
>> 311 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
>> 312 faces_vec[iface][2] = ixyz + 2;
>> 313 faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
>> 314 }
>> 315 }
>> 316 ++iface;
>> 317 ++ixyz;
>> 318 }
>> 319 phi += dPhi;
>> 320 }
242 321
243 // Check bounding box (bbox) << 322 // Last corners...
244 // <<
245 BoundingLimits(bmin,bmax); <<
246 G4BoundingEnvelope bbox(bmin,bmax); <<
247 #ifdef G4BBOX_EXTENT <<
248 return bbox.CalculateExtent(pAxis,pVoxelLimi <<
249 #endif <<
250 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox <<
251 { <<
252 return exist = pMin < pMax; <<
253 } <<
254 <<
255 // To find the extent, RZ contour of the pol <<
256 // in triangles. The extent is calculated as <<
257 // all sub-polycones formed by rotation of t <<
258 // <<
259 G4TwoVectorList contourRZ; <<
260 G4TwoVectorList triangles; <<
261 G4double eminlim = pVoxelLimit.GetMinExtent( <<
262 G4double emaxlim = pVoxelLimit.GetMaxExtent( <<
263 <<
264 // get RZ contour, ensure anticlockwise orde <<
265 for (G4int i=0; i<GetNumRZCorner(); ++i) <<
266 { <<
267 G4PolyconeSideRZ corner = GetCorner(i); <<
268 contourRZ.emplace_back(corner.r,corner.z); <<
269 } <<
270 G4double area = G4GeomTools::PolygonArea(con <<
271 if (area < 0.) std::reverse(contourRZ.begin( <<
272 <<
273 // triangulate RZ countour <<
274 if (!G4GeomTools::TriangulatePolygon(contour <<
275 { <<
276 std::ostringstream message; <<
277 message << "Triangulation of RZ contour ha <<
278 << GetName() << " !" <<
279 << "\nExtent has been calculated u <<
280 G4Exception("G4UGenericPolycone::Calculate <<
281 "GeomMgt1002", JustWarning, me <<
282 return bbox.CalculateExtent(pAxis,pVoxelLi <<
283 } <<
284 <<
285 // set trigonometric values <<
286 const G4int NSTEPS = 24; // numbe <<
287 G4double astep = twopi/NSTEPS; // max a <<
288 <<
289 G4double sphi = GetStartPhi(); <<
290 G4double ephi = GetEndPhi(); <<
291 G4double dphi = IsOpen() ? ephi-sphi : two <<
292 G4int ksteps = (dphi <= astep) ? 1 : (G4i <<
293 G4double ang = dphi/ksteps; <<
294 <<
295 G4double sinHalf = std::sin(0.5*ang); <<
296 G4double cosHalf = std::cos(0.5*ang); <<
297 G4double sinStep = 2.*sinHalf*cosHalf; <<
298 G4double cosStep = 1. - 2.*sinHalf*sinHalf; <<
299 <<
300 G4double sinStart = GetSinStartPhi(); <<
301 G4double cosStart = GetCosStartPhi(); <<
302 G4double sinEnd = GetSinEndPhi(); <<
303 G4double cosEnd = GetCosEndPhi(); <<
304 <<
305 // define vectors and arrays <<
306 std::vector<const G4ThreeVectorList *> polyg <<
307 polygons.resize(ksteps+2); <<
308 G4ThreeVectorList pols[NSTEPS+2]; <<
309 for (G4int k=0; k<ksteps+2; ++k) pols[k].res <<
310 for (G4int k=0; k<ksteps+2; ++k) polygons[k] <<
311 G4double r0[6],z0[6]; // contour with origin <<
312 G4double r1[6]; // shifted radii of ex <<
313 <<
314 // main loop along triangles <<
315 pMin = kInfinity; <<
316 pMax =-kInfinity; <<
317 G4int ntria = triangles.size()/3; <<
318 for (G4int i=0; i<ntria; ++i) <<
319 { <<
320 G4int i3 = i*3; <<
321 for (G4int k=0; k<3; ++k) <<
322 { <<
323 G4int e0 = i3+k, e1 = (k<2) ? e0+1 : i3; <<
324 G4int k2 = k*2; <<
325 // set contour with original edges of tr <<
326 r0[k2+0] = triangles[e0].x(); z0[k2+0] = <<
327 r0[k2+1] = triangles[e1].x(); z0[k2+1] = <<
328 // set shifted radii <<
329 r1[k2+0] = r0[k2+0]; <<
330 r1[k2+1] = r0[k2+1]; <<
331 if (z0[k2+1] - z0[k2+0] <= 0) continue; <<
332 r1[k2+0] /= cosHalf; <<
333 r1[k2+1] /= cosHalf; <<
334 } <<
335 323
336 // rotate countour, set sequence of 6-side << 324 for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
337 G4double sinCur = sinStart*cosHalf + cosSt << 325 {
338 G4double cosCur = cosStart*cosHalf - sinSt << 326 xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
339 for (G4int j=0; j<6; ++j) << 327 xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
340 { << 328 xyz[ixyz][2] = GetCorner(iCorner).z;
341 pols[0][j].set(r0[j]*cosStart,r0[j]*sinS << 329 ++ixyz;
>> 330 }
342 } 331 }
343 for (G4int k=1; k<ksteps+1; ++k) << 332 else // !phiIsOpen - i.e., a complete 360 degrees.
344 { 333 {
345 for (G4int j=0; j<6; ++j) << 334 nNodes = numSide * GetNumRZCorner();
>> 335 nFaces = numSide * GetNumRZCorner();;
>> 336 xyz = new double3[nNodes];
>> 337 faces_vec = new int4[nFaces];
>> 338 const G4double dPhi = (GetEndPhi() - GetStartPhi()) / numSide;
>> 339 G4double phi = GetStartPhi();
>> 340 G4int ixyz = 0, iface = 0;
>> 341 for (G4int iSide = 0; iSide < numSide; ++iSide)
346 { 342 {
347 pols[k][j].set(r1[j]*cosCur,r1[j]*sinC << 343 for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
>> 344 {
>> 345 xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
>> 346 xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
>> 347 xyz[ixyz][2] = GetCorner(iCorner).z;
>> 348
>> 349 if (iSide < numSide - 1)
>> 350 {
>> 351 if (iCorner < GetNumRZCorner() - 1)
>> 352 {
>> 353 faces_vec[iface][0] = ixyz + 1;
>> 354 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
>> 355 faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
>> 356 faces_vec[iface][3] = -(ixyz + 2);
>> 357 }
>> 358 else
>> 359 {
>> 360 faces_vec[iface][0] = ixyz + 1;
>> 361 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
>> 362 faces_vec[iface][2] = ixyz + 2;
>> 363 faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
>> 364 }
>> 365 }
>> 366 else // Last side joins ends...
>> 367 {
>> 368 if (iCorner < GetNumRZCorner() - 1)
>> 369 {
>> 370 faces_vec[iface][0] = ixyz + 1;
>> 371 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() - nFaces + 1);
>> 372 faces_vec[iface][2] = ixyz + GetNumRZCorner() - nFaces + 2;
>> 373 faces_vec[iface][3] = -(ixyz + 2);
>> 374 }
>> 375 else
>> 376 {
>> 377 faces_vec[iface][0] = ixyz + 1;
>> 378 faces_vec[iface][1] = -(ixyz - nFaces + GetNumRZCorner() + 1);
>> 379 faces_vec[iface][2] = ixyz - nFaces + 2;
>> 380 faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
>> 381 }
>> 382 }
>> 383 ++ixyz;
>> 384 ++iface;
>> 385 }
>> 386 phi += dPhi;
348 } 387 }
349 G4double sinTmp = sinCur; <<
350 sinCur = sinCur*cosStep + cosCur*sinStep <<
351 cosCur = cosCur*cosStep - sinTmp*sinStep <<
352 } 388 }
353 for (G4int j=0; j<6; ++j) << 389 G4Polyhedron* polyhedron = new G4Polyhedron;
>> 390 G4int problem = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
>> 391 delete [] faces_vec;
>> 392 delete [] xyz;
>> 393 if (problem)
354 { 394 {
355 pols[ksteps+1][j].set(r0[j]*cosEnd,r0[j] << 395 std::ostringstream message;
>> 396 message << "Problem creating G4Polyhedron for: " << GetName();
>> 397 G4Exception("G4GenericPolycone::CreatePolyhedron()", "GeomSolids1002",
>> 398 JustWarning, message);
>> 399 delete polyhedron;
>> 400 return 0;
>> 401 }
>> 402 else
>> 403 {
>> 404 return polyhedron;
356 } 405 }
357 <<
358 // set sub-envelope and adjust extent <<
359 G4double emin,emax; <<
360 G4BoundingEnvelope benv(polygons); <<
361 if (!benv.CalculateExtent(pAxis,pVoxelLimi <<
362 if (emin < pMin) pMin = emin; <<
363 if (emax > pMax) pMax = emax; <<
364 if (eminlim > pMin && emaxlim < pMax) retu <<
365 } <<
366 return (pMin < pMax); <<
367 } <<
368 <<
369 ////////////////////////////////////////////// <<
370 // <<
371 // CreatePolyhedron <<
372 <<
373 G4Polyhedron* G4UGenericPolycone::CreatePolyhe <<
374 { <<
375 return new G4PolyhedronPcon(wrStart, wrDelta <<
376 } 406 }
377 407
378 #endif // G4GEOM_USE_USOLIDS 408 #endif // G4GEOM_USE_USOLIDS
379 409