G4SubtractionSolid.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 // Implementation of methods for the class G4IntersectionSolid
00030 //
00031 // History:
00032 //
00033 // 14.10.98 V.Grichine: implementation of the first version 
00034 // 19.10.98 V.Grichine: new algorithm of DistanceToIn(p,v)
00035 // 02.08.99 V.Grichine: bugs fixed in DistanceToOut(p,v,...)
00036 //                      while -> do-while & surfaceA limitations
00037 // 13.09.00 V.Grichine: bug fixed in SurfaceNormal(p), p can be inside
00038 // 22.07.11 T.Nikitina: add detection of Infinite Loop in DistanceToIn(p,v)
00039 //
00040 // --------------------------------------------------------------------
00041 
00042 #include <sstream>
00043 
00044 #include "G4SubtractionSolid.hh"
00045 
00046 #include "G4SystemOfUnits.hh"
00047 #include "G4VoxelLimits.hh"
00048 #include "G4VPVParameterisation.hh"
00049 #include "G4GeometryTolerance.hh"
00050 
00051 #include "G4VGraphicsScene.hh"
00052 #include "G4Polyhedron.hh"
00053 #include "HepPolyhedronProcessor.h"
00054 #include "G4NURBS.hh"
00055 // #include "G4NURBSbox.hh"
00056 
00058 //
00059 // Transfer all data members to G4BooleanSolid which is responsible
00060 // for them. pName will be in turn sent to G4VSolid
00061 
00062 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
00063                                               G4VSolid* pSolidA ,
00064                                               G4VSolid* pSolidB   )
00065   : G4BooleanSolid(pName,pSolidA,pSolidB)
00066 {
00067 }
00068 
00070 //
00071 // Constructor
00072  
00073 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
00074                                               G4VSolid* pSolidA ,
00075                                               G4VSolid* pSolidB ,
00076                                               G4RotationMatrix* rotMatrix,
00077                                         const G4ThreeVector& transVector )
00078   : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
00079 {
00080 } 
00081 
00083 //
00084 // Constructor
00085 
00086 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
00087                                               G4VSolid* pSolidA ,
00088                                               G4VSolid* pSolidB ,
00089                                         const G4Transform3D& transform )
00090   : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
00091 {
00092 }
00093 
00095 //
00096 // Fake default constructor - sets only member data and allocates memory
00097 //                            for usage restricted to object persistency.
00098 
00099 G4SubtractionSolid::G4SubtractionSolid( __void__& a )
00100   : G4BooleanSolid(a)
00101 {
00102 }
00103 
00105 //
00106 // Destructor
00107 
00108 G4SubtractionSolid::~G4SubtractionSolid()
00109 {
00110 }
00111 
00113 //
00114 // Copy constructor
00115 
00116 G4SubtractionSolid::G4SubtractionSolid(const G4SubtractionSolid& rhs)
00117   : G4BooleanSolid (rhs)
00118 {
00119 }
00120 
00122 //
00123 // Assignment operator
00124 
00125 G4SubtractionSolid&
00126 G4SubtractionSolid::operator = (const G4SubtractionSolid& rhs) 
00127 {
00128   // Check assignment to self
00129   //
00130   if (this == &rhs)  { return *this; }
00131 
00132   // Copy base class data
00133   //
00134   G4BooleanSolid::operator=(rhs);
00135 
00136   return *this;
00137 }  
00138 
00140 //
00141 // CalculateExtent
00142      
00143 G4bool 
00144 G4SubtractionSolid::CalculateExtent( const EAxis pAxis,
00145                                      const G4VoxelLimits& pVoxelLimit,
00146                                      const G4AffineTransform& pTransform,
00147                                            G4double& pMin,
00148                                            G4double& pMax ) const 
00149 {
00150   // Since we cannot be sure how much the second solid subtracts 
00151   // from the first,    we must use the first solid's extent!
00152 
00153   return fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit, 
00154                                       pTransform, pMin, pMax );
00155 }
00156  
00158 //
00159 // Touching ? Empty subtraction ?
00160 
00161 EInside G4SubtractionSolid::Inside( const G4ThreeVector& p ) const
00162 {
00163   EInside positionA = fPtrSolidA->Inside(p);
00164   if (positionA == kOutside) return kOutside;
00165 
00166   EInside positionB = fPtrSolidB->Inside(p);
00167   
00168   if(positionA == kInside && positionB == kOutside)
00169   {
00170     return kInside ;
00171   }
00172   else
00173   {
00174     if(( positionA == kInside && positionB == kSurface) ||
00175        ( positionB == kOutside && positionA == kSurface) ||
00176        ( positionA == kSurface && positionB == kSurface &&
00177          ( fPtrSolidA->SurfaceNormal(p) - 
00178            fPtrSolidB->SurfaceNormal(p) ).mag2() > 
00179             1000*G4GeometryTolerance::GetInstance()->GetRadialTolerance() ) )
00180     {
00181       return kSurface;
00182     }
00183     else
00184     {
00185       return kOutside;
00186     }
00187   }
00188 }
00189 
00191 //
00192 // SurfaceNormal
00193 
00194 G4ThreeVector 
00195 G4SubtractionSolid::SurfaceNormal( const G4ThreeVector& p ) const 
00196 {
00197   G4ThreeVector normal;
00198   if( Inside(p) == kOutside )
00199   {
00200 #ifdef G4BOOLDEBUG
00201     G4cout << "WARNING - Invalid call [1] in "
00202            << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
00203            << "  Point p is inside !" << G4endl;
00204     G4cout << "          p = " << p << G4endl;
00205     G4cerr << "WARNING - Invalid call [1] in "
00206            << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
00207            << "  Point p is inside !" << G4endl;
00208     G4cerr << "          p = " << p << G4endl;
00209 #endif
00210   }
00211   else
00212   { 
00213     if( fPtrSolidA->Inside(p) == kSurface && 
00214         fPtrSolidB->Inside(p) != kInside      ) 
00215     {
00216       normal = fPtrSolidA->SurfaceNormal(p) ;
00217     }
00218     else if( fPtrSolidA->Inside(p) == kInside && 
00219              fPtrSolidB->Inside(p) != kOutside    )
00220     {
00221       normal = -fPtrSolidB->SurfaceNormal(p) ;
00222     }
00223     else 
00224     {
00225       if ( fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToIn(p) )
00226       {
00227         normal = fPtrSolidA->SurfaceNormal(p) ;
00228       }
00229       else
00230       {
00231         normal = -fPtrSolidB->SurfaceNormal(p) ;
00232       }
00233 #ifdef G4BOOLDEBUG
00234       if(Inside(p) == kInside)
00235       {
00236         G4cout << "WARNING - Invalid call [2] in "
00237              << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
00238              << "  Point p is inside !" << G4endl;
00239         G4cout << "          p = " << p << G4endl;
00240         G4cerr << "WARNING - Invalid call [2] in "
00241              << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
00242              << "  Point p is inside !" << G4endl;
00243         G4cerr << "          p = " << p << G4endl;
00244       }
00245 #endif
00246     }
00247   }
00248   return normal;
00249 }
00250 
00252 //
00253 // The same algorithm as in DistanceToIn(p)
00254 
00255 G4double 
00256 G4SubtractionSolid::DistanceToIn(  const G4ThreeVector& p,
00257                                    const G4ThreeVector& v  ) const 
00258 {
00259   G4double dist = 0.0,disTmp = 0.0 ;
00260     
00261 #ifdef G4BOOLDEBUG
00262   if( Inside(p) == kInside )
00263   {
00264     G4cout << "WARNING - Invalid call in "
00265            << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
00266            << "  Point p is inside !" << G4endl;
00267     G4cout << "          p = " << p << G4endl;
00268     G4cout << "          v = " << v << G4endl;
00269     G4cerr << "WARNING - Invalid call in "
00270            << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
00271            << "  Point p is inside !" << G4endl;
00272     G4cerr << "          p = " << p << G4endl;
00273     G4cerr << "          v = " << v << G4endl;
00274   }
00275 #endif
00276 
00277     // if( // ( fPtrSolidA->Inside(p) != kOutside) &&  // case1:p in both A&B 
00278     if ( fPtrSolidB->Inside(p) != kOutside )   // start: out of B
00279     {
00280       dist = fPtrSolidB->DistanceToOut(p,v) ; // ,calcNorm,validNorm,n) ;
00281       
00282       if( fPtrSolidA->Inside(p+dist*v) != kInside )
00283       {
00284         G4int count1=0;
00285         do
00286         {
00287           disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
00288 
00289           if(disTmp == kInfinity)
00290           {  
00291             return kInfinity ;
00292           }
00293           dist += disTmp ;
00294 
00295           if( Inside(p+dist*v) == kOutside )
00296           {
00297             disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ; 
00298             dist += disTmp ;
00299             count1++;
00300             if( count1 > 1000 )  // Infinite loop detected
00301             {
00302               G4String nameB = fPtrSolidB->GetName();
00303               if(fPtrSolidB->GetEntityType()=="G4DisplacedSolid")
00304               {
00305                 nameB = (dynamic_cast<G4DisplacedSolid*>(fPtrSolidB))
00306                         ->GetConstituentMovedSolid()->GetName();
00307               }
00308               std::ostringstream message;
00309               message << "Illegal condition caused by solids: "
00310                       << fPtrSolidA->GetName() << " and " << nameB << G4endl;
00311               message.precision(16);
00312               message << "Looping detected in point " << p+dist*v
00313                       << ", from original point " << p
00314                       << " and direction " << v << G4endl
00315                       << "Computed candidate distance: " << dist << "*mm.";
00316               message.precision(6);
00317               DumpInfo();
00318               G4Exception("G4SubtractionSolid::DistanceToIn(p,v)",
00319                           "GeomSolids1001", JustWarning, message,
00320                           "Returning candidate distance.");
00321               return dist;
00322             }
00323           }    
00324         }
00325         while( Inside(p+dist*v) == kOutside ) ;
00326       }
00327     }
00328     else // p outside A, start in A
00329     {
00330       dist = fPtrSolidA->DistanceToIn(p,v) ;
00331 
00332       if( dist == kInfinity ) // past A, hence past A\B
00333       {
00334         return kInfinity ;
00335       }
00336       else
00337       {
00338         G4int count2=0;
00339         while( Inside(p+dist*v) == kOutside )  // pushing loop
00340         {
00341           disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
00342           dist += disTmp ;
00343 
00344           if( Inside(p+dist*v) == kOutside )
00345           { 
00346             disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
00347 
00348             if(disTmp == kInfinity) // past A, hence past A\B
00349             {  
00350               return kInfinity ;
00351             }                 
00352             dist += disTmp ;
00353             count2++;
00354             if( count2 > 1000 )  // Infinite loop detected
00355             {
00356               G4String nameB = fPtrSolidB->GetName();
00357               if(fPtrSolidB->GetEntityType()=="G4DisplacedSolid")
00358               {
00359                 nameB = (dynamic_cast<G4DisplacedSolid*>(fPtrSolidB))
00360                         ->GetConstituentMovedSolid()->GetName();
00361               }
00362               std::ostringstream message;
00363               message << "Illegal condition caused by solids: "
00364                       << fPtrSolidA->GetName() << " and " << nameB << G4endl;
00365               message.precision(16);
00366               message << "Looping detected in point " << p+dist*v
00367                       << ", from original point " << p
00368                       << " and direction " << v << G4endl
00369                       << "Computed candidate distance: " << dist << "*mm.";
00370               message.precision(6);
00371               DumpInfo();
00372               G4Exception("G4SubtractionSolid::DistanceToIn(p,v)",
00373                           "GeomSolids1001", JustWarning, message);
00374               return dist;         
00375 
00376             }
00377           }
00378         }
00379       }
00380     }
00381   
00382   return dist ;
00383 }
00384 
00386 //
00387 // Approximate nearest distance from the point p to the intersection of
00388 // two solids. It is usually underestimated from the point of view of
00389 // isotropic safety
00390 
00391 G4double 
00392 G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p ) const 
00393 {
00394   G4double dist=0.0;
00395 
00396 #ifdef G4BOOLDEBUG
00397   if( Inside(p) == kInside )
00398   {
00399     G4cout << "WARNING - Invalid call in "
00400            << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
00401            << "  Point p is inside !" << G4endl;
00402     G4cout << "          p = " << p << G4endl;
00403     G4cerr << "WARNING - Invalid call in "
00404            << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
00405            << "  Point p is inside !" << G4endl;
00406     G4cerr << "          p = " << p << G4endl;
00407   }
00408 #endif
00409 
00410   if( ( fPtrSolidA->Inside(p) != kOutside) &&   // case 1
00411       ( fPtrSolidB->Inside(p) != kOutside)    )
00412   {
00413       dist= fPtrSolidB->DistanceToOut(p)  ;
00414   }
00415   else
00416   {
00417       dist= fPtrSolidA->DistanceToIn(p) ; 
00418   }
00419   
00420   return dist; 
00421 }
00422 
00424 //
00425 // The same algorithm as DistanceToOut(p)
00426 
00427 G4double 
00428 G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p,
00429                  const G4ThreeVector& v,
00430                  const G4bool calcNorm,
00431                        G4bool *validNorm,
00432                        G4ThreeVector *n ) const 
00433 {
00434 #ifdef G4BOOLDEBUG
00435     if( Inside(p) == kOutside )
00436     {
00437       G4cout << "Position:"  << G4endl << G4endl;
00438       G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl;
00439       G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl;
00440       G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl;
00441       G4cout << "Direction:" << G4endl << G4endl;
00442       G4cout << "v.x() = "   << v.x() << G4endl;
00443       G4cout << "v.y() = "   << v.y() << G4endl;
00444       G4cout << "v.z() = "   << v.z() << G4endl << G4endl;
00445       G4cout << "WARNING - Invalid call in "
00446              << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
00447              << "  Point p is outside !" << G4endl;
00448       G4cout << "          p = " << p << G4endl;
00449       G4cout << "          v = " << v << G4endl;
00450       G4cerr << "WARNING - Invalid call in "
00451              << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
00452              << "  Point p is outside !" << G4endl;
00453       G4cerr << "          p = " << p << G4endl;
00454       G4cerr << "          v = " << v << G4endl;
00455     }
00456 #endif
00457 
00458     G4double distout;
00459     G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,validNorm,n) ;
00460     G4double distB = fPtrSolidB->DistanceToIn(p,v) ;
00461     if(distB < distA)
00462     {
00463       if(calcNorm)
00464       {
00465         *n = -(fPtrSolidB->SurfaceNormal(p+distB*v)) ;
00466         *validNorm = false ;
00467       }
00468       distout= distB ;
00469     }
00470     else
00471     {
00472       distout= distA ; 
00473     } 
00474     return distout;
00475 }
00476 
00478 //
00479 // Inverted algorithm of DistanceToIn(p)
00480 
00481 G4double 
00482 G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p ) const 
00483 {
00484   G4double dist=0.0;
00485 
00486   if( Inside(p) == kOutside )
00487   { 
00488 #ifdef G4BOOLDEBUG
00489     G4cout << "WARNING - Invalid call in "
00490            << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
00491            << "  Point p is outside" << G4endl;
00492     G4cout << "          p = " << p << G4endl;
00493     G4cerr << "WARNING - Invalid call in "
00494            << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
00495            << "  Point p is outside" << G4endl;
00496     G4cerr << "          p = " << p << G4endl;
00497 #endif
00498   }
00499   else
00500   {
00501      dist= std::min(fPtrSolidA->DistanceToOut(p),
00502                       fPtrSolidB->DistanceToIn(p) ) ; 
00503   }
00504   return dist; 
00505 }
00506 
00508 //
00509 //
00510 
00511 G4GeometryType G4SubtractionSolid::GetEntityType() const 
00512 {
00513   return G4String("G4SubtractionSolid");
00514 }
00515 
00517 //
00518 // Make a clone of the object
00519 
00520 G4VSolid* G4SubtractionSolid::Clone() const
00521 {
00522   return new G4SubtractionSolid(*this);
00523 }
00524 
00526 //
00527 //
00528 
00529 void 
00530 G4SubtractionSolid::ComputeDimensions(       G4VPVParameterisation*,
00531                                        const G4int,
00532                                        const G4VPhysicalVolume* ) 
00533 {
00534 }
00535 
00537 //
00538 //                    
00539 
00540 void 
00541 G4SubtractionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 
00542 {
00543   scene.AddSolid (*this);
00544 }
00545 
00547 //
00548 //
00549 
00550 G4Polyhedron* 
00551 G4SubtractionSolid::CreatePolyhedron () const 
00552 {
00553   HepPolyhedronProcessor processor;
00554   // Stack components and components of components recursively
00555   // See G4BooleanSolid::StackPolyhedron
00556   G4Polyhedron* top = StackPolyhedron(processor, this);
00557   G4Polyhedron* result = new G4Polyhedron(*top);
00558   if (processor.execute(*result)) { return result; }
00559   else { return 0; }
00560 }
00561 
00563 //
00564 //
00565 
00566 G4NURBS*      
00567 G4SubtractionSolid::CreateNURBS () const 
00568 {
00569   // Take into account boolean operation - see CreatePolyhedron.
00570   // return new G4NURBSbox (1.0, 1.0, 1.0);
00571   return 0;
00572 }

Generated on Mon May 27 17:49:55 2013 for Geant4 by  doxygen 1.4.7