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