Geant4  10.01.p02
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 66356 2012-12-18 09:02:32Z 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  EInside positionB = fPtrSolidB->Inside(p);
177 
178  if( positionB == kInside ||
179  ( positionA == kSurface && positionB == kSurface &&
180  ( fPtrSolidA->SurfaceNormal(p) +
181  fPtrSolidB->SurfaceNormal(p) ).mag2() <
183  {
184  return kInside;
185  }
186  else
187  {
188  if( ( positionB == kSurface ) || ( positionA == kSurface ) )
189  { return kSurface; }
190  else
191  { return kOutside; }
192  }
193 }
194 
196 //
197 //
198 
201 {
203 
204 #ifdef G4BOOLDEBUG
205  if( Inside(p) == kOutside )
206  {
207  G4cout << "WARNING - Invalid call in "
208  << "G4UnionSolid::SurfaceNormal(p)" << G4endl
209  << " Point p is outside !" << G4endl;
210  G4cout << " p = " << p << G4endl;
211  G4cerr << "WARNING - Invalid call in "
212  << "G4UnionSolid::SurfaceNormal(p)" << G4endl
213  << " Point p is outside !" << G4endl;
214  G4cerr << " p = " << p << G4endl;
215  }
216 #endif
217 
218  if(fPtrSolidA->Inside(p) == kSurface && fPtrSolidB->Inside(p) != kInside)
219  {
220  normal= fPtrSolidA->SurfaceNormal(p) ;
221  }
222  else if(fPtrSolidB->Inside(p) == kSurface &&
223  fPtrSolidA->Inside(p) != kInside)
224  {
225  normal= fPtrSolidB->SurfaceNormal(p) ;
226  }
227  else
228  {
229  normal= fPtrSolidA->SurfaceNormal(p) ;
230 #ifdef G4BOOLDEBUG
231  if(Inside(p)==kInside)
232  {
233  G4cout << "WARNING - Invalid call in "
234  << "G4UnionSolid::SurfaceNormal(p)" << G4endl
235  << " Point p is inside !" << G4endl;
236  G4cout << " p = " << p << G4endl;
237  G4cerr << "WARNING - Invalid call in "
238  << "G4UnionSolid::SurfaceNormal(p)" << G4endl
239  << " Point p is inside !" << G4endl;
240  G4cerr << " p = " << p << G4endl;
241  }
242 #endif
243  }
244  return normal;
245 }
246 
248 //
249 // The same algorithm as in DistanceToIn(p)
250 
251 G4double
253  const G4ThreeVector& v ) const
254 {
255 #ifdef G4BOOLDEBUG
256  if( Inside(p) == kInside )
257  {
258  G4cout << "WARNING - Invalid call in "
259  << "G4UnionSolid::DistanceToIn(p,v)" << G4endl
260  << " Point p is inside !" << G4endl;
261  G4cout << " p = " << p << G4endl;
262  G4cout << " v = " << v << G4endl;
263  G4cerr << "WARNING - Invalid call in "
264  << "G4UnionSolid::DistanceToIn(p,v)" << G4endl
265  << " Point p is inside !" << G4endl;
266  G4cerr << " p = " << p << G4endl;
267  G4cerr << " v = " << v << G4endl;
268  }
269 #endif
270 
271  return std::min(fPtrSolidA->DistanceToIn(p,v),
272  fPtrSolidB->DistanceToIn(p,v) ) ;
273 }
274 
276 //
277 // Approximate nearest distance from the point p to the union of
278 // two solids
279 
280 G4double
282 {
283 #ifdef G4BOOLDEBUG
284  if( Inside(p) == kInside )
285  {
286  G4cout << "WARNING - Invalid call in "
287  << "G4UnionSolid::DistanceToIn(p)" << G4endl
288  << " Point p is inside !" << G4endl;
289  G4cout << " p = " << p << G4endl;
290  G4cerr << "WARNING - Invalid call in "
291  << "G4UnionSolid::DistanceToIn(p)" << G4endl
292  << " Point p is inside !" << G4endl;
293  G4cerr << " p = " << p << G4endl;
294  }
295 #endif
296  G4double distA = fPtrSolidA->DistanceToIn(p) ;
297  G4double distB = fPtrSolidB->DistanceToIn(p) ;
298  G4double safety = std::min(distA,distB) ;
299  if(safety < 0.0) safety = 0.0 ;
300  return safety ;
301 }
302 
304 //
305 // The same algorithm as DistanceToOut(p)
306 
307 G4double
309  const G4ThreeVector& v,
310  const G4bool calcNorm,
311  G4bool *validNorm,
312  G4ThreeVector *n ) const
313 {
314  G4double dist = 0.0, disTmp = 0.0 ;
315  G4ThreeVector normTmp;
316  G4ThreeVector* nTmp= &normTmp;
317 
318  if( Inside(p) == kOutside )
319  {
320 #ifdef G4BOOLDEBUG
321  G4cout << "Position:" << G4endl << G4endl;
322  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
323  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
324  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
325  G4cout << "Direction:" << G4endl << G4endl;
326  G4cout << "v.x() = " << v.x() << G4endl;
327  G4cout << "v.y() = " << v.y() << G4endl;
328  G4cout << "v.z() = " << v.z() << G4endl << G4endl;
329  G4cout << "WARNING - Invalid call in "
330  << "G4UnionSolid::DistanceToOut(p,v)" << G4endl
331  << " Point p is outside !" << G4endl;
332  G4cout << " p = " << p << G4endl;
333  G4cout << " v = " << v << G4endl;
334  G4cerr << "WARNING - Invalid call in "
335  << "G4UnionSolid::DistanceToOut(p,v)" << G4endl
336  << " Point p is outside !" << G4endl;
337  G4cerr << " p = " << p << G4endl;
338  G4cerr << " v = " << v << G4endl;
339 #endif
340  }
341  else
342  {
343  EInside positionA = fPtrSolidA->Inside(p) ;
344  // EInside positionB = fPtrSolidB->Inside(p) ;
345 
346  if( positionA != kOutside )
347  {
348  do
349  {
350  disTmp = fPtrSolidA->DistanceToOut(p+dist*v,v,calcNorm,
351  validNorm,nTmp) ;
352  dist += disTmp ;
353 
354  if(fPtrSolidB->Inside(p+dist*v) != kOutside)
355  {
356  disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v,calcNorm,
357  validNorm,nTmp) ;
358  dist += disTmp ;
359  }
360  }
361  // while( Inside(p+dist*v) == kInside ) ;
362  while( fPtrSolidA->Inside(p+dist*v) != kOutside &&
363  disTmp > 0.5*kCarTolerance ) ;
364  }
365  else // if( positionB != kOutside )
366  {
367  do
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( Inside(p+dist*v) == kInside ) ;
381  while( (fPtrSolidB->Inside(p+dist*v) != kOutside)
382  && (disTmp > 0.5*kCarTolerance) ) ;
383  }
384  }
385  if( calcNorm )
386  {
387  *validNorm = false ;
388  *n = *nTmp ;
389  }
390  return dist ;
391 }
392 
394 //
395 // Inverted algorithm of DistanceToIn(p)
396 
397 G4double
399 {
400  G4double distout = 0.0;
401  if( Inside(p) == kOutside )
402  {
403 #ifdef G4BOOLDEBUG
404  G4cout << "WARNING - Invalid call in "
405  << "G4UnionSolid::DistanceToOut(p)" << G4endl
406  << " Point p is outside !" << G4endl;
407  G4cout << " p = " << p << G4endl;
408  G4cerr << "WARNING - Invalid call in "
409  << "G4UnionSolid::DistanceToOut(p)" << G4endl
410  << " Point p is outside !" << G4endl;
411  G4cerr << " p = " << p << G4endl;
412 #endif
413  }
414  else
415  {
416  EInside positionA = fPtrSolidA->Inside(p) ;
417  EInside positionB = fPtrSolidB->Inside(p) ;
418 
419  // Is this equivalent ??
420  // if( ! ( (positionA == kOutside)) &&
421  // (positionB == kOutside)) )
422  if((positionA == kInside && positionB == kInside ) ||
423  (positionA == kInside && positionB == kSurface ) ||
424  (positionA == kSurface && positionB == kInside ) )
425  {
426  distout= std::max(fPtrSolidA->DistanceToOut(p),
427  fPtrSolidB->DistanceToOut(p) ) ;
428  }
429  else
430  {
431  if(positionA == kOutside)
432  {
433  distout= fPtrSolidB->DistanceToOut(p) ;
434  }
435  else
436  {
437  distout= fPtrSolidA->DistanceToOut(p) ;
438  }
439  }
440  }
441  return distout;
442 }
443 
445 //
446 //
447 
449 {
450  return G4String("G4UnionSolid");
451 }
452 
454 //
455 // Make a clone of the object
456 
458 {
459  return new G4UnionSolid(*this);
460 }
461 
463 //
464 //
465 
466 void
468  const G4int,
469  const G4VPhysicalVolume* )
470 {
471 }
472 
474 //
475 //
476 
477 void
479 {
480  scene.AddSolid (*this);
481 }
482 
484 //
485 //
486 
487 G4Polyhedron*
489 {
490  HepPolyhedronProcessor processor;
491  // Stack components and components of components recursively
492  // See G4BooleanSolid::StackPolyhedron
493  G4Polyhedron* top = StackPolyhedron(processor, this);
494  G4Polyhedron* result = new G4Polyhedron(*top);
495  if (processor.execute(*result)) { return result; }
496  else { return 0; }
497 }
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:600
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:305
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:102
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