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 8.1.p1)


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