| 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 methods for the class G4IntersectionSolid
27 //
28 // 22.07.11 T.Nikitina: added detection of infinite loop in DistanceToIn(p,v)
29 // 19.10.98 V.Grichine: new algorithm of DistanceToIn(p,v)
30 // 14.10.98 V.Grichine: implementation of the first version
31 // --------------------------------------------------------------------
32
33 #include "G4SubtractionSolid.hh"
34
35 #include "G4SystemOfUnits.hh"
36 #include "G4VoxelLimits.hh"
37 #include "G4VPVParameterisation.hh"
38 #include "G4GeometryTolerance.hh"
39
40 #include "G4VGraphicsScene.hh"
41 #include "G4Polyhedron.hh"
42 #include "G4PolyhedronArbitrary.hh"
43 #include "HepPolyhedronProcessor.h"
44
45 #include "G4IntersectionSolid.hh"
46
47 #include <sstream>
48
49 //////////////////////////////////////////////////////////////////////////
50 //
51 // Transfer all data members to G4BooleanSolid which is responsible
52 // for them. pName will be in turn sent to G4VSolid
53
54 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
55 G4VSolid* pSolidA ,
56 G4VSolid* pSolidB )
57 : G4BooleanSolid(pName,pSolidA,pSolidB)
58 {
59 }
60
61 //////////////////////////////////////////////////////////////////////////
62 //
63 // Constructor
64
65 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
66 G4VSolid* pSolidA ,
67 G4VSolid* pSolidB ,
68 G4RotationMatrix* rotMatrix,
69 const G4ThreeVector& transVector )
70 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
71 {
72 }
73
74 //////////////////////////////////////////////////////////////////////////
75 //
76 // Constructor
77
78 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
79 G4VSolid* pSolidA ,
80 G4VSolid* pSolidB ,
81 const G4Transform3D& transform )
82 : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
83 {
84 }
85
86 //////////////////////////////////////////////////////////////////////////
87 //
88 // Fake default constructor - sets only member data and allocates memory
89 // for usage restricted to object persistency.
90
91 G4SubtractionSolid::G4SubtractionSolid( __void__& a )
92 : G4BooleanSolid(a)
93 {
94 }
95
96 //////////////////////////////////////////////////////////////////////////
97 //
98 // Destructor
99
100 G4SubtractionSolid::~G4SubtractionSolid() = default;
101
102 //////////////////////////////////////////////////////////////////////////
103 //
104 // Copy constructor
105
106 G4SubtractionSolid::G4SubtractionSolid(const G4SubtractionSolid&) = default;
107
108 //////////////////////////////////////////////////////////////////////////
109 //
110 // Assignment operator
111
112 G4SubtractionSolid&
113 G4SubtractionSolid::operator = (const G4SubtractionSolid& rhs)
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(G4ThreeVector& pMin,
132 G4ThreeVector& pMax) const
133 {
134 // Since it is unclear how the shape of the first solid will be changed
135 // after subtraction, just return its original bounding box.
136 //
137 fPtrSolidA->BoundingLimits(pMin,pMax);
138
139 // Check correctness of the bounding box
140 //
141 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
142 {
143 std::ostringstream message;
144 message << "Bad bounding box (min >= max) for solid: "
145 << GetName() << " !"
146 << "\npMin = " << pMin
147 << "\npMax = " << pMax;
148 G4Exception("G4SubtractionSolid::BoundingLimits()", "GeomMgt0001",
149 JustWarning, message);
150 DumpInfo();
151 }
152 }
153
154 //////////////////////////////////////////////////////////////////////////
155 //
156 // Calculate extent under transform and specified limit
157
158 G4bool
159 G4SubtractionSolid::CalculateExtent( const EAxis pAxis,
160 const G4VoxelLimits& pVoxelLimit,
161 const G4AffineTransform& pTransform,
162 G4double& pMin,
163 G4double& pMax ) const
164 {
165 // Since we cannot be sure how much the second solid subtracts
166 // from the first, we must use the first solid's extent!
167
168 return fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
169 pTransform, pMin, pMax );
170 }
171
172 //////////////////////////////////////////////////////////////////////////
173 //
174 // Touching ? Empty subtraction ?
175
176 EInside G4SubtractionSolid::Inside( const G4ThreeVector& p ) const
177 {
178 EInside positionA = fPtrSolidA->Inside(p);
179 if (positionA == kOutside) return positionA; // outside A
180
181 EInside positionB = fPtrSolidB->Inside(p);
182 if (positionB == kOutside) return positionA;
183
184 if (positionB == kInside) return kOutside;
185 if (positionA == kInside) return kSurface; // surface B
186
187 // Point is on both surfaces
188 //
189 static const G4double rtol = 1000*kCarTolerance;
190
191 return ((fPtrSolidA->SurfaceNormal(p) -
192 fPtrSolidB->SurfaceNormal(p)).mag2() > rtol) ? kSurface : kOutside;
193 }
194
195 //////////////////////////////////////////////////////////////////////////
196 //
197 // SurfaceNormal
198
199 G4ThreeVector
200 G4SubtractionSolid::SurfaceNormal( const G4ThreeVector& p ) const
201 {
202 G4ThreeVector normal;
203
204 EInside InsideA = fPtrSolidA->Inside(p);
205 EInside InsideB = fPtrSolidB->Inside(p);
206
207 if( InsideA == kOutside )
208 {
209 #ifdef G4BOOLDEBUG
210 G4cout << "WARNING - Invalid call [1] in "
211 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
212 << " Point p is outside !" << G4endl;
213 G4cout << " p = " << p << G4endl;
214 G4cerr << "WARNING - Invalid call [1] in "
215 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
216 << " Point p is outside !" << G4endl;
217 G4cerr << " p = " << p << G4endl;
218 #endif
219 normal = fPtrSolidA->SurfaceNormal(p) ;
220 }
221 else if( InsideA == kSurface &&
222 InsideB != kInside )
223 {
224 normal = fPtrSolidA->SurfaceNormal(p) ;
225 }
226 else if( InsideA == kInside &&
227 InsideB != kOutside )
228 {
229 normal = -fPtrSolidB->SurfaceNormal(p) ;
230 }
231 else
232 {
233 if ( fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToIn(p) )
234 {
235 normal = fPtrSolidA->SurfaceNormal(p) ;
236 }
237 else
238 {
239 normal = -fPtrSolidB->SurfaceNormal(p) ;
240 }
241 #ifdef G4BOOLDEBUG
242 if(Inside(p) == kInside)
243 {
244 G4cout << "WARNING - Invalid call [2] in "
245 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
246 << " Point p is inside !" << G4endl;
247 G4cout << " p = " << p << G4endl;
248 G4cerr << "WARNING - Invalid call [2] in "
249 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
250 << " Point p is inside !" << G4endl;
251 G4cerr << " p = " << p << G4endl;
252 }
253 #endif
254 }
255 return normal;
256 }
257
258 //////////////////////////////////////////////////////////////////////////
259 //
260 // The same algorithm as in DistanceToIn(p)
261
262 G4double
263 G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p,
264 const G4ThreeVector& v ) const
265 {
266 G4double dist = 0.0, dist2 = 0.0, disTmp = 0.0;
267
268 #ifdef G4BOOLDEBUG
269 if( Inside(p) == kInside )
270 {
271 G4cout << "WARNING - Invalid call in "
272 << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
273 << " Point p is inside !" << G4endl;
274 G4cout << " p = " << p << G4endl;
275 G4cout << " v = " << v << G4endl;
276 G4cerr << "WARNING - Invalid call in "
277 << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
278 << " Point p is inside !" << G4endl;
279 G4cerr << " p = " << p << G4endl;
280 G4cerr << " v = " << v << G4endl;
281 }
282 #endif
283
284 // if( // ( fPtrSolidA->Inside(p) != kOutside) && // case1:p in both A&B
285 if ( fPtrSolidB->Inside(p) != kOutside ) // start: out of B
286 {
287 dist = fPtrSolidB->DistanceToOut(p,v) ; // ,calcNorm,validNorm,n) ;
288
289 if( fPtrSolidA->Inside(p+dist*v) != kInside )
290 {
291 G4int count1=0;
292 do // Loop checking, 13.08.2015, G.Cosmo
293 {
294 disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
295
296 if(disTmp == kInfinity)
297 {
298 return kInfinity ;
299 }
300 dist += disTmp ;
301
302 if( Inside(p+dist*v) == kOutside )
303 {
304 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
305 dist2 = dist+disTmp;
306 if (dist == dist2) { return dist; } // no progress
307 dist = dist2 ;
308 ++count1;
309 if( count1 > 1000 ) // Infinite loop detected
310 {
311 G4String nameB = fPtrSolidB->GetName();
312 if(fPtrSolidB->GetEntityType()=="G4DisplacedSolid")
313 {
314 nameB = (dynamic_cast<G4DisplacedSolid*>(fPtrSolidB))
315 ->GetConstituentMovedSolid()->GetName();
316 }
317 std::ostringstream message;
318 message << "Illegal condition caused by solids: "
319 << fPtrSolidA->GetName() << " and " << nameB << G4endl;
320 message.precision(16);
321 message << "Looping detected in point " << p+dist*v
322 << ", from original point " << p
323 << " and direction " << v << G4endl
324 << "Computed candidate distance: " << dist << "*mm. ";
325 message.precision(6);
326 DumpInfo();
327 G4Exception("G4SubtractionSolid::DistanceToIn(p,v)",
328 "GeomSolids1001", JustWarning, message,
329 "Returning candidate distance.");
330 return dist;
331 }
332 }
333 }
334 while( Inside(p+dist*v) == kOutside ) ;
335 }
336 }
337 else // p outside A, start in A
338 {
339 dist = fPtrSolidA->DistanceToIn(p,v) ;
340
341 if( dist == kInfinity ) // past A, hence past A\B
342 {
343 return kInfinity ;
344 }
345 else
346 {
347 G4int count2=0;
348 while( Inside(p+dist*v) == kOutside ) // pushing loop
349 {
350 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
351 dist += disTmp ;
352
353 if( Inside(p+dist*v) == kOutside )
354 {
355 disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
356
357 if(disTmp == kInfinity) // past A, hence past A\B
358 {
359 return kInfinity ;
360 }
361 dist2 = dist+disTmp;
362 if (dist == dist2) { return dist; } // no progress
363 dist = dist2 ;
364 ++count2;
365 if( count2 > 1000 ) // Infinite loop detected
366 {
367 G4String nameB = fPtrSolidB->GetName();
368 if(fPtrSolidB->GetEntityType()=="G4DisplacedSolid")
369 {
370 nameB = (dynamic_cast<G4DisplacedSolid*>(fPtrSolidB))
371 ->GetConstituentMovedSolid()->GetName();
372 }
373 std::ostringstream message;
374 message << "Illegal condition caused by solids: "
375 << fPtrSolidA->GetName() << " and " << nameB << G4endl;
376 message.precision(16);
377 message << "Looping detected in point " << p+dist*v
378 << ", from original point " << p
379 << " and direction " << v << G4endl
380 << "Computed candidate distance: " << dist << "*mm. ";
381 message.precision(6);
382 DumpInfo();
383 G4Exception("G4SubtractionSolid::DistanceToIn(p,v)",
384 "GeomSolids1001", JustWarning, message,
385 "Returning candidate distance.");
386 return dist;
387 }
388 }
389 } // Loop checking, 13.08.2015, G.Cosmo
390 }
391 }
392
393 return dist ;
394 }
395
396 //////////////////////////////////////////////////////////////////////////
397 //
398 // Approximate nearest distance from the point p to the intersection of
399 // two solids. It is usually underestimated from the point of view of
400 // isotropic safety
401
402 G4double
403 G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p ) const
404 {
405 G4double dist = 0.0;
406
407 #ifdef G4BOOLDEBUG
408 if( Inside(p) == kInside )
409 {
410 G4cout << "WARNING - Invalid call in "
411 << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
412 << " Point p is inside !" << G4endl;
413 G4cout << " p = " << p << G4endl;
414 G4cerr << "WARNING - Invalid call in "
415 << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
416 << " Point p is inside !" << G4endl;
417 G4cerr << " p = " << p << G4endl;
418 }
419 #endif
420
421 if( ( fPtrSolidA->Inside(p) != kOutside) && // case 1
422 ( fPtrSolidB->Inside(p) != kOutside) )
423 {
424 dist = fPtrSolidB->DistanceToOut(p);
425 }
426 else
427 {
428 dist = fPtrSolidA->DistanceToIn(p);
429 }
430
431 return dist;
432 }
433
434 //////////////////////////////////////////////////////////////////////////
435 //
436 // The same algorithm as DistanceToOut(p)
437
438 G4double
439 G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p,
440 const G4ThreeVector& v,
441 const G4bool calcNorm,
442 G4bool* validNorm,
443 G4ThreeVector* n ) const
444 {
445 #ifdef G4BOOLDEBUG
446 if( Inside(p) == kOutside )
447 {
448 G4cout << "Position:" << G4endl << G4endl;
449 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
450 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
451 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
452 G4cout << "Direction:" << G4endl << G4endl;
453 G4cout << "v.x() = " << v.x() << G4endl;
454 G4cout << "v.y() = " << v.y() << G4endl;
455 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
456 G4cout << "WARNING - Invalid call in "
457 << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
458 << " Point p is outside !" << G4endl;
459 G4cout << " p = " << p << G4endl;
460 G4cout << " v = " << v << G4endl;
461 G4cerr << "WARNING - Invalid call in "
462 << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
463 << " Point p is outside !" << G4endl;
464 G4cerr << " p = " << p << G4endl;
465 G4cerr << " v = " << v << G4endl;
466 }
467 #endif
468
469 G4double distout;
470 G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,validNorm,n) ;
471 G4double distB = fPtrSolidB->DistanceToIn(p,v) ;
472 if(distB < distA)
473 {
474 if(calcNorm)
475 {
476 *n = -(fPtrSolidB->SurfaceNormal(p+distB*v)) ;
477 *validNorm = false ;
478 }
479 distout= distB ;
480 }
481 else
482 {
483 distout= distA ;
484 }
485 return distout;
486 }
487
488 //////////////////////////////////////////////////////////////////////////
489 //
490 // Inverted algorithm of DistanceToIn(p)
491
492 G4double
493 G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p ) const
494 {
495 G4double dist=0.0;
496
497 if( Inside(p) == kOutside )
498 {
499 #ifdef G4BOOLDEBUG
500 G4cout << "WARNING - Invalid call in "
501 << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
502 << " Point p is outside" << G4endl;
503 G4cout << " p = " << p << G4endl;
504 G4cerr << "WARNING - Invalid call in "
505 << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
506 << " Point p is outside" << G4endl;
507 G4cerr << " p = " << p << G4endl;
508 #endif
509 }
510 else
511 {
512 dist= std::min(fPtrSolidA->DistanceToOut(p),
513 fPtrSolidB->DistanceToIn(p) ) ;
514 }
515 return dist;
516 }
517
518 //////////////////////////////////////////////////////////////////////////
519 //
520 //
521
522 G4GeometryType G4SubtractionSolid::GetEntityType() const
523 {
524 return {"G4SubtractionSolid"};
525 }
526
527 //////////////////////////////////////////////////////////////////////////
528 //
529 // Make a clone of the object
530
531 G4VSolid* G4SubtractionSolid::Clone() const
532 {
533 return new G4SubtractionSolid(*this);
534 }
535
536 //////////////////////////////////////////////////////////////////////////
537 //
538 // ComputeDimensions
539
540 void
541 G4SubtractionSolid::ComputeDimensions( G4VPVParameterisation*,
542 const G4int,
543 const G4VPhysicalVolume* )
544 {
545 }
546
547 //////////////////////////////////////////////////////////////////////////
548 //
549 // DescribeYourselfTo
550
551 void
552 G4SubtractionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
553 {
554 scene.AddSolid (*this);
555 }
556
557 //////////////////////////////////////////////////////////////////////////
558 //
559 // CreatePolyhedron
560
561 G4Polyhedron* G4SubtractionSolid::CreatePolyhedron () const
562 {
563 if (fExternalBoolProcessor == nullptr)
564 {
565 HepPolyhedronProcessor processor;
566 // Stack components and components of components recursively
567 // See G4BooleanSolid::StackPolyhedron
568 G4Polyhedron* top = StackPolyhedron(processor, this);
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(this);
582 }
583 }
584
585 //////////////////////////////////////////////////////////////////////////
586 //
587 // GetCubicVolume
588 //
589
590 G4double G4SubtractionSolid::GetCubicVolume()
591 {
592 if( fCubicVolume >= 0. )
593 {
594 return fCubicVolume;
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() >= bmaxB.y() || bminA.z() >= bmaxB.z() ||
601 bminB.x() >= bmaxA.x() || bminB.y() >= bmaxA.y() || bminB.z() >= bmaxA.z();
602
603 if (noIntersection)
604 {
605 fCubicVolume = fPtrSolidA->GetCubicVolume();
606 }
607 else
608 {
609 if (GetNumOfConstituents() > 10)
610 {
611 fCubicVolume = G4BooleanSolid::GetCubicVolume();
612 }
613 else
614 {
615 G4IntersectionSolid intersectVol("Temporary-Intersection-for-Subtraction",
616 fPtrSolidA, fPtrSolidB);
617 intersectVol.SetCubVolStatistics(GetCubVolStatistics());
618 intersectVol.SetCubVolEpsilon(GetCubVolEpsilon());
619
620 G4double cubVolumeA = fPtrSolidA->GetCubicVolume();
621 fCubicVolume = cubVolumeA - intersectVol.GetCubicVolume();
622 if (fCubicVolume < 0.01*cubVolumeA) fCubicVolume = G4BooleanSolid::GetCubicVolume();
623 }
624 }
625 return fCubicVolume;
626 }
627