Geant4  10.02.p01
G4UnionSolid.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: G4UnionSolid.cc 95422 2016-02-10 14:35:47Z gcosmo $
28 //
29 // Implementation of methods for the class G4IntersectionSolid
30 //
31 // History:
32 //
33 // 12.09.98 V.Grichine: first implementation
34 // 28.11.98 V.Grichine: fix while loops in DistToIn/Out
35 // 27.07.99 V.Grichine: modifications in DistToOut(p,v,...), while -> do-while
36 // 16.03.01 V.Grichine: modifications in CalculateExtent()
37 //
38 // --------------------------------------------------------------------
39 
40 #include <sstream>
41 
42 #include "G4UnionSolid.hh"
43 
44 #include "G4SystemOfUnits.hh"
45 #include "G4VoxelLimits.hh"
46 #include "G4VPVParameterisation.hh"
47 #include "G4GeometryTolerance.hh"
48 
49 #include "G4VGraphicsScene.hh"
50 #include "G4Polyhedron.hh"
51 #include "HepPolyhedronProcessor.h"
52 
54 //
55 // Transfer all data members to G4BooleanSolid which is responsible
56 // for them. pName will be in turn sent to G4VSolid
57 
59  G4VSolid* pSolidA ,
60  G4VSolid* pSolidB )
61  : G4BooleanSolid(pName,pSolidA,pSolidB)
62 {
63 }
64 
66 //
67 // Constructor
68 
70  G4VSolid* pSolidA ,
71  G4VSolid* pSolidB ,
72  G4RotationMatrix* rotMatrix,
73  const G4ThreeVector& transVector )
74  : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
75 
76 {
77 }
78 
80 //
81 // Constructor
82 
84  G4VSolid* pSolidA ,
85  G4VSolid* pSolidB ,
86  const G4Transform3D& transform )
87  : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
88 {
89 }
90 
92 //
93 // Fake default constructor - sets only member data and allocates memory
94 // for usage restricted to object persistency.
95 
97  : G4BooleanSolid(a)
98 {
99 }
100 
102 //
103 // Destructor
104 
106 {
107 }
108 
110 //
111 // Copy constructor
112 
114  : G4BooleanSolid (rhs)
115 {
116 }
117 
119 //
120 // Assignment operator
121 
123 {
124  // Check assignment to self
125  //
126  if (this == &rhs) { return *this; }
127 
128  // Copy base class data
129  //
131 
132  return *this;
133 }
134 
136 //
137 //
138 
139 G4bool
141  const G4VoxelLimits& pVoxelLimit,
142  const G4AffineTransform& pTransform,
143  G4double& pMin,
144  G4double& pMax ) const
145 {
146  G4bool touchesA, touchesB, out ;
147  G4double minA = kInfinity, minB = kInfinity,
148  maxA = -kInfinity, maxB = -kInfinity;
149 
150  touchesA = fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
151  pTransform, minA, maxA);
152  touchesB= fPtrSolidB->CalculateExtent( pAxis, pVoxelLimit,
153  pTransform, minB, maxB);
154  if( touchesA || touchesB )
155  {
156  pMin = std::min( minA, minB );
157  pMax = std::max( maxA, maxB );
158  out = true ;
159  }
160  else out = false ;
161 
162  return out ; // It exists in this slice if either one does.
163 }
164 
166 //
167 // Important comment: When solids A and B touch together along flat
168 // surface the surface points will be considered as kSurface, while points
169 // located around will correspond to kInside
170 
172 {
173  EInside positionA = fPtrSolidA->Inside(p);
174  if (positionA == kInside) { return kInside; }
175 
176  static const G4double rtol
178  EInside positionB = fPtrSolidB->Inside(p);
179 
180  if( positionB == kInside ||
181  ( positionA == kSurface && positionB == kSurface &&
182  ( fPtrSolidA->SurfaceNormal(p) +
183  fPtrSolidB->SurfaceNormal(p) ).mag2() < rtol ) )
184  {
185  return kInside;
186  }
187  else
188  {
189  if( ( positionB == kSurface ) || ( positionA == kSurface ) )
190  { return kSurface; }
191  else
192  { return kOutside; }
193  }
194 }
195 
197 //
198 //
199 
202 {
204 
205 #ifdef G4BOOLDEBUG
206  if( Inside(p) == kOutside )
207  {
208  G4cout << "WARNING - Invalid call in "
209  << "G4UnionSolid::SurfaceNormal(p)" << G4endl
210  << " Point p is outside !" << G4endl;
211  G4cout << " p = " << p << G4endl;
212  G4cerr << "WARNING - Invalid call in "
213  << "G4UnionSolid::SurfaceNormal(p)" << G4endl
214  << " Point p is outside !" << G4endl;
215  G4cerr << " p = " << p << G4endl;
216  }
217 #endif
218 
219  if(fPtrSolidA->Inside(p) == kSurface && fPtrSolidB->Inside(p) != kInside)
220  {
221  normal= fPtrSolidA->SurfaceNormal(p) ;
222  }
223  else if(fPtrSolidB->Inside(p) == kSurface &&
224  fPtrSolidA->Inside(p) != kInside)
225  {
226  normal= fPtrSolidB->SurfaceNormal(p) ;
227  }
228  else
229  {
230  normal= fPtrSolidA->SurfaceNormal(p) ;
231 #ifdef G4BOOLDEBUG
232  if(Inside(p)==kInside)
233  {
234  G4cout << "WARNING - Invalid call in "
235  << "G4UnionSolid::SurfaceNormal(p)" << G4endl
236  << " Point p is inside !" << G4endl;
237  G4cout << " p = " << p << G4endl;
238  G4cerr << "WARNING - Invalid call in "
239  << "G4UnionSolid::SurfaceNormal(p)" << G4endl
240  << " Point p is inside !" << G4endl;
241  G4cerr << " p = " << p << G4endl;
242  }
243 #endif
244  }
245  return normal;
246 }
247 
249 //
250 // The same algorithm as in DistanceToIn(p)
251 
252 G4double
254  const G4ThreeVector& v ) const
255 {
256 #ifdef G4BOOLDEBUG
257  if( Inside(p) == kInside )
258  {
259  G4cout << "WARNING - Invalid call in "
260  << "G4UnionSolid::DistanceToIn(p,v)" << G4endl
261  << " Point p is inside !" << G4endl;
262  G4cout << " p = " << p << G4endl;
263  G4cout << " v = " << v << G4endl;
264  G4cerr << "WARNING - Invalid call in "
265  << "G4UnionSolid::DistanceToIn(p,v)" << G4endl
266  << " Point p is inside !" << G4endl;
267  G4cerr << " p = " << p << G4endl;
268  G4cerr << " v = " << v << G4endl;
269  }
270 #endif
271 
272  return std::min(fPtrSolidA->DistanceToIn(p,v),
273  fPtrSolidB->DistanceToIn(p,v) ) ;
274 }
275 
277 //
278 // Approximate nearest distance from the point p to the union of
279 // two solids
280 
281 G4double
283 {
284 #ifdef G4BOOLDEBUG
285  if( Inside(p) == kInside )
286  {
287  G4cout << "WARNING - Invalid call in "
288  << "G4UnionSolid::DistanceToIn(p)" << G4endl
289  << " Point p is inside !" << G4endl;
290  G4cout << " p = " << p << G4endl;
291  G4cerr << "WARNING - Invalid call in "
292  << "G4UnionSolid::DistanceToIn(p)" << G4endl
293  << " Point p is inside !" << G4endl;
294  G4cerr << " p = " << p << G4endl;
295  }
296 #endif
297  G4double distA = fPtrSolidA->DistanceToIn(p) ;
298  G4double distB = fPtrSolidB->DistanceToIn(p) ;
299  G4double safety = std::min(distA,distB) ;
300  if(safety < 0.0) safety = 0.0 ;
301  return safety ;
302 }
303 
305 //
306 // The same algorithm as DistanceToOut(p)
307 
308 G4double
310  const G4ThreeVector& v,
311  const G4bool calcNorm,
312  G4bool *validNorm,
313  G4ThreeVector *n ) const
314 {
315  G4double dist = 0.0, disTmp = 0.0 ;
316  G4ThreeVector normTmp;
317  G4ThreeVector* nTmp= &normTmp;
318 
319  if( Inside(p) == kOutside )
320  {
321 #ifdef G4BOOLDEBUG
322  G4cout << "Position:" << G4endl << G4endl;
323  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
324  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
325  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
326  G4cout << "Direction:" << G4endl << G4endl;
327  G4cout << "v.x() = " << v.x() << G4endl;
328  G4cout << "v.y() = " << v.y() << G4endl;
329  G4cout << "v.z() = " << v.z() << G4endl << G4endl;
330  G4cout << "WARNING - Invalid call in "
331  << "G4UnionSolid::DistanceToOut(p,v)" << G4endl
332  << " Point p is outside !" << G4endl;
333  G4cout << " p = " << p << G4endl;
334  G4cout << " v = " << v << G4endl;
335  G4cerr << "WARNING - Invalid call in "
336  << "G4UnionSolid::DistanceToOut(p,v)" << G4endl
337  << " Point p is outside !" << G4endl;
338  G4cerr << " p = " << p << G4endl;
339  G4cerr << " v = " << v << G4endl;
340 #endif
341  }
342  else
343  {
344  EInside positionA = fPtrSolidA->Inside(p) ;
345  // EInside positionB = fPtrSolidB->Inside(p) ;
346 
347  if( positionA != kOutside )
348  {
349  do // Loop checking, 13.08.2015, G.Cosmo
350  {
351  disTmp = fPtrSolidA->DistanceToOut(p+dist*v,v,calcNorm,
352  validNorm,nTmp);
353  dist += disTmp ;
354 
355  if(fPtrSolidB->Inside(p+dist*v) != kOutside)
356  {
357  disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v,calcNorm,
358  validNorm,nTmp);
359  dist += disTmp ;
360  }
361  }
362  while( (fPtrSolidA->Inside(p+dist*v) != kOutside)
363  && (disTmp > 0.5*kCarTolerance) );
364  }
365  else // if( positionB != kOutside )
366  {
367  do // Loop checking, 13.08.2015, G.Cosmo
368  {
369  disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v,calcNorm,
370  validNorm,nTmp);
371  dist += disTmp ;
372 
373  if(fPtrSolidA->Inside(p+dist*v) != kOutside)
374  {
375  disTmp = fPtrSolidA->DistanceToOut(p+dist*v,v,calcNorm,
376  validNorm,nTmp);
377  dist += disTmp ;
378  }
379  }
380  while( (fPtrSolidB->Inside(p+dist*v) != kOutside)
381  && (disTmp > 0.5*kCarTolerance) );
382  }
383  }
384  if( calcNorm )
385  {
386  *validNorm = false ;
387  *n = *nTmp ;
388  }
389  return dist ;
390 }
391 
393 //
394 // Inverted algorithm of DistanceToIn(p)
395 
396 G4double
398 {
399  G4double distout = 0.0;
400  if( Inside(p) == kOutside )
401  {
402 #ifdef G4BOOLDEBUG
403  G4cout << "WARNING - Invalid call in "
404  << "G4UnionSolid::DistanceToOut(p)" << G4endl
405  << " Point p is outside !" << G4endl;
406  G4cout << " p = " << p << G4endl;
407  G4cerr << "WARNING - Invalid call in "
408  << "G4UnionSolid::DistanceToOut(p)" << G4endl
409  << " Point p is outside !" << G4endl;
410  G4cerr << " p = " << p << G4endl;
411 #endif
412  }
413  else
414  {
415  EInside positionA = fPtrSolidA->Inside(p) ;
416  EInside positionB = fPtrSolidB->Inside(p) ;
417 
418  // Is this equivalent ??
419  // if( ! ( (positionA == kOutside)) &&
420  // (positionB == kOutside)) )
421  if((positionA == kInside && positionB == kInside ) ||
422  (positionA == kInside && positionB == kSurface ) ||
423  (positionA == kSurface && positionB == kInside ) )
424  {
425  distout= std::max(fPtrSolidA->DistanceToOut(p),
426  fPtrSolidB->DistanceToOut(p) ) ;
427  }
428  else
429  {
430  if(positionA == kOutside)
431  {
432  distout= fPtrSolidB->DistanceToOut(p) ;
433  }
434  else
435  {
436  distout= fPtrSolidA->DistanceToOut(p) ;
437  }
438  }
439  }
440  return distout;
441 }
442 
444 //
445 //
446 
448 {
449  return G4String("G4UnionSolid");
450 }
451 
453 //
454 // Make a clone of the object
455 
457 {
458  return new G4UnionSolid(*this);
459 }
460 
462 //
463 //
464 
465 void
467  const G4int,
468  const G4VPhysicalVolume* )
469 {
470 }
471 
473 //
474 //
475 
476 void
478 {
479  scene.AddSolid (*this);
480 }
481 
483 //
484 //
485 
486 G4Polyhedron*
488 {
489  HepPolyhedronProcessor processor;
490  // Stack components and components of components recursively
491  // See G4BooleanSolid::StackPolyhedron
492  G4Polyhedron* top = StackPolyhedron(processor, this);
493  G4Polyhedron* result = new G4Polyhedron(*top);
494  if (processor.execute(*result)) { return result; }
495  else { return 0; }
496 }
G4VSolid * fPtrSolidB
virtual G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const =0
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
EInside Inside(const G4ThreeVector &p) const
CLHEP::HepRotation G4RotationMatrix
virtual ~G4UnionSolid()
G4double a
Definition: TRTMaterials.hh:39
virtual void AddSolid(const G4Box &)=0
void DescribeYourselfTo(G4VGraphicsScene &scene) const
int G4int
Definition: G4Types.hh:78
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4GLOB_DLL std::ostream G4cout
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
virtual EInside Inside(const G4ThreeVector &p) const =0
G4double GetRadialTolerance() const
bool G4bool
Definition: G4Types.hh:79
G4Polyhedron * CreatePolyhedron() const
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const =0
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const =0
HepGeom::Transform3D G4Transform3D
const G4int n
#define processor
Definition: xmlparse.cc:617
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
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
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4VSolid * fPtrSolidA
#define G4endl
Definition: G4ios.hh:61
G4double kCarTolerance
Definition: G4VSolid.hh:304
G4BooleanSolid & operator=(const G4BooleanSolid &rhs)
double G4double
Definition: G4Types.hh:76
G4VSolid * Clone() const
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
G4GeometryType GetEntityType() const
static const double mm
Definition: G4SIunits.hh:114
G4Polyhedron * StackPolyhedron(HepPolyhedronProcessor &, const G4VSolid *) const
static G4GeometryTolerance * GetInstance()
G4GLOB_DLL std::ostream G4cerr
G4UnionSolid & operator=(const G4UnionSolid &rhs)
G4UnionSolid(const G4String &pName, G4VSolid *pSolidA, G4VSolid *pSolidB)
Definition: G4UnionSolid.cc:58