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