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