Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4SubtractionSolid.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: G4SubtractionSolid.cc 66356 2012-12-18 09:02:32Z gcosmo $
28 //
29 // Implementation of methods for the class G4IntersectionSolid
30 //
31 // History:
32 //
33 // 14.10.98 V.Grichine: implementation of the first version
34 // 19.10.98 V.Grichine: new algorithm of DistanceToIn(p,v)
35 // 02.08.99 V.Grichine: bugs fixed in DistanceToOut(p,v,...)
36 // while -> do-while & surfaceA limitations
37 // 13.09.00 V.Grichine: bug fixed in SurfaceNormal(p), p can be inside
38 // 22.07.11 T.Nikitina: add detection of Infinite Loop in DistanceToIn(p,v)
39 //
40 // --------------------------------------------------------------------
41 
42 #include "G4SubtractionSolid.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 
53 #include <sstream>
54 
55 ///////////////////////////////////////////////////////////////////
56 //
57 // Transfer all data members to G4BooleanSolid which is responsible
58 // for them. pName will be in turn sent to G4VSolid
59 
61  G4VSolid* pSolidA ,
62  G4VSolid* pSolidB )
63  : G4BooleanSolid(pName,pSolidA,pSolidB)
64 {
65 }
66 
67 ///////////////////////////////////////////////////////////////
68 //
69 // Constructor
70 
72  G4VSolid* pSolidA ,
73  G4VSolid* pSolidB ,
74  G4RotationMatrix* rotMatrix,
75  const G4ThreeVector& transVector )
76  : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
77 {
78 }
79 
80 ///////////////////////////////////////////////////////////////
81 //
82 // Constructor
83 
85  G4VSolid* pSolidA ,
86  G4VSolid* pSolidB ,
87  const G4Transform3D& transform )
88  : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
89 {
90 }
91 
92 //////////////////////////////////////////////////////////////////
93 //
94 // Fake default constructor - sets only member data and allocates memory
95 // for usage restricted to object persistency.
96 
98  : G4BooleanSolid(a)
99 {
100 }
101 
102 ///////////////////////////////////////////////////////////////
103 //
104 // Destructor
105 
107 {
108 }
109 
110 ///////////////////////////////////////////////////////////////
111 //
112 // Copy constructor
113 
115  : G4BooleanSolid (rhs)
116 {
117 }
118 
119 ///////////////////////////////////////////////////////////////
120 //
121 // Assignment operator
122 
125 {
126  // Check assignment to self
127  //
128  if (this == &rhs) { return *this; }
129 
130  // Copy base class data
131  //
133 
134  return *this;
135 }
136 
137 ///////////////////////////////////////////////////////////////
138 //
139 // CalculateExtent
140 
141 G4bool
143  const G4VoxelLimits& pVoxelLimit,
144  const G4AffineTransform& pTransform,
145  G4double& pMin,
146  G4double& pMax ) const
147 {
148  // Since we cannot be sure how much the second solid subtracts
149  // from the first, we must use the first solid's extent!
150 
151  return fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
152  pTransform, pMin, pMax );
153 }
154 
155 /////////////////////////////////////////////////////
156 //
157 // Touching ? Empty subtraction ?
158 
160 {
161  EInside positionA = fPtrSolidA->Inside(p);
162  if (positionA == kOutside) return kOutside;
163 
164  EInside positionB = fPtrSolidB->Inside(p);
165 
166  if(positionA == kInside && positionB == kOutside)
167  {
168  return kInside ;
169  }
170  else
171  {
172  if(( positionA == kInside && positionB == kSurface) ||
173  ( positionB == kOutside && positionA == kSurface) ||
174  ( positionA == kSurface && positionB == kSurface &&
175  ( fPtrSolidA->SurfaceNormal(p) -
176  fPtrSolidB->SurfaceNormal(p) ).mag2() >
178  {
179  return kSurface;
180  }
181  else
182  {
183  return kOutside;
184  }
185  }
186 }
187 
188 //////////////////////////////////////////////////////////////
189 //
190 // SurfaceNormal
191 
194 {
195  G4ThreeVector normal;
196  EInside insideThis= Inside(p);
197  if( insideThis == kOutside )
198  {
199 #ifdef G4BOOLDEBUG
200  G4cout << "WARNING - Invalid call [1] in "
201  << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
202  << " Point p is outside !" << G4endl;
203  G4cout << " p = " << p << G4endl;
204  G4cerr << "WARNING - Invalid call [1] in "
205  << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
206  << " Point p is outside !" << G4endl;
207  G4cerr << " p = " << p << G4endl;
208 #endif
209  }
210  else
211  {
212  EInside InsideA = fPtrSolidA->Inside(p);
213  EInside InsideB = fPtrSolidB->Inside(p);
214 
215  if( InsideA == kSurface &&
216  InsideB != kInside )
217  {
218  normal = fPtrSolidA->SurfaceNormal(p) ;
219  }
220  else if( InsideA == kInside &&
221  InsideB != kOutside )
222  {
223  normal = -fPtrSolidB->SurfaceNormal(p) ;
224  }
225  else
226  {
228  {
229  normal = fPtrSolidA->SurfaceNormal(p) ;
230  }
231  else
232  {
233  normal = -fPtrSolidB->SurfaceNormal(p) ;
234  }
235 #ifdef G4BOOLDEBUG
236  if(insideThis == kInside)
237  {
238  G4cout << "WARNING - Invalid call [2] in "
239  << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
240  << " Point p is inside !" << G4endl;
241  G4cout << " p = " << p << G4endl;
242  G4cerr << "WARNING - Invalid call [2] in "
243  << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
244  << " Point p is inside !" << G4endl;
245  G4cerr << " p = " << p << G4endl;
246  }
247 #endif
248  }
249  }
250  return normal;
251 }
252 
253 /////////////////////////////////////////////////////////////
254 //
255 // The same algorithm as in DistanceToIn(p)
256 
257 G4double
259  const G4ThreeVector& v ) const
260 {
261  G4double dist = 0.0,disTmp = 0.0 ;
262 
263 #ifdef G4BOOLDEBUG
264  if( Inside(p) == kInside )
265  {
266  G4cout << "WARNING - Invalid call in "
267  << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
268  << " Point p is inside !" << G4endl;
269  G4cout << " p = " << p << G4endl;
270  G4cout << " v = " << v << G4endl;
271  G4cerr << "WARNING - Invalid call in "
272  << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
273  << " Point p is inside !" << G4endl;
274  G4cerr << " p = " << p << G4endl;
275  G4cerr << " v = " << v << G4endl;
276  }
277 #endif
278 
279  // if( // ( fPtrSolidA->Inside(p) != kOutside) && // case1:p in both A&B
280  if ( fPtrSolidB->Inside(p) != kOutside ) // start: out of B
281  {
282  dist = fPtrSolidB->DistanceToOut(p,v) ; // ,calcNorm,validNorm,n) ;
283 
284  if( fPtrSolidA->Inside(p+dist*v) != kInside )
285  {
286  G4int count1=0;
287  do
288  {
289  disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
290 
291  if(disTmp == kInfinity)
292  {
293  return kInfinity ;
294  }
295  dist += disTmp ;
296 
297  if( Inside(p+dist*v) == kOutside )
298  {
299  disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
300  dist += disTmp ;
301  count1++;
302  if( count1 > 1000 ) // Infinite loop detected
303  {
304  G4String nameB = fPtrSolidB->GetName();
305  if(fPtrSolidB->GetEntityType()=="G4DisplacedSolid")
306  {
307  nameB = (dynamic_cast<G4DisplacedSolid*>(fPtrSolidB))
308  ->GetConstituentMovedSolid()->GetName();
309  }
310  std::ostringstream message;
311  message << "Illegal condition caused by solids: "
312  << fPtrSolidA->GetName() << " and " << nameB << G4endl;
313  message.precision(16);
314  message << "Looping detected in point " << p+dist*v
315  << ", from original point " << p
316  << " and direction " << v << G4endl
317  << "Computed candidate distance: " << dist << "*mm. ";
318  message.precision(6);
319  DumpInfo();
320  G4Exception("G4SubtractionSolid::DistanceToIn(p,v)",
321  "GeomSolids1001", JustWarning, message,
322  "Returning candidate distance.");
323  return dist;
324  }
325  }
326  }
327  while( Inside(p+dist*v) == kOutside ) ;
328  }
329  }
330  else // p outside A, start in A
331  {
332  dist = fPtrSolidA->DistanceToIn(p,v) ;
333 
334  if( dist == kInfinity ) // past A, hence past A\B
335  {
336  return kInfinity ;
337  }
338  else
339  {
340  G4int count2=0;
341  while( Inside(p+dist*v) == kOutside ) // pushing loop
342  {
343  disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
344  dist += disTmp ;
345 
346  if( Inside(p+dist*v) == kOutside )
347  {
348  disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
349 
350  if(disTmp == kInfinity) // past A, hence past A\B
351  {
352  return kInfinity ;
353  }
354  dist += disTmp ;
355  count2++;
356  if( count2 > 1000 ) // Infinite loop detected
357  {
358  G4String nameB = fPtrSolidB->GetName();
359  if(fPtrSolidB->GetEntityType()=="G4DisplacedSolid")
360  {
361  nameB = (dynamic_cast<G4DisplacedSolid*>(fPtrSolidB))
362  ->GetConstituentMovedSolid()->GetName();
363  }
364  std::ostringstream message;
365  message << "Illegal condition caused by solids: "
366  << fPtrSolidA->GetName() << " and " << nameB << G4endl;
367  message.precision(16);
368  message << "Looping detected in point " << p+dist*v
369  << ", from original point " << p
370  << " and direction " << v << G4endl
371  << "Computed candidate distance: " << dist << "*mm. ";
372  message.precision(6);
373  DumpInfo();
374  G4Exception("G4SubtractionSolid::DistanceToIn(p,v)",
375  "GeomSolids1001", JustWarning, message,
376  "Returning candidate distance.");
377  return dist;
378  }
379  }
380  }
381  }
382  }
383 
384  return dist ;
385 }
386 
387 ////////////////////////////////////////////////////////
388 //
389 // Approximate nearest distance from the point p to the intersection of
390 // two solids. It is usually underestimated from the point of view of
391 // isotropic safety
392 
393 G4double
395 {
396  G4double dist=0.0;
397 
398 #ifdef G4BOOLDEBUG
399  if( Inside(p) == kInside )
400  {
401  G4cout << "WARNING - Invalid call in "
402  << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
403  << " Point p is inside !" << G4endl;
404  G4cout << " p = " << p << G4endl;
405  G4cerr << "WARNING - Invalid call in "
406  << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
407  << " Point p is inside !" << G4endl;
408  G4cerr << " p = " << p << G4endl;
409  }
410 #endif
411 
412  if( ( fPtrSolidA->Inside(p) != kOutside) && // case 1
413  ( fPtrSolidB->Inside(p) != kOutside) )
414  {
415  dist= fPtrSolidB->DistanceToOut(p) ;
416  }
417  else
418  {
419  dist= fPtrSolidA->DistanceToIn(p) ;
420  }
421 
422  return dist;
423 }
424 
425 //////////////////////////////////////////////////////////
426 //
427 // The same algorithm as DistanceToOut(p)
428 
429 G4double
431  const G4ThreeVector& v,
432  const G4bool calcNorm,
433  G4bool *validNorm,
434  G4ThreeVector *n ) const
435 {
436 #ifdef G4BOOLDEBUG
437  if( Inside(p) == kOutside )
438  {
439  G4cout << "Position:" << G4endl << G4endl;
440  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
441  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
442  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
443  G4cout << "Direction:" << G4endl << G4endl;
444  G4cout << "v.x() = " << v.x() << G4endl;
445  G4cout << "v.y() = " << v.y() << G4endl;
446  G4cout << "v.z() = " << v.z() << G4endl << G4endl;
447  G4cout << "WARNING - Invalid call in "
448  << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
449  << " Point p is outside !" << G4endl;
450  G4cout << " p = " << p << G4endl;
451  G4cout << " v = " << v << G4endl;
452  G4cerr << "WARNING - Invalid call in "
453  << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
454  << " Point p is outside !" << G4endl;
455  G4cerr << " p = " << p << G4endl;
456  G4cerr << " v = " << v << G4endl;
457  }
458 #endif
459 
460  G4double distout;
461  G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,validNorm,n) ;
462  G4double distB = fPtrSolidB->DistanceToIn(p,v) ;
463  if(distB < distA)
464  {
465  if(calcNorm)
466  {
467  *n = -(fPtrSolidB->SurfaceNormal(p+distB*v)) ;
468  *validNorm = false ;
469  }
470  distout= distB ;
471  }
472  else
473  {
474  distout= distA ;
475  }
476  return distout;
477 }
478 
479 //////////////////////////////////////////////////////////////
480 //
481 // Inverted algorithm of DistanceToIn(p)
482 
483 G4double
485 {
486  G4double dist=0.0;
487 
488  if( Inside(p) == kOutside )
489  {
490 #ifdef G4BOOLDEBUG
491  G4cout << "WARNING - Invalid call in "
492  << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
493  << " Point p is outside" << G4endl;
494  G4cout << " p = " << p << G4endl;
495  G4cerr << "WARNING - Invalid call in "
496  << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
497  << " Point p is outside" << G4endl;
498  G4cerr << " p = " << p << G4endl;
499 #endif
500  }
501  else
502  {
504  fPtrSolidB->DistanceToIn(p) ) ;
505  }
506  return dist;
507 }
508 
509 //////////////////////////////////////////////////////////////
510 //
511 //
512 
514 {
515  return G4String("G4SubtractionSolid");
516 }
517 
518 //////////////////////////////////////////////////////////////////////////
519 //
520 // Make a clone of the object
521 
523 {
524  return new G4SubtractionSolid(*this);
525 }
526 
527 //////////////////////////////////////////////////////////////
528 //
529 //
530 
531 void
533  const G4int,
534  const G4VPhysicalVolume* )
535 {
536 }
537 
538 /////////////////////////////////////////////////
539 //
540 //
541 
542 void
544 {
545  scene.AddSolid (*this);
546 }
547 
548 ////////////////////////////////////////////////////
549 //
550 //
551 
552 G4Polyhedron*
554 {
556  // Stack components and components of components recursively
557  // See G4BooleanSolid::StackPolyhedron
558  G4Polyhedron* top = StackPolyhedron(processor, this);
559  G4Polyhedron* result = new G4Polyhedron(*top);
560  if (processor.execute(*result)) { return result; }
561  else { return 0; }
562 }
G4String GetName() const
G4VSolid * fPtrSolidB
G4VSolid * Clone() const
virtual G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const =0
double x() const
const char * p
Definition: xmltok.h:285
virtual G4GeometryType GetEntityType() const =0
virtual void AddSolid(const G4Box &)=0
int G4int
Definition: G4Types.hh:78
double z() const
G4Polyhedron * CreatePolyhedron() const
void DumpInfo() const
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
G4GLOB_DLL std::ostream G4cout
virtual EInside Inside(const G4ThreeVector &p) const =0
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
G4double GetRadialTolerance() const
bool G4bool
Definition: G4Types.hh:79
G4SubtractionSolid & operator=(const G4SubtractionSolid &rhs)
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const =0
EInside Inside(const G4ThreeVector &p) const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const =0
G4SubtractionSolid(const G4String &pName, G4VSolid *pSolidA, G4VSolid *pSolidB)
const G4int n
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void DescribeYourselfTo(G4VGraphicsScene &scene) const
#define processor
Definition: xmlparse.cc:600
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
EInside
Definition: geomdefs.hh:58
EAxis
Definition: geomdefs.hh:54
double y() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4VSolid * fPtrSolidA
#define G4endl
Definition: G4ios.hh:61
G4GeometryType GetEntityType() const
G4BooleanSolid & operator=(const G4BooleanSolid &rhs)
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 =0
G4Polyhedron * StackPolyhedron(HepPolyhedronProcessor &, const G4VSolid *) const
bool execute(HepPolyhedron &)
static G4GeometryTolerance * GetInstance()
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4GLOB_DLL std::ostream G4cerr