Geant4  10.02.p02
G4Tet.cc
Go to the documentation of this file.
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 // $Id: G4Tet.cc 83572 2014-09-01 15:23:27Z gcosmo $
31 //
32 // class G4Tet
33 //
34 // Implementation for G4Tet class
35 //
36 // History:
37 //
38 // 20040903 - Marcus Mendenhall, created G4Tet
39 // 20041101 - Marcus Mendenhall, optimized constant dot products with
40 // fCdotNijk values
41 // 20041101 - MHM removed tracking error by clipping DistanceToOut to 0
42 // for surface cases
43 // 20041101 - MHM many speed optimizations in if statements
44 // 20041101 - MHM changed vdotn comparisons to 1e-12 instead of 0.0 to
45 // avoid nearly-parallel problems
46 // 20041102 - MHM Added extra distance into solid to DistanceToIn(p,v)
47 // hit testing
48 // 20041102 - MHM added ability to check for degeneracy without throwing
49 // G4Exception
50 // 20041103 - MHM removed many unused variables from class
51 // 20040803 - Dionysios Anninos, added GetPointOnSurface() method
52 // 20061112 - MHM added code for G4VSolid GetSurfaceArea()
53 // 20100920 - Gabriele Cosmo added copy-ctor and operator=()
54 //
55 // --------------------------------------------------------------------
56 
57 #include "G4Tet.hh"
58 
59 #if !defined(G4GEOM_USE_UTET)
60 
61 const char G4Tet::CVSVers[]="$Id: G4Tet.cc 83572 2014-09-01 15:23:27Z gcosmo $";
62 
63 #include "G4VoxelLimits.hh"
64 #include "G4AffineTransform.hh"
65 
66 #include "G4VPVParameterisation.hh"
67 
68 #include "Randomize.hh"
69 
70 #include "G4VGraphicsScene.hh"
71 #include "G4Polyhedron.hh"
72 #include "G4VisExtent.hh"
73 
74 #include "G4ThreeVector.hh"
75 
76 #include <cmath>
77 
78 #include "G4AutoLock.hh"
79 
80 namespace
81 {
82  G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
83 }
84 
85 using namespace CLHEP;
86 
88 //
89 // Constructor - create a tetrahedron
90 // This class is implemented separately from general polyhedra,
91 // because the simplex geometry can be computed very quickly,
92 // which may become important in situations imported from mesh generators,
93 // in which a very large number of G4Tets are created.
94 // A Tet has all of its geometrical information precomputed
95 
96 G4Tet::G4Tet(const G4String& pName,
97  G4ThreeVector anchor,
98  G4ThreeVector p2,
99  G4ThreeVector p3,
100  G4ThreeVector p4, G4bool *degeneracyFlag)
101  : G4VSolid(pName), fRebuildPolyhedron(false), fpPolyhedron(0), warningFlag(0)
102 {
103  // fV<x><y> is vector from vertex <y> to vertex <x>
104  //
105  G4ThreeVector fV21=p2-anchor;
106  G4ThreeVector fV31=p3-anchor;
107  G4ThreeVector fV41=p4-anchor;
108 
109  // make sure this is a correctly oriented set of points for the tetrahedron
110  //
111  G4double signed_vol=fV21.cross(fV31).dot(fV41);
112  if(signed_vol<0.0)
113  {
114  G4ThreeVector temp(p4);
115  p4=p3;
116  p3=temp;
117  temp=fV41;
118  fV41=fV31;
119  fV31=temp;
120  }
121  fCubicVolume = std::fabs(signed_vol) / 6.;
122 
123  G4ThreeVector fV24=p2-p4;
124  G4ThreeVector fV43=p4-p3;
125  G4ThreeVector fV32=p3-p2;
126 
127  fXMin=std::min(std::min(std::min(anchor.x(), p2.x()),p3.x()),p4.x());
128  fXMax=std::max(std::max(std::max(anchor.x(), p2.x()),p3.x()),p4.x());
129  fYMin=std::min(std::min(std::min(anchor.y(), p2.y()),p3.y()),p4.y());
130  fYMax=std::max(std::max(std::max(anchor.y(), p2.y()),p3.y()),p4.y());
131  fZMin=std::min(std::min(std::min(anchor.z(), p2.z()),p3.z()),p4.z());
132  fZMax=std::max(std::max(std::max(anchor.z(), p2.z()),p3.z()),p4.z());
133 
134  fDx=(fXMax-fXMin)*0.5; fDy=(fYMax-fYMin)*0.5; fDz=(fZMax-fZMin)*0.5;
135 
136  fMiddle=G4ThreeVector(fXMax+fXMin, fYMax+fYMin, fZMax+fZMin)*0.5;
137  fMaxSize=std::max(std::max(std::max((anchor-fMiddle).mag(),
138  (p2-fMiddle).mag()),
139  (p3-fMiddle).mag()),
140  (p4-fMiddle).mag());
141 
142  G4bool degenerate=std::fabs(signed_vol) < 1e-9*fMaxSize*fMaxSize*fMaxSize;
143 
144  if(degeneracyFlag) *degeneracyFlag=degenerate;
145  else if (degenerate)
146  {
147  G4Exception("G4Tet::G4Tet()", "GeomSolids0002", FatalException,
148  "Degenerate tetrahedron not allowed.");
149  }
150 
151  fTol=1e-9*(std::fabs(fXMin)+std::fabs(fXMax)+std::fabs(fYMin)
152  +std::fabs(fYMax)+std::fabs(fZMin)+std::fabs(fZMax));
153  //fTol=kCarTolerance;
154 
155  fAnchor=anchor;
156  fP2=p2;
157  fP3=p3;
158  fP4=p4;
159 
160  G4ThreeVector fCenter123=(anchor+p2+p3)*(1.0/3.0); // face center
161  G4ThreeVector fCenter134=(anchor+p4+p3)*(1.0/3.0);
162  G4ThreeVector fCenter142=(anchor+p4+p2)*(1.0/3.0);
163  G4ThreeVector fCenter234=(p2+p3+p4)*(1.0/3.0);
164 
165  // compute area of each triangular face by cross product
166  // and sum for total surface area
167 
168  G4ThreeVector normal123=fV31.cross(fV21);
169  G4ThreeVector normal134=fV41.cross(fV31);
170  G4ThreeVector normal142=fV21.cross(fV41);
171  G4ThreeVector normal234=fV32.cross(fV43);
172 
173  fSurfaceArea=(
174  normal123.mag()+
175  normal134.mag()+
176  normal142.mag()+
177  normal234.mag()
178  )/2.0;
179 
180  fNormal123=normal123.unit();
181  fNormal134=normal134.unit();
182  fNormal142=normal142.unit();
183  fNormal234=normal234.unit();
184 
185  fCdotN123=fCenter123.dot(fNormal123);
186  fCdotN134=fCenter134.dot(fNormal134);
187  fCdotN142=fCenter142.dot(fNormal142);
188  fCdotN234=fCenter234.dot(fNormal234);
189 }
190 
192 //
193 // Fake default constructor - sets only member data and allocates memory
194 // for usage restricted to object persistency.
195 //
196 G4Tet::G4Tet( __void__& a )
197  : G4VSolid(a), fCubicVolume(0.), fSurfaceArea(0.),
198  fRebuildPolyhedron(false), fpPolyhedron(0),
199  fAnchor(0,0,0), fP2(0,0,0), fP3(0,0,0), fP4(0,0,0), fMiddle(0,0,0),
200  fNormal123(0,0,0), fNormal142(0,0,0), fNormal134(0,0,0),
201  fNormal234(0,0,0), warningFlag(0),
202  fCdotN123(0.), fCdotN142(0.), fCdotN134(0.), fCdotN234(0.),
203  fXMin(0.), fXMax(0.), fYMin(0.), fYMax(0.), fZMin(0.), fZMax(0.),
204  fDx(0.), fDy(0.), fDz(0.), fTol(0.), fMaxSize(0.)
205 {
206 }
207 
209 //
210 // Destructor
211 
213 {
214  delete fpPolyhedron; fpPolyhedron = 0;
215 }
216 
218 //
219 // Copy constructor
220 
221 G4Tet::G4Tet(const G4Tet& rhs)
222  : G4VSolid(rhs),
223  fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
224  fRebuildPolyhedron(false), fpPolyhedron(0), fAnchor(rhs.fAnchor),
225  fP2(rhs.fP2), fP3(rhs.fP3), fP4(rhs.fP4), fMiddle(rhs.fMiddle),
226  fNormal123(rhs.fNormal123), fNormal142(rhs.fNormal142),
227  fNormal134(rhs.fNormal134), fNormal234(rhs.fNormal234),
228  warningFlag(rhs.warningFlag), fCdotN123(rhs.fCdotN123),
229  fCdotN142(rhs.fCdotN142), fCdotN134(rhs.fCdotN134),
230  fCdotN234(rhs.fCdotN234), fXMin(rhs.fXMin), fXMax(rhs.fXMax),
231  fYMin(rhs.fYMin), fYMax(rhs.fYMax), fZMin(rhs.fZMin), fZMax(rhs.fZMax),
232  fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz), fTol(rhs.fTol),
233  fMaxSize(rhs.fMaxSize)
234 {
235 }
236 
237 
239 //
240 // Assignment operator
241 
243 {
244  // Check assignment to self
245  //
246  if (this == &rhs) { return *this; }
247 
248  // Copy base class data
249  //
250  G4VSolid::operator=(rhs);
251 
252  // Copy data
253  //
255  fAnchor = rhs.fAnchor;
256  fP2 = rhs.fP2; fP3 = rhs.fP3; fP4 = rhs.fP4; fMiddle = rhs.fMiddle;
260  fCdotN142 = rhs.fCdotN142; fCdotN134 = rhs.fCdotN134;
261  fCdotN234 = rhs.fCdotN234; fXMin = rhs.fXMin; fXMax = rhs.fXMax;
262  fYMin = rhs.fYMin; fYMax = rhs.fYMax; fZMin = rhs.fZMin; fZMax = rhs.fZMax;
263  fDx = rhs.fDx; fDy = rhs.fDy; fDz = rhs.fDz; fTol = rhs.fTol;
264  fMaxSize = rhs.fMaxSize;
265  fRebuildPolyhedron = false;
266  delete fpPolyhedron; fpPolyhedron = 0;
267 
268  return *this;
269 }
270 
272 //
273 // CheckDegeneracy
274 
276  G4ThreeVector p2,
277  G4ThreeVector p3,
278  G4ThreeVector p4 )
279 {
280  G4bool result;
281  G4Tet *object=new G4Tet("temp",anchor,p2,p3,p4,&result);
282  delete object;
283  return result;
284 }
285 
287 //
288 // Dispatch to parameterisation for replication mechanism dimension
289 // computation & modification.
290 
292  const G4int ,
293  const G4VPhysicalVolume* )
294 {
295 }
296 
298 //
299 // Calculate extent under transform and specified limit
300 
302  const G4VoxelLimits& pVoxelLimit,
303  const G4AffineTransform& pTransform,
304  G4double& pMin, G4double& pMax) const
305 {
306  G4double xMin,xMax;
307  G4double yMin,yMax;
308  G4double zMin,zMax;
309 
310  if (pTransform.IsRotated())
311  {
312  G4ThreeVector pp0=pTransform.TransformPoint(fAnchor);
313  G4ThreeVector pp1=pTransform.TransformPoint(fP2);
314  G4ThreeVector pp2=pTransform.TransformPoint(fP3);
315  G4ThreeVector pp3=pTransform.TransformPoint(fP4);
316 
317  xMin = std::min(std::min(std::min(pp0.x(), pp1.x()),pp2.x()),pp3.x());
318  xMax = std::max(std::max(std::max(pp0.x(), pp1.x()),pp2.x()),pp3.x());
319  yMin = std::min(std::min(std::min(pp0.y(), pp1.y()),pp2.y()),pp3.y());
320  yMax = std::max(std::max(std::max(pp0.y(), pp1.y()),pp2.y()),pp3.y());
321  zMin = std::min(std::min(std::min(pp0.z(), pp1.z()),pp2.z()),pp3.z());
322  zMax = std::max(std::max(std::max(pp0.z(), pp1.z()),pp2.z()),pp3.z());
323 
324  }
325  else
326  {
327  G4double xoffset = pTransform.NetTranslation().x() ;
328  xMin = xoffset + fXMin;
329  xMax = xoffset + fXMax;
330  G4double yoffset = pTransform.NetTranslation().y() ;
331  yMin = yoffset + fYMin;
332  yMax = yoffset + fYMax;
333  G4double zoffset = pTransform.NetTranslation().z() ;
334  zMin = zoffset + fZMin;
335  zMax = zoffset + fZMax;
336  }
337 
338  if (pVoxelLimit.IsXLimited())
339  {
340  if ( (xMin > pVoxelLimit.GetMaxXExtent()+fTol) ||
341  (xMax < pVoxelLimit.GetMinXExtent()-fTol) ) { return false; }
342  else
343  {
344  xMin = std::max(xMin, pVoxelLimit.GetMinXExtent());
345  xMax = std::min(xMax, pVoxelLimit.GetMaxXExtent());
346  }
347  }
348 
349  if (pVoxelLimit.IsYLimited())
350  {
351  if ( (yMin > pVoxelLimit.GetMaxYExtent()+fTol) ||
352  (yMax < pVoxelLimit.GetMinYExtent()-fTol) ) { return false; }
353  else
354  {
355  yMin = std::max(yMin, pVoxelLimit.GetMinYExtent());
356  yMax = std::min(yMax, pVoxelLimit.GetMaxYExtent());
357  }
358  }
359 
360  if (pVoxelLimit.IsZLimited())
361  {
362  if ( (zMin > pVoxelLimit.GetMaxZExtent()+fTol) ||
363  (zMax < pVoxelLimit.GetMinZExtent()-fTol) ) { return false; }
364  else
365  {
366  zMin = std::max(zMin, pVoxelLimit.GetMinZExtent());
367  zMax = std::min(zMax, pVoxelLimit.GetMaxZExtent());
368  }
369  }
370 
371  switch (pAxis)
372  {
373  case kXAxis:
374  pMin=xMin;
375  pMax=xMax;
376  break;
377  case kYAxis:
378  pMin=yMin;
379  pMax=yMax;
380  break;
381  case kZAxis:
382  pMin=zMin;
383  pMax=zMax;
384  break;
385  default:
386  break;
387  }
388 
389  return true;
390 }
391 
393 //
394 // Return whether point inside/outside/on surface, using tolerance
395 
397 {
398  G4double r123, r134, r142, r234;
399 
400  // this is written to allow if-statement truncation so the outside test
401  // (where most of the world is) can fail very quickly and efficiently
402 
403  if ( (r123=p.dot(fNormal123)-fCdotN123) > fTol ||
404  (r134=p.dot(fNormal134)-fCdotN134) > fTol ||
405  (r142=p.dot(fNormal142)-fCdotN142) > fTol ||
406  (r234=p.dot(fNormal234)-fCdotN234) > fTol )
407  {
408  return kOutside; // at least one is out!
409  }
410  else if( (r123 < -fTol)&&(r134 < -fTol)&&(r142 < -fTol)&&(r234 < -fTol) )
411  {
412  return kInside; // all are definitively inside
413  }
414  else
415  {
416  return kSurface; // too close to tell
417  }
418 }
419 
421 //
422 // Calculate side nearest to p, and return normal
423 // If two sides are equidistant, normal of first side (x/y/z)
424 // encountered returned.
425 // This assumes that we are looking from the inside!
426 
428 {
429  G4double r123=std::fabs(p.dot(fNormal123)-fCdotN123);
430  G4double r134=std::fabs(p.dot(fNormal134)-fCdotN134);
431  G4double r142=std::fabs(p.dot(fNormal142)-fCdotN142);
432  G4double r234=std::fabs(p.dot(fNormal234)-fCdotN234);
433 
434  const G4double delta = 0.5*kCarTolerance;
435  G4ThreeVector sumnorm(0., 0., 0.);
436  G4int noSurfaces=0;
437 
438  if (r123 <= delta)
439  {
440  noSurfaces ++;
441  sumnorm= fNormal123;
442  }
443 
444  if (r134 <= delta)
445  {
446  noSurfaces ++;
447  sumnorm += fNormal134;
448  }
449 
450  if (r142 <= delta)
451  {
452  noSurfaces ++;
453  sumnorm += fNormal142;
454  }
455  if (r234 <= delta)
456  {
457  noSurfaces ++;
458  sumnorm += fNormal234;
459  }
460 
461  if( noSurfaces > 0 )
462  {
463  if( noSurfaces == 1 )
464  {
465  return sumnorm;
466  }
467  else
468  {
469  return sumnorm.unit();
470  }
471  }
472  else // Approximative Surface Normal
473  {
474 
475  if( (r123<=r134) && (r123<=r142) && (r123<=r234) ) { return fNormal123; }
476  else if ( (r134<=r142) && (r134<=r234) ) { return fNormal134; }
477  else if (r142 <= r234) { return fNormal142; }
478  return fNormal234;
479  }
480 }
482 //
483 // Calculate distance to box from an outside point
484 // - return kInfinity if no intersection.
485 // All this is very unrolled, for speed.
486 
488  const G4ThreeVector& v) const
489 {
490  G4ThreeVector vu(v.unit()), hp;
491  G4double vdotn, t, tmin=kInfinity;
492 
493  G4double extraDistance=10.0*fTol; // a little ways into the solid
494 
495  vdotn=-vu.dot(fNormal123);
496  if(vdotn > 1e-12)
497  { // this is a candidate face, since it is pointing at us
498  t=(p.dot(fNormal123)-fCdotN123)/vdotn; // # distance to intersection
499  if( (t>=-fTol) && (t<tmin) )
500  { // if not true, we're going away from this face or it's not close
501  hp=p+vu*(t+extraDistance); // a little beyond point of intersection
502  if ( ( hp.dot(fNormal134)-fCdotN134 < 0.0 ) &&
503  ( hp.dot(fNormal142)-fCdotN142 < 0.0 ) &&
504  ( hp.dot(fNormal234)-fCdotN234 < 0.0 ) )
505  {
506  tmin=t;
507  }
508  }
509  }
510 
511  vdotn=-vu.dot(fNormal134);
512  if(vdotn > 1e-12)
513  { // # this is a candidate face, since it is pointing at us
514  t=(p.dot(fNormal134)-fCdotN134)/vdotn; // # distance to intersection
515  if( (t>=-fTol) && (t<tmin) )
516  { // if not true, we're going away from this face
517  hp=p+vu*(t+extraDistance); // a little beyond point of intersection
518  if ( ( hp.dot(fNormal123)-fCdotN123 < 0.0 ) &&
519  ( hp.dot(fNormal142)-fCdotN142 < 0.0 ) &&
520  ( hp.dot(fNormal234)-fCdotN234 < 0.0 ) )
521  {
522  tmin=t;
523  }
524  }
525  }
526 
527  vdotn=-vu.dot(fNormal142);
528  if(vdotn > 1e-12)
529  { // # this is a candidate face, since it is pointing at us
530  t=(p.dot(fNormal142)-fCdotN142)/vdotn; // # distance to intersection
531  if( (t>=-fTol) && (t<tmin) )
532  { // if not true, we're going away from this face
533  hp=p+vu*(t+extraDistance); // a little beyond point of intersection
534  if ( ( hp.dot(fNormal123)-fCdotN123 < 0.0 ) &&
535  ( hp.dot(fNormal134)-fCdotN134 < 0.0 ) &&
536  ( hp.dot(fNormal234)-fCdotN234 < 0.0 ) )
537  {
538  tmin=t;
539  }
540  }
541  }
542 
543  vdotn=-vu.dot(fNormal234);
544  if(vdotn > 1e-12)
545  { // # this is a candidate face, since it is pointing at us
546  t=(p.dot(fNormal234)-fCdotN234)/vdotn; // # distance to intersection
547  if( (t>=-fTol) && (t<tmin) )
548  { // if not true, we're going away from this face
549  hp=p+vu*(t+extraDistance); // a little beyond point of intersection
550  if ( ( hp.dot(fNormal123)-fCdotN123 < 0.0 ) &&
551  ( hp.dot(fNormal134)-fCdotN134 < 0.0 ) &&
552  ( hp.dot(fNormal142)-fCdotN142 < 0.0 ) )
553  {
554  tmin=t;
555  }
556  }
557  }
558 
559  return std::max(0.0,tmin);
560 }
561 
563 //
564 // Approximate distance to tet.
565 // returns distance to sphere centered on bounding box
566 // - If inside return 0
567 
569 {
570  G4double dd=(p-fMiddle).mag() - fMaxSize - fTol;
571  return std::max(0.0, dd);
572 }
573 
575 //
576 // Calcluate distance to surface of box from inside
577 // by calculating distances to box's x/y/z planes.
578 // Smallest distance is exact distance to exiting.
579 
581  const G4bool calcNorm,
582  G4bool *validNorm, G4ThreeVector *n) const
583 {
584  G4ThreeVector vu(v.unit());
585  G4double t1=kInfinity,t2=kInfinity,t3=kInfinity,t4=kInfinity, vdotn, tt;
586 
587  vdotn=vu.dot(fNormal123);
588  if(vdotn > 1e-12) // #we're heading towards this face, so it is a candidate
589  {
590  t1=(fCdotN123-p.dot(fNormal123))/vdotn; // # distance to intersection
591  }
592 
593  vdotn=vu.dot(fNormal134);
594  if(vdotn > 1e-12) // #we're heading towards this face, so it is a candidate
595  {
596  t2=(fCdotN134-p.dot(fNormal134))/vdotn; // # distance to intersection
597  }
598 
599  vdotn=vu.dot(fNormal142);
600  if(vdotn > 1e-12) // #we're heading towards this face, so it is a candidate
601  {
602  t3=(fCdotN142-p.dot(fNormal142))/vdotn; // # distance to intersection
603  }
604 
605  vdotn=vu.dot(fNormal234);
606  if(vdotn > 1e-12) // #we're heading towards this face, so it is a candidate
607  {
608  t4=(fCdotN234-p.dot(fNormal234))/vdotn; // # distance to intersection
609  }
610 
611  tt=std::min(std::min(std::min(t1,t2),t3),t4);
612 
613  if (warningFlag && (tt == kInfinity || tt < -fTol))
614  {
615  DumpInfo();
616  std::ostringstream message;
617  message << "No good intersection found or already outside!?" << G4endl
618  << "p = " << p / mm << "mm" << G4endl
619  << "v = " << v << G4endl
620  << "t1, t2, t3, t4 (mm) "
621  << t1/mm << ", " << t2/mm << ", " << t3/mm << ", " << t4/mm;
622  G4Exception("G4Tet::DistanceToOut(p,v,...)", "GeomSolids1002",
623  JustWarning, message);
624  if(validNorm)
625  {
626  *validNorm=false; // flag normal as meaningless
627  }
628  }
629  else if(calcNorm && n)
630  {
632  if(tt==t1) { normal=fNormal123; }
633  else if (tt==t2) { normal=fNormal134; }
634  else if (tt==t3) { normal=fNormal142; }
635  else if (tt==t4) { normal=fNormal234; }
636  *n=normal;
637  if(validNorm) { *validNorm=true; }
638  }
639 
640  return std::max(tt,0.0); // avoid tt<0.0 by a tiny bit
641  // if we are right on a face
642 }
643 
645 //
646 // Calculate exact shortest distance to any boundary from inside
647 // - If outside return 0
648 
650 {
651  G4double t1,t2,t3,t4;
652  t1=fCdotN123-p.dot(fNormal123); // distance to plane, positive if inside
653  t2=fCdotN134-p.dot(fNormal134); // distance to plane
654  t3=fCdotN142-p.dot(fNormal142); // distance to plane
655  t4=fCdotN234-p.dot(fNormal234); // distance to plane
656 
657  // if any one of these is negative, we are outside,
658  // so return zero in that case
659 
660  G4double tmin=std::min(std::min(std::min(t1,t2),t3),t4);
661  return (tmin < fTol)? 0:tmin;
662 }
663 
665 //
666 // Create a List containing the transformed vertices
667 // Note: Caller has deletion responsibility
668 
671 {
672  G4ThreeVectorList* vertices = new G4ThreeVectorList();
673 
674  if (vertices)
675  {
676  vertices->reserve(4);
677  G4ThreeVector vertex0(fAnchor);
678  G4ThreeVector vertex1(fP2);
679  G4ThreeVector vertex2(fP3);
680  G4ThreeVector vertex3(fP4);
681 
682  vertices->push_back(pTransform.TransformPoint(vertex0));
683  vertices->push_back(pTransform.TransformPoint(vertex1));
684  vertices->push_back(pTransform.TransformPoint(vertex2));
685  vertices->push_back(pTransform.TransformPoint(vertex3));
686  }
687  else
688  {
689  DumpInfo();
690  G4Exception("G4Tet::CreateRotatedVertices()",
691  "GeomSolids0003", FatalException,
692  "Error in allocation of vertices. Out of memory !");
693  }
694  return vertices;
695 }
696 
698 //
699 // GetEntityType
700 
702 {
703  return G4String("G4Tet");
704 }
705 
707 //
708 // Make a clone of the object
709 
711 {
712  return new G4Tet(*this);
713 }
714 
716 //
717 // Stream object contents to an output stream
718 
719 std::ostream& G4Tet::StreamInfo(std::ostream& os) const
720 {
721  G4int oldprc = os.precision(16);
722  os << "-----------------------------------------------------------\n"
723  << " *** Dump for solid - " << GetName() << " ***\n"
724  << " ===================================================\n"
725  << " Solid type: G4Tet\n"
726  << " Parameters: \n"
727  << " anchor: " << fAnchor/mm << " mm \n"
728  << " p2: " << fP2/mm << " mm \n"
729  << " p3: " << fP3/mm << " mm \n"
730  << " p4: " << fP4/mm << " mm \n"
731  << " normal123: " << fNormal123 << " \n"
732  << " normal134: " << fNormal134 << " \n"
733  << " normal142: " << fNormal142 << " \n"
734  << " normal234: " << fNormal234 << " \n"
735  << "-----------------------------------------------------------\n";
736  os.precision(oldprc);
737 
738  return os;
739 }
740 
741 
743 //
744 // GetPointOnFace
745 //
746 // Auxiliary method for get point on surface
747 
749  G4ThreeVector p3, G4double& area) const
750 {
751  G4double lambda1,lambda2;
752  G4ThreeVector v, w;
753 
754  v = p3 - p1;
755  w = p1 - p2;
756 
757  lambda1 = RandFlat::shoot(0.,1.);
758  lambda2 = RandFlat::shoot(0.,lambda1);
759 
760  area = 0.5*(v.cross(w)).mag();
761 
762  return (p2 + lambda1*w + lambda2*v);
763 }
764 
766 //
767 // GetPointOnSurface
768 
770 {
771  G4double chose,aOne,aTwo,aThree,aFour;
772  G4ThreeVector p1, p2, p3, p4;
773 
774  p1 = GetPointOnFace(fAnchor,fP2,fP3,aOne);
775  p2 = GetPointOnFace(fAnchor,fP4,fP3,aTwo);
776  p3 = GetPointOnFace(fAnchor,fP4,fP2,aThree);
777  p4 = GetPointOnFace(fP4,fP3,fP2,aFour);
778 
779  chose = RandFlat::shoot(0.,aOne+aTwo+aThree+aFour);
780  if( (chose>=0.) && (chose <aOne) ) {return p1;}
781  else if( (chose>=aOne) && (chose < aOne+aTwo) ) {return p2;}
782  else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThree) ) {return p3;}
783  return p4;
784 }
785 
787 //
788 // GetVertices
789 
790 std::vector<G4ThreeVector> G4Tet::GetVertices() const
791 {
792  std::vector<G4ThreeVector> vertices(4);
793  vertices[0] = fAnchor;
794  vertices[1] = fP2;
795  vertices[2] = fP3;
796  vertices[3] = fP4;
797 
798  return vertices;
799 }
800 
802 //
803 // GetCubicVolume
804 
806 {
807  return fCubicVolume;
808 }
809 
811 //
812 // GetSurfaceArea
813 
815 {
816  return fSurfaceArea;
817 }
818 
819 // Methods for visualisation
820 
822 //
823 // DescribeYourselfTo
824 
826 {
827  scene.AddSolid (*this);
828 }
829 
831 //
832 // GetExtent
833 
835 {
836  return G4VisExtent (fXMin, fXMax, fYMin, fYMax, fZMin, fZMax);
837 }
838 
840 //
841 // CreatePolyhedron
842 
844 {
845  G4Polyhedron *ph=new G4Polyhedron;
846  G4double xyz[4][3];
847  const G4int faces[4][4]={{1,3,2,0},{1,4,3,0},{1,2,4,0},{2,3,4,0}};
848  xyz[0][0]=fAnchor.x(); xyz[0][1]=fAnchor.y(); xyz[0][2]=fAnchor.z();
849  xyz[1][0]=fP2.x(); xyz[1][1]=fP2.y(); xyz[1][2]=fP2.z();
850  xyz[2][0]=fP3.x(); xyz[2][1]=fP3.y(); xyz[2][2]=fP3.z();
851  xyz[3][0]=fP4.x(); xyz[3][1]=fP4.y(); xyz[3][2]=fP4.z();
852 
853  ph->createPolyhedron(4,4,xyz,faces);
854 
855  return ph;
856 }
857 
859 //
860 // GetPolyhedron
861 
863 {
864  if (!fpPolyhedron ||
867  fpPolyhedron->GetNumberOfRotationSteps())
868  {
869  G4AutoLock l(&polyhedronMutex);
870  delete fpPolyhedron;
872  fRebuildPolyhedron = false;
873  l.unlock();
874  }
875  return fpPolyhedron;
876 }
877 
878 #endif
G4ThreeVector fP4
Definition: G4Tet.hh:165
G4double fTol
Definition: G4Tet.hh:172
G4double fYMin
Definition: G4Tet.hh:171
G4String GetName() const
G4double fXMax
Definition: G4Tet.hh:171
ThreeVector shoot(const G4int Ap, const G4int Af)
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition: G4Tet.cc:580
G4ThreeVector fNormal142
Definition: G4Tet.hh:166
G4bool fRebuildPolyhedron
Definition: G4Tet.hh:156
G4double fDz
Definition: G4Tet.hh:172
G4ThreeVector fP2
Definition: G4Tet.hh:165
static const G4double kInfinity
Definition: geomdefs.hh:42
G4double GetMinYExtent() const
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector GetPointOnSurface() const
Definition: G4Tet.cc:769
G4double fYMax
Definition: G4Tet.hh:171
G4bool IsYLimited() const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Tet.cc:487
G4GeometryType GetEntityType() const
Definition: G4Tet.cc:701
G4bool IsRotated() const
const G4double w[NPOINTSGL]
G4ThreeVector fNormal123
Definition: G4Tet.hh:166
G4double fCdotN142
Definition: G4Tet.hh:170
G4ThreeVector fP3
Definition: G4Tet.hh:165
G4ThreeVector NetTranslation() const
G4bool IsXLimited() const
G4double a
Definition: TRTMaterials.hh:39
G4double fCdotN134
Definition: G4Tet.hh:170
virtual void AddSolid(const G4Box &)=0
int G4int
Definition: G4Types.hh:78
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
Definition: G4Tet.cc:301
void DumpInfo() const
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4double fMaxSize
Definition: G4Tet.hh:172
std::vector< G4ThreeVector > GetVertices() const
Definition: G4Tet.cc:790
G4double GetMaxXExtent() const
G4ThreeVector fAnchor
Definition: G4Tet.hh:165
G4double GetMinZExtent() const
G4double fZMin
Definition: G4Tet.hh:171
G4double fXMin
Definition: G4Tet.hh:171
Definition: G4Tet.hh:65
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Tet.cc:825
bool G4bool
Definition: G4Types.hh:79
G4ThreeVector GetPointOnFace(G4ThreeVector p1, G4ThreeVector p2, G4ThreeVector p3, G4double &area) const
Definition: G4Tet.cc:748
static G4bool CheckDegeneracy(G4ThreeVector anchor, G4ThreeVector p2, G4ThreeVector p3, G4ThreeVector p4)
Definition: G4Tet.cc:275
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:79
G4ThreeVector fMiddle
Definition: G4Tet.hh:165
const G4int n
G4double fSurfaceArea
Definition: G4Tet.hh:154
G4Polyhedron * CreatePolyhedron() const
Definition: G4Tet.cc:843
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double fCdotN234
Definition: G4Tet.hh:170
G4double GetMinXExtent() const
G4int G4Mutex
Definition: G4Threading.hh:173
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4double GetMaxZExtent() const
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Tet.cc:719
G4double fDy
Definition: G4Tet.hh:172
T max(const T t1, const T t2)
brief Return the largest of the two arguments
EInside
Definition: geomdefs.hh:58
G4ThreeVector fNormal134
Definition: G4Tet.hh:166
EAxis
Definition: geomdefs.hh:54
G4double GetCubicVolume()
Definition: G4Tet.cc:805
G4double fCdotN123
Definition: G4Tet.hh:170
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double fCubicVolume
Definition: G4Tet.hh:154
G4double fDx
Definition: G4Tet.hh:172
G4double GetSurfaceArea()
Definition: G4Tet.cc:814
EInside Inside(const G4ThreeVector &p) const
Definition: G4Tet.cc:396
G4VisExtent GetExtent() const
Definition: G4Tet.cc:834
#define G4endl
Definition: G4ios.hh:61
G4double GetMaxYExtent() const
static const char CVSVers[]
Definition: G4Tet.hh:161
G4VSolid & operator=(const G4VSolid &rhs)
Definition: G4VSolid.cc:110
G4double kCarTolerance
Definition: G4VSolid.hh:304
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
G4VSolid * Clone() const
Definition: G4Tet.cc:710
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform) const
Definition: G4Tet.cc:670
G4ThreeVector fNormal234
Definition: G4Tet.hh:166
G4Tet & operator=(const G4Tet &rhs)
Definition: G4Tet.cc:242
double G4double
Definition: G4Types.hh:76
virtual ~G4Tet()
Definition: G4Tet.cc:212
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Tet.cc:291
G4bool IsZLimited() const
static const double mm
Definition: G4SIunits.hh:114
G4Polyhedron * fpPolyhedron
Definition: G4Tet.hh:157
G4Polyhedron * GetPolyhedron() const
Definition: G4Tet.cc:862
G4bool warningFlag
Definition: G4Tet.hh:168
G4Tet(const G4String &pName, G4ThreeVector anchor, G4ThreeVector p2, G4ThreeVector p3, G4ThreeVector p4, G4bool *degeneracyFlag=0)
Definition: G4Tet.cc:96
G4double fZMax
Definition: G4Tet.hh:171
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Tet.cc:427