00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "G4SubtractionSolid.hh"
00043
00044 #include "G4VoxelLimits.hh"
00045 #include "G4VPVParameterisation.hh"
00046 #include "G4GeometryTolerance.hh"
00047
00048 #include "G4VGraphicsScene.hh"
00049 #include "G4Polyhedron.hh"
00050 #include "G4NURBS.hh"
00051
00052
00053 #include <sstream>
00054
00056
00057
00058
00059
00060 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
00061 G4VSolid* pSolidA ,
00062 G4VSolid* pSolidB )
00063 : G4BooleanSolid(pName,pSolidA,pSolidB)
00064 {
00065 }
00066
00068
00069
00070
00071 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
00072 G4VSolid* pSolidA ,
00073 G4VSolid* pSolidB ,
00074 G4RotationMatrix* rotMatrix,
00075 const G4ThreeVector& transVector )
00076 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
00077 {
00078 }
00079
00081
00082
00083
00084 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
00085 G4VSolid* pSolidA ,
00086 G4VSolid* pSolidB ,
00087 const G4Transform3D& transform )
00088 : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
00089 {
00090 }
00091
00093
00094
00095
00096
00097 G4SubtractionSolid::G4SubtractionSolid( __void__& a )
00098 : G4BooleanSolid(a)
00099 {
00100 }
00101
00103
00104
00105
00106 G4SubtractionSolid::~G4SubtractionSolid()
00107 {
00108 }
00109
00111
00112
00113
00114 G4bool
00115 G4SubtractionSolid::CalculateExtent( const EAxis pAxis,
00116 const G4VoxelLimits& pVoxelLimit,
00117 const G4AffineTransform& pTransform,
00118 G4double& pMin,
00119 G4double& pMax ) const
00120 {
00121
00122
00123
00124 return fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
00125 pTransform, pMin, pMax );
00126 }
00127
00129
00130
00131
00132 EInside G4SubtractionSolid::Inside( const G4ThreeVector& p ) const
00133 {
00134 EInside positionA = fPtrSolidA->Inside(p);
00135 if (positionA == kOutside) return kOutside;
00136
00137 EInside positionB = fPtrSolidB->Inside(p);
00138
00139 if(positionA == kInside && positionB == kOutside)
00140 {
00141 return kInside ;
00142 }
00143 else
00144 {
00145 if(( positionA == kInside && positionB == kSurface) ||
00146 ( positionB == kOutside && positionA == kSurface) ||
00147 ( positionA == kSurface && positionB == kSurface &&
00148 ( fPtrSolidA->SurfaceNormal(p) -
00149 fPtrSolidB->SurfaceNormal(p) ).mag2() >
00150 1000*G4GeometryTolerance::GetInstance()->GetRadialTolerance() ) )
00151 {
00152 return kSurface;
00153 }
00154 else
00155 {
00156 return kOutside;
00157 }
00158 }
00159 }
00160
00162
00163
00164
00165 G4ThreeVector
00166 G4SubtractionSolid::SurfaceNormal( const G4ThreeVector& p ) const
00167 {
00168 G4ThreeVector normal;
00169 if( Inside(p) == kOutside )
00170 {
00171 #ifdef G4BOOLDEBUG
00172 G4cout << "WARNING - Invalid call [1] in "
00173 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
00174 << " Point p is inside !" << G4endl;
00175 G4cout << " p = " << p << G4endl;
00176 G4cerr << "WARNING - Invalid call [1] in "
00177 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
00178 << " Point p is inside !" << G4endl;
00179 G4cerr << " p = " << p << G4endl;
00180 #endif
00181 }
00182 else
00183 {
00184 if( fPtrSolidA->Inside(p) == kSurface &&
00185 fPtrSolidB->Inside(p) != kInside )
00186 {
00187 normal = fPtrSolidA->SurfaceNormal(p) ;
00188 }
00189 else if( fPtrSolidA->Inside(p) == kInside &&
00190 fPtrSolidB->Inside(p) != kOutside )
00191 {
00192 normal = -fPtrSolidB->SurfaceNormal(p) ;
00193 }
00194 else
00195 {
00196 if ( fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToIn(p) )
00197 {
00198 normal = fPtrSolidA->SurfaceNormal(p) ;
00199 }
00200 else
00201 {
00202 normal = -fPtrSolidB->SurfaceNormal(p) ;
00203 }
00204 #ifdef G4BOOLDEBUG
00205 if(Inside(p) == kInside)
00206 {
00207 G4cout << "WARNING - Invalid call [2] in "
00208 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
00209 << " Point p is inside !" << G4endl;
00210 G4cout << " p = " << p << G4endl;
00211 G4cerr << "WARNING - Invalid call [2] in "
00212 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
00213 << " Point p is inside !" << G4endl;
00214 G4cerr << " p = " << p << G4endl;
00215 }
00216 #endif
00217 }
00218 }
00219 return normal;
00220 }
00221
00223
00224
00225
00226 G4double
00227 G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p,
00228 const G4ThreeVector& v ) const
00229 {
00230 G4double dist = 0.0,disTmp = 0.0 ;
00231
00232 #ifdef G4BOOLDEBUG
00233 if( Inside(p) == kInside )
00234 {
00235 G4cout << "WARNING - Invalid call in "
00236 << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
00237 << " Point p is inside !" << G4endl;
00238 G4cout << " p = " << p << G4endl;
00239 G4cout << " v = " << v << G4endl;
00240 G4cerr << "WARNING - Invalid call in "
00241 << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
00242 << " Point p is inside !" << G4endl;
00243 G4cerr << " p = " << p << G4endl;
00244 G4cerr << " v = " << v << G4endl;
00245 }
00246 #endif
00247
00248
00249 if ( fPtrSolidB->Inside(p) != kOutside )
00250 {
00251 dist = fPtrSolidB->DistanceToOut(p,v) ;
00252
00253 if( fPtrSolidA->Inside(p+dist*v) != kInside )
00254 {
00255 do
00256 {
00257 disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
00258
00259 if(disTmp == kInfinity)
00260 {
00261 return kInfinity ;
00262 }
00263 dist += disTmp ;
00264
00265 if( Inside(p+dist*v) == kOutside )
00266 {
00267 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
00268 dist += disTmp ;
00269 }
00270 }
00271 while( Inside(p+dist*v) == kOutside ) ;
00272 }
00273 }
00274 else
00275 {
00276 dist = fPtrSolidA->DistanceToIn(p,v) ;
00277
00278 if( dist == kInfinity )
00279 {
00280 return kInfinity ;
00281 }
00282 else
00283 {
00284 int loopflag=0;
00285 while( Inside(p+dist*v) == kOutside )
00286 {
00287 loopflag++;
00288 if(loopflag>10000) {
00289 std::cout<<"G4SubtractionSolid::DistanceToIn loop >10000"<<std::endl;
00290 return kInfinity ;
00291 }
00292 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
00293 dist += disTmp ;
00294
00295 if( Inside(p+dist*v) == kOutside )
00296 {
00297 disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
00298
00299 if(disTmp == kInfinity)
00300 {
00301 return kInfinity ;
00302 }
00303 dist += disTmp ;
00304 }
00305 }
00306 }
00307 }
00308
00309 return dist ;
00310 }
00311
00313
00314
00315
00316
00317
00318 G4double
00319 G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p ) const
00320 {
00321 G4double dist=0.0;
00322
00323 #ifdef G4BOOLDEBUG
00324 if( Inside(p) == kInside )
00325 {
00326 G4cout << "WARNING - Invalid call in "
00327 << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
00328 << " Point p is inside !" << G4endl;
00329 G4cout << " p = " << p << G4endl;
00330 G4cerr << "WARNING - Invalid call in "
00331 << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
00332 << " Point p is inside !" << G4endl;
00333 G4cerr << " p = " << p << G4endl;
00334 }
00335 #endif
00336
00337 if( ( fPtrSolidA->Inside(p) != kOutside) &&
00338 ( fPtrSolidB->Inside(p) != kOutside) )
00339 {
00340 dist= fPtrSolidB->DistanceToOut(p) ;
00341 }
00342 else
00343 {
00344 dist= fPtrSolidA->DistanceToIn(p) ;
00345 }
00346
00347 return dist;
00348 }
00349
00351
00352
00353
00354 G4double
00355 G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p,
00356 const G4ThreeVector& v,
00357 const G4bool calcNorm,
00358 G4bool *validNorm,
00359 G4ThreeVector *n ) const
00360 {
00361 #ifdef G4BOOLDEBUG
00362 if( Inside(p) == kOutside )
00363 {
00364 G4cout << "Position:" << G4endl << G4endl;
00365 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
00366 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
00367 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
00368 G4cout << "Direction:" << G4endl << G4endl;
00369 G4cout << "v.x() = " << v.x() << G4endl;
00370 G4cout << "v.y() = " << v.y() << G4endl;
00371 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
00372 G4cout << "WARNING - Invalid call in "
00373 << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
00374 << " Point p is outside !" << G4endl;
00375 G4cout << " p = " << p << G4endl;
00376 G4cout << " v = " << v << G4endl;
00377 G4cerr << "WARNING - Invalid call in "
00378 << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
00379 << " Point p is outside !" << G4endl;
00380 G4cerr << " p = " << p << G4endl;
00381 G4cerr << " v = " << v << G4endl;
00382 }
00383 #endif
00384
00385 G4double distout;
00386 G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,validNorm,n) ;
00387 G4double distB = fPtrSolidB->DistanceToIn(p,v) ;
00388 if(distB < distA)
00389 {
00390 if(calcNorm)
00391 {
00392 *n = -(fPtrSolidB->SurfaceNormal(p+distB*v)) ;
00393 *validNorm = false ;
00394 }
00395 distout= distB ;
00396 }
00397 else
00398 {
00399 distout= distA ;
00400 }
00401 return distout;
00402 }
00403
00405
00406
00407
00408 G4double
00409 G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p ) const
00410 {
00411 G4double dist=0.0;
00412
00413 if( Inside(p) == kOutside )
00414 {
00415 #ifdef G4BOOLDEBUG
00416 G4cout << "WARNING - Invalid call in "
00417 << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
00418 << " Point p is outside" << G4endl;
00419 G4cout << " p = " << p << G4endl;
00420 G4cerr << "WARNING - Invalid call in "
00421 << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
00422 << " Point p is outside" << G4endl;
00423 G4cerr << " p = " << p << G4endl;
00424 #endif
00425 }
00426 else
00427 {
00428 dist= std::min(fPtrSolidA->DistanceToOut(p),
00429 fPtrSolidB->DistanceToIn(p) ) ;
00430 }
00431 return dist;
00432 }
00433
00435
00436
00437
00438 G4GeometryType G4SubtractionSolid::GetEntityType() const
00439 {
00440 return G4String("G4SubtractionSolid");
00441 }
00442
00444
00445
00446
00447 void
00448 G4SubtractionSolid::ComputeDimensions( G4VPVParameterisation*,
00449 const G4int,
00450 const G4VPhysicalVolume* )
00451 {
00452 }
00453
00455
00456
00457
00458 void
00459 G4SubtractionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
00460 {
00461 scene.AddSolid (*this);
00462 }
00463
00465
00466
00467
00468 G4Polyhedron*
00469 G4SubtractionSolid::CreatePolyhedron () const
00470 {
00471 G4Polyhedron* pA = fPtrSolidA->GetPolyhedron();
00472 G4Polyhedron* pB = fPtrSolidB->GetPolyhedron();
00473 if (pA && pB)
00474 {
00475 G4Polyhedron* resultant = new G4Polyhedron (pA->subtract(*pB));
00476 return resultant;
00477 }
00478 else
00479 {
00480 std::ostringstream oss;
00481 oss << "Solid - " << GetName()
00482 << " - one of the Boolean components has no" << G4endl
00483 << " corresponding polyhedron. Returning NULL !";
00484 G4Exception("G4SubtractionSolid::CreatePolyhedron()", "InvalidSetup",
00485 JustWarning, oss.str().c_str());
00486 return 0;
00487 }
00488 }
00489
00491
00492
00493
00494 G4NURBS*
00495 G4SubtractionSolid::CreateNURBS () const
00496 {
00497
00498
00499 return 0;
00500 }