| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // G4VSolid implementation for solid base clas <<
27 // 23 //
28 // 10.10.18 E.Tcherniaev, more robust Estimate << 24 // $Id: G4VSolid.cc,v 1.7.2.1 2001/06/28 19:08:32 gunter Exp $
29 // 30.06.95 P.Kent, Created. << 25 // GEANT4 tag $Name: $
30 // ------------------------------------------- << 26 //
>> 27 // class G4VSolid
>> 28 //
>> 29 // Implementation for solid base class
>> 30 //
>> 31 //
>> 32 // History:
>> 33 // 10.07.95 P.Kent Added == operator, solid Store entry
>> 34 // 30.06.95 P.Kent
>> 35 // 15.11.00 D.Williams, V.Grichine change in CalculateClippedPolygonExtent:
>> 36 // else if(component>pMax) ---> if
31 37
32 #include "G4VSolid.hh" 38 #include "G4VSolid.hh"
33 #include "G4SolidStore.hh" 39 #include "G4SolidStore.hh"
34 #include "globals.hh" <<
35 #include "G4QuickRand.hh" <<
36 #include "G4GeometryTolerance.hh" <<
37 40
38 #include "G4VoxelLimits.hh" 41 #include "G4VoxelLimits.hh"
39 #include "G4AffineTransform.hh" 42 #include "G4AffineTransform.hh"
40 #include "G4VisExtent.hh" 43 #include "G4VisExtent.hh"
41 44
42 ////////////////////////////////////////////// <<
43 // <<
44 // Streaming operator dumping solid contents <<
45 <<
46 std::ostream& operator<< ( std::ostream& os, c <<
47 { <<
48 return e.StreamInfo(os); <<
49 } <<
50 <<
51 ////////////////////////////////////////////// <<
52 // <<
53 // Constructor 45 // Constructor
54 // - Copies name 46 // - Copies name
55 // - Add ourselves to solid Store 47 // - Add ourselves to solid Store
56 << 48 G4VSolid::G4VSolid(const G4String& name) :
57 G4VSolid::G4VSolid(const G4String& name) << 49 fshapeName(name)
58 : fshapeName(name) <<
59 { <<
60 kCarTolerance = G4GeometryTolerance::GetIn <<
61 <<
62 // Register to store <<
63 // <<
64 G4SolidStore::GetInstance()->Register(this <<
65 } <<
66 <<
67 ////////////////////////////////////////////// <<
68 // <<
69 // Copy constructor <<
70 // <<
71 <<
72 G4VSolid::G4VSolid(const G4VSolid& rhs) <<
73 : kCarTolerance(rhs.kCarTolerance), fshapeNa <<
74 { <<
75 // Register to store <<
76 // <<
77 G4SolidStore::GetInstance()->Register(this <<
78 } <<
79 <<
80 ////////////////////////////////////////////// <<
81 // <<
82 // Fake default constructor - sets only member <<
83 // for usage restri <<
84 // <<
85 G4VSolid::G4VSolid( __void__& ) <<
86 : fshapeName("") <<
87 { 50 {
88 // Register to store << 51 G4SolidStore::GetInstance()->push_back(this);
89 // <<
90 G4SolidStore::GetInstance()->Register(this <<
91 } 52 }
92 53
93 ////////////////////////////////////////////// <<
94 // <<
95 // Destructor (virtual) 54 // Destructor (virtual)
96 // - Remove ourselves from solid Store 55 // - Remove ourselves from solid Store
97 <<
98 G4VSolid::~G4VSolid() 56 G4VSolid::~G4VSolid()
99 { 57 {
100 G4SolidStore::GetInstance()->DeRegister(th 58 G4SolidStore::GetInstance()->DeRegister(this);
101 } 59 }
102 << 60
103 ////////////////////////////////////////////// << 61 // Returns name by value
104 // << 62 G4String G4VSolid::GetName() const
105 // Assignment operator <<
106 <<
107 G4VSolid& G4VSolid::operator = (const G4VSolid <<
108 { 63 {
109 // Check assignment to self << 64 return fshapeName;
110 // << 65 }
111 if (this == &rhs) { return *this; } <<
112 <<
113 // Copy data <<
114 // <<
115 kCarTolerance = rhs.kCarTolerance; <<
116 fshapeName = rhs.fshapeName; <<
117 <<
118 return *this; <<
119 } <<
120 <<
121 <<
122 <<
123 ////////////////////////////////////////////// <<
124 // <<
125 // Set solid name and notify store of the chan <<
126 66
127 void G4VSolid::SetName(const G4String& name) 67 void G4VSolid::SetName(const G4String& name)
128 { 68 {
129 fshapeName = name; << 69 fshapeName=name;
130 G4SolidStore::GetInstance()->SetMapValid(fal << 70 }
131 } <<
132 71
133 ////////////////////////////////////////////// <<
134 // <<
135 // Throw exception if ComputeDimensions called 72 // Throw exception if ComputeDimensions called for illegal derived class
136 << 73 void G4VSolid::ComputeDimensions(G4VPVParameterisation* p,
137 void G4VSolid::ComputeDimensions(G4VPVParamete << 74 const G4int n,
138 const G4int, << 75 const G4VPhysicalVolume* pRep)
139 const G4VPhys <<
140 { <<
141 std::ostringstream message; <<
142 message << "Illegal call to G4VSolid::Comp <<
143 << "Method not overloaded by deriv <<
144 G4Exception("G4VSolid::ComputeDimensions() <<
145 FatalException, message); <<
146 } <<
147 <<
148 ////////////////////////////////////////////// <<
149 // <<
150 // Throw exception (warning) for solids not im <<
151 <<
152 G4ThreeVector G4VSolid::GetPointOnSurface() co <<
153 { 76 {
154 std::ostringstream message; << 77 G4Exception("G4VSolid::ComputeDimensions called illegally: not overloaded by derived class");
155 message << "Not implemented for solid: " <<
156 << GetEntityType() << " !" << G4en <<
157 << "Returning origin."; <<
158 G4Exception("G4VSolid::GetPointOnSurface() <<
159 JustWarning, message); <<
160 return {0,0,0}; <<
161 } 78 }
162 79
163 ////////////////////////////////////////////// << 80 ///////////////////////////////////////////////////////////////////////////////
164 // <<
165 // Returns total number of constituents that w <<
166 // of the solid. For non-Boolean solids the re <<
167 <<
168 G4int G4VSolid::GetNumOfConstituents() const <<
169 { return 1; } <<
170 <<
171 ////////////////////////////////////////////// <<
172 // <<
173 // Returns true if the solid has only planar f <<
174 <<
175 G4bool G4VSolid::IsFaceted() const <<
176 { return false; } <<
177 <<
178 ////////////////////////////////////////////// <<
179 // <<
180 // Dummy implementations ... <<
181 <<
182 const G4VSolid* G4VSolid::GetConstituentSolid( <<
183 { return nullptr; } <<
184 <<
185 G4VSolid* G4VSolid::GetConstituentSolid(G4int) <<
186 { return nullptr; } <<
187 <<
188 const G4DisplacedSolid* G4VSolid::GetDisplaced <<
189 { return nullptr; } <<
190 <<
191 G4DisplacedSolid* G4VSolid::GetDisplacedSolidP <<
192 { return nullptr; } <<
193 <<
194 ////////////////////////////////////////////// <<
195 // 81 //
196 // Returns an estimation of the solid volume i << 82 // Calculate the maximum and minimum extents of the convex polygon pPolygon
197 // The number of statistics and error accuracy << 83 // along the axis pAxis, within the limits pVoxelLimit
198 // This method may be overloaded by derived cl <<
199 // exact geometrical quantity for solids where <<
200 // or anyway to cache the computed value. <<
201 // This implementation does NOT cache the comp <<
202 <<
203 G4double G4VSolid::GetCubicVolume() <<
204 { <<
205 G4int cubVolStatistics = 1000000; <<
206 G4double cubVolEpsilon = 0.001; <<
207 return EstimateCubicVolume(cubVolStatistics, <<
208 } <<
209 <<
210 ////////////////////////////////////////////// <<
211 // 84 //
212 // Calculate cubic volume based on Inside() me <<
213 // Accuracy is limited by the second argument <<
214 // expressed by the first argument. <<
215 // Implementation is courtesy of Vasiliki Desp <<
216 // University of Athens. <<
217 85
218 G4double G4VSolid::EstimateCubicVolume(G4int n << 86 void G4VSolid::CalculateClippedPolygonExtent(G4ThreeVectorList& pPolygon,
>> 87 const G4VoxelLimits& pVoxelLimit,
>> 88 const EAxis pAxis,
>> 89 G4double& pMin, G4double& pMax) const
219 { 90 {
220 G4int iInside=0; << 91 G4int noLeft,i;
221 G4double px,py,pz,minX,maxX,minY,maxY,minZ,m << 92 G4double component;
222 G4ThreeVector p; << 93 ClipPolygon(pPolygon,pVoxelLimit);
223 EInside in; << 94 noLeft = pPolygon.size();
224 <<
225 // values needed for CalculateExtent signatu <<
226 <<
227 G4VoxelLimits limit; // Unlim <<
228 G4AffineTransform origin; <<
229 <<
230 // min max extents of pSolid along X,Y,Z <<
231 <<
232 CalculateExtent(kXAxis,limit,origin,minX,max <<
233 CalculateExtent(kYAxis,limit,origin,minY,max <<
234 CalculateExtent(kZAxis,limit,origin,minZ,max <<
235 <<
236 // limits <<
237 <<
238 if(nStat < 100) nStat = 100; <<
239 if(epsilon > 0.01) epsilon = 0.01; <<
240 halfepsilon = 0.5*epsilon; <<
241 <<
242 for(auto i = 0; i < nStat; ++i ) <<
243 { <<
244 px = minX-halfepsilon+(maxX-minX+epsilon)* <<
245 py = minY-halfepsilon+(maxY-minY+epsilon)* <<
246 pz = minZ-halfepsilon+(maxZ-minZ+epsilon)* <<
247 p = G4ThreeVector(px,py,pz); <<
248 in = Inside(p); <<
249 if(in != kOutside) ++iInside; <<
250 } <<
251 volume = (maxX-minX+epsilon)*(maxY-minY+epsi <<
252 * (maxZ-minZ+epsilon)*iInside/nStat; <<
253 return volume; <<
254 } <<
255 <<
256 ////////////////////////////////////////////// <<
257 // <<
258 // Returns an estimation of the solid surface <<
259 // The number of statistics and error accuracy <<
260 // This method may be overloaded by derived cl <<
261 // exact geometrical quantity for solids where <<
262 // or anyway to cache the computed value. <<
263 // This implementation does NOT cache the comp <<
264 <<
265 G4double G4VSolid::GetSurfaceArea() <<
266 { <<
267 G4int stat = 1000000; <<
268 G4double ell = -1.; <<
269 return EstimateSurfaceArea(stat,ell); <<
270 } <<
271 95
272 ////////////////////////////////////////////// << 96 if (noLeft)
273 // <<
274 // Calculate surface area by estimating volume <<
275 // surrounding the surface using Monte-Carlo m <<
276 // Input parameters: <<
277 // nstat - statistics (number of random poi <<
278 // eps - shell thinkness <<
279 <<
280 G4double G4VSolid::EstimateSurfaceArea(G4int n <<
281 { <<
282 static const G4double s2 = 1./std::sqrt(2.); <<
283 static const G4double s3 = 1./std::sqrt(3.); <<
284 static const G4ThreeVector directions[64] = <<
285 { 97 {
286 G4ThreeVector( 0, 0, 0), G4ThreeVector( << 98 for (i=0;i<noLeft;i++)
287 G4ThreeVector( 1, 0, 0), G4ThreeVector( <<
288 G4ThreeVector( 0, -1, 0), G4ThreeVector( <<
289 G4ThreeVector( s2, -s2, 0), G4ThreeVector( <<
290 <<
291 G4ThreeVector( 0, 1, 0), G4ThreeVector( <<
292 G4ThreeVector( s2, s2, 0), G4ThreeVector( <<
293 G4ThreeVector( 0, -1, 0), G4ThreeVector( <<
294 G4ThreeVector( 1, 0, 0), G4ThreeVector( <<
295 <<
296 G4ThreeVector( 0, 0, -1), G4ThreeVector( <<
297 G4ThreeVector( s2, 0,-s2), G4ThreeVector( <<
298 G4ThreeVector( 0,-s2,-s2), G4ThreeVector( <<
299 G4ThreeVector( s3,-s3,-s3), G4ThreeVector( <<
300 <<
301 G4ThreeVector( 0, s2,-s2), G4ThreeVector( <<
302 G4ThreeVector( s3, s3,-s3), G4ThreeVector( <<
303 G4ThreeVector( 0, 0, -1), G4ThreeVector( <<
304 G4ThreeVector( s2, 0,-s2), G4ThreeVector( <<
305 <<
306 G4ThreeVector( 0, 0, 1), G4ThreeVector( <<
307 G4ThreeVector( s2, 0, s2), G4ThreeVector( <<
308 G4ThreeVector( 0,-s2, s2), G4ThreeVector( <<
309 G4ThreeVector( s3,-s3, s3), G4ThreeVector( <<
310 <<
311 G4ThreeVector( 0, s2, s2), G4ThreeVector( <<
312 G4ThreeVector( s3, s3, s3), G4ThreeVector( <<
313 G4ThreeVector( 0, 0, 1), G4ThreeVector( <<
314 G4ThreeVector( s2, 0, s2), G4ThreeVector( <<
315 <<
316 G4ThreeVector( 0, 0, -1), G4ThreeVector( <<
317 G4ThreeVector( 1, 0, 0), G4ThreeVector( <<
318 G4ThreeVector( 0, -1, 0), G4ThreeVector( <<
319 G4ThreeVector( s2, -s2, 0), G4ThreeVector( <<
320 <<
321 G4ThreeVector( 0, 1, 0), G4ThreeVector( <<
322 G4ThreeVector( s2, s2, 0), G4ThreeVector( <<
323 G4ThreeVector( 0, -1, 0), G4ThreeVector( <<
324 G4ThreeVector( 1, 0, 0), G4ThreeVector( <<
325 }; <<
326 <<
327 G4ThreeVector bmin, bmax; <<
328 BoundingLimits(bmin, bmax); <<
329 <<
330 G4double dX = bmax.x() - bmin.x(); <<
331 G4double dY = bmax.y() - bmin.y(); <<
332 G4double dZ = bmax.z() - bmin.z(); <<
333 <<
334 // Define statistics and shell thickness <<
335 // <<
336 G4int npoints = (nstat < 1000) ? 1000 : nsta <<
337 G4double coeff = 0.5 / std::cbrt(G4double(np <<
338 G4double eps = (ell > 0) ? ell : coeff * std <<
339 G4double del = 1.8 * eps; // shold be more t <<
340 <<
341 G4double minX = bmin.x() - eps; <<
342 G4double minY = bmin.y() - eps; <<
343 G4double minZ = bmin.z() - eps; <<
344 <<
345 G4double dd = 2. * eps; <<
346 dX += dd; <<
347 dY += dd; <<
348 dZ += dd; <<
349 <<
350 // Calculate surface area <<
351 // <<
352 G4int icount = 0; <<
353 for(auto i = 0; i < npoints; ++i) <<
354 { <<
355 G4double px = minX + dX*G4QuickRand(); <<
356 G4double py = minY + dY*G4QuickRand(); <<
357 G4double pz = minZ + dZ*G4QuickRand(); <<
358 G4ThreeVector p = G4ThreeVector(px, py, p <<
359 EInside in = Inside(p); <<
360 G4double dist = 0; <<
361 if (in == kInside) <<
362 { <<
363 if (DistanceToOut(p) >= eps) continue; <<
364 G4int icase = 0; <<
365 if (Inside(G4ThreeVector(px-del, py, pz) <<
366 if (Inside(G4ThreeVector(px+del, py, pz) <<
367 if (Inside(G4ThreeVector(px, py-del, pz) <<
368 if (Inside(G4ThreeVector(px, py+del, pz) <<
369 if (Inside(G4ThreeVector(px, py, pz-del) <<
370 if (Inside(G4ThreeVector(px, py, pz+del) <<
371 if (icase == 0) continue; <<
372 G4ThreeVector v = directions[icase]; <<
373 dist = DistanceToOut(p, v); <<
374 G4ThreeVector n = SurfaceNormal(p + v*di <<
375 dist *= v.dot(n); <<
376 } <<
377 else if (in == kOutside) <<
378 { 99 {
379 if (DistanceToIn(p) >= eps) continue; << 100 component = pPolygon[i].operator()(pAxis);
380 G4int icase = 0; << 101
381 if (Inside(G4ThreeVector(px-del, py, pz) << 102 if (component < pMin)
382 if (Inside(G4ThreeVector(px+del, py, pz) << 103 {
383 if (Inside(G4ThreeVector(px, py-del, pz) << 104 pMin = component;
384 if (Inside(G4ThreeVector(px, py+del, pz) << 105 }
385 if (Inside(G4ThreeVector(px, py, pz-del) << 106 // else
386 if (Inside(G4ThreeVector(px, py, pz+del) << 107 if (component > pMax)
387 if (icase == 0) continue; << 108 {
388 G4ThreeVector v = directions[icase]; << 109 pMax = component;
389 dist = DistanceToIn(p, v); << 110 }
390 if (dist == kInfinity) continue; <<
391 G4ThreeVector n = SurfaceNormal(p + v*di <<
392 dist *= -(v.dot(n)); <<
393 } 111 }
394 if (dist < eps) ++icount; <<
395 } 112 }
396 return dX*dY*dZ*icount/npoints/dd; <<
397 } 113 }
398 114
399 ////////////////////////////////////////////// 115 ///////////////////////////////////////////////////////////////////////////
400 // << 116 //
401 // Returns a pointer of a dynamically allocate <<
402 // Returns NULL pointer with warning in case t <<
403 // implement this method. The caller has respo <<
404 // <<
405 <<
406 G4VSolid* G4VSolid::Clone() const <<
407 { <<
408 std::ostringstream message; <<
409 message << "Clone() method not implemented f <<
410 << GetEntityType() << "!" << G4endl <<
411 << "Returning NULL pointer!"; <<
412 G4Exception("G4VSolid::Clone()", "GeomMgt100 <<
413 return nullptr; <<
414 } <<
415 <<
416 ////////////////////////////////////////////// <<
417 // <<
418 // Calculate the maximum and minimum extents o 117 // Calculate the maximum and minimum extents of the polygon described
419 // by the vertices: pSectionIndex->pSectionInd 118 // by the vertices: pSectionIndex->pSectionIndex+1->
420 // pSectionIndex+2->pSection 119 // pSectionIndex+2->pSectionIndex+3->pSectionIndex
421 // in the List pVertices 120 // in the List pVertices
422 // 121 //
423 // If the minimum is <pMin pMin is set to the 122 // If the minimum is <pMin pMin is set to the new minimum
424 // If the maximum is >pMax pMax is set to the 123 // If the maximum is >pMax pMax is set to the new maximum
425 // 124 //
426 // No modifications are made to pVertices 125 // No modifications are made to pVertices
427 // << 126 void G4VSolid::ClipCrossSection(G4ThreeVectorList* pVertices,
428 << 127 const G4int pSectionIndex,
429 void G4VSolid::ClipCrossSection( G4Three << 128 const G4VoxelLimits& pVoxelLimit,
430 const G4int p << 129 const EAxis pAxis,
431 const G4Voxel << 130 G4double& pMin, G4double& pMax) const
432 const EAxis p <<
433 G4doubl <<
434 { 131 {
435 132
436 G4ThreeVectorList polygon; << 133 G4ThreeVectorList polygon;
437 polygon.reserve(4); << 134 polygon.push_back((*pVertices)[pSectionIndex]);
438 polygon.push_back((*pVertices)[pSectionIndex << 135 polygon.push_back((*pVertices)[pSectionIndex+1]);
439 polygon.push_back((*pVertices)[pSectionIndex << 136 polygon.push_back((*pVertices)[pSectionIndex+2]);
440 polygon.push_back((*pVertices)[pSectionIndex << 137 polygon.push_back((*pVertices)[pSectionIndex+3]);
441 polygon.push_back((*pVertices)[pSectionIndex << 138 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
442 CalculateClippedPolygonExtent(polygon,pVoxel << 139 return;
443 return; <<
444 } 140 }
445 141
446 ////////////////////////////////////////////// <<
447 // <<
448 // Calculate the maximum and minimum extents o 142 // Calculate the maximum and minimum extents of the polygons
449 // joining the CrossSections at pSectionIndex- 143 // joining the CrossSections at pSectionIndex->pSectionIndex+3 and
450 // pSectionIndex+ 144 // pSectionIndex+4->pSectionIndex7
451 // 145 //
452 // in the List pVertices, within the boundarie 146 // in the List pVertices, within the boundaries of the voxel limits pVoxelLimit
453 // 147 //
454 // If the minimum is <pMin pMin is set to the 148 // If the minimum is <pMin pMin is set to the new minimum
455 // If the maximum is >pMax pMax is set to the 149 // If the maximum is >pMax pMax is set to the new maximum
456 // 150 //
457 // No modifications are made to pVertices 151 // No modifications are made to pVertices
458 152
459 void G4VSolid::ClipBetweenSections( G4Thr << 153 void G4VSolid::ClipBetweenSections(G4ThreeVectorList* pVertices,
460 const G4int << 154 const G4int pSectionIndex,
461 const G4Vox << 155 const G4VoxelLimits& pVoxelLimit,
462 const EAxis << 156 const EAxis pAxis,
463 G4dou << 157 G4double& pMin, G4double& pMax) const
464 { << 158 {
465 G4ThreeVectorList polygon; << 159 G4ThreeVectorList polygon;
466 polygon.reserve(4); << 160 polygon.push_back((*pVertices)[pSectionIndex]);
467 polygon.push_back((*pVertices)[pSectionIndex << 161 polygon.push_back((*pVertices)[pSectionIndex+4]);
468 polygon.push_back((*pVertices)[pSectionIndex << 162 polygon.push_back((*pVertices)[pSectionIndex+5]);
469 polygon.push_back((*pVertices)[pSectionIndex << 163 polygon.push_back((*pVertices)[pSectionIndex+1]);
470 polygon.push_back((*pVertices)[pSectionIndex << 164 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
471 CalculateClippedPolygonExtent(polygon,pVoxel << 165 polygon.clear();
472 polygon.clear(); << 166 polygon.push_back((*pVertices)[pSectionIndex+1]);
473 << 167 polygon.push_back((*pVertices)[pSectionIndex+5]);
474 polygon.push_back((*pVertices)[pSectionIndex << 168 polygon.push_back((*pVertices)[pSectionIndex+6]);
475 polygon.push_back((*pVertices)[pSectionIndex << 169 polygon.push_back((*pVertices)[pSectionIndex+2]);
476 polygon.push_back((*pVertices)[pSectionIndex << 170 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
477 polygon.push_back((*pVertices)[pSectionIndex << 171 polygon.clear();
478 CalculateClippedPolygonExtent(polygon,pVoxel << 172 polygon.push_back((*pVertices)[pSectionIndex+2]);
479 polygon.clear(); << 173 polygon.push_back((*pVertices)[pSectionIndex+6]);
480 << 174 polygon.push_back((*pVertices)[pSectionIndex+7]);
481 polygon.push_back((*pVertices)[pSectionIndex << 175 polygon.push_back((*pVertices)[pSectionIndex+3]);
482 polygon.push_back((*pVertices)[pSectionIndex << 176 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
483 polygon.push_back((*pVertices)[pSectionIndex << 177 polygon.clear();
484 polygon.push_back((*pVertices)[pSectionIndex << 178 polygon.push_back((*pVertices)[pSectionIndex+3]);
485 CalculateClippedPolygonExtent(polygon,pVoxel << 179 polygon.push_back((*pVertices)[pSectionIndex+7]);
486 polygon.clear(); << 180 polygon.push_back((*pVertices)[pSectionIndex+4]);
487 << 181 polygon.push_back((*pVertices)[pSectionIndex]);
488 polygon.push_back((*pVertices)[pSectionIndex << 182 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
489 polygon.push_back((*pVertices)[pSectionIndex << 183 return;
490 polygon.push_back((*pVertices)[pSectionIndex <<
491 polygon.push_back((*pVertices)[pSectionIndex <<
492 CalculateClippedPolygonExtent(polygon,pVoxel <<
493 return; <<
494 } 184 }
495 185
496 <<
497 ////////////////////////////////////////////// <<
498 // <<
499 // Calculate the maximum and minimum extents o <<
500 // along the axis pAxis, within the limits pVo <<
501 // <<
502 <<
503 void <<
504 G4VSolid::CalculateClippedPolygonExtent(G4Thre <<
505 const G4Voxe <<
506 const EAxis <<
507 G4doub <<
508 G4doub <<
509 { <<
510 G4int noLeft,i; <<
511 G4double component; <<
512 <<
513 ClipPolygon(pPolygon,pVoxelLimit,pAxis); <<
514 noLeft = (G4int)pPolygon.size(); <<
515 <<
516 if ( noLeft != 0 ) <<
517 { <<
518 for (i=0; i<noLeft; ++i) <<
519 { <<
520 component = pPolygon[i].operator()(pAxis <<
521 <<
522 if (component < pMin) <<
523 { <<
524 pMin = component; <<
525 } <<
526 if (component > pMax) <<
527 { <<
528 pMax = component; <<
529 } <<
530 } <<
531 } <<
532 } <<
533 <<
534 ////////////////////////////////////////////// <<
535 // <<
536 // Clip the convex polygon described by the ve 186 // Clip the convex polygon described by the vertices at
537 // pSectionIndex ->pSectionIndex+3 within pVer 187 // pSectionIndex ->pSectionIndex+3 within pVertices to the limits pVoxelLimit
538 // 188 //
539 // Set pMin to the smallest 189 // Set pMin to the smallest
540 // 190 //
541 // Calculate the extent of the polygon along p 191 // Calculate the extent of the polygon along pAxis, when clipped to the
542 // limits pVoxelLimit. If the polygon exists a 192 // limits pVoxelLimit. If the polygon exists after clippin, set pMin to
543 // the polygon's minimum extent along the axis 193 // the polygon's minimum extent along the axis if <pMin, and set pMax to
544 // the polygon's maximum extent along the axis 194 // the polygon's maximum extent along the axis if >pMax.
545 // 195 //
546 // The polygon is described by a set of vector 196 // The polygon is described by a set of vectors, where each vector represents
547 // a vertex, so that the polygon is described 197 // a vertex, so that the polygon is described by the vertex sequence:
548 // 0th->1st 1st->2nd 2nd->... nth->0th 198 // 0th->1st 1st->2nd 2nd->... nth->0th
549 // 199 //
550 // Modifications to the polygon are made 200 // Modifications to the polygon are made
551 // 201 //
552 // NOTE: Execessive copying during clipping 202 // NOTE: Execessive copying during clipping
553 << 203 void G4VSolid::ClipPolygon(G4ThreeVectorList& pPolygon,
554 void G4VSolid::ClipPolygon( G4ThreeVector << 204 const G4VoxelLimits& pVoxelLimit) const
555 const G4VoxelLimits <<
556 const EAxis <<
557 { 205 {
558 G4ThreeVectorList outputPolygon; << 206 G4ThreeVectorList outputPolygon;
559 << 207 if (pVoxelLimit.IsLimited())
560 if ( pVoxelLimit.IsLimited() ) << 208 {
561 { << 209 if (pVoxelLimit.IsXLimited())
562 if (pVoxelLimit.IsXLimited() ) // && pAxis << 210 {
563 { << 211 G4VoxelLimits simpleLimit1;
564 G4VoxelLimits simpleLimit1; << 212 simpleLimit1.AddLimit(kXAxis,pVoxelLimit.GetMinXExtent(),kInfinity);
565 simpleLimit1.AddLimit(kXAxis,pVoxelLimit << 213 ClipPolygonToSimpleLimits(pPolygon,outputPolygon,
566 ClipPolygonToSimpleLimits(pPolygon,outpu << 214 simpleLimit1);
567 << 215 pPolygon.clear();
568 pPolygon.clear(); << 216 if (!outputPolygon.size())
569 << 217 {
570 if ( outputPolygon.empty() ) return; << 218 return;
571 << 219 }
572 G4VoxelLimits simpleLimit2; << 220
573 simpleLimit2.AddLimit(kXAxis,-kInfinity, << 221 G4VoxelLimits simpleLimit2;
574 ClipPolygonToSimpleLimits(outputPolygon, << 222 simpleLimit2.AddLimit(kXAxis,-kInfinity,pVoxelLimit.GetMaxXExtent());
575 << 223 ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
576 if ( pPolygon.empty() ) return; << 224 if (!pPolygon.size())
577 else outputPoly << 225 {
578 } << 226 return;
579 if ( pVoxelLimit.IsYLimited() ) // && pAxi << 227 }
580 { << 228 else
581 G4VoxelLimits simpleLimit1; << 229 {
582 simpleLimit1.AddLimit(kYAxis,pVoxelLimit << 230 outputPolygon.clear();
583 ClipPolygonToSimpleLimits(pPolygon,outpu << 231 }
584 << 232 }
585 // Must always clear pPolygon - for clip << 233
586 // early exit << 234 if (pVoxelLimit.IsYLimited())
587 << 235 {
588 pPolygon.clear(); << 236 G4VoxelLimits simpleLimit1;
589 << 237 simpleLimit1.AddLimit(kYAxis,pVoxelLimit.GetMinYExtent(),kInfinity);
590 if ( outputPolygon.empty() ) return; << 238 ClipPolygonToSimpleLimits(pPolygon,outputPolygon,
591 << 239 simpleLimit1);
592 G4VoxelLimits simpleLimit2; << 240 // Must always clear pPolygon - for clip to simpleLimit2 and incase of
593 simpleLimit2.AddLimit(kYAxis,-kInfinity, << 241 // early exit
594 ClipPolygonToSimpleLimits(outputPolygon, << 242 pPolygon.clear();
595 << 243 if (!outputPolygon.size())
596 if ( pPolygon.empty() ) return; << 244 {
597 else outputPoly << 245 return;
598 } << 246 }
599 if ( pVoxelLimit.IsZLimited() ) // && pAxi << 247
600 { << 248 G4VoxelLimits simpleLimit2;
601 G4VoxelLimits simpleLimit1; << 249 simpleLimit2.AddLimit(kYAxis,-kInfinity,pVoxelLimit.GetMaxYExtent());
602 simpleLimit1.AddLimit(kZAxis,pVoxelLimit << 250 ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
603 ClipPolygonToSimpleLimits(pPolygon,outpu << 251 if (!pPolygon.size())
604 << 252 {
605 // Must always clear pPolygon - for clip << 253 return;
606 // early exit << 254 }
607 << 255 else
608 pPolygon.clear(); << 256 {
609 << 257 outputPolygon.clear();
610 if ( outputPolygon.empty() ) return; << 258 }
611 << 259 }
612 G4VoxelLimits simpleLimit2; << 260
613 simpleLimit2.AddLimit(kZAxis,-kInfinity, << 261 if (pVoxelLimit.IsZLimited())
614 ClipPolygonToSimpleLimits(outputPolygon, << 262 {
615 << 263 G4VoxelLimits simpleLimit1;
616 // Return after final clip - no cleanup << 264 simpleLimit1.AddLimit(kZAxis,pVoxelLimit.GetMinZExtent(),kInfinity);
617 } << 265 ClipPolygonToSimpleLimits(pPolygon,outputPolygon,
618 } << 266 simpleLimit1);
>> 267 // Must always clear pPolygon - for clip to simpleLimit2 and incase of
>> 268 // early exit
>> 269 pPolygon.clear();
>> 270 if (!outputPolygon.size())
>> 271 {
>> 272 return;
>> 273 }
>> 274
>> 275 G4VoxelLimits simpleLimit2;
>> 276 simpleLimit2.AddLimit(kZAxis,-kInfinity,pVoxelLimit.GetMaxZExtent());
>> 277 ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
>> 278
>> 279 // Return after final clip - no cleanup
>> 280 }
>> 281 }
619 } 282 }
620 283
621 ////////////////////////////////////////////// <<
622 // <<
623 // pVoxelLimits must be only limited along one 284 // pVoxelLimits must be only limited along one axis, and either the maximum
624 // along the axis must be +kInfinity, or the m 285 // along the axis must be +kInfinity, or the minimum -kInfinity
>> 286 void G4VSolid::ClipPolygonToSimpleLimits(G4ThreeVectorList& pPolygon,
>> 287 G4ThreeVectorList& outputPolygon,
>> 288 const G4VoxelLimits& pVoxelLimit) const
>> 289 {
>> 290 G4int i;
>> 291 G4int noVertices=pPolygon.size();
>> 292 G4ThreeVector vEnd,vStart;
>> 293 for (i=0;i<noVertices;i++)
>> 294 {
>> 295 vStart=pPolygon[i];
>> 296 if (i==noVertices-1)
>> 297 {
>> 298 vEnd=pPolygon[0];
>> 299 }
>> 300 else
>> 301 {
>> 302 vEnd=pPolygon[i+1];
>> 303 }
>> 304
>> 305 if (pVoxelLimit.Inside(vStart))
>> 306 {
>> 307 if (pVoxelLimit.Inside(vEnd))
>> 308 {
>> 309 // vStart and vEnd inside -> output end point
>> 310 outputPolygon.push_back(vEnd);
>> 311 }
>> 312 else
>> 313 {
>> 314 // vStart inside, vEnd outside -> output crossing point
>> 315 pVoxelLimit.ClipToLimits(vStart,vEnd);
>> 316 outputPolygon.push_back(vEnd);
>> 317 }
>> 318
>> 319 }
>> 320 else
>> 321 {
>> 322 if (pVoxelLimit.Inside(vEnd))
>> 323 {
>> 324 // vStart outside, vEnd inside -> output inside section
>> 325 pVoxelLimit.ClipToLimits(vStart,vEnd);
>> 326 outputPolygon.push_back(vStart);
>> 327 outputPolygon.push_back(vEnd);
>> 328 }
>> 329 else
>> 330 // Both point outside -> no output
>> 331 {
>> 332 }
>> 333
>> 334 }
>> 335 }
>> 336
>> 337 }
>> 338
>> 339 const G4VSolid* G4VSolid::GetConstituentSolid(G4int no) const
>> 340 { return 0; }
>> 341
>> 342 G4VSolid* G4VSolid::GetConstituentSolid(G4int no)
>> 343 { return 0; }
625 344
626 void << 345 const G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr() const
627 G4VSolid::ClipPolygonToSimpleLimits( G4ThreeVe << 346 { return 0; }
628 G4ThreeVe <<
629 const G4VoxelLi <<
630 { <<
631 G4int i; <<
632 auto noVertices = (G4int)pPolygon.size(); <<
633 G4ThreeVector vEnd,vStart; <<
634 <<
635 for (i = 0 ; i < noVertices ; ++i ) <<
636 { <<
637 vStart = pPolygon[i]; <<
638 if ( i == noVertices-1 ) vEnd = pPolygo <<
639 else vEnd = pPolygo <<
640 <<
641 if ( pVoxelLimit.Inside(vStart) ) <<
642 { <<
643 if (pVoxelLimit.Inside(vEnd)) <<
644 { <<
645 // vStart and vEnd inside -> output en <<
646 // <<
647 outputPolygon.push_back(vEnd); <<
648 } <<
649 else <<
650 { <<
651 // vStart inside, vEnd outside -> outp <<
652 // <<
653 pVoxelLimit.ClipToLimits(vStart,vEnd); <<
654 outputPolygon.push_back(vEnd); <<
655 } <<
656 } <<
657 else <<
658 { <<
659 if (pVoxelLimit.Inside(vEnd)) <<
660 { <<
661 // vStart outside, vEnd inside -> outp <<
662 // <<
663 pVoxelLimit.ClipToLimits(vStart,vEnd); <<
664 outputPolygon.push_back(vStart); <<
665 outputPolygon.push_back(vEnd); <<
666 } <<
667 else // Both point outside -> no output <<
668 { <<
669 // outputPolygon.push_back(vStart); <<
670 // outputPolygon.push_back(vEnd); <<
671 } <<
672 } <<
673 } <<
674 } <<
675 347
676 ////////////////////////////////////////////// << 348 G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr()
677 // << 349 { return 0; }
678 // Throw exception (warning) for solids not im <<
679 <<
680 void G4VSolid::BoundingLimits(G4ThreeVector& p <<
681 { <<
682 std::ostringstream message; <<
683 message << "Not implemented for solid: " <<
684 << GetEntityType() << " !" <<
685 << "\nReturning infinite boundinx bo <<
686 G4Exception("G4VSolid::BoundingLimits()", "G <<
687 JustWarning, message); <<
688 350
689 pMin.set(-kInfinity,-kInfinity,-kInfinity); << 351 G4VisExtent G4VSolid::GetExtent () const {
690 pMax.set( kInfinity, kInfinity, kInfinity); <<
691 } <<
692 <<
693 ////////////////////////////////////////////// <<
694 // <<
695 // Get G4VisExtent - bounding box for graphics <<
696 <<
697 G4VisExtent G4VSolid::GetExtent () const <<
698 { <<
699 G4VisExtent extent; 352 G4VisExtent extent;
700 G4VoxelLimits voxelLimits; // Defaults to " 353 G4VoxelLimits voxelLimits; // Defaults to "infinite" limits.
701 G4AffineTransform affineTransform; 354 G4AffineTransform affineTransform;
702 G4double vmin, vmax; 355 G4double vmin, vmax;
703 CalculateExtent(kXAxis,voxelLimits,affineTra 356 CalculateExtent(kXAxis,voxelLimits,affineTransform,vmin,vmax);
704 extent.SetXmin (vmin); 357 extent.SetXmin (vmin);
705 extent.SetXmax (vmax); 358 extent.SetXmax (vmax);
706 CalculateExtent(kYAxis,voxelLimits,affineTra 359 CalculateExtent(kYAxis,voxelLimits,affineTransform,vmin,vmax);
707 extent.SetYmin (vmin); 360 extent.SetYmin (vmin);
708 extent.SetYmax (vmax); 361 extent.SetYmax (vmax);
709 CalculateExtent(kZAxis,voxelLimits,affineTra 362 CalculateExtent(kZAxis,voxelLimits,affineTransform,vmin,vmax);
710 extent.SetZmin (vmin); 363 extent.SetZmin (vmin);
711 extent.SetZmax (vmax); 364 extent.SetZmax (vmax);
712 return extent; 365 return extent;
713 } 366 }
714 367
715 G4Polyhedron* G4VSolid::CreatePolyhedron () co 368 G4Polyhedron* G4VSolid::CreatePolyhedron () const
716 { 369 {
717 return nullptr; << 370 return 0;
718 } 371 }
719 372
720 G4Polyhedron* G4VSolid::GetPolyhedron () const << 373 G4NURBS* G4VSolid::CreateNURBS () const
721 { 374 {
722 return nullptr; << 375 return 0;
723 } 376 }
724 377