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 <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
00056
00058
00059
00060
00061
00062 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
00063 G4VSolid* pSolidA ,
00064 G4VSolid* pSolidB )
00065 : G4BooleanSolid(pName,pSolidA,pSolidB)
00066 {
00067 }
00068
00070
00071
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
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
00097
00098
00099 G4SubtractionSolid::G4SubtractionSolid( __void__& a )
00100 : G4BooleanSolid(a)
00101 {
00102 }
00103
00105
00106
00107
00108 G4SubtractionSolid::~G4SubtractionSolid()
00109 {
00110 }
00111
00113
00114
00115
00116 G4SubtractionSolid::G4SubtractionSolid(const G4SubtractionSolid& rhs)
00117 : G4BooleanSolid (rhs)
00118 {
00119 }
00120
00122
00123
00124
00125 G4SubtractionSolid&
00126 G4SubtractionSolid::operator = (const G4SubtractionSolid& rhs)
00127 {
00128
00129
00130 if (this == &rhs) { return *this; }
00131
00132
00133
00134 G4BooleanSolid::operator=(rhs);
00135
00136 return *this;
00137 }
00138
00140
00141
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
00151
00152
00153 return fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
00154 pTransform, pMin, pMax );
00155 }
00156
00158
00159
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
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
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
00278 if ( fPtrSolidB->Inside(p) != kOutside )
00279 {
00280 dist = fPtrSolidB->DistanceToOut(p,v) ;
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 )
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
00329 {
00330 dist = fPtrSolidA->DistanceToIn(p,v) ;
00331
00332 if( dist == kInfinity )
00333 {
00334 return kInfinity ;
00335 }
00336 else
00337 {
00338 G4int count2=0;
00339 while( Inside(p+dist*v) == kOutside )
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)
00349 {
00350 return kInfinity ;
00351 }
00352 dist += disTmp ;
00353 count2++;
00354 if( count2 > 1000 )
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
00388
00389
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) &&
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
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
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
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
00555
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
00570
00571 return 0;
00572 }