| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Implementation of the base class for solids created by Boolean
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_INITIALIZER;
46 }
47
48 G4VBooleanProcessor* G4BooleanSolid::fExternalBoolProcessor = nullptr;
49
50 //////////////////////////////////////////////////////////////////
51 //
52 // Constructor
53
54 G4BooleanSolid::G4BooleanSolid( const G4String& pName,
55 G4VSolid* pSolidA ,
56 G4VSolid* pSolidB )
57 : G4VSolid(pName), fPtrSolidA(pSolidA), fPtrSolidB(pSolidB)
58 {
59 }
60
61 //////////////////////////////////////////////////////////////////
62 //
63 // Constructor
64
65 G4BooleanSolid::G4BooleanSolid( const G4String& pName,
66 G4VSolid* pSolidA ,
67 G4VSolid* pSolidB ,
68 G4RotationMatrix* rotMatrix,
69 const G4ThreeVector& transVector )
70 : G4VSolid(pName), createdDisplacedSolid(true)
71 {
72 fPtrSolidA = pSolidA ;
73 fPtrSolidB = new G4DisplacedSolid("placedB",pSolidB,rotMatrix,transVector) ;
74 }
75
76 //////////////////////////////////////////////////////////////////
77 //
78 // Constructor
79
80 G4BooleanSolid::G4BooleanSolid( const G4String& pName,
81 G4VSolid* pSolidA ,
82 G4VSolid* pSolidB ,
83 const G4Transform3D& transform )
84 : G4VSolid(pName), createdDisplacedSolid(true)
85 {
86 fPtrSolidA = pSolidA ;
87 fPtrSolidB = new G4DisplacedSolid("placedB",pSolidB,transform) ;
88 }
89
90 ///////////////////////////////////////////////////////////////
91 //
92 // Fake default constructor - sets only member data and allocates memory
93 // for usage restricted to object persistency.
94
95 G4BooleanSolid::G4BooleanSolid( __void__& a )
96 : G4VSolid(a)
97 {
98 }
99
100 ///////////////////////////////////////////////////////////////
101 //
102 // Destructor deletes transformation contents of the created displaced solid
103
104 G4BooleanSolid::~G4BooleanSolid()
105 {
106 if(createdDisplacedSolid)
107 {
108 ((G4DisplacedSolid*)fPtrSolidB)->CleanTransformations();
109 }
110 delete fpPolyhedron; fpPolyhedron = nullptr;
111 }
112
113 ///////////////////////////////////////////////////////////////
114 //
115 // Copy constructor
116
117 G4BooleanSolid::G4BooleanSolid(const G4BooleanSolid& rhs)
118 : G4VSolid (rhs), fPtrSolidA(rhs.fPtrSolidA), fPtrSolidB(rhs.fPtrSolidB),
119 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
120 fCubVolStatistics(rhs.fCubVolStatistics),
121 fAreaStatistics(rhs.fAreaStatistics),
122 fCubVolEpsilon(rhs.fCubVolEpsilon),
123 fAreaAccuracy(rhs.fAreaAccuracy),
124 createdDisplacedSolid(rhs.createdDisplacedSolid)
125 {
126 fPrimitives.resize(0); fPrimitivesSurfaceArea = 0.;
127 }
128
129 ///////////////////////////////////////////////////////////////
130 //
131 // Assignment operator
132
133 G4BooleanSolid& G4BooleanSolid::operator = (const G4BooleanSolid& rhs)
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.fPtrSolidB;
146 fCubicVolume= rhs.fCubicVolume; fSurfaceArea= rhs.fSurfaceArea;
147 fCubVolStatistics = rhs.fCubVolStatistics; fCubVolEpsilon = rhs.fCubVolEpsilon;
148 fAreaStatistics = rhs.fAreaStatistics; fAreaAccuracy = rhs.fAreaAccuracy;
149 createdDisplacedSolid= rhs.createdDisplacedSolid;
150
151 fRebuildPolyhedron = false;
152 delete fpPolyhedron; fpPolyhedron = nullptr;
153 fPrimitives.resize(0); fPrimitivesSurfaceArea = 0.;
154
155 return *this;
156 }
157
158 ///////////////////////////////////////////////////////////////
159 //
160 // If solid is made up from a Boolean operation of two solids,
161 // return the corresponding solid (for no=0 and 1)
162 // If the solid is not a "Boolean", return 0
163
164 const G4VSolid* G4BooleanSolid::GetConstituentSolid(G4int no) const
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::GetConstituentSolid()",
175 "GeomSolids0002", FatalException, "Invalid solid index.");
176 }
177 return subSolid;
178 }
179
180 ///////////////////////////////////////////////////////////////
181 //
182 // If solid is made up from a Boolean operation of two solids,
183 // return the corresponding solid (for no=0 and 1)
184 // If the solid is not a "Boolean", return 0
185
186 G4VSolid* G4BooleanSolid::GetConstituentSolid(G4int no)
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::GetConstituentSolid()",
197 "GeomSolids0002", FatalException, "Invalid solid index.");
198 }
199 return subSolid;
200 }
201
202 //////////////////////////////////////////////////////////////////////////
203 //
204 // Returns entity type
205
206 G4GeometryType G4BooleanSolid::GetEntityType() const
207 {
208 return {"G4BooleanSolid"};
209 }
210
211 //////////////////////////////////////////////////////////////////////////
212 //
213 // Set number of random points to be used for computing cubic volume
214
215 void G4BooleanSolid::SetCubVolStatistics(G4int st)
216 {
217 if (st != fCubVolStatistics) { fCubicVolume = -1.; }
218 fCubVolStatistics = st;
219
220 // Propagate st to all components of the 1st solid
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)->GetConstituentMovedSolid();
230 continue;
231 }
232 if (type == "G4ReflectedSolid")
233 {
234 ptr = ((G4ReflectedSolid*)ptr)->GetConstituentMovedSolid();
235 continue;
236 }
237 if (type == "G4ScaledSolid")
238 {
239 ptr = ((G4ScaledSolid*)ptr)->GetUnscaledSolid();
240 continue;
241 }
242 if (type != "G4MultiUnion") // G4MultiUnion doesn't have SetCubVolStatistics()
243 {
244 ((G4BooleanSolid*)ptr)->SetCubVolStatistics(st);
245 }
246 break;
247 }
248 }
249
250 // Propagate st to all components of the 2nd solid
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)->GetConstituentMovedSolid();
260 continue;
261 }
262 if (type == "G4ReflectedSolid")
263 {
264 ptr = ((G4ReflectedSolid*)ptr)->GetConstituentMovedSolid();
265 continue;
266 }
267 if (type == "G4ScaledSolid")
268 {
269 ptr = ((G4ScaledSolid*)ptr)->GetUnscaledSolid();
270 continue;
271 }
272 if (type != "G4MultiUnion") // G4MultiUnion doesn't have SetCubVolStatistics()
273 {
274 ((G4BooleanSolid*)ptr)->SetCubVolStatistics(st);
275 }
276 break;
277 }
278 }
279 }
280
281 //////////////////////////////////////////////////////////////////////////
282 //
283 // Set epsilon for computing cubic volume
284
285 void G4BooleanSolid::SetCubVolEpsilon(G4double ep)
286 {
287 if (ep != fCubVolEpsilon) { fCubicVolume = -1.; }
288 fCubVolEpsilon = ep;
289
290 // Propagate ep to all components of the 1st solid
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)->GetConstituentMovedSolid();
300 continue;
301 }
302 if (type == "G4ReflectedSolid")
303 {
304 ptr = ((G4ReflectedSolid*)ptr)->GetConstituentMovedSolid();
305 continue;
306 }
307 if (type == "G4ScaledSolid")
308 {
309 ptr = ((G4ScaledSolid*)ptr)->GetUnscaledSolid();
310 continue;
311 }
312 if (type != "G4MultiUnion") // G4MultiUnion doesn't have SetCubVolEpsilon()
313 {
314 ((G4BooleanSolid*)ptr)->SetCubVolEpsilon(ep);
315 }
316 break;
317 }
318 }
319
320 // Propagate ep to all components of the 2nd solid
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)->GetConstituentMovedSolid();
330 continue;
331 }
332 if (type == "G4ReflectedSolid")
333 {
334 ptr = ((G4ReflectedSolid*)ptr)->GetConstituentMovedSolid();
335 continue;
336 }
337 if (type == "G4ScaledSolid")
338 {
339 ptr = ((G4ScaledSolid*)ptr)->GetUnscaledSolid();
340 continue;
341 }
342 if (type != "G4MultiUnion") // G4MultiUnion doesn't have SetCubVolEpsilon()
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::ostream& os) const
356 {
357 os << "-----------------------------------------------------------\n"
358 << " *** Dump for Boolean solid - " << GetName() << " ***\n"
359 << " ===================================================\n"
360 << " Solid type: " << GetEntityType() << "\n"
361 << " Parameters of constituent solids: \n"
362 << "===========================================================\n";
363 fPtrSolidA->StreamInfo(os);
364 fPtrSolidB->StreamInfo(os);
365 os << "===========================================================\n";
366
367 return os;
368 }
369
370 //////////////////////////////////////////////////////////////////////////
371 //
372 // Creates list of constituent primitives of and their placements
373
374 void G4BooleanSolid::GetListOfPrimitives(
375 std::vector<std::pair<G4VSolid*,G4Transform3D>>& primitives,
376 const G4Transform3D& curPlacement) const
377 {
378 G4Transform3D transform;
379 G4VSolid* solid;
380 G4String type;
381
382 // Repeat two times, first time for fPtrSolidA and then for fPtrSolidB
383 //
384 for (auto i=0; i<2; ++i)
385 {
386 transform = curPlacement;
387 solid = (i == 0) ? fPtrSolidA : fPtrSolidB;
388 type = solid->GetEntityType();
389
390 // While current solid is a trasformed solid just modify transform
391 //
392 while (type == "G4DisplacedSolid" ||
393 type == "G4ReflectedSolid" ||
394 type == "G4ScaledSolid")
395 {
396 if (type == "G4DisplacedSolid")
397 {
398 transform = transform * G4Transform3D(
399 ((G4DisplacedSolid*)solid)->GetObjectRotation(),
400 ((G4DisplacedSolid*)solid)->GetObjectTranslation());
401 solid = ((G4DisplacedSolid*)solid)->GetConstituentMovedSolid();
402 }
403 else if (type == "G4ReflectedSolid")
404 {
405 transform= transform*((G4ReflectedSolid*)solid)->GetDirectTransform3D();
406 solid = ((G4ReflectedSolid*)solid)->GetConstituentMovedSolid();
407 }
408 else if (type == "G4ScaledSolid")
409 {
410 transform = transform * ((G4ScaledSolid*)solid)->GetScaleTransform();
411 solid = ((G4ScaledSolid*)solid)->GetUnscaledSolid();
412 }
413 type = solid->GetEntityType();
414 }
415
416 // If current solid is a Boolean solid then continue recursion,
417 // otherwise add it to the list of primitives
418 //
419 if (type == "G4UnionSolid" ||
420 type == "G4SubtractionSolid" ||
421 type == "G4IntersectionSolid" ||
422 type == "G4BooleanSolid")
423 {
424 ((G4BooleanSolid *)solid)->GetListOfPrimitives(primitives,transform);
425 }
426 else
427 {
428 primitives.emplace_back(solid,transform);
429 }
430 }
431 }
432
433 //////////////////////////////////////////////////////////////////////////
434 //
435 // Returns a point (G4ThreeVector) randomly and uniformly selected
436 // on the surface of the solid
437
438 G4ThreeVector G4BooleanSolid::GetPointOnSurface() const
439 {
440 std::size_t nprims = fPrimitives.size();
441 std::pair<G4VSolid *, G4Transform3D> prim;
442
443 // Get list of primitives and find the total area of their surfaces
444 //
445 if (nprims == 0)
446 {
447 GetListOfPrimitives(fPrimitives, G4Transform3D());
448 nprims = fPrimitives.size();
449 fPrimitivesSurfaceArea = 0.;
450 for (std::size_t i=0; i<nprims; ++i)
451 {
452 fPrimitivesSurfaceArea += fPrimitives[i].first->GetSurfaceArea();
453 }
454 }
455
456 // Select random primitive, get random point on its surface and
457 // check that the point belongs to the surface of the solid
458 //
459 G4ThreeVector p;
460 for (std::size_t k=0; k<100000; ++k) // try 100k times
461 {
462 G4double rand = fPrimitivesSurfaceArea * G4QuickRand();
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 point on the surface have failed!\n"
477 << "The solid created may be an invalid Boolean construct!";
478 G4Exception("G4BooleanSolid::GetPointOnSurface()",
479 "GeomSolids1001", JustWarning, message);
480 return p;
481 }
482
483 //////////////////////////////////////////////////////////////////////////
484 //
485 // Return total number of constituents used for construction of the solid
486
487 G4int G4BooleanSolid::GetNumOfConstituents() const
488 {
489 return (fPtrSolidA->GetNumOfConstituents() + fPtrSolidB->GetNumOfConstituents());
490 }
491
492 //////////////////////////////////////////////////////////////////////////
493 //
494 // Return true if the resulting solid has only planar faces
495
496 G4bool G4BooleanSolid::IsFaceted() const
497 {
498 return (fPtrSolidA->IsFaceted() && fPtrSolidB->IsFaceted());
499 }
500
501 //////////////////////////////////////////////////////////////////////////
502 //
503 // Returns polyhedron for visualization
504
505 G4Polyhedron* G4BooleanSolid::GetPolyhedron () const
506 {
507 if (fpPolyhedron == nullptr ||
508 fRebuildPolyhedron ||
509 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
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 top polyhedron.
524
525 G4Polyhedron*
526 G4BooleanSolid::StackPolyhedron(HepPolyhedronProcessor& processor,
527 const G4VSolid* solid) const
528 {
529 HepPolyhedronProcessor::Operation operation;
530 const G4String& type = solid->GetEntityType();
531 if (type == "G4UnionSolid")
532 { operation = HepPolyhedronProcessor::UNION; }
533 else if (type == "G4IntersectionSolid")
534 { operation = HepPolyhedronProcessor::INTERSECTION; }
535 else if (type == "G4SubtractionSolid")
536 { operation = HepPolyhedronProcessor::SUBTRACTION; }
537 else
538 {
539 std::ostringstream message;
540 message << "Solid - " << solid->GetName()
541 << " - Unrecognised composite solid" << G4endl
542 << " Returning NULL !";
543 G4Exception("StackPolyhedron()", "GeomSolids1001", JustWarning, message);
544 return nullptr;
545 }
546
547 G4Polyhedron* top = nullptr;
548 const G4VSolid* solidA = solid->GetConstituentSolid(0);
549 const G4VSolid* solidB = solid->GetConstituentSolid(1);
550
551 if (solidA->GetConstituentSolid(0) != nullptr)
552 {
553 top = StackPolyhedron(processor, solidA);
554 }
555 else
556 {
557 top = solidA->GetPolyhedron();
558 }
559 G4Polyhedron* operand = solidB->GetPolyhedron();
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 component";
569 G4Exception("G4BooleanSolid::StackPolyhedron()",
570 "GeomSolids2001", JustWarning, message);
571 }
572
573 return top;
574 }
575
576
577 //////////////////////////////////////////////////////////////////////////
578 //
579 // Estimate Cubic Volume (capacity) and cache it for reuse.
580
581 G4double G4BooleanSolid::GetCubicVolume()
582 {
583 if(fCubicVolume < 0.)
584 {
585 fCubicVolume = EstimateCubicVolume(fCubVolStatistics, fCubVolEpsilon);
586 }
587 return fCubicVolume;
588 }
589
590 //////////////////////////////////////////////////////////////////////////
591 //
592 // Set external Boolean processor.
593
594 void
595 G4BooleanSolid::SetExternalBooleanProcessor(G4VBooleanProcessor* extProcessor)
596 {
597 fExternalBoolProcessor = extProcessor;
598 }
599
600 //////////////////////////////////////////////////////////////////////////
601 //
602 // Get external Boolean processor.
603
604 G4VBooleanProcessor* G4BooleanSolid::GetExternalBooleanProcessor()
605 {
606 return fExternalBoolProcessor;
607 }
608