| Geant4 Cross Reference |
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 // G4VCSGfaceted implementation; a virtual class of a CSG type shape
27 // that is built entirely out of G4VCSGface faces.
28 //
29 // Author: David C. Williams (davidw@scipp.ucsc.edu)
30 // --------------------------------------------------------------------
31
32 #include "G4VCSGfaceted.hh"
33 #include "G4VCSGface.hh"
34 #include "G4SolidExtentList.hh"
35
36 #include "G4VoxelLimits.hh"
37 #include "G4AffineTransform.hh"
38
39 #include "Randomize.hh"
40
41 #include "G4Polyhedron.hh"
42 #include "G4VGraphicsScene.hh"
43 #include "G4VisExtent.hh"
44
45 #include "G4AutoLock.hh"
46
47 namespace
48 {
49 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
50 }
51
52 //
53 // Constructor
54 //
55 G4VCSGfaceted::G4VCSGfaceted( const G4String& name )
56 : G4VSolid(name),
57 fStatistics(1000000), fCubVolEpsilon(0.001), fAreaAccuracy(-1.)
58 {
59 }
60
61
62 //
63 // Fake default constructor - sets only member data and allocates memory
64 // for usage restricted to object persistency.
65 //
66 G4VCSGfaceted::G4VCSGfaceted( __void__& a )
67 : G4VSolid(a),
68 fStatistics(1000000), fCubVolEpsilon(0.001), fAreaAccuracy(-1.)
69 {
70 }
71
72 //
73 // Destructor
74 //
75 G4VCSGfaceted::~G4VCSGfaceted()
76 {
77 DeleteStuff();
78 delete fpPolyhedron; fpPolyhedron = nullptr;
79 }
80
81
82 //
83 // Copy constructor
84 //
85 G4VCSGfaceted::G4VCSGfaceted( const G4VCSGfaceted& source )
86 : G4VSolid( source )
87 {
88 fStatistics = source.fStatistics;
89 fCubVolEpsilon = source.fCubVolEpsilon;
90 fAreaAccuracy = source.fAreaAccuracy;
91
92 CopyStuff( source );
93 }
94
95
96 //
97 // Assignment operator
98 //
99 G4VCSGfaceted& G4VCSGfaceted::operator=( const G4VCSGfaceted& source )
100 {
101 if (&source == this) { return *this; }
102
103 // Copy base class data
104 //
105 G4VSolid::operator=(source);
106
107 // Copy data
108 //
109 fStatistics = source.fStatistics;
110 fCubVolEpsilon = source.fCubVolEpsilon;
111 fAreaAccuracy = source.fAreaAccuracy;
112
113 DeleteStuff();
114 CopyStuff( source );
115
116 return *this;
117 }
118
119
120 //
121 // CopyStuff (protected)
122 //
123 // Copy the contents of source
124 //
125 void G4VCSGfaceted::CopyStuff( const G4VCSGfaceted& source )
126 {
127 numFace = source.numFace;
128 if (numFace == 0) { return; } // odd, but permissable?
129
130 faces = new G4VCSGface*[numFace];
131
132 G4VCSGface **face = faces,
133 **sourceFace = source.faces;
134 do // Loop checking, 13.08.2015, G.Cosmo
135 {
136 *face = (*sourceFace)->Clone();
137 } while( ++sourceFace, ++face < faces+numFace );
138 fCubicVolume = source.fCubicVolume;
139 fSurfaceArea = source.fSurfaceArea;
140 fRebuildPolyhedron = false;
141 fpPolyhedron = nullptr;
142 }
143
144
145 //
146 // DeleteStuff (protected)
147 //
148 // Delete all allocated objects
149 //
150 void G4VCSGfaceted::DeleteStuff()
151 {
152 if (numFace != 0)
153 {
154 G4VCSGface **face = faces;
155 do // Loop checking, 13.08.2015, G.Cosmo
156 {
157 delete *face;
158 } while( ++face < faces + numFace );
159
160 delete [] faces;
161 }
162 delete fpPolyhedron; fpPolyhedron = nullptr;
163 }
164
165
166 //
167 // CalculateExtent
168 //
169 G4bool G4VCSGfaceted::CalculateExtent( const EAxis axis,
170 const G4VoxelLimits& voxelLimit,
171 const G4AffineTransform& transform,
172 G4double& min,
173 G4double& max ) const
174 {
175 G4SolidExtentList extentList( axis, voxelLimit );
176
177 //
178 // Loop over all faces, checking min/max extent as we go.
179 //
180 G4VCSGface **face = faces;
181 do // Loop checking, 13.08.2015, G.Cosmo
182 {
183 (*face)->CalculateExtent( axis, voxelLimit, transform, extentList );
184 } while( ++face < faces + numFace );
185
186 //
187 // Return min/max value
188 //
189 return extentList.GetExtent( min, max );
190 }
191
192
193 //
194 // Inside
195 //
196 // It could be a good idea to override this virtual
197 // member to add first a simple test (such as spherical
198 // test or whatnot) and to call this version only if
199 // the simplier test fails.
200 //
201 EInside G4VCSGfaceted::Inside( const G4ThreeVector& p ) const
202 {
203 EInside answer=kOutside;
204 G4VCSGface **face = faces;
205 G4double best = kInfinity;
206 do // Loop checking, 13.08.2015, G.Cosmo
207 {
208 G4double distance;
209 EInside result = (*face)->Inside( p, kCarTolerance/2, &distance );
210 if (result == kSurface) { return kSurface; }
211 if (distance < best)
212 {
213 best = distance;
214 answer = result;
215 }
216 } while( ++face < faces + numFace );
217
218 return answer;
219 }
220
221
222 //
223 // SurfaceNormal
224 //
225 G4ThreeVector G4VCSGfaceted::SurfaceNormal( const G4ThreeVector& p ) const
226 {
227 G4ThreeVector answer;
228 G4VCSGface **face = faces;
229 G4double best = kInfinity;
230 do // Loop checking, 13.08.2015, G.Cosmo
231 {
232 G4double distance = kInfinity;
233 G4ThreeVector normal = (*face)->Normal( p, &distance );
234 if (distance < best)
235 {
236 best = distance;
237 answer = normal;
238 }
239 } while( ++face < faces + numFace );
240
241 return answer;
242 }
243
244
245 //
246 // DistanceToIn(p,v)
247 //
248 G4double G4VCSGfaceted::DistanceToIn( const G4ThreeVector& p,
249 const G4ThreeVector& v ) const
250 {
251 G4double distance = kInfinity;
252 G4double distFromSurface = kInfinity;
253 G4VCSGface **face = faces;
254 G4VCSGface *bestFace = *face;
255 do // Loop checking, 13.08.2015, G.Cosmo
256 {
257 G4double faceDistance,
258 faceDistFromSurface;
259 G4ThreeVector faceNormal;
260 G4bool faceAllBehind;
261 if ((*face)->Intersect( p, v, false, kCarTolerance/2,
262 faceDistance, faceDistFromSurface,
263 faceNormal, faceAllBehind ) )
264 {
265 //
266 // Intersecting face
267 //
268 if (faceDistance < distance)
269 {
270 distance = faceDistance;
271 distFromSurface = faceDistFromSurface;
272 bestFace = *face;
273 if (distFromSurface <= 0) { return 0; }
274 }
275 }
276 } while( ++face < faces + numFace );
277
278 if (distance < kInfinity && distFromSurface<kCarTolerance/2)
279 {
280 if (bestFace->Distance(p,false) < kCarTolerance/2) { distance = 0; }
281 }
282
283 return distance;
284 }
285
286
287 //
288 // DistanceToIn(p)
289 //
290 G4double G4VCSGfaceted::DistanceToIn( const G4ThreeVector& p ) const
291 {
292 return DistanceTo( p, false );
293 }
294
295
296 //
297 // DistanceToOut(p,v)
298 //
299 G4double G4VCSGfaceted::DistanceToOut( const G4ThreeVector& p,
300 const G4ThreeVector& v,
301 const G4bool calcNorm,
302 G4bool* validNorm,
303 G4ThreeVector* n ) const
304 {
305 G4bool allBehind = true;
306 G4double distance = kInfinity;
307 G4double distFromSurface = kInfinity;
308 G4ThreeVector normal;
309
310 G4VCSGface **face = faces;
311 G4VCSGface *bestFace = *face;
312 do // Loop checking, 13.08.2015, G.Cosmo
313 {
314 G4double faceDistance,
315 faceDistFromSurface;
316 G4ThreeVector faceNormal;
317 G4bool faceAllBehind;
318 if ((*face)->Intersect( p, v, true, kCarTolerance/2,
319 faceDistance, faceDistFromSurface,
320 faceNormal, faceAllBehind ) )
321 {
322 //
323 // Intersecting face
324 //
325 if ( (distance < kInfinity) || (!faceAllBehind) ) { allBehind = false; }
326 if (faceDistance < distance)
327 {
328 distance = faceDistance;
329 distFromSurface = faceDistFromSurface;
330 normal = faceNormal;
331 bestFace = *face;
332 if (distFromSurface <= 0.) { break; }
333 }
334 }
335 } while( ++face < faces + numFace );
336
337 if (distance < kInfinity)
338 {
339 if (distFromSurface <= 0.)
340 {
341 distance = 0.;
342 }
343 else if (distFromSurface<kCarTolerance/2)
344 {
345 if (bestFace->Distance(p,true) < kCarTolerance/2) { distance = 0.; }
346 }
347
348 if (calcNorm)
349 {
350 *validNorm = allBehind;
351 *n = normal;
352 }
353 }
354 else
355 {
356 if (Inside(p) == kSurface) { distance = 0.; }
357 if (calcNorm) { *validNorm = false; }
358 }
359
360 return distance;
361 }
362
363
364 //
365 // DistanceToOut(p)
366 //
367 G4double G4VCSGfaceted::DistanceToOut( const G4ThreeVector& p ) const
368 {
369 return DistanceTo( p, true );
370 }
371
372
373 //
374 // DistanceTo
375 //
376 // Protected routine called by DistanceToIn and DistanceToOut
377 //
378 G4double G4VCSGfaceted::DistanceTo( const G4ThreeVector& p,
379 const G4bool outgoing ) const
380 {
381 G4VCSGface **face = faces;
382 G4double best = kInfinity;
383 do // Loop checking, 13.08.2015, G.Cosmo
384 {
385 G4double distance = (*face)->Distance( p, outgoing );
386 if (distance < best) { best = distance; }
387 } while( ++face < faces + numFace );
388
389 return (best < 0.5*kCarTolerance) ? 0. : best;
390 }
391
392
393 //
394 // DescribeYourselfTo
395 //
396 void G4VCSGfaceted::DescribeYourselfTo( G4VGraphicsScene& scene ) const
397 {
398 scene.AddSolid( *this );
399 }
400
401
402 //
403 // GetExtent
404 //
405 // Define the sides of the box into which our solid instance would fit.
406 //
407 G4VisExtent G4VCSGfaceted::GetExtent() const
408 {
409 static const G4ThreeVector xMax(1,0,0), xMin(-1,0,0),
410 yMax(0,1,0), yMin(0,-1,0),
411 zMax(0,0,1), zMin(0,0,-1);
412 static const G4ThreeVector *axes[6] =
413 { &xMin, &xMax, &yMin, &yMax, &zMin, &zMax };
414
415 G4double answers[6] =
416 {-kInfinity, -kInfinity, -kInfinity, -kInfinity, -kInfinity, -kInfinity};
417
418 G4VCSGface **face = faces;
419 do // Loop checking, 13.08.2015, G.Cosmo
420 {
421 const G4ThreeVector **axis = axes+5 ;
422 G4double* answer = answers+5;
423 do // Loop checking, 13.08.2015, G.Cosmo
424 {
425 G4double testFace = (*face)->Extent( **axis );
426 if (testFace > *answer) { *answer = testFace; }
427 }
428 while( --axis, --answer >= answers );
429
430 } while( ++face < faces + numFace );
431
432 return { -answers[0], answers[1],
433 -answers[2], answers[3],
434 -answers[4], answers[5] };
435 }
436
437
438 //
439 // GetEntityType
440 //
441 G4GeometryType G4VCSGfaceted::GetEntityType() const
442 {
443 return {"G4CSGfaceted"};
444 }
445
446
447 //
448 // Stream object contents to an output stream
449 //
450 std::ostream& G4VCSGfaceted::StreamInfo( std::ostream& os ) const
451 {
452 os << "-----------------------------------------------------------\n"
453 << " *** Dump for solid - " << GetName() << " ***\n"
454 << " ===================================================\n"
455 << " Solid type: G4VCSGfaceted\n"
456 << " Parameters: \n"
457 << " number of faces: " << numFace << "\n"
458 << "-----------------------------------------------------------\n";
459
460 return os;
461 }
462
463
464 //
465 // GetCubVolStatistics
466 //
467 G4int G4VCSGfaceted::GetCubVolStatistics() const
468 {
469 return fStatistics;
470 }
471
472
473 //
474 // GetCubVolEpsilon
475 //
476 G4double G4VCSGfaceted::GetCubVolEpsilon() const
477 {
478 return fCubVolEpsilon;
479 }
480
481
482 //
483 // SetCubVolStatistics
484 //
485 void G4VCSGfaceted::SetCubVolStatistics(G4int st)
486 {
487 fCubicVolume=0.;
488 fStatistics=st;
489 }
490
491
492 //
493 // SetCubVolEpsilon
494 //
495 void G4VCSGfaceted::SetCubVolEpsilon(G4double ep)
496 {
497 fCubicVolume=0.;
498 fCubVolEpsilon=ep;
499 }
500
501
502 //
503 // GetAreaStatistics
504 //
505 G4int G4VCSGfaceted::GetAreaStatistics() const
506 {
507 return fStatistics;
508 }
509
510
511 //
512 // GetAreaAccuracy
513 //
514 G4double G4VCSGfaceted::GetAreaAccuracy() const
515 {
516 return fAreaAccuracy;
517 }
518
519
520 //
521 // SetAreaStatistics
522 //
523 void G4VCSGfaceted::SetAreaStatistics(G4int st)
524 {
525 fSurfaceArea=0.;
526 fStatistics=st;
527 }
528
529
530 //
531 // SetAreaAccuracy
532 //
533 void G4VCSGfaceted::SetAreaAccuracy(G4double ep)
534 {
535 fSurfaceArea=0.;
536 fAreaAccuracy=ep;
537 }
538
539
540 //
541 // GetCubicVolume
542 //
543 G4double G4VCSGfaceted::GetCubicVolume()
544 {
545 if(fCubicVolume != 0.) {;}
546 else { fCubicVolume = EstimateCubicVolume(fStatistics,fCubVolEpsilon); }
547 return fCubicVolume;
548 }
549
550
551 //
552 // GetSurfaceArea
553 //
554 G4double G4VCSGfaceted::GetSurfaceArea()
555 {
556 if(fSurfaceArea != 0.) {;}
557 else { fSurfaceArea = EstimateSurfaceArea(fStatistics,fAreaAccuracy); }
558 return fSurfaceArea;
559 }
560
561
562 //
563 // GetPolyhedron
564 //
565 G4Polyhedron* G4VCSGfaceted::GetPolyhedron () const
566 {
567 if (fpPolyhedron == nullptr ||
568 fRebuildPolyhedron ||
569 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
570 fpPolyhedron->GetNumberOfRotationSteps())
571 {
572 G4AutoLock l(&polyhedronMutex);
573 delete fpPolyhedron;
574 fpPolyhedron = CreatePolyhedron();
575 fRebuildPolyhedron = false;
576 l.unlock();
577 }
578 return fpPolyhedron;
579 }
580
581
582 //
583 // GetPointOnSurfaceGeneric proportional to Areas of faces
584 // in case of GenericPolycone or GenericPolyhedra
585 //
586 G4ThreeVector G4VCSGfaceted::GetPointOnSurfaceGeneric( ) const
587 {
588 // Preparing variables
589 //
590 G4ThreeVector answer=G4ThreeVector(0.,0.,0.);
591 G4VCSGface **face = faces;
592 G4double area = 0.;
593 G4int i;
594 std::vector<G4double> areas;
595
596 // First step: calculate surface areas
597 //
598 do // Loop checking, 13.08.2015, G.Cosmo
599 {
600 G4double result = (*face)->SurfaceArea( );
601 areas.push_back(result);
602 area=area+result;
603 } while( ++face < faces + numFace );
604
605 // Second Step: choose randomly one surface
606 //
607 G4VCSGface **face1 = faces;
608 G4double chose = area*G4UniformRand();
609 G4double Achose1, Achose2;
610 Achose1=0.; Achose2=0.;
611 i=0;
612
613 do
614 {
615 Achose2+=areas[i];
616 if(chose>=Achose1 && chose<Achose2)
617 {
618 G4ThreeVector point;
619 point= (*face1)->GetPointOnFace();
620 return point;
621 }
622 ++i;
623 Achose1=Achose2;
624 } while( ++face1 < faces + numFace );
625
626 return answer;
627 }
628