Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/Boolean/src/G4MultiUnion.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /geometry/solids/Boolean/src/G4MultiUnion.cc (Version 11.3.0) and /geometry/solids/Boolean/src/G4MultiUnion.cc (Version 10.3.p3)


  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 G4MultiUnion class           
 27 //                                                
 28 // 19.10.12 M.Gayer - Original implementation     
 29 // 06.04.17 G.Cosmo - Adapted implementation i    
 30 // -------------------------------------------    
 31                                                   
 32 #include <iostream>                               
 33 #include <sstream>                                
 34                                                   
 35 #include "G4MultiUnion.hh"                        
 36 #include "Randomize.hh"                           
 37 #include "G4GeometryTolerance.hh"                 
 38 #include "G4BoundingEnvelope.hh"                  
 39 #include "G4AffineTransform.hh"                   
 40 #include "G4DisplacedSolid.hh"                    
 41                                                   
 42 #include "G4VGraphicsScene.hh"                    
 43 #include "G4Polyhedron.hh"                        
 44 #include "G4PolyhedronArbitrary.hh"               
 45 #include "HepPolyhedronProcessor.h"               
 46                                                   
 47 #include "G4BooleanSolid.hh"                      
 48                                                   
 49 #include "G4AutoLock.hh"                          
 50                                                   
 51 namespace                                         
 52 {                                                 
 53   G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE    
 54 }                                                 
 55                                                   
 56 //____________________________________________    
 57 G4MultiUnion::G4MultiUnion(const G4String& nam    
 58   : G4VSolid(name)                                
 59 {                                                 
 60   SetName(name);                                  
 61   fSolids.clear();                                
 62   fTransformObjs.clear();                         
 63   kRadTolerance = G4GeometryTolerance::GetInst    
 64 }                                                 
 65                                                   
 66 //____________________________________________    
 67 G4MultiUnion::~G4MultiUnion()                     
 68 = default;                                        
 69                                                   
 70 //____________________________________________    
 71 void G4MultiUnion::AddNode(G4VSolid& solid, co    
 72 {                                                 
 73   fSolids.push_back(&solid);                      
 74   fTransformObjs.push_back(trans);  // Store a    
 75 }                                                 
 76                                                   
 77 //____________________________________________    
 78 void G4MultiUnion::AddNode(G4VSolid* solid, co    
 79 {                                                 
 80   fSolids.push_back(solid);                       
 81   fTransformObjs.push_back(trans);  // Store a    
 82 }                                                 
 83                                                   
 84 //____________________________________________    
 85 G4VSolid* G4MultiUnion::Clone() const             
 86 {                                                 
 87   return new G4MultiUnion(*this);                 
 88 }                                                 
 89                                                   
 90 // Copy constructor                               
 91 //____________________________________________    
 92 G4MultiUnion::G4MultiUnion(const G4MultiUnion&    
 93   : G4VSolid(rhs), fCubicVolume(rhs.fCubicVolu    
 94     fSurfaceArea(rhs.fSurfaceArea),               
 95     kRadTolerance(rhs.kRadTolerance), fAccurat    
 96 {                                                 
 97 }                                                 
 98                                                   
 99 // Fake default constructor for persistency       
100 //____________________________________________    
101 G4MultiUnion::G4MultiUnion( __void__& a )         
102   : G4VSolid(a)                                   
103 {                                                 
104 }                                                 
105                                                   
106 // Assignment operator                            
107 //____________________________________________    
108 G4MultiUnion& G4MultiUnion::operator = (const     
109 {                                                 
110   // Check assignment to self                     
111   //                                              
112   if (this == &rhs)                               
113   {                                               
114     return *this;                                 
115   }                                               
116                                                   
117   // Copy base class data                         
118   //                                              
119   G4VSolid::operator=(rhs);                       
120                                                   
121   return *this;                                   
122 }                                                 
123                                                   
124 //____________________________________________    
125 G4double G4MultiUnion::GetCubicVolume()           
126 {                                                 
127   if (fCubicVolume == 0.0)                        
128   {                                               
129     fCubicVolume = EstimateCubicVolume(1000000    
130   }                                               
131   return fCubicVolume;                            
132 }                                                 
133                                                   
134 //____________________________________________    
135 G4double                                          
136 G4MultiUnion::DistanceToInNoVoxels(const G4Thr    
137                                    const G4Thr    
138 {                                                 
139   G4ThreeVector direction = aDirection.unit();    
140   G4ThreeVector localPoint, localDirection;       
141   G4double minDistance = kInfinity;               
142                                                   
143   std::size_t numNodes = fSolids.size();          
144   for (std::size_t i = 0 ; i < numNodes ; ++i)    
145   {                                               
146     G4VSolid& solid = *fSolids[i];                
147     const G4Transform3D& transform = fTransfor    
148                                                   
149     localPoint = GetLocalPoint(transform, aPoi    
150     localDirection = GetLocalVector(transform,    
151                                                   
152     G4double distance = solid.DistanceToIn(loc    
153     if (minDistance > distance) minDistance =     
154   }                                               
155   return minDistance;                             
156 }                                                 
157                                                   
158 //____________________________________________    
159 G4double G4MultiUnion::DistanceToInCandidates(    
160                                                   
161                                                   
162                                                   
163 {                                                 
164   std::size_t candidatesCount = candidates.siz    
165   G4ThreeVector localPoint, localDirection;       
166                                                   
167   G4double minDistance = kInfinity;               
168   for (std::size_t i = 0 ; i < candidatesCount    
169   {                                               
170     G4int candidate = candidates[i];              
171     G4VSolid& solid = *fSolids[candidate];        
172     const G4Transform3D& transform = fTransfor    
173                                                   
174     localPoint = GetLocalPoint(transform, aPoi    
175     localDirection = GetLocalVector(transform,    
176     G4double distance = solid.DistanceToIn(loc    
177     if (minDistance > distance) minDistance =     
178     bits.SetBitNumber(candidate);                 
179     if (minDistance == 0) break;                  
180   }                                               
181   return minDistance;                             
182 }                                                 
183                                                   
184 // Algorithm note: we have to look also for al    
185 // if the distance is not shorter ... we have     
186 // for example for objects which starts in fir    
187 // do not collide with direction line, but in     
188 // The idea of crossing voxels would be still     
189 // because this way we could exclude from the     
190 // which were found that obviously are not goo    
191 // they would return infinity                     
192 // But if distance is smaller than the shift t    
193 // it immediately                                 
194 //____________________________________________    
195 G4double G4MultiUnion::DistanceToIn(const G4Th    
196                                     const G4Th    
197 {                                                 
198   G4double minDistance = kInfinity;               
199   G4ThreeVector direction = aDirection.unit();    
200   G4double shift = fVoxels.DistanceToFirst(aPo    
201   if (shift == kInfinity) return shift;           
202                                                   
203   G4ThreeVector currentPoint = aPoint;            
204   if (shift != 0.0) currentPoint += direction     
205                                                   
206   G4SurfBits exclusion(fVoxels.GetBitsPerSlice    
207   std::vector<G4int> candidates, curVoxel(3);     
208   fVoxels.GetVoxel(curVoxel, currentPoint);       
209                                                   
210   do                                              
211   {                                               
212     {                                             
213       if (fVoxels.GetCandidatesVoxelArray(curV    
214       {                                           
215         G4double distance = DistanceToInCandid    
216                                                   
217         if (minDistance > distance) minDistanc    
218         if (distance < shift) break;              
219       }                                           
220     }                                             
221     shift = fVoxels.DistanceToNext(aPoint, dir    
222   }                                               
223   while (minDistance > shift);                    
224                                                   
225   return minDistance;                             
226 }                                                 
227                                                   
228 //____________________________________________    
229 G4double G4MultiUnion::DistanceToOutNoVoxels(c    
230                                              c    
231                                              G    
232 {                                                 
233   // Computes distance from a point presumably    
234   // surface. Ignores first surface if the poi    
235   // Early return infinity in case the safety     
236   // than the proposed step aPstep.               
237   // The normal vector to the crossed surface     
238   // is crossed, otherwise aNormal->IsNull() i    
239                                                   
240   // algorithm:                                   
241   G4ThreeVector direction = aDirection.unit();    
242   G4ThreeVector localPoint, localDirection;       
243   G4int ignoredSolid = -1;                        
244   G4double resultDistToOut = 0;                   
245   G4ThreeVector currentPoint = aPoint;            
246                                                   
247   auto numNodes = (G4int)fSolids.size();          
248   for (auto i = 0; i < numNodes; ++i)             
249   {                                               
250     if (i != ignoredSolid)                        
251     {                                             
252       G4VSolid& solid = *fSolids[i];              
253       const G4Transform3D& transform = fTransf    
254       localPoint = GetLocalPoint(transform, cu    
255       localDirection = GetLocalVector(transfor    
256       EInside location = solid.Inside(localPoi    
257       if (location != EInside::kOutside)          
258       {                                           
259         G4double distance = solid.DistanceToOu    
260                                                   
261         if (distance < kInfinity)                 
262         {                                         
263           if (resultDistToOut == kInfinity) re    
264           if (distance > 0)                       
265           {                                       
266             currentPoint = GetGlobalPoint(tran    
267                                           + di    
268             resultDistToOut += distance;          
269             ignoredSolid = i; // skip the soli    
270             i = -1; // force the loop to conti    
271           }                                       
272         }                                         
273       }                                           
274     }                                             
275   }                                               
276   return resultDistToOut;                         
277 }                                                 
278                                                   
279 //____________________________________________    
280 G4double G4MultiUnion::DistanceToOut(const G4T    
281                                      const G4T    
282                                      const G4b    
283                                      G4bool* /    
284                                      G4ThreeVe    
285 {                                                 
286   return DistanceToOutVoxels(aPoint, aDirectio    
287 }                                                 
288                                                   
289 //____________________________________________    
290 G4double G4MultiUnion::DistanceToOutVoxels(con    
291                                            con    
292                                            G4T    
293 {                                                 
294   // Computes distance from a point presumably    
295   // surface. Ignores first surface along each    
296   // inside or outside. Early returns zero in     
297   // the starting point.                          
298   // o The proposed step is ignored.              
299   // o The normal vector to the crossed surfac    
300                                                   
301   // In the case the considered point is locat    
302   // structure, the treatments are as follows:    
303   //      - investigation of the candidates fo    
304   //      - progressive moving of the point to    
305   //        passed direction                      
306   //      - processing of the normal              
307                                                   
308   G4ThreeVector direction = aDirection.unit();    
309   std::vector<G4int> candidates;                  
310   G4double distance = 0;                          
311   std::size_t numNodes = 2*fSolids.size();        
312   std::size_t count=0;                            
313                                                   
314   if (fVoxels.GetCandidatesVoxelArray(aPoint,     
315   {                                               
316     // For normal case for which we presume th    
317     G4ThreeVector localPoint, localDirection,     
318     G4ThreeVector currentPoint = aPoint;          
319     G4SurfBits exclusion(fVoxels.GetBitsPerSli    
320     G4bool notOutside;                            
321     G4ThreeVector maxNormal;                      
322                                                   
323     do                                            
324     {                                             
325       notOutside = false;                         
326                                                   
327       G4double maxDistance = -kInfinity;          
328       G4int maxCandidate = 0;                     
329       G4ThreeVector maxLocalPoint;                
330                                                   
331       std::size_t limit = candidates.size();      
332       for (std::size_t i = 0 ; i < limit ; ++i    
333       {                                           
334         G4int candidate = candidates[i];          
335         // ignore the current component (that     
336         // numerically the propagated point wi    
337                                                   
338         G4VSolid& solid = *fSolids[candidate];    
339         const G4Transform3D& transform = fTran    
340                                                   
341         // The coordinates of the point are mo    
342         // intrinsic solid local frame:           
343         localPoint = GetLocalPoint(transform,     
344                                                   
345         // DistanceToOut at least for Trd some    
346         // even from points that are outside.     
347         // must currently be here, otherwise i    
348         // But it means it would be slower.       
349                                                   
350         if (solid.Inside(localPoint) != EInsid    
351         {                                         
352           notOutside = true;                      
353                                                   
354           localDirection = GetLocalVector(tran    
355                                                   
356           // propagate with solid.DistanceToOu    
357           G4double shift = solid.DistanceToOut    
358                                                   
359           if (maxDistance < shift)                
360           {                                       
361             maxDistance = shift;                  
362             maxCandidate = candidate;             
363             maxNormal = localNormal;              
364           }                                       
365         }                                         
366       }                                           
367                                                   
368       if (notOutside)                             
369       {                                           
370         const G4Transform3D& transform = fTran    
371                                                   
372         // convert from local normal              
373         if (aNormal != nullptr) *aNormal = Get    
374                                                   
375         distance += maxDistance;                  
376         currentPoint += maxDistance * directio    
377         if(maxDistance == 0.) ++count;            
378                                                   
379         // the current component will be ignor    
380         exclusion.SetBitNumber(maxCandidate);     
381         EInside location = InsideWithExclusion    
382                                                   
383         // perform a Inside                       
384         // it should be excluded current solid    
385         // we have to collect the maximum dist    
386         // such "maximum" candidate should be     
387         // candidates                             
388         if (location == EInside::kOutside)        
389         {                                         
390           // else return cumulated distances t    
391           // components                           
392           break;                                  
393         }                                         
394         // if inside another component, redo 1    
395         // DistanceToOut on top of the previou    
396                                                   
397         // and fill the candidates for the cor    
398         // exiting current component along dir    
399         candidates.clear();                       
400                                                   
401         fVoxels.GetCandidatesVoxelArray(curren    
402         exclusion.ResetBitNumber(maxCandidate)    
403       }                                           
404     }                                             
405     while  ((notOutside) && (count < numNodes)    
406   }                                               
407                                                   
408   return distance;                                
409 }                                                 
410                                                   
411 //____________________________________________    
412 EInside G4MultiUnion::InsideWithExclusion(cons    
413                                           G4Su    
414 {                                                 
415   // Classify point location with respect to s    
416   //  o eInside       - inside the solid          
417   //  o eSurface      - close to surface withi    
418   //  o eOutside      - outside the solid         
419                                                   
420   // Hitherto, it is considered that only para    
421   // can be added to the container                
422                                                   
423   // Implementation using voxelisation techniq    
424   // -----------------------------------------    
425                                                   
426   G4ThreeVector localPoint;                       
427   EInside location = EInside::kOutside;           
428                                                   
429   std::vector<G4int> candidates;                  
430   std::vector<G4MultiUnionSurface> surfaces;      
431                                                   
432   // TODO: test if it works well and if so mea    
433   // TODO: getPointIndex should not be used, i    
434   //       should be used                         
435   // TODO: than pass result to GetVoxel furthe    
436   // TODO: eventually GetVoxel should be inlin    
437   //       binary search is -1                    
438                                                   
439   G4int limit = fVoxels.GetCandidatesVoxelArra    
440   for (G4int i = 0 ; i < limit ; ++i)             
441   {                                               
442     G4int candidate = candidates[i];              
443     G4VSolid& solid = *fSolids[candidate];        
444     const G4Transform3D& transform = fTransfor    
445                                                   
446     // The coordinates of the point are modifi    
447     // solid local frame:                         
448     localPoint = GetLocalPoint(transform, aPoi    
449     location = solid.Inside(localPoint);          
450     if (location == EInside::kInside) return E    
451     else if (location == EInside::kSurface)       
452     {                                             
453       G4MultiUnionSurface surface;                
454       surface.point = localPoint;                 
455       surface.solid = &solid;                     
456       surfaces.push_back(surface);                
457     }                                             
458   }                                               
459                                                   
460   ////////////////////////////////////////////    
461   // Important comment: When two solids touch     
462   // surface, the surface points will be consi    
463   // located around will correspond to kInside    
464                                                   
465   std::size_t size = surfaces.size();             
466                                                   
467   if (size == 0)                                  
468   {                                               
469     return EInside::kOutside;                     
470   }                                               
471                                                   
472   for (std::size_t i = 0; i < size - 1; ++i)      
473   {                                               
474     G4MultiUnionSurface& left = surfaces[i];      
475     for (std::size_t j = i + 1; j < size; ++j)    
476     {                                             
477       G4MultiUnionSurface& right = surfaces[j]    
478       G4ThreeVector n, n2;                        
479       n = left.solid->SurfaceNormal(left.point    
480       n2 = right.solid->SurfaceNormal(right.po    
481       if ((n +  n2).mag2() < 1000 * kRadTolera    
482         return EInside::kInside;                  
483     }                                             
484   }                                               
485                                                   
486   return EInside::kSurface;                       
487 }                                                 
488                                                   
489 //____________________________________________    
490 EInside G4MultiUnion::Inside(const G4ThreeVect    
491 {                                                 
492   // Classify point location with respect to s    
493   //  o eInside       - inside the solid          
494   //  o eSurface      - close to surface withi    
495   //  o eOutside      - outside the solid         
496                                                   
497   // Hitherto, it is considered that only para    
498   // added to the container                       
499                                                   
500   // Implementation using voxelisation techniq    
501   // -----------------------------------------    
502                                                   
503   //  return InsideIterator(aPoint);              
504                                                   
505   EInside location = InsideWithExclusion(aPoin    
506   return location;                                
507 }                                                 
508                                                   
509 //____________________________________________    
510 EInside G4MultiUnion::InsideNoVoxels(const G4T    
511 {                                                 
512   G4ThreeVector localPoint;                       
513   EInside location = EInside::kOutside;           
514   G4int countSurface = 0;                         
515                                                   
516   auto numNodes = (G4int)fSolids.size();          
517   for (auto i = 0 ; i < numNodes ; ++i)           
518   {                                               
519     G4VSolid& solid = *fSolids[i];                
520     G4Transform3D transform = GetTransformatio    
521                                                   
522     // The coordinates of the point are modifi    
523     // intrinsic solid local frame:               
524     localPoint = GetLocalPoint(transform, aPoi    
525                                                   
526     location = solid.Inside(localPoint);          
527                                                   
528     if (location == EInside::kSurface)            
529       ++countSurface;                             
530                                                   
531     if (location == EInside::kInside) return E    
532   }                                               
533   if (countSurface != 0) return EInside::kSurf    
534   return EInside::kOutside;                       
535 }                                                 
536                                                   
537 //____________________________________________    
538 void G4MultiUnion::Extent(EAxis aAxis, G4doubl    
539 {                                                 
540   // Determines the bounding box for the consi    
541   G4ThreeVector min, max;                         
542                                                   
543   auto numNodes = (G4int)fSolids.size();          
544   for (auto i = 0 ; i < numNodes ; ++i)           
545   {                                               
546     G4VSolid& solid = *fSolids[i];                
547     G4Transform3D transform = GetTransformatio    
548     solid.BoundingLimits(min, max);               
549                                                   
550     TransformLimits(min, max, transform);         
551                                                   
552     if (i == 0)                                   
553     {                                             
554       switch (aAxis)                              
555       {                                           
556         case kXAxis:                              
557           aMin = min.x();                         
558           aMax = max.x();                         
559           break;                                  
560         case kYAxis:                              
561           aMin = min.y();                         
562           aMax = max.y();                         
563           break;                                  
564         case kZAxis:                              
565           aMin = min.z();                         
566           aMax = max.z();                         
567           break;                                  
568         default:                                  
569           break;                                  
570       }                                           
571     }                                             
572     else                                          
573     {                                             
574       // Determine the min/max on the consider    
575       switch (aAxis)                              
576       {                                           
577         case kXAxis:                              
578           if (min.x() < aMin)                     
579             aMin = min.x();                       
580           if (max.x() > aMax)                     
581             aMax = max.x();                       
582           break;                                  
583         case kYAxis:                              
584           if (min.y() < aMin)                     
585             aMin = min.y();                       
586           if (max.y() > aMax)                     
587             aMax = max.y();                       
588           break;                                  
589         case kZAxis:                              
590           if (min.z() < aMin)                     
591             aMin = min.z();                       
592           if (max.z() > aMax)                     
593             aMax = max.z();                       
594           break;                                  
595         default:                                  
596           break;                                  
597       }                                           
598     }                                             
599   }                                               
600 }                                                 
601                                                   
602 //____________________________________________    
603 void G4MultiUnion::BoundingLimits(G4ThreeVecto    
604                                   G4ThreeVecto    
605 {                                                 
606    Extent(kXAxis, aMin[0], aMax[0]);              
607    Extent(kYAxis, aMin[1], aMax[1]);              
608    Extent(kZAxis, aMin[2], aMax[2]);              
609 }                                                 
610                                                   
611 //____________________________________________    
612 G4bool                                            
613 G4MultiUnion::CalculateExtent(const EAxis pAxi    
614                               const G4VoxelLim    
615                               const G4AffineTr    
616                                     G4double&     
617 {                                                 
618   G4ThreeVector bmin, bmax;                       
619                                                   
620   // Get bounding box                             
621   BoundingLimits(bmin,bmax);                      
622                                                   
623   // Find extent                                  
624   G4BoundingEnvelope bbox(bmin,bmax);             
625   return bbox.CalculateExtent(pAxis,pVoxelLimi    
626 }                                                 
627                                                   
628 //____________________________________________    
629 G4ThreeVector G4MultiUnion::SurfaceNormal(cons    
630 {                                                 
631   // Computes the localNormal on a surface and    
632   // Must return a valid vector. (even if the     
633   //                                              
634   //   On an edge or corner, provide an averag    
635   //   tolerance                                  
636   //   NOTE: the tolerance value used in here     
637   //         tolerance - we will have to revis    
638                                                   
639   std::vector<G4int> candidates;                  
640   G4ThreeVector localPoint, normal, localNorma    
641   G4double safety = kInfinity;                    
642   G4int node = 0;                                 
643                                                   
644   ////////////////////////////////////////////    
645   // Important comment: Cases for which the po    
646   // on a vertice remain to be treated            
647                                                   
648   // determine weather we are in voxel area       
649   if (fVoxels.GetCandidatesVoxelArray(aPoint,     
650   {                                               
651     std::size_t limit = candidates.size();        
652     for (std::size_t i = 0 ; i < limit ; ++i)     
653     {                                             
654       G4int candidate = candidates[i];            
655       const G4Transform3D& transform = fTransf    
656                                                   
657       // The coordinates of the point are modi    
658       // solid local frame:                       
659       localPoint = GetLocalPoint(transform, aP    
660       G4VSolid& solid = *fSolids[candidate];      
661       EInside location = solid.Inside(localPoi    
662                                                   
663       if (location == EInside::kSurface)          
664       {                                           
665         // normal case when point is on surfac    
666         normal = GetGlobalVector(transform, so    
667         return normal.unit();                     
668       }                                           
669       else                                        
670       {                                           
671         // collect the smallest safety and rem    
672         G4double s = (location == EInside::kIn    
673                    ? solid.DistanceToOut(local    
674                    : solid.DistanceToIn(localP    
675         if (s < safety)                           
676         {                                         
677           safety = s;                             
678           node = candidate;                       
679         }                                         
680       }                                           
681     }                                             
682     // on none of the solids, the point was no    
683     G4VSolid& solid = *fSolids[node];             
684     const G4Transform3D& transform = fTransfor    
685     localPoint = GetLocalPoint(transform, aPoi    
686                                                   
687     normal = GetGlobalVector(transform, solid.    
688     return normal.unit();                         
689   }                                               
690   else                                            
691   {                                               
692     // for the case when point is certainly ou    
693                                                   
694     // find a solid in union with the smallest    
695     node = SafetyFromOutsideNumberNode(aPoint,    
696     G4VSolid& solid = *fSolids[node];             
697                                                   
698     const G4Transform3D& transform = fTransfor    
699     localPoint = GetLocalPoint(transform, aPoi    
700                                                   
701     // evaluate normal for point at this found    
702     // and transform multi-union coordinates      
703     normal = GetGlobalVector(transform, solid.    
704                                                   
705     return normal.unit();                         
706   }                                               
707 }                                                 
708                                                   
709 //____________________________________________    
710 G4double G4MultiUnion::DistanceToOut(const G4T    
711 {                                                 
712   // Estimates isotropic distance to the surfa    
713   // be either accurate or an underestimate.      
714   // Two modes: - default/fast mode, sacrifici    
715   //            - "precise" mode,  requests ac    
716                                                   
717   std::vector<G4int> candidates;                  
718   G4ThreeVector localPoint;                       
719   G4double safetyMin = kInfinity;                 
720                                                   
721   // In general, the value return by DistanceT    
722   // but only an undervalue (cf. overlaps)        
723   fVoxels.GetCandidatesVoxelArray(point, candi    
724                                                   
725   std::size_t limit = candidates.size();          
726   for (std::size_t i = 0; i < limit; ++i)         
727   {                                               
728     G4int candidate = candidates[i];              
729                                                   
730     // The coordinates of the point are modifi    
731     // solid local frame:                         
732     const G4Transform3D& transform = fTransfor    
733     localPoint = GetLocalPoint(transform, poin    
734     G4VSolid& solid = *fSolids[candidate];        
735     if (solid.Inside(localPoint) == EInside::k    
736     {                                             
737       G4double safety = solid.DistanceToOut(lo    
738       if (safetyMin > safety) safetyMin = safe    
739     }                                             
740   }                                               
741   if (safetyMin == kInfinity) safetyMin = 0; /    
742                                                   
743   return safetyMin;                               
744 }                                                 
745                                                   
746 //____________________________________________    
747 G4double G4MultiUnion::DistanceToIn(const G4Th    
748 {                                                 
749   // Estimates the isotropic safety from a poi    
750   // any of its surfaces. The algorithm may be    
751   // underestimate.                               
752                                                   
753   if (!fAccurate)  { return fVoxels.DistanceTo    
754                                                   
755   const std::vector<G4VoxelBox>& boxes = fVoxe    
756   G4double safetyMin = kInfinity;                 
757   G4ThreeVector localPoint;                       
758                                                   
759   std::size_t numNodes = fSolids.size();          
760   for (std::size_t j = 0; j < numNodes; ++j)      
761   {                                               
762     G4ThreeVector dxyz;                           
763     if (j > 0)                                    
764     {                                             
765       const G4ThreeVector& pos = boxes[j].pos;    
766       const G4ThreeVector& hlen = boxes[j].hle    
767       for (auto i = 0; i <= 2; ++i)               
768         // distance to middle point - hlength     
769         // of x,y,z                               
770         if ((dxyz[i] = std::abs(point[i] - pos    
771           continue;                               
772                                                   
773       G4double d2xyz = 0.;                        
774       for (auto i = 0; i <= 2; ++i)               
775         if (dxyz[i] > 0) d2xyz += dxyz[i] * dx    
776                                                   
777       // minimal distance is at least this, bu    
778       // we can stop if previous was already l    
779       // chance to be better tha previous valu    
780       if (d2xyz >= safetyMin * safetyMin)         
781       {                                           
782         continue;                                 
783       }                                           
784     }                                             
785     const G4Transform3D& transform = fTransfor    
786     localPoint = GetLocalPoint(transform, poin    
787     G4VSolid& solid = *fSolids[j];                
788                                                   
789     G4double safety = solid.DistanceToIn(local    
790     if (safety <= 0) return safety;               
791       // it was detected, that the point is no    
792     if (safetyMin > safety) safetyMin = safety    
793   }                                               
794   return safetyMin;                               
795 }                                                 
796                                                   
797 //____________________________________________    
798 G4double G4MultiUnion::GetSurfaceArea()           
799 {                                                 
800   if (fSurfaceArea == 0.0)                        
801   {                                               
802     fSurfaceArea = EstimateSurfaceArea(1000000    
803   }                                               
804   return fSurfaceArea;                            
805 }                                                 
806                                                   
807 //____________________________________________    
808 G4int G4MultiUnion::GetNumOfConstituents() con    
809 {                                                 
810   G4int num = 0;                                  
811   for (const auto solid : fSolids) { num += so    
812   return num;                                     
813 }                                                 
814                                                   
815 //____________________________________________    
816 G4bool G4MultiUnion::IsFaceted() const            
817 {                                                 
818   for (const auto solid : fSolids) { if (!soli    
819   return true;                                    
820 }                                                 
821                                                   
822 //____________________________________________    
823 void G4MultiUnion::Voxelize()                     
824 {                                                 
825   fVoxels.Voxelize(fSolids, fTransformObjs);      
826 }                                                 
827                                                   
828 //____________________________________________    
829 G4int G4MultiUnion::SafetyFromOutsideNumberNod    
830                                                   
831 {                                                 
832   // Method returning the closest node from a     
833   // G4MultiUnion.                                
834   // This is used to compute the normal in the    
835                                                   
836   const std::vector<G4VoxelBox>& boxes = fVoxe    
837   safetyMin = kInfinity;                          
838   std::size_t safetyNode = 0;                     
839   G4ThreeVector localPoint;                       
840                                                   
841   std::size_t numNodes = fSolids.size();          
842   for (std::size_t i = 0; i < numNodes; ++i)      
843   {                                               
844     G4double d2xyz = 0.;                          
845     G4double dxyz0 = std::abs(aPoint.x() - box    
846     if (dxyz0 > safetyMin) continue;              
847     G4double dxyz1 = std::abs(aPoint.y() - box    
848     if (dxyz1 > safetyMin) continue;              
849     G4double dxyz2 = std::abs(aPoint.z() - box    
850     if (dxyz2 > safetyMin) continue;              
851                                                   
852     if (dxyz0 > 0) d2xyz += dxyz0 * dxyz0;        
853     if (dxyz1 > 0) d2xyz += dxyz1 * dxyz1;        
854     if (dxyz2 > 0) d2xyz += dxyz2 * dxyz2;        
855     if (d2xyz >= safetyMin * safetyMin) contin    
856                                                   
857     G4VSolid& solid = *fSolids[i];                
858     const G4Transform3D& transform = fTransfor    
859     localPoint = GetLocalPoint(transform, aPoi    
860     fAccurate = true;                             
861     G4double safety = solid.DistanceToIn(local    
862     fAccurate = false;                            
863     if (safetyMin > safety)                       
864     {                                             
865       safetyMin = safety;                         
866       safetyNode = i;                             
867     }                                             
868   }                                               
869   return (G4int)safetyNode;                       
870 }                                                 
871                                                   
872 //____________________________________________    
873 void G4MultiUnion::TransformLimits(G4ThreeVect    
874                                    const G4Tra    
875 {                                                 
876   // The goal of this method is to convert the    
877   // (representing the bounding box of a given    
878   // to the main frame, using "transformation"    
879                                                   
880   G4ThreeVector vertices[8] =     // Deteminat    
881   {                               // the exten    
882     G4ThreeVector(min.x(), min.y(), min.z()),     
883     G4ThreeVector(min.x(), max.y(), min.z()),     
884     G4ThreeVector(max.x(), max.y(), min.z()),     
885     G4ThreeVector(max.x(), min.y(), min.z()),     
886     G4ThreeVector(min.x(), min.y(), max.z()),     
887     G4ThreeVector(min.x(), max.y(), max.z()),     
888     G4ThreeVector(max.x(), max.y(), max.z()),     
889     G4ThreeVector(max.x(), min.y(), max.z())      
890   };                                              
891                                                   
892   min.set(kInfinity,kInfinity,kInfinity);         
893   max.set(-kInfinity,-kInfinity,-kInfinity);      
894                                                   
895   // Loop on th vertices                          
896   G4int limit = sizeof(vertices) / sizeof(G4Th    
897   for (G4int i = 0 ; i < limit; ++i)              
898   {                                               
899     // From local frame to the global one:        
900     // Current positions on the three axis:       
901     G4ThreeVector current = GetGlobalPoint(tra    
902                                                   
903     // If need be, replacement of the min & ma    
904     if (current.x() > max.x()) max.setX(curren    
905     if (current.x() < min.x()) min.setX(curren    
906                                                   
907     if (current.y() > max.y()) max.setY(curren    
908     if (current.y() < min.y()) min.setY(curren    
909                                                   
910     if (current.z() > max.z()) max.setZ(curren    
911     if (current.z() < min.z()) min.setZ(curren    
912   }                                               
913 }                                                 
914                                                   
915 // Stream object contents to an output stream     
916 //____________________________________________    
917 std::ostream& G4MultiUnion::StreamInfo(std::os    
918 {                                                 
919   G4long oldprc = os.precision(16);               
920   os << "-------------------------------------    
921      << "                *** Dump for solid -     
922      << "                =====================    
923      << " Solid type: G4MultiUnion\n"             
924      << " Parameters: \n";                        
925      std::size_t numNodes = fSolids.size();       
926      for (std::size_t i = 0 ; i < numNodes ; +    
927      {                                            
928       G4VSolid& solid = *fSolids[i];              
929       solid.StreamInfo(os);                       
930       const G4Transform3D& transform = fTransf    
931       os << " Translation is " << transform.ge    
932       os << " Rotation is :" << " \n";            
933       os << " " << transform.getRotation() <<     
934      }                                            
935   os << "             \n"                         
936      << "-------------------------------------    
937   os.precision(oldprc);                           
938                                                   
939   return os;                                      
940 }                                                 
941                                                   
942 //____________________________________________    
943 G4ThreeVector G4MultiUnion::GetPointOnSurface(    
944 {                                                 
945   G4ThreeVector point;                            
946                                                   
947   G4long size = fSolids.size();                   
948                                                   
949   do                                              
950   {                                               
951     G4long rnd = G4RandFlat::shootInt(G4long(0    
952     G4VSolid& solid = *fSolids[rnd];              
953     point = solid.GetPointOnSurface();            
954     const G4Transform3D& transform = fTransfor    
955     point = GetGlobalPoint(transform, point);     
956   }                                               
957   while (Inside(point) != EInside::kSurface);     
958                                                   
959   return point;                                   
960 }                                                 
961                                                   
962 //____________________________________________    
963 void                                              
964 G4MultiUnion::DescribeYourselfTo ( G4VGraphics    
965 {                                                 
966   scene.AddSolid (*this);                         
967 }                                                 
968                                                   
969 //____________________________________________    
970 G4Polyhedron* G4MultiUnion::CreatePolyhedron()    
971 {                                                 
972   if (G4BooleanSolid::GetExternalBooleanProces    
973   {                                               
974     HepPolyhedronProcessor processor;             
975     HepPolyhedronProcessor::Operation operatio    
976                                                   
977     G4VSolid* solidA = GetSolid(0);               
978     const G4Transform3D transform0 = GetTransf    
979     G4DisplacedSolid dispSolidA("placedA", sol    
980                                                   
981     auto top = new G4Polyhedron(*dispSolidA.Ge    
982                                                   
983     for (G4int i = 1; i < GetNumberOfSolids();    
984     {                                             
985       G4VSolid* solidB = GetSolid(i);             
986       const G4Transform3D transform = GetTrans    
987       G4DisplacedSolid dispSolidB("placedB", s    
988       G4Polyhedron* operand = dispSolidB.GetPo    
989       processor.push_back(operation, *operand)    
990     }                                             
991                                                   
992     if (processor.execute(*top))                  
993     {                                             
994       return top;                                 
995     }                                             
996     else                                          
997     {                                             
998       return nullptr;                             
999     }                                             
1000   }                                              
1001   else                                           
1002   {                                              
1003     return G4BooleanSolid::GetExternalBoolean    
1004   }                                              
1005 }                                                
1006                                                  
1007 //___________________________________________    
1008 G4Polyhedron* G4MultiUnion::GetPolyhedron() c    
1009 {                                                
1010   if (fpPolyhedron == nullptr ||                 
1011       fRebuildPolyhedron ||                      
1012       fpPolyhedron->GetNumberOfRotationStepsA    
1013       fpPolyhedron->GetNumberOfRotationSteps(    
1014   {                                              
1015     G4AutoLock l(&polyhedronMutex);              
1016     delete fpPolyhedron;                         
1017     fpPolyhedron = CreatePolyhedron();           
1018     fRebuildPolyhedron = false;                  
1019     l.unlock();                                  
1020   }                                              
1021   return fpPolyhedron;                           
1022 }                                                
1023