ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Tet.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Tet.cc
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 intellectual property of the *
19 // * Vanderbilt University Free Electron Laser Center *
20 // * Vanderbilt University, Nashville, TN, USA *
21 // * Development supported by: *
22 // * United States MFEL program under grant FA9550-04-1-0045 *
23 // * and NASA under contract number NNG04CT05P *
24 // * Written by Marcus H. Mendenhall and Robert A. Weller. *
25 // * *
26 // * Contributed to the Geant4 Core, January, 2005. *
27 // * *
28 // ********************************************************************
29 //
30 // Implementation for G4Tet class
31 //
32 // 03.09.2004 - Marcus Mendenhall, created
33 // 08.01.2020 - Evgueni Tcherniaev, complete revision, speed up
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_INITIALIZER;
57 }
58 
59 using namespace CLHEP;
60 
62 //
63 // Constructor - create a tetrahedron
64 // A Tet has all of its geometrical information precomputed
65 //
66 G4Tet::G4Tet(const G4String& pName,
67  const G4ThreeVector& p0,
68  const G4ThreeVector& p1,
69  const G4ThreeVector& p2,
70  const G4ThreeVector& p3, G4bool* degeneracyFlag)
71  : G4VSolid(pName)
72 {
73  // Check for degeneracy
74  G4bool degenerate = CheckDegeneracy(p0, p1, p2, p3);
75  if (degeneracyFlag)
76  {
77  *degeneracyFlag = degenerate;
78  }
79  else if (degenerate)
80  {
81  std::ostringstream message;
82  message << "Degenerate tetrahedron: " << GetName() << " !\n"
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 - p0).dot(p3 - p0))/6.;
89  G4Exception("G4Tet::G4Tet()", "GeomSolids0002", FatalException, message);
90  }
91 
92  // Define surface thickness
94 
95  // Set data members
96  Initialize(p0, p1, p2, p3);
97 }
98 
100 //
101 // Fake default constructor - sets only member data and allocates memory
102 // for usage restricted to object persistency.
103 //
104 G4Tet::G4Tet( __void__& a )
105  : G4VSolid(a)
106 {
107 }
108 
110 //
111 // Destructor
112 //
114 {
115  delete fpPolyhedron; fpPolyhedron = nullptr;
116 }
117 
119 //
120 // Copy constructor
121 //
122 G4Tet::G4Tet(const G4Tet& rhs)
123  : G4VSolid(rhs)
124 {
128  for (G4int i = 0; i < 4; ++i) { fVertex[i] = rhs.fVertex[i]; }
129  for (G4int i = 0; i < 4; ++i) { fNormal[i] = rhs.fNormal[i]; }
130  for (G4int i = 0; i < 4; ++i) { fDist[i] = rhs.fDist[i]; }
131  for (G4int i = 0; i < 4; ++i) { fArea[i] = rhs.fArea[i]; }
132  fBmin = rhs.fBmin;
133  fBmax = rhs.fBmax;
134 }
135 
137 //
138 // Assignment operator
139 //
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  //
155  for (G4int i = 0; i < 4; ++i) { fVertex[i] = rhs.fVertex[i]; }
156  for (G4int i = 0; i < 4; ++i) { fNormal[i] = rhs.fNormal[i]; }
157  for (G4int i = 0; i < 4; ++i) { fDist[i] = rhs.fDist[i]; }
158  for (G4int i = 0; i < 4; ++i) { fArea[i] = rhs.fArea[i]; }
159  fBmin = rhs.fBmin;
160  fBmax = rhs.fBmax;
161  fRebuildPolyhedron = false;
162  delete fpPolyhedron; fpPolyhedron = nullptr;
163 
164  return *this;
165 }
166 
168 //
169 // Return true if tetrahedron is degenerate
170 // Tetrahedron is concidered as degenerate in case if its minimal
171 // height is less than degeneracy tolerance
172 //
174  const G4ThreeVector& p1,
175  const G4ThreeVector& p2,
176  const G4ThreeVector& p3) const
177 {
178  G4double hmin = 4. * kCarTolerance; // degeneracy tolerance
179 
180  // Calculate volume
181  G4double vol = std::abs((p1 - p0).cross(p2 - p0).dot(p3 - p0));
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] > ss[k]) k = i; }
193 
194  // Check: vol^2 / s^2 <= hmin^2
195  return (vol*vol <= ss[k]*hmin*hmin);
196 }
197 
199 //
200 // Set data members
201 //
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 (G4int i = 0; i < 4; ++i) { norm[i] = -norm[i]; }
222  }
223 
224  // Set normals to face planes
225  for (G4int i = 0; i < 4; ++i) { fNormal[i] = norm[i].unit(); }
226 
227  // Set distances to planes
228  for (G4int i = 0; i < 3; ++i) { fDist[i] = fNormal[i].dot(p0); }
229  fDist[3] = fNormal[3].dot(p1);
230 
231  // Set face areas
232  for (G4int i = 0; i < 4; ++i) { fArea[i] = 0.5*norm[i].mag(); }
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], p1[i]), p2[i]), p3[i]);
238  fBmax[i] = std::max(std::max(std::max(p0[i], p1[i]), p2[i]), p3[i]);
239  }
240 
241  // Set volume and surface area
242  fCubicVolume = std::abs(volume)/6.;
243  fSurfaceArea = fArea[0] + fArea[1] + fArea[2] + fArea[3];
244 }
245 
247 //
248 // Set vertices
249 //
251  const G4ThreeVector& p1,
252  const G4ThreeVector& p2,
253  const G4ThreeVector& p3)
254 {
255  // Check for degeneracy
256  G4bool degenerate = CheckDegeneracy(p0, p1, p2, p3);
257  if (degenerate)
258  {
259  std::ostringstream message;
260  message << "Degenerate tetrahedron is not permitted: " << GetName() << " !\n"
261  << " anchor: " << p0 << "\n"
262  << " p1 : " << p1 << "\n"
263  << " p2 : " << p2 << "\n"
264  << " p3 : " << p3 << "\n"
265  << " volume: "
266  << std::abs((p1 - p0).cross(p2 - p0).dot(p3 - p0))/6.;
267  G4Exception("G4Tet::G4SetVertices()", "GeomSolids0002",
268  FatalException, message);
269  }
270 
271  // Set data members
272  Initialize(p0, p1, p2, p3);
273 
274  // Set flag to rebuild polyhedron
275  fRebuildPolyhedron = true;
276 }
277 
279 //
280 // Return four vertices
281 //
283  G4ThreeVector& p1,
284  G4ThreeVector& p2,
285  G4ThreeVector& p3) const
286 {
287  p0 = fVertex[0];
288  p1 = fVertex[1];
289  p2 = fVertex[2];
290  p3 = fVertex[3];
291 }
292 
294 //
295 // Return std::vector of vertices
296 //
297 std::vector<G4ThreeVector> G4Tet::GetVertices() const
298 {
299  std::vector<G4ThreeVector> vertices(4);
300  for (G4int i = 0; i < 4; ++i) { vertices[i] = fVertex[i]; }
301  return vertices;
302 }
303 
305 //
306 // Dispatch to parameterisation for replication mechanism dimension
307 // computation & modification.
308 //
310  const G4int ,
311  const G4VPhysicalVolume* )
312 {
313 }
314 
316 //
317 // Return bounding box
318 //
320 {
321  pMin = fBmin;
322  pMax = fBmax;
323 }
324 
326 //
327 // Calculate extent under transform and specified limit
328 //
330  const G4VoxelLimits& pVoxelLimit,
331  const G4AffineTransform& pTransform,
332  G4double& pMin, G4double& pMax) const
333 {
334  G4ThreeVector bmin, bmax;
335 
336  // Check bounding box (bbox)
337  //
338  BoundingLimits(bmin,bmax);
339  G4BoundingEnvelope bbox(bmin,bmax);
340 
341  // Use simple bounding-box to help in the case of complex 3D meshes
342  //
343  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
344 
345 #if 0
346  // Precise extent computation (disabled by default for this shape)
347  //
348  G4bool exist;
349  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
350  {
351  return exist = (pMin < pMax) ? true : false;
352  }
353 
354  // Set bounding envelope (benv) and calculate extent
355  //
356  std::vector<G4ThreeVector> vec = GetVertices();
357 
358  G4ThreeVectorList anchor(1);
359  anchor[0].set(vec[0].x(),vec[0].y(),vec[0].z());
360 
362  base[0].set(vec[1].x(),vec[1].y(),vec[1].z());
363  base[1].set(vec[2].x(),vec[2].y(),vec[2].z());
364  base[2].set(vec[3].x(),vec[3].y(),vec[3].z());
365 
366  std::vector<const G4ThreeVectorList *> polygons(2);
367  polygons[0] = &anchor;
368  polygons[1] = &base;
369 
370  G4BoundingEnvelope benv(bmin,bmax,polygons);
371  return exists = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
372 #endif
373 }
374 
376 //
377 // Return whether point inside/outside/on surface
378 //
380 {
381  G4double dd[4];
382  for (G4int i = 0; i < 4; ++i) { dd[i] = fNormal[i].dot(p) - fDist[i]; }
383 
384  G4double dist = std::max(std::max(std::max(dd[0], dd[1]), dd[2]), dd[3]);
385  return (dist <= -halfTolerance) ?
386  kInside : ((dist <= halfTolerance) ? kSurface : kOutside);
387 }
388 
390 //
391 // Return unit normal to surface at p
392 //
394 {
395  G4double k[4];
396  for (G4int i = 0; i < 4; ++i)
397  {
398  k[i] = std::abs(fNormal[i].dot(p) - fDist[i]) <= halfTolerance;
399  }
400  G4double nsurf = k[0] + k[1] + k[2] + k[3];
402  k[0]*fNormal[0] + k[1]*fNormal[1] + k[2]*fNormal[2] + k[3]*fNormal[3];
403 
404  if (nsurf == 1.) return norm;
405  else if (nsurf > 1.) return norm.unit(); // edge or vertex
406  {
407 #ifdef G4SPECSDEBUG
408  std::ostringstream message;
409  G4int oldprc = message.precision(16);
410  message << "Point p is not on surface (!?) of solid: "
411  << GetName() << "\n";
412  message << "Position:\n";
413  message << " p.x() = " << p.x()/mm << " mm\n";
414  message << " p.y() = " << p.y()/mm << " mm\n";
415  message << " p.z() = " << p.z()/mm << " mm";
416  G4cout.precision(oldprc);
417  G4Exception("G4Tet::SurfaceNormal(p)", "GeomSolids1002",
418  JustWarning, message );
419  DumpInfo();
420 #endif
421  return ApproxSurfaceNormal(p);
422  }
423 }
424 
426 //
427 // Find surface nearest to point and return corresponding normal
428 // This method normally should not be called
429 //
431 {
432  G4double dist = -DBL_MAX;
433  G4int iside = 0;
434  for (G4int i = 0; i < 4; ++i)
435  {
436  G4double d = fNormal[i].dot(p) - fDist[i];
437  if (d > dist) { dist = d; iside = i; }
438  }
439  return fNormal[iside];
440 }
441 
443 //
444 // Calculate distance to surface from outside,
445 // return kInfinity if no intersection
446 //
448  const G4ThreeVector& v) const
449 {
450  G4double tin = -DBL_MAX, tout = DBL_MAX;
451  for (G4int i = 0; i < 4; ++i)
452  {
453  G4double cosa = fNormal[i].dot(v);
454  G4double dist = fNormal[i].dot(p) - fDist[i];
455  if (dist >= -halfTolerance)
456  {
457  if (cosa >= 0.) { return kInfinity; }
458  tin = std::max(tin, -dist/cosa);
459  }
460  else if (cosa > 0.)
461  {
462  tout = std::min(tout, -dist/cosa);
463  }
464  }
465 
466  return (tout - tin <= halfTolerance) ?
467  kInfinity : ((tin < halfTolerance) ? 0. : tin);
468 }
469 
471 //
472 // Estimate safety distance to surface from outside
473 //
475 {
476  G4double dd[4];
477  for (G4int i = 0; i < 4; ++i) { dd[i] = fNormal[i].dot(p) - fDist[i]; }
478 
479  G4double dist = std::max(std::max(std::max(dd[0], dd[1]), dd[2]), dd[3]);
480  return (dist > 0.) ? dist : 0.;
481 }
482 
484 //
485 // Calcluate distance to surface from inside
486 //
488  const G4ThreeVector& v,
489  const G4bool calcNorm,
490  G4bool* validNorm,
491  G4ThreeVector* n) const
492 {
493  // Calculate distances and cosines
494  G4double cosa[4], dist[4];
495  G4int ind[4] = {0}, nside = 0;
496  for (G4int i = 0; i < 4; ++i)
497  {
498  G4double tmp = fNormal[i].dot(v);
499  cosa[i] = tmp;
500  ind[nside] = (tmp > 0) * i;
501  nside += (tmp > 0);
502  dist[i] = fNormal[i].dot(p) - fDist[i];
503  }
504 
505  // Find intersection (in most of cases nside == 1)
506  G4double tout = DBL_MAX;
507  G4int iside = 0;
508  for (G4int i = 0; i < nside; ++i)
509  {
510  G4int k = ind[i];
511  // Check: leaving the surface
512  if (dist[k] >= -halfTolerance) { tout = 0.; iside = k; break; }
513  // Compute distance to intersection
514  G4double tmp = -dist[k]/cosa[k];
515  if (tmp < tout) { tout = tmp; iside = k; }
516  }
517 
518  // Set normal, if required, and return distance to out
519  if (calcNorm)
520  {
521  *validNorm = true;
522  *n = fNormal[iside];
523  }
524  return tout;
525 }
526 
528 //
529 // Calculate safety distance to surface from inside
530 //
532 {
533  G4double dd[4];
534  for (G4int i = 0; i < 4; ++i) { dd[i] = fDist[i] - fNormal[i].dot(p); }
535 
536  G4double dist = std::min(std::min(std::min(dd[0], dd[1]), dd[2]), dd[3]);
537  return (dist > 0.) ? dist : 0.;
538 }
539 
541 //
542 // GetEntityType
543 //
545 {
546  return G4String("G4Tet");
547 }
548 
550 //
551 // Make a clone of the object
552 //
554 {
555  return new G4Tet(*this);
556 }
557 
559 //
560 // Stream object contents to an output stream
561 //
562 std::ostream& G4Tet::StreamInfo(std::ostream& os) const
563 {
564  G4int oldprc = os.precision(16);
565  os << "-----------------------------------------------------------\n"
566  << " *** Dump for solid - " << GetName() << " ***\n"
567  << " ===================================================\n"
568  << " Solid type: " << GetEntityType() << "\n"
569  << " Parameters: \n"
570  << " anchor: " << fVertex[0]/mm << " mm\n"
571  << " p1 : " << fVertex[1]/mm << " mm\n"
572  << " p2 : " << fVertex[2]/mm << " mm\n"
573  << " p3 : " << fVertex[3]/mm << " mm\n"
574  << "-----------------------------------------------------------\n";
575  os.precision(oldprc);
576  return os;
577 }
578 
580 //
581 // Return random point on the surface
582 //
584 {
585  constexpr G4int iface[4][3] = { {0,1,2}, {0,2,3}, {0,3,1}, {1,2,3} };
586 
587  // Select face
588  G4double select = fSurfaceArea*G4QuickRand();
589  G4int i = 0;
590  for ( ; i < 4; ++i) { if ((select -= fArea[i]) <= 0.) break; }
591 
592  // Set selected triangle
593  G4ThreeVector p0 = fVertex[iface[i][0]];
594  G4ThreeVector e1 = fVertex[iface[i][1]] - p0;
595  G4ThreeVector e2 = fVertex[iface[i][2]] - p0;
596 
597  // Return random point
598  G4double r1 = G4QuickRand();
599  G4double r2 = G4QuickRand();
600  return (r1 + r2 > 1.) ?
601  p0 + e1*(1. - r1) + e2*(1. - r2) : p0 + e1*r1 + e2*r2;
602 }
603 
605 //
606 // Return volume of the tetrahedron
607 //
609 {
610  return fCubicVolume;
611 }
612 
614 //
615 // Return surface area of the tetrahedron
616 //
618 {
619  return fSurfaceArea;
620 }
621 
623 //
624 // Methods for visualisation
625 //
627 {
628  scene.AddSolid (*this);
629 }
630 
632 //
633 // Return VisExtent
634 //
636 {
637  return G4VisExtent(fBmin.x(), fBmax.x(),
638  fBmin.y(), fBmax.y(),
639  fBmin.z(), fBmax.z());
640 }
641 
643 //
644 // CreatePolyhedron
645 //
647 {
648  // Check orientation of vertices
649  G4ThreeVector v1 = fVertex[1] - fVertex[0];
650  G4ThreeVector v2 = fVertex[2] - fVertex[0];
651  G4ThreeVector v3 = fVertex[3] - fVertex[0];
652  G4bool invert = v1.cross(v2).dot(v3) < 0.;
653  G4int k2 = (invert) ? 3 : 2;
654  G4int k3 = (invert) ? 2 : 3;
655 
656  // Set coordinates of vertices
657  G4double xyz[4][3];
658  for (G4int i = 0; i < 3; ++i)
659  {
660  xyz[0][i] = fVertex[0][i];
661  xyz[1][i] = fVertex[1][i];
662  xyz[2][i] = fVertex[k2][i];
663  xyz[3][i] = fVertex[k3][i];
664  }
665 
666  // Create polyhedron
667  G4int faces[4][4] = { {1,3,2,0}, {1,4,3,0}, {1,2,4,0}, {2,3,4,0} };
668  G4Polyhedron* ph = new G4Polyhedron;
669  ph->createPolyhedron(4,4,xyz,faces);
670 
671  return ph;
672 }
673 
675 //
676 // GetPolyhedron
677 //
679 {
680  if (fpPolyhedron == nullptr ||
684  {
685  G4AutoLock l(&polyhedronMutex);
686  delete fpPolyhedron;
688  fRebuildPolyhedron = false;
689  l.unlock();
690  }
691  return fpPolyhedron;
692 }
693 
694 #endif