Geant4 Cross Reference |
>> 1 // This code implementation is the intellectual property of >> 2 // the GEANT4 collaboration. 1 // 3 // 2 // ******************************************* << 4 // By copying, distributing or modifying the Program (or any work 3 // * License and Disclaimer << 5 // based on the Program) you indicate your acceptance of this statement, 4 // * << 6 // and all its terms. 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 // 7 // 26 // Implementation of the base class for solids << 8 // $Id: G4BooleanSolid.cc,v 1.4 2000/11/20 17:56:01 gcosmo Exp $ >> 9 // GEANT4 tag $Name: geant4-03-00 $ >> 10 // >> 11 // Implementation for the abstract base class for solids created by boolean 27 // operations between other solids 12 // operations between other solids 28 // 13 // 29 // 1998.09.10 V.Grichine - created << 14 // History: 30 // ------------------------------------------- << 15 // >> 16 // 10.09.98 V.Grichine, creation according J. Apostolakis's recommendations 31 17 32 #include "G4BooleanSolid.hh" 18 #include "G4BooleanSolid.hh" 33 #include "G4VSolid.hh" << 34 #include "G4DisplacedSolid.hh" 19 #include "G4DisplacedSolid.hh" 35 #include "G4ReflectedSolid.hh" << 36 #include "G4ScaledSolid.hh" << 37 #include "G4Polyhedron.hh" << 38 #include "HepPolyhedronProcessor.h" << 39 #include "G4QuickRand.hh" << 40 20 41 #include "G4AutoLock.hh" << 21 #include "G4RotationMatrix.hh" 42 << 22 #include "G4ThreeVector.hh" 43 namespace << 23 #include "G4Transform3D.hh" 44 { << 24 #include "G4AffineTransform.hh" 45 G4RecursiveMutex polyhedronMutex = G4MUTEX_I << 46 } << 47 << 48 G4VBooleanProcessor* G4BooleanSolid::fExternal << 49 25 50 ////////////////////////////////////////////// 26 ////////////////////////////////////////////////////////////////// 51 // 27 // 52 // Constructor << 28 // 53 29 54 G4BooleanSolid::G4BooleanSolid( const G4String 30 G4BooleanSolid::G4BooleanSolid( const G4String& pName, 55 G4VSolid* pSol 31 G4VSolid* pSolidA , 56 G4VSolid* pSol << 32 G4VSolid* pSolidB ) : 57 : G4VSolid(pName), fPtrSolidA(pSolidA), fPtr << 33 G4VSolid(pName), >> 34 createdDisplacedSolid(false) 58 { 35 { >> 36 fPtrSolidA = pSolidA ; >> 37 fPtrSolidB = pSolidB ; 59 } 38 } 60 39 61 ////////////////////////////////////////////// 40 ////////////////////////////////////////////////////////////////// 62 // 41 // 63 // Constructor << 42 // 64 43 65 G4BooleanSolid::G4BooleanSolid( const G4String 44 G4BooleanSolid::G4BooleanSolid( const G4String& pName, 66 G4VSolid 45 G4VSolid* pSolidA , 67 G4VSolid 46 G4VSolid* pSolidB , 68 G4Rotati 47 G4RotationMatrix* rotMatrix, 69 const G4ThreeV << 48 const G4ThreeVector& transVector ) : 70 : G4VSolid(pName), createdDisplacedSolid(tru << 49 G4VSolid(pName), >> 50 createdDisplacedSolid(true) 71 { 51 { 72 fPtrSolidA = pSolidA ; 52 fPtrSolidA = pSolidA ; 73 fPtrSolidB = new G4DisplacedSolid("placedB", 53 fPtrSolidB = new G4DisplacedSolid("placedB",pSolidB,rotMatrix,transVector) ; 74 } 54 } 75 55 76 ////////////////////////////////////////////// 56 ////////////////////////////////////////////////////////////////// 77 // 57 // 78 // Constructor << 58 // 79 59 80 G4BooleanSolid::G4BooleanSolid( const G4String 60 G4BooleanSolid::G4BooleanSolid( const G4String& pName, 81 G4VSolid 61 G4VSolid* pSolidA , 82 G4VSolid 62 G4VSolid* pSolidB , 83 const G4Transf << 63 const G4Transform3D& transform ) : 84 : G4VSolid(pName), createdDisplacedSolid(tru << 64 G4VSolid(pName), >> 65 createdDisplacedSolid(true) 85 { 66 { 86 fPtrSolidA = pSolidA ; 67 fPtrSolidA = pSolidA ; 87 fPtrSolidB = new G4DisplacedSolid("placedB", 68 fPtrSolidB = new G4DisplacedSolid("placedB",pSolidB,transform) ; 88 } 69 } 89 70 90 ////////////////////////////////////////////// 71 /////////////////////////////////////////////////////////////// 91 // 72 // 92 // Fake default constructor - sets only member << 73 // Destructor deletes second pointer created by 'new' 93 // for usage restri << 94 << 95 G4BooleanSolid::G4BooleanSolid( __void__& a ) << 96 : G4VSolid(a) << 97 { << 98 } << 99 << 100 ////////////////////////////////////////////// << 101 // << 102 // Destructor deletes transformation contents << 103 74 104 G4BooleanSolid::~G4BooleanSolid() 75 G4BooleanSolid::~G4BooleanSolid() 105 { 76 { 106 if(createdDisplacedSolid) << 77 if(createdDisplacedSolid) delete fPtrSolidB ; 107 { << 108 ((G4DisplacedSolid*)fPtrSolidB)->CleanTran << 109 } << 110 delete fpPolyhedron; fpPolyhedron = nullptr; << 111 } 78 } 112 79 113 ////////////////////////////////////////////// 80 /////////////////////////////////////////////////////////////// 114 // 81 // 115 // Copy constructor << 82 // If Solid is made up from a Boolean operation of two solids, 116 << 83 // return the corresponding solid (for no=0 and 1) 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 84 // If the solid is not a "Boolean", return 0 163 << 164 const G4VSolid* G4BooleanSolid::GetConstituent 85 const G4VSolid* G4BooleanSolid::GetConstituentSolid(G4int no) const 165 { 86 { 166 const G4VSolid* subSolid = nullptr; << 87 const G4VSolid* subSolid=0; 167 if( no == 0 ) 88 if( no == 0 ) 168 subSolid = fPtrSolidA; 89 subSolid = fPtrSolidA; 169 else if( no == 1 ) 90 else if( no == 1 ) 170 subSolid = fPtrSolidB; 91 subSolid = fPtrSolidB; 171 else 92 else 172 { << 93 G4Exception("G4BooleanSolid::GetConstituentSolid()const invalid subsolid index"); 173 DumpInfo(); << 94 174 G4Exception("G4BooleanSolid::GetConstituen << 175 "GeomSolids0002", FatalExcepti << 176 } << 177 return subSolid; 95 return subSolid; 178 } 96 } 179 97 180 ////////////////////////////////////////////// << 98 G4VSolid* G4BooleanSolid::GetConstituentSolid(G4int no) 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 { 99 { 188 G4VSolid* subSolid = nullptr; << 100 G4VSolid* subSolid=0; 189 if( no == 0 ) 101 if( no == 0 ) 190 subSolid = fPtrSolidA; 102 subSolid = fPtrSolidA; 191 else if( no == 1 ) 103 else if( no == 1 ) 192 subSolid = fPtrSolidB; 104 subSolid = fPtrSolidB; 193 else 105 else 194 { << 106 G4Exception("G4BooleanSolid::GetConstituentSolid invalid subsolid index"); 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 107 351 ////////////////////////////////////////////// << 108 return subSolid; 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 } 109 } 575 110 576 111 577 ////////////////////////////////////////////// << 578 // << 579 // Estimate Cubic Volume (capacity) and cache << 580 112 581 G4double G4BooleanSolid::GetCubicVolume() << 582 { << 583 if(fCubicVolume < 0.) << 584 { << 585 fCubicVolume = EstimateCubicVolume(fCubVol << 586 } << 587 return fCubicVolume; << 588 } << 589 113 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 114