| 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 G4MultiUnion class
27 //
28 // 19.10.12 M.Gayer - Original implementation
29 // 06.04.17 G.Cosmo - Adapted implementation i
30 // -------------------------------------------
31
32 #include <iostream>
33 #include <sstream>
34
35 #include "G4MultiUnion.hh"
36 #include "Randomize.hh"
37 #include "G4GeometryTolerance.hh"
38 #include "G4BoundingEnvelope.hh"
39 #include "G4AffineTransform.hh"
40 #include "G4DisplacedSolid.hh"
41
42 #include "G4VGraphicsScene.hh"
43 #include "G4Polyhedron.hh"
44 #include "G4PolyhedronArbitrary.hh"
45 #include "HepPolyhedronProcessor.h"
46
47 #include "G4BooleanSolid.hh"
48
49 #include "G4AutoLock.hh"
50
51 namespace
52 {
53 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE
54 }
55
56 //____________________________________________
57 G4MultiUnion::G4MultiUnion(const G4String& nam
58 : G4VSolid(name)
59 {
60 SetName(name);
61 fSolids.clear();
62 fTransformObjs.clear();
63 kRadTolerance = G4GeometryTolerance::GetInst
64 }
65
66 //____________________________________________
67 G4MultiUnion::~G4MultiUnion()
68 = default;
69
70 //____________________________________________
71 void G4MultiUnion::AddNode(G4VSolid& solid, co
72 {
73 fSolids.push_back(&solid);
74 fTransformObjs.push_back(trans); // Store a
75 }
76
77 //____________________________________________
78 void G4MultiUnion::AddNode(G4VSolid* solid, co
79 {
80 fSolids.push_back(solid);
81 fTransformObjs.push_back(trans); // Store a
82 }
83
84 //____________________________________________
85 G4VSolid* G4MultiUnion::Clone() const
86 {
87 return new G4MultiUnion(*this);
88 }
89
90 // Copy constructor
91 //____________________________________________
92 G4MultiUnion::G4MultiUnion(const G4MultiUnion&
93 : G4VSolid(rhs), fCubicVolume(rhs.fCubicVolu
94 fSurfaceArea(rhs.fSurfaceArea),
95 kRadTolerance(rhs.kRadTolerance), fAccurat
96 {
97 }
98
99 // Fake default constructor for persistency
100 //____________________________________________
101 G4MultiUnion::G4MultiUnion( __void__& a )
102 : G4VSolid(a)
103 {
104 }
105
106 // Assignment operator
107 //____________________________________________
108 G4MultiUnion& G4MultiUnion::operator = (const
109 {
110 // Check assignment to self
111 //
112 if (this == &rhs)
113 {
114 return *this;
115 }
116
117 // Copy base class data
118 //
119 G4VSolid::operator=(rhs);
120
121 return *this;
122 }
123
124 //____________________________________________
125 G4double G4MultiUnion::GetCubicVolume()
126 {
127 if (fCubicVolume == 0.0)
128 {
129 fCubicVolume = EstimateCubicVolume(1000000
130 }
131 return fCubicVolume;
132 }
133
134 //____________________________________________
135 G4double
136 G4MultiUnion::DistanceToInNoVoxels(const G4Thr
137 const G4Thr
138 {
139 G4ThreeVector direction = aDirection.unit();
140 G4ThreeVector localPoint, localDirection;
141 G4double minDistance = kInfinity;
142
143 std::size_t numNodes = fSolids.size();
144 for (std::size_t i = 0 ; i < numNodes ; ++i)
145 {
146 G4VSolid& solid = *fSolids[i];
147 const G4Transform3D& transform = fTransfor
148
149 localPoint = GetLocalPoint(transform, aPoi
150 localDirection = GetLocalVector(transform,
151
152 G4double distance = solid.DistanceToIn(loc
153 if (minDistance > distance) minDistance =
154 }
155 return minDistance;
156 }
157
158 //____________________________________________
159 G4double G4MultiUnion::DistanceToInCandidates(
160
161
162
163 {
164 std::size_t candidatesCount = candidates.siz
165 G4ThreeVector localPoint, localDirection;
166
167 G4double minDistance = kInfinity;
168 for (std::size_t i = 0 ; i < candidatesCount
169 {
170 G4int candidate = candidates[i];
171 G4VSolid& solid = *fSolids[candidate];
172 const G4Transform3D& transform = fTransfor
173
174 localPoint = GetLocalPoint(transform, aPoi
175 localDirection = GetLocalVector(transform,
176 G4double distance = solid.DistanceToIn(loc
177 if (minDistance > distance) minDistance =
178 bits.SetBitNumber(candidate);
179 if (minDistance == 0) break;
180 }
181 return minDistance;
182 }
183
184 // Algorithm note: we have to look also for al
185 // if the distance is not shorter ... we have
186 // for example for objects which starts in fir
187 // do not collide with direction line, but in
188 // The idea of crossing voxels would be still
189 // because this way we could exclude from the
190 // which were found that obviously are not goo
191 // they would return infinity
192 // But if distance is smaller than the shift t
193 // it immediately
194 //____________________________________________
195 G4double G4MultiUnion::DistanceToIn(const G4Th
196 const G4Th
197 {
198 G4double minDistance = kInfinity;
199 G4ThreeVector direction = aDirection.unit();
200 G4double shift = fVoxels.DistanceToFirst(aPo
201 if (shift == kInfinity) return shift;
202
203 G4ThreeVector currentPoint = aPoint;
204 if (shift != 0.0) currentPoint += direction
205
206 G4SurfBits exclusion(fVoxels.GetBitsPerSlice
207 std::vector<G4int> candidates, curVoxel(3);
208 fVoxels.GetVoxel(curVoxel, currentPoint);
209
210 do
211 {
212 {
213 if (fVoxels.GetCandidatesVoxelArray(curV
214 {
215 G4double distance = DistanceToInCandid
216
217 if (minDistance > distance) minDistanc
218 if (distance < shift) break;
219 }
220 }
221 shift = fVoxels.DistanceToNext(aPoint, dir
222 }
223 while (minDistance > shift);
224
225 return minDistance;
226 }
227
228 //____________________________________________
229 G4double G4MultiUnion::DistanceToOutNoVoxels(c
230 c
231 G
232 {
233 // Computes distance from a point presumably
234 // surface. Ignores first surface if the poi
235 // Early return infinity in case the safety
236 // than the proposed step aPstep.
237 // The normal vector to the crossed surface
238 // is crossed, otherwise aNormal->IsNull() i
239
240 // algorithm:
241 G4ThreeVector direction = aDirection.unit();
242 G4ThreeVector localPoint, localDirection;
243 G4int ignoredSolid = -1;
244 G4double resultDistToOut = 0;
245 G4ThreeVector currentPoint = aPoint;
246
247 auto numNodes = (G4int)fSolids.size();
248 for (auto i = 0; i < numNodes; ++i)
249 {
250 if (i != ignoredSolid)
251 {
252 G4VSolid& solid = *fSolids[i];
253 const G4Transform3D& transform = fTransf
254 localPoint = GetLocalPoint(transform, cu
255 localDirection = GetLocalVector(transfor
256 EInside location = solid.Inside(localPoi
257 if (location != EInside::kOutside)
258 {
259 G4double distance = solid.DistanceToOu
260
261 if (distance < kInfinity)
262 {
263 if (resultDistToOut == kInfinity) re
264 if (distance > 0)
265 {
266 currentPoint = GetGlobalPoint(tran
267 + di
268 resultDistToOut += distance;
269 ignoredSolid = i; // skip the soli
270 i = -1; // force the loop to conti
271 }
272 }
273 }
274 }
275 }
276 return resultDistToOut;
277 }
278
279 //____________________________________________
280 G4double G4MultiUnion::DistanceToOut(const G4T
281 const G4T
282 const G4b
283 G4bool* /
284 G4ThreeVe
285 {
286 return DistanceToOutVoxels(aPoint, aDirectio
287 }
288
289 //____________________________________________
290 G4double G4MultiUnion::DistanceToOutVoxels(con
291 con
292 G4T
293 {
294 // Computes distance from a point presumably
295 // surface. Ignores first surface along each
296 // inside or outside. Early returns zero in
297 // the starting point.
298 // o The proposed step is ignored.
299 // o The normal vector to the crossed surfac
300
301 // In the case the considered point is locat
302 // structure, the treatments are as follows:
303 // - investigation of the candidates fo
304 // - progressive moving of the point to
305 // passed direction
306 // - processing of the normal
307
308 G4ThreeVector direction = aDirection.unit();
309 std::vector<G4int> candidates;
310 G4double distance = 0;
311 std::size_t numNodes = 2*fSolids.size();
312 std::size_t count=0;
313
314 if (fVoxels.GetCandidatesVoxelArray(aPoint,
315 {
316 // For normal case for which we presume th
317 G4ThreeVector localPoint, localDirection,
318 G4ThreeVector currentPoint = aPoint;
319 G4SurfBits exclusion(fVoxels.GetBitsPerSli
320 G4bool notOutside;
321 G4ThreeVector maxNormal;
322
323 do
324 {
325 notOutside = false;
326
327 G4double maxDistance = -kInfinity;
328 G4int maxCandidate = 0;
329 G4ThreeVector maxLocalPoint;
330
331 std::size_t limit = candidates.size();
332 for (std::size_t i = 0 ; i < limit ; ++i
333 {
334 G4int candidate = candidates[i];
335 // ignore the current component (that
336 // numerically the propagated point wi
337
338 G4VSolid& solid = *fSolids[candidate];
339 const G4Transform3D& transform = fTran
340
341 // The coordinates of the point are mo
342 // intrinsic solid local frame:
343 localPoint = GetLocalPoint(transform,
344
345 // DistanceToOut at least for Trd some
346 // even from points that are outside.
347 // must currently be here, otherwise i
348 // But it means it would be slower.
349
350 if (solid.Inside(localPoint) != EInsid
351 {
352 notOutside = true;
353
354 localDirection = GetLocalVector(tran
355
356 // propagate with solid.DistanceToOu
357 G4double shift = solid.DistanceToOut
358
359 if (maxDistance < shift)
360 {
361 maxDistance = shift;
362 maxCandidate = candidate;
363 maxNormal = localNormal;
364 }
365 }
366 }
367
368 if (notOutside)
369 {
370 const G4Transform3D& transform = fTran
371
372 // convert from local normal
373 if (aNormal != nullptr) *aNormal = Get
374
375 distance += maxDistance;
376 currentPoint += maxDistance * directio
377 if(maxDistance == 0.) ++count;
378
379 // the current component will be ignor
380 exclusion.SetBitNumber(maxCandidate);
381 EInside location = InsideWithExclusion
382
383 // perform a Inside
384 // it should be excluded current solid
385 // we have to collect the maximum dist
386 // such "maximum" candidate should be
387 // candidates
388 if (location == EInside::kOutside)
389 {
390 // else return cumulated distances t
391 // components
392 break;
393 }
394 // if inside another component, redo 1
395 // DistanceToOut on top of the previou
396
397 // and fill the candidates for the cor
398 // exiting current component along dir
399 candidates.clear();
400
401 fVoxels.GetCandidatesVoxelArray(curren
402 exclusion.ResetBitNumber(maxCandidate)
403 }
404 }
405 while ((notOutside) && (count < numNodes)
406 }
407
408 return distance;
409 }
410
411 //____________________________________________
412 EInside G4MultiUnion::InsideWithExclusion(cons
413 G4Su
414 {
415 // Classify point location with respect to s
416 // o eInside - inside the solid
417 // o eSurface - close to surface withi
418 // o eOutside - outside the solid
419
420 // Hitherto, it is considered that only para
421 // can be added to the container
422
423 // Implementation using voxelisation techniq
424 // -----------------------------------------
425
426 G4ThreeVector localPoint;
427 EInside location = EInside::kOutside;
428
429 std::vector<G4int> candidates;
430 std::vector<G4MultiUnionSurface> surfaces;
431
432 // TODO: test if it works well and if so mea
433 // TODO: getPointIndex should not be used, i
434 // should be used
435 // TODO: than pass result to GetVoxel furthe
436 // TODO: eventually GetVoxel should be inlin
437 // binary search is -1
438
439 G4int limit = fVoxels.GetCandidatesVoxelArra
440 for (G4int i = 0 ; i < limit ; ++i)
441 {
442 G4int candidate = candidates[i];
443 G4VSolid& solid = *fSolids[candidate];
444 const G4Transform3D& transform = fTransfor
445
446 // The coordinates of the point are modifi
447 // solid local frame:
448 localPoint = GetLocalPoint(transform, aPoi
449 location = solid.Inside(localPoint);
450 if (location == EInside::kInside) return E
451 else if (location == EInside::kSurface)
452 {
453 G4MultiUnionSurface surface;
454 surface.point = localPoint;
455 surface.solid = &solid;
456 surfaces.push_back(surface);
457 }
458 }
459
460 ////////////////////////////////////////////
461 // Important comment: When two solids touch
462 // surface, the surface points will be consi
463 // located around will correspond to kInside
464
465 std::size_t size = surfaces.size();
466
467 if (size == 0)
468 {
469 return EInside::kOutside;
470 }
471
472 for (std::size_t i = 0; i < size - 1; ++i)
473 {
474 G4MultiUnionSurface& left = surfaces[i];
475 for (std::size_t j = i + 1; j < size; ++j)
476 {
477 G4MultiUnionSurface& right = surfaces[j]
478 G4ThreeVector n, n2;
479 n = left.solid->SurfaceNormal(left.point
480 n2 = right.solid->SurfaceNormal(right.po
481 if ((n + n2).mag2() < 1000 * kRadTolera
482 return EInside::kInside;
483 }
484 }
485
486 return EInside::kSurface;
487 }
488
489 //____________________________________________
490 EInside G4MultiUnion::Inside(const G4ThreeVect
491 {
492 // Classify point location with respect to s
493 // o eInside - inside the solid
494 // o eSurface - close to surface withi
495 // o eOutside - outside the solid
496
497 // Hitherto, it is considered that only para
498 // added to the container
499
500 // Implementation using voxelisation techniq
501 // -----------------------------------------
502
503 // return InsideIterator(aPoint);
504
505 EInside location = InsideWithExclusion(aPoin
506 return location;
507 }
508
509 //____________________________________________
510 EInside G4MultiUnion::InsideNoVoxels(const G4T
511 {
512 G4ThreeVector localPoint;
513 EInside location = EInside::kOutside;
514 G4int countSurface = 0;
515
516 auto numNodes = (G4int)fSolids.size();
517 for (auto i = 0 ; i < numNodes ; ++i)
518 {
519 G4VSolid& solid = *fSolids[i];
520 G4Transform3D transform = GetTransformatio
521
522 // The coordinates of the point are modifi
523 // intrinsic solid local frame:
524 localPoint = GetLocalPoint(transform, aPoi
525
526 location = solid.Inside(localPoint);
527
528 if (location == EInside::kSurface)
529 ++countSurface;
530
531 if (location == EInside::kInside) return E
532 }
533 if (countSurface != 0) return EInside::kSurf
534 return EInside::kOutside;
535 }
536
537 //____________________________________________
538 void G4MultiUnion::Extent(EAxis aAxis, G4doubl
539 {
540 // Determines the bounding box for the consi
541 G4ThreeVector min, max;
542
543 auto numNodes = (G4int)fSolids.size();
544 for (auto i = 0 ; i < numNodes ; ++i)
545 {
546 G4VSolid& solid = *fSolids[i];
547 G4Transform3D transform = GetTransformatio
548 solid.BoundingLimits(min, max);
549
550 TransformLimits(min, max, transform);
551
552 if (i == 0)
553 {
554 switch (aAxis)
555 {
556 case kXAxis:
557 aMin = min.x();
558 aMax = max.x();
559 break;
560 case kYAxis:
561 aMin = min.y();
562 aMax = max.y();
563 break;
564 case kZAxis:
565 aMin = min.z();
566 aMax = max.z();
567 break;
568 default:
569 break;
570 }
571 }
572 else
573 {
574 // Determine the min/max on the consider
575 switch (aAxis)
576 {
577 case kXAxis:
578 if (min.x() < aMin)
579 aMin = min.x();
580 if (max.x() > aMax)
581 aMax = max.x();
582 break;
583 case kYAxis:
584 if (min.y() < aMin)
585 aMin = min.y();
586 if (max.y() > aMax)
587 aMax = max.y();
588 break;
589 case kZAxis:
590 if (min.z() < aMin)
591 aMin = min.z();
592 if (max.z() > aMax)
593 aMax = max.z();
594 break;
595 default:
596 break;
597 }
598 }
599 }
600 }
601
602 //____________________________________________
603 void G4MultiUnion::BoundingLimits(G4ThreeVecto
604 G4ThreeVecto
605 {
606 Extent(kXAxis, aMin[0], aMax[0]);
607 Extent(kYAxis, aMin[1], aMax[1]);
608 Extent(kZAxis, aMin[2], aMax[2]);
609 }
610
611 //____________________________________________
612 G4bool
613 G4MultiUnion::CalculateExtent(const EAxis pAxi
614 const G4VoxelLim
615 const G4AffineTr
616 G4double&
617 {
618 G4ThreeVector bmin, bmax;
619
620 // Get bounding box
621 BoundingLimits(bmin,bmax);
622
623 // Find extent
624 G4BoundingEnvelope bbox(bmin,bmax);
625 return bbox.CalculateExtent(pAxis,pVoxelLimi
626 }
627
628 //____________________________________________
629 G4ThreeVector G4MultiUnion::SurfaceNormal(cons
630 {
631 // Computes the localNormal on a surface and
632 // Must return a valid vector. (even if the
633 //
634 // On an edge or corner, provide an averag
635 // tolerance
636 // NOTE: the tolerance value used in here
637 // tolerance - we will have to revis
638
639 std::vector<G4int> candidates;
640 G4ThreeVector localPoint, normal, localNorma
641 G4double safety = kInfinity;
642 G4int node = 0;
643
644 ////////////////////////////////////////////
645 // Important comment: Cases for which the po
646 // on a vertice remain to be treated
647
648 // determine weather we are in voxel area
649 if (fVoxels.GetCandidatesVoxelArray(aPoint,
650 {
651 std::size_t limit = candidates.size();
652 for (std::size_t i = 0 ; i < limit ; ++i)
653 {
654 G4int candidate = candidates[i];
655 const G4Transform3D& transform = fTransf
656
657 // The coordinates of the point are modi
658 // solid local frame:
659 localPoint = GetLocalPoint(transform, aP
660 G4VSolid& solid = *fSolids[candidate];
661 EInside location = solid.Inside(localPoi
662
663 if (location == EInside::kSurface)
664 {
665 // normal case when point is on surfac
666 normal = GetGlobalVector(transform, so
667 return normal.unit();
668 }
669 else
670 {
671 // collect the smallest safety and rem
672 G4double s = (location == EInside::kIn
673 ? solid.DistanceToOut(local
674 : solid.DistanceToIn(localP
675 if (s < safety)
676 {
677 safety = s;
678 node = candidate;
679 }
680 }
681 }
682 // on none of the solids, the point was no
683 G4VSolid& solid = *fSolids[node];
684 const G4Transform3D& transform = fTransfor
685 localPoint = GetLocalPoint(transform, aPoi
686
687 normal = GetGlobalVector(transform, solid.
688 return normal.unit();
689 }
690 else
691 {
692 // for the case when point is certainly ou
693
694 // find a solid in union with the smallest
695 node = SafetyFromOutsideNumberNode(aPoint,
696 G4VSolid& solid = *fSolids[node];
697
698 const G4Transform3D& transform = fTransfor
699 localPoint = GetLocalPoint(transform, aPoi
700
701 // evaluate normal for point at this found
702 // and transform multi-union coordinates
703 normal = GetGlobalVector(transform, solid.
704
705 return normal.unit();
706 }
707 }
708
709 //____________________________________________
710 G4double G4MultiUnion::DistanceToOut(const G4T
711 {
712 // Estimates isotropic distance to the surfa
713 // be either accurate or an underestimate.
714 // Two modes: - default/fast mode, sacrifici
715 // - "precise" mode, requests ac
716
717 std::vector<G4int> candidates;
718 G4ThreeVector localPoint;
719 G4double safetyMin = kInfinity;
720
721 // In general, the value return by DistanceT
722 // but only an undervalue (cf. overlaps)
723 fVoxels.GetCandidatesVoxelArray(point, candi
724
725 std::size_t limit = candidates.size();
726 for (std::size_t i = 0; i < limit; ++i)
727 {
728 G4int candidate = candidates[i];
729
730 // The coordinates of the point are modifi
731 // solid local frame:
732 const G4Transform3D& transform = fTransfor
733 localPoint = GetLocalPoint(transform, poin
734 G4VSolid& solid = *fSolids[candidate];
735 if (solid.Inside(localPoint) == EInside::k
736 {
737 G4double safety = solid.DistanceToOut(lo
738 if (safetyMin > safety) safetyMin = safe
739 }
740 }
741 if (safetyMin == kInfinity) safetyMin = 0; /
742
743 return safetyMin;
744 }
745
746 //____________________________________________
747 G4double G4MultiUnion::DistanceToIn(const G4Th
748 {
749 // Estimates the isotropic safety from a poi
750 // any of its surfaces. The algorithm may be
751 // underestimate.
752
753 if (!fAccurate) { return fVoxels.DistanceTo
754
755 const std::vector<G4VoxelBox>& boxes = fVoxe
756 G4double safetyMin = kInfinity;
757 G4ThreeVector localPoint;
758
759 std::size_t numNodes = fSolids.size();
760 for (std::size_t j = 0; j < numNodes; ++j)
761 {
762 G4ThreeVector dxyz;
763 if (j > 0)
764 {
765 const G4ThreeVector& pos = boxes[j].pos;
766 const G4ThreeVector& hlen = boxes[j].hle
767 for (auto i = 0; i <= 2; ++i)
768 // distance to middle point - hlength
769 // of x,y,z
770 if ((dxyz[i] = std::abs(point[i] - pos
771 continue;
772
773 G4double d2xyz = 0.;
774 for (auto i = 0; i <= 2; ++i)
775 if (dxyz[i] > 0) d2xyz += dxyz[i] * dx
776
777 // minimal distance is at least this, bu
778 // we can stop if previous was already l
779 // chance to be better tha previous valu
780 if (d2xyz >= safetyMin * safetyMin)
781 {
782 continue;
783 }
784 }
785 const G4Transform3D& transform = fTransfor
786 localPoint = GetLocalPoint(transform, poin
787 G4VSolid& solid = *fSolids[j];
788
789 G4double safety = solid.DistanceToIn(local
790 if (safety <= 0) return safety;
791 // it was detected, that the point is no
792 if (safetyMin > safety) safetyMin = safety
793 }
794 return safetyMin;
795 }
796
797 //____________________________________________
798 G4double G4MultiUnion::GetSurfaceArea()
799 {
800 if (fSurfaceArea == 0.0)
801 {
802 fSurfaceArea = EstimateSurfaceArea(1000000
803 }
804 return fSurfaceArea;
805 }
806
807 //____________________________________________
808 G4int G4MultiUnion::GetNumOfConstituents() con
809 {
810 G4int num = 0;
811 for (const auto solid : fSolids) { num += so
812 return num;
813 }
814
815 //____________________________________________
816 G4bool G4MultiUnion::IsFaceted() const
817 {
818 for (const auto solid : fSolids) { if (!soli
819 return true;
820 }
821
822 //____________________________________________
823 void G4MultiUnion::Voxelize()
824 {
825 fVoxels.Voxelize(fSolids, fTransformObjs);
826 }
827
828 //____________________________________________
829 G4int G4MultiUnion::SafetyFromOutsideNumberNod
830
831 {
832 // Method returning the closest node from a
833 // G4MultiUnion.
834 // This is used to compute the normal in the
835
836 const std::vector<G4VoxelBox>& boxes = fVoxe
837 safetyMin = kInfinity;
838 std::size_t safetyNode = 0;
839 G4ThreeVector localPoint;
840
841 std::size_t numNodes = fSolids.size();
842 for (std::size_t i = 0; i < numNodes; ++i)
843 {
844 G4double d2xyz = 0.;
845 G4double dxyz0 = std::abs(aPoint.x() - box
846 if (dxyz0 > safetyMin) continue;
847 G4double dxyz1 = std::abs(aPoint.y() - box
848 if (dxyz1 > safetyMin) continue;
849 G4double dxyz2 = std::abs(aPoint.z() - box
850 if (dxyz2 > safetyMin) continue;
851
852 if (dxyz0 > 0) d2xyz += dxyz0 * dxyz0;
853 if (dxyz1 > 0) d2xyz += dxyz1 * dxyz1;
854 if (dxyz2 > 0) d2xyz += dxyz2 * dxyz2;
855 if (d2xyz >= safetyMin * safetyMin) contin
856
857 G4VSolid& solid = *fSolids[i];
858 const G4Transform3D& transform = fTransfor
859 localPoint = GetLocalPoint(transform, aPoi
860 fAccurate = true;
861 G4double safety = solid.DistanceToIn(local
862 fAccurate = false;
863 if (safetyMin > safety)
864 {
865 safetyMin = safety;
866 safetyNode = i;
867 }
868 }
869 return (G4int)safetyNode;
870 }
871
872 //____________________________________________
873 void G4MultiUnion::TransformLimits(G4ThreeVect
874 const G4Tra
875 {
876 // The goal of this method is to convert the
877 // (representing the bounding box of a given
878 // to the main frame, using "transformation"
879
880 G4ThreeVector vertices[8] = // Deteminat
881 { // the exten
882 G4ThreeVector(min.x(), min.y(), min.z()),
883 G4ThreeVector(min.x(), max.y(), min.z()),
884 G4ThreeVector(max.x(), max.y(), min.z()),
885 G4ThreeVector(max.x(), min.y(), min.z()),
886 G4ThreeVector(min.x(), min.y(), max.z()),
887 G4ThreeVector(min.x(), max.y(), max.z()),
888 G4ThreeVector(max.x(), max.y(), max.z()),
889 G4ThreeVector(max.x(), min.y(), max.z())
890 };
891
892 min.set(kInfinity,kInfinity,kInfinity);
893 max.set(-kInfinity,-kInfinity,-kInfinity);
894
895 // Loop on th vertices
896 G4int limit = sizeof(vertices) / sizeof(G4Th
897 for (G4int i = 0 ; i < limit; ++i)
898 {
899 // From local frame to the global one:
900 // Current positions on the three axis:
901 G4ThreeVector current = GetGlobalPoint(tra
902
903 // If need be, replacement of the min & ma
904 if (current.x() > max.x()) max.setX(curren
905 if (current.x() < min.x()) min.setX(curren
906
907 if (current.y() > max.y()) max.setY(curren
908 if (current.y() < min.y()) min.setY(curren
909
910 if (current.z() > max.z()) max.setZ(curren
911 if (current.z() < min.z()) min.setZ(curren
912 }
913 }
914
915 // Stream object contents to an output stream
916 //____________________________________________
917 std::ostream& G4MultiUnion::StreamInfo(std::os
918 {
919 G4long oldprc = os.precision(16);
920 os << "-------------------------------------
921 << " *** Dump for solid -
922 << " =====================
923 << " Solid type: G4MultiUnion\n"
924 << " Parameters: \n";
925 std::size_t numNodes = fSolids.size();
926 for (std::size_t i = 0 ; i < numNodes ; +
927 {
928 G4VSolid& solid = *fSolids[i];
929 solid.StreamInfo(os);
930 const G4Transform3D& transform = fTransf
931 os << " Translation is " << transform.ge
932 os << " Rotation is :" << " \n";
933 os << " " << transform.getRotation() <<
934 }
935 os << " \n"
936 << "-------------------------------------
937 os.precision(oldprc);
938
939 return os;
940 }
941
942 //____________________________________________
943 G4ThreeVector G4MultiUnion::GetPointOnSurface(
944 {
945 G4ThreeVector point;
946
947 G4long size = fSolids.size();
948
949 do
950 {
951 G4long rnd = G4RandFlat::shootInt(G4long(0
952 G4VSolid& solid = *fSolids[rnd];
953 point = solid.GetPointOnSurface();
954 const G4Transform3D& transform = fTransfor
955 point = GetGlobalPoint(transform, point);
956 }
957 while (Inside(point) != EInside::kSurface);
958
959 return point;
960 }
961
962 //____________________________________________
963 void
964 G4MultiUnion::DescribeYourselfTo ( G4VGraphics
965 {
966 scene.AddSolid (*this);
967 }
968
969 //____________________________________________
970 G4Polyhedron* G4MultiUnion::CreatePolyhedron()
971 {
972 if (G4BooleanSolid::GetExternalBooleanProces
973 {
974 HepPolyhedronProcessor processor;
975 HepPolyhedronProcessor::Operation operatio
976
977 G4VSolid* solidA = GetSolid(0);
978 const G4Transform3D transform0 = GetTransf
979 G4DisplacedSolid dispSolidA("placedA", sol
980
981 auto top = new G4Polyhedron(*dispSolidA.Ge
982
983 for (G4int i = 1; i < GetNumberOfSolids();
984 {
985 G4VSolid* solidB = GetSolid(i);
986 const G4Transform3D transform = GetTrans
987 G4DisplacedSolid dispSolidB("placedB", s
988 G4Polyhedron* operand = dispSolidB.GetPo
989 processor.push_back(operation, *operand)
990 }
991
992 if (processor.execute(*top))
993 {
994 return top;
995 }
996 else
997 {
998 return nullptr;
999 }
1000 }
1001 else
1002 {
1003 return G4BooleanSolid::GetExternalBoolean
1004 }
1005 }
1006
1007 //___________________________________________
1008 G4Polyhedron* G4MultiUnion::GetPolyhedron() c
1009 {
1010 if (fpPolyhedron == nullptr ||
1011 fRebuildPolyhedron ||
1012 fpPolyhedron->GetNumberOfRotationStepsA
1013 fpPolyhedron->GetNumberOfRotationSteps(
1014 {
1015 G4AutoLock l(&polyhedronMutex);
1016 delete fpPolyhedron;
1017 fpPolyhedron = CreatePolyhedron();
1018 fRebuildPolyhedron = false;
1019 l.unlock();
1020 }
1021 return fpPolyhedron;
1022 }
1023