| 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 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
61 ////////////////////////////////////////////////////////////////////////
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 != nullptr)
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
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 data and allocates memory
102 // for usage restricted to object persistency.
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] = 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
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] = 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
167 ////////////////////////////////////////////////////////////////////////
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 //
173 G4bool G4Tet::CheckDegeneracy(const G4ThreeVector& p0,
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
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] = 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
246 ////////////////////////////////////////////////////////////////////////
247 //
248 // Set vertices
249 //
250 void G4Tet::SetVertices(const G4ThreeVector& p0,
251 const G4ThreeVector& p1,
252 const G4ThreeVector& p2,
253 const G4ThreeVector& p3, G4bool* degeneracyFlag)
254 {
255 // Check for degeneracy
256 G4bool degenerate = CheckDegeneracy(p0, p1, p2, p3);
257 if (degeneracyFlag != nullptr)
258 {
259 *degeneracyFlag = degenerate;
260 }
261 else if (degenerate)
262 {
263 std::ostringstream message;
264 message << "Degenerate tetrahedron is not permitted: " << GetName() << " !\n"
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 - p0).dot(p3 - p0))/6.;
271 G4Exception("G4Tet::SetVertices()", "GeomSolids0002",
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) const
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() const
302 {
303 std::vector<G4ThreeVector> vertices(4);
304 for (G4int i = 0; i < 4; ++i) { vertices[i] = fVertex[i]; }
305 return vertices;
306 }
307
308 ////////////////////////////////////////////////////////////////////////
309 //
310 // Dispatch to parameterisation for replication mechanism dimension
311 // computation & modification.
312 //
313 void G4Tet::ComputeDimensions(G4VPVParameterisation* ,
314 const G4int ,
315 const G4VPhysicalVolume* )
316 {
317 }
318
319 ////////////////////////////////////////////////////////////////////////
320 //
321 // Set bounding box
322 //
323 void G4Tet::SetBoundingLimits(const G4ThreeVector& pMin,
324 const G4ThreeVector& pMax)
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] != 0)
337 {
338 std::ostringstream message;
339 message << "Attempt to set bounding box that does not encapsulate solid: "
340 << GetName() << " !\n"
341 << " Specified bounding box limits:\n"
342 << " pmin: " << pMin << "\n"
343 << " pmax: " << pMax << "\n"
344 << " Tetrahedron vertices:\n"
345 << " anchor " << fVertex[0] << ((iout[0]) != 0 ? " is outside\n" : "\n")
346 << " p1 " << fVertex[1] << ((iout[1]) != 0 ? " is outside\n" : "\n")
347 << " p2 " << fVertex[2] << ((iout[2]) != 0 ? " is outside\n" : "\n")
348 << " p3 " << fVertex[3] << ((iout[3]) != 0 ? " is outside" : "");
349 G4Exception("G4Tet::SetBoundingLimits()", "GeomSolids0002",
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, G4ThreeVector& pMax) const
361 {
362 pMin = fBmin;
363 pMax = fBmax;
364 }
365
366 ////////////////////////////////////////////////////////////////////////
367 //
368 // Calculate extent under transform and specified limit
369 //
370 G4bool G4Tet::CalculateExtent(const EAxis pAxis,
371 const G4VoxelLimits& pVoxelLimit,
372 const G4AffineTransform& pTransform,
373 G4double& pMin, G4double& pMax) const
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 case of complex 3D meshes
383 //
384 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
385
386 #if 0
387 // Precise extent computation (disabled by default for this shape)
388 //
389 G4bool exist;
390 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
391 {
392 return exist = (pMin < pMax) ? true : false;
393 }
394
395 // Set bounding envelope (benv) and calculate extent
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 *> polygons(2);
408 polygons[0] = &anchor;
409 polygons[1] = &base;
410
411 G4BoundingEnvelope benv(bmin,bmax,polygons);
412 return exists = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
413 #endif
414 }
415
416 ////////////////////////////////////////////////////////////////////////
417 //
418 // Return whether point inside/outside/on surface
419 //
420 EInside G4Tet::Inside(const G4ThreeVector& p) const
421 {
422 G4double dd[4];
423 for (G4int i = 0; i < 4; ++i) { dd[i] = fNormal[i].dot(p) - fDist[i]; }
424
425 G4double dist = std::max(std::max(std::max(dd[0], dd[1]), dd[2]), dd[3]);
426 return (dist <= -halfTolerance) ?
427 kInside : ((dist <= halfTolerance) ? kSurface : kOutside);
428 }
429
430 ////////////////////////////////////////////////////////////////////////
431 //
432 // Return unit normal to surface at p
433 //
434 G4ThreeVector G4Tet::SurfaceNormal( const G4ThreeVector& p) const
435 {
436 G4double k[4];
437 for (G4int i = 0; i < 4; ++i)
438 {
439 k[i] = (G4double)(std::abs(fNormal[i].dot(p) - fDist[i]) <= halfTolerance);
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]*fNormal[2] + k[3]*fNormal[3];
444
445 if (nsurf == 1.) return norm;
446 else if (nsurf > 1.) return norm.unit(); // edge or vertex
447 {
448 #ifdef G4SPECSDEBUG
449 std::ostringstream message;
450 G4long oldprc = message.precision(16);
451 message << "Point p is not on surface (!?) of solid: "
452 << GetName() << "\n";
453 message << "Position:\n";
454 message << " p.x() = " << p.x()/mm << " mm\n";
455 message << " p.y() = " << p.y()/mm << " mm\n";
456 message << " p.z() = " << p.z()/mm << " mm";
457 G4cout.precision(oldprc);
458 G4Exception("G4Tet::SurfaceNormal(p)", "GeomSolids1002",
459 JustWarning, message );
460 DumpInfo();
461 #endif
462 return ApproxSurfaceNormal(p);
463 }
464 }
465
466 ////////////////////////////////////////////////////////////////////////
467 //
468 // Find surface nearest to point and return corresponding normal
469 // This method normally should not be called
470 //
471 G4ThreeVector G4Tet::ApproxSurfaceNormal(const G4ThreeVector& p) 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 G4ThreeVector& p,
489 const G4ThreeVector& v) const
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[i];
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. : tin);
509 }
510
511 ////////////////////////////////////////////////////////////////////////
512 //
513 // Estimate safety distance to surface from outside
514 //
515 G4double G4Tet::DistanceToIn(const G4ThreeVector& p) const
516 {
517 G4double dd[4];
518 for (G4int i = 0; i < 4; ++i) { dd[i] = fNormal[i].dot(p) - fDist[i]; }
519
520 G4double dist = std::max(std::max(std::max(dd[0], dd[1]), dd[2]), dd[3]);
521 return (dist > 0.) ? dist : 0.;
522 }
523
524 ////////////////////////////////////////////////////////////////////////
525 //
526 // Calcluate distance to surface from inside
527 //
528 G4double G4Tet::DistanceToOut(const G4ThreeVector& p,
529 const G4ThreeVector& v,
530 const G4bool calcNorm,
531 G4bool* validNorm,
532 G4ThreeVector* n) const
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 == 1)
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.; iside = k; break; }
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 distance to out
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 inside
571 //
572 G4double G4Tet::DistanceToOut(const G4ThreeVector& p) const
573 {
574 G4double dd[4];
575 for (G4int i = 0; i < 4; ++i) { dd[i] = fDist[i] - fNormal[i].dot(p); }
576
577 G4double dist = std::min(std::min(std::min(dd[0], dd[1]), dd[2]), dd[3]);
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& os) const
613 {
614 G4long oldprc = os.precision(16);
615 os << "-----------------------------------------------------------\n"
616 << " *** Dump for solid - " << GetName() << " ***\n"
617 << " ===================================================\n"
618 << " Solid type: " << GetEntityType() << "\n"
619 << " Parameters: \n"
620 << " anchor: " << fVertex[0]/mm << " mm\n"
621 << " p1 : " << fVertex[1]/mm << " mm\n"
622 << " p2 : " << fVertex[2]/mm << " mm\n"
623 << " p3 : " << fVertex[3]/mm << " mm\n"
624 << "-----------------------------------------------------------\n";
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,2,3}, {0,3,1}, {1,2,3} };
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] + fArea[2]);
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*r1 + e2*r2;
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 (G4VGraphicsScene& scene) const
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}, {1,2,4,0}, {2,3,4,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->GetNumberOfRotationStepsAtTimeOfCreation() !=
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