| Geant4 Cross Reference |
>> 1 // This code implementation is the intellectual property of
>> 2 // the GEANT4 collaboration.
1 // 3 //
2 // ******************************************* << 4 // By copying, distributing or modifying the Program (or any work
3 // * License and Disclaimer << 5 // based on the Program) you indicate your acceptance of this statement,
4 // * << 6 // and all its terms.
5 // * The Geant4 software is copyright of th << 7 //
6 // * the Geant4 Collaboration. It is provided << 8 // $Id: G4IntersectionSolid.cc,v 1.11 2001/03/19 10:30:18 gcosmo Exp $
7 // * conditions of the Geant4 Software License << 9 // GEANT4 tag $Name: geant4-03-01 $
8 // * LICENSE and available at http://cern.ch/ <<
9 // * include a list of copyright holders. <<
10 // * <<
11 // * Neither the authors of this software syst <<
12 // * institutes,nor the agencies providing fin <<
13 // * work make any representation or warran <<
14 // * regarding this software system or assum <<
15 // * use. Please see the license in the file <<
16 // * for the full disclaimer and the limitatio <<
17 // * <<
18 // * This code implementation is the result <<
19 // * technical work of the GEANT4 collaboratio <<
20 // * By using, copying, modifying or distri <<
21 // * any work based on the software) you ag <<
22 // * use in resulting scientific publicati <<
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* <<
25 // 10 //
26 // Implementation of methods for the class G4I 11 // Implementation of methods for the class G4IntersectionSolid
27 // 12 //
28 // 12.09.98 V.Grichine: first implementation << 13 // History:
29 // ------------------------------------------- << 14 //
30 << 15 // 12.09.98 V.Grichine
31 #include <sstream> << 16 // 29.07.99 V.Grichine, modifications in DistanceToIn(p,v)
>> 17 // 16.03.01 V.Grichine, modifications in CalculateExtent() and Inside()
>> 18 // based on D.Williams proposal
32 19
33 #include "G4IntersectionSolid.hh" 20 #include "G4IntersectionSolid.hh"
>> 21 // #include "G4DisplacedSolid.hh"
>> 22
>> 23 #include "G4RotationMatrix.hh"
>> 24 #include "G4ThreeVector.hh"
>> 25 #include "G4Transform3D.hh"
>> 26 #include "G4AffineTransform.hh"
34 27
35 #include "G4SystemOfUnits.hh" <<
36 #include "G4VoxelLimits.hh" 28 #include "G4VoxelLimits.hh"
37 #include "G4VPVParameterisation.hh" 29 #include "G4VPVParameterisation.hh"
38 30
39 #include "G4VGraphicsScene.hh" 31 #include "G4VGraphicsScene.hh"
40 #include "G4Polyhedron.hh" 32 #include "G4Polyhedron.hh"
41 #include "G4PolyhedronArbitrary.hh" << 33 #include "G4NURBS.hh"
42 #include "HepPolyhedronProcessor.h" << 34 #include "G4NURBSbox.hh"
43 35
44 ////////////////////////////////////////////// << 36 /////////////////////////////////////////////////////////////////////
45 // 37 //
46 // Transfer all data members to G4BooleanSolid 38 // Transfer all data members to G4BooleanSolid which is responsible
47 // for them. pName will be in turn sent to G4V 39 // for them. pName will be in turn sent to G4VSolid
48 // 40 //
49 41
50 G4IntersectionSolid::G4IntersectionSolid( cons << 42 G4IntersectionSolid:: G4IntersectionSolid( const G4String& pName,
51 << 43 G4VSolid* pSolidA ,
52 << 44 G4VSolid* pSolidB ):
53 : G4BooleanSolid(pName,pSolidA,pSolidB) << 45 G4BooleanSolid(pName,pSolidA,pSolidB)
54 { 46 {
>> 47 ;
55 } 48 }
56 49
57 ////////////////////////////////////////////// << 50
>> 51 ///////////////////////////////////////////////////////////////////
58 // 52 //
59 53
60 G4IntersectionSolid::G4IntersectionSolid( cons << 54 G4IntersectionSolid::
61 << 55 G4IntersectionSolid( const G4String& pName,
62 << 56 G4VSolid* pSolidA ,
63 << 57 G4VSolid* pSolidB,
64 cons << 58 G4RotationMatrix* rotMatrix,
65 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMa << 59 const G4ThreeVector& transVector ):
>> 60 G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
66 { 61 {
>> 62 ;
67 } 63 }
68 64
69 ////////////////////////////////////////////// << 65 //////////////////////////////////////////////////////////////////
70 // 66 //
71 // 67 //
72 68
73 G4IntersectionSolid::G4IntersectionSolid( cons << 69 G4IntersectionSolid::
74 << 70 G4IntersectionSolid( const G4String& pName,
75 << 71 G4VSolid* pSolidA ,
76 cons << 72 G4VSolid* pSolidB ,
77 : G4BooleanSolid(pName,pSolidA,pSolidB,trans << 73 const G4Transform3D& transform ):
>> 74 G4BooleanSolid(pName,pSolidA,pSolidB,transform)
78 { 75 {
>> 76 ;
79 } 77 }
80 78
81 ////////////////////////////////////////////// <<
82 // <<
83 // Fake default constructor - sets only member <<
84 // for usage restri <<
85 79
86 G4IntersectionSolid::G4IntersectionSolid( __vo << 80 G4IntersectionSolid::~G4IntersectionSolid()
87 : G4BooleanSolid(a) <<
88 { 81 {
>> 82 ;
89 } 83 }
90 84
91 ////////////////////////////////////////////// << 85 ///////////////////////////////////////////////////////////////
92 // <<
93 // <<
94 <<
95 G4IntersectionSolid::~G4IntersectionSolid() = <<
96 <<
97 ////////////////////////////////////////////// <<
98 // 86 //
99 // Copy constructor <<
100 <<
101 G4IntersectionSolid::G4IntersectionSolid(const <<
102 <<
103 ////////////////////////////////////////////// <<
104 // <<
105 // Assignment operator <<
106 <<
107 G4IntersectionSolid& <<
108 G4IntersectionSolid::operator = (const G4Inter <<
109 { <<
110 // Check assignment to self <<
111 // <<
112 if (this == &rhs) { return *this; } <<
113 <<
114 // Copy base class data <<
115 // <<
116 G4BooleanSolid::operator=(rhs); <<
117 <<
118 return *this; <<
119 } <<
120 <<
121 ////////////////////////////////////////////// <<
122 // 87 //
123 // Get bounding box <<
124 <<
125 void <<
126 G4IntersectionSolid::BoundingLimits(G4ThreeVec <<
127 G4ThreeVec <<
128 { <<
129 G4ThreeVector minA,maxA, minB,maxB; <<
130 fPtrSolidA->BoundingLimits(minA,maxA); <<
131 fPtrSolidB->BoundingLimits(minB,maxB); <<
132 <<
133 pMin.set(std::max(minA.x(),minB.x()), <<
134 std::max(minA.y(),minB.y()), <<
135 std::max(minA.z(),minB.z())); <<
136 <<
137 pMax.set(std::min(maxA.x(),maxB.x()), <<
138 std::min(maxA.y(),maxB.y()), <<
139 std::min(maxA.z(),maxB.z())); <<
140 <<
141 // Check correctness of the bounding box <<
142 // <<
143 if (pMin.x() >= pMax.x() || pMin.y() >= pMax <<
144 { <<
145 std::ostringstream message; <<
146 message << "Bad bounding box (min >= max) <<
147 << GetName() << " !" <<
148 << "\npMin = " << pMin <<
149 << "\npMax = " << pMax; <<
150 G4Exception("G4IntersectionSolid::Bounding <<
151 JustWarning, message); <<
152 DumpInfo(); <<
153 } <<
154 } <<
155 <<
156 ////////////////////////////////////////////// <<
157 // <<
158 // Calculate extent under transform and specif <<
159 88
160 G4bool 89 G4bool
161 G4IntersectionSolid::CalculateExtent(const EAx 90 G4IntersectionSolid::CalculateExtent(const EAxis pAxis,
162 const G4V << 91 const G4VoxelLimits& pVoxelLimit,
163 const G4A << 92 const G4AffineTransform& pTransform,
164 G4d << 93 G4double& pMin, G4double& pMax) const
165 G4d << 94 {
166 { << 95 G4bool retA, retB, out ;
167 G4bool retA, retB, out; << 96 G4double minA, minB, maxA, maxB ;
168 G4double minA, minB, maxA, maxB; << 97
169 << 98 retA= fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit, pTransform, minA, maxA);
170 retA = fPtrSolidA << 99 retB= fPtrSolidB->CalculateExtent( pAxis, pVoxelLimit, pTransform, minB, maxB);
171 ->CalculateExtent( pAxis, pVoxelLimi << 100 if(retA && retB)
172 retB = fPtrSolidB <<
173 ->CalculateExtent( pAxis, pVoxelLimi <<
174 <<
175 if( retA && retB ) <<
176 { <<
177 pMin = std::max( minA, minB ); <<
178 pMax = std::min( maxA, maxB ); <<
179 out = (pMax > pMin); // true; <<
180 } <<
181 else <<
182 { 101 {
183 out = false; << 102 pMin = G4std::max( minA, minB ) ;
>> 103 pMax = G4std::min( maxA, maxB ) ;
>> 104 out = true ;
184 } 105 }
>> 106 else out = false ;
185 107
186 return out; // It exists in this slice only << 108 return out ; // It exists in this slice only if both exist in it.
187 } 109 }
188 110
189 ////////////////////////////////////////////// << 111 /////////////////////////////////////////////////////
190 // 112 //
191 // Touching ? Empty intersection ? 113 // Touching ? Empty intersection ?
192 114
193 EInside G4IntersectionSolid::Inside(const G4Th 115 EInside G4IntersectionSolid::Inside(const G4ThreeVector& p) const
194 { 116 {
195 EInside positionA = fPtrSolidA->Inside(p); << 117 EInside positionA = fPtrSolidA->Inside(p) ;
196 if(positionA == kOutside) return positionA; <<
197 118
198 EInside positionB = fPtrSolidB->Inside(p); << 119 if( positionA == kOutside ) return kOutside ;
199 if(positionA == kInside) return positionB; <<
200 120
201 if(positionB == kOutside) return positionB; << 121 EInside positionB = fPtrSolidB->Inside(p) ;
202 return kSurface; << 122
>> 123 if(positionA == kInside && positionB == kInside)
>> 124 {
>> 125 return kInside ;
>> 126 }
>> 127 else
>> 128 {
>> 129 if((positionA == kInside && positionB == kSurface) ||
>> 130 (positionB == kInside && positionA == kSurface) ||
>> 131 (positionA == kSurface && positionB == kSurface) )
>> 132 {
>> 133 return kSurface ;
>> 134 }
>> 135 else
>> 136 {
>> 137 return kOutside ;
>> 138 }
>> 139 }
203 } 140 }
204 141
205 ////////////////////////////////////////////// << 142 //////////////////////////////////////////////////////////////
206 // 143 //
207 144
208 G4ThreeVector 145 G4ThreeVector
209 G4IntersectionSolid::SurfaceNormal( const G4Th 146 G4IntersectionSolid::SurfaceNormal( const G4ThreeVector& p ) const
210 { 147 {
211 G4ThreeVector normal; 148 G4ThreeVector normal;
212 EInside insideA, insideB; << 149 G4bool insideA, insideB;
213 150
214 insideA = fPtrSolidA->Inside(p); << 151 insideA= fPtrSolidA->Inside(p);
215 insideB = fPtrSolidB->Inside(p); << 152 insideB= fPtrSolidB->Inside(p);
216 153
217 #ifdef G4BOOLDEBUG << 154 // if( Inside(p) == kOutside )
218 if( (insideA == kOutside) || (insideB == kOu 155 if( (insideA == kOutside) || (insideB == kOutside) )
219 { 156 {
220 G4cout << "WARNING - Invalid call in " << 157 G4Exception("Invalid call in G4IntersectionSolid::SurfaceNormal(p), point p is outside") ;
221 << "G4IntersectionSolid::SurfaceNor <<
222 << " Point p is outside !" << G4en <<
223 G4cout << " p = " << p << G4endl; <<
224 G4cerr << "WARNING - Invalid call in " <<
225 << "G4IntersectionSolid::SurfaceNor <<
226 << " Point p is outside !" << G4en <<
227 G4cerr << " p = " << p << G4endl; <<
228 } 158 }
229 #endif << 159
>> 160 // OLD: if(fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToOut(p) )
230 161
231 // On the surface of both is difficult ... t 162 // On the surface of both is difficult ... treat it like on A now!
232 // 163 //
>> 164 // if( (insideA == kSurface) && (insideB == kSurface) )
>> 165 // normal= fPtrSolidA->SurfaceNormal(p) ;
>> 166 // else
233 if( insideA == kSurface ) 167 if( insideA == kSurface )
234 { <<
235 normal = fPtrSolidA->SurfaceNormal(p) ; <<
236 } <<
237 else if( insideB == kSurface ) <<
238 { <<
239 normal = fPtrSolidB->SurfaceNormal(p) ; <<
240 } <<
241 else // We are on neither surface, so we sh <<
242 { <<
243 if(fPtrSolidA->DistanceToOut(p) <= fPtrSol <<
244 { <<
245 normal= fPtrSolidA->SurfaceNormal(p) ; <<
246 } <<
247 else <<
248 { 168 {
249 normal= fPtrSolidB->SurfaceNormal(p) ; << 169 normal= fPtrSolidA->SurfaceNormal(p) ;
250 } 170 }
251 #ifdef G4BOOLDEBUG << 171 else if( insideB == kSurface )
252 G4cout << "WARNING - Invalid call in " << 172 {
253 << "G4IntersectionSolid::SurfaceNor << 173 normal= fPtrSolidB->SurfaceNormal(p) ;
254 << " Point p is out of surface !" << 174 }
255 G4cout << " p = " << p << G4endl; << 175 // We are on neither surface, so we should generate an exception
256 G4cerr << "WARNING - Invalid call in " << 176 else
257 << "G4IntersectionSolid::SurfaceNor << 177 {
258 << " Point p is out of surface !" << 178 G4Exception("Invalid call in G4IntersectionSolid::SurfaceNormal(p), point p is not on the surface of the volume.") ;
259 G4cerr << " p = " << p << G4endl; << 179
260 #endif << 180 // Or else
>> 181 if(fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToOut(p) )
>> 182 normal= fPtrSolidA->SurfaceNormal(p) ;
>> 183 else
>> 184 normal= fPtrSolidB->SurfaceNormal(p) ;
261 } 185 }
262 186
263 return normal; 187 return normal;
264 } 188 }
265 189
266 ////////////////////////////////////////////// << 190 /////////////////////////////////////////////////////////////
267 // 191 //
268 // The same algorithm as in DistanceToIn(p) 192 // The same algorithm as in DistanceToIn(p)
269 193
270 G4double 194 G4double
271 G4IntersectionSolid::DistanceToIn( const G4Thr 195 G4IntersectionSolid::DistanceToIn( const G4ThreeVector& p,
272 const G4Thr 196 const G4ThreeVector& v ) const
273 { 197 {
274 G4double dist = 0.0; << 198 G4double dist = 0.0, disTmp = 0.0 ;
275 if( Inside(p) == kInside ) 199 if( Inside(p) == kInside )
276 { 200 {
277 #ifdef G4BOOLDEBUG << 201 G4Exception("Invalid call in G4IntersectionSolid::DistanceToIn(p,v), point p is inside") ;
278 G4cout << "WARNING - Invalid call in " << 202 }
279 << "G4IntersectionSolid::DistanceTo << 203 else
280 << " Point p is inside !" << G4end << 204 {
281 G4cout << " p = " << p << G4endl; << 205 if( fPtrSolidA->Inside(p) != kOutside )
282 G4cout << " v = " << v << G4endl; <<
283 G4cerr << "WARNING - Invalid call in " <<
284 << "G4IntersectionSolid::DistanceTo <<
285 << " Point p is inside !" << G4end <<
286 G4cerr << " p = " << p << G4endl; <<
287 G4cerr << " v = " << v << G4endl; <<
288 #endif <<
289 } <<
290 else // if( Inside(p) == kSurface ) <<
291 { <<
292 EInside wA = fPtrSolidA->Inside(p); <<
293 EInside wB = fPtrSolidB->Inside(p); <<
294 <<
295 G4ThreeVector pA = p, pB = p; <<
296 G4double dA = 0., dA1=0., dA2=0.; <<
297 G4double dB = 0., dB1=0., dB2=0.; <<
298 G4bool doA = true, doB = true; <<
299 <<
300 static const std::size_t max_trials=10000; <<
301 for (std::size_t trial=0; trial<max_trials <<
302 { 206 {
303 if(doA) << 207 do
304 { 208 {
305 // find next valid range for A << 209 disTmp = fPtrSolidB->DistanceToIn(p+dist*v,v) ;
306 << 210
307 dA1 = 0.; << 211 if( disTmp != kInfinity )
308 <<
309 if( wA != kInside ) <<
310 { 212 {
311 dA1 = fPtrSolidA->DistanceToIn(pA, v << 213 dist += disTmp ;
312 << 214 if(Inside(p+dist*v) == kOutside )
313 if( dA1 == kInfinity ) return kInf << 215 {
314 << 216 disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
315 pA += dA1*v; << 217 if( disTmp != kInfinity )
>> 218 {
>> 219 dist += disTmp ;
>> 220 }
>> 221 else
>> 222 {
>> 223 return kInfinity ;
>> 224 }
>> 225 }
>> 226 else
>> 227 {
>> 228 break ;
>> 229 }
316 } 230 }
317 dA2 = dA1 + fPtrSolidA->DistanceToOut( << 231 else
>> 232 {
>> 233 return kInfinity ;
>> 234 }
318 } 235 }
319 dA1 += dA; << 236 while( Inside(p+dist*v) == kOutside ) ;
320 dA2 += dA; << 237 // while( fPtrSolidB->Inside(p) != kOutside ) ;
321 << 238 }
322 if(doB) << 239 else if( fPtrSolidB->Inside(p) != kOutside )
>> 240 {
>> 241 do
323 { 242 {
324 // find next valid range for B << 243 disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
325 << 244
326 dB1 = 0.; << 245 if( disTmp != kInfinity )
327 if(wB != kInside) <<
328 { 246 {
329 dB1 = fPtrSolidB->DistanceToIn(pB, v << 247 dist += disTmp ;
330 << 248 if(Inside(p+dist*v) == kOutside )
331 if(dB1 == kInfinity) return kInfin << 249 {
332 << 250 disTmp = fPtrSolidB->DistanceToIn(p+dist*v,v) ;
333 pB += dB1*v; << 251 if( disTmp != kInfinity )
>> 252 {
>> 253 dist += disTmp ;
>> 254 }
>> 255 else
>> 256 {
>> 257 return kInfinity ;
>> 258 }
>> 259 }
>> 260 else
>> 261 {
>> 262 break ;
>> 263 }
334 } 264 }
335 dB2 = dB1 + fPtrSolidB->DistanceToOut( << 265 else
336 } << 266 {
337 dB1 += dB; << 267 return kInfinity ;
338 dB2 += dB; << 268 }
339 <<
340 // check if they overlap <<
341 <<
342 if( dA1 < dB1 ) <<
343 { <<
344 if( dB1 < dA2 ) return dB1; <<
345 <<
346 dA = dA2; <<
347 pA = p + dA*v; // continue from her <<
348 wA = kSurface; <<
349 doA = true; <<
350 doB = false; <<
351 } 269 }
352 else << 270 while( Inside(p+dist*v) == kOutside ) ;
>> 271 }
>> 272 else
>> 273 {
>> 274 do
353 { 275 {
354 if( dA1 < dB2 ) return dA1; << 276 disTmp = fPtrSolidB->DistanceToIn(p+dist*v,v) ;
355 << 277
356 dB = dB2; << 278 if( disTmp != kInfinity )
357 pB = p + dB*v; // continue from her << 279 {
358 wB = kSurface; << 280 dist += disTmp ;
359 doB = true; << 281 if(Inside(p+dist*v) == kOutside )
360 doA = false; << 282 {
>> 283 disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
>> 284 if( disTmp != kInfinity )
>> 285 {
>> 286 dist += disTmp ;
>> 287 }
>> 288 else
>> 289 {
>> 290 return kInfinity ;
>> 291 }
>> 292 }
>> 293 else
>> 294 {
>> 295 break ;
>> 296 }
>> 297 }
>> 298 else
>> 299 {
>> 300 return kInfinity ;
>> 301 }
361 } 302 }
>> 303 while( Inside(p+dist*v) == kOutside ) ;
362 } 304 }
363 } 305 }
364 #ifdef G4BOOLDEBUG <<
365 G4Exception("G4IntersectionSolid::DistanceTo <<
366 "GeomSolids0001", JustWarning, <<
367 "Reached maximum number of itera <<
368 #endif <<
369 return dist ; 306 return dist ;
370 } 307 }
371 308
372 ////////////////////////////////////////////// << 309 ////////////////////////////////////////////////////////
373 // 310 //
374 // Approximate nearest distance from the point 311 // Approximate nearest distance from the point p to the intersection of
375 // two solids 312 // two solids
376 313
377 G4double 314 G4double
378 G4IntersectionSolid::DistanceToIn( const G4Thr 315 G4IntersectionSolid::DistanceToIn( const G4ThreeVector& p) const
379 { 316 {
380 #ifdef G4BOOLDEBUG <<
381 if( Inside(p) == kInside ) 317 if( Inside(p) == kInside )
382 { 318 {
383 G4cout << "WARNING - Invalid call in " << 319 G4Exception("Invalid call in G4IntersectionSolid::DistanceToIn(p), point p is inside") ;
384 << "G4IntersectionSolid::DistanceTo <<
385 << " Point p is inside !" << G4end <<
386 G4cout << " p = " << p << G4endl; <<
387 G4cerr << "WARNING - Invalid call in " <<
388 << "G4IntersectionSolid::DistanceTo <<
389 << " Point p is inside !" << G4end <<
390 G4cerr << " p = " << p << G4endl; <<
391 } 320 }
392 #endif <<
393 EInside sideA = fPtrSolidA->Inside(p) ; 321 EInside sideA = fPtrSolidA->Inside(p) ;
394 EInside sideB = fPtrSolidB->Inside(p) ; 322 EInside sideB = fPtrSolidB->Inside(p) ;
395 G4double dist=0.0 ; << 323 G4double dist ;
396 324
397 if( sideA != kInside && sideB != kOutside ) << 325 if( sideA != kInside && sideB != kOutside )
398 { 326 {
399 dist = fPtrSolidA->DistanceToIn(p) ; 327 dist = fPtrSolidA->DistanceToIn(p) ;
400 } 328 }
401 else 329 else
402 { 330 {
403 if( sideB != kInside && sideA != kOutside << 331 if( sideB != kInside && sideA != kOutside )
404 { 332 {
405 dist = fPtrSolidB->DistanceToIn(p) ; 333 dist = fPtrSolidB->DistanceToIn(p) ;
406 } 334 }
407 else 335 else
408 { 336 {
409 dist = std::min(fPtrSolidA->DistanceToI << 337 dist = G4std::min(fPtrSolidA->DistanceToIn(p),
410 fPtrSolidB->DistanceToI << 338 fPtrSolidB->DistanceToIn(p) ) ;
411 } 339 }
412 } 340 }
413 return dist ; 341 return dist ;
414 } 342 }
415 343
416 ////////////////////////////////////////////// << 344 //////////////////////////////////////////////////////////
417 // 345 //
418 // The same algorithm as DistanceToOut(p) 346 // The same algorithm as DistanceToOut(p)
419 347
420 G4double 348 G4double
421 G4IntersectionSolid::DistanceToOut( const G4Th 349 G4IntersectionSolid::DistanceToOut( const G4ThreeVector& p,
422 const G4Th << 350 const G4ThreeVector& v,
423 const G4bo << 351 const G4bool calcNorm,
424 G4bo << 352 G4bool *validNorm,
425 G4Th << 353 G4ThreeVector *n ) const
426 { 354 {
427 G4bool validNormA, validNormB; 355 G4bool validNormA, validNormB;
428 G4ThreeVector nA, nB; 356 G4ThreeVector nA, nB;
429 357
430 #ifdef G4BOOLDEBUG <<
431 if( Inside(p) == kOutside ) 358 if( Inside(p) == kOutside )
432 { 359 {
433 G4cout << "Position:" << G4endl << G4endl 360 G4cout << "Position:" << G4endl << G4endl;
434 G4cout << "p.x() = " << p.x()/mm << " mm 361 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
435 G4cout << "p.y() = " << p.y()/mm << " mm 362 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
436 G4cout << "p.z() = " << p.z()/mm << " mm 363 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
437 G4cout << "Direction:" << G4endl << G4endl 364 G4cout << "Direction:" << G4endl << G4endl;
438 G4cout << "v.x() = " << v.x() << G4endl; 365 G4cout << "v.x() = " << v.x() << G4endl;
439 G4cout << "v.y() = " << v.y() << G4endl; 366 G4cout << "v.y() = " << v.y() << G4endl;
440 G4cout << "v.z() = " << v.z() << G4endl 367 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
441 G4cout << "WARNING - Invalid call in " << 368 G4Exception("Invalid call in G4IntersectionSolid::DistanceToOut(p,v), point p is outside") ;
442 << "G4IntersectionSolid::DistanceTo <<
443 << " Point p is outside !" << G4en <<
444 G4cout << " p = " << p << G4endl; <<
445 G4cout << " v = " << v << G4endl; <<
446 G4cerr << "WARNING - Invalid call in " <<
447 << "G4IntersectionSolid::DistanceTo <<
448 << " Point p is outside !" << G4en <<
449 G4cerr << " p = " << p << G4endl; <<
450 G4cerr << " v = " << v << G4endl; <<
451 } 369 }
452 #endif <<
453 G4double distA = fPtrSolidA->DistanceToOut(p 370 G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,&validNormA,&nA) ;
454 G4double distB = fPtrSolidB->DistanceToOut(p 371 G4double distB = fPtrSolidB->DistanceToOut(p,v,calcNorm,&validNormB,&nB) ;
455 372
456 G4double dist = std::min(distA,distB) ; << 373 G4double dist = G4std::min(distA,distB) ;
457 374
458 if( calcNorm ) << 375 if( calcNorm ){
459 { << 376 if (distA < distB ) {
460 if ( distA < distB ) << 377 *validNorm = validNormA;
461 { << 378 *n = nA;
462 *validNorm = validNormA; << 379 }else{
463 *n = nA; << 380 *validNorm = validNormB;
464 } << 381 *n = nB;
465 else << 382 }
466 { <<
467 *validNorm = validNormB; <<
468 *n = nB; <<
469 } <<
470 } 383 }
471 384
472 return dist ; 385 return dist ;
>> 386 // return G4std::min(fPtrSolidA->DistanceToOut(p,v,calcNorm,validNorm,n),
>> 387 // fPtrSolidB->DistanceToOut(p,v,calcNorm,validNorm,n) ) ;
473 } 388 }
474 389
475 ////////////////////////////////////////////// << 390 //////////////////////////////////////////////////////////////
476 // 391 //
477 // Inverted algorithm of DistanceToIn(p) 392 // Inverted algorithm of DistanceToIn(p)
478 393
479 G4double 394 G4double
480 G4IntersectionSolid::DistanceToOut( const G4Th 395 G4IntersectionSolid::DistanceToOut( const G4ThreeVector& p ) const
481 { 396 {
482 #ifdef G4BOOLDEBUG <<
483 if( Inside(p) == kOutside ) 397 if( Inside(p) == kOutside )
484 { 398 {
485 G4cout << "WARNING - Invalid call in " << 399 G4Exception("Invalid call in G4IntersectionSolid::DistanceToOut(p), point p is outside") ;
486 << "G4IntersectionSolid::DistanceTo <<
487 << " Point p is outside !" << G4en <<
488 G4cout << " p = " << p << G4endl; <<
489 G4cerr << "WARNING - Invalid call in " <<
490 << "G4IntersectionSolid::DistanceTo <<
491 << " Point p is outside !" << G4en <<
492 G4cerr << " p = " << p << G4endl; <<
493 } 400 }
494 #endif <<
495 401
496 return std::min(fPtrSolidA->DistanceToOut(p) << 402
497 fPtrSolidB->DistanceToOut(p) << 403 return G4std::min(fPtrSolidA->DistanceToOut(p),
>> 404 fPtrSolidB->DistanceToOut(p) ) ;
498 405
499 } 406 }
500 407
501 ////////////////////////////////////////////// << 408 //////////////////////////////////////////////////////////////
>> 409 //
502 // 410 //
503 // ComputeDimensions <<
504 411
505 void 412 void
506 G4IntersectionSolid::ComputeDimensions( G4VPVP << 413 G4IntersectionSolid::ComputeDimensions( G4VPVParameterisation* p,
507 const << 414 const G4int n,
508 const << 415 const G4VPhysicalVolume* pRep )
509 { 416 {
>> 417 return ;
510 } 418 }
511 419
512 ////////////////////////////////////////////// << 420 /////////////////////////////////////////////////
513 // 421 //
514 // GetEntityType << 422 //
515 423
516 G4GeometryType G4IntersectionSolid::GetEntityT 424 G4GeometryType G4IntersectionSolid::GetEntityType() const
517 { 425 {
518 return {"G4IntersectionSolid"}; << 426 return G4String("G4IntersectionSolid");
519 } 427 }
520 428
521 ////////////////////////////////////////////// << 429 /////////////////////////////////////////////////
522 // 430 //
523 // Make a clone of the object << 431 //
524 432
525 G4VSolid* G4IntersectionSolid::Clone() const << 433 void
>> 434 G4IntersectionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
526 { 435 {
527 return new G4IntersectionSolid(*this); << 436 scene.AddThis (*this);
528 } 437 }
529 438
530 ////////////////////////////////////////////// << 439 ////////////////////////////////////////////////////
>> 440 //
531 // 441 //
532 // DescribeYourselfTo <<
533 442
534 void << 443 G4Polyhedron*
535 G4IntersectionSolid::DescribeYourselfTo ( G4VG << 444 G4IntersectionSolid::CreatePolyhedron () const
536 { 445 {
537 scene.AddSolid (*this); << 446 G4Polyhedron* pA = fPtrSolidA->CreatePolyhedron();
>> 447 G4Polyhedron* pB = fPtrSolidB->CreatePolyhedron();
>> 448 G4Polyhedron* resultant = new G4Polyhedron (pA->intersect(*pB));
>> 449 delete pB;
>> 450 delete pA;
>> 451 return resultant;
538 } 452 }
539 453
540 ////////////////////////////////////////////// << 454 /////////////////////////////////////////////////////////
>> 455 //
541 // 456 //
542 // CreatePolyhedron <<
543 457
544 G4Polyhedron* << 458 G4NURBS*
545 G4IntersectionSolid::CreatePolyhedron () const << 459 G4IntersectionSolid::CreateNURBS () const
546 { 460 {
547 if (fExternalBoolProcessor == nullptr) << 461 // Take into account boolean operation - see CreatePolyhedron.
548 { << 462 // return new G4NURBSbox (1.0, 1.0, 1.0);
549 HepPolyhedronProcessor processor; << 463 return 0;
550 // Stack components and components of comp <<
551 // See G4BooleanSolid::StackPolyhedron <<
552 G4Polyhedron* top = StackPolyhedron(proces <<
553 auto result = new G4Polyhedron(*top); <<
554 if (processor.execute(*result)) <<
555 { <<
556 return result; <<
557 } <<
558 else <<
559 { <<
560 return nullptr; <<
561 } <<
562 } <<
563 else <<
564 { <<
565 return fExternalBoolProcessor->Process(thi <<
566 } <<
567 } 464 }
568 465