source: trunk/source/geometry/management/src/G4VSolid.cc@ 989

Last change on this file since 989 was 921, checked in by garnier, 17 years ago

en test de gl2ps. Problemes de libraries
File size: 20.1 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//
27// $Id: G4VSolid.cc,v 1.39 2008/09/23 13:07:41 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-02-cand-01 $
29//
30// class G4VSolid
31//
32// Implementation for solid base class
33//
34// History:
35//
36// 06.12.02 V.Grichine, restored original conditions in ClipPolygon()
37// 10.05.02 V.Grichine, ClipPolygon(): clip only other axis and limited voxels
38// 15.04.02 V.Grichine, bug fixed in ClipPolygon(): clip only one axis
39// 13.03.02 V.Grichine, cosmetics of voxel limit functions
40// 15.11.00 D.Williams, V.Grichine, fix in CalculateClippedPolygonExtent()
41// 10.07.95 P.Kent, Added == operator, solid Store entry
42// 30.06.95 P.Kent, Created.
43// --------------------------------------------------------------------
44
45#include "G4VSolid.hh"
46#include "G4SolidStore.hh"
47#include "globals.hh"
48#include "Randomize.hh"
49#include "G4GeometryTolerance.hh"
50
51#include "G4VoxelLimits.hh"
52#include "G4AffineTransform.hh"
53#include "G4VisExtent.hh"
54
55//////////////////////////////////////////////////////////////////////////
56//
57// Constructor
58// - Copies name
59// - Add ourselves to solid Store
60
61G4VSolid::G4VSolid(const G4String& name)
62 : fshapeName(name)
63{
64 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
65
66 // Register to store
67 //
68 G4SolidStore::GetInstance()->Register(this);
69}
70
71//////////////////////////////////////////////////////////////////////////
72//
73// Copy constructor
74//
75
76G4VSolid::G4VSolid(const G4VSolid& rhs)
77 : kCarTolerance(rhs.kCarTolerance), fshapeName(rhs.fshapeName)
78{
79 // Register to store
80 //
81 G4SolidStore::GetInstance()->Register(this);
82}
83
84//////////////////////////////////////////////////////////////////////////
85//
86// Fake default constructor - sets only member data and allocates memory
87// for usage restricted to object persistency.
88//
89G4VSolid::G4VSolid( __void__& )
90 : fshapeName("")
91{
92 // Register to store
93 //
94 G4SolidStore::GetInstance()->Register(this);
95}
96
97//////////////////////////////////////////////////////////////////////////
98//
99// Destructor (virtual)
100// - Remove ourselves from solid Store
101
102G4VSolid::~G4VSolid()
103{
104 G4SolidStore::GetInstance()->DeRegister(this);
105}
106
107//////////////////////////////////////////////////////////////////////////
108//
109// Assignment operator
110
111G4VSolid& G4VSolid::operator = (const G4VSolid& rhs)
112{
113 // Check assignment to self
114 //
115 if (this == &rhs) { return *this; }
116
117 // Copy data
118 //
119 kCarTolerance = rhs.kCarTolerance;
120 fshapeName = rhs.fshapeName;
121
122 return *this;
123}
124
125//////////////////////////////////////////////////////////////////////////
126//
127// Streaming operator dumping solid contents
128
129std::ostream& operator<< ( std::ostream& os, const G4VSolid& e )
130{
131 return e.StreamInfo(os);
132}
133
134//////////////////////////////////////////////////////////////////////////
135//
136// Throw exception if ComputeDimensions called for illegal derived class
137
138void G4VSolid::ComputeDimensions(G4VPVParameterisation*,
139 const G4int,
140 const G4VPhysicalVolume*)
141{
142 G4cerr << "ERROR - Illegal call to G4VSolid::ComputeDimensions()" << G4endl
143 << " Method not overloaded by derived class !" << G4endl;
144 G4Exception("G4VSolid::ComputeDimensions()", "NotApplicable",
145 FatalException, "Illegal call to case class.");
146}
147
148//////////////////////////////////////////////////////////////////////////
149//
150// Throw exception (warning) for solids not implementing the method
151
152G4ThreeVector G4VSolid::GetPointOnSurface() const
153{
154 G4cerr << "WARNING - G4VSolid::GetPointOnSurface()" << G4endl
155 << " Not implemented for solid: "
156 << this->GetEntityType() << " !" << G4endl;
157 G4Exception("G4VSolid::GetPointOnSurface()", "NotImplemented",
158 JustWarning, "Not implemented for this solid ! Returning origin.");
159 return G4ThreeVector(0,0,0);
160}
161
162///////////////////////////////////////////////////////////////////////////
163//
164// Calculate the maximum and minimum extents of the polygon described
165// by the vertices: pSectionIndex->pSectionIndex+1->
166// pSectionIndex+2->pSectionIndex+3->pSectionIndex
167// in the List pVertices
168//
169// If the minimum is <pMin pMin is set to the new minimum
170// If the maximum is >pMax pMax is set to the new maximum
171//
172// No modifications are made to pVertices
173//
174
175void G4VSolid::ClipCrossSection( G4ThreeVectorList* pVertices,
176 const G4int pSectionIndex,
177 const G4VoxelLimits& pVoxelLimit,
178 const EAxis pAxis,
179 G4double& pMin, G4double& pMax) const
180{
181
182 G4ThreeVectorList polygon;
183 polygon.reserve(4);
184 polygon.push_back((*pVertices)[pSectionIndex]);
185 polygon.push_back((*pVertices)[pSectionIndex+1]);
186 polygon.push_back((*pVertices)[pSectionIndex+2]);
187 polygon.push_back((*pVertices)[pSectionIndex+3]);
188 // G4cout<<"ClipCrossSection: 0-1-2-3"<<G4endl;
189 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
190 return;
191}
192
193//////////////////////////////////////////////////////////////////////////////////
194//
195// Calculate the maximum and minimum extents of the polygons
196// joining the CrossSections at pSectionIndex->pSectionIndex+3 and
197// pSectionIndex+4->pSectionIndex7
198//
199// in the List pVertices, within the boundaries of the voxel limits pVoxelLimit
200//
201// If the minimum is <pMin pMin is set to the new minimum
202// If the maximum is >pMax pMax is set to the new maximum
203//
204// No modifications are made to pVertices
205
206void G4VSolid::ClipBetweenSections( G4ThreeVectorList* pVertices,
207 const G4int pSectionIndex,
208 const G4VoxelLimits& pVoxelLimit,
209 const EAxis pAxis,
210 G4double& pMin, G4double& pMax) const
211{
212 G4ThreeVectorList polygon;
213 polygon.reserve(4);
214 polygon.push_back((*pVertices)[pSectionIndex]);
215 polygon.push_back((*pVertices)[pSectionIndex+4]);
216 polygon.push_back((*pVertices)[pSectionIndex+5]);
217 polygon.push_back((*pVertices)[pSectionIndex+1]);
218 // G4cout<<"ClipBetweenSections: 0-4-5-1"<<G4endl;
219 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
220 polygon.clear();
221
222 polygon.push_back((*pVertices)[pSectionIndex+1]);
223 polygon.push_back((*pVertices)[pSectionIndex+5]);
224 polygon.push_back((*pVertices)[pSectionIndex+6]);
225 polygon.push_back((*pVertices)[pSectionIndex+2]);
226 // G4cout<<"ClipBetweenSections: 1-5-6-2"<<G4endl;
227 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
228 polygon.clear();
229
230 polygon.push_back((*pVertices)[pSectionIndex+2]);
231 polygon.push_back((*pVertices)[pSectionIndex+6]);
232 polygon.push_back((*pVertices)[pSectionIndex+7]);
233 polygon.push_back((*pVertices)[pSectionIndex+3]);
234 // G4cout<<"ClipBetweenSections: 2-6-7-3"<<G4endl;
235 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
236 polygon.clear();
237
238 polygon.push_back((*pVertices)[pSectionIndex+3]);
239 polygon.push_back((*pVertices)[pSectionIndex+7]);
240 polygon.push_back((*pVertices)[pSectionIndex+4]);
241 polygon.push_back((*pVertices)[pSectionIndex]);
242 // G4cout<<"ClipBetweenSections: 3-7-4-0"<<G4endl;
243 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
244 return;
245}
246
247
248///////////////////////////////////////////////////////////////////////////////
249//
250// Calculate the maximum and minimum extents of the convex polygon pPolygon
251// along the axis pAxis, within the limits pVoxelLimit
252//
253
254void
255G4VSolid::CalculateClippedPolygonExtent(G4ThreeVectorList& pPolygon,
256 const G4VoxelLimits& pVoxelLimit,
257 const EAxis pAxis,
258 G4double& pMin,
259 G4double& pMax) const
260{
261 G4int noLeft,i;
262 G4double component;
263 /*
264 G4cout<<G4endl;
265 for(i = 0 ; i < pPolygon.size() ; i++ )
266 {
267 G4cout << i << "\t"
268 << "p.x = " << pPolygon[i].operator()(pAxis) << "\t"
269 // << "p.y = " << pPolygon[i].y() << "\t"
270 // << "p.z = " << pPolygon[i].z() << "\t"
271 << G4endl;
272 }
273 G4cout<<G4endl;
274 */
275 ClipPolygon(pPolygon,pVoxelLimit,pAxis);
276 noLeft = pPolygon.size();
277
278 if ( noLeft )
279 {
280 // G4cout<<G4endl;
281 for (i=0;i<noLeft;i++)
282 {
283 component = pPolygon[i].operator()(pAxis);
284 // G4cout <<i<<"\t"<<component<<G4endl;
285
286 if (component < pMin)
287 {
288 // G4cout <<i<<"\t"<<"Pmin = "<<component<<G4endl;
289 pMin = component;
290 }
291 if (component > pMax)
292 {
293 // G4cout <<i<<"\t"<<"PMax = "<<component<<G4endl;
294 pMax = component;
295 }
296 }
297 // G4cout<<G4endl;
298 }
299 // G4cout<<"pMin = "<<pMin<<"\t"<<"pMax = "<<pMax<<G4endl;
300}
301
302/////////////////////////////////////////////////////////////////////////////
303//
304// Clip the convex polygon described by the vertices at
305// pSectionIndex ->pSectionIndex+3 within pVertices to the limits pVoxelLimit
306//
307// Set pMin to the smallest
308//
309// Calculate the extent of the polygon along pAxis, when clipped to the
310// limits pVoxelLimit. If the polygon exists after clippin, set pMin to
311// the polygon's minimum extent along the axis if <pMin, and set pMax to
312// the polygon's maximum extent along the axis if >pMax.
313//
314// The polygon is described by a set of vectors, where each vector represents
315// a vertex, so that the polygon is described by the vertex sequence:
316// 0th->1st 1st->2nd 2nd->... nth->0th
317//
318// Modifications to the polygon are made
319//
320// NOTE: Execessive copying during clipping
321
322void G4VSolid::ClipPolygon( G4ThreeVectorList& pPolygon,
323 const G4VoxelLimits& pVoxelLimit,
324 const EAxis ) const
325{
326 G4ThreeVectorList outputPolygon;
327
328 if ( pVoxelLimit.IsLimited() )
329 {
330 if (pVoxelLimit.IsXLimited() ) // && pAxis != kXAxis)
331 {
332 G4VoxelLimits simpleLimit1;
333 simpleLimit1.AddLimit(kXAxis,pVoxelLimit.GetMinXExtent(),kInfinity);
334 // G4cout<<"MinXExtent()"<<G4endl;
335 ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1);
336
337 pPolygon.clear();
338
339 if ( !outputPolygon.size() ) return;
340
341 G4VoxelLimits simpleLimit2;
342 // G4cout<<"MaxXExtent()"<<G4endl;
343 simpleLimit2.AddLimit(kXAxis,-kInfinity,pVoxelLimit.GetMaxXExtent());
344 ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
345
346 if ( !pPolygon.size() ) return;
347 else outputPolygon.clear();
348 }
349 if ( pVoxelLimit.IsYLimited() ) // && pAxis != kYAxis)
350 {
351 G4VoxelLimits simpleLimit1;
352 simpleLimit1.AddLimit(kYAxis,pVoxelLimit.GetMinYExtent(),kInfinity);
353 ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1);
354
355 // Must always clear pPolygon - for clip to simpleLimit2 and in case of
356 // early exit
357
358 pPolygon.clear();
359
360 if ( !outputPolygon.size() ) return;
361
362 G4VoxelLimits simpleLimit2;
363 simpleLimit2.AddLimit(kYAxis,-kInfinity,pVoxelLimit.GetMaxYExtent());
364 ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
365
366 if ( !pPolygon.size() ) return;
367 else outputPolygon.clear();
368 }
369 if ( pVoxelLimit.IsZLimited() ) // && pAxis != kZAxis)
370 {
371 G4VoxelLimits simpleLimit1;
372 simpleLimit1.AddLimit(kZAxis,pVoxelLimit.GetMinZExtent(),kInfinity);
373 ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1);
374
375 // Must always clear pPolygon - for clip to simpleLimit2 and in case of
376 // early exit
377
378 pPolygon.clear();
379
380 if ( !outputPolygon.size() ) return;
381
382 G4VoxelLimits simpleLimit2;
383 simpleLimit2.AddLimit(kZAxis,-kInfinity,pVoxelLimit.GetMaxZExtent());
384 ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
385
386 // Return after final clip - no cleanup
387 }
388 }
389}
390
391////////////////////////////////////////////////////////////////////////////
392//
393// pVoxelLimits must be only limited along one axis, and either the maximum
394// along the axis must be +kInfinity, or the minimum -kInfinity
395
396void
397G4VSolid::ClipPolygonToSimpleLimits( G4ThreeVectorList& pPolygon,
398 G4ThreeVectorList& outputPolygon,
399 const G4VoxelLimits& pVoxelLimit ) const
400{
401 G4int i;
402 G4int noVertices=pPolygon.size();
403 G4ThreeVector vEnd,vStart;
404
405 for (i = 0 ; i < noVertices ; i++ )
406 {
407 vStart = pPolygon[i];
408 // G4cout << "i = " << i << G4endl;
409 if ( i == noVertices-1 ) vEnd = pPolygon[0];
410 else vEnd = pPolygon[i+1];
411
412 if ( pVoxelLimit.Inside(vStart) )
413 {
414 if (pVoxelLimit.Inside(vEnd))
415 {
416 // vStart and vEnd inside -> output end point
417 //
418 outputPolygon.push_back(vEnd);
419 }
420 else
421 {
422 // vStart inside, vEnd outside -> output crossing point
423 //
424 // G4cout << "vStart inside, vEnd outside" << G4endl;
425 pVoxelLimit.ClipToLimits(vStart,vEnd);
426 outputPolygon.push_back(vEnd);
427 }
428 }
429 else
430 {
431 if (pVoxelLimit.Inside(vEnd))
432 {
433 // vStart outside, vEnd inside -> output inside section
434 //
435 // G4cout << "vStart outside, vEnd inside" << G4endl;
436 pVoxelLimit.ClipToLimits(vStart,vEnd);
437 outputPolygon.push_back(vStart);
438 outputPolygon.push_back(vEnd);
439 }
440 else // Both point outside -> no output
441 {
442 // outputPolygon.push_back(vStart);
443 // outputPolygon.push_back(vEnd);
444 }
445 }
446 }
447}
448
449const G4VSolid* G4VSolid::GetConstituentSolid(G4int) const
450{ return 0; }
451
452G4VSolid* G4VSolid::GetConstituentSolid(G4int)
453{ return 0; }
454
455const G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr() const
456{ return 0; }
457
458G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr()
459{ return 0; }
460
461G4VisExtent G4VSolid::GetExtent () const
462{
463 G4VisExtent extent;
464 G4VoxelLimits voxelLimits; // Defaults to "infinite" limits.
465 G4AffineTransform affineTransform;
466 G4double vmin, vmax;
467 CalculateExtent(kXAxis,voxelLimits,affineTransform,vmin,vmax);
468 extent.SetXmin (vmin);
469 extent.SetXmax (vmax);
470 CalculateExtent(kYAxis,voxelLimits,affineTransform,vmin,vmax);
471 extent.SetYmin (vmin);
472 extent.SetYmax (vmax);
473 CalculateExtent(kZAxis,voxelLimits,affineTransform,vmin,vmax);
474 extent.SetZmin (vmin);
475 extent.SetZmax (vmax);
476 return extent;
477}
478
479G4Polyhedron* G4VSolid::CreatePolyhedron () const
480{
481 return 0;
482}
483
484G4NURBS* G4VSolid::CreateNURBS () const
485{
486 return 0;
487}
488
489G4Polyhedron* G4VSolid::GetPolyhedron () const
490{
491 return 0;
492}
493
494////////////////////////////////////////////////////////////////
495//
496// Returns an estimation of the solid volume in internal units.
497// The number of statistics and error accuracy is fixed.
498// This method may be overloaded by derived classes to compute the
499// exact geometrical quantity for solids where this is possible.
500// or anyway to cache the computed value.
501// This implementation does NOT cache the computed value.
502
503G4double G4VSolid::GetCubicVolume()
504{
505 G4int cubVolStatistics = 1000000;
506 G4double cubVolEpsilon = 0.001;
507 return EstimateCubicVolume(cubVolStatistics, cubVolEpsilon);
508}
509
510////////////////////////////////////////////////////////////////
511//
512// Calculate cubic volume based on Inside() method.
513// Accuracy is limited by the second argument or the statistics
514// expressed by the first argument.
515// Implementation is courtesy of Vasiliki Despoina Mitsou,
516// University of Athens.
517
518G4double G4VSolid::EstimateCubicVolume(G4int nStat, G4double epsilon) const
519{
520 G4int iInside=0;
521 G4double px,py,pz,minX,maxX,minY,maxY,minZ,maxZ,volume;
522 G4bool yesno;
523 G4ThreeVector p;
524 EInside in;
525
526 // values needed for CalculateExtent signature
527
528 G4VoxelLimits limit; // Unlimited
529 G4AffineTransform origin;
530
531 // min max extents of pSolid along X,Y,Z
532
533 yesno = this->CalculateExtent(kXAxis,limit,origin,minX,maxX);
534 yesno = this->CalculateExtent(kYAxis,limit,origin,minY,maxY);
535 yesno = this->CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
536
537 // limits
538
539 if(nStat < 100) nStat = 100;
540 if(epsilon > 0.01) epsilon = 0.01;
541
542 for(G4int i = 0; i < nStat; i++ )
543 {
544 px = minX+(maxX-minX)*G4UniformRand();
545 py = minY+(maxY-minY)*G4UniformRand();
546 pz = minZ+(maxZ-minZ)*G4UniformRand();
547 p = G4ThreeVector(px,py,pz);
548 in = this->Inside(p);
549 if(in != kOutside) iInside++;
550 }
551 volume = (maxX-minX)*(maxY-minY)*(maxZ-minZ)*iInside/nStat;
552 return volume;
553}
554
555////////////////////////////////////////////////////////////////
556//
557// Returns an estimation of the solid surface area in internal units.
558// The number of statistics and error accuracy is fixed.
559// This method may be overloaded by derived classes to compute the
560// exact geometrical quantity for solids where this is possible.
561// or anyway to cache the computed value.
562// This implementation does NOT cache the computed value.
563
564G4double G4VSolid::GetSurfaceArea()
565{
566 G4int stat = 1000000;
567 G4double ell = -1.;
568 return EstimateSurfaceArea(stat,ell);
569}
570
571////////////////////////////////////////////////////////////////
572//
573// Estimate surface area based on Inside(), DistanceToIn(), and
574// DistanceToOut() methods. Accuracy is limited by the statistics
575// defined by the first argument. Implemented by Mikhail Kosov.
576
577G4double G4VSolid::EstimateSurfaceArea(G4int nStat, G4double ell) const
578{
579 G4int inside=0;
580 G4double px,py,pz,minX,maxX,minY,maxY,minZ,maxZ,surf;
581 G4bool yesno;
582 G4ThreeVector p;
583 EInside in;
584
585 // values needed for CalculateExtent signature
586
587 G4VoxelLimits limit; // Unlimited
588 G4AffineTransform origin;
589
590 // min max extents of pSolid along X,Y,Z
591
592 yesno = this->CalculateExtent(kXAxis,limit,origin,minX,maxX);
593 yesno = this->CalculateExtent(kYAxis,limit,origin,minY,maxY);
594 yesno = this->CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
595
596 // limits
597
598 if(nStat < 100) { nStat = 100; }
599
600 G4double dX=maxX-minX;
601 G4double dY=maxY-minY;
602 G4double dZ=maxZ-minZ;
603 if(ell<=0.) // Automatic definition of skin thickness
604 {
605 G4double minval=dX;
606 if(dY<dX) { minval=dY; }
607 if(dZ<minval) { minval=dZ; }
608 ell=.01*minval;
609 }
610
611 G4double dd=2*ell;
612 minX-=ell; minY-=ell; minZ-=ell; dX+=dd; dY+=dd; dZ+=dd;
613
614 for(G4int i = 0; i < nStat; i++ )
615 {
616 px = minX+dX*G4UniformRand();
617 py = minY+dY*G4UniformRand();
618 pz = minZ+dZ*G4UniformRand();
619 p = G4ThreeVector(px,py,pz);
620 in = this->Inside(p);
621 if(in != kOutside)
622 {
623 if (DistanceToOut(p)<ell) { inside++; }
624 }
625 else if(DistanceToIn(p)<ell) { inside++; }
626 }
627 // @@ The conformal correction can be upgraded
628 surf = dX*dY*dZ*inside/dd/nStat;
629 return surf;
630}
Note: See TracBrowser for help on using the repository browser.