Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // >> 23 // >> 24 // $Id: G4SubtractionSolid.cc,v 1.21 2004/02/27 08:38:09 grichine Exp $ >> 25 // GEANT4 tag $Name: geant4-06-01 $ >> 26 // 26 // Implementation of methods for the class G4I 27 // Implementation of methods for the class G4IntersectionSolid 27 // 28 // 28 // 22.07.11 T.Nikitina: added detection of inf << 29 // History: 29 // 19.10.98 V.Grichine: new algorithm of Dista << 30 // 30 // 14.10.98 V.Grichine: implementation of the 31 // 14.10.98 V.Grichine: implementation of the first version >> 32 // 19.10.98 V.Grichine: new algorithm of DistanceToIn(p,v) >> 33 // 02.08.99 V.Grichine: bugs fixed in DistanceToOut(p,v,...) >> 34 // while -> do-while & surfaceA limitations >> 35 // 13.09.00 V.Grichine: bug fixed in SurfaceNormal(p), p can be inside >> 36 // 31 // ------------------------------------------- 37 // -------------------------------------------------------------------- 32 38 33 #include "G4SubtractionSolid.hh" 39 #include "G4SubtractionSolid.hh" 34 40 35 #include "G4SystemOfUnits.hh" << 36 #include "G4VoxelLimits.hh" 41 #include "G4VoxelLimits.hh" 37 #include "G4VPVParameterisation.hh" 42 #include "G4VPVParameterisation.hh" 38 #include "G4GeometryTolerance.hh" << 39 43 40 #include "G4VGraphicsScene.hh" 44 #include "G4VGraphicsScene.hh" 41 #include "G4Polyhedron.hh" 45 #include "G4Polyhedron.hh" 42 #include "G4PolyhedronArbitrary.hh" << 46 #include "G4NURBS.hh" 43 #include "HepPolyhedronProcessor.h" << 47 // #include "G4NURBSbox.hh" 44 << 45 #include "G4IntersectionSolid.hh" << 46 48 47 #include <sstream> << 49 /////////////////////////////////////////////////////////////////// 48 << 49 ////////////////////////////////////////////// << 50 // 50 // 51 // Transfer all data members to G4BooleanSolid 51 // Transfer all data members to G4BooleanSolid which is responsible 52 // for them. pName will be in turn sent to G4V 52 // for them. pName will be in turn sent to G4VSolid 53 53 54 G4SubtractionSolid::G4SubtractionSolid( const 54 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName, 55 55 G4VSolid* pSolidA , 56 56 G4VSolid* pSolidB ) 57 : G4BooleanSolid(pName,pSolidA,pSolidB) 57 : G4BooleanSolid(pName,pSolidA,pSolidB) 58 { 58 { 59 } 59 } 60 60 61 ////////////////////////////////////////////// << 61 /////////////////////////////////////////////////////////////// 62 // 62 // 63 // Constructor 63 // Constructor 64 64 65 G4SubtractionSolid::G4SubtractionSolid( const 65 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName, 66 66 G4VSolid* pSolidA , 67 67 G4VSolid* pSolidB , 68 68 G4RotationMatrix* rotMatrix, 69 const 69 const G4ThreeVector& transVector ) 70 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMa 70 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector) 71 { 71 { 72 } 72 } 73 73 74 ////////////////////////////////////////////// << 74 /////////////////////////////////////////////////////////////// 75 // 75 // 76 // Constructor 76 // Constructor 77 77 78 G4SubtractionSolid::G4SubtractionSolid( const 78 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName, 79 79 G4VSolid* pSolidA , 80 80 G4VSolid* pSolidB , 81 const 81 const G4Transform3D& transform ) 82 : G4BooleanSolid(pName,pSolidA,pSolidB,trans 82 : G4BooleanSolid(pName,pSolidA,pSolidB,transform) 83 { 83 { 84 } 84 } 85 85 86 ////////////////////////////////////////////// << 87 // << 88 // Fake default constructor - sets only member << 89 // for usage restri << 90 << 91 G4SubtractionSolid::G4SubtractionSolid( __void << 92 : G4BooleanSolid(a) << 93 { << 94 } << 95 86 96 ////////////////////////////////////////////// << 87 /////////////////////////////////////////////////////////////// 97 // 88 // 98 // Destructor 89 // Destructor 99 90 100 G4SubtractionSolid::~G4SubtractionSolid() = de << 91 G4SubtractionSolid::~G4SubtractionSolid() 101 << 102 ////////////////////////////////////////////// << 103 // << 104 // Copy constructor << 105 << 106 G4SubtractionSolid::G4SubtractionSolid(const G << 107 << 108 ////////////////////////////////////////////// << 109 // << 110 // Assignment operator << 111 << 112 G4SubtractionSolid& << 113 G4SubtractionSolid::operator = (const G4Subtra << 114 { << 115 // Check assignment to self << 116 // << 117 if (this == &rhs) { return *this; } << 118 << 119 // Copy base class data << 120 // << 121 G4BooleanSolid::operator=(rhs); << 122 << 123 return *this; << 124 } << 125 << 126 ////////////////////////////////////////////// << 127 // << 128 // Get bounding box << 129 << 130 void << 131 G4SubtractionSolid::BoundingLimits(G4ThreeVect << 132 G4ThreeVect << 133 { 92 { 134 // Since it is unclear how the shape of the << 135 // after subtraction, just return its origin << 136 // << 137 fPtrSolidA->BoundingLimits(pMin,pMax); << 138 << 139 // Check correctness of the bounding box << 140 // << 141 if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 142 { << 143 std::ostringstream message; << 144 message << "Bad bounding box (min >= max) << 145 << GetName() << " !" << 146 << "\npMin = " << pMin << 147 << "\npMax = " << pMax; << 148 G4Exception("G4SubtractionSolid::BoundingL << 149 JustWarning, message); << 150 DumpInfo(); << 151 } << 152 } 93 } 153 94 154 ////////////////////////////////////////////// << 95 /////////////////////////////////////////////////////////////// 155 // 96 // 156 // Calculate extent under transform and specif << 97 // CalculateExtent 157 98 158 G4bool 99 G4bool 159 G4SubtractionSolid::CalculateExtent( const EAx 100 G4SubtractionSolid::CalculateExtent( const EAxis pAxis, 160 const G4V 101 const G4VoxelLimits& pVoxelLimit, 161 const G4A 102 const G4AffineTransform& pTransform, 162 G4d 103 G4double& pMin, 163 G4d 104 G4double& pMax ) const 164 { 105 { 165 // Since we cannot be sure how much the seco 106 // Since we cannot be sure how much the second solid subtracts 166 // from the first, we must use the first sol << 107 // from the first, we must use the first solid's extent! 167 108 168 return fPtrSolidA->CalculateExtent( pAxis, p 109 return fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit, 169 pTransfo 110 pTransform, pMin, pMax ); 170 } 111 } 171 112 172 ////////////////////////////////////////////// << 113 ///////////////////////////////////////////////////// 173 // 114 // 174 // Touching ? Empty subtraction ? 115 // Touching ? Empty subtraction ? 175 116 176 EInside G4SubtractionSolid::Inside( const G4Th 117 EInside G4SubtractionSolid::Inside( const G4ThreeVector& p ) const 177 { 118 { 178 EInside positionA = fPtrSolidA->Inside(p); << 119 EInside positionA = fPtrSolidA->Inside(p) ; 179 if (positionA == kOutside) return positionA; << 120 EInside positionB = fPtrSolidB->Inside(p) ; 180 << 121 181 EInside positionB = fPtrSolidB->Inside(p); << 122 if(positionA == kInside && positionB == kOutside) 182 if (positionB == kOutside) return positionA; << 123 { 183 << 124 return kInside ; 184 if (positionB == kInside) return kOutside; << 125 } 185 if (positionA == kInside) return kSurface; / << 126 else 186 << 127 { 187 // Point is on both surfaces << 128 if((positionA == kInside && positionB == kSurface) || 188 // << 129 (positionB == kOutside && positionA == kSurface) || 189 static const G4double rtol = 1000*kCarTolera << 130 (positionA == kSurface && positionB == kSurface) ) 190 << 131 { 191 return ((fPtrSolidA->SurfaceNormal(p) - << 132 return kSurface ; 192 fPtrSolidB->SurfaceNormal(p)).mag2( << 133 } >> 134 else >> 135 { >> 136 return kOutside ; >> 137 } >> 138 } 193 } 139 } 194 140 195 ////////////////////////////////////////////// << 141 ////////////////////////////////////////////////////////////// 196 // 142 // 197 // SurfaceNormal 143 // SurfaceNormal 198 144 199 G4ThreeVector 145 G4ThreeVector 200 G4SubtractionSolid::SurfaceNormal( const G4Thr 146 G4SubtractionSolid::SurfaceNormal( const G4ThreeVector& p ) const 201 { 147 { 202 G4ThreeVector normal; 148 G4ThreeVector normal; 203 << 149 if( Inside(p) == kOutside ) 204 EInside InsideA = fPtrSolidA->Inside(p); << 205 EInside InsideB = fPtrSolidB->Inside(p); << 206 << 207 if( InsideA == kOutside ) << 208 { 150 { 209 #ifdef G4BOOLDEBUG 151 #ifdef G4BOOLDEBUG 210 G4cout << "WARNING - Invalid call [1] in " 152 G4cout << "WARNING - Invalid call [1] in " 211 << "G4SubtractionSolid::SurfaceNorm 153 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl 212 << " Point p is outside !" << G4en << 154 << " Point p is inside !" << G4endl; 213 G4cout << " p = " << p << G4endl; 155 G4cout << " p = " << p << G4endl; 214 G4cerr << "WARNING - Invalid call [1] in " 156 G4cerr << "WARNING - Invalid call [1] in " 215 << "G4SubtractionSolid::SurfaceNorm 157 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl 216 << " Point p is outside !" << G4en << 158 << " Point p is inside !" << G4endl; 217 G4cerr << " p = " << p << G4endl; 159 G4cerr << " p = " << p << G4endl; 218 #endif 160 #endif 219 normal = fPtrSolidA->SurfaceNormal(p) ; << 220 } << 221 else if( InsideA == kSurface && << 222 InsideB != kInside ) << 223 { << 224 normal = fPtrSolidA->SurfaceNormal(p) ; << 225 } 161 } 226 else if( InsideA == kInside && << 162 else 227 InsideB != kOutside ) << 163 { 228 { << 164 if( fPtrSolidA->Inside(p) == kSurface && 229 normal = -fPtrSolidB->SurfaceNormal(p) ; << 165 fPtrSolidB->Inside(p) != kInside ) 230 } << 231 else << 232 { << 233 if ( fPtrSolidA->DistanceToOut(p) <= fPtrS << 234 { 166 { 235 normal = fPtrSolidA->SurfaceNormal(p) ; 167 normal = fPtrSolidA->SurfaceNormal(p) ; 236 } 168 } 237 else << 169 else if( fPtrSolidA->Inside(p) == kInside && >> 170 fPtrSolidB->Inside(p) != kOutside ) 238 { 171 { 239 normal = -fPtrSolidB->SurfaceNormal(p) ; 172 normal = -fPtrSolidB->SurfaceNormal(p) ; 240 } 173 } 241 #ifdef G4BOOLDEBUG << 174 else 242 if(Inside(p) == kInside) << 243 { 175 { 244 G4cout << "WARNING - Invalid call [2] in << 176 if ( fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToIn(p) ) >> 177 { >> 178 normal = fPtrSolidA->SurfaceNormal(p) ; >> 179 } >> 180 else >> 181 { >> 182 normal = -fPtrSolidB->SurfaceNormal(p) ; >> 183 } >> 184 #ifdef G4BOOLDEBUG >> 185 if(Inside(p) == kInside) >> 186 { >> 187 G4cout << "WARNING - Invalid call [2] in " 245 << "G4SubtractionSolid::SurfaceNo 188 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl 246 << " Point p is inside !" << G4e 189 << " Point p is inside !" << G4endl; 247 G4cout << " p = " << p << G4end << 190 G4cout << " p = " << p << G4endl; 248 G4cerr << "WARNING - Invalid call [2] in << 191 G4cerr << "WARNING - Invalid call [2] in " 249 << "G4SubtractionSolid::SurfaceNo 192 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl 250 << " Point p is inside !" << G4e 193 << " Point p is inside !" << G4endl; 251 G4cerr << " p = " << p << G4end << 194 G4cerr << " p = " << p << G4endl; 252 } << 195 } 253 #endif 196 #endif >> 197 } 254 } 198 } 255 return normal; 199 return normal; 256 } 200 } 257 201 258 ////////////////////////////////////////////// << 202 ///////////////////////////////////////////////////////////// 259 // 203 // 260 // The same algorithm as in DistanceToIn(p) 204 // The same algorithm as in DistanceToIn(p) 261 205 262 G4double 206 G4double 263 G4SubtractionSolid::DistanceToIn( const G4Thre << 207 G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p, 264 const G4Thre << 208 const G4ThreeVector& v ) const 265 { 209 { 266 G4double dist = 0.0, dist2 = 0.0, disTmp = 0 << 210 G4double dist = 0.0,disTmp = 0.0 ; 267 211 268 #ifdef G4BOOLDEBUG 212 #ifdef G4BOOLDEBUG 269 if( Inside(p) == kInside ) 213 if( Inside(p) == kInside ) 270 { 214 { 271 G4cout << "WARNING - Invalid call in " 215 G4cout << "WARNING - Invalid call in " 272 << "G4SubtractionSolid::DistanceToI 216 << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl 273 << " Point p is inside !" << G4end 217 << " Point p is inside !" << G4endl; 274 G4cout << " p = " << p << G4endl; 218 G4cout << " p = " << p << G4endl; 275 G4cout << " v = " << v << G4endl; 219 G4cout << " v = " << v << G4endl; 276 G4cerr << "WARNING - Invalid call in " 220 G4cerr << "WARNING - Invalid call in " 277 << "G4SubtractionSolid::DistanceToI 221 << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl 278 << " Point p is inside !" << G4end 222 << " Point p is inside !" << G4endl; 279 G4cerr << " p = " << p << G4endl; 223 G4cerr << " p = " << p << G4endl; 280 G4cerr << " v = " << v << G4endl; 224 G4cerr << " v = " << v << G4endl; 281 } 225 } 282 #endif 226 #endif 283 227 284 // if( // ( fPtrSolidA->Inside(p) != kOuts 228 // if( // ( fPtrSolidA->Inside(p) != kOutside) && // case1:p in both A&B 285 if ( fPtrSolidB->Inside(p) != kOutside ) 229 if ( fPtrSolidB->Inside(p) != kOutside ) // start: out of B 286 { 230 { 287 dist = fPtrSolidB->DistanceToOut(p,v) ; 231 dist = fPtrSolidB->DistanceToOut(p,v) ; // ,calcNorm,validNorm,n) ; 288 232 289 if( fPtrSolidA->Inside(p+dist*v) != kIns 233 if( fPtrSolidA->Inside(p+dist*v) != kInside ) 290 { 234 { 291 G4int count1=0; << 235 do 292 do // Loop checking, 13.08.2015, G.C << 293 { 236 { 294 disTmp = fPtrSolidA->DistanceToIn(p+ 237 disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ; 295 238 296 if(disTmp == kInfinity) 239 if(disTmp == kInfinity) 297 { 240 { 298 return kInfinity ; 241 return kInfinity ; 299 } 242 } 300 dist += disTmp ; 243 dist += disTmp ; 301 244 302 if( Inside(p+dist*v) == kOutside ) 245 if( Inside(p+dist*v) == kOutside ) 303 { 246 { 304 disTmp = fPtrSolidB->DistanceToOut << 247 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ; 305 dist2 = dist+disTmp; << 248 dist += disTmp ; 306 if (dist == dist2) { return dist; << 249 } 307 dist = dist2 ; << 308 ++count1; << 309 if( count1 > 1000 ) // Infinite l << 310 { << 311 G4String nameB = fPtrSolidB->Get << 312 if(fPtrSolidB->GetEntityType()== << 313 { << 314 nameB = (dynamic_cast<G4Displa << 315 ->GetConstituentMovedS << 316 } << 317 std::ostringstream message; << 318 message << "Illegal condition ca << 319 << fPtrSolidA->GetName() << 320 message.precision(16); << 321 message << "Looping detected in << 322 << ", from original poin << 323 << " and direction " << << 324 << "Computed candidate d << 325 message.precision(6); << 326 DumpInfo(); << 327 G4Exception("G4SubtractionSolid: << 328 "GeomSolids1001", Ju << 329 "Returning candidate << 330 return dist; << 331 } << 332 } << 333 } 250 } 334 while( Inside(p+dist*v) == kOutside ) 251 while( Inside(p+dist*v) == kOutside ) ; 335 } 252 } 336 } 253 } 337 else // p outside A, start in A 254 else // p outside A, start in A 338 { 255 { 339 dist = fPtrSolidA->DistanceToIn(p,v) ; 256 dist = fPtrSolidA->DistanceToIn(p,v) ; 340 257 341 if( dist == kInfinity ) // past A, hence 258 if( dist == kInfinity ) // past A, hence past A\B 342 { 259 { 343 return kInfinity ; 260 return kInfinity ; 344 } 261 } 345 else 262 else 346 { 263 { 347 G4int count2=0; << 348 while( Inside(p+dist*v) == kOutside ) 264 while( Inside(p+dist*v) == kOutside ) // pushing loop 349 { 265 { 350 disTmp = fPtrSolidB->DistanceToOut(p 266 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ; 351 dist += disTmp ; 267 dist += disTmp ; 352 268 353 if( Inside(p+dist*v) == kOutside ) 269 if( Inside(p+dist*v) == kOutside ) 354 { 270 { 355 disTmp = fPtrSolidA->DistanceToIn( 271 disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ; 356 272 357 if(disTmp == kInfinity) // past A, 273 if(disTmp == kInfinity) // past A, hence past A\B 358 { 274 { 359 return kInfinity ; 275 return kInfinity ; 360 } 276 } 361 dist2 = dist+disTmp; << 277 dist += disTmp ; 362 if (dist == dist2) { return dist; << 363 dist = dist2 ; << 364 ++count2; << 365 if( count2 > 1000 ) // Infinite l << 366 { << 367 G4String nameB = fPtrSolidB->Get << 368 if(fPtrSolidB->GetEntityType()== << 369 { << 370 nameB = (dynamic_cast<G4Displa << 371 ->GetConstituentMovedS << 372 } << 373 std::ostringstream message; << 374 message << "Illegal condition ca << 375 << fPtrSolidA->GetName() << 376 message.precision(16); << 377 message << "Looping detected in << 378 << ", from original poin << 379 << " and direction " << << 380 << "Computed candidate d << 381 message.precision(6); << 382 DumpInfo(); << 383 G4Exception("G4SubtractionSolid: << 384 "GeomSolids1001", Ju << 385 "Returning candidate << 386 return dist; << 387 } << 388 } 278 } 389 } // Loop checking, 13.08.2015, G.C << 279 } 390 } 280 } 391 } 281 } 392 282 393 return dist ; 283 return dist ; 394 } 284 } 395 285 396 ////////////////////////////////////////////// << 286 //////////////////////////////////////////////////////// 397 // 287 // 398 // Approximate nearest distance from the point 288 // Approximate nearest distance from the point p to the intersection of 399 // two solids. It is usually underestimated fr 289 // two solids. It is usually underestimated from the point of view of 400 // isotropic safety 290 // isotropic safety 401 291 402 G4double 292 G4double 403 G4SubtractionSolid::DistanceToIn( const G4Thre 293 G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p ) const 404 { 294 { 405 G4double dist = 0.0; << 295 G4double dist=0.0; 406 296 407 #ifdef G4BOOLDEBUG 297 #ifdef G4BOOLDEBUG 408 if( Inside(p) == kInside ) 298 if( Inside(p) == kInside ) 409 { 299 { 410 G4cout << "WARNING - Invalid call in " 300 G4cout << "WARNING - Invalid call in " 411 << "G4SubtractionSolid::DistanceToI 301 << "G4SubtractionSolid::DistanceToIn(p)" << G4endl 412 << " Point p is inside !" << G4end 302 << " Point p is inside !" << G4endl; 413 G4cout << " p = " << p << G4endl; 303 G4cout << " p = " << p << G4endl; 414 G4cerr << "WARNING - Invalid call in " 304 G4cerr << "WARNING - Invalid call in " 415 << "G4SubtractionSolid::DistanceToI 305 << "G4SubtractionSolid::DistanceToIn(p)" << G4endl 416 << " Point p is inside !" << G4end 306 << " Point p is inside !" << G4endl; 417 G4cerr << " p = " << p << G4endl; 307 G4cerr << " p = " << p << G4endl; 418 } 308 } 419 #endif 309 #endif 420 310 421 if( ( fPtrSolidA->Inside(p) != kOutside) && 311 if( ( fPtrSolidA->Inside(p) != kOutside) && // case 1 422 ( fPtrSolidB->Inside(p) != kOutside) 312 ( fPtrSolidB->Inside(p) != kOutside) ) 423 { 313 { 424 dist = fPtrSolidB->DistanceToOut(p); << 314 dist= fPtrSolidB->DistanceToOut(p) ; 425 } 315 } 426 else 316 else 427 { 317 { 428 dist = fPtrSolidA->DistanceToIn(p); << 318 dist= fPtrSolidA->DistanceToIn(p) ; 429 } 319 } 430 320 431 return dist; 321 return dist; 432 } 322 } 433 323 434 ////////////////////////////////////////////// << 324 ////////////////////////////////////////////////////////// 435 // 325 // 436 // The same algorithm as DistanceToOut(p) 326 // The same algorithm as DistanceToOut(p) 437 327 438 G4double 328 G4double 439 G4SubtractionSolid::DistanceToOut( const G4Thr 329 G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p, 440 const G4Thr << 330 const G4ThreeVector& v, 441 const G4boo << 331 const G4bool calcNorm, 442 G4boo << 332 G4bool *validNorm, 443 G4Thr << 333 G4ThreeVector *n ) const 444 { 334 { 445 #ifdef G4BOOLDEBUG 335 #ifdef G4BOOLDEBUG 446 if( Inside(p) == kOutside ) 336 if( Inside(p) == kOutside ) 447 { 337 { 448 G4cout << "Position:" << G4endl << G4en 338 G4cout << "Position:" << G4endl << G4endl; 449 G4cout << "p.x() = " << p.x()/mm << " 339 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl; 450 G4cout << "p.y() = " << p.y()/mm << " 340 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl; 451 G4cout << "p.z() = " << p.z()/mm << " 341 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl; 452 G4cout << "Direction:" << G4endl << G4en 342 G4cout << "Direction:" << G4endl << G4endl; 453 G4cout << "v.x() = " << v.x() << G4end 343 G4cout << "v.x() = " << v.x() << G4endl; 454 G4cout << "v.y() = " << v.y() << G4end 344 G4cout << "v.y() = " << v.y() << G4endl; 455 G4cout << "v.z() = " << v.z() << G4end 345 G4cout << "v.z() = " << v.z() << G4endl << G4endl; 456 G4cout << "WARNING - Invalid call in " 346 G4cout << "WARNING - Invalid call in " 457 << "G4SubtractionSolid::DistanceT 347 << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl 458 << " Point p is outside !" << G4 348 << " Point p is outside !" << G4endl; 459 G4cout << " p = " << p << G4end 349 G4cout << " p = " << p << G4endl; 460 G4cout << " v = " << v << G4end 350 G4cout << " v = " << v << G4endl; 461 G4cerr << "WARNING - Invalid call in " 351 G4cerr << "WARNING - Invalid call in " 462 << "G4SubtractionSolid::DistanceT 352 << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl 463 << " Point p is outside !" << G4 353 << " Point p is outside !" << G4endl; 464 G4cerr << " p = " << p << G4end 354 G4cerr << " p = " << p << G4endl; 465 G4cerr << " v = " << v << G4end 355 G4cerr << " v = " << v << G4endl; 466 } 356 } 467 #endif 357 #endif 468 358 469 G4double distout; 359 G4double distout; 470 G4double distA = fPtrSolidA->DistanceToOut 360 G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,validNorm,n) ; 471 G4double distB = fPtrSolidB->DistanceToIn( 361 G4double distB = fPtrSolidB->DistanceToIn(p,v) ; 472 if(distB < distA) 362 if(distB < distA) 473 { 363 { 474 if(calcNorm) 364 if(calcNorm) 475 { 365 { 476 *n = -(fPtrSolidB->SurfaceNormal(p+dis 366 *n = -(fPtrSolidB->SurfaceNormal(p+distB*v)) ; 477 *validNorm = false ; 367 *validNorm = false ; 478 } 368 } 479 distout= distB ; 369 distout= distB ; 480 } 370 } 481 else 371 else 482 { 372 { 483 distout= distA ; 373 distout= distA ; 484 } 374 } 485 return distout; 375 return distout; 486 } 376 } 487 377 488 ////////////////////////////////////////////// << 378 ////////////////////////////////////////////////////////////// 489 // 379 // 490 // Inverted algorithm of DistanceToIn(p) 380 // Inverted algorithm of DistanceToIn(p) 491 381 492 G4double 382 G4double 493 G4SubtractionSolid::DistanceToOut( const G4Thr 383 G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p ) const 494 { 384 { 495 G4double dist=0.0; 385 G4double dist=0.0; 496 386 497 if( Inside(p) == kOutside ) 387 if( Inside(p) == kOutside ) 498 { 388 { 499 #ifdef G4BOOLDEBUG 389 #ifdef G4BOOLDEBUG 500 G4cout << "WARNING - Invalid call in " 390 G4cout << "WARNING - Invalid call in " 501 << "G4SubtractionSolid::DistanceToO 391 << "G4SubtractionSolid::DistanceToOut(p)" << G4endl 502 << " Point p is outside" << G4endl 392 << " Point p is outside" << G4endl; 503 G4cout << " p = " << p << G4endl; 393 G4cout << " p = " << p << G4endl; 504 G4cerr << "WARNING - Invalid call in " 394 G4cerr << "WARNING - Invalid call in " 505 << "G4SubtractionSolid::DistanceToO 395 << "G4SubtractionSolid::DistanceToOut(p)" << G4endl 506 << " Point p is outside" << G4endl 396 << " Point p is outside" << G4endl; 507 G4cerr << " p = " << p << G4endl; 397 G4cerr << " p = " << p << G4endl; 508 #endif 398 #endif 509 } 399 } 510 else 400 else 511 { 401 { 512 dist= std::min(fPtrSolidA->DistanceToOut( 402 dist= std::min(fPtrSolidA->DistanceToOut(p), 513 fPtrSolidB->DistanceToIn 403 fPtrSolidB->DistanceToIn(p) ) ; 514 } 404 } 515 return dist; 405 return dist; 516 } 406 } 517 407 518 ////////////////////////////////////////////// << 408 ////////////////////////////////////////////////////////////// 519 // 409 // 520 // 410 // 521 411 522 G4GeometryType G4SubtractionSolid::GetEntityTy 412 G4GeometryType G4SubtractionSolid::GetEntityType() const 523 { 413 { 524 return {"G4SubtractionSolid"}; << 414 return G4String("G4SubtractionSolid"); 525 } 415 } 526 416 527 ////////////////////////////////////////////// << 417 ////////////////////////////////////////////////////////////// 528 // 418 // 529 // Make a clone of the object << 530 << 531 G4VSolid* G4SubtractionSolid::Clone() const << 532 { << 533 return new G4SubtractionSolid(*this); << 534 } << 535 << 536 ////////////////////////////////////////////// << 537 // 419 // 538 // ComputeDimensions << 539 420 540 void 421 void 541 G4SubtractionSolid::ComputeDimensions( G 422 G4SubtractionSolid::ComputeDimensions( G4VPVParameterisation*, 542 const G 423 const G4int, 543 const G 424 const G4VPhysicalVolume* ) 544 { 425 { 545 } 426 } 546 427 547 ////////////////////////////////////////////// << 428 ///////////////////////////////////////////////// 548 // 429 // 549 // DescribeYourselfTo << 430 // 550 431 551 void 432 void 552 G4SubtractionSolid::DescribeYourselfTo ( G4VGr 433 G4SubtractionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 553 { 434 { 554 scene.AddSolid (*this); << 435 scene.AddThis (*this); 555 } 436 } 556 437 557 ////////////////////////////////////////////// << 438 //////////////////////////////////////////////////// >> 439 // 558 // 440 // 559 // CreatePolyhedron << 560 441 561 G4Polyhedron* G4SubtractionSolid::CreatePolyhe << 442 G4Polyhedron* >> 443 G4SubtractionSolid::CreatePolyhedron () const 562 { 444 { 563 if (fExternalBoolProcessor == nullptr) << 445 G4Polyhedron* pA = fPtrSolidA->CreatePolyhedron(); 564 { << 446 G4Polyhedron* pB = fPtrSolidB->CreatePolyhedron(); 565 HepPolyhedronProcessor processor; << 447 G4Polyhedron* resultant = new G4Polyhedron (pA->subtract(*pB)); 566 // Stack components and components of comp << 448 delete pB; 567 // See G4BooleanSolid::StackPolyhedron << 449 delete pA; 568 G4Polyhedron* top = StackPolyhedron(proces << 450 return resultant; 569 auto result = new G4Polyhedron(*top); << 570 if (processor.execute(*result)) << 571 { << 572 return result; << 573 } << 574 else << 575 { << 576 return nullptr; << 577 } << 578 } << 579 else << 580 { << 581 return fExternalBoolProcessor->Process(thi << 582 } << 583 } 451 } 584 452 585 ////////////////////////////////////////////// << 453 ///////////////////////////////////////////////////////// 586 // 454 // 587 // GetCubicVolume << 588 // 455 // 589 456 590 G4double G4SubtractionSolid::GetCubicVolume() << 457 G4NURBS* >> 458 G4SubtractionSolid::CreateNURBS () const 591 { 459 { 592 if( fCubicVolume >= 0. ) << 460 // Take into account boolean operation - see CreatePolyhedron. 593 { << 461 // return new G4NURBSbox (1.0, 1.0, 1.0); 594 return fCubicVolume; << 462 return 0; 595 } << 596 G4ThreeVector bminA, bmaxA, bminB, bmaxB; << 597 fPtrSolidA->BoundingLimits(bminA, bmaxA); << 598 fPtrSolidB->BoundingLimits(bminB, bmaxB); << 599 G4bool noIntersection = << 600 bminA.x() >= bmaxB.x() || bminA.y() >= bm << 601 bminB.x() >= bmaxA.x() || bminB.y() >= bm << 602 << 603 if (noIntersection) << 604 { << 605 fCubicVolume = fPtrSolidA->GetCubicVolume( << 606 } << 607 else << 608 { << 609 if (GetNumOfConstituents() > 10) << 610 { << 611 fCubicVolume = G4BooleanSolid::GetCubicV << 612 } << 613 else << 614 { << 615 G4IntersectionSolid intersectVol("Tempor << 616 fPtrSo << 617 intersectVol.SetCubVolStatistics(GetCubV << 618 intersectVol.SetCubVolEpsilon(GetCubVolE << 619 << 620 G4double cubVolumeA = fPtrSolidA->GetCub << 621 fCubicVolume = cubVolumeA - intersectVol << 622 if (fCubicVolume < 0.01*cubVolumeA) fCub << 623 } << 624 } << 625 return fCubicVolume; << 626 } 463 } 627 464