| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Implementation of G4MultiUnion class
27 //
28 // 19.10.12 M.Gayer - Original implementation from USolids module
29 // 06.04.17 G.Cosmo - Adapted implementation in Geant4 for VecGeom migration
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_INITIALIZER;
54 }
55
56 //______________________________________________________________________________
57 G4MultiUnion::G4MultiUnion(const G4String& name)
58 : G4VSolid(name)
59 {
60 SetName(name);
61 fSolids.clear();
62 fTransformObjs.clear();
63 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
64 }
65
66 //______________________________________________________________________________
67 G4MultiUnion::~G4MultiUnion()
68 = default;
69
70 //______________________________________________________________________________
71 void G4MultiUnion::AddNode(G4VSolid& solid, const G4Transform3D& trans)
72 {
73 fSolids.push_back(&solid);
74 fTransformObjs.push_back(trans); // Store a local copy of transformations
75 }
76
77 //______________________________________________________________________________
78 void G4MultiUnion::AddNode(G4VSolid* solid, const G4Transform3D& trans)
79 {
80 fSolids.push_back(solid);
81 fTransformObjs.push_back(trans); // Store a local copy of transformations
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& rhs)
93 : G4VSolid(rhs), fCubicVolume(rhs.fCubicVolume),
94 fSurfaceArea(rhs.fSurfaceArea),
95 kRadTolerance(rhs.kRadTolerance), fAccurate(rhs.fAccurate)
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 G4MultiUnion& rhs)
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, 0.001);
130 }
131 return fCubicVolume;
132 }
133
134 //______________________________________________________________________________
135 G4double
136 G4MultiUnion::DistanceToInNoVoxels(const G4ThreeVector& aPoint,
137 const G4ThreeVector& aDirection) const
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 = fTransformObjs[i];
148
149 localPoint = GetLocalPoint(transform, aPoint);
150 localDirection = GetLocalVector(transform, direction);
151
152 G4double distance = solid.DistanceToIn(localPoint, localDirection);
153 if (minDistance > distance) minDistance = distance;
154 }
155 return minDistance;
156 }
157
158 //______________________________________________________________________________
159 G4double G4MultiUnion::DistanceToInCandidates(const G4ThreeVector& aPoint,
160 const G4ThreeVector& direction,
161 std::vector<G4int>& candidates,
162 G4SurfBits& bits) const
163 {
164 std::size_t candidatesCount = candidates.size();
165 G4ThreeVector localPoint, localDirection;
166
167 G4double minDistance = kInfinity;
168 for (std::size_t i = 0 ; i < candidatesCount; ++i)
169 {
170 G4int candidate = candidates[i];
171 G4VSolid& solid = *fSolids[candidate];
172 const G4Transform3D& transform = fTransformObjs[candidate];
173
174 localPoint = GetLocalPoint(transform, aPoint);
175 localDirection = GetLocalVector(transform, direction);
176 G4double distance = solid.DistanceToIn(localPoint, localDirection);
177 if (minDistance > distance) minDistance = distance;
178 bits.SetBitNumber(candidate);
179 if (minDistance == 0) break;
180 }
181 return minDistance;
182 }
183
184 // Algorithm note: we have to look also for all other objects in next voxels,
185 // if the distance is not shorter ... we have to do it because,
186 // for example for objects which starts in first voxel in which they
187 // do not collide with direction line, but in second it collides...
188 // The idea of crossing voxels would be still applicable,
189 // because this way we could exclude from the testing such solids,
190 // which were found that obviously are not good candidates, because
191 // they would return infinity
192 // But if distance is smaller than the shift to next voxel, we can return
193 // it immediately
194 //______________________________________________________________________________
195 G4double G4MultiUnion::DistanceToIn(const G4ThreeVector& aPoint,
196 const G4ThreeVector& aDirection) const
197 {
198 G4double minDistance = kInfinity;
199 G4ThreeVector direction = aDirection.unit();
200 G4double shift = fVoxels.DistanceToFirst(aPoint, direction);
201 if (shift == kInfinity) return shift;
202
203 G4ThreeVector currentPoint = aPoint;
204 if (shift != 0.0) currentPoint += direction * shift;
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(curVoxel, candidates, &exclusion) != 0)
214 {
215 G4double distance = DistanceToInCandidates(aPoint, direction,
216 candidates, exclusion);
217 if (minDistance > distance) minDistance = distance;
218 if (distance < shift) break;
219 }
220 }
221 shift = fVoxels.DistanceToNext(aPoint, direction, curVoxel);
222 }
223 while (minDistance > shift);
224
225 return minDistance;
226 }
227
228 //______________________________________________________________________________
229 G4double G4MultiUnion::DistanceToOutNoVoxels(const G4ThreeVector& aPoint,
230 const G4ThreeVector& aDirection,
231 G4ThreeVector* aNormal) const
232 {
233 // Computes distance from a point presumably outside the solid to the solid
234 // surface. Ignores first surface if the point is actually inside.
235 // Early return infinity in case the safety to any surface is found greater
236 // than the proposed step aPstep.
237 // The normal vector to the crossed surface is filled only in case the box
238 // is crossed, otherwise aNormal->IsNull() is true.
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 = fTransformObjs[i];
254 localPoint = GetLocalPoint(transform, currentPoint);
255 localDirection = GetLocalVector(transform, direction);
256 EInside location = solid.Inside(localPoint);
257 if (location != EInside::kOutside)
258 {
259 G4double distance = solid.DistanceToOut(localPoint, localDirection,
260 false, nullptr, aNormal);
261 if (distance < kInfinity)
262 {
263 if (resultDistToOut == kInfinity) resultDistToOut = 0;
264 if (distance > 0)
265 {
266 currentPoint = GetGlobalPoint(transform, localPoint
267 + distance*localDirection);
268 resultDistToOut += distance;
269 ignoredSolid = i; // skip the solid which we have just left
270 i = -1; // force the loop to continue from 0
271 }
272 }
273 }
274 }
275 }
276 return resultDistToOut;
277 }
278
279 //______________________________________________________________________________
280 G4double G4MultiUnion::DistanceToOut(const G4ThreeVector& aPoint,
281 const G4ThreeVector& aDirection,
282 const G4bool /* calcNorm */,
283 G4bool* /* validNorm */,
284 G4ThreeVector* aNormal) const
285 {
286 return DistanceToOutVoxels(aPoint, aDirection, aNormal);
287 }
288
289 //______________________________________________________________________________
290 G4double G4MultiUnion::DistanceToOutVoxels(const G4ThreeVector& aPoint,
291 const G4ThreeVector& aDirection,
292 G4ThreeVector* aNormal) const
293 {
294 // Computes distance from a point presumably inside the solid to the solid
295 // surface. Ignores first surface along each axis systematically (for points
296 // inside or outside. Early returns zero in case the second surface is behind
297 // the starting point.
298 // o The proposed step is ignored.
299 // o The normal vector to the crossed surface is always filled.
300
301 // In the case the considered point is located inside the G4MultiUnion
302 // structure, the treatments are as follows:
303 // - investigation of the candidates for the passed point
304 // - progressive moving of the point towards the surface, along the
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, candidates) != 0)
315 {
316 // For normal case for which we presume the point is inside
317 G4ThreeVector localPoint, localDirection, localNormal;
318 G4ThreeVector currentPoint = aPoint;
319 G4SurfBits exclusion(fVoxels.GetBitsPerSlice());
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 you just got out of) since
336 // numerically the propagated point will be on its surface
337
338 G4VSolid& solid = *fSolids[candidate];
339 const G4Transform3D& transform = fTransformObjs[candidate];
340
341 // The coordinates of the point are modified so as to fit the
342 // intrinsic solid local frame:
343 localPoint = GetLocalPoint(transform, currentPoint);
344
345 // DistanceToOut at least for Trd sometimes return non-zero value
346 // even from points that are outside. Therefore, this condition
347 // must currently be here, otherwise it would not work.
348 // But it means it would be slower.
349
350 if (solid.Inside(localPoint) != EInside::kOutside)
351 {
352 notOutside = true;
353
354 localDirection = GetLocalVector(transform, direction);
355
356 // propagate with solid.DistanceToOut
357 G4double shift = solid.DistanceToOut(localPoint, localDirection,
358 false, nullptr, &localNormal);
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 = fTransformObjs[maxCandidate];
371
372 // convert from local normal
373 if (aNormal != nullptr) *aNormal = GetGlobalVector(transform, maxNormal);
374
375 distance += maxDistance;
376 currentPoint += maxDistance * direction;
377 if(maxDistance == 0.) ++count;
378
379 // the current component will be ignored
380 exclusion.SetBitNumber(maxCandidate);
381 EInside location = InsideWithExclusion(currentPoint, &exclusion);
382
383 // perform a Inside
384 // it should be excluded current solid from checking
385 // we have to collect the maximum distance from all given candidates.
386 // such "maximum" candidate should be then used for finding next
387 // candidates
388 if (location == EInside::kOutside)
389 {
390 // else return cumulated distances to outside of the traversed
391 // components
392 break;
393 }
394 // if inside another component, redo 1 to 3 but add the next
395 // DistanceToOut on top of the previous.
396
397 // and fill the candidates for the corresponding voxel (just
398 // exiting current component along direction)
399 candidates.clear();
400
401 fVoxels.GetCandidatesVoxelArray(currentPoint, candidates, &exclusion);
402 exclusion.ResetBitNumber(maxCandidate);
403 }
404 }
405 while ((notOutside) && (count < numNodes));
406 }
407
408 return distance;
409 }
410
411 //______________________________________________________________________________
412 EInside G4MultiUnion::InsideWithExclusion(const G4ThreeVector& aPoint,
413 G4SurfBits* exclusion) const
414 {
415 // Classify point location with respect to solid:
416 // o eInside - inside the solid
417 // o eSurface - close to surface within tolerance
418 // o eOutside - outside the solid
419
420 // Hitherto, it is considered that only parallelepipedic nodes
421 // can be added to the container
422
423 // Implementation using voxelisation techniques:
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 measure performance
433 // TODO: getPointIndex should not be used, instead GetVoxel + GetVoxelsIndex
434 // should be used
435 // TODO: than pass result to GetVoxel further to GetCandidatesVoxelArray
436 // TODO: eventually GetVoxel should be inlined here, early exit if any
437 // binary search is -1
438
439 G4int limit = fVoxels.GetCandidatesVoxelArray(aPoint, candidates, exclusion);
440 for (G4int i = 0 ; i < limit ; ++i)
441 {
442 G4int candidate = candidates[i];
443 G4VSolid& solid = *fSolids[candidate];
444 const G4Transform3D& transform = fTransformObjs[candidate];
445
446 // The coordinates of the point are modified so as to fit the intrinsic
447 // solid local frame:
448 localPoint = GetLocalPoint(transform, aPoint);
449 location = solid.Inside(localPoint);
450 if (location == EInside::kInside) return EInside::kInside;
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 each other along a flat
462 // surface, the surface points will be considered as kSurface, while points
463 // located around will correspond to kInside (cf. G4UnionSolid)
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.point);
481 if ((n + n2).mag2() < 1000 * kRadTolerance)
482 return EInside::kInside;
483 }
484 }
485
486 return EInside::kSurface;
487 }
488
489 //______________________________________________________________________________
490 EInside G4MultiUnion::Inside(const G4ThreeVector& aPoint) const
491 {
492 // Classify point location with respect to solid:
493 // o eInside - inside the solid
494 // o eSurface - close to surface within tolerance
495 // o eOutside - outside the solid
496
497 // Hitherto, it is considered that only parallelepipedic nodes can be
498 // added to the container
499
500 // Implementation using voxelisation techniques:
501 // ---------------------------------------------
502
503 // return InsideIterator(aPoint);
504
505 EInside location = InsideWithExclusion(aPoint);
506 return location;
507 }
508
509 //______________________________________________________________________________
510 EInside G4MultiUnion::InsideNoVoxels(const G4ThreeVector& aPoint) const
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 = GetTransformation(i);
521
522 // The coordinates of the point are modified so as to fit the
523 // intrinsic solid local frame:
524 localPoint = GetLocalPoint(transform, aPoint);
525
526 location = solid.Inside(localPoint);
527
528 if (location == EInside::kSurface)
529 ++countSurface;
530
531 if (location == EInside::kInside) return EInside::kInside;
532 }
533 if (countSurface != 0) return EInside::kSurface;
534 return EInside::kOutside;
535 }
536
537 //______________________________________________________________________________
538 void G4MultiUnion::Extent(EAxis aAxis, G4double& aMin, G4double& aMax) const
539 {
540 // Determines the bounding box for the considered instance of "UMultipleUnion"
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 = GetTransformation(i);
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 considered axis:
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(G4ThreeVector& aMin,
604 G4ThreeVector& aMax) const
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 pAxis,
614 const G4VoxelLimits& pVoxelLimit,
615 const G4AffineTransform& pTransform,
616 G4double& pMin, G4double& pMax) const
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,pVoxelLimit,pTransform,pMin,pMax);
626 }
627
628 //______________________________________________________________________________
629 G4ThreeVector G4MultiUnion::SurfaceNormal(const G4ThreeVector& aPoint) const
630 {
631 // Computes the localNormal on a surface and returns it as a unit vector.
632 // Must return a valid vector. (even if the point is not on the surface).
633 //
634 // On an edge or corner, provide an average localNormal of all facets within
635 // tolerance
636 // NOTE: the tolerance value used in here is not yet the global surface
637 // tolerance - we will have to revise this value - TODO
638
639 std::vector<G4int> candidates;
640 G4ThreeVector localPoint, normal, localNormal;
641 G4double safety = kInfinity;
642 G4int node = 0;
643
644 ///////////////////////////////////////////////////////////////////////////
645 // Important comment: Cases for which the point is located on an edge or
646 // on a vertice remain to be treated
647
648 // determine weather we are in voxel area
649 if (fVoxels.GetCandidatesVoxelArray(aPoint, candidates) != 0)
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 = fTransformObjs[candidate];
656
657 // The coordinates of the point are modified so as to fit the intrinsic
658 // solid local frame:
659 localPoint = GetLocalPoint(transform, aPoint);
660 G4VSolid& solid = *fSolids[candidate];
661 EInside location = solid.Inside(localPoint);
662
663 if (location == EInside::kSurface)
664 {
665 // normal case when point is on surface, we pick first solid
666 normal = GetGlobalVector(transform, solid.SurfaceNormal(localPoint));
667 return normal.unit();
668 }
669 else
670 {
671 // collect the smallest safety and remember solid node
672 G4double s = (location == EInside::kInside)
673 ? solid.DistanceToOut(localPoint)
674 : solid.DistanceToIn(localPoint);
675 if (s < safety)
676 {
677 safety = s;
678 node = candidate;
679 }
680 }
681 }
682 // on none of the solids, the point was not on the surface
683 G4VSolid& solid = *fSolids[node];
684 const G4Transform3D& transform = fTransformObjs[node];
685 localPoint = GetLocalPoint(transform, aPoint);
686
687 normal = GetGlobalVector(transform, solid.SurfaceNormal(localPoint));
688 return normal.unit();
689 }
690 else
691 {
692 // for the case when point is certainly outside:
693
694 // find a solid in union with the smallest safety
695 node = SafetyFromOutsideNumberNode(aPoint, safety);
696 G4VSolid& solid = *fSolids[node];
697
698 const G4Transform3D& transform = fTransformObjs[node];
699 localPoint = GetLocalPoint(transform, aPoint);
700
701 // evaluate normal for point at this found solid
702 // and transform multi-union coordinates
703 normal = GetGlobalVector(transform, solid.SurfaceNormal(localPoint));
704
705 return normal.unit();
706 }
707 }
708
709 //______________________________________________________________________________
710 G4double G4MultiUnion::DistanceToOut(const G4ThreeVector& point) const
711 {
712 // Estimates isotropic distance to the surface of the solid. This must
713 // be either accurate or an underestimate.
714 // Two modes: - default/fast mode, sacrificing accuracy for speed
715 // - "precise" mode, requests accurate value if available.
716
717 std::vector<G4int> candidates;
718 G4ThreeVector localPoint;
719 G4double safetyMin = kInfinity;
720
721 // In general, the value return by DistanceToIn(p) will not be the exact
722 // but only an undervalue (cf. overlaps)
723 fVoxels.GetCandidatesVoxelArray(point, candidates);
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 modified so as to fit the intrinsic
731 // solid local frame:
732 const G4Transform3D& transform = fTransformObjs[candidate];
733 localPoint = GetLocalPoint(transform, point);
734 G4VSolid& solid = *fSolids[candidate];
735 if (solid.Inside(localPoint) == EInside::kInside)
736 {
737 G4double safety = solid.DistanceToOut(localPoint);
738 if (safetyMin > safety) safetyMin = safety;
739 }
740 }
741 if (safetyMin == kInfinity) safetyMin = 0; // we are not inside
742
743 return safetyMin;
744 }
745
746 //______________________________________________________________________________
747 G4double G4MultiUnion::DistanceToIn(const G4ThreeVector& point) const
748 {
749 // Estimates the isotropic safety from a point outside the current solid to
750 // any of its surfaces. The algorithm may be accurate or should provide a fast
751 // underestimate.
752
753 if (!fAccurate) { return fVoxels.DistanceToBoundingBox(point); }
754
755 const std::vector<G4VoxelBox>& boxes = fVoxels.GetBoxes();
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].hlen;
767 for (auto i = 0; i <= 2; ++i)
768 // distance to middle point - hlength => distance from point to border
769 // of x,y,z
770 if ((dxyz[i] = std::abs(point[i] - pos[i]) - hlen[i]) > safetyMin)
771 continue;
772
773 G4double d2xyz = 0.;
774 for (auto i = 0; i <= 2; ++i)
775 if (dxyz[i] > 0) d2xyz += dxyz[i] * dxyz[i];
776
777 // minimal distance is at least this, but could be even higher. therefore,
778 // we can stop if previous was already lower, let us check if it does any
779 // chance to be better tha previous values...
780 if (d2xyz >= safetyMin * safetyMin)
781 {
782 continue;
783 }
784 }
785 const G4Transform3D& transform = fTransformObjs[j];
786 localPoint = GetLocalPoint(transform, point);
787 G4VSolid& solid = *fSolids[j];
788
789 G4double safety = solid.DistanceToIn(localPoint);
790 if (safety <= 0) return safety;
791 // it was detected, that the point is not located outside
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, 0.001);
803 }
804 return fSurfaceArea;
805 }
806
807 //______________________________________________________________________________
808 G4int G4MultiUnion::GetNumOfConstituents() const
809 {
810 G4int num = 0;
811 for (const auto solid : fSolids) { num += solid->GetNumOfConstituents(); }
812 return num;
813 }
814
815 //______________________________________________________________________________
816 G4bool G4MultiUnion::IsFaceted() const
817 {
818 for (const auto solid : fSolids) { if (!solid->IsFaceted()) return false; }
819 return true;
820 }
821
822 //______________________________________________________________________________
823 void G4MultiUnion::Voxelize()
824 {
825 fVoxels.Voxelize(fSolids, fTransformObjs);
826 }
827
828 //______________________________________________________________________________
829 G4int G4MultiUnion::SafetyFromOutsideNumberNode(const G4ThreeVector& aPoint,
830 G4double& safetyMin) const
831 {
832 // Method returning the closest node from a point located outside a
833 // G4MultiUnion.
834 // This is used to compute the normal in the case no candidate has been found.
835
836 const std::vector<G4VoxelBox>& boxes = fVoxels.GetBoxes();
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() - boxes[i].pos.x()) - boxes[i].hlen.x();
846 if (dxyz0 > safetyMin) continue;
847 G4double dxyz1 = std::abs(aPoint.y() - boxes[i].pos.y()) - boxes[i].hlen.y();
848 if (dxyz1 > safetyMin) continue;
849 G4double dxyz2 = std::abs(aPoint.z() - boxes[i].pos.z()) - boxes[i].hlen.z();
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) continue;
856
857 G4VSolid& solid = *fSolids[i];
858 const G4Transform3D& transform = fTransformObjs[i];
859 localPoint = GetLocalPoint(transform, aPoint);
860 fAccurate = true;
861 G4double safety = solid.DistanceToIn(localPoint);
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(G4ThreeVector& min, G4ThreeVector& max,
874 const G4Transform3D& transformation) const
875 {
876 // The goal of this method is to convert the quantities min and max
877 // (representing the bounding box of a given solid in its local frame)
878 // to the main frame, using "transformation"
879
880 G4ThreeVector vertices[8] = // Detemination of the vertices thanks to
881 { // the extension of each solid:
882 G4ThreeVector(min.x(), min.y(), min.z()), // 1st vertice:
883 G4ThreeVector(min.x(), max.y(), min.z()), // 2nd vertice:
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(G4ThreeVector);
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(transformation, vertices[i]);
902
903 // If need be, replacement of the min & max values:
904 if (current.x() > max.x()) max.setX(current.x());
905 if (current.x() < min.x()) min.setX(current.x());
906
907 if (current.y() > max.y()) max.setY(current.y());
908 if (current.y() < min.y()) min.setY(current.y());
909
910 if (current.z() > max.z()) max.setZ(current.z());
911 if (current.z() < min.z()) min.setZ(current.z());
912 }
913 }
914
915 // Stream object contents to an output stream
916 //______________________________________________________________________________
917 std::ostream& G4MultiUnion::StreamInfo(std::ostream& os) const
918 {
919 G4long oldprc = os.precision(16);
920 os << "-----------------------------------------------------------\n"
921 << " *** Dump for solid - " << GetName() << " ***\n"
922 << " ===================================================\n"
923 << " Solid type: G4MultiUnion\n"
924 << " Parameters: \n";
925 std::size_t numNodes = fSolids.size();
926 for (std::size_t i = 0 ; i < numNodes ; ++i)
927 {
928 G4VSolid& solid = *fSolids[i];
929 solid.StreamInfo(os);
930 const G4Transform3D& transform = fTransformObjs[i];
931 os << " Translation is " << transform.getTranslation() << " \n";
932 os << " Rotation is :" << " \n";
933 os << " " << transform.getRotation() << "\n";
934 }
935 os << " \n"
936 << "-----------------------------------------------------------\n";
937 os.precision(oldprc);
938
939 return os;
940 }
941
942 //______________________________________________________________________________
943 G4ThreeVector G4MultiUnion::GetPointOnSurface() const
944 {
945 G4ThreeVector point;
946
947 G4long size = fSolids.size();
948
949 do
950 {
951 G4long rnd = G4RandFlat::shootInt(G4long(0), size);
952 G4VSolid& solid = *fSolids[rnd];
953 point = solid.GetPointOnSurface();
954 const G4Transform3D& transform = fTransformObjs[rnd];
955 point = GetGlobalPoint(transform, point);
956 }
957 while (Inside(point) != EInside::kSurface);
958
959 return point;
960 }
961
962 //______________________________________________________________________________
963 void
964 G4MultiUnion::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
965 {
966 scene.AddSolid (*this);
967 }
968
969 //______________________________________________________________________________
970 G4Polyhedron* G4MultiUnion::CreatePolyhedron() const
971 {
972 if (G4BooleanSolid::GetExternalBooleanProcessor() == nullptr)
973 {
974 HepPolyhedronProcessor processor;
975 HepPolyhedronProcessor::Operation operation = HepPolyhedronProcessor::UNION;
976
977 G4VSolid* solidA = GetSolid(0);
978 const G4Transform3D transform0 = GetTransformation(0);
979 G4DisplacedSolid dispSolidA("placedA", solidA, transform0);
980
981 auto top = new G4Polyhedron(*dispSolidA.GetPolyhedron());
982
983 for (G4int i = 1; i < GetNumberOfSolids(); ++i)
984 {
985 G4VSolid* solidB = GetSolid(i);
986 const G4Transform3D transform = GetTransformation(i);
987 G4DisplacedSolid dispSolidB("placedB", solidB, transform);
988 G4Polyhedron* operand = dispSolidB.GetPolyhedron();
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::GetExternalBooleanProcessor()->Process(this);
1004 }
1005 }
1006
1007 //______________________________________________________________________________
1008 G4Polyhedron* G4MultiUnion::GetPolyhedron() const
1009 {
1010 if (fpPolyhedron == nullptr ||
1011 fRebuildPolyhedron ||
1012 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
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