Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/Boolean/src/G4UnionSolid.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/G4UnionSolid.cc (Version 11.3.0) and /geometry/solids/Boolean/src/G4UnionSolid.cc (Version 7.1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
 26 // Implementation of methods for the class G4U << 
 27 //                                                 23 //
 28 // 23.04.18 E.Tcherniaev: added extended BBox, <<  24 // $Id: G4UnionSolid.cc,v 1.29 2005/05/09 13:44:58 gcosmo Exp $
 29 // 17.03.17 E.Tcherniaev: revision of SurfaceN <<  25 // GEANT4 tag $Name: geant4-07-01 $
                                                   >>  26 //
                                                   >>  27 // Implementation of methods for the class G4IntersectionSolid
                                                   >>  28 //
                                                   >>  29 // History:
                                                   >>  30 //
 30 // 12.09.98 V.Grichine: first implementation       31 // 12.09.98 V.Grichine: first implementation
                                                   >>  32 // 28.11.98 V.Grichine: fix while loops in DistToIn/Out 
                                                   >>  33 // 27.07.99 V.Grichine: modifications in DistToOut(p,v,...), while -> do-while
                                                   >>  34 // 16.03.01 V.Grichine: modifications in CalculateExtent()
                                                   >>  35 //
 31 // -------------------------------------------     36 // --------------------------------------------------------------------
 32                                                    37 
 33 #include <sstream>                             << 
 34                                                << 
 35 #include "G4UnionSolid.hh"                         38 #include "G4UnionSolid.hh"
 36                                                    39 
 37 #include "G4SystemOfUnits.hh"                  << 
 38 #include "G4VoxelLimits.hh"                        40 #include "G4VoxelLimits.hh"
 39 #include "G4VPVParameterisation.hh"                41 #include "G4VPVParameterisation.hh"
 40 #include "G4GeometryTolerance.hh"              << 
 41                                                    42 
 42 #include "G4VGraphicsScene.hh"                     43 #include "G4VGraphicsScene.hh"
 43 #include "G4Polyhedron.hh"                         44 #include "G4Polyhedron.hh"
 44 #include "G4PolyhedronArbitrary.hh"            <<  45 #include "G4NURBS.hh"
 45 #include "HepPolyhedronProcessor.h"            <<  46 // #include "G4NURBSbox.hh"
 46                                                << 
 47 #include "G4IntersectionSolid.hh"              << 
 48                                                    47 
 49 ////////////////////////////////////////////// <<  48 ///////////////////////////////////////////////////////////////////
 50 //                                                 49 //
 51 // Transfer all data members to G4BooleanSolid     50 // Transfer all data members to G4BooleanSolid which is responsible
 52 // for them. pName will be in turn sent to G4V     51 // for them. pName will be in turn sent to G4VSolid
 53                                                    52 
 54 G4UnionSolid:: G4UnionSolid( const G4String& p     53 G4UnionSolid:: G4UnionSolid( const G4String& pName,
 55                                    G4VSolid* p     54                                    G4VSolid* pSolidA ,
 56                                    G4VSolid* p     55                                    G4VSolid* pSolidB )
 57   : G4BooleanSolid(pName,pSolidA,pSolidB)          56   : G4BooleanSolid(pName,pSolidA,pSolidB)
 58 {                                                  57 {
 59   Init();                                      << 
 60 }                                                  58 }
 61                                                    59 
 62 ////////////////////////////////////////////// <<  60 /////////////////////////////////////////////////////////////////////
 63 //                                                 61 //
 64 // Constructor                                     62 // Constructor
 65                                                    63  
 66 G4UnionSolid::G4UnionSolid( const G4String& pN     64 G4UnionSolid::G4UnionSolid( const G4String& pName,
 67                                   G4VSolid* pS     65                                   G4VSolid* pSolidA ,
 68                                   G4VSolid* pS     66                                   G4VSolid* pSolidB ,
 69                                   G4RotationMa     67                                   G4RotationMatrix* rotMatrix,
 70                             const G4ThreeVecto     68                             const G4ThreeVector& transVector )
 71   : G4BooleanSolid(pName,pSolidA,pSolidB,rotMa     69   : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
 72                                                    70 
 73 {                                                  71 {
 74   Init();                                      << 
 75 }                                                  72 }
 76                                                    73 
 77 ////////////////////////////////////////////// <<  74 ///////////////////////////////////////////////////////////
 78 //                                                 75 //
 79 // Constructor                                     76 // Constructor
 80                                                    77  
 81 G4UnionSolid::G4UnionSolid( const G4String& pN     78 G4UnionSolid::G4UnionSolid( const G4String& pName,
 82                                   G4VSolid* pS     79                                   G4VSolid* pSolidA ,
 83                                   G4VSolid* pS     80                                   G4VSolid* pSolidB ,
 84                             const G4Transform3     81                             const G4Transform3D& transform )
 85   : G4BooleanSolid(pName,pSolidA,pSolidB,trans     82   : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
 86 {                                                  83 {
 87   Init();                                      << 
 88 }                                                  84 } 
 89                                                    85 
 90 ////////////////////////////////////////////// << 
 91 //                                             << 
 92 // Fake default constructor - sets only member << 
 93 //                            for usage restri << 
 94                                                << 
 95 G4UnionSolid::G4UnionSolid( __void__& a )      << 
 96   : G4BooleanSolid(a)                          << 
 97 {                                              << 
 98 }                                              << 
 99                                                    86 
100 ////////////////////////////////////////////// <<  87 ///////////////////////////////////////////////////////////
101 //                                                 88 //
102 // Destructor                                      89 // Destructor
103                                                    90 
104 G4UnionSolid::~G4UnionSolid()                      91 G4UnionSolid::~G4UnionSolid()
105 = default;                                     << 
106                                                << 
107 ////////////////////////////////////////////// << 
108 //                                             << 
109 // Copy constructor                            << 
110                                                << 
111 G4UnionSolid::G4UnionSolid(const G4UnionSolid& << 
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  << 
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 {                                                  92 {
146   G4ThreeVector pdelta(kCarTolerance,kCarToler << 
147   G4ThreeVector pmin, pmax;                    << 
148   BoundingLimits(pmin, pmax);                  << 
149   fPMin = pmin - pdelta;                       << 
150   fPMax = pmax + pdelta;                       << 
151   halfCarTolerance=0.5*kCarTolerance;          << 
152 }                                                  93 }
153                                                    94 
154 ////////////////////////////////////////////// <<  95 ///////////////////////////////////////////////////////////////
155 //                                                 96 //
156 // Get bounding box                            << 
157                                                << 
158 void G4UnionSolid::BoundingLimits(G4ThreeVecto << 
159                                   G4ThreeVecto << 
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 << 
176   {                                            << 
177     std::ostringstream message;                << 
178     message << "Bad bounding box (min >= max)  << 
179             << GetName() << " !"               << 
180             << "\npMin = " << pMin             << 
181             << "\npMax = " << pMax;            << 
182     G4Exception("G4UnionSolid::BoundingLimits( << 
183                 JustWarning, message);         << 
184     DumpInfo();                                << 
185   }                                            << 
186 }                                              << 
187                                                << 
188 ////////////////////////////////////////////// << 
189 //                                                 97 //
190 // Calculate extent under transform and specif << 
191                                                    98      
192 G4bool                                             99 G4bool 
193 G4UnionSolid::CalculateExtent( const EAxis pAx    100 G4UnionSolid::CalculateExtent( const EAxis pAxis,
194                                const G4VoxelLi    101                                const G4VoxelLimits& pVoxelLimit,
195                                const G4AffineT    102                                const G4AffineTransform& pTransform,
196                                      G4double&    103                                      G4double& pMin,
197                                      G4double&    104                                      G4double& pMax ) const 
198 {                                                 105 {
199   G4bool   touchesA, touchesB, out ;              106   G4bool   touchesA, touchesB, out ;
200   G4double minA =  kInfinity, minB =  kInfinit    107   G4double minA =  kInfinity, minB =  kInfinity, 
201            maxA = -kInfinity, maxB = -kInfinit    108            maxA = -kInfinity, maxB = -kInfinity; 
202                                                   109 
203   touchesA = fPtrSolidA->CalculateExtent( pAxi    110   touchesA = fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit, 
204                                           pTra    111                                           pTransform, minA, maxA);
205   touchesB = fPtrSolidB->CalculateExtent( pAxi << 112   touchesB= fPtrSolidB->CalculateExtent( pAxis, pVoxelLimit, 
206                                           pTra << 113                                          pTransform, minB, maxB);
207   if( touchesA || touchesB )                      114   if( touchesA || touchesB )
208   {                                               115   {
209     pMin = std::min( minA, minB );                116     pMin = std::min( minA, minB ); 
210     pMax = std::max( maxA, maxB );                117     pMax = std::max( maxA, maxB );
211     out  = true ;                                 118     out  = true ; 
212   }                                               119   }
213   else                                         << 120   else out = false ;
214   {                                            << 
215     out = false ;                              << 
216   }                                            << 
217                                                   121 
218   return out ;  // It exists in this slice if     122   return out ;  // It exists in this slice if either one does.
219 }                                                 123 }
220                                                   124  
221 ////////////////////////////////////////////// << 125 /////////////////////////////////////////////////////
222 //                                                126 //
223 // Important comment: When solids A and B touc    127 // Important comment: When solids A and B touch together along flat
224 // surface the surface points will be consider    128 // surface the surface points will be considered as kSurface, while points 
225 // located around will correspond to kInside      129 // located around will correspond to kInside
226                                                   130 
227 EInside G4UnionSolid::Inside( const G4ThreeVec    131 EInside G4UnionSolid::Inside( const G4ThreeVector& p ) const
228 {                                                 132 {
229   if (std::max(p.z()-fPMax.z(), fPMin.z()-p.z( << 
230                                                << 
231   EInside positionA = fPtrSolidA->Inside(p);      133   EInside positionA = fPtrSolidA->Inside(p);
232   if (positionA == kInside)  { return position << 
233   EInside positionB = fPtrSolidB->Inside(p);      134   EInside positionB = fPtrSolidB->Inside(p);
234   if (positionA == kOutside) { return position << 
235                                                << 
236   if (positionB == kInside)  { return position << 
237   if (positionB == kOutside) { return position << 
238                                                   135 
239   // Both points are on surface                << 136   if( positionA == kInside  || positionB == kInside   ||
240   //                                           << 137     ( positionA == kSurface && positionB == kSurface &&
241   static const G4double rtol                   << 138         ( fPtrSolidA->SurfaceNormal(p) + 
242     = 1000*G4GeometryTolerance::GetInstance()- << 139           fPtrSolidB->SurfaceNormal(p) ).mag2() < 
243                                                << 140           1000*kRadTolerance ) )                              return kInside;
244   return ((fPtrSolidA->SurfaceNormal(p) +      << 141   else
245            fPtrSolidB->SurfaceNormal(p)).mag2( << 142   {
                                                   >> 143     if( ( positionA != kInside  && positionB == kSurface ) ||
                                                   >> 144         ( positionB != kInside  && positionA == kSurface ) ||
                                                   >> 145         ( positionA == kSurface && positionB == kSurface )    ) return kSurface;
                                                   >> 146     else                                                        return kOutside;
                                                   >> 147   }
246 }                                                 148 }
247                                                   149 
248 ////////////////////////////////////////////// << 150 //////////////////////////////////////////////////////////////
                                                   >> 151 //
249 //                                                152 //
250 // Get surface normal                          << 
251                                                   153 
252 G4ThreeVector                                     154 G4ThreeVector 
253 G4UnionSolid::SurfaceNormal( const G4ThreeVect    155 G4UnionSolid::SurfaceNormal( const G4ThreeVector& p ) const 
254 {                                                 156 {
255   EInside positionA = fPtrSolidA->Inside(p);   << 157     G4ThreeVector normal;
256   EInside positionB = fPtrSolidB->Inside(p);   << 
257                                                << 
258   if (positionA == kSurface &&                 << 
259       positionB == kOutside) return fPtrSolidA << 
260                                                   158 
261   if (positionA == kOutside &&                 << 159 #ifdef G4BOOLDEBUG
262       positionB == kSurface) return fPtrSolidB << 160     if( Inside(p) == kOutside )
                                                   >> 161     {
                                                   >> 162       G4cout << "WARNING - Invalid call in "
                                                   >> 163              << "G4UnionSolid::SurfaceNormal(p)" << G4endl
                                                   >> 164              << "  Point p is outside !" << G4endl;
                                                   >> 165       G4cout << "          p = " << p << G4endl;
                                                   >> 166       G4cerr << "WARNING - Invalid call in "
                                                   >> 167              << "G4UnionSolid::SurfaceNormal(p)" << G4endl
                                                   >> 168              << "  Point p is outside !" << G4endl;
                                                   >> 169       G4cerr << "          p = " << p << G4endl;
                                                   >> 170     }
                                                   >> 171 #endif
263                                                   172 
264   if (positionA == kSurface &&                 << 173     if(fPtrSolidA->Inside(p) == kSurface && fPtrSolidB->Inside(p) != kInside) 
265       positionB == kSurface)                   << 
266   {                                            << 
267     if (Inside(p) == kSurface)                 << 
268     {                                             174     {
269       G4ThreeVector normalA = fPtrSolidA->Surf << 175        normal= fPtrSolidA->SurfaceNormal(p) ;
270       G4ThreeVector normalB = fPtrSolidB->Surf << 
271       return (normalA + normalB).unit();       << 
272     }                                             176     }
273   }                                            << 177     else if(fPtrSolidB->Inside(p) == kSurface && 
                                                   >> 178             fPtrSolidA->Inside(p) != kInside)
                                                   >> 179     {
                                                   >> 180        normal= fPtrSolidB->SurfaceNormal(p) ;
                                                   >> 181     }
                                                   >> 182     else 
                                                   >> 183     {
                                                   >> 184       normal= fPtrSolidA->SurfaceNormal(p) ;
274 #ifdef G4BOOLDEBUG                                185 #ifdef G4BOOLDEBUG
275   G4String surf[3] = { "OUTSIDE", "SURFACE", " << 186       if(Inside(p)==kInside)
276   std::ostringstream message;                  << 187       {
277   G4int oldprc = message.precision(16);        << 188         G4cout << "WARNING - Invalid call in "
278   message << "Invalid call of SurfaceNormal(p) << 189              << "G4UnionSolid::SurfaceNormal(p)" << G4endl
279           << GetName() << " !"                 << 190              << "  Point p is inside !" << G4endl;
280           << "\nPoint p" << p << " is " << sur << 191         G4cout << "          p = " << p << G4endl;
281   message.precision(oldprc);                   << 192         G4cerr << "WARNING - Invalid call in "
282   G4Exception("G4UnionSolid::SurfaceNormal()", << 193              << "G4UnionSolid::SurfaceNormal(p)" << G4endl
283               JustWarning, message);           << 194              << "  Point p is inside !" << G4endl;
                                                   >> 195         G4cerr << "          p = " << p << G4endl;
                                                   >> 196       }
284 #endif                                            197 #endif
285   return fPtrSolidA->SurfaceNormal(p);         << 198     }
                                                   >> 199     return normal;
286 }                                                 200 }
287                                                   201 
288 ////////////////////////////////////////////// << 202 /////////////////////////////////////////////////////////////
289 //                                                203 //
290 // The same algorithm as in DistanceToIn(p)       204 // The same algorithm as in DistanceToIn(p)
291                                                   205 
292 G4double                                          206 G4double 
293 G4UnionSolid::DistanceToIn( const G4ThreeVecto    207 G4UnionSolid::DistanceToIn( const G4ThreeVector& p,
294                             const G4ThreeVecto << 208                                    const G4ThreeVector& v  ) const 
295 {                                                 209 {
296 #ifdef G4BOOLDEBUG                                210 #ifdef G4BOOLDEBUG
297   if( Inside(p) == kInside )                      211   if( Inside(p) == kInside )
298   {                                               212   {
299     G4cout << "WARNING - Invalid call in "        213     G4cout << "WARNING - Invalid call in "
300            << "G4UnionSolid::DistanceToIn(p,v)    214            << "G4UnionSolid::DistanceToIn(p,v)" << G4endl
301            << "  Point p is inside !" << G4end    215            << "  Point p is inside !" << G4endl;
302     G4cout << "          p = " << p << G4endl;    216     G4cout << "          p = " << p << G4endl;
303     G4cout << "          v = " << v << G4endl;    217     G4cout << "          v = " << v << G4endl;
304     G4cerr << "WARNING - Invalid call in "        218     G4cerr << "WARNING - Invalid call in "
305            << "G4UnionSolid::DistanceToIn(p,v)    219            << "G4UnionSolid::DistanceToIn(p,v)" << G4endl
306            << "  Point p is inside !" << G4end    220            << "  Point p is inside !" << G4endl;
307     G4cerr << "          p = " << p << G4endl;    221     G4cerr << "          p = " << p << G4endl;
308     G4cerr << "          v = " << v << G4endl;    222     G4cerr << "          v = " << v << G4endl;
309   }                                               223   }
310 #endif                                            224 #endif
311                                                   225 
312   return std::min(fPtrSolidA->DistanceToIn(p,v    226   return std::min(fPtrSolidA->DistanceToIn(p,v),
313                   fPtrSolidB->DistanceToIn(p,v << 227                     fPtrSolidB->DistanceToIn(p,v) ) ;
314 }                                                 228 }
315                                                   229 
316 ////////////////////////////////////////////// << 230 ////////////////////////////////////////////////////////
317 //                                                231 //
318 // Approximate nearest distance from the point    232 // Approximate nearest distance from the point p to the union of
319 // two solids                                     233 // two solids
320                                                   234 
321 G4double                                          235 G4double 
322 G4UnionSolid::DistanceToIn( const G4ThreeVecto << 236 G4UnionSolid::DistanceToIn( const G4ThreeVector& p) const 
323 {                                                 237 {
324 #ifdef G4BOOLDEBUG                                238 #ifdef G4BOOLDEBUG
325   if( Inside(p) == kInside )                      239   if( Inside(p) == kInside )
326   {                                               240   {
327     G4cout << "WARNING - Invalid call in "        241     G4cout << "WARNING - Invalid call in "
328            << "G4UnionSolid::DistanceToIn(p)"     242            << "G4UnionSolid::DistanceToIn(p)" << G4endl
329            << "  Point p is inside !" << G4end    243            << "  Point p is inside !" << G4endl;
330     G4cout << "          p = " << p << G4endl;    244     G4cout << "          p = " << p << G4endl;
331     G4cerr << "WARNING - Invalid call in "        245     G4cerr << "WARNING - Invalid call in "
332            << "G4UnionSolid::DistanceToIn(p)"     246            << "G4UnionSolid::DistanceToIn(p)" << G4endl
333            << "  Point p is inside !" << G4end    247            << "  Point p is inside !" << G4endl;
334     G4cerr << "          p = " << p << G4endl;    248     G4cerr << "          p = " << p << G4endl;
335   }                                               249   }
336 #endif                                            250 #endif
337   G4double distA = fPtrSolidA->DistanceToIn(p)    251   G4double distA = fPtrSolidA->DistanceToIn(p) ;
338   G4double distB = fPtrSolidB->DistanceToIn(p)    252   G4double distB = fPtrSolidB->DistanceToIn(p) ;
339   G4double safety = std::min(distA,distB) ;       253   G4double safety = std::min(distA,distB) ;
340   if(safety < 0.0) safety = 0.0 ;                 254   if(safety < 0.0) safety = 0.0 ;
341   return safety ;                                 255   return safety ;
342 }                                                 256 }
343                                                   257 
344 ////////////////////////////////////////////// << 258 //////////////////////////////////////////////////////////
345 //                                                259 //
346 // The same algorithm as DistanceToOut(p)         260 // The same algorithm as DistanceToOut(p)
347                                                   261 
348 G4double                                          262 G4double 
349 G4UnionSolid::DistanceToOut( const G4ThreeVect    263 G4UnionSolid::DistanceToOut( const G4ThreeVector& p,
350                              const G4ThreeVect << 264            const G4ThreeVector& v,
351                              const G4bool calc << 265            const G4bool calcNorm,
352                                    G4bool* val << 266                  G4bool *validNorm,
353                                    G4ThreeVect << 267                  G4ThreeVector *n      ) const 
354 {                                                 268 {
355   G4double  dist = 0.0, disTmp = 0.0 ;            269   G4double  dist = 0.0, disTmp = 0.0 ;
356   G4ThreeVector normTmp;                          270   G4ThreeVector normTmp;
357   G4ThreeVector* nTmp = &normTmp;              << 271   G4ThreeVector* nTmp= &normTmp;
358                                                   272 
359   if( Inside(p) == kOutside )                     273   if( Inside(p) == kOutside )
360   {                                               274   {
361 #ifdef G4BOOLDEBUG                                275 #ifdef G4BOOLDEBUG
362       G4cout << "Position:"  << G4endl << G4en    276       G4cout << "Position:"  << G4endl << G4endl;
363       G4cout << "p.x() = "   << p.x()/mm << "     277       G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl;
364       G4cout << "p.y() = "   << p.y()/mm << "     278       G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl;
365       G4cout << "p.z() = "   << p.z()/mm << "     279       G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl;
366       G4cout << "Direction:" << G4endl << G4en    280       G4cout << "Direction:" << G4endl << G4endl;
367       G4cout << "v.x() = "   << v.x() << G4end    281       G4cout << "v.x() = "   << v.x() << G4endl;
368       G4cout << "v.y() = "   << v.y() << G4end    282       G4cout << "v.y() = "   << v.y() << G4endl;
369       G4cout << "v.z() = "   << v.z() << G4end    283       G4cout << "v.z() = "   << v.z() << G4endl << G4endl;
370       G4cout << "WARNING - Invalid call in "      284       G4cout << "WARNING - Invalid call in "
371              << "G4UnionSolid::DistanceToOut(p    285              << "G4UnionSolid::DistanceToOut(p,v)" << G4endl
372              << "  Point p is outside !" << G4    286              << "  Point p is outside !" << G4endl;
373       G4cout << "          p = " << p << G4end    287       G4cout << "          p = " << p << G4endl;
374       G4cout << "          v = " << v << G4end    288       G4cout << "          v = " << v << G4endl;
375       G4cerr << "WARNING - Invalid call in "      289       G4cerr << "WARNING - Invalid call in "
376              << "G4UnionSolid::DistanceToOut(p    290              << "G4UnionSolid::DistanceToOut(p,v)" << G4endl
377              << "  Point p is outside !" << G4    291              << "  Point p is outside !" << G4endl;
378       G4cerr << "          p = " << p << G4end    292       G4cerr << "          p = " << p << G4endl;
379       G4cerr << "          v = " << v << G4end    293       G4cerr << "          v = " << v << G4endl;
380 #endif                                            294 #endif
381   }                                               295   }
382   else                                            296   else
383   {                                               297   {
384     EInside positionA = fPtrSolidA->Inside(p)     298     EInside positionA = fPtrSolidA->Inside(p) ;
                                                   >> 299     // EInside positionB = fPtrSolidB->Inside(p) ;
385                                                   300 
386     if( positionA != kOutside )                   301     if( positionA != kOutside )
387     {                                             302     { 
388       do  // Loop checking, 13.08.2015, G.Cosm << 303       do
389       {                                           304       {
390         disTmp = fPtrSolidA->DistanceToOut(p+d    305         disTmp = fPtrSolidA->DistanceToOut(p+dist*v,v,calcNorm,
391                                            val << 306                                            validNorm,nTmp)        ;
392         dist += disTmp ;                          307         dist += disTmp ;
393                                                   308 
394         if(fPtrSolidB->Inside(p+dist*v) != kOu    309         if(fPtrSolidB->Inside(p+dist*v) != kOutside)
395         {                                         310         { 
396           disTmp = fPtrSolidB->DistanceToOut(p    311           disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v,calcNorm,
397                                              v << 312                                             validNorm,nTmp)         ;
398           dist += disTmp ;                        313           dist += disTmp ;
399         }                                         314         }
400       }                                           315       }
401       while( (fPtrSolidA->Inside(p+dist*v) !=  << 316       //     while( Inside(p+dist*v) == kInside ) ;
402           && (disTmp > halfCarTolerance) );    << 317            while( fPtrSolidA->Inside(p+dist*v) != kOutside && 
                                                   >> 318                   disTmp > 0.5*kCarTolerance ) ;
403     }                                             319     }
404     else // if( positionB != kOutside )           320     else // if( positionB != kOutside )
405     {                                             321     {
406       do  // Loop checking, 13.08.2015, G.Cosm << 322       do
407       {                                           323       {
408         disTmp = fPtrSolidB->DistanceToOut(p+d    324         disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v,calcNorm,
409                                            val << 325                                            validNorm,nTmp)        ; 
410         dist += disTmp ;                          326         dist += disTmp ;
411                                                   327 
412         if(fPtrSolidA->Inside(p+dist*v) != kOu    328         if(fPtrSolidA->Inside(p+dist*v) != kOutside)
413         {                                         329         { 
414           disTmp = fPtrSolidA->DistanceToOut(p    330           disTmp = fPtrSolidA->DistanceToOut(p+dist*v,v,calcNorm,
415                                              v << 331                                             validNorm,nTmp)         ;
416           dist += disTmp ;                        332           dist += disTmp ;
417         }                                         333         }
418       }                                           334       }
419       while( (fPtrSolidB->Inside(p+dist*v) !=  << 335       //  while( Inside(p+dist*v) == kInside ) ;
420           && (disTmp > halfCarTolerance) );    << 336         while( (fPtrSolidB->Inside(p+dist*v) != kOutside)
                                                   >> 337             && (disTmp > 0.5*kCarTolerance) ) ;
421     }                                             338     }
422   }                                               339   }
423   if( calcNorm )                                  340   if( calcNorm )
424   {                                               341   { 
425      *validNorm = false ;                         342      *validNorm = false ;
426      *n         = *nTmp ;                         343      *n         = *nTmp ;   
427   }                                               344   }
428   return dist ;                                   345   return dist ;
429 }                                                 346 }
430                                                   347 
431 ////////////////////////////////////////////// << 348 //////////////////////////////////////////////////////////////
432 //                                                349 //
433 // Inverted algorithm of DistanceToIn(p)          350 // Inverted algorithm of DistanceToIn(p)
434                                                   351 
435 G4double                                          352 G4double 
436 G4UnionSolid::DistanceToOut( const G4ThreeVect    353 G4UnionSolid::DistanceToOut( const G4ThreeVector& p ) const 
437 {                                                 354 {
438   G4double distout = 0.0;                         355   G4double distout = 0.0;
439   if( Inside(p) == kOutside )                     356   if( Inside(p) == kOutside )
440   {                                               357   {
441 #ifdef G4BOOLDEBUG                                358 #ifdef G4BOOLDEBUG
442     G4cout << "WARNING - Invalid call in "        359     G4cout << "WARNING - Invalid call in "
443            << "G4UnionSolid::DistanceToOut(p)"    360            << "G4UnionSolid::DistanceToOut(p)" << G4endl
444            << "  Point p is outside !" << G4en    361            << "  Point p is outside !" << G4endl;
445     G4cout << "          p = " << p << G4endl;    362     G4cout << "          p = " << p << G4endl;
446     G4cerr << "WARNING - Invalid call in "        363     G4cerr << "WARNING - Invalid call in "
447            << "G4UnionSolid::DistanceToOut(p)"    364            << "G4UnionSolid::DistanceToOut(p)" << G4endl
448            << "  Point p is outside !" << G4en    365            << "  Point p is outside !" << G4endl;
449     G4cerr << "          p = " << p << G4endl;    366     G4cerr << "          p = " << p << G4endl;
450 #endif                                            367 #endif
451   }                                               368   }
452   else                                            369   else
453   {                                               370   {
454     EInside positionA = fPtrSolidA->Inside(p)     371     EInside positionA = fPtrSolidA->Inside(p) ;
455     EInside positionB = fPtrSolidB->Inside(p)     372     EInside positionB = fPtrSolidB->Inside(p) ;
456                                                   373   
457     //  Is this equivalent ??                     374     //  Is this equivalent ??
458     //    if( ! (  (positionA == kOutside)) &&    375     //    if( ! (  (positionA == kOutside)) && 
459     //             (positionB == kOutside))  )    376     //             (positionB == kOutside))  ) 
460     if((positionA == kInside  && positionB ==     377     if((positionA == kInside  && positionB == kInside  ) ||
461        (positionA == kInside  && positionB ==     378        (positionA == kInside  && positionB == kSurface ) ||
462        (positionA == kSurface && positionB ==     379        (positionA == kSurface && positionB == kInside  )     )
463     {                                             380     {     
464       distout= std::max(fPtrSolidA->DistanceTo    381       distout= std::max(fPtrSolidA->DistanceToOut(p),
465                         fPtrSolidB->DistanceTo << 382                           fPtrSolidB->DistanceToOut(p) ) ;
466     }                                             383     }
467     else                                          384     else
468     {                                             385     {
469       if(positionA == kOutside)                   386       if(positionA == kOutside)
470       {                                           387       {
471         distout= fPtrSolidB->DistanceToOut(p)     388         distout= fPtrSolidB->DistanceToOut(p) ;
472       }                                           389       }
473       else                                        390       else
474       {                                           391       {
475         distout= fPtrSolidA->DistanceToOut(p)     392         distout= fPtrSolidA->DistanceToOut(p) ;
476       }                                           393       }
477     }                                             394     }
478   }                                               395   }
479   return distout;                                 396   return distout;
480 }                                                 397 }
481                                                   398 
482 ////////////////////////////////////////////// << 399 //////////////////////////////////////////////////////////////
                                                   >> 400 //
483 //                                                401 //
484 // GetEntityType                               << 
485                                                   402 
486 G4GeometryType G4UnionSolid::GetEntityType() c    403 G4GeometryType G4UnionSolid::GetEntityType() const 
487 {                                                 404 {
488   return {"G4UnionSolid"};                     << 405   return G4String("G4UnionSolid");
489 }                                                 406 }
490                                                   407 
491 ////////////////////////////////////////////// << 408 //////////////////////////////////////////////////////////////
492 //                                                409 //
493 // Make a clone of the object                  << 
494                                                << 
495 G4VSolid* G4UnionSolid::Clone() const          << 
496 {                                              << 
497   return new G4UnionSolid(*this);              << 
498 }                                              << 
499                                                << 
500 ////////////////////////////////////////////// << 
501 //                                                410 //
502 // ComputeDimensions                           << 
503                                                   411 
504 void                                              412 void 
505 G4UnionSolid::ComputeDimensions(       G4VPVPa    413 G4UnionSolid::ComputeDimensions(       G4VPVParameterisation*,
506                                  const G4int,     414                                  const G4int,
507                                  const G4VPhys    415                                  const G4VPhysicalVolume* ) 
508 {                                                 416 {
509 }                                                 417 }
510                                                   418 
511 ////////////////////////////////////////////// << 419 /////////////////////////////////////////////////
512 //                                                420 //
513 // DescribeYourselfTo                          << 421 //                    
514                                                   422 
515 void                                              423 void 
516 G4UnionSolid::DescribeYourselfTo ( G4VGraphics    424 G4UnionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 
517 {                                                 425 {
518   scene.AddSolid (*this);                         426   scene.AddSolid (*this);
519 }                                                 427 }
520                                                   428 
521 ////////////////////////////////////////////// << 429 ////////////////////////////////////////////////////
                                                   >> 430 //
522 //                                                431 //
523 // CreatePolyhedron                            << 
524                                                   432 
525 G4Polyhedron*                                     433 G4Polyhedron* 
526 G4UnionSolid::CreatePolyhedron () const        << 434 G4UnionSolid::CreatePolyhedron () const 
527 {                                                 435 {
528   if (fExternalBoolProcessor == nullptr)       << 436   G4Polyhedron* pA = fPtrSolidA->GetPolyhedron();
529   {                                            << 437   G4Polyhedron* pB = fPtrSolidB->GetPolyhedron();
530     HepPolyhedronProcessor processor;          << 438   G4Polyhedron* resultant = new G4Polyhedron (pA->add(*pB));
531     // Stack components and components of comp << 439   return resultant;
532     // See G4BooleanSolid::StackPolyhedron     << 
533     G4Polyhedron* top = StackPolyhedron(proces << 
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(thi << 
547   }                                            << 
548 }                                                 440 }
549                                                   441 
550 ////////////////////////////////////////////// << 442 /////////////////////////////////////////////////////////
                                                   >> 443 //
551 //                                                444 //
552 // GetCubicVolume                              << 
553                                                   445 
554 G4double G4UnionSolid::GetCubicVolume()        << 446 G4NURBS*      
                                                   >> 447 G4UnionSolid::CreateNURBS      () const 
555 {                                                 448 {
556   if( fCubicVolume >= 0. )                     << 449   // Take into account boolean operation - see CreatePolyhedron.
557   {                                            << 450   // return new G4NURBSbox (1.0, 1.0, 1.0);
558     return fCubicVolume;                       << 451   return 0;
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() >= bm << 
565      bminB.x() >= bmaxA.x() || bminB.y() >= bm << 
566                                                << 
567   if (noIntersection)                          << 
568   {                                            << 
569     fCubicVolume = fPtrSolidA->GetCubicVolume( << 
570   }                                            << 
571   else                                         << 
572   {                                            << 
573     if (GetNumOfConstituents() > 10)           << 
574     {                                          << 
575       fCubicVolume = G4BooleanSolid::GetCubicV << 
576     }                                          << 
577     else                                       << 
578     {                                          << 
579       G4IntersectionSolid intersectVol("Tempor << 
580                                        fPtrSol << 
581       intersectVol.SetCubVolStatistics(GetCubV << 
582       intersectVol.SetCubVolEpsilon(GetCubVolE << 
583                                                << 
584       fCubicVolume = fPtrSolidA->GetCubicVolum << 
585         - intersectVol.GetCubicVolume();       << 
586     }                                          << 
587   }                                            << 
588   return fCubicVolume;                         << 
589 }                                                 452 }
590                                                   453