Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // Implementation of the base class for solids 27 // operations between other solids 28 // 29 // 1998.09.10 V.Grichine - created 30 // ------------------------------------------- 31 32 #include "G4BooleanSolid.hh" 33 #include "G4VSolid.hh" 34 #include "G4DisplacedSolid.hh" 35 #include "G4ReflectedSolid.hh" 36 #include "G4ScaledSolid.hh" 37 #include "G4Polyhedron.hh" 38 #include "HepPolyhedronProcessor.h" 39 #include "G4QuickRand.hh" 40 41 #include "G4AutoLock.hh" 42 43 namespace 44 { 45 G4RecursiveMutex polyhedronMutex = G4MUTEX_I 46 } 47 48 G4VBooleanProcessor* G4BooleanSolid::fExternal 49 50 ////////////////////////////////////////////// 51 // 52 // Constructor 53 54 G4BooleanSolid::G4BooleanSolid( const G4String 55 G4VSolid* pSol 56 G4VSolid* pSol 57 : G4VSolid(pName), fPtrSolidA(pSolidA), fPtr 58 { 59 } 60 61 ////////////////////////////////////////////// 62 // 63 // Constructor 64 65 G4BooleanSolid::G4BooleanSolid( const G4String 66 G4VSolid 67 G4VSolid 68 G4Rotati 69 const G4ThreeV 70 : G4VSolid(pName), createdDisplacedSolid(tru 71 { 72 fPtrSolidA = pSolidA ; 73 fPtrSolidB = new G4DisplacedSolid("placedB", 74 } 75 76 ////////////////////////////////////////////// 77 // 78 // Constructor 79 80 G4BooleanSolid::G4BooleanSolid( const G4String 81 G4VSolid 82 G4VSolid 83 const G4Transf 84 : G4VSolid(pName), createdDisplacedSolid(tru 85 { 86 fPtrSolidA = pSolidA ; 87 fPtrSolidB = new G4DisplacedSolid("placedB", 88 } 89 90 ////////////////////////////////////////////// 91 // 92 // Fake default constructor - sets only member 93 // for usage restri 94 95 G4BooleanSolid::G4BooleanSolid( __void__& a ) 96 : G4VSolid(a) 97 { 98 } 99 100 ////////////////////////////////////////////// 101 // 102 // Destructor deletes transformation contents 103 104 G4BooleanSolid::~G4BooleanSolid() 105 { 106 if(createdDisplacedSolid) 107 { 108 ((G4DisplacedSolid*)fPtrSolidB)->CleanTran 109 } 110 delete fpPolyhedron; fpPolyhedron = nullptr; 111 } 112 113 ////////////////////////////////////////////// 114 // 115 // Copy constructor 116 117 G4BooleanSolid::G4BooleanSolid(const G4Boolean 118 : G4VSolid (rhs), fPtrSolidA(rhs.fPtrSolidA) 119 fCubicVolume(rhs.fCubicVolume), fSurfaceAr 120 fCubVolStatistics(rhs.fCubVolStatistics), 121 fAreaStatistics(rhs.fAreaStatistics), 122 fCubVolEpsilon(rhs.fCubVolEpsilon), 123 fAreaAccuracy(rhs.fAreaAccuracy), 124 createdDisplacedSolid(rhs.createdDisplaced 125 { 126 fPrimitives.resize(0); fPrimitivesSurfaceAre 127 } 128 129 ////////////////////////////////////////////// 130 // 131 // Assignment operator 132 133 G4BooleanSolid& G4BooleanSolid::operator = (co 134 { 135 // Check assignment to self 136 // 137 if (this == &rhs) { return *this; } 138 139 // Copy base class data 140 // 141 G4VSolid::operator=(rhs); 142 143 // Copy data 144 // 145 fPtrSolidA= rhs.fPtrSolidA; fPtrSolidB= rhs. 146 fCubicVolume= rhs.fCubicVolume; fSurfaceArea 147 fCubVolStatistics = rhs.fCubVolStatistics; f 148 fAreaStatistics = rhs.fAreaStatistics; fArea 149 createdDisplacedSolid= rhs.createdDisplacedS 150 151 fRebuildPolyhedron = false; 152 delete fpPolyhedron; fpPolyhedron = nullptr; 153 fPrimitives.resize(0); fPrimitivesSurfaceAre 154 155 return *this; 156 } 157 158 ////////////////////////////////////////////// 159 // 160 // If solid is made up from a Boolean operatio 161 // return the corresponding solid (for no=0 an 162 // If the solid is not a "Boolean", return 0 163 164 const G4VSolid* G4BooleanSolid::GetConstituent 165 { 166 const G4VSolid* subSolid = nullptr; 167 if( no == 0 ) 168 subSolid = fPtrSolidA; 169 else if( no == 1 ) 170 subSolid = fPtrSolidB; 171 else 172 { 173 DumpInfo(); 174 G4Exception("G4BooleanSolid::GetConstituen 175 "GeomSolids0002", FatalExcepti 176 } 177 return subSolid; 178 } 179 180 ////////////////////////////////////////////// 181 // 182 // If solid is made up from a Boolean operatio 183 // return the corresponding solid (for no=0 an 184 // If the solid is not a "Boolean", return 0 185 186 G4VSolid* G4BooleanSolid::GetConstituentSolid( 187 { 188 G4VSolid* subSolid = nullptr; 189 if( no == 0 ) 190 subSolid = fPtrSolidA; 191 else if( no == 1 ) 192 subSolid = fPtrSolidB; 193 else 194 { 195 DumpInfo(); 196 G4Exception("G4BooleanSolid::GetConstituen 197 "GeomSolids0002", FatalExcepti 198 } 199 return subSolid; 200 } 201 202 ////////////////////////////////////////////// 203 // 204 // Returns entity type 205 206 G4GeometryType G4BooleanSolid::GetEntityType() 207 { 208 return {"G4BooleanSolid"}; 209 } 210 211 ////////////////////////////////////////////// 212 // 213 // Set number of random points to be used for 214 215 void G4BooleanSolid::SetCubVolStatistics(G4int 216 { 217 if (st != fCubVolStatistics) { fCubicVolume 218 fCubVolStatistics = st; 219 220 // Propagate st to all components of the 1st 221 if (fPtrSolidA->GetNumOfConstituents() > 1) 222 { 223 G4VSolid* ptr = fPtrSolidA; 224 while(true) 225 { 226 G4String type = ptr->GetEntityType(); 227 if (type == "G4DisplacedSolid") 228 { 229 ptr = ((G4DisplacedSolid*)ptr)->GetCon 230 continue; 231 } 232 if (type == "G4ReflectedSolid") 233 { 234 ptr = ((G4ReflectedSolid*)ptr)->GetCon 235 continue; 236 } 237 if (type == "G4ScaledSolid") 238 { 239 ptr = ((G4ScaledSolid*)ptr)->GetUnscal 240 continue; 241 } 242 if (type != "G4MultiUnion") // G4MultiUn 243 { 244 ((G4BooleanSolid*)ptr)->SetCubVolStatistics( 245 } 246 break; 247 } 248 } 249 250 // Propagate st to all components of the 2nd 251 if (fPtrSolidB->GetNumOfConstituents() > 1) 252 { 253 G4VSolid* ptr = fPtrSolidB; 254 while(true) 255 { 256 G4String type = ptr->GetEntityType(); 257 if (type == "G4DisplacedSolid") 258 { 259 ptr = ((G4DisplacedSolid*)ptr)->GetCon 260 continue; 261 } 262 if (type == "G4ReflectedSolid") 263 { 264 ptr = ((G4ReflectedSolid*)ptr)->GetCon 265 continue; 266 } 267 if (type == "G4ScaledSolid") 268 { 269 ptr = ((G4ScaledSolid*)ptr)->GetUnscal 270 continue; 271 } 272 if (type != "G4MultiUnion") // G4MultiUn 273 { 274 ((G4BooleanSolid*)ptr)->SetCubVolStatistics( 275 } 276 break; 277 } 278 } 279 } 280 281 ////////////////////////////////////////////// 282 // 283 // Set epsilon for computing cubic volume 284 285 void G4BooleanSolid::SetCubVolEpsilon(G4double 286 { 287 if (ep != fCubVolEpsilon) { fCubicVolume = - 288 fCubVolEpsilon = ep; 289 290 // Propagate ep to all components of the 1st 291 if (fPtrSolidA->GetNumOfConstituents() > 1) 292 { 293 G4VSolid* ptr = fPtrSolidA; 294 while(true) 295 { 296 G4String type = ptr->GetEntityType(); 297 if (type == "G4DisplacedSolid") 298 { 299 ptr = ((G4DisplacedSolid*)ptr)->GetCon 300 continue; 301 } 302 if (type == "G4ReflectedSolid") 303 { 304 ptr = ((G4ReflectedSolid*)ptr)->GetCon 305 continue; 306 } 307 if (type == "G4ScaledSolid") 308 { 309 ptr = ((G4ScaledSolid*)ptr)->GetUnscal 310 continue; 311 } 312 if (type != "G4MultiUnion") // G4MultiUn 313 { 314 ((G4BooleanSolid*)ptr)->SetCubVolEpsilon(ep) 315 } 316 break; 317 } 318 } 319 320 // Propagate ep to all components of the 2nd 321 if (fPtrSolidB->GetNumOfConstituents() > 1) 322 { 323 G4VSolid* ptr = fPtrSolidB; 324 while(true) 325 { 326 G4String type = ptr->GetEntityType(); 327 if (type == "G4DisplacedSolid") 328 { 329 ptr = ((G4DisplacedSolid*)ptr)->GetCon 330 continue; 331 } 332 if (type == "G4ReflectedSolid") 333 { 334 ptr = ((G4ReflectedSolid*)ptr)->GetCon 335 continue; 336 } 337 if (type == "G4ScaledSolid") 338 { 339 ptr = ((G4ScaledSolid*)ptr)->GetUnscal 340 continue; 341 } 342 if (type != "G4MultiUnion") // G4MultiUn 343 { 344 ((G4BooleanSolid*)ptr)->SetCubVolEpsilon(ep) 345 } 346 break; 347 } 348 } 349 } 350 351 ////////////////////////////////////////////// 352 // 353 // Stream object contents to an output stream 354 355 std::ostream& G4BooleanSolid::StreamInfo(std:: 356 { 357 os << "------------------------------------- 358 << " *** Dump for Boolean solid - " << 359 << " ================================= 360 << " Solid type: " << GetEntityType() << 361 << " Parameters of constituent solids: \n 362 << "===================================== 363 fPtrSolidA->StreamInfo(os); 364 fPtrSolidB->StreamInfo(os); 365 os << "===================================== 366 367 return os; 368 } 369 370 ////////////////////////////////////////////// 371 // 372 // Creates list of constituent primitives of a 373 374 void G4BooleanSolid::GetListOfPrimitives( 375 std::vector<std::pair<G4VSolid*,G4Trans 376 const G4Transform3D& curPlacement) cons 377 { 378 G4Transform3D transform; 379 G4VSolid* solid; 380 G4String type; 381 382 // Repeat two times, first time for fPtrSoli 383 // 384 for (auto i=0; i<2; ++i) 385 { 386 transform = curPlacement; 387 solid = (i == 0) ? fPtrSolidA : fPtrSo 388 type = solid->GetEntityType(); 389 390 // While current solid is a trasformed sol 391 // 392 while (type == "G4DisplacedSolid" || 393 type == "G4ReflectedSolid" || 394 type == "G4ScaledSolid") 395 { 396 if (type == "G4DisplacedSolid") 397 { 398 transform = transform * G4Transform3D( 399 ((G4DisplacedSolid*)solid) 400 ((G4DisplacedSolid*)solid) 401 solid = ((G4DisplacedSolid*)solid) 402 } 403 else if (type == "G4ReflectedSolid") 404 { 405 transform= transform*((G4ReflectedSoli 406 solid = ((G4ReflectedSolid*)solid)- 407 } 408 else if (type == "G4ScaledSolid") 409 { 410 transform = transform * ((G4ScaledSoli 411 solid = ((G4ScaledSolid*)solid)->G 412 } 413 type = solid->GetEntityType(); 414 } 415 416 // If current solid is a Boolean solid the 417 // otherwise add it to the list of primiti 418 // 419 if (type == "G4UnionSolid" || 420 type == "G4SubtractionSolid" || 421 type == "G4IntersectionSolid" || 422 type == "G4BooleanSolid") 423 { 424 ((G4BooleanSolid *)solid)->GetListOfPrim 425 } 426 else 427 { 428 primitives.emplace_back(solid,transform) 429 } 430 } 431 } 432 433 ////////////////////////////////////////////// 434 // 435 // Returns a point (G4ThreeVector) randomly an 436 // on the surface of the solid 437 438 G4ThreeVector G4BooleanSolid::GetPointOnSurfac 439 { 440 std::size_t nprims = fPrimitives.size(); 441 std::pair<G4VSolid *, G4Transform3D> prim; 442 443 // Get list of primitives and find the total 444 // 445 if (nprims == 0) 446 { 447 GetListOfPrimitives(fPrimitives, G4Transfo 448 nprims = fPrimitives.size(); 449 fPrimitivesSurfaceArea = 0.; 450 for (std::size_t i=0; i<nprims; ++i) 451 { 452 fPrimitivesSurfaceArea += fPrimitives[i] 453 } 454 } 455 456 // Select random primitive, get random point 457 // check that the point belongs to the surfa 458 // 459 G4ThreeVector p; 460 for (std::size_t k=0; k<100000; ++k) // try 461 { 462 G4double rand = fPrimitivesSurfaceArea * 463 G4double area = 0.; 464 for (std::size_t i=0; i<nprims; ++i) 465 { 466 prim = fPrimitives[i]; 467 area += prim.first->GetSurfaceArea(); 468 if (rand < area) break; 469 } 470 p = prim.first->GetPointOnSurface(); 471 p = prim.second * G4Point3D(p); 472 if (Inside(p) == kSurface) return p; 473 } 474 std::ostringstream message; 475 message << "Solid - " << GetName() << "\n" 476 << "All 100k attempts to generate a 477 << "The solid created may be an inva 478 G4Exception("G4BooleanSolid::GetPointOnSurfa 479 "GeomSolids1001", JustWarning, m 480 return p; 481 } 482 483 ////////////////////////////////////////////// 484 // 485 // Return total number of constituents used fo 486 487 G4int G4BooleanSolid::GetNumOfConstituents() c 488 { 489 return (fPtrSolidA->GetNumOfConstituents() + 490 } 491 492 ////////////////////////////////////////////// 493 // 494 // Return true if the resulting solid has only 495 496 G4bool G4BooleanSolid::IsFaceted() const 497 { 498 return (fPtrSolidA->IsFaceted() && fPtrSolid 499 } 500 501 ////////////////////////////////////////////// 502 // 503 // Returns polyhedron for visualization 504 505 G4Polyhedron* G4BooleanSolid::GetPolyhedron () 506 { 507 if (fpPolyhedron == nullptr || 508 fRebuildPolyhedron || 509 fpPolyhedron->GetNumberOfRotationStepsAt 510 fpPolyhedron->GetNumberOfRotationSteps() 511 { 512 G4RecursiveAutoLock l(&polyhedronMutex); 513 delete fpPolyhedron; 514 fpPolyhedron = CreatePolyhedron(); 515 fRebuildPolyhedron = false; 516 l.unlock(); 517 } 518 return fpPolyhedron; 519 } 520 521 ////////////////////////////////////////////// 522 // 523 // Stacks polyhedra for processing. Returns to 524 525 G4Polyhedron* 526 G4BooleanSolid::StackPolyhedron(HepPolyhedronP 527 const G4VSolid 528 { 529 HepPolyhedronProcessor::Operation operation; 530 const G4String& type = solid->GetEntityType( 531 if (type == "G4UnionSolid") 532 { operation = HepPolyhedronProcessor::UNIO 533 else if (type == "G4IntersectionSolid") 534 { operation = HepPolyhedronProcessor::INTE 535 else if (type == "G4SubtractionSolid") 536 { operation = HepPolyhedronProcessor::SUBT 537 else 538 { 539 std::ostringstream message; 540 message << "Solid - " << solid->GetName() 541 << " - Unrecognised composite soli 542 << " Returning NULL !"; 543 G4Exception("StackPolyhedron()", "GeomSoli 544 return nullptr; 545 } 546 547 G4Polyhedron* top = nullptr; 548 const G4VSolid* solidA = solid->GetConstitue 549 const G4VSolid* solidB = solid->GetConstitue 550 551 if (solidA->GetConstituentSolid(0) != nullpt 552 { 553 top = StackPolyhedron(processor, solidA); 554 } 555 else 556 { 557 top = solidA->GetPolyhedron(); 558 } 559 G4Polyhedron* operand = solidB->GetPolyhedro 560 if (operand != nullptr) 561 { 562 processor.push_back (operation, *operand); 563 } 564 else 565 { 566 std::ostringstream message; 567 message << "Solid - " << solid->GetName() 568 << " - No G4Polyhedron for Boolean 569 G4Exception("G4BooleanSolid::StackPolyhedr 570 "GeomSolids2001", JustWarning, 571 } 572 573 return top; 574 } 575 576 577 ////////////////////////////////////////////// 578 // 579 // Estimate Cubic Volume (capacity) and cache 580 581 G4double G4BooleanSolid::GetCubicVolume() 582 { 583 if(fCubicVolume < 0.) 584 { 585 fCubicVolume = EstimateCubicVolume(fCubVol 586 } 587 return fCubicVolume; 588 } 589 590 ////////////////////////////////////////////// 591 // 592 // Set external Boolean processor. 593 594 void 595 G4BooleanSolid::SetExternalBooleanProcessor(G4 596 { 597 fExternalBoolProcessor = extProcessor; 598 } 599 600 ////////////////////////////////////////////// 601 // 602 // Get external Boolean processor. 603 604 G4VBooleanProcessor* G4BooleanSolid::GetExtern 605 { 606 return fExternalBoolProcessor; 607 } 608