Geant4  9.6.p02
 All Classes 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$
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 <sstream>
43 
44 #include "G4SubtractionSolid.hh"
45 
46 #include "G4SystemOfUnits.hh"
47 #include "G4VoxelLimits.hh"
48 #include "G4VPVParameterisation.hh"
49 #include "G4GeometryTolerance.hh"
50 
51 #include "G4VGraphicsScene.hh"
52 #include "G4Polyhedron.hh"
53 #include "HepPolyhedronProcessor.h"
54 #include "G4NURBS.hh"
55 // #include "G4NURBSbox.hh"
56 
58 //
59 // Transfer all data members to G4BooleanSolid which is responsible
60 // for them. pName will be in turn sent to G4VSolid
61 
62 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
63  G4VSolid* pSolidA ,
64  G4VSolid* pSolidB )
65  : G4BooleanSolid(pName,pSolidA,pSolidB)
66 {
67 }
68 
70 //
71 // Constructor
72 
73 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
74  G4VSolid* pSolidA ,
75  G4VSolid* pSolidB ,
76  G4RotationMatrix* rotMatrix,
77  const G4ThreeVector& transVector )
78  : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
79 {
80 }
81 
83 //
84 // Constructor
85 
86 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
87  G4VSolid* pSolidA ,
88  G4VSolid* pSolidB ,
89  const G4Transform3D& transform )
90  : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
91 {
92 }
93 
95 //
96 // Fake default constructor - sets only member data and allocates memory
97 // for usage restricted to object persistency.
98 
99 G4SubtractionSolid::G4SubtractionSolid( __void__& a )
100  : G4BooleanSolid(a)
101 {
102 }
103 
105 //
106 // Destructor
107 
108 G4SubtractionSolid::~G4SubtractionSolid()
109 {
110 }
111 
113 //
114 // Copy constructor
115 
116 G4SubtractionSolid::G4SubtractionSolid(const G4SubtractionSolid& rhs)
117  : G4BooleanSolid (rhs)
118 {
119 }
120 
122 //
123 // Assignment operator
124 
125 G4SubtractionSolid&
126 G4SubtractionSolid::operator = (const G4SubtractionSolid& rhs)
127 {
128  // Check assignment to self
129  //
130  if (this == &rhs) { return *this; }
131 
132  // Copy base class data
133  //
134  G4BooleanSolid::operator=(rhs);
135 
136  return *this;
137 }
138 
140 //
141 // CalculateExtent
142 
143 G4bool
144 G4SubtractionSolid::CalculateExtent( const EAxis pAxis,
145  const G4VoxelLimits& pVoxelLimit,
146  const G4AffineTransform& pTransform,
147  G4double& pMin,
148  G4double& pMax ) const
149 {
150  // Since we cannot be sure how much the second solid subtracts
151  // from the first, we must use the first solid's extent!
152 
153  return fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
154  pTransform, pMin, pMax );
155 }
156 
158 //
159 // Touching ? Empty subtraction ?
160 
161 EInside G4SubtractionSolid::Inside( const G4ThreeVector& p ) const
162 {
163  EInside positionA = fPtrSolidA->Inside(p);
164  if (positionA == kOutside) return kOutside;
165 
166  EInside positionB = fPtrSolidB->Inside(p);
167 
168  if(positionA == kInside && positionB == kOutside)
169  {
170  return kInside ;
171  }
172  else
173  {
174  if(( positionA == kInside && positionB == kSurface) ||
175  ( positionB == kOutside && positionA == kSurface) ||
176  ( positionA == kSurface && positionB == kSurface &&
177  ( fPtrSolidA->SurfaceNormal(p) -
178  fPtrSolidB->SurfaceNormal(p) ).mag2() >
179  1000*G4GeometryTolerance::GetInstance()->GetRadialTolerance() ) )
180  {
181  return kSurface;
182  }
183  else
184  {
185  return kOutside;
186  }
187  }
188 }
189 
191 //
192 // SurfaceNormal
193 
194 G4ThreeVector
195 G4SubtractionSolid::SurfaceNormal( const G4ThreeVector& p ) const
196 {
197  G4ThreeVector normal;
198  if( Inside(p) == kOutside )
199  {
200 #ifdef G4BOOLDEBUG
201  G4cout << "WARNING - Invalid call [1] in "
202  << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
203  << " Point p is inside !" << G4endl;
204  G4cout << " p = " << p << G4endl;
205  G4cerr << "WARNING - Invalid call [1] in "
206  << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
207  << " Point p is inside !" << G4endl;
208  G4cerr << " p = " << p << G4endl;
209 #endif
210  }
211  else
212  {
213  if( fPtrSolidA->Inside(p) == kSurface &&
214  fPtrSolidB->Inside(p) != kInside )
215  {
216  normal = fPtrSolidA->SurfaceNormal(p) ;
217  }
218  else if( fPtrSolidA->Inside(p) == kInside &&
219  fPtrSolidB->Inside(p) != kOutside )
220  {
221  normal = -fPtrSolidB->SurfaceNormal(p) ;
222  }
223  else
224  {
225  if ( fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToIn(p) )
226  {
227  normal = fPtrSolidA->SurfaceNormal(p) ;
228  }
229  else
230  {
231  normal = -fPtrSolidB->SurfaceNormal(p) ;
232  }
233 #ifdef G4BOOLDEBUG
234  if(Inside(p) == kInside)
235  {
236  G4cout << "WARNING - Invalid call [2] in "
237  << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
238  << " Point p is inside !" << G4endl;
239  G4cout << " p = " << p << G4endl;
240  G4cerr << "WARNING - Invalid call [2] in "
241  << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
242  << " Point p is inside !" << G4endl;
243  G4cerr << " p = " << p << G4endl;
244  }
245 #endif
246  }
247  }
248  return normal;
249 }
250 
252 //
253 // The same algorithm as in DistanceToIn(p)
254 
255 G4double
256 G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p,
257  const G4ThreeVector& v ) const
258 {
259  G4double dist = 0.0,disTmp = 0.0 ;
260 
261 #ifdef G4BOOLDEBUG
262  if( Inside(p) == kInside )
263  {
264  G4cout << "WARNING - Invalid call in "
265  << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
266  << " Point p is inside !" << G4endl;
267  G4cout << " p = " << p << G4endl;
268  G4cout << " v = " << v << G4endl;
269  G4cerr << "WARNING - Invalid call in "
270  << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
271  << " Point p is inside !" << G4endl;
272  G4cerr << " p = " << p << G4endl;
273  G4cerr << " v = " << v << G4endl;
274  }
275 #endif
276 
277  // if( // ( fPtrSolidA->Inside(p) != kOutside) && // case1:p in both A&B
278  if ( fPtrSolidB->Inside(p) != kOutside ) // start: out of B
279  {
280  dist = fPtrSolidB->DistanceToOut(p,v) ; // ,calcNorm,validNorm,n) ;
281 
282  if( fPtrSolidA->Inside(p+dist*v) != kInside )
283  {
284  G4int count1=0;
285  do
286  {
287  disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
288 
289  if(disTmp == kInfinity)
290  {
291  return kInfinity ;
292  }
293  dist += disTmp ;
294 
295  if( Inside(p+dist*v) == kOutside )
296  {
297  disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
298  dist += disTmp ;
299  count1++;
300  if( count1 > 1000 ) // Infinite loop detected
301  {
302  G4String nameB = fPtrSolidB->GetName();
303  if(fPtrSolidB->GetEntityType()=="G4DisplacedSolid")
304  {
305  nameB = (dynamic_cast<G4DisplacedSolid*>(fPtrSolidB))
306  ->GetConstituentMovedSolid()->GetName();
307  }
308  std::ostringstream message;
309  message << "Illegal condition caused by solids: "
310  << fPtrSolidA->GetName() << " and " << nameB << G4endl;
311  message.precision(16);
312  message << "Looping detected in point " << p+dist*v
313  << ", from original point " << p
314  << " and direction " << v << G4endl
315  << "Computed candidate distance: " << dist << "*mm.";
316  message.precision(6);
317  DumpInfo();
318  G4Exception("G4SubtractionSolid::DistanceToIn(p,v)",
319  "GeomSolids1001", JustWarning, message,
320  "Returning candidate distance.");
321  return dist;
322  }
323  }
324  }
325  while( Inside(p+dist*v) == kOutside ) ;
326  }
327  }
328  else // p outside A, start in A
329  {
330  dist = fPtrSolidA->DistanceToIn(p,v) ;
331 
332  if( dist == kInfinity ) // past A, hence past A\B
333  {
334  return kInfinity ;
335  }
336  else
337  {
338  G4int count2=0;
339  while( Inside(p+dist*v) == kOutside ) // pushing loop
340  {
341  disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
342  dist += disTmp ;
343 
344  if( Inside(p+dist*v) == kOutside )
345  {
346  disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
347 
348  if(disTmp == kInfinity) // past A, hence past A\B
349  {
350  return kInfinity ;
351  }
352  dist += disTmp ;
353  count2++;
354  if( count2 > 1000 ) // Infinite loop detected
355  {
356  G4String nameB = fPtrSolidB->GetName();
357  if(fPtrSolidB->GetEntityType()=="G4DisplacedSolid")
358  {
359  nameB = (dynamic_cast<G4DisplacedSolid*>(fPtrSolidB))
360  ->GetConstituentMovedSolid()->GetName();
361  }
362  std::ostringstream message;
363  message << "Illegal condition caused by solids: "
364  << fPtrSolidA->GetName() << " and " << nameB << G4endl;
365  message.precision(16);
366  message << "Looping detected in point " << p+dist*v
367  << ", from original point " << p
368  << " and direction " << v << G4endl
369  << "Computed candidate distance: " << dist << "*mm.";
370  message.precision(6);
371  DumpInfo();
372  G4Exception("G4SubtractionSolid::DistanceToIn(p,v)",
373  "GeomSolids1001", JustWarning, message);
374  return dist;
375 
376  }
377  }
378  }
379  }
380  }
381 
382  return dist ;
383 }
384 
386 //
387 // Approximate nearest distance from the point p to the intersection of
388 // two solids. It is usually underestimated from the point of view of
389 // isotropic safety
390 
391 G4double
392 G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p ) const
393 {
394  G4double dist=0.0;
395 
396 #ifdef G4BOOLDEBUG
397  if( Inside(p) == kInside )
398  {
399  G4cout << "WARNING - Invalid call in "
400  << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
401  << " Point p is inside !" << G4endl;
402  G4cout << " p = " << p << G4endl;
403  G4cerr << "WARNING - Invalid call in "
404  << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
405  << " Point p is inside !" << G4endl;
406  G4cerr << " p = " << p << G4endl;
407  }
408 #endif
409 
410  if( ( fPtrSolidA->Inside(p) != kOutside) && // case 1
411  ( fPtrSolidB->Inside(p) != kOutside) )
412  {
413  dist= fPtrSolidB->DistanceToOut(p) ;
414  }
415  else
416  {
417  dist= fPtrSolidA->DistanceToIn(p) ;
418  }
419 
420  return dist;
421 }
422 
424 //
425 // The same algorithm as DistanceToOut(p)
426 
427 G4double
428 G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p,
429  const G4ThreeVector& v,
430  const G4bool calcNorm,
431  G4bool *validNorm,
432  G4ThreeVector *n ) const
433 {
434 #ifdef G4BOOLDEBUG
435  if( Inside(p) == kOutside )
436  {
437  G4cout << "Position:" << G4endl << G4endl;
438  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
439  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
440  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
441  G4cout << "Direction:" << G4endl << G4endl;
442  G4cout << "v.x() = " << v.x() << G4endl;
443  G4cout << "v.y() = " << v.y() << G4endl;
444  G4cout << "v.z() = " << v.z() << G4endl << G4endl;
445  G4cout << "WARNING - Invalid call in "
446  << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
447  << " Point p is outside !" << G4endl;
448  G4cout << " p = " << p << G4endl;
449  G4cout << " v = " << v << G4endl;
450  G4cerr << "WARNING - Invalid call in "
451  << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
452  << " Point p is outside !" << G4endl;
453  G4cerr << " p = " << p << G4endl;
454  G4cerr << " v = " << v << G4endl;
455  }
456 #endif
457 
458  G4double distout;
459  G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,validNorm,n) ;
460  G4double distB = fPtrSolidB->DistanceToIn(p,v) ;
461  if(distB < distA)
462  {
463  if(calcNorm)
464  {
465  *n = -(fPtrSolidB->SurfaceNormal(p+distB*v)) ;
466  *validNorm = false ;
467  }
468  distout= distB ;
469  }
470  else
471  {
472  distout= distA ;
473  }
474  return distout;
475 }
476 
478 //
479 // Inverted algorithm of DistanceToIn(p)
480 
481 G4double
482 G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p ) const
483 {
484  G4double dist=0.0;
485 
486  if( Inside(p) == kOutside )
487  {
488 #ifdef G4BOOLDEBUG
489  G4cout << "WARNING - Invalid call in "
490  << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
491  << " Point p is outside" << G4endl;
492  G4cout << " p = " << p << G4endl;
493  G4cerr << "WARNING - Invalid call in "
494  << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
495  << " Point p is outside" << G4endl;
496  G4cerr << " p = " << p << G4endl;
497 #endif
498  }
499  else
500  {
501  dist= std::min(fPtrSolidA->DistanceToOut(p),
502  fPtrSolidB->DistanceToIn(p) ) ;
503  }
504  return dist;
505 }
506 
508 //
509 //
510 
511 G4GeometryType G4SubtractionSolid::GetEntityType() const
512 {
513  return G4String("G4SubtractionSolid");
514 }
515 
517 //
518 // Make a clone of the object
519 
520 G4VSolid* G4SubtractionSolid::Clone() const
521 {
522  return new G4SubtractionSolid(*this);
523 }
524 
526 //
527 //
528 
529 void
530 G4SubtractionSolid::ComputeDimensions( G4VPVParameterisation*,
531  const G4int,
532  const G4VPhysicalVolume* )
533 {
534 }
535 
537 //
538 //
539 
540 void
541 G4SubtractionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
542 {
543  scene.AddSolid (*this);
544 }
545 
547 //
548 //
549 
550 G4Polyhedron*
551 G4SubtractionSolid::CreatePolyhedron () const
552 {
553  HepPolyhedronProcessor processor;
554  // Stack components and components of components recursively
555  // See G4BooleanSolid::StackPolyhedron
556  G4Polyhedron* top = StackPolyhedron(processor, this);
557  G4Polyhedron* result = new G4Polyhedron(*top);
558  if (processor.execute(*result)) { return result; }
559  else { return 0; }
560 }
561 
563 //
564 //
565 
566 G4NURBS*
567 G4SubtractionSolid::CreateNURBS () const
568 {
569  // Take into account boolean operation - see CreatePolyhedron.
570  // return new G4NURBSbox (1.0, 1.0, 1.0);
571  return 0;
572 }
Generated on Sat May 25 2013 14:33:12 for Geant4 by   doxygen 1.8.4