source: trunk/source/geometry/solids/specific/src/G4Tet.cc@ 1127

Last change on this file since 1127 was 1058, checked in by garnier, 17 years ago

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