| 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 intell
19 // * Vanderbilt University Free Electron Laser
20 // * Vanderbilt University, Nashville, TN, USA
21 // * Development supported by:
22 // * United States MFEL program under grant F
23 // * and NASA under contract number NNG04CT05P
24 // * Written by Marcus H. Mendenhall and Rober
25 // *
26 // * Contributed to the Geant4 Core, January,
27 // *
28 // *******************************************
29 //
30 // Implementation for G4Tet class
31 //
32 // 03.09.2004 - Marcus Mendenhall, created
33 // 08.01.2020 - Evgueni Tcherniaev, complete r
34 // -------------------------------------------
35
36 #include "G4Tet.hh"
37
38 #if !defined(G4GEOM_USE_UTET)
39
40 #include "G4VoxelLimits.hh"
41 #include "G4AffineTransform.hh"
42 #include "G4BoundingEnvelope.hh"
43
44 #include "G4VPVParameterisation.hh"
45
46 #include "G4QuickRand.hh"
47
48 #include "G4VGraphicsScene.hh"
49 #include "G4Polyhedron.hh"
50 #include "G4VisExtent.hh"
51
52 #include "G4AutoLock.hh"
53
54 namespace
55 {
56 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE
57 }
58
59 using namespace CLHEP;
60
61 //////////////////////////////////////////////
62 //
63 // Constructor - create a tetrahedron
64 // A Tet has all of its geometrical informatio
65 //
66 G4Tet::G4Tet(const G4String& pName,
67 const G4ThreeVector& p0,
68 const G4ThreeVector& p1,
69 const G4ThreeVector& p2,
70 const G4ThreeVector& p3, G4bool*
71 : G4VSolid(pName)
72 {
73 // Check for degeneracy
74 G4bool degenerate = CheckDegeneracy(p0, p1,
75 if (degeneracyFlag != nullptr)
76 {
77 *degeneracyFlag = degenerate;
78 }
79 else if (degenerate)
80 {
81 std::ostringstream message;
82 message << "Degenerate tetrahedron: " << G
83 << " anchor: " << p0 << "\n"
84 << " p1 : " << p1 << "\n"
85 << " p2 : " << p2 << "\n"
86 << " p3 : " << p3 << "\n"
87 << " volume: "
88 << std::abs((p1 - p0).cross(p2 - p
89 G4Exception("G4Tet::G4Tet()", "GeomSolids0
90 }
91
92 // Define surface thickness
93 halfTolerance = 0.5 * kCarTolerance;
94
95 // Set data members
96 Initialize(p0, p1, p2, p3);
97 }
98
99 //////////////////////////////////////////////
100 //
101 // Fake default constructor - sets only member
102 // for usage restri
103 //
104 G4Tet::G4Tet( __void__& a )
105 : G4VSolid(a)
106 {
107 }
108
109 //////////////////////////////////////////////
110 //
111 // Destructor
112 //
113 G4Tet::~G4Tet()
114 {
115 delete fpPolyhedron; fpPolyhedron = nullptr;
116 }
117
118 //////////////////////////////////////////////
119 //
120 // Copy constructor
121 //
122 G4Tet::G4Tet(const G4Tet& rhs)
123 : G4VSolid(rhs)
124 {
125 halfTolerance = rhs.halfTolerance;
126 fCubicVolume = rhs.fCubicVolume;
127 fSurfaceArea = rhs.fSurfaceArea;
128 for (G4int i = 0; i < 4; ++i) { fVertex[i]
129 for (G4int i = 0; i < 4; ++i) { fNormal[i]
130 for (G4int i = 0; i < 4; ++i) { fDist[i] =
131 for (G4int i = 0; i < 4; ++i) { fArea[i] =
132 fBmin = rhs.fBmin;
133 fBmax = rhs.fBmax;
134 }
135
136 //////////////////////////////////////////////
137 //
138 // Assignment operator
139 //
140 G4Tet& G4Tet::operator = (const G4Tet& rhs)
141 {
142 // Check assignment to self
143 //
144 if (this == &rhs) { return *this; }
145
146 // Copy base class data
147 //
148 G4VSolid::operator=(rhs);
149
150 // Copy data
151 //
152 halfTolerance = rhs.halfTolerance;
153 fCubicVolume = rhs.fCubicVolume;
154 fSurfaceArea = rhs.fSurfaceArea;
155 for (G4int i = 0; i < 4; ++i) { fVertex[i]
156 for (G4int i = 0; i < 4; ++i) { fNormal[i]
157 for (G4int i = 0; i < 4; ++i) { fDist[i] =
158 for (G4int i = 0; i < 4; ++i) { fArea[i] =
159 fBmin = rhs.fBmin;
160 fBmax = rhs.fBmax;
161 fRebuildPolyhedron = false;
162 delete fpPolyhedron; fpPolyhedron = nullptr
163
164 return *this;
165 }
166
167 //////////////////////////////////////////////
168 //
169 // Return true if tetrahedron is degenerate
170 // Tetrahedron is concidered as degenerate in
171 // height is less than degeneracy tolerance
172 //
173 G4bool G4Tet::CheckDegeneracy(const G4ThreeVec
174 const G4ThreeVec
175 const G4ThreeVec
176 const G4ThreeVec
177 {
178 G4double hmin = 4. * kCarTolerance; // degen
179
180 // Calculate volume
181 G4double vol = std::abs((p1 - p0).cross(p2 -
182
183 // Calculate face areas squared
184 G4double ss[4];
185 ss[0] = ((p1 - p0).cross(p2 - p0)).mag2();
186 ss[1] = ((p2 - p0).cross(p3 - p0)).mag2();
187 ss[2] = ((p3 - p0).cross(p1 - p0)).mag2();
188 ss[3] = ((p2 - p1).cross(p3 - p1)).mag2();
189
190 // Find face with max area
191 G4int k = 0;
192 for (G4int i = 1; i < 4; ++i) { if (ss[i] >
193
194 // Check: vol^2 / s^2 <= hmin^2
195 return (vol*vol <= ss[k]*hmin*hmin);
196 }
197
198 //////////////////////////////////////////////
199 //
200 // Set data members
201 //
202 void G4Tet::Initialize(const G4ThreeVector& p0
203 const G4ThreeVector& p1
204 const G4ThreeVector& p2
205 const G4ThreeVector& p3
206 {
207 // Set vertices
208 fVertex[0] = p0;
209 fVertex[1] = p1;
210 fVertex[2] = p2;
211 fVertex[3] = p3;
212
213 G4ThreeVector norm[4];
214 norm[0] = (p2 - p0).cross(p1 - p0);
215 norm[1] = (p3 - p0).cross(p2 - p0);
216 norm[2] = (p1 - p0).cross(p3 - p0);
217 norm[3] = (p2 - p1).cross(p3 - p1);
218 G4double volume = norm[0].dot(p3 - p0);
219 if (volume > 0.)
220 {
221 for (auto & i : norm) { i = -i; }
222 }
223
224 // Set normals to face planes
225 for (G4int i = 0; i < 4; ++i) { fNormal[i] =
226
227 // Set distances to planes
228 for (G4int i = 0; i < 3; ++i) { fDist[i] = f
229 fDist[3] = fNormal[3].dot(p1);
230
231 // Set face areas
232 for (G4int i = 0; i < 4; ++i) { fArea[i] = 0
233
234 // Set bounding box
235 for (G4int i = 0; i < 3; ++i)
236 {
237 fBmin[i] = std::min(std::min(std::min(p0[i
238 fBmax[i] = std::max(std::max(std::max(p0[i
239 }
240
241 // Set volume and surface area
242 fCubicVolume = std::abs(volume)/6.;
243 fSurfaceArea = fArea[0] + fArea[1] + fArea[2
244 }
245
246 //////////////////////////////////////////////
247 //
248 // Set vertices
249 //
250 void G4Tet::SetVertices(const G4ThreeVector& p
251 const G4ThreeVector& p
252 const G4ThreeVector& p
253 const G4ThreeVector& p
254 {
255 // Check for degeneracy
256 G4bool degenerate = CheckDegeneracy(p0, p1,
257 if (degeneracyFlag != nullptr)
258 {
259 *degeneracyFlag = degenerate;
260 }
261 else if (degenerate)
262 {
263 std::ostringstream message;
264 message << "Degenerate tetrahedron is not
265 << " anchor: " << p0 << "\n"
266 << " p1 : " << p1 << "\n"
267 << " p2 : " << p2 << "\n"
268 << " p3 : " << p3 << "\n"
269 << " volume: "
270 << std::abs((p1 - p0).cross(p2 - p
271 G4Exception("G4Tet::SetVertices()", "GeomS
272 FatalException, message);
273 }
274
275 // Set data members
276 Initialize(p0, p1, p2, p3);
277
278 // Set flag to rebuild polyhedron
279 fRebuildPolyhedron = true;
280 }
281
282 //////////////////////////////////////////////
283 //
284 // Return four vertices
285 //
286 void G4Tet::GetVertices(G4ThreeVector& p0,
287 G4ThreeVector& p1,
288 G4ThreeVector& p2,
289 G4ThreeVector& p3) con
290 {
291 p0 = fVertex[0];
292 p1 = fVertex[1];
293 p2 = fVertex[2];
294 p3 = fVertex[3];
295 }
296
297 //////////////////////////////////////////////
298 //
299 // Return std::vector of vertices
300 //
301 std::vector<G4ThreeVector> G4Tet::GetVertices(
302 {
303 std::vector<G4ThreeVector> vertices(4);
304 for (G4int i = 0; i < 4; ++i) { vertices[i]
305 return vertices;
306 }
307
308 //////////////////////////////////////////////
309 //
310 // Dispatch to parameterisation for replicatio
311 // computation & modification.
312 //
313 void G4Tet::ComputeDimensions(G4VPVParameteris
314 const G4int ,
315 const G4VPhysica
316 {
317 }
318
319 //////////////////////////////////////////////
320 //
321 // Set bounding box
322 //
323 void G4Tet::SetBoundingLimits(const G4ThreeVec
324 const G4ThreeVec
325 {
326 G4int iout[4] = { 0, 0, 0, 0 };
327 for (G4int i = 0; i < 4; ++i)
328 {
329 iout[i] = (G4int)(fVertex[i].x() < pMin.x(
330 fVertex[i].y() < pMin.y(
331 fVertex[i].z() < pMin.z(
332 fVertex[i].x() > pMax.x(
333 fVertex[i].y() > pMax.y(
334 fVertex[i].z() > pMax.z(
335 }
336 if (iout[0] + iout[1] + iout[2] + iout[3] !=
337 {
338 std::ostringstream message;
339 message << "Attempt to set bounding box th
340 << GetName() << " !\n"
341 << " Specified bounding box limit
342 << " pmin: " << pMin << "\n"
343 << " pmax: " << pMax << "\n"
344 << " Tetrahedron vertices:\n"
345 << " anchor " << fVertex[0] <<
346 << " p1 " << fVertex[1] <<
347 << " p2 " << fVertex[2] <<
348 << " p3 " << fVertex[3] <<
349 G4Exception("G4Tet::SetBoundingLimits()",
350 FatalException, message);
351 }
352 fBmin = pMin;
353 fBmax = pMax;
354 }
355
356 //////////////////////////////////////////////
357 //
358 // Return bounding box
359 //
360 void G4Tet::BoundingLimits(G4ThreeVector& pMin
361 {
362 pMin = fBmin;
363 pMax = fBmax;
364 }
365
366 //////////////////////////////////////////////
367 //
368 // Calculate extent under transform and specif
369 //
370 G4bool G4Tet::CalculateExtent(const EAxis pAxi
371 const G4VoxelLim
372 const G4AffineTr
373 G4double&
374 {
375 G4ThreeVector bmin, bmax;
376
377 // Check bounding box (bbox)
378 //
379 BoundingLimits(bmin,bmax);
380 G4BoundingEnvelope bbox(bmin,bmax);
381
382 // Use simple bounding-box to help in the ca
383 //
384 return bbox.CalculateExtent(pAxis,pVoxelLimi
385
386 #if 0
387 // Precise extent computation (disabled by d
388 //
389 G4bool exist;
390 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox
391 {
392 return exist = (pMin < pMax) ? true : fals
393 }
394
395 // Set bounding envelope (benv) and calculat
396 //
397 std::vector<G4ThreeVector> vec = GetVertices
398
399 G4ThreeVectorList anchor(1);
400 anchor[0].set(vec[0].x(),vec[0].y(),vec[0].z
401
402 G4ThreeVectorList base(3);
403 base[0].set(vec[1].x(),vec[1].y(),vec[1].z()
404 base[1].set(vec[2].x(),vec[2].y(),vec[2].z()
405 base[2].set(vec[3].x(),vec[3].y(),vec[3].z()
406
407 std::vector<const G4ThreeVectorList *> polyg
408 polygons[0] = &anchor;
409 polygons[1] = &base;
410
411 G4BoundingEnvelope benv(bmin,bmax,polygons);
412 return exists = benv.CalculateExtent(pAxis,p
413 #endif
414 }
415
416 //////////////////////////////////////////////
417 //
418 // Return whether point inside/outside/on surf
419 //
420 EInside G4Tet::Inside(const G4ThreeVector& p)
421 {
422 G4double dd[4];
423 for (G4int i = 0; i < 4; ++i) { dd[i] = fNor
424
425 G4double dist = std::max(std::max(std::max(d
426 return (dist <= -halfTolerance) ?
427 kInside : ((dist <= halfTolerance) ? kSurf
428 }
429
430 //////////////////////////////////////////////
431 //
432 // Return unit normal to surface at p
433 //
434 G4ThreeVector G4Tet::SurfaceNormal( const G4Th
435 {
436 G4double k[4];
437 for (G4int i = 0; i < 4; ++i)
438 {
439 k[i] = (G4double)(std::abs(fNormal[i].dot(
440 }
441 G4double nsurf = k[0] + k[1] + k[2] + k[3];
442 G4ThreeVector norm =
443 k[0]*fNormal[0] + k[1]*fNormal[1] + k[2]*f
444
445 if (nsurf == 1.) return norm;
446 else if (nsurf > 1.) return norm.unit(); //
447 {
448 #ifdef G4SPECSDEBUG
449 std::ostringstream message;
450 G4long oldprc = message.precision(16);
451 message << "Point p is not on surface (!?)
452 << GetName() << "\n";
453 message << "Position:\n";
454 message << " p.x() = " << p.x()/mm << "
455 message << " p.y() = " << p.y()/mm << "
456 message << " p.z() = " << p.z()/mm << "
457 G4cout.precision(oldprc);
458 G4Exception("G4Tet::SurfaceNormal(p)", "Ge
459 JustWarning, message );
460 DumpInfo();
461 #endif
462 return ApproxSurfaceNormal(p);
463 }
464 }
465
466 //////////////////////////////////////////////
467 //
468 // Find surface nearest to point and return co
469 // This method normally should not be called
470 //
471 G4ThreeVector G4Tet::ApproxSurfaceNormal(const
472 {
473 G4double dist = -DBL_MAX;
474 G4int iside = 0;
475 for (G4int i = 0; i < 4; ++i)
476 {
477 G4double d = fNormal[i].dot(p) - fDist[i];
478 if (d > dist) { dist = d; iside = i; }
479 }
480 return fNormal[iside];
481 }
482
483 //////////////////////////////////////////////
484 //
485 // Calculate distance to surface from outside,
486 // return kInfinity if no intersection
487 //
488 G4double G4Tet::DistanceToIn(const G4ThreeVect
489 const G4ThreeVect
490 {
491 G4double tin = -DBL_MAX, tout = DBL_MAX;
492 for (G4int i = 0; i < 4; ++i)
493 {
494 G4double cosa = fNormal[i].dot(v);
495 G4double dist = fNormal[i].dot(p) - fDist[
496 if (dist >= -halfTolerance)
497 {
498 if (cosa >= 0.) { return kInfinity; }
499 tin = std::max(tin, -dist/cosa);
500 }
501 else if (cosa > 0.)
502 {
503 tout = std::min(tout, -dist/cosa);
504 }
505 }
506
507 return (tout - tin <= halfTolerance) ?
508 kInfinity : ((tin < halfTolerance) ? 0. :
509 }
510
511 //////////////////////////////////////////////
512 //
513 // Estimate safety distance to surface from ou
514 //
515 G4double G4Tet::DistanceToIn(const G4ThreeVect
516 {
517 G4double dd[4];
518 for (G4int i = 0; i < 4; ++i) { dd[i] = fNor
519
520 G4double dist = std::max(std::max(std::max(d
521 return (dist > 0.) ? dist : 0.;
522 }
523
524 //////////////////////////////////////////////
525 //
526 // Calcluate distance to surface from inside
527 //
528 G4double G4Tet::DistanceToOut(const G4ThreeVec
529 const G4ThreeVec
530 const G4bool cal
531 G4bool* va
532 G4ThreeVec
533 {
534 // Calculate distances and cosines
535 G4double cosa[4], dist[4];
536 G4int ind[4] = {0}, nside = 0;
537 for (G4int i = 0; i < 4; ++i)
538 {
539 G4double tmp = fNormal[i].dot(v);
540 cosa[i] = tmp;
541 ind[nside] = (G4int)(tmp > 0) * i;
542 nside += (G4int)(tmp > 0);
543 dist[i] = fNormal[i].dot(p) - fDist[i];
544 }
545
546 // Find intersection (in most of cases nside
547 G4double tout = DBL_MAX;
548 G4int iside = 0;
549 for (G4int i = 0; i < nside; ++i)
550 {
551 G4int k = ind[i];
552 // Check: leaving the surface
553 if (dist[k] >= -halfTolerance) { tout = 0.
554 // Compute distance to intersection
555 G4double tmp = -dist[k]/cosa[k];
556 if (tmp < tout) { tout = tmp; iside = k; }
557 }
558
559 // Set normal, if required, and return dista
560 if (calcNorm)
561 {
562 *validNorm = true;
563 *n = fNormal[iside];
564 }
565 return tout;
566 }
567
568 //////////////////////////////////////////////
569 //
570 // Calculate safety distance to surface from i
571 //
572 G4double G4Tet::DistanceToOut(const G4ThreeVec
573 {
574 G4double dd[4];
575 for (G4int i = 0; i < 4; ++i) { dd[i] = fDis
576
577 G4double dist = std::min(std::min(std::min(d
578 return (dist > 0.) ? dist : 0.;
579 }
580
581 //////////////////////////////////////////////
582 //
583 // GetEntityType
584 //
585 G4GeometryType G4Tet::GetEntityType() const
586 {
587 return {"G4Tet"};
588 }
589
590 //////////////////////////////////////////////
591 //
592 // IsFaceted
593 //
594 G4bool G4Tet::IsFaceted() const
595 {
596 return true;
597 }
598
599 //////////////////////////////////////////////
600 //
601 // Make a clone of the object
602 //
603 G4VSolid* G4Tet::Clone() const
604 {
605 return new G4Tet(*this);
606 }
607
608 //////////////////////////////////////////////
609 //
610 // Stream object contents to an output stream
611 //
612 std::ostream& G4Tet::StreamInfo(std::ostream&
613 {
614 G4long oldprc = os.precision(16);
615 os << "-------------------------------------
616 << " *** Dump for solid - " << GetName
617 << " =================================
618 << " Solid type: " << GetEntityType() <<
619 << " Parameters: \n"
620 << " anchor: " << fVertex[0]/mm << " m
621 << " p1 : " << fVertex[1]/mm << " m
622 << " p2 : " << fVertex[2]/mm << " m
623 << " p3 : " << fVertex[3]/mm << " m
624 << "-------------------------------------
625 os.precision(oldprc);
626 return os;
627 }
628
629 //////////////////////////////////////////////
630 //
631 // Return random point on the surface
632 //
633 G4ThreeVector G4Tet::GetPointOnSurface() const
634 {
635 constexpr G4int iface[4][3] = { {0,1,2}, {0,
636
637 // Select face
638 G4double select = fSurfaceArea*G4QuickRand()
639 G4int i = 0;
640 i += (G4int)(select > fArea[0]);
641 i += (G4int)(select > fArea[0] + fArea[1]);
642 i += (G4int)(select > fArea[0] + fArea[1] +
643
644 // Set selected triangle
645 G4ThreeVector p0 = fVertex[iface[i][0]];
646 G4ThreeVector e1 = fVertex[iface[i][1]] - p0
647 G4ThreeVector e2 = fVertex[iface[i][2]] - p0
648
649 // Return random point
650 G4double r1 = G4QuickRand();
651 G4double r2 = G4QuickRand();
652 return (r1 + r2 > 1.) ?
653 p0 + e1*(1. - r1) + e2*(1. - r2) : p0 + e1
654 }
655
656 //////////////////////////////////////////////
657 //
658 // Return volume of the tetrahedron
659 //
660 G4double G4Tet::GetCubicVolume()
661 {
662 return fCubicVolume;
663 }
664
665 //////////////////////////////////////////////
666 //
667 // Return surface area of the tetrahedron
668 //
669 G4double G4Tet::GetSurfaceArea()
670 {
671 return fSurfaceArea;
672 }
673
674 //////////////////////////////////////////////
675 //
676 // Methods for visualisation
677 //
678 void G4Tet::DescribeYourselfTo (G4VGraphicsSce
679 {
680 scene.AddSolid (*this);
681 }
682
683 //////////////////////////////////////////////
684 //
685 // Return VisExtent
686 //
687 G4VisExtent G4Tet::GetExtent() const
688 {
689 return { fBmin.x(), fBmax.x(),
690 fBmin.y(), fBmax.y(),
691 fBmin.z(), fBmax.z() };
692 }
693
694 //////////////////////////////////////////////
695 //
696 // CreatePolyhedron
697 //
698 G4Polyhedron* G4Tet::CreatePolyhedron() const
699 {
700 // Check orientation of vertices
701 G4ThreeVector v1 = fVertex[1] - fVertex[0];
702 G4ThreeVector v2 = fVertex[2] - fVertex[0];
703 G4ThreeVector v3 = fVertex[3] - fVertex[0];
704 G4bool invert = v1.cross(v2).dot(v3) < 0.;
705 G4int k2 = (invert) ? 3 : 2;
706 G4int k3 = (invert) ? 2 : 3;
707
708 // Set coordinates of vertices
709 G4double xyz[4][3];
710 for (G4int i = 0; i < 3; ++i)
711 {
712 xyz[0][i] = fVertex[0][i];
713 xyz[1][i] = fVertex[1][i];
714 xyz[2][i] = fVertex[k2][i];
715 xyz[3][i] = fVertex[k3][i];
716 }
717
718 // Create polyhedron
719 G4int faces[4][4] = { {1,3,2,0}, {1,4,3,0},
720 auto ph = new G4Polyhedron;
721 ph->createPolyhedron(4,4,xyz,faces);
722
723 return ph;
724 }
725
726 //////////////////////////////////////////////
727 //
728 // GetPolyhedron
729 //
730 G4Polyhedron* G4Tet::GetPolyhedron() const
731 {
732 if (fpPolyhedron == nullptr ||
733 fRebuildPolyhedron ||
734 fpPolyhedron->GetNumberOfRotationStepsAt
735 fpPolyhedron->GetNumberOfRotationSteps()
736 {
737 G4AutoLock l(&polyhedronMutex);
738 delete fpPolyhedron;
739 fpPolyhedron = CreatePolyhedron();
740 fRebuildPolyhedron = false;
741 l.unlock();
742 }
743 return fpPolyhedron;
744 }
745
746 #endif
747