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