00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 #include "G4VSolid.hh"
00045 #include "G4SolidStore.hh"
00046 #include "globals.hh"
00047 #include "Randomize.hh"
00048 #include "G4GeometryTolerance.hh"
00049
00050 #include "G4VoxelLimits.hh"
00051 #include "G4AffineTransform.hh"
00052 #include "G4VisExtent.hh"
00053
00055
00056
00057
00058
00059
00060 G4VSolid::G4VSolid(const G4String& name)
00061 : fshapeName(name)
00062 {
00063 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00064
00065
00066
00067 G4SolidStore::GetInstance()->Register(this);
00068 }
00069
00071
00072
00073
00074
00075 G4VSolid::G4VSolid(const G4VSolid& rhs)
00076 : kCarTolerance(rhs.kCarTolerance), fshapeName(rhs.fshapeName)
00077 {
00078
00079
00080 G4SolidStore::GetInstance()->Register(this);
00081 }
00082
00084
00085
00086
00087
00088 G4VSolid::G4VSolid( __void__& )
00089 : fshapeName("")
00090 {
00091
00092
00093 G4SolidStore::GetInstance()->Register(this);
00094 }
00095
00097
00098
00099
00100
00101 G4VSolid::~G4VSolid()
00102 {
00103 G4SolidStore::GetInstance()->DeRegister(this);
00104 }
00105
00107
00108
00109
00110 G4VSolid& G4VSolid::operator = (const G4VSolid& rhs)
00111 {
00112
00113
00114 if (this == &rhs) { return *this; }
00115
00116
00117
00118 kCarTolerance = rhs.kCarTolerance;
00119 fshapeName = rhs.fshapeName;
00120
00121 return *this;
00122 }
00123
00125
00126
00127
00128 std::ostream& operator<< ( std::ostream& os, const G4VSolid& e )
00129 {
00130 return e.StreamInfo(os);
00131 }
00132
00134
00135
00136
00137 void G4VSolid::ComputeDimensions(G4VPVParameterisation*,
00138 const G4int,
00139 const G4VPhysicalVolume*)
00140 {
00141 std::ostringstream message;
00142 message << "Illegal call to G4VSolid::ComputeDimensions()" << G4endl
00143 << "Method not overloaded by derived class !";
00144 G4Exception("G4VSolid::ComputeDimensions()", "GeomMgt0003",
00145 FatalException, message);
00146 }
00147
00149
00150
00151
00152 G4ThreeVector G4VSolid::GetPointOnSurface() const
00153 {
00154 std::ostringstream message;
00155 message << "Not implemented for solid: "
00156 << this->GetEntityType() << " !" << G4endl
00157 << "Returning origin.";
00158 G4Exception("G4VSolid::GetPointOnSurface()", "GeomMgt1001",
00159 JustWarning, message);
00160 return G4ThreeVector(0,0,0);
00161 }
00162
00164
00165
00166
00167 const G4VSolid* G4VSolid::GetConstituentSolid(G4int) const
00168 { return 0; }
00169
00170 G4VSolid* G4VSolid::GetConstituentSolid(G4int)
00171 { return 0; }
00172
00173 const G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr() const
00174 { return 0; }
00175
00176 G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr()
00177 { return 0; }
00178
00180
00181
00182
00183
00184
00185
00186
00187
00188 G4double G4VSolid::GetCubicVolume()
00189 {
00190 G4int cubVolStatistics = 1000000;
00191 G4double cubVolEpsilon = 0.001;
00192 return EstimateCubicVolume(cubVolStatistics, cubVolEpsilon);
00193 }
00194
00196
00197
00198
00199
00200
00201
00202
00203 G4double G4VSolid::EstimateCubicVolume(G4int nStat, G4double epsilon) const
00204 {
00205 G4int iInside=0;
00206 G4double px,py,pz,minX,maxX,minY,maxY,minZ,maxZ,volume;
00207 G4ThreeVector p;
00208 EInside in;
00209
00210
00211
00212 G4VoxelLimits limit;
00213 G4AffineTransform origin;
00214
00215
00216
00217 this->CalculateExtent(kXAxis,limit,origin,minX,maxX);
00218 this->CalculateExtent(kYAxis,limit,origin,minY,maxY);
00219 this->CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
00220
00221
00222
00223 if(nStat < 100) nStat = 100;
00224 if(epsilon > 0.01) epsilon = 0.01;
00225
00226 for(G4int i = 0; i < nStat; i++ )
00227 {
00228 px = minX+(maxX-minX)*G4UniformRand();
00229 py = minY+(maxY-minY)*G4UniformRand();
00230 pz = minZ+(maxZ-minZ)*G4UniformRand();
00231 p = G4ThreeVector(px,py,pz);
00232 in = this->Inside(p);
00233 if(in != kOutside) iInside++;
00234 }
00235 volume = (maxX-minX)*(maxY-minY)*(maxZ-minZ)*iInside/nStat;
00236 return volume;
00237 }
00238
00240
00241
00242
00243
00244
00245
00246
00247
00248 G4double G4VSolid::GetSurfaceArea()
00249 {
00250 G4int stat = 1000000;
00251 G4double ell = -1.;
00252 return EstimateSurfaceArea(stat,ell);
00253 }
00254
00256
00257
00258
00259
00260
00261 G4double G4VSolid::EstimateSurfaceArea(G4int nStat, G4double ell) const
00262 {
00263 G4int inside=0;
00264 G4double px,py,pz,minX,maxX,minY,maxY,minZ,maxZ,surf;
00265 G4ThreeVector p;
00266 EInside in;
00267
00268
00269
00270 G4VoxelLimits limit;
00271 G4AffineTransform origin;
00272
00273
00274
00275 this->CalculateExtent(kXAxis,limit,origin,minX,maxX);
00276 this->CalculateExtent(kYAxis,limit,origin,minY,maxY);
00277 this->CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
00278
00279
00280
00281 if(nStat < 100) { nStat = 100; }
00282
00283 G4double dX=maxX-minX;
00284 G4double dY=maxY-minY;
00285 G4double dZ=maxZ-minZ;
00286 if(ell<=0.)
00287 {
00288 G4double minval=dX;
00289 if(dY<dX) { minval=dY; }
00290 if(dZ<minval) { minval=dZ; }
00291 ell=.01*minval;
00292 }
00293
00294 G4double dd=2*ell;
00295 minX-=ell; minY-=ell; minZ-=ell; dX+=dd; dY+=dd; dZ+=dd;
00296
00297 for(G4int i = 0; i < nStat; i++ )
00298 {
00299 px = minX+dX*G4UniformRand();
00300 py = minY+dY*G4UniformRand();
00301 pz = minZ+dZ*G4UniformRand();
00302 p = G4ThreeVector(px,py,pz);
00303 in = this->Inside(p);
00304 if(in != kOutside)
00305 {
00306 if (DistanceToOut(p)<ell) { inside++; }
00307 }
00308 else if(DistanceToIn(p)<ell) { inside++; }
00309 }
00310
00311 surf = dX*dY*dZ*inside/dd/nStat;
00312 return surf;
00313 }
00314
00316
00317
00318
00319
00320
00321
00322 G4VSolid* G4VSolid::Clone() const
00323 {
00324 std::ostringstream message;
00325 message << "Clone() method not implemented for type: "
00326 << GetEntityType() << "!" << G4endl
00327 << "Returning NULL pointer!";
00328 G4Exception("G4VSolid::Clone()", "GeomMgt1001", JustWarning, message);
00329 return 0;
00330 }
00331
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345 void G4VSolid::ClipCrossSection( G4ThreeVectorList* pVertices,
00346 const G4int pSectionIndex,
00347 const G4VoxelLimits& pVoxelLimit,
00348 const EAxis pAxis,
00349 G4double& pMin, G4double& pMax) const
00350 {
00351
00352 G4ThreeVectorList polygon;
00353 polygon.reserve(4);
00354 polygon.push_back((*pVertices)[pSectionIndex]);
00355 polygon.push_back((*pVertices)[pSectionIndex+1]);
00356 polygon.push_back((*pVertices)[pSectionIndex+2]);
00357 polygon.push_back((*pVertices)[pSectionIndex+3]);
00358
00359 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
00360 return;
00361 }
00362
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376 void G4VSolid::ClipBetweenSections( G4ThreeVectorList* pVertices,
00377 const G4int pSectionIndex,
00378 const G4VoxelLimits& pVoxelLimit,
00379 const EAxis pAxis,
00380 G4double& pMin, G4double& pMax) const
00381 {
00382 G4ThreeVectorList polygon;
00383 polygon.reserve(4);
00384 polygon.push_back((*pVertices)[pSectionIndex]);
00385 polygon.push_back((*pVertices)[pSectionIndex+4]);
00386 polygon.push_back((*pVertices)[pSectionIndex+5]);
00387 polygon.push_back((*pVertices)[pSectionIndex+1]);
00388
00389 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
00390 polygon.clear();
00391
00392 polygon.push_back((*pVertices)[pSectionIndex+1]);
00393 polygon.push_back((*pVertices)[pSectionIndex+5]);
00394 polygon.push_back((*pVertices)[pSectionIndex+6]);
00395 polygon.push_back((*pVertices)[pSectionIndex+2]);
00396
00397 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
00398 polygon.clear();
00399
00400 polygon.push_back((*pVertices)[pSectionIndex+2]);
00401 polygon.push_back((*pVertices)[pSectionIndex+6]);
00402 polygon.push_back((*pVertices)[pSectionIndex+7]);
00403 polygon.push_back((*pVertices)[pSectionIndex+3]);
00404
00405 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
00406 polygon.clear();
00407
00408 polygon.push_back((*pVertices)[pSectionIndex+3]);
00409 polygon.push_back((*pVertices)[pSectionIndex+7]);
00410 polygon.push_back((*pVertices)[pSectionIndex+4]);
00411 polygon.push_back((*pVertices)[pSectionIndex]);
00412
00413 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
00414 return;
00415 }
00416
00417
00419
00420
00421
00422
00423
00424 void
00425 G4VSolid::CalculateClippedPolygonExtent(G4ThreeVectorList& pPolygon,
00426 const G4VoxelLimits& pVoxelLimit,
00427 const EAxis pAxis,
00428 G4double& pMin,
00429 G4double& pMax) const
00430 {
00431 G4int noLeft,i;
00432 G4double component;
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445 ClipPolygon(pPolygon,pVoxelLimit,pAxis);
00446 noLeft = pPolygon.size();
00447
00448 if ( noLeft )
00449 {
00450
00451 for (i=0;i<noLeft;i++)
00452 {
00453 component = pPolygon[i].operator()(pAxis);
00454
00455
00456 if (component < pMin)
00457 {
00458
00459 pMin = component;
00460 }
00461 if (component > pMax)
00462 {
00463
00464 pMax = component;
00465 }
00466 }
00467
00468 }
00469
00470 }
00471
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492 void G4VSolid::ClipPolygon( G4ThreeVectorList& pPolygon,
00493 const G4VoxelLimits& pVoxelLimit,
00494 const EAxis ) const
00495 {
00496 G4ThreeVectorList outputPolygon;
00497
00498 if ( pVoxelLimit.IsLimited() )
00499 {
00500 if (pVoxelLimit.IsXLimited() )
00501 {
00502 G4VoxelLimits simpleLimit1;
00503 simpleLimit1.AddLimit(kXAxis,pVoxelLimit.GetMinXExtent(),kInfinity);
00504
00505 ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1);
00506
00507 pPolygon.clear();
00508
00509 if ( !outputPolygon.size() ) return;
00510
00511 G4VoxelLimits simpleLimit2;
00512
00513 simpleLimit2.AddLimit(kXAxis,-kInfinity,pVoxelLimit.GetMaxXExtent());
00514 ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
00515
00516 if ( !pPolygon.size() ) return;
00517 else outputPolygon.clear();
00518 }
00519 if ( pVoxelLimit.IsYLimited() )
00520 {
00521 G4VoxelLimits simpleLimit1;
00522 simpleLimit1.AddLimit(kYAxis,pVoxelLimit.GetMinYExtent(),kInfinity);
00523 ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1);
00524
00525
00526
00527
00528 pPolygon.clear();
00529
00530 if ( !outputPolygon.size() ) return;
00531
00532 G4VoxelLimits simpleLimit2;
00533 simpleLimit2.AddLimit(kYAxis,-kInfinity,pVoxelLimit.GetMaxYExtent());
00534 ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
00535
00536 if ( !pPolygon.size() ) return;
00537 else outputPolygon.clear();
00538 }
00539 if ( pVoxelLimit.IsZLimited() )
00540 {
00541 G4VoxelLimits simpleLimit1;
00542 simpleLimit1.AddLimit(kZAxis,pVoxelLimit.GetMinZExtent(),kInfinity);
00543 ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1);
00544
00545
00546
00547
00548 pPolygon.clear();
00549
00550 if ( !outputPolygon.size() ) return;
00551
00552 G4VoxelLimits simpleLimit2;
00553 simpleLimit2.AddLimit(kZAxis,-kInfinity,pVoxelLimit.GetMaxZExtent());
00554 ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
00555
00556
00557 }
00558 }
00559 }
00560
00562
00563
00564
00565
00566 void
00567 G4VSolid::ClipPolygonToSimpleLimits( G4ThreeVectorList& pPolygon,
00568 G4ThreeVectorList& outputPolygon,
00569 const G4VoxelLimits& pVoxelLimit ) const
00570 {
00571 G4int i;
00572 G4int noVertices=pPolygon.size();
00573 G4ThreeVector vEnd,vStart;
00574
00575 for (i = 0 ; i < noVertices ; i++ )
00576 {
00577 vStart = pPolygon[i];
00578
00579 if ( i == noVertices-1 ) vEnd = pPolygon[0];
00580 else vEnd = pPolygon[i+1];
00581
00582 if ( pVoxelLimit.Inside(vStart) )
00583 {
00584 if (pVoxelLimit.Inside(vEnd))
00585 {
00586
00587
00588 outputPolygon.push_back(vEnd);
00589 }
00590 else
00591 {
00592
00593
00594
00595 pVoxelLimit.ClipToLimits(vStart,vEnd);
00596 outputPolygon.push_back(vEnd);
00597 }
00598 }
00599 else
00600 {
00601 if (pVoxelLimit.Inside(vEnd))
00602 {
00603
00604
00605
00606 pVoxelLimit.ClipToLimits(vStart,vEnd);
00607 outputPolygon.push_back(vStart);
00608 outputPolygon.push_back(vEnd);
00609 }
00610 else
00611 {
00612
00613
00614 }
00615 }
00616 }
00617 }
00618
00619 G4VisExtent G4VSolid::GetExtent () const
00620 {
00621 G4VisExtent extent;
00622 G4VoxelLimits voxelLimits;
00623 G4AffineTransform affineTransform;
00624 G4double vmin, vmax;
00625 CalculateExtent(kXAxis,voxelLimits,affineTransform,vmin,vmax);
00626 extent.SetXmin (vmin);
00627 extent.SetXmax (vmax);
00628 CalculateExtent(kYAxis,voxelLimits,affineTransform,vmin,vmax);
00629 extent.SetYmin (vmin);
00630 extent.SetYmax (vmax);
00631 CalculateExtent(kZAxis,voxelLimits,affineTransform,vmin,vmax);
00632 extent.SetZmin (vmin);
00633 extent.SetZmax (vmax);
00634 return extent;
00635 }
00636
00637 G4Polyhedron* G4VSolid::CreatePolyhedron () const
00638 {
00639 return 0;
00640 }
00641
00642 G4NURBS* G4VSolid::CreateNURBS () const
00643 {
00644 return 0;
00645 }
00646
00647 G4Polyhedron* G4VSolid::GetPolyhedron () const
00648 {
00649 return 0;
00650 }