| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
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 //
26 // Implementation of methods for the class G4I
27 //
28 // 12.09.98 V.Grichine: first implementation
29 // -------------------------------------------
30
31 #include <sstream>
32
33 #include "G4IntersectionSolid.hh"
34
35 #include "G4SystemOfUnits.hh"
36 #include "G4VoxelLimits.hh"
37 #include "G4VPVParameterisation.hh"
38
39 #include "G4VGraphicsScene.hh"
40 #include "G4Polyhedron.hh"
41 #include "G4PolyhedronArbitrary.hh"
42 #include "HepPolyhedronProcessor.h"
43
44 //////////////////////////////////////////////
45 //
46 // Transfer all data members to G4BooleanSolid
47 // for them. pName will be in turn sent to G4V
48 //
49
50 G4IntersectionSolid::G4IntersectionSolid( cons
51
52
53 : G4BooleanSolid(pName,pSolidA,pSolidB)
54 {
55 }
56
57 //////////////////////////////////////////////
58 //
59
60 G4IntersectionSolid::G4IntersectionSolid( cons
61
62
63
64 cons
65 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMa
66 {
67 }
68
69 //////////////////////////////////////////////
70 //
71 //
72
73 G4IntersectionSolid::G4IntersectionSolid( cons
74
75
76 cons
77 : G4BooleanSolid(pName,pSolidA,pSolidB,trans
78 {
79 }
80
81 //////////////////////////////////////////////
82 //
83 // Fake default constructor - sets only member
84 // for usage restri
85
86 G4IntersectionSolid::G4IntersectionSolid( __vo
87 : G4BooleanSolid(a)
88 {
89 }
90
91 //////////////////////////////////////////////
92 //
93 //
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
160 G4bool
161 G4IntersectionSolid::CalculateExtent(const EAx
162 const G4V
163 const G4A
164 G4d
165 G4d
166 {
167 G4bool retA, retB, out;
168 G4double minA, minB, maxA, maxB;
169
170 retA = fPtrSolidA
171 ->CalculateExtent( pAxis, pVoxelLimi
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 {
183 out = false;
184 }
185
186 return out; // It exists in this slice only
187 }
188
189 //////////////////////////////////////////////
190 //
191 // Touching ? Empty intersection ?
192
193 EInside G4IntersectionSolid::Inside(const G4Th
194 {
195 EInside positionA = fPtrSolidA->Inside(p);
196 if(positionA == kOutside) return positionA;
197
198 EInside positionB = fPtrSolidB->Inside(p);
199 if(positionA == kInside) return positionB;
200
201 if(positionB == kOutside) return positionB;
202 return kSurface;
203 }
204
205 //////////////////////////////////////////////
206 //
207
208 G4ThreeVector
209 G4IntersectionSolid::SurfaceNormal( const G4Th
210 {
211 G4ThreeVector normal;
212 EInside insideA, insideB;
213
214 insideA = fPtrSolidA->Inside(p);
215 insideB = fPtrSolidB->Inside(p);
216
217 #ifdef G4BOOLDEBUG
218 if( (insideA == kOutside) || (insideB == kOu
219 {
220 G4cout << "WARNING - Invalid call in "
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 }
229 #endif
230
231 // On the surface of both is difficult ... t
232 //
233 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 {
249 normal= fPtrSolidB->SurfaceNormal(p) ;
250 }
251 #ifdef G4BOOLDEBUG
252 G4cout << "WARNING - Invalid call in "
253 << "G4IntersectionSolid::SurfaceNor
254 << " Point p is out of surface !"
255 G4cout << " p = " << p << G4endl;
256 G4cerr << "WARNING - Invalid call in "
257 << "G4IntersectionSolid::SurfaceNor
258 << " Point p is out of surface !"
259 G4cerr << " p = " << p << G4endl;
260 #endif
261 }
262
263 return normal;
264 }
265
266 //////////////////////////////////////////////
267 //
268 // The same algorithm as in DistanceToIn(p)
269
270 G4double
271 G4IntersectionSolid::DistanceToIn( const G4Thr
272 const G4Thr
273 {
274 G4double dist = 0.0;
275 if( Inside(p) == kInside )
276 {
277 #ifdef G4BOOLDEBUG
278 G4cout << "WARNING - Invalid call in "
279 << "G4IntersectionSolid::DistanceTo
280 << " Point p is inside !" << G4end
281 G4cout << " p = " << p << G4endl;
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 {
303 if(doA)
304 {
305 // find next valid range for A
306
307 dA1 = 0.;
308
309 if( wA != kInside )
310 {
311 dA1 = fPtrSolidA->DistanceToIn(pA, v
312
313 if( dA1 == kInfinity ) return kInf
314
315 pA += dA1*v;
316 }
317 dA2 = dA1 + fPtrSolidA->DistanceToOut(
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 {
329 dB1 = fPtrSolidB->DistanceToIn(pB, v
330
331 if(dB1 == kInfinity) return kInfin
332
333 pB += dB1*v;
334 }
335 dB2 = dB1 + fPtrSolidB->DistanceToOut(
336 }
337 dB1 += dB;
338 dB2 += dB;
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 }
352 else
353 {
354 if( dA1 < dB2 ) return dA1;
355
356 dB = dB2;
357 pB = p + dB*v; // continue from her
358 wB = kSurface;
359 doB = true;
360 doA = false;
361 }
362 }
363 }
364 #ifdef G4BOOLDEBUG
365 G4Exception("G4IntersectionSolid::DistanceTo
366 "GeomSolids0001", JustWarning,
367 "Reached maximum number of itera
368 #endif
369 return dist ;
370 }
371
372 //////////////////////////////////////////////
373 //
374 // Approximate nearest distance from the point
375 // two solids
376
377 G4double
378 G4IntersectionSolid::DistanceToIn( const G4Thr
379 {
380 #ifdef G4BOOLDEBUG
381 if( Inside(p) == kInside )
382 {
383 G4cout << "WARNING - Invalid call in "
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 }
392 #endif
393 EInside sideA = fPtrSolidA->Inside(p) ;
394 EInside sideB = fPtrSolidB->Inside(p) ;
395 G4double dist=0.0 ;
396
397 if( sideA != kInside && sideB != kOutside )
398 {
399 dist = fPtrSolidA->DistanceToIn(p) ;
400 }
401 else
402 {
403 if( sideB != kInside && sideA != kOutside
404 {
405 dist = fPtrSolidB->DistanceToIn(p) ;
406 }
407 else
408 {
409 dist = std::min(fPtrSolidA->DistanceToI
410 fPtrSolidB->DistanceToI
411 }
412 }
413 return dist ;
414 }
415
416 //////////////////////////////////////////////
417 //
418 // The same algorithm as DistanceToOut(p)
419
420 G4double
421 G4IntersectionSolid::DistanceToOut( const G4Th
422 const G4Th
423 const G4bo
424 G4bo
425 G4Th
426 {
427 G4bool validNormA, validNormB;
428 G4ThreeVector nA, nB;
429
430 #ifdef G4BOOLDEBUG
431 if( Inside(p) == kOutside )
432 {
433 G4cout << "Position:" << G4endl << G4endl
434 G4cout << "p.x() = " << p.x()/mm << " mm
435 G4cout << "p.y() = " << p.y()/mm << " mm
436 G4cout << "p.z() = " << p.z()/mm << " mm
437 G4cout << "Direction:" << G4endl << G4endl
438 G4cout << "v.x() = " << v.x() << G4endl;
439 G4cout << "v.y() = " << v.y() << G4endl;
440 G4cout << "v.z() = " << v.z() << G4endl
441 G4cout << "WARNING - Invalid call in "
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 }
452 #endif
453 G4double distA = fPtrSolidA->DistanceToOut(p
454 G4double distB = fPtrSolidB->DistanceToOut(p
455
456 G4double dist = std::min(distA,distB) ;
457
458 if( calcNorm )
459 {
460 if ( distA < distB )
461 {
462 *validNorm = validNormA;
463 *n = nA;
464 }
465 else
466 {
467 *validNorm = validNormB;
468 *n = nB;
469 }
470 }
471
472 return dist ;
473 }
474
475 //////////////////////////////////////////////
476 //
477 // Inverted algorithm of DistanceToIn(p)
478
479 G4double
480 G4IntersectionSolid::DistanceToOut( const G4Th
481 {
482 #ifdef G4BOOLDEBUG
483 if( Inside(p) == kOutside )
484 {
485 G4cout << "WARNING - Invalid call in "
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 }
494 #endif
495
496 return std::min(fPtrSolidA->DistanceToOut(p)
497 fPtrSolidB->DistanceToOut(p)
498
499 }
500
501 //////////////////////////////////////////////
502 //
503 // ComputeDimensions
504
505 void
506 G4IntersectionSolid::ComputeDimensions( G4VPVP
507 const
508 const
509 {
510 }
511
512 //////////////////////////////////////////////
513 //
514 // GetEntityType
515
516 G4GeometryType G4IntersectionSolid::GetEntityT
517 {
518 return {"G4IntersectionSolid"};
519 }
520
521 //////////////////////////////////////////////
522 //
523 // Make a clone of the object
524
525 G4VSolid* G4IntersectionSolid::Clone() const
526 {
527 return new G4IntersectionSolid(*this);
528 }
529
530 //////////////////////////////////////////////
531 //
532 // DescribeYourselfTo
533
534 void
535 G4IntersectionSolid::DescribeYourselfTo ( G4VG
536 {
537 scene.AddSolid (*this);
538 }
539
540 //////////////////////////////////////////////
541 //
542 // CreatePolyhedron
543
544 G4Polyhedron*
545 G4IntersectionSolid::CreatePolyhedron () const
546 {
547 if (fExternalBoolProcessor == nullptr)
548 {
549 HepPolyhedronProcessor processor;
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 }
568