| 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 G4UnionSolid
27 //
28 // 23.04.18 E.Tcherniaev: added extended BBox, yearly return in Inside()
29 // 17.03.17 E.Tcherniaev: revision of SurfaceNormal()
30 // 12.09.98 V.Grichine: first implementation
31 // --------------------------------------------------------------------
32
33 #include <sstream>
34
35 #include "G4UnionSolid.hh"
36
37 #include "G4SystemOfUnits.hh"
38 #include "G4VoxelLimits.hh"
39 #include "G4VPVParameterisation.hh"
40 #include "G4GeometryTolerance.hh"
41
42 #include "G4VGraphicsScene.hh"
43 #include "G4Polyhedron.hh"
44 #include "G4PolyhedronArbitrary.hh"
45 #include "HepPolyhedronProcessor.h"
46
47 #include "G4IntersectionSolid.hh"
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 G4UnionSolid:: G4UnionSolid( const G4String& pName,
55 G4VSolid* pSolidA ,
56 G4VSolid* pSolidB )
57 : G4BooleanSolid(pName,pSolidA,pSolidB)
58 {
59 Init();
60 }
61
62 //////////////////////////////////////////////////////////////////////////
63 //
64 // Constructor
65
66 G4UnionSolid::G4UnionSolid( const G4String& pName,
67 G4VSolid* pSolidA ,
68 G4VSolid* pSolidB ,
69 G4RotationMatrix* rotMatrix,
70 const G4ThreeVector& transVector )
71 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
72
73 {
74 Init();
75 }
76
77 //////////////////////////////////////////////////////////////////////////
78 //
79 // Constructor
80
81 G4UnionSolid::G4UnionSolid( const G4String& pName,
82 G4VSolid* pSolidA ,
83 G4VSolid* pSolidB ,
84 const G4Transform3D& transform )
85 : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
86 {
87 Init();
88 }
89
90 //////////////////////////////////////////////////////////////////////////
91 //
92 // Fake default constructor - sets only member data and allocates memory
93 // for usage restricted to object persistency.
94
95 G4UnionSolid::G4UnionSolid( __void__& a )
96 : G4BooleanSolid(a)
97 {
98 }
99
100 //////////////////////////////////////////////////////////////////////////
101 //
102 // Destructor
103
104 G4UnionSolid::~G4UnionSolid()
105 = default;
106
107 //////////////////////////////////////////////////////////////////////////
108 //
109 // Copy constructor
110
111 G4UnionSolid::G4UnionSolid(const G4UnionSolid& rhs)
112 : G4BooleanSolid (rhs)
113 {
114 fPMin = rhs.fPMin;
115 fPMax = rhs.fPMax;
116 halfCarTolerance=0.5*kCarTolerance;
117 }
118
119 //////////////////////////////////////////////////////////////////////////
120 //
121 // Assignment operator
122
123 G4UnionSolid& G4UnionSolid::operator = (const G4UnionSolid& rhs)
124 {
125 // Check assignment to self
126 //
127 if (this == &rhs) { return *this; }
128
129 // Copy base class data
130 //
131 G4BooleanSolid::operator=(rhs);
132
133 fPMin = rhs.fPMin;
134 fPMax = rhs.fPMax;
135 halfCarTolerance = rhs.halfCarTolerance;
136
137 return *this;
138 }
139
140 //////////////////////////////////////////////////////////////////////////
141 //
142 // Initialisation
143
144 void G4UnionSolid::Init()
145 {
146 G4ThreeVector pdelta(kCarTolerance,kCarTolerance,kCarTolerance);
147 G4ThreeVector pmin, pmax;
148 BoundingLimits(pmin, pmax);
149 fPMin = pmin - pdelta;
150 fPMax = pmax + pdelta;
151 halfCarTolerance=0.5*kCarTolerance;
152 }
153
154 //////////////////////////////////////////////////////////////////////////
155 //
156 // Get bounding box
157
158 void G4UnionSolid::BoundingLimits(G4ThreeVector& pMin,
159 G4ThreeVector& pMax) const
160 {
161 G4ThreeVector minA,maxA, minB,maxB;
162 fPtrSolidA->BoundingLimits(minA,maxA);
163 fPtrSolidB->BoundingLimits(minB,maxB);
164
165 pMin.set(std::min(minA.x(),minB.x()),
166 std::min(minA.y(),minB.y()),
167 std::min(minA.z(),minB.z()));
168
169 pMax.set(std::max(maxA.x(),maxB.x()),
170 std::max(maxA.y(),maxB.y()),
171 std::max(maxA.z(),maxB.z()));
172
173 // Check correctness of the bounding box
174 //
175 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
176 {
177 std::ostringstream message;
178 message << "Bad bounding box (min >= max) for solid: "
179 << GetName() << " !"
180 << "\npMin = " << pMin
181 << "\npMax = " << pMax;
182 G4Exception("G4UnionSolid::BoundingLimits()", "GeomMgt0001",
183 JustWarning, message);
184 DumpInfo();
185 }
186 }
187
188 //////////////////////////////////////////////////////////////////////////
189 //
190 // Calculate extent under transform and specified limit
191
192 G4bool
193 G4UnionSolid::CalculateExtent( const EAxis pAxis,
194 const G4VoxelLimits& pVoxelLimit,
195 const G4AffineTransform& pTransform,
196 G4double& pMin,
197 G4double& pMax ) const
198 {
199 G4bool touchesA, touchesB, out ;
200 G4double minA = kInfinity, minB = kInfinity,
201 maxA = -kInfinity, maxB = -kInfinity;
202
203 touchesA = fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
204 pTransform, minA, maxA);
205 touchesB = fPtrSolidB->CalculateExtent( pAxis, pVoxelLimit,
206 pTransform, minB, maxB);
207 if( touchesA || touchesB )
208 {
209 pMin = std::min( minA, minB );
210 pMax = std::max( maxA, maxB );
211 out = true ;
212 }
213 else
214 {
215 out = false ;
216 }
217
218 return out ; // It exists in this slice if either one does.
219 }
220
221 //////////////////////////////////////////////////////////////////////////
222 //
223 // Important comment: When solids A and B touch together along flat
224 // surface the surface points will be considered as kSurface, while points
225 // located around will correspond to kInside
226
227 EInside G4UnionSolid::Inside( const G4ThreeVector& p ) const
228 {
229 if (std::max(p.z()-fPMax.z(), fPMin.z()-p.z()) > 0) { return kOutside; }
230
231 EInside positionA = fPtrSolidA->Inside(p);
232 if (positionA == kInside) { return positionA; } // inside A
233 EInside positionB = fPtrSolidB->Inside(p);
234 if (positionA == kOutside) { return positionB; }
235
236 if (positionB == kInside) { return positionB; } // inside B
237 if (positionB == kOutside) { return positionA; } // surface A
238
239 // Both points are on surface
240 //
241 static const G4double rtol
242 = 1000*G4GeometryTolerance::GetInstance()->GetRadialTolerance();
243
244 return ((fPtrSolidA->SurfaceNormal(p) +
245 fPtrSolidB->SurfaceNormal(p)).mag2() < rtol) ? kInside : kSurface;
246 }
247
248 //////////////////////////////////////////////////////////////////////////
249 //
250 // Get surface normal
251
252 G4ThreeVector
253 G4UnionSolid::SurfaceNormal( const G4ThreeVector& p ) const
254 {
255 EInside positionA = fPtrSolidA->Inside(p);
256 EInside positionB = fPtrSolidB->Inside(p);
257
258 if (positionA == kSurface &&
259 positionB == kOutside) return fPtrSolidA->SurfaceNormal(p);
260
261 if (positionA == kOutside &&
262 positionB == kSurface) return fPtrSolidB->SurfaceNormal(p);
263
264 if (positionA == kSurface &&
265 positionB == kSurface)
266 {
267 if (Inside(p) == kSurface)
268 {
269 G4ThreeVector normalA = fPtrSolidA->SurfaceNormal(p);
270 G4ThreeVector normalB = fPtrSolidB->SurfaceNormal(p);
271 return (normalA + normalB).unit();
272 }
273 }
274 #ifdef G4BOOLDEBUG
275 G4String surf[3] = { "OUTSIDE", "SURFACE", "INSIDE" };
276 std::ostringstream message;
277 G4int oldprc = message.precision(16);
278 message << "Invalid call of SurfaceNormal(p) for union solid: "
279 << GetName() << " !"
280 << "\nPoint p" << p << " is " << surf[Inside(p)] << " !!!";
281 message.precision(oldprc);
282 G4Exception("G4UnionSolid::SurfaceNormal()", "GeomMgt0001",
283 JustWarning, message);
284 #endif
285 return fPtrSolidA->SurfaceNormal(p);
286 }
287
288 //////////////////////////////////////////////////////////////////////////
289 //
290 // The same algorithm as in DistanceToIn(p)
291
292 G4double
293 G4UnionSolid::DistanceToIn( const G4ThreeVector& p,
294 const G4ThreeVector& v ) const
295 {
296 #ifdef G4BOOLDEBUG
297 if( Inside(p) == kInside )
298 {
299 G4cout << "WARNING - Invalid call in "
300 << "G4UnionSolid::DistanceToIn(p,v)" << G4endl
301 << " Point p is inside !" << G4endl;
302 G4cout << " p = " << p << G4endl;
303 G4cout << " v = " << v << G4endl;
304 G4cerr << "WARNING - Invalid call in "
305 << "G4UnionSolid::DistanceToIn(p,v)" << G4endl
306 << " Point p is inside !" << G4endl;
307 G4cerr << " p = " << p << G4endl;
308 G4cerr << " v = " << v << G4endl;
309 }
310 #endif
311
312 return std::min(fPtrSolidA->DistanceToIn(p,v),
313 fPtrSolidB->DistanceToIn(p,v) ) ;
314 }
315
316 //////////////////////////////////////////////////////////////////////////
317 //
318 // Approximate nearest distance from the point p to the union of
319 // two solids
320
321 G4double
322 G4UnionSolid::DistanceToIn( const G4ThreeVector& p ) const
323 {
324 #ifdef G4BOOLDEBUG
325 if( Inside(p) == kInside )
326 {
327 G4cout << "WARNING - Invalid call in "
328 << "G4UnionSolid::DistanceToIn(p)" << G4endl
329 << " Point p is inside !" << G4endl;
330 G4cout << " p = " << p << G4endl;
331 G4cerr << "WARNING - Invalid call in "
332 << "G4UnionSolid::DistanceToIn(p)" << G4endl
333 << " Point p is inside !" << G4endl;
334 G4cerr << " p = " << p << G4endl;
335 }
336 #endif
337 G4double distA = fPtrSolidA->DistanceToIn(p) ;
338 G4double distB = fPtrSolidB->DistanceToIn(p) ;
339 G4double safety = std::min(distA,distB) ;
340 if(safety < 0.0) safety = 0.0 ;
341 return safety ;
342 }
343
344 //////////////////////////////////////////////////////////////////////////
345 //
346 // The same algorithm as DistanceToOut(p)
347
348 G4double
349 G4UnionSolid::DistanceToOut( const G4ThreeVector& p,
350 const G4ThreeVector& v,
351 const G4bool calcNorm,
352 G4bool* validNorm,
353 G4ThreeVector* n ) const
354 {
355 G4double dist = 0.0, disTmp = 0.0 ;
356 G4ThreeVector normTmp;
357 G4ThreeVector* nTmp = &normTmp;
358
359 if( Inside(p) == kOutside )
360 {
361 #ifdef G4BOOLDEBUG
362 G4cout << "Position:" << G4endl << G4endl;
363 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
364 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
365 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
366 G4cout << "Direction:" << G4endl << G4endl;
367 G4cout << "v.x() = " << v.x() << G4endl;
368 G4cout << "v.y() = " << v.y() << G4endl;
369 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
370 G4cout << "WARNING - Invalid call in "
371 << "G4UnionSolid::DistanceToOut(p,v)" << G4endl
372 << " Point p is outside !" << G4endl;
373 G4cout << " p = " << p << G4endl;
374 G4cout << " v = " << v << G4endl;
375 G4cerr << "WARNING - Invalid call in "
376 << "G4UnionSolid::DistanceToOut(p,v)" << G4endl
377 << " Point p is outside !" << G4endl;
378 G4cerr << " p = " << p << G4endl;
379 G4cerr << " v = " << v << G4endl;
380 #endif
381 }
382 else
383 {
384 EInside positionA = fPtrSolidA->Inside(p) ;
385
386 if( positionA != kOutside )
387 {
388 do // Loop checking, 13.08.2015, G.Cosmo
389 {
390 disTmp = fPtrSolidA->DistanceToOut(p+dist*v,v,calcNorm,
391 validNorm,nTmp);
392 dist += disTmp ;
393
394 if(fPtrSolidB->Inside(p+dist*v) != kOutside)
395 {
396 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v,calcNorm,
397 validNorm,nTmp);
398 dist += disTmp ;
399 }
400 }
401 while( (fPtrSolidA->Inside(p+dist*v) != kOutside)
402 && (disTmp > halfCarTolerance) );
403 }
404 else // if( positionB != kOutside )
405 {
406 do // Loop checking, 13.08.2015, G.Cosmo
407 {
408 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v,calcNorm,
409 validNorm,nTmp);
410 dist += disTmp ;
411
412 if(fPtrSolidA->Inside(p+dist*v) != kOutside)
413 {
414 disTmp = fPtrSolidA->DistanceToOut(p+dist*v,v,calcNorm,
415 validNorm,nTmp);
416 dist += disTmp ;
417 }
418 }
419 while( (fPtrSolidB->Inside(p+dist*v) != kOutside)
420 && (disTmp > halfCarTolerance) );
421 }
422 }
423 if( calcNorm )
424 {
425 *validNorm = false ;
426 *n = *nTmp ;
427 }
428 return dist ;
429 }
430
431 //////////////////////////////////////////////////////////////////////////
432 //
433 // Inverted algorithm of DistanceToIn(p)
434
435 G4double
436 G4UnionSolid::DistanceToOut( const G4ThreeVector& p ) const
437 {
438 G4double distout = 0.0;
439 if( Inside(p) == kOutside )
440 {
441 #ifdef G4BOOLDEBUG
442 G4cout << "WARNING - Invalid call in "
443 << "G4UnionSolid::DistanceToOut(p)" << G4endl
444 << " Point p is outside !" << G4endl;
445 G4cout << " p = " << p << G4endl;
446 G4cerr << "WARNING - Invalid call in "
447 << "G4UnionSolid::DistanceToOut(p)" << G4endl
448 << " Point p is outside !" << G4endl;
449 G4cerr << " p = " << p << G4endl;
450 #endif
451 }
452 else
453 {
454 EInside positionA = fPtrSolidA->Inside(p) ;
455 EInside positionB = fPtrSolidB->Inside(p) ;
456
457 // Is this equivalent ??
458 // if( ! ( (positionA == kOutside)) &&
459 // (positionB == kOutside)) )
460 if((positionA == kInside && positionB == kInside ) ||
461 (positionA == kInside && positionB == kSurface ) ||
462 (positionA == kSurface && positionB == kInside ) )
463 {
464 distout= std::max(fPtrSolidA->DistanceToOut(p),
465 fPtrSolidB->DistanceToOut(p) ) ;
466 }
467 else
468 {
469 if(positionA == kOutside)
470 {
471 distout= fPtrSolidB->DistanceToOut(p) ;
472 }
473 else
474 {
475 distout= fPtrSolidA->DistanceToOut(p) ;
476 }
477 }
478 }
479 return distout;
480 }
481
482 //////////////////////////////////////////////////////////////////////////
483 //
484 // GetEntityType
485
486 G4GeometryType G4UnionSolid::GetEntityType() const
487 {
488 return {"G4UnionSolid"};
489 }
490
491 //////////////////////////////////////////////////////////////////////////
492 //
493 // Make a clone of the object
494
495 G4VSolid* G4UnionSolid::Clone() const
496 {
497 return new G4UnionSolid(*this);
498 }
499
500 //////////////////////////////////////////////////////////////////////////
501 //
502 // ComputeDimensions
503
504 void
505 G4UnionSolid::ComputeDimensions( G4VPVParameterisation*,
506 const G4int,
507 const G4VPhysicalVolume* )
508 {
509 }
510
511 //////////////////////////////////////////////////////////////////////////
512 //
513 // DescribeYourselfTo
514
515 void
516 G4UnionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
517 {
518 scene.AddSolid (*this);
519 }
520
521 //////////////////////////////////////////////////////////////////////////
522 //
523 // CreatePolyhedron
524
525 G4Polyhedron*
526 G4UnionSolid::CreatePolyhedron () const
527 {
528 if (fExternalBoolProcessor == nullptr)
529 {
530 HepPolyhedronProcessor processor;
531 // Stack components and components of components recursively
532 // See G4BooleanSolid::StackPolyhedron
533 G4Polyhedron* top = StackPolyhedron(processor, this);
534 auto result = new G4Polyhedron(*top);
535 if (processor.execute(*result))
536 {
537 return result;
538 }
539 else
540 {
541 return nullptr;
542 }
543 }
544 else
545 {
546 return fExternalBoolProcessor->Process(this);
547 }
548 }
549
550 //////////////////////////////////////////////////////////////////////////
551 //
552 // GetCubicVolume
553
554 G4double G4UnionSolid::GetCubicVolume()
555 {
556 if( fCubicVolume >= 0. )
557 {
558 return fCubicVolume;
559 }
560 G4ThreeVector bminA, bmaxA, bminB, bmaxB;
561 fPtrSolidA->BoundingLimits(bminA, bmaxA);
562 fPtrSolidB->BoundingLimits(bminB, bmaxB);
563 G4bool noIntersection =
564 bminA.x() >= bmaxB.x() || bminA.y() >= bmaxB.y() || bminA.z() >= bmaxB.z() ||
565 bminB.x() >= bmaxA.x() || bminB.y() >= bmaxA.y() || bminB.z() >= bmaxA.z();
566
567 if (noIntersection)
568 {
569 fCubicVolume = fPtrSolidA->GetCubicVolume() + fPtrSolidB->GetCubicVolume();
570 }
571 else
572 {
573 if (GetNumOfConstituents() > 10)
574 {
575 fCubicVolume = G4BooleanSolid::GetCubicVolume();
576 }
577 else
578 {
579 G4IntersectionSolid intersectVol("Temporary-Intersection-for-Union",
580 fPtrSolidA, fPtrSolidB);
581 intersectVol.SetCubVolStatistics(GetCubVolStatistics());
582 intersectVol.SetCubVolEpsilon(GetCubVolEpsilon());
583
584 fCubicVolume = fPtrSolidA->GetCubicVolume() + fPtrSolidB->GetCubicVolume()
585 - intersectVol.GetCubicVolume();
586 }
587 }
588 return fCubicVolume;
589 }
590