Geant4  10.02.p01
G4USolid.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 // GEANT4 tag $Name:$
29 //
30 //
31 // G4USolid implementation
32 //
33 // --------------------------------------------------------------------
34 
35 #include "G4USolid.hh"
36 
37 #if defined(G4GEOM_USE_USOLIDS)
38 
39 #include "G4AffineTransform.hh"
40 #include "G4VoxelLimits.hh"
41 #include "G4VGraphicsScene.hh"
42 #include "G4Polyhedron.hh"
43 #include "G4PolyhedronArbitrary.hh"
44 #include "G4VisExtent.hh"
45 #include "G4PhysicalConstants.hh"
46 
47 #include "G4AutoLock.hh"
48 
49 namespace
50 {
51  G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
52 }
53 
54 G4USolid::G4USolid(const G4String& name, VUSolid* s) :
55  G4VSolid(name), fShape(s), fRebuildPolyhedron(false), fPolyhedron(0)
56 {
57 }
58 
59 G4USolid::G4USolid(__void__& a)
60  : G4VSolid(a), fShape(0), fRebuildPolyhedron(false), fPolyhedron(0)
61 {
62 }
63 
64 G4USolid::~G4USolid()
65 {
66  delete fPolyhedron; fPolyhedron = 0;
67 }
68 
69 G4bool G4USolid::operator==(const G4USolid& s) const
70 {
71  return (this == &s) ? true : false;
72 }
73 
74 EInside G4USolid::Inside(const G4ThreeVector& p) const
75 {
76  UVector3 pt;
77  VUSolid::EnumInside in_temp;
78  EInside in = kOutside;
79  pt.x() = p.x();
80  pt.y() = p.y();
81  pt.z() = p.z(); // better assign at construction
82 
83  in_temp = fShape->Inside(pt);
84 
85 #ifndef G4USE_STD11
86  if (in_temp == VUSolid::eSurface)return kSurface;
87  if (in_temp == VUSolid::eInside)return kInside;
88 #else
89  if (in_temp == VUSolid::EnumInside::eSurface)return kSurface;
90  if (in_temp == VUSolid::EnumInside::eInside)return kInside;
91 #endif
92 
93  return in;
94 }
95 
96 G4ThreeVector G4USolid::SurfaceNormal(const G4ThreeVector& pt) const
97 {
98  UVector3 p;
99  p.x() = pt.x();
100  p.y() = pt.y();
101  p.z() = pt.z();
102  UVector3 n;
103  fShape->Normal(p, n);
104  return G4ThreeVector(n.x(), n.y(), n.z());
105 }
106 
107 G4double G4USolid::DistanceToIn(const G4ThreeVector& pt,
108  const G4ThreeVector& d)const
109 {
110  UVector3 p;
111  p.x() = pt.x();
112  p.y() = pt.y();
113  p.z() = pt.z(); // better assign at construction
114  UVector3 v;
115  v.x() = d.x();
116  v.y() = d.y();
117  v.z() = d.z(); // better assign at construction
118  G4double dist = fShape->DistanceToIn(p, v);
119  if (dist > kInfinity) dist = kInfinity;
120  return dist;
121 }
122 
123 G4double G4USolid::DistanceToIn(const G4ThreeVector& pt) const
124 {
125  UVector3 p;
126  p.x() = pt.x();
127  p.y() = pt.y();
128  p.z() = pt.z(); // better assign at construction
129  G4double dist = fShape->SafetyFromOutside(p); // true?
130  if (dist > kInfinity) dist = kInfinity;
131  return dist;
132 }
133 
134 G4double G4USolid::DistanceToOut(const G4ThreeVector& pt,
135  const G4ThreeVector& d,
136  const G4bool calcNorm,
137  G4bool* validNorm,
138  G4ThreeVector* norm) const
139 {
140  UVector3 p;
141  p.x() = pt.x();
142  p.y() = pt.y();
143  p.z() = pt.z(); // better assign at construction
144  UVector3 v;
145  v.x() = d.x();
146  v.y() = d.y();
147  v.z() = d.z(); // better assign at construction
148  UVector3 n;
149  bool valid;
150  G4double dist = fShape->DistanceToOut(p, v, n,valid); // should use local variable
151  if(calcNorm)
152  {
153  if(valid){ *validNorm = true;}
154  else {* validNorm =false;}
155  if(*validNorm)
156  { norm->setX(n.x());
157  norm->setY(n.y());
158  norm->setZ(n.z());
159  } // *norm = n, but only after calcNorm check
160  }
161  if (dist > kInfinity) dist = kInfinity;
162  return dist;
163 }
164 
165 G4double G4USolid::DistanceToOut(const G4ThreeVector& pt) const
166 {
167  UVector3 p;
168  p.x() = pt.x();
169  p.y() = pt.y();
170  p.z() = pt.z(); // better assign at construction
171  return fShape->SafetyFromInside(p); // true?
172 }
173 
174 G4double G4USolid::GetCubicVolume()
175 {
176  return fShape->Capacity();
177 }
178 
179 G4double G4USolid::GetSurfaceArea()
180 {
181  return fShape->SurfaceArea();
182 }
183 
184 G4ThreeVector G4USolid::GetPointOnSurface() const
185 {
186  UVector3 p;
187  p = fShape->GetPointOnSurface();
188  return G4ThreeVector(p.x(), p.y(), p.z());
189 }
190 
191 G4bool G4USolid::CalculateExtent(const EAxis pAxis,
192  const G4VoxelLimits& pVoxelLimit,
193  const G4AffineTransform& pTransform,
194  G4double& pMin, G4double& pMax) const
195 {
196  if (!pTransform.IsRotated())
197  {
198  VUSolid::EAxisType eAxis = VUSolid::eXaxis;
199  G4double offset = pTransform.NetTranslation().x();
200  if (pAxis == kYAxis)
201  {
202  eAxis = VUSolid::eYaxis;
203  offset = pTransform.NetTranslation().y();
204  }
205  if (pAxis == kZAxis)
206  {
207  eAxis = VUSolid::eZaxis;
208  offset = pTransform.NetTranslation().z();
209  }
210  fShape->ExtentAxis(eAxis, pMin, pMax);
211 
212  pMin += offset;
213  pMax += offset;
214 
215  if (pVoxelLimit.IsLimited())
216  {
217  switch (pAxis)
218  {
219  case kXAxis:
220  if ((pMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance) ||
221  (pMax < pVoxelLimit.GetMinXExtent() - kCarTolerance))
222  {
223  return false;
224  }
225  else
226  {
227  pMin = std::max(pMin, pVoxelLimit.GetMinXExtent());
228  pMax = std::min(pMax, pVoxelLimit.GetMaxXExtent());
229  }
230  break;
231  case kYAxis:
232  if ((pMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance) ||
233  (pMax < pVoxelLimit.GetMinYExtent() - kCarTolerance))
234  {
235  return false;
236  }
237  else
238  {
239  pMin = std::max(pMin, pVoxelLimit.GetMinYExtent());
240  pMax = std::min(pMax, pVoxelLimit.GetMaxYExtent());
241  }
242  break;
243  case kZAxis:
244  if ((pMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance) ||
245  (pMax < pVoxelLimit.GetMinZExtent() - kCarTolerance))
246  {
247  return false;
248  }
249  else
250  {
251  pMin = std::max(pMin, pVoxelLimit.GetMinZExtent());
252  pMax = std::min(pMax, pVoxelLimit.GetMaxZExtent());
253  }
254  break;
255  default:
256  break;
257  }
258  pMin -= kCarTolerance ;
259  pMax += kCarTolerance ;
260  }
261  return true;
262  }
263  else // General rotated case - create and clip mesh to boundaries
264  {
265  // Rotate BoundingBox and Calculate Extent as for BREPS
266 
267  G4bool existsAfterClip = false ;
268  G4ThreeVectorList* vertices ;
269 
270  pMin = +kInfinity ;
271  pMax = -kInfinity ;
272 
273  // Calculate rotated vertex coordinates
274 
275  vertices = CreateRotatedVertices(pTransform) ;
276  ClipCrossSection(vertices, 0, pVoxelLimit, pAxis, pMin, pMax) ;
277  ClipCrossSection(vertices, 4, pVoxelLimit, pAxis, pMin, pMax) ;
278  ClipBetweenSections(vertices, 0, pVoxelLimit, pAxis, pMin, pMax) ;
279 
280  if (pVoxelLimit.IsLimited(pAxis) == false)
281  {
282  if ((pMin != kInfinity) || (pMax != -kInfinity))
283  {
284  existsAfterClip = true ;
285 
286  // Add 2*tolerance to avoid precision troubles
287 
288  pMin -= kCarTolerance;
289  pMax += kCarTolerance;
290  }
291  }
292  else
293  {
294  G4ThreeVector clipCentre(
295  (pVoxelLimit.GetMinXExtent() + pVoxelLimit.GetMaxXExtent()) * 0.5,
296  (pVoxelLimit.GetMinYExtent() + pVoxelLimit.GetMaxYExtent()) * 0.5,
297  (pVoxelLimit.GetMinZExtent() + pVoxelLimit.GetMaxZExtent()) * 0.5);
298 
299  if ((pMin != kInfinity) || (pMax != -kInfinity))
300  {
301  existsAfterClip = true ;
302 
303 
304  // Check to see if endpoints are in the solid
305 
306  clipCentre(pAxis) = pVoxelLimit.GetMinExtent(pAxis);
307 
308  if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside)
309  {
310  pMin = pVoxelLimit.GetMinExtent(pAxis);
311  }
312  else
313  {
314  pMin -= kCarTolerance;
315  }
316  clipCentre(pAxis) = pVoxelLimit.GetMaxExtent(pAxis);
317 
318  if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside)
319  {
320  pMax = pVoxelLimit.GetMaxExtent(pAxis);
321  }
322  else
323  {
324  pMax += kCarTolerance;
325  }
326  }
327 
328  // Check for case where completely enveloping clipping volume
329  // If point inside then we are confident that the solid completely
330  // envelopes the clipping volume. Hence set min/max extents according
331  // to clipping volume extents along the specified axis.
332 
333  else if (Inside(pTransform.Inverse().TransformPoint(clipCentre))
334  != kOutside)
335  {
336  existsAfterClip = true ;
337  pMin = pVoxelLimit.GetMinExtent(pAxis) ;
338  pMax = pVoxelLimit.GetMaxExtent(pAxis) ;
339  }
340  }
341  delete vertices;
342  return existsAfterClip;
343  }
344 }
345 
346 void G4USolid::ComputeDimensions(G4VPVParameterisation*,
347  const G4int,
348  const G4VPhysicalVolume*)
349 {
350  std::ostringstream message;
351  message << "Illegal call to G4USolid::ComputeDimensions()" << G4endl
352  << "Method not overloaded by derived class !";
353  G4Exception("G4USolid::ComputeDimensions()", "GeomSolids0003",
354  FatalException, message);
355 }
356 
357 void G4USolid::DescribeYourselfTo(G4VGraphicsScene& scene) const
358 {
359  scene.AddSolid(*this);
360 }
361 
362 G4GeometryType G4USolid::GetEntityType() const
363 {
364 
365  G4String string = fShape->GetEntityType();
366  return "G4" + string;
367 }
368 
369 std::ostream& G4USolid::StreamInfo(std::ostream& os) const
370 {
371  return fShape->StreamInfo(os);
372 }
373 
374 G4USolid::G4USolid(const G4USolid& rhs)
375  : G4VSolid(rhs), fRebuildPolyhedron(false), fPolyhedron(0)
376 {
377  fShape = rhs.fShape->Clone();
378 }
379 
380 G4USolid& G4USolid::operator=(const G4USolid& rhs)
381 {
382  // Check assignment to self
383  //
384  if (this == &rhs)
385  {
386  return *this;
387  }
388 
389  // Copy base class data
390  //
391  G4VSolid::operator=(rhs);
392 
393  // Copy data
394  //
395  fShape = rhs.fShape->Clone();
396  fRebuildPolyhedron = false;
397  delete fPolyhedron; fPolyhedron = 0;
398 
399  return *this;
400 }
401 
402 G4VSolid* G4USolid::Clone() const
403 {
404  std::ostringstream message;
405  message << "Clone() method not implemented for type: "
406  << GetEntityType() << "!" << G4endl
407  << "Returning NULL pointer!";
408  G4Exception("G4USolid::Clone()", "GeomSolids1001", JustWarning, message);
409  return 0;
410 }
411 
413 G4USolid::CreateRotatedVertices(const G4AffineTransform& pTransform) const
414 {
415  G4double xMin, xMax, yMin, yMax, zMin, zMax;
416 
417  fShape->ExtentAxis(VUSolid::eXaxis, xMin, xMax);
418  fShape->ExtentAxis(VUSolid::eYaxis, yMin, yMax);
419  fShape->ExtentAxis(VUSolid::eZaxis, zMin, zMax);
420 
421  G4ThreeVectorList* vertices;
422  vertices = new G4ThreeVectorList();
423 
424  if (vertices)
425  {
426  vertices->reserve(8);
427  G4ThreeVector vertex0(xMin, yMin, zMin);
428  G4ThreeVector vertex1(xMax, yMin, zMin);
429  G4ThreeVector vertex2(xMax, yMax, zMin);
430  G4ThreeVector vertex3(xMin, yMax, zMin);
431  G4ThreeVector vertex4(xMin, yMin, zMax);
432  G4ThreeVector vertex5(xMax, yMin, zMax);
433  G4ThreeVector vertex6(xMax, yMax, zMax);
434  G4ThreeVector vertex7(xMin, yMax, zMax);
435 
436  vertices->push_back(pTransform.TransformPoint(vertex0));
437  vertices->push_back(pTransform.TransformPoint(vertex1));
438  vertices->push_back(pTransform.TransformPoint(vertex2));
439  vertices->push_back(pTransform.TransformPoint(vertex3));
440  vertices->push_back(pTransform.TransformPoint(vertex4));
441  vertices->push_back(pTransform.TransformPoint(vertex5));
442  vertices->push_back(pTransform.TransformPoint(vertex6));
443  vertices->push_back(pTransform.TransformPoint(vertex7));
444  }
445  else
446  {
447  G4Exception("G4VUSolid::CreateRotatedVertices()", "FatalError",
448  FatalException, "Out of memory - Cannot allocate vertices!");
449  }
450  return vertices;
451 }
452 
453 G4Polyhedron* G4USolid::CreatePolyhedron() const
454 {
455  G4int index = 0;
456  if (fShape->GetEntityType() == "Box")
457  {
458  double array[3];
459  fShape->GetParametersList(index, array);
460  return new G4PolyhedronBox(array[0], array[1], array[2]);
461  }
462  if (fShape->GetEntityType() == "Tubs")
463  {
464  double array[5];
465  fShape->GetParametersList(index, array);
466  return new G4PolyhedronTubs(array[0], array[1], array[2], array[3], array[4]);
467  }
468  if (fShape->GetEntityType() == "Cons")
469  {
470  double array[7];
471  fShape->GetParametersList(index, array);
472  return new G4PolyhedronCons(array[0], array[1], array[2], array[3], array[4], array[5], array[6]);
473  }
474  if (fShape->GetEntityType() == "Orb")
475  {
476  double array[1];
477  fShape->GetParametersList(index, array);
478  return new G4PolyhedronSphere(0., array[0], 0., 2 * pi, 0., pi);
479  }
480  if (fShape->GetEntityType() == "Sphere")
481  {
482  double array[6];
483  fShape->GetParametersList(index, array);
484  return new G4PolyhedronSphere(array[0], array[1], array[2], array[3], array[4], array[5]);
485  }
486  if (fShape->GetEntityType() == "Tet")
487  {
488  double array[12];
489  fShape->GetParametersList(index, array);
490  G4Polyhedron* ph = new G4Polyhedron;
491  double xyz[4][3];
492  static int faces[4][4] = {{1, 3, 2, 0}, {1, 4, 3, 0}, {1, 2, 4, 0}, {2, 3, 4, 0}};
493  xyz[0][0] = array[0];
494  xyz[0][1] = array[1];
495  xyz[0][2] = array[2];
496  xyz[1][0] = array[3];
497  xyz[1][1] = array[4];
498  xyz[1][2] = array[5];
499  xyz[2][0] = array[6];
500  xyz[2][1] = array[7];
501  xyz[2][2] = array[8];
502  xyz[3][0] = array[9];
503  xyz[3][1] = array[10];
504  xyz[3][2] = array[11];
505 
506  ph->createPolyhedron(4, 4, xyz, faces);
507  return ph;
508  }
509  if (fShape->GetEntityType() == "Trd")
510  {
511  double array[5];
512  fShape->GetParametersList(index, array);
513  return new G4PolyhedronTrd2(array[0], array[1], array[2], array[3], array[4]);
514  }
515  if (fShape->GetEntityType() == "Trap")
516  {
517  double array[12];
518  fShape->GetParametersList(index, array);
519  double phi = (array[11] != 1.0) ? (std::atan(array[10] / array[9])) : (0.0);
520  double alpha1 = std::atan(array[4]);
521  double alpha2 = std::atan(array[8]);
522  double theta = std::acos(array[11]);
523 
524  return new G4PolyhedronTrap(array[0], theta, phi,
525  array[1], array[2], array[3], alpha1,
526  array[5], array[6], array[7], alpha2);
527  }
528 
529  /*
530  if(fShape->GetEntityType()=="TessellatedSolid"){
531 
532  G4Polyhedron *uPolyhedron=fShape->GetPolyhedron();
533  std::size_t nVertices = (*uPolyhedron).vertices.size();
534  std::size_t nFacets = (*uPolyhedron).facets.size();
535 
536  G4PolyhedronArbitrary *polyhedron =
537  new G4PolyhedronArbitrary (nVertices, nFacets);
538 
539  for (std::vector<UVector3>::const_iterator v = (*uPolyhedron).vertices.begin();
540  v!=(*uPolyhedron).vertices.end(); v++)
541  {
542  UVector3 p=(*v);
543  G4ThreeVector pt(p.x(),p.y(),p.z());
544 
545  polyhedron->AddVertex(pt);
546  }
547  for (std::vector<UFacet>::const_iterator f=(*uPolyhedron).facets.begin();
548  f != (*uPolyhedron).facets.end(); f++)
549  {
550  polyhedron->AddFacet((*f).f1,(*f).f2,(*f).f3,(*f).f4);
551  }
552 
553  return (G4Polyhedron*) polyhedron;
554  }
555  */
556 
557  return 0;
558 }
559 
560 G4Polyhedron* G4USolid::GetPolyhedron() const
561 {
562  if (!fPolyhedron ||
563  fRebuildPolyhedron ||
564  fPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
565  fPolyhedron->GetNumberOfRotationSteps())
566  {
567  G4AutoLock l(&polyhedronMutex);
568  delete fPolyhedron;
569  fPolyhedron = CreatePolyhedron();
570  fRebuildPolyhedron = false;
571  l.unlock();
572  }
573  return fPolyhedron;
574 }
575 
576 G4VisExtent G4USolid:: GetExtent() const
577 {
578  G4VisExtent extent;
579  G4VoxelLimits voxelLimits; // Defaults to "infinite" limits.
580  G4AffineTransform affineTransform;
581  G4double vmin, vmax;
582  CalculateExtent(kXAxis, voxelLimits, affineTransform, vmin, vmax);
583  extent.SetXmin(vmin);
584  extent.SetXmax(vmax);
585  CalculateExtent(kYAxis, voxelLimits, affineTransform, vmin, vmax);
586  extent.SetYmin(vmin);
587  extent.SetYmax(vmax);
588  CalculateExtent(kZAxis, voxelLimits, affineTransform, vmin, vmax);
589  extent.SetZmin(vmin);
590  extent.SetZmax(vmax);
591  return extent;
592 }
593 
594 #endif // G4GEOM_USE_USOLIDS
static const G4double kInfinity
Definition: geomdefs.hh:42
G4double GetMinYExtent() const
CLHEP::Hep3Vector G4ThreeVector
G4AffineTransform Inverse() const
G4String name
Definition: TRTMaterials.hh:40
G4bool IsRotated() const
void SetXmax(G4double xmax)
Definition: G4VisExtent.hh:102
G4ThreeVector NetTranslation() const
G4double a
Definition: TRTMaterials.hh:39
virtual void AddSolid(const G4Box &)=0
int G4int
Definition: G4Types.hh:78
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
void SetYmax(G4double ymax)
Definition: G4VisExtent.hh:106
static const double s
Definition: G4SIunits.hh:168
G4double GetMaxXExtent() const
G4double GetMinZExtent() const
G4bool IsLimited() const
bool G4bool
Definition: G4Types.hh:79
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:79
const G4int n
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double GetMinXExtent() const
G4int G4Mutex
Definition: G4Threading.hh:173
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4double GetMaxZExtent() const
void SetZmin(G4double zmin)
Definition: G4VisExtent.hh:108
static const double pi
Definition: G4SIunits.hh:74
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
G4bool operator==(const G4GeometryCell &k1, const G4GeometryCell &k2)
#define G4endl
Definition: G4ios.hh:61
G4double GetMaxYExtent() const
G4VSolid & operator=(const G4VSolid &rhs)
Definition: G4VSolid.cc:110
virtual G4VSolid * Clone() const
Definition: G4VSolid.cc:324
void SetZmax(G4double zmax)
Definition: G4VisExtent.hh:110
void SetXmin(G4double xmin)
Definition: G4VisExtent.hh:100
double G4double
Definition: G4Types.hh:76
G4double GetMaxExtent(const EAxis pAxis) const
void SetYmin(G4double ymin)
Definition: G4VisExtent.hh:104
G4double GetMinExtent(const EAxis pAxis) const