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