source: trunk/source/geometry/solids/specific/src/G4VCSGfaceted.cc@ 1202

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

file release beta
File size: 13.7 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 result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// the GEANT4 collaboration.
27//
28// By copying, distributing or modifying the Program (or any work
29// based on the Program) you indicate your acceptance of this statement,
30// and all its terms.
31//
32// $Id: G4VCSGfaceted.cc,v 1.25 2008/05/22 10:22:52 gcosmo Exp $
33// GEANT4 tag $Name: geant4-09-02-ref-02 $
34//
35//
36// --------------------------------------------------------------------
37// GEANT 4 class source file
38//
39//
40// G4VCSGfaceted.cc
41//
42// Implementation of the virtual class of a CSG type shape that is built
43// entirely out of G4VCSGface faces.
44//
45// --------------------------------------------------------------------
46
47#include "G4VCSGfaceted.hh"
48#include "G4VCSGface.hh"
49#include "G4SolidExtentList.hh"
50
51#include "G4VoxelLimits.hh"
52#include "G4AffineTransform.hh"
53
54#include "Randomize.hh"
55
56#include "G4Polyhedron.hh"
57#include "G4VGraphicsScene.hh"
58#include "G4NURBS.hh"
59#include "G4NURBSbox.hh"
60#include "G4VisExtent.hh"
61
62//
63// Constructor
64//
65G4VCSGfaceted::G4VCSGfaceted( const G4String& name )
66 : G4VSolid(name),
67 numFace(0), faces(0), fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0),
68 fStatistics(1000000), fCubVolEpsilon(0.001), fAreaAccuracy(-1.)
69{
70}
71
72
73//
74// Fake default constructor - sets only member data and allocates memory
75// for usage restricted to object persistency.
76//
77G4VCSGfaceted::G4VCSGfaceted( __void__& a )
78 : G4VSolid(a),
79 numFace(0), faces(0), fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0),
80 fStatistics(1000000), fCubVolEpsilon(0.001), fAreaAccuracy(-1.)
81{
82}
83
84//
85// Destructor
86//
87G4VCSGfaceted::~G4VCSGfaceted()
88{
89 DeleteStuff();
90 delete fpPolyhedron;
91}
92
93
94//
95// Copy constructor
96//
97G4VCSGfaceted::G4VCSGfaceted( const G4VCSGfaceted &source )
98 : G4VSolid( source )
99{
100 CopyStuff( source );
101}
102
103
104//
105// Assignment operator
106//
107const G4VCSGfaceted &G4VCSGfaceted::operator=( const G4VCSGfaceted &source )
108{
109 if (&source == this) { return *this; }
110
111 DeleteStuff();
112 CopyStuff( source );
113
114 return *this;
115}
116
117
118//
119// CopyStuff (protected)
120//
121// Copy the contents of source
122//
123void G4VCSGfaceted::CopyStuff( const G4VCSGfaceted &source )
124{
125 numFace = source.numFace;
126 if (numFace == 0) { return; } // odd, but permissable?
127
128 faces = new G4VCSGface*[numFace];
129
130 G4VCSGface **face = faces,
131 **sourceFace = source.faces;
132 do
133 {
134 *face = (*sourceFace)->Clone();
135 } while( ++sourceFace, ++face < faces+numFace );
136 fCubicVolume = source.fCubicVolume;
137 fpPolyhedron = source.fpPolyhedron;
138}
139
140
141//
142// DeleteStuff (protected)
143//
144// Delete all allocated objects
145//
146void G4VCSGfaceted::DeleteStuff()
147{
148 if (numFace)
149 {
150 G4VCSGface **face = faces;
151 do
152 {
153 delete *face;
154 } while( ++face < faces + numFace );
155
156 delete [] faces;
157 }
158}
159
160
161//
162// CalculateExtent
163//
164G4bool G4VCSGfaceted::CalculateExtent( const EAxis axis,
165 const G4VoxelLimits &voxelLimit,
166 const G4AffineTransform &transform,
167 G4double &min,
168 G4double &max ) const
169{
170 G4SolidExtentList extentList( axis, voxelLimit );
171
172 //
173 // Loop over all faces, checking min/max extent as we go.
174 //
175 G4VCSGface **face = faces;
176 do
177 {
178 (*face)->CalculateExtent( axis, voxelLimit, transform, extentList );
179 } while( ++face < faces + numFace );
180
181 //
182 // Return min/max value
183 //
184 return extentList.GetExtent( min, max );
185}
186
187
188//
189// Inside
190//
191// It could be a good idea to override this virtual
192// member to add first a simple test (such as spherical
193// test or whatnot) and to call this version only if
194// the simplier test fails.
195//
196EInside G4VCSGfaceted::Inside( const G4ThreeVector &p ) const
197{
198 EInside answer=kOutside;
199 G4VCSGface **face = faces;
200 G4double best = kInfinity;
201 do
202 {
203 G4double distance;
204 EInside result = (*face)->Inside( p, kCarTolerance/2, &distance );
205 if (result == kSurface) { return kSurface; }
206 if (distance < best)
207 {
208 best = distance;
209 answer = result;
210 }
211 } while( ++face < faces + numFace );
212
213 return answer;
214}
215
216
217//
218// SurfaceNormal
219//
220G4ThreeVector G4VCSGfaceted::SurfaceNormal( const G4ThreeVector& p ) const
221{
222 G4ThreeVector answer;
223 G4VCSGface **face = faces;
224 G4double best = kInfinity;
225 do
226 {
227 G4double distance;
228 G4ThreeVector normal = (*face)->Normal( p, &distance );
229 if (distance < best)
230 {
231 best = distance;
232 answer = normal;
233 }
234 } while( ++face < faces + numFace );
235
236 return answer;
237}
238
239
240//
241// DistanceToIn(p,v)
242//
243G4double G4VCSGfaceted::DistanceToIn( const G4ThreeVector &p,
244 const G4ThreeVector &v ) const
245{
246 G4double distance = kInfinity;
247 G4double distFromSurface = kInfinity;
248 G4VCSGface *bestFace=0;
249 G4VCSGface **face = faces;
250 do
251 {
252 G4double faceDistance,
253 faceDistFromSurface;
254 G4ThreeVector faceNormal;
255 G4bool faceAllBehind;
256 if ((*face)->Intersect( p, v, false, kCarTolerance/2,
257 faceDistance, faceDistFromSurface,
258 faceNormal, faceAllBehind ) )
259 {
260 //
261 // Intersecting face
262 //
263 if (faceDistance < distance)
264 {
265 distance = faceDistance;
266 distFromSurface = faceDistFromSurface;
267 bestFace = *face;
268 if (distFromSurface <= 0) { return 0; }
269 }
270 }
271 } while( ++face < faces + numFace );
272
273 if (distance < kInfinity && distFromSurface<kCarTolerance/2)
274 {
275 if (bestFace->Distance(p,false) < kCarTolerance/2) { distance = 0; }
276 }
277
278 return distance;
279}
280
281
282//
283// DistanceToIn(p)
284//
285G4double G4VCSGfaceted::DistanceToIn( const G4ThreeVector &p ) const
286{
287 return DistanceTo( p, false );
288}
289
290
291//
292// DistanceToOut(p,v)
293//
294G4double G4VCSGfaceted::DistanceToOut( const G4ThreeVector &p,
295 const G4ThreeVector &v,
296 const G4bool calcNorm,
297 G4bool *validNorm,
298 G4ThreeVector *n ) const
299{
300 G4bool allBehind = true;
301 G4double distance = kInfinity;
302 G4double distFromSurface = kInfinity;
303 G4ThreeVector normal;
304 G4VCSGface *bestFace=0;
305
306 G4VCSGface **face = faces;
307 do
308 {
309 G4double faceDistance,
310 faceDistFromSurface;
311 G4ThreeVector faceNormal;
312 G4bool faceAllBehind;
313 if ((*face)->Intersect( p, v, true, kCarTolerance/2,
314 faceDistance, faceDistFromSurface,
315 faceNormal, faceAllBehind ) )
316 {
317 //
318 // Intersecting face
319 //
320 if ( (distance < kInfinity) || (!faceAllBehind) ) { allBehind = false; }
321 if (faceDistance < distance)
322 {
323 distance = faceDistance;
324 distFromSurface = faceDistFromSurface;
325 normal = faceNormal;
326 bestFace = *face;
327 if (distFromSurface <= 0) { break; }
328 }
329 }
330 } while( ++face < faces + numFace );
331
332 if (distance < kInfinity)
333 {
334 if (distFromSurface <= 0)
335 {
336 distance = 0;
337 }
338 else if (distFromSurface<kCarTolerance/2)
339 {
340 if (bestFace->Distance(p,true) < kCarTolerance/2) { distance = 0; }
341 }
342
343 if (calcNorm)
344 {
345 *validNorm = allBehind;
346 *n = normal;
347 }
348 }
349 else
350 {
351 if (Inside(p) == kSurface) { distance = 0; }
352 if (calcNorm) { *validNorm = false; }
353 }
354
355 return distance;
356}
357
358
359//
360// DistanceToOut(p)
361//
362G4double G4VCSGfaceted::DistanceToOut( const G4ThreeVector &p ) const
363{
364 return DistanceTo( p, true );
365}
366
367
368//
369// DistanceTo
370//
371// Protected routine called by DistanceToIn and DistanceToOut
372//
373G4double G4VCSGfaceted::DistanceTo( const G4ThreeVector &p,
374 const G4bool outgoing ) const
375{
376 G4VCSGface **face = faces;
377 G4double best = kInfinity;
378 do
379 {
380 G4double distance = (*face)->Distance( p, outgoing );
381 if (distance < best) { best = distance; }
382 } while( ++face < faces + numFace );
383
384 return (best < 0.5*kCarTolerance) ? 0 : best;
385}
386
387
388//
389// DescribeYourselfTo
390//
391void G4VCSGfaceted::DescribeYourselfTo( G4VGraphicsScene& scene ) const
392{
393 scene.AddSolid( *this );
394}
395
396
397//
398// GetExtent
399//
400// Define the sides of the box into which our solid instance would fit.
401//
402G4VisExtent G4VCSGfaceted::GetExtent() const
403{
404 static const G4ThreeVector xMax(1,0,0), xMin(-1,0,0),
405 yMax(0,1,0), yMin(0,-1,0),
406 zMax(0,0,1), zMin(0,0,-1);
407 static const G4ThreeVector *axes[6] =
408 { &xMin, &xMax, &yMin, &yMax, &zMin, &zMax };
409
410 G4double answers[6] =
411 {-kInfinity, -kInfinity, -kInfinity, -kInfinity, -kInfinity, -kInfinity};
412
413 G4VCSGface **face = faces;
414 do
415 {
416 const G4ThreeVector **axis = axes+5 ;
417 G4double *answer = answers+5;
418 do
419 {
420 G4double testFace = (*face)->Extent( **axis );
421 if (testFace > *answer) { *answer = testFace; }
422 }
423 while( --axis, --answer >= answers );
424
425 } while( ++face < faces + numFace );
426
427 return G4VisExtent( -answers[0], answers[1],
428 -answers[2], answers[3],
429 -answers[4], answers[5] );
430}
431
432
433//
434// GetEntityType
435//
436G4GeometryType G4VCSGfaceted::GetEntityType() const
437{
438 return G4String("G4CSGfaceted");
439}
440
441
442//
443// Stream object contents to an output stream
444//
445std::ostream& G4VCSGfaceted::StreamInfo( std::ostream& os ) const
446{
447 os << "-----------------------------------------------------------\n"
448 << " *** Dump for solid - " << GetName() << " ***\n"
449 << " ===================================================\n"
450 << " Solid type: G4VCSGfaceted\n"
451 << " Parameters: \n"
452 << " number of faces: " << numFace << "\n"
453 << "-----------------------------------------------------------\n";
454
455 return os;
456}
457
458
459//
460// GetCubVolStatistics
461//
462G4int G4VCSGfaceted::GetCubVolStatistics() const
463{
464 return fStatistics;
465}
466
467
468//
469// GetCubVolEpsilon
470//
471G4double G4VCSGfaceted::GetCubVolEpsilon() const
472{
473 return fCubVolEpsilon;
474}
475
476
477//
478// SetCubVolStatistics
479//
480void G4VCSGfaceted::SetCubVolStatistics(G4int st)
481{
482 fCubicVolume=0.;
483 fStatistics=st;
484}
485
486
487//
488// SetCubVolEpsilon
489//
490void G4VCSGfaceted::SetCubVolEpsilon(G4double ep)
491{
492 fCubicVolume=0.;
493 fCubVolEpsilon=ep;
494}
495
496
497//
498// GetAreaStatistics
499//
500G4int G4VCSGfaceted::GetAreaStatistics() const
501{
502 return fStatistics;
503}
504
505
506//
507// GetAreaAccuracy
508//
509G4double G4VCSGfaceted::GetAreaAccuracy() const
510{
511 return fAreaAccuracy;
512}
513
514
515//
516// SetAreaStatistics
517//
518void G4VCSGfaceted::SetAreaStatistics(G4int st)
519{
520 fSurfaceArea=0.;
521 fStatistics=st;
522}
523
524
525//
526// SetAreaAccuracy
527//
528void G4VCSGfaceted::SetAreaAccuracy(G4double ep)
529{
530 fSurfaceArea=0.;
531 fAreaAccuracy=ep;
532}
533
534
535//
536// GetCubicVolume
537//
538G4double G4VCSGfaceted::GetCubicVolume()
539{
540 if(fCubicVolume != 0.) {;}
541 else { fCubicVolume = EstimateCubicVolume(fStatistics,fCubVolEpsilon); }
542 return fCubicVolume;
543}
544
545
546//
547// GetSurfaceArea
548//
549G4double G4VCSGfaceted::GetSurfaceArea()
550{
551 if(fSurfaceArea != 0.) {;}
552 else { fSurfaceArea = EstimateCubicVolume(fStatistics,fAreaAccuracy); }
553 return fSurfaceArea;
554}
555
556
557//
558// GetPolyhedron
559//
560G4Polyhedron* G4VCSGfaceted::GetPolyhedron () const
561{
562 if (!fpPolyhedron ||
563 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
564 fpPolyhedron->GetNumberOfRotationSteps())
565 {
566 delete fpPolyhedron;
567 fpPolyhedron = CreatePolyhedron();
568 }
569 return fpPolyhedron;
570}
571
572
573//
574// GetPointOnSurfaceGeneric proportional to Areas of faces
575// in case of GenericPolycone or GenericPolyhedra
576//
577G4ThreeVector G4VCSGfaceted::GetPointOnSurfaceGeneric( ) const
578{
579 // Preparing variables
580 //
581 G4ThreeVector answer=G4ThreeVector(0.,0.,0.);
582 G4VCSGface **face = faces;
583 G4double area = 0;
584 G4int i;
585 std::vector<G4double> areas;
586
587 // First step: calculate surface areas
588 //
589 do
590 {
591 G4double result = (*face)->SurfaceArea( );
592 areas.push_back(result);
593 area=area+result;
594 } while( ++face < faces + numFace );
595
596 // Second Step: choose randomly one surface
597 //
598 G4VCSGface **face1 = faces;
599 G4double chose = area*G4UniformRand();
600 G4double Achose1, Achose2;
601 Achose1=0; Achose2=0.;
602 i=0;
603
604 do
605 {
606 Achose2+=areas[i];
607 if(chose>=Achose1 && chose<Achose2)
608 {
609 G4ThreeVector point;
610 point= (*face1)->GetPointOnFace();
611 return point;
612 }
613 i++;
614 Achose1=Achose2;
615 } while( ++face1 < faces + numFace );
616
617 return answer;
618}
Note: See TracBrowser for help on using the repository browser.