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
00045
00046 #include "G4VCSGfaceted.hh"
00047 #include "G4VCSGface.hh"
00048 #include "G4SolidExtentList.hh"
00049
00050 #include "G4VoxelLimits.hh"
00051 #include "G4AffineTransform.hh"
00052
00053 #include "Randomize.hh"
00054
00055 #include "G4Polyhedron.hh"
00056 #include "G4VGraphicsScene.hh"
00057 #include "G4NURBS.hh"
00058 #include "G4NURBSbox.hh"
00059 #include "G4VisExtent.hh"
00060
00061
00062
00063
00064 G4VCSGfaceted::G4VCSGfaceted( const G4String& name )
00065 : G4VSolid(name),
00066 numFace(0), faces(0), fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0),
00067 fStatistics(1000000), fCubVolEpsilon(0.001), fAreaAccuracy(-1.)
00068 {
00069 }
00070
00071
00072
00073
00074
00075
00076 G4VCSGfaceted::G4VCSGfaceted( __void__& a )
00077 : G4VSolid(a),
00078 numFace(0), faces(0), fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0),
00079 fStatistics(1000000), fCubVolEpsilon(0.001), fAreaAccuracy(-1.)
00080 {
00081 }
00082
00083
00084
00085
00086 G4VCSGfaceted::~G4VCSGfaceted()
00087 {
00088 DeleteStuff();
00089 delete fpPolyhedron;
00090 }
00091
00092
00093
00094
00095
00096 G4VCSGfaceted::G4VCSGfaceted( const G4VCSGfaceted &source )
00097 : G4VSolid( source )
00098 {
00099 fStatistics = source.fStatistics;
00100 fCubVolEpsilon = source.fCubVolEpsilon;
00101 fAreaAccuracy = source.fAreaAccuracy;
00102
00103 CopyStuff( source );
00104 }
00105
00106
00107
00108
00109
00110 const G4VCSGfaceted &G4VCSGfaceted::operator=( const G4VCSGfaceted &source )
00111 {
00112 if (&source == this) { return *this; }
00113
00114
00115
00116 G4VSolid::operator=(source);
00117
00118
00119
00120 fStatistics = source.fStatistics;
00121 fCubVolEpsilon = source.fCubVolEpsilon;
00122 fAreaAccuracy = source.fAreaAccuracy;
00123
00124 DeleteStuff();
00125 CopyStuff( source );
00126
00127 return *this;
00128 }
00129
00130
00131
00132
00133
00134
00135
00136 void G4VCSGfaceted::CopyStuff( const G4VCSGfaceted &source )
00137 {
00138 numFace = source.numFace;
00139 if (numFace == 0) { return; }
00140
00141 faces = new G4VCSGface*[numFace];
00142
00143 G4VCSGface **face = faces,
00144 **sourceFace = source.faces;
00145 do
00146 {
00147 *face = (*sourceFace)->Clone();
00148 } while( ++sourceFace, ++face < faces+numFace );
00149 fCubicVolume = source.fCubicVolume;
00150 fSurfaceArea = source.fSurfaceArea;
00151 fpPolyhedron = 0;
00152 }
00153
00154
00155
00156
00157
00158
00159
00160 void G4VCSGfaceted::DeleteStuff()
00161 {
00162 if (numFace)
00163 {
00164 G4VCSGface **face = faces;
00165 do
00166 {
00167 delete *face;
00168 } while( ++face < faces + numFace );
00169
00170 delete [] faces;
00171 }
00172 }
00173
00174
00175
00176
00177
00178 G4bool G4VCSGfaceted::CalculateExtent( const EAxis axis,
00179 const G4VoxelLimits &voxelLimit,
00180 const G4AffineTransform &transform,
00181 G4double &min,
00182 G4double &max ) const
00183 {
00184 G4SolidExtentList extentList( axis, voxelLimit );
00185
00186
00187
00188
00189 G4VCSGface **face = faces;
00190 do
00191 {
00192 (*face)->CalculateExtent( axis, voxelLimit, transform, extentList );
00193 } while( ++face < faces + numFace );
00194
00195
00196
00197
00198 return extentList.GetExtent( min, max );
00199 }
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210 EInside G4VCSGfaceted::Inside( const G4ThreeVector &p ) const
00211 {
00212 EInside answer=kOutside;
00213 G4VCSGface **face = faces;
00214 G4double best = kInfinity;
00215 do
00216 {
00217 G4double distance;
00218 EInside result = (*face)->Inside( p, kCarTolerance/2, &distance );
00219 if (result == kSurface) { return kSurface; }
00220 if (distance < best)
00221 {
00222 best = distance;
00223 answer = result;
00224 }
00225 } while( ++face < faces + numFace );
00226
00227 return answer;
00228 }
00229
00230
00231
00232
00233
00234 G4ThreeVector G4VCSGfaceted::SurfaceNormal( const G4ThreeVector& p ) const
00235 {
00236 G4ThreeVector answer;
00237 G4VCSGface **face = faces;
00238 G4double best = kInfinity;
00239 do
00240 {
00241 G4double distance;
00242 G4ThreeVector normal = (*face)->Normal( p, &distance );
00243 if (distance < best)
00244 {
00245 best = distance;
00246 answer = normal;
00247 }
00248 } while( ++face < faces + numFace );
00249
00250 return answer;
00251 }
00252
00253
00254
00255
00256
00257 G4double G4VCSGfaceted::DistanceToIn( const G4ThreeVector &p,
00258 const G4ThreeVector &v ) const
00259 {
00260 G4double distance = kInfinity;
00261 G4double distFromSurface = kInfinity;
00262 G4VCSGface **face = faces;
00263 G4VCSGface *bestFace = *face;
00264 do
00265 {
00266 G4double faceDistance,
00267 faceDistFromSurface;
00268 G4ThreeVector faceNormal;
00269 G4bool faceAllBehind;
00270 if ((*face)->Intersect( p, v, false, kCarTolerance/2,
00271 faceDistance, faceDistFromSurface,
00272 faceNormal, faceAllBehind ) )
00273 {
00274
00275
00276
00277 if (faceDistance < distance)
00278 {
00279 distance = faceDistance;
00280 distFromSurface = faceDistFromSurface;
00281 bestFace = *face;
00282 if (distFromSurface <= 0) { return 0; }
00283 }
00284 }
00285 } while( ++face < faces + numFace );
00286
00287 if (distance < kInfinity && distFromSurface<kCarTolerance/2)
00288 {
00289 if (bestFace->Distance(p,false) < kCarTolerance/2) { distance = 0; }
00290 }
00291
00292 return distance;
00293 }
00294
00295
00296
00297
00298
00299 G4double G4VCSGfaceted::DistanceToIn( const G4ThreeVector &p ) const
00300 {
00301 return DistanceTo( p, false );
00302 }
00303
00304
00305
00306
00307
00308 G4double G4VCSGfaceted::DistanceToOut( const G4ThreeVector &p,
00309 const G4ThreeVector &v,
00310 const G4bool calcNorm,
00311 G4bool *validNorm,
00312 G4ThreeVector *n ) const
00313 {
00314 G4bool allBehind = true;
00315 G4double distance = kInfinity;
00316 G4double distFromSurface = kInfinity;
00317 G4ThreeVector normal;
00318
00319 G4VCSGface **face = faces;
00320 G4VCSGface *bestFace = *face;
00321 do
00322 {
00323 G4double faceDistance,
00324 faceDistFromSurface;
00325 G4ThreeVector faceNormal;
00326 G4bool faceAllBehind;
00327 if ((*face)->Intersect( p, v, true, kCarTolerance/2,
00328 faceDistance, faceDistFromSurface,
00329 faceNormal, faceAllBehind ) )
00330 {
00331
00332
00333
00334 if ( (distance < kInfinity) || (!faceAllBehind) ) { allBehind = false; }
00335 if (faceDistance < distance)
00336 {
00337 distance = faceDistance;
00338 distFromSurface = faceDistFromSurface;
00339 normal = faceNormal;
00340 bestFace = *face;
00341 if (distFromSurface <= 0) { break; }
00342 }
00343 }
00344 } while( ++face < faces + numFace );
00345
00346 if (distance < kInfinity)
00347 {
00348 if (distFromSurface <= 0)
00349 {
00350 distance = 0;
00351 }
00352 else if (distFromSurface<kCarTolerance/2)
00353 {
00354 if (bestFace->Distance(p,true) < kCarTolerance/2) { distance = 0; }
00355 }
00356
00357 if (calcNorm)
00358 {
00359 *validNorm = allBehind;
00360 *n = normal;
00361 }
00362 }
00363 else
00364 {
00365 if (Inside(p) == kSurface) { distance = 0; }
00366 if (calcNorm) { *validNorm = false; }
00367 }
00368
00369 return distance;
00370 }
00371
00372
00373
00374
00375
00376 G4double G4VCSGfaceted::DistanceToOut( const G4ThreeVector &p ) const
00377 {
00378 return DistanceTo( p, true );
00379 }
00380
00381
00382
00383
00384
00385
00386
00387 G4double G4VCSGfaceted::DistanceTo( const G4ThreeVector &p,
00388 const G4bool outgoing ) const
00389 {
00390 G4VCSGface **face = faces;
00391 G4double best = kInfinity;
00392 do
00393 {
00394 G4double distance = (*face)->Distance( p, outgoing );
00395 if (distance < best) { best = distance; }
00396 } while( ++face < faces + numFace );
00397
00398 return (best < 0.5*kCarTolerance) ? 0 : best;
00399 }
00400
00401
00402
00403
00404
00405 void G4VCSGfaceted::DescribeYourselfTo( G4VGraphicsScene& scene ) const
00406 {
00407 scene.AddSolid( *this );
00408 }
00409
00410
00411
00412
00413
00414
00415
00416 G4VisExtent G4VCSGfaceted::GetExtent() const
00417 {
00418 static const G4ThreeVector xMax(1,0,0), xMin(-1,0,0),
00419 yMax(0,1,0), yMin(0,-1,0),
00420 zMax(0,0,1), zMin(0,0,-1);
00421 static const G4ThreeVector *axes[6] =
00422 { &xMin, &xMax, &yMin, &yMax, &zMin, &zMax };
00423
00424 G4double answers[6] =
00425 {-kInfinity, -kInfinity, -kInfinity, -kInfinity, -kInfinity, -kInfinity};
00426
00427 G4VCSGface **face = faces;
00428 do
00429 {
00430 const G4ThreeVector **axis = axes+5 ;
00431 G4double *answer = answers+5;
00432 do
00433 {
00434 G4double testFace = (*face)->Extent( **axis );
00435 if (testFace > *answer) { *answer = testFace; }
00436 }
00437 while( --axis, --answer >= answers );
00438
00439 } while( ++face < faces + numFace );
00440
00441 return G4VisExtent( -answers[0], answers[1],
00442 -answers[2], answers[3],
00443 -answers[4], answers[5] );
00444 }
00445
00446
00447
00448
00449
00450 G4GeometryType G4VCSGfaceted::GetEntityType() const
00451 {
00452 return G4String("G4CSGfaceted");
00453 }
00454
00455
00456
00457
00458
00459 std::ostream& G4VCSGfaceted::StreamInfo( std::ostream& os ) const
00460 {
00461 os << "-----------------------------------------------------------\n"
00462 << " *** Dump for solid - " << GetName() << " ***\n"
00463 << " ===================================================\n"
00464 << " Solid type: G4VCSGfaceted\n"
00465 << " Parameters: \n"
00466 << " number of faces: " << numFace << "\n"
00467 << "-----------------------------------------------------------\n";
00468
00469 return os;
00470 }
00471
00472
00473
00474
00475
00476 G4int G4VCSGfaceted::GetCubVolStatistics() const
00477 {
00478 return fStatistics;
00479 }
00480
00481
00482
00483
00484
00485 G4double G4VCSGfaceted::GetCubVolEpsilon() const
00486 {
00487 return fCubVolEpsilon;
00488 }
00489
00490
00491
00492
00493
00494 void G4VCSGfaceted::SetCubVolStatistics(G4int st)
00495 {
00496 fCubicVolume=0.;
00497 fStatistics=st;
00498 }
00499
00500
00501
00502
00503
00504 void G4VCSGfaceted::SetCubVolEpsilon(G4double ep)
00505 {
00506 fCubicVolume=0.;
00507 fCubVolEpsilon=ep;
00508 }
00509
00510
00511
00512
00513
00514 G4int G4VCSGfaceted::GetAreaStatistics() const
00515 {
00516 return fStatistics;
00517 }
00518
00519
00520
00521
00522
00523 G4double G4VCSGfaceted::GetAreaAccuracy() const
00524 {
00525 return fAreaAccuracy;
00526 }
00527
00528
00529
00530
00531
00532 void G4VCSGfaceted::SetAreaStatistics(G4int st)
00533 {
00534 fSurfaceArea=0.;
00535 fStatistics=st;
00536 }
00537
00538
00539
00540
00541
00542 void G4VCSGfaceted::SetAreaAccuracy(G4double ep)
00543 {
00544 fSurfaceArea=0.;
00545 fAreaAccuracy=ep;
00546 }
00547
00548
00549
00550
00551
00552 G4double G4VCSGfaceted::GetCubicVolume()
00553 {
00554 if(fCubicVolume != 0.) {;}
00555 else { fCubicVolume = EstimateCubicVolume(fStatistics,fCubVolEpsilon); }
00556 return fCubicVolume;
00557 }
00558
00559
00560
00561
00562
00563 G4double G4VCSGfaceted::GetSurfaceArea()
00564 {
00565 if(fSurfaceArea != 0.) {;}
00566 else { fSurfaceArea = EstimateSurfaceArea(fStatistics,fAreaAccuracy); }
00567 return fSurfaceArea;
00568 }
00569
00570
00571
00572
00573
00574 G4Polyhedron* G4VCSGfaceted::GetPolyhedron () const
00575 {
00576 if (!fpPolyhedron ||
00577 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
00578 fpPolyhedron->GetNumberOfRotationSteps())
00579 {
00580 delete fpPolyhedron;
00581 fpPolyhedron = CreatePolyhedron();
00582 }
00583 return fpPolyhedron;
00584 }
00585
00586
00587
00588
00589
00590
00591 G4ThreeVector G4VCSGfaceted::GetPointOnSurfaceGeneric( ) const
00592 {
00593
00594
00595 G4ThreeVector answer=G4ThreeVector(0.,0.,0.);
00596 G4VCSGface **face = faces;
00597 G4double area = 0;
00598 G4int i;
00599 std::vector<G4double> areas;
00600
00601
00602
00603 do
00604 {
00605 G4double result = (*face)->SurfaceArea( );
00606 areas.push_back(result);
00607 area=area+result;
00608 } while( ++face < faces + numFace );
00609
00610
00611
00612 G4VCSGface **face1 = faces;
00613 G4double chose = area*G4UniformRand();
00614 G4double Achose1, Achose2;
00615 Achose1=0; Achose2=0.;
00616 i=0;
00617
00618 do
00619 {
00620 Achose2+=areas[i];
00621 if(chose>=Achose1 && chose<Achose2)
00622 {
00623 G4ThreeVector point;
00624 point= (*face1)->GetPointOnFace();
00625 return point;
00626 }
00627 i++;
00628 Achose1=Achose2;
00629 } while( ++face1 < faces + numFace );
00630
00631 return answer;
00632 }