G4IntersectionSolid.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 // 17.02.05 V.Grichine: bug was fixed in DistanceToIn(p,v) based on algorithm
00034 //                      proposed by Dino Bazzacco <dino.bazzacco@pd.infn.it>
00035 // 29.05.01 V.Grichine: bug was fixed in DistanceToIn(p,v)
00036 // 16.03.01 V.Grichine: modifications in CalculateExtent() and Inside()
00037 // 29.07.99 V.Grichine: modifications in DistanceToIn(p,v)
00038 // 12.09.98 V.Grichine: first implementation
00039 //
00040 // --------------------------------------------------------------------
00041 
00042 
00043 #include <sstream>
00044 
00045 #include "G4IntersectionSolid.hh"
00046 
00047 #include "G4SystemOfUnits.hh"
00048 #include "G4VoxelLimits.hh"
00049 #include "G4VPVParameterisation.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 
00063 G4IntersectionSolid::G4IntersectionSolid( const G4String& pName,
00064                                                 G4VSolid* pSolidA ,
00065                                                 G4VSolid* pSolidB   )
00066   : G4BooleanSolid(pName,pSolidA,pSolidB)
00067 {
00068 } 
00069 
00071 //
00072 
00073 G4IntersectionSolid::G4IntersectionSolid( 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 // 
00085  
00086 G4IntersectionSolid::G4IntersectionSolid( 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 G4IntersectionSolid::G4IntersectionSolid( __void__& a )
00100   : G4BooleanSolid(a)
00101 {
00102 }
00103 
00105 //
00106 //
00107 
00108 G4IntersectionSolid::~G4IntersectionSolid()
00109 {
00110 }
00111 
00113 //
00114 // Copy constructor
00115 
00116 G4IntersectionSolid::G4IntersectionSolid(const G4IntersectionSolid& rhs)
00117   : G4BooleanSolid (rhs)
00118 {
00119 }
00120 
00122 //
00123 // Assignment operator
00124 
00125 G4IntersectionSolid&
00126 G4IntersectionSolid::operator = (const G4IntersectionSolid& 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 //
00142      
00143 G4bool 
00144 G4IntersectionSolid::CalculateExtent(const EAxis pAxis,
00145                                      const G4VoxelLimits& pVoxelLimit,
00146                                      const G4AffineTransform& pTransform,
00147                                            G4double& pMin,
00148                                            G4double& pMax) const 
00149 {
00150   G4bool   retA, retB, out;
00151   G4double minA, minB, maxA, maxB; 
00152 
00153   retA = fPtrSolidA
00154           ->CalculateExtent( pAxis, pVoxelLimit, pTransform, minA, maxA);
00155   retB = fPtrSolidB
00156           ->CalculateExtent( pAxis, pVoxelLimit, pTransform, minB, maxB);
00157 
00158   if( retA && retB )
00159   {
00160     pMin = std::max( minA, minB ); 
00161     pMax = std::min( maxA, maxB );
00162     out  = (pMax > pMin); // true;
00163 #ifdef G4BOOLDEBUG
00164     // G4cout.precision(16);
00165     // G4cout<<"pMin = "<<pMin<<"; pMax = "<<pMax<<G4endl;
00166 #endif
00167   }
00168   else out = false;
00169 
00170   return out; // It exists in this slice only if both exist in it.
00171 }
00172  
00174 //
00175 // Touching ? Empty intersection ?
00176 
00177 EInside G4IntersectionSolid::Inside(const G4ThreeVector& p) const
00178 {
00179   EInside positionA = fPtrSolidA->Inside(p) ;
00180 
00181   if( positionA == kOutside ) return kOutside ;
00182 
00183   EInside positionB = fPtrSolidB->Inside(p) ;
00184   
00185   if(positionA == kInside && positionB == kInside)
00186   {
00187     return kInside ;
00188   }
00189   else
00190   {
00191     if((positionA == kInside && positionB == kSurface) ||
00192        (positionB == kInside && positionA == kSurface) ||
00193        (positionA == kSurface && positionB == kSurface)   )
00194     {
00195       return kSurface ;
00196     }
00197     else
00198     {
00199       return kOutside ;
00200     }
00201   }
00202 }
00203 
00205 //
00206 
00207 G4ThreeVector 
00208 G4IntersectionSolid::SurfaceNormal( const G4ThreeVector& p ) const 
00209 {
00210   G4ThreeVector normal;
00211   EInside insideA, insideB;
00212   
00213   insideA= fPtrSolidA->Inside(p);
00214   insideB= fPtrSolidB->Inside(p);
00215 
00216 #ifdef G4BOOLDEBUG
00217   if( (insideA == kOutside) || (insideB == kOutside) )
00218   {
00219     G4cout << "WARNING - Invalid call in "
00220            << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
00221            << "  Point p is outside !" << G4endl;
00222     G4cout << "          p = " << p << G4endl;
00223     G4cerr << "WARNING - Invalid call in "
00224            << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
00225            << "  Point p is outside !" << G4endl;
00226     G4cerr << "          p = " << p << G4endl;
00227   }
00228 #endif
00229 
00230   // OLD: if(fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToOut(p) ) 
00231 
00232   // On the surface of both is difficult ... treat it like on A now!
00233   //
00234   // if( (insideA == kSurface) && (insideB == kSurface) )
00235   //    normal= fPtrSolidA->SurfaceNormal(p) ;
00236   // else 
00237   if( insideA == kSurface )
00238     {
00239       normal= fPtrSolidA->SurfaceNormal(p) ;
00240     }
00241   else if( insideB == kSurface )
00242     {
00243       normal= fPtrSolidB->SurfaceNormal(p) ;
00244     } 
00245     // We are on neither surface, so we should generate an exception
00246   else
00247     {
00248       if(fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToOut(p) ) 
00249    normal= fPtrSolidA->SurfaceNormal(p) ;   
00250       else
00251    normal= fPtrSolidB->SurfaceNormal(p) ;   
00252 #ifdef G4BOOLDEBUG
00253       G4cout << "WARNING - Invalid call in "
00254              << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
00255              << "  Point p is out of surface !" << G4endl;
00256       G4cout << "          p = " << p << G4endl;
00257       G4cerr << "WARNING - Invalid call in "
00258              << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
00259              << "  Point p is out of surface !" << G4endl;
00260       G4cerr << "          p = " << p << G4endl;
00261 #endif
00262     }
00263 
00264   return normal;
00265 }
00266 
00268 //
00269 // The same algorithm as in DistanceToIn(p)
00270 
00271 G4double 
00272 G4IntersectionSolid::DistanceToIn( const G4ThreeVector& p,
00273                                    const G4ThreeVector& v  ) const 
00274 {
00275   G4double dist = 0.0;
00276   if( Inside(p) == kInside )
00277   {
00278 #ifdef G4BOOLDEBUG
00279     G4cout << "WARNING - Invalid call in "
00280            << "G4IntersectionSolid::DistanceToIn(p,v)" << G4endl
00281            << "  Point p is inside !" << G4endl;
00282     G4cout << "          p = " << p << G4endl;
00283     G4cout << "          v = " << v << G4endl;
00284     G4cerr << "WARNING - Invalid call in "
00285            << "G4IntersectionSolid::DistanceToIn(p,v)" << G4endl
00286            << "  Point p is inside !" << G4endl;
00287     G4cerr << "          p = " << p << G4endl;
00288     G4cerr << "          v = " << v << G4endl;
00289 #endif
00290   }
00291   else // if( Inside(p) == kSurface ) 
00292   {
00293     EInside wA = fPtrSolidA->Inside(p);
00294     EInside wB = fPtrSolidB->Inside(p);
00295 
00296     G4ThreeVector pA = p,  pB = p;
00297     G4double      dA = 0., dA1=0., dA2=0.;
00298     G4double      dB = 0., dB1=0., dB2=0.;
00299     G4bool        doA = true, doB = true;
00300 
00301     while(true) 
00302     {
00303       if(doA) 
00304       {
00305         // find next valid range for A
00306 
00307         dA1 = 0.;
00308 
00309         if( wA != kInside ) 
00310         {
00311           dA1 = fPtrSolidA->DistanceToIn(pA, v);
00312 
00313           if( dA1 == kInfinity )   return kInfinity;
00314         
00315           pA += dA1*v;
00316         }
00317         dA2 = dA1 + fPtrSolidA->DistanceToOut(pA, v);
00318       }
00319       dA1 += dA;
00320       dA2 += dA;
00321 
00322       if(doB) 
00323       {
00324         // find next valid range for B
00325 
00326         dB1 = 0.;
00327         if(wB != kInside) 
00328         {
00329           dB1 = fPtrSolidB->DistanceToIn(pB, v);
00330 
00331           if(dB1 == kInfinity)   return kInfinity;
00332         
00333           pB += dB1*v;
00334         }
00335         dB2 = dB1 + fPtrSolidB->DistanceToOut(pB, v);
00336       }
00337       dB1 += dB;
00338       dB2 += dB;
00339 
00340        // check if they overlap
00341 
00342       if( dA1 < dB1 ) 
00343       {
00344         if( dB1 < dA2 )  return dB1;
00345 
00346         dA   = dA2;
00347         pA   = p + dA*v;  // continue from here
00348         wA   = kSurface;
00349         doA  = true;
00350         doB  = false;
00351       }
00352       else 
00353       {
00354         if( dA1 < dB2 )  return dA1;
00355 
00356         dB   = dB2;
00357         pB   = p + dB*v;  // continue from here
00358         wB   = kSurface;
00359         doB  = true;
00360         doA  = false;
00361       }
00362     }
00363   }
00364   return dist ;  
00365 }
00366 
00368 //
00369 // Approximate nearest distance from the point p to the intersection of
00370 // two solids
00371 
00372 G4double 
00373 G4IntersectionSolid::DistanceToIn( const G4ThreeVector& p) const 
00374 {
00375 #ifdef G4BOOLDEBUG
00376   if( Inside(p) == kInside )
00377   {
00378     G4cout << "WARNING - Invalid call in "
00379            << "G4IntersectionSolid::DistanceToIn(p)" << G4endl
00380            << "  Point p is inside !" << G4endl;
00381     G4cout << "          p = " << p << G4endl;
00382     G4cerr << "WARNING - Invalid call in "
00383            << "G4IntersectionSolid::DistanceToIn(p)" << G4endl
00384            << "  Point p is inside !" << G4endl;
00385     G4cerr << "          p = " << p << G4endl;
00386   }
00387 #endif
00388   EInside sideA = fPtrSolidA->Inside(p) ;
00389   EInside sideB = fPtrSolidB->Inside(p) ;
00390   G4double dist=0.0 ;
00391 
00392   if( sideA != kInside && sideB  != kOutside )
00393   {
00394     dist = fPtrSolidA->DistanceToIn(p) ;
00395   }
00396   else
00397   {
00398     if( sideB != kInside  && sideA != kOutside )
00399     {
00400       dist = fPtrSolidB->DistanceToIn(p) ;
00401     }
00402     else
00403     {
00404       dist =  std::min(fPtrSolidA->DistanceToIn(p),
00405                     fPtrSolidB->DistanceToIn(p) ) ; 
00406     }
00407   }
00408   return dist ;
00409 }
00410 
00412 //
00413 // The same algorithm as DistanceToOut(p)
00414 
00415 G4double 
00416 G4IntersectionSolid::DistanceToOut( const G4ThreeVector& p,
00417                                     const G4ThreeVector& v,
00418                                     const G4bool calcNorm,
00419                                           G4bool *validNorm,
00420                                           G4ThreeVector *n      ) const 
00421 {
00422   G4bool         validNormA, validNormB;
00423   G4ThreeVector  nA, nB;
00424 
00425 #ifdef G4BOOLDEBUG
00426   if( Inside(p) == kOutside )
00427   {
00428     G4cout << "Position:"  << G4endl << G4endl;
00429     G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl;
00430     G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl;
00431     G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl;
00432     G4cout << "Direction:" << G4endl << G4endl;
00433     G4cout << "v.x() = "   << v.x() << G4endl;
00434     G4cout << "v.y() = "   << v.y() << G4endl;
00435     G4cout << "v.z() = "   << v.z() << G4endl << G4endl;
00436     G4cout << "WARNING - Invalid call in "
00437            << "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl
00438            << "  Point p is outside !" << G4endl;
00439     G4cout << "          p = " << p << G4endl;
00440     G4cout << "          v = " << v << G4endl;
00441     G4cerr << "WARNING - Invalid call in "
00442            << "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl
00443            << "  Point p is outside !" << G4endl;
00444     G4cerr << "          p = " << p << G4endl;
00445     G4cerr << "          v = " << v << G4endl;
00446   }
00447 #endif
00448   G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,&validNormA,&nA) ;
00449   G4double distB = fPtrSolidB->DistanceToOut(p,v,calcNorm,&validNormB,&nB) ;
00450 
00451   G4double dist = std::min(distA,distB) ; 
00452 
00453   if( calcNorm )
00454   {
00455     if ( distA < distB )
00456     {
00457        *validNorm = validNormA;
00458        *n =         nA;
00459     }
00460     else
00461     {   
00462        *validNorm = validNormB;
00463        *n =         nB;
00464     }
00465   }
00466 
00467   return dist ; 
00468 }
00469 
00471 //
00472 // Inverted algorithm of DistanceToIn(p)
00473 
00474 G4double 
00475 G4IntersectionSolid::DistanceToOut( const G4ThreeVector& p ) const 
00476 {
00477 #ifdef G4BOOLDEBUG
00478   if( Inside(p) == kOutside )
00479   {
00480     G4cout << "WARNING - Invalid call in "
00481            << "G4IntersectionSolid::DistanceToOut(p)" << G4endl
00482            << "  Point p is outside !" << G4endl;
00483     G4cout << "          p = " << p << G4endl;
00484     G4cerr << "WARNING - Invalid call in "
00485            << "G4IntersectionSolid::DistanceToOut(p)" << G4endl
00486            << "  Point p is outside !" << G4endl;
00487     G4cerr << "          p = " << p << G4endl;
00488   }
00489 #endif
00490 
00491   return std::min(fPtrSolidA->DistanceToOut(p),
00492                   fPtrSolidB->DistanceToOut(p) ) ; 
00493 
00494 }
00495 
00497 //
00498 //
00499 
00500 void 
00501 G4IntersectionSolid::ComputeDimensions( G4VPVParameterisation*,
00502                                   const G4int,
00503                                         const G4VPhysicalVolume* ) 
00504 {
00505 }
00506 
00508 //
00509 //                    
00510 
00511 G4GeometryType G4IntersectionSolid::GetEntityType() const 
00512 {
00513   return G4String("G4IntersectionSolid");
00514 }
00515 
00517 //
00518 // Make a clone of the object
00519 
00520 G4VSolid* G4IntersectionSolid::Clone() const
00521 {
00522   return new G4IntersectionSolid(*this);
00523 }
00524 
00526 //
00527 //                    
00528 
00529 void 
00530 G4IntersectionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 
00531 {
00532   scene.AddSolid (*this);
00533 }
00534 
00536 //
00537 //
00538 
00539 G4Polyhedron* 
00540 G4IntersectionSolid::CreatePolyhedron () const 
00541 {
00542   HepPolyhedronProcessor processor;
00543   // Stack components and components of components recursively
00544   // See G4BooleanSolid::StackPolyhedron
00545   G4Polyhedron* top = StackPolyhedron(processor, this);
00546   G4Polyhedron* result = new G4Polyhedron(*top);
00547   if (processor.execute(*result)) { return result; }
00548   else { return 0; }
00549 }
00550 
00552 //
00553 //
00554 
00555 G4NURBS*      
00556 G4IntersectionSolid::CreateNURBS      () const 
00557 {
00558   // Take into account boolean operation - see CreatePolyhedron.
00559   // return new G4NURBSbox (1.0, 1.0, 1.0);
00560   return 0;
00561 }

Generated on Mon May 27 17:48:38 2013 for Geant4 by  doxygen 1.4.7