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 #include "G4BREPSolid.hh"
00037 #include "G4VoxelLimits.hh"
00038 #include "G4AffineTransform.hh"
00039 #include "G4VGraphicsScene.hh"
00040 #include "G4Polyhedron.hh"
00041 #include "G4NURBSbox.hh"
00042 #include "G4BoundingBox3D.hh"
00043 #include "G4FPlane.hh"
00044 #include "G4BSplineSurface.hh"
00045 #include "G4ToroidalSurface.hh"
00046 #include "G4SphericalSurface.hh"
00047
00048 G4Ray G4BREPSolid::Track;
00049 G4double G4BREPSolid::ShortestDistance= kInfinity;
00050 G4int G4BREPSolid::NumberOfSolids=0;
00051
00052 G4BREPSolid::G4BREPSolid(const G4String& name)
00053 : G4VSolid(name),
00054 Box(0), Convex(0), AxisBox(0), PlaneSolid(0), place(0), bbox(0),
00055 intersectionDistance(kInfinity), active(1), startInside(0),
00056 nb_of_surfaces(0), SurfaceVec(0), RealDist(0.), solidname(name), Id(0),
00057 fStatistics(1000000), fCubVolEpsilon(0.001), fAreaAccuracy(-1.),
00058 fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
00059 {
00060 static G4bool warn=true;
00061 if (warn)
00062 {
00063 G4cout
00064 << "--------------------------------------------------------" << G4endl
00065 << "WARNING: BREPS classes are being dismissed. They will |" << G4endl
00066 << " be removed starting from next Geant4 major |" << G4endl
00067 << " release. Please, consider switching to adopt |" << G4endl
00068 << " correspondent CSG or specific primitives. |" << G4endl
00069 << "--------------------------------------------------------"
00070 << G4endl;
00071 warn = false;
00072 }
00073 }
00074
00075 G4BREPSolid::G4BREPSolid( const G4String& name ,
00076 G4Surface** srfVec ,
00077 G4int numberOfSrfs )
00078 : G4VSolid(name),
00079 Box(0), Convex(0), AxisBox(0), PlaneSolid(0), place(0), bbox(0),
00080 intersectionDistance(kInfinity), active(1), startInside(0),
00081 nb_of_surfaces(numberOfSrfs), SurfaceVec(srfVec), RealDist(0.),
00082 solidname(name), Id(0),
00083 fStatistics(1000000), fCubVolEpsilon(0.001), fAreaAccuracy(-1.),
00084 fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
00085 {
00086 static G4bool warn=true;
00087 if (warn)
00088 {
00089 G4cout
00090 << "--------------------------------------------------------" << G4endl
00091 << "WARNING: BREPS classes are being dismissed. They will |" << G4endl
00092 << " be removed starting from next Geant4 major |" << G4endl
00093 << " release. Please, consider switching to adopt |" << G4endl
00094 << " correspondent CSG or specific primitives. |" << G4endl
00095 << "--------------------------------------------------------"
00096 << G4endl;
00097 warn = false;
00098 }
00099 Initialize();
00100 }
00101
00102 G4BREPSolid::G4BREPSolid( __void__& a )
00103 : G4VSolid(a),
00104 Box(0), Convex(0), AxisBox(0), PlaneSolid(0), place(0), bbox(0),
00105 intersectionDistance(kInfinity), active(1), startInside(0),
00106 nb_of_surfaces(0), SurfaceVec(0), RealDist(0.), solidname(""), Id(0),
00107 fStatistics(1000000), fCubVolEpsilon(0.001), fAreaAccuracy(-1.),
00108 fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
00109 {
00110 }
00111
00112 G4BREPSolid::~G4BREPSolid()
00113 {
00114 delete place;
00115 delete bbox;
00116
00117 for(G4int a=0;a<nb_of_surfaces;++a)
00118 { delete SurfaceVec[a]; }
00119
00120 delete [] SurfaceVec;
00121 delete fpPolyhedron;
00122 }
00123
00124 G4BREPSolid::G4BREPSolid(const G4BREPSolid& rhs)
00125 : G4VSolid(rhs), Box(rhs.Box), Convex(rhs.Convex), AxisBox(rhs.AxisBox),
00126 PlaneSolid(rhs.PlaneSolid), place(0), bbox(0),
00127 intersectionDistance(rhs.intersectionDistance), active(rhs.active),
00128 startInside(rhs.startInside), nb_of_surfaces(rhs.nb_of_surfaces),
00129 intersection_point(rhs.intersection_point), SurfaceVec(0),
00130 RealDist(rhs.RealDist), solidname(rhs.solidname), Id(rhs.Id),
00131 fStatistics(rhs.fStatistics), fCubVolEpsilon(rhs.fCubVolEpsilon),
00132 fAreaAccuracy(rhs.fAreaAccuracy), fCubicVolume(rhs.fCubicVolume),
00133 fSurfaceArea(rhs.fSurfaceArea), fpPolyhedron(0)
00134 {
00135 Initialize();
00136 }
00137
00138 G4BREPSolid& G4BREPSolid::operator = (const G4BREPSolid& rhs)
00139 {
00140
00141
00142 if (this == &rhs) { return *this; }
00143
00144
00145
00146 G4VSolid::operator=(rhs);
00147
00148
00149
00150 Box= rhs.Box; Convex= rhs.Convex; AxisBox= rhs.AxisBox;
00151 PlaneSolid= rhs.PlaneSolid;
00152 intersectionDistance= rhs.intersectionDistance; active= rhs.active;
00153 startInside= rhs.startInside; nb_of_surfaces= rhs.nb_of_surfaces;
00154 intersection_point= rhs.intersection_point;
00155 RealDist= rhs.RealDist; solidname= rhs.solidname; Id= rhs.Id;
00156 fStatistics= rhs.fStatistics; fCubVolEpsilon= rhs.fCubVolEpsilon;
00157 fAreaAccuracy= rhs.fAreaAccuracy; fCubicVolume= rhs.fCubicVolume;
00158 fSurfaceArea= rhs.fSurfaceArea;
00159
00160 delete place; delete bbox; place= 0; bbox= 0;
00161 for(G4int a=0;a<nb_of_surfaces;++a)
00162 { delete SurfaceVec[a]; }
00163 delete [] SurfaceVec; SurfaceVec= 0;
00164 delete fpPolyhedron; fpPolyhedron= 0;
00165
00166 Initialize();
00167
00168 return *this;
00169 }
00170
00171 void G4BREPSolid::Initialize()
00172 {
00173 if(active)
00174 {
00175
00176
00177
00178 ShortestDistance= kInfinity;
00179 if (!SurfaceVec) { return; }
00180
00181 IsBox();
00182 CheckSurfaceNormals();
00183
00184 if(!Box || !AxisBox)
00185 IsConvex();
00186
00187 CalcBBoxes();
00188 }
00189 }
00190
00191 G4String G4BREPSolid::GetEntityType() const
00192 {
00193 return "Closed_Shell";
00194 }
00195
00196 G4VSolid* G4BREPSolid::Clone() const
00197 {
00198 return new G4BREPSolid(*this);
00199 }
00200
00201 void G4BREPSolid::Reset() const
00202 {
00203 ((G4BREPSolid*)this)->active=1;
00204 ((G4BREPSolid*)this)->intersectionDistance=kInfinity;
00205 ((G4BREPSolid*)this)->startInside=0;
00206
00207 for(register G4int a=0;a<nb_of_surfaces;a++)
00208 SurfaceVec[a]->Reset();
00209
00210 ShortestDistance = kInfinity;
00211 }
00212
00213 void G4BREPSolid::CheckSurfaceNormals()
00214 {
00215 if(!PlaneSolid)
00216 return;
00217
00218 Convex=1;
00219
00220
00221
00222
00223
00224
00225 G4Surface* srf;
00226 G4Point3D V;
00227
00228 G4int SrfNum = 0;
00229 G4double YValue=0;
00230 G4Point3D Pt;
00231
00232 G4int a, b;
00233 for(a=0; a<nb_of_surfaces; a++)
00234 {
00235
00236
00237 srf = SurfaceVec[a];
00238 G4int Points = srf->GetNumberOfPoints();
00239
00240 for(b =0; b<Points; b++)
00241 {
00242 Pt = (G4Point3D)srf->GetPoint(b);
00243 if(YValue < Pt.y())
00244 {
00245 YValue = Pt.y();
00246 SrfNum = a;
00247 }
00248 }
00249 }
00250
00251
00252
00253 srf = SurfaceVec[SrfNum];
00254
00255
00256
00257
00258
00259 G4Point3D Pt1;
00260 G4Point3D Pt2;
00261 G4Point3D Pt3;
00262 G4Point3D Pt4;
00263
00264 G4Vector3D N1;
00265 G4Vector3D N2;
00266 G4Vector3D N3;
00267 G4Vector3D N4;
00268
00269 G4int* ConnectedList = new G4int[nb_of_surfaces];
00270
00271 for(a=0; a<nb_of_surfaces; a++)
00272 ConnectedList[a]=0;
00273
00274 G4Surface* ConnectedSrf;
00275
00276 for(a=0; a<nb_of_surfaces-1; a++)
00277 {
00278 if(ConnectedList[a] == 0)
00279 break;
00280 else
00281 ConnectedList[a]=1;
00282
00283 srf = SurfaceVec[a];
00284 G4int SrfPoints = srf->GetNumberOfPoints();
00285 N1 = (srf->Norm())->GetDir();
00286
00287 for(b=a+1; b<nb_of_surfaces; b++)
00288 {
00289 if(ConnectedList[b] == 1)
00290 break;
00291 else
00292 ConnectedList[b]=1;
00293
00294
00295
00296 ConnectedSrf = SurfaceVec[b];
00297
00298
00299
00300 G4int ConnSrfPoints = ConnectedSrf->GetNumberOfPoints();
00301
00302 for(G4int c=0;c<SrfPoints;c++)
00303 {
00304 Pt1 = srf->GetPoint(c);
00305
00306 for(G4int d=0;d<ConnSrfPoints;d++)
00307 {
00308
00309
00310 Pt2 = (ConnectedSrf)->GetPoint(d);
00311
00312 if( Pt1 == Pt2 )
00313 {
00314
00315
00316 N2 = ((ConnectedSrf)->Norm())->GetDir();
00317
00318
00319
00320 G4Vector3D CP1 = G4Vector3D( N1.cross(N2) );
00321 G4double CrossProd1 = CP1.x()+CP1.y()+CP1.z();
00322
00323
00324
00325 if(c==0)
00326 Pt3 = srf->GetPoint(c+1);
00327 else
00328 Pt3 = srf->GetPoint(0);
00329
00330 N3 = (Pt1-Pt3);
00331
00332 if(d==0)
00333 Pt4 = (ConnectedSrf)->GetPoint(d+1);
00334 else
00335 Pt4 = (ConnectedSrf)->GetPoint(0);
00336
00337 N4 = (Pt1-Pt4);
00338
00339 G4Vector3D CP2 = G4Vector3D( N3.cross(N4) );
00340 G4double CrossProd2 = CP2.x()+CP2.y()+CP2.z();
00341
00342 G4cout << "\nCroosProd2: " << CrossProd2;
00343
00344 if( (CrossProd1 < 0 && CrossProd2 < 0) ||
00345 (CrossProd1 > 0 && CrossProd2 > 0) )
00346 {
00347
00348
00349 (ConnectedSrf)->Norm()
00350 ->SetDir(-1 * (ConnectedSrf)->Norm()->GetDir());
00351
00352
00353
00354 CP1 = N1.cross(N2);
00355 CrossProd1 = CP1.x()+CP1.y()+CP1.z();
00356 }
00357 if(CrossProd1 > 0)
00358 Convex=0;
00359 }
00360 }
00361 }
00362 }
00363 }
00364
00365 delete []ConnectedList;
00366 }
00367
00368 G4int G4BREPSolid::IsBox()
00369 {
00370
00371
00372
00373
00374
00375
00376
00377 Box=0;
00378 G4Surface* srf1, *srf2;
00379 register G4int a;
00380
00381
00382
00383 for(a=0; a < nb_of_surfaces;a++)
00384 {
00385 srf1 = SurfaceVec[a];
00386
00387 if(srf1->MyType()==1)
00388 (srf1)->Project();
00389 else
00390 {
00391 PlaneSolid=0;
00392 return 0;
00393 }
00394 }
00395
00396
00397
00398 for(a=0; a < nb_of_surfaces;a++)
00399 {
00400 srf1 = SurfaceVec[a];
00401
00402 if (srf1->MyType()!=1)
00403 return 0;
00404 }
00405
00406 PlaneSolid = 1;
00407
00408
00409
00410 if(nb_of_surfaces!=6) return 0;
00411
00412 G4Point3D Pt;
00413 G4int Points;
00414 G4int Sides=0;
00415 G4int Opposite=0;
00416
00417 srf1 = SurfaceVec[0];
00418 Points = (srf1)->GetNumberOfPoints();
00419
00420 if(Points!=4)
00421 return 0;
00422
00423 G4Vector3D Normal1 = (srf1->Norm())->GetDir();
00424 G4double Result;
00425
00426 for(G4int b=1; b < nb_of_surfaces;b++)
00427 {
00428 srf2 = SurfaceVec[b];
00429 G4Vector3D Normal2 = ((srf2)->Norm())->GetDir();
00430 Result = std::fabs(Normal1 * Normal2);
00431
00432 if((Result != 0) && (Result != 1))
00433 return 0;
00434 else
00435 {
00436 if(!(G4int)Result)
00437 Sides++;
00438 else
00439 if(((G4int)Result) == 1)
00440 Opposite++;
00441 }
00442 }
00443
00444 if((Opposite != 1) && (Sides != nb_of_surfaces-2))
00445 return 0;
00446
00447 G4Vector3D x_axis(1,0,0);
00448 G4Vector3D y_axis(0,1,0);
00449
00450 if(((std::fabs(x_axis * Normal1) == 1) && (std::fabs(y_axis * Normal1) == 0)) ||
00451 ((std::fabs(x_axis * Normal1) == 0) && (std::fabs(y_axis * Normal1) == 1)) ||
00452 ((std::fabs(x_axis * Normal1) == 0) && (std::fabs(y_axis * Normal1) == 0)))
00453 AxisBox=1;
00454 else
00455 Box=1;
00456
00457 return 1;
00458 }
00459
00460 G4bool G4BREPSolid::IsConvex()
00461 {
00462 if (!PlaneSolid) { return 0; }
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472 G4Surface* Srf=0;
00473 G4Surface* ConnectedSrf=0;
00474 G4int Result;
00475 Convex = 1;
00476
00477 G4int a, b, c, d;
00478 for(a=0;a<nb_of_surfaces;a++)
00479 {
00480 Srf = SurfaceVec[a];
00481
00482
00483
00484
00485
00486 Result = Srf->IsConvex();
00487
00488 if (Result != -1)
00489 {
00490 Convex = 0;
00491 return 0;
00492 }
00493 }
00494
00495 Srf = SurfaceVec[0];
00496
00497 G4int ConnectingPoints=0;
00498 G4Point3D Pt1, Pt2;
00499 G4Vector3D N1, N2;
00500
00501
00502
00503
00504
00505
00506 const G4int maxCNum = (nb_of_surfaces-1)*nb_of_surfaces;
00507 G4int* ConnectedList = new G4int[maxCNum];
00508
00509 for(a=0; a<maxCNum; a++)
00510 {
00511 ConnectedList[a]=0;
00512 }
00513
00514 G4int Connections=0;
00515
00516 for(a=0; a<nb_of_surfaces-1; a++)
00517 {
00518 Srf = SurfaceVec[a];
00519 G4int SrfPoints = Srf->GetNumberOfPoints();
00520 Result=0;
00521
00522 for(b=0; b<nb_of_surfaces; b++)
00523 {
00524 if (b==a) { continue; }
00525
00526
00527
00528 ConnectedSrf = SurfaceVec[b];
00529
00530
00531
00532 G4int ConnSrfPoints = ConnectedSrf->GetNumberOfPoints();
00533
00534 for(c=0; c<SrfPoints; c++)
00535 {
00536 const G4Point3D& Pts1 =Srf->GetPoint(c);
00537
00538 for(d=0; d<ConnSrfPoints; d++)
00539 {
00540
00541
00542 const G4Point3D& Pts2 = ConnectedSrf->GetPoint(d);
00543 if (Pts1 == Pts2) { ConnectingPoints++; }
00544 }
00545 if (ConnectingPoints > 0) { break; }
00546 }
00547
00548 if (ConnectingPoints > 0)
00549 {
00550 Connections++;
00551 ConnectedList[Connections]=b;
00552 }
00553 ConnectingPoints=0;
00554 }
00555 }
00556
00557
00558
00559
00560 for(c=0; c<Connections; c++)
00561 {
00562 G4int Left=0;
00563 G4int Right =0;
00564 G4int tmp = ConnectedList[c];
00565
00566 Srf = SurfaceVec[tmp];
00567 ConnectedSrf = SurfaceVec[tmp+1];
00568
00569
00570
00571 N1 = Srf->Norm()->GetDir();
00572 N2 = ConnectedSrf->Norm()->GetDir();
00573
00574
00575
00576 G4Vector3D CP = G4Vector3D( N1.cross(N2) );
00577 G4double CrossProd = CP.x()+CP.y()+CP.z();
00578 if (CrossProd > 0) { Left++; }
00579 if (CrossProd < 0) { Right++; }
00580 if (Left&&Right)
00581 {
00582 Convex = 0;
00583 delete [] ConnectedList;
00584 return 0;
00585 }
00586 Connections=0;
00587 }
00588
00589 Convex=1;
00590
00591 delete [] ConnectedList;
00592
00593 return 1;
00594 }
00595
00596 G4bool G4BREPSolid::CalculateExtent(const EAxis pAxis,
00597 const G4VoxelLimits& pVoxelLimit,
00598 const G4AffineTransform& pTransform,
00599 G4double& pMin, G4double& pMax) const
00600 {
00601 G4Point3D Min = bbox->GetBoxMin();
00602 G4Point3D Max = bbox->GetBoxMax();
00603
00604 if (!pTransform.IsRotated())
00605 {
00606
00607
00608
00609
00610 G4double xoffset,xMin,xMax;
00611 G4double yoffset,yMin,yMax;
00612 G4double zoffset,zMin,zMax;
00613
00614 xoffset=pTransform.NetTranslation().x();
00615 xMin=xoffset+Min.x();
00616 xMax=xoffset+Max.x();
00617 if (pVoxelLimit.IsXLimited())
00618 {
00619 if (xMin>pVoxelLimit.GetMaxXExtent()
00620 ||xMax<pVoxelLimit.GetMinXExtent())
00621 {
00622 return false;
00623 }
00624 else
00625 {
00626 if (xMin<pVoxelLimit.GetMinXExtent())
00627 {
00628 xMin=pVoxelLimit.GetMinXExtent();
00629 }
00630 if (xMax>pVoxelLimit.GetMaxXExtent())
00631 {
00632 xMax=pVoxelLimit.GetMaxXExtent();
00633 }
00634 }
00635 }
00636
00637 yoffset=pTransform.NetTranslation().y();
00638 yMin=yoffset+Min.y();
00639 yMax=yoffset+Max.y();
00640 if (pVoxelLimit.IsYLimited())
00641 {
00642 if (yMin>pVoxelLimit.GetMaxYExtent()
00643 ||yMax<pVoxelLimit.GetMinYExtent())
00644 {
00645 return false;
00646 }
00647 else
00648 {
00649 if (yMin<pVoxelLimit.GetMinYExtent())
00650 {
00651 yMin=pVoxelLimit.GetMinYExtent();
00652 }
00653 if (yMax>pVoxelLimit.GetMaxYExtent())
00654 {
00655 yMax=pVoxelLimit.GetMaxYExtent();
00656 }
00657 }
00658 }
00659
00660 zoffset=pTransform.NetTranslation().z();
00661 zMin=zoffset+Min.z();
00662 zMax=zoffset+Max.z();
00663 if (pVoxelLimit.IsZLimited())
00664 {
00665 if (zMin>pVoxelLimit.GetMaxZExtent()
00666 ||zMax<pVoxelLimit.GetMinZExtent())
00667 {
00668 return false;
00669 }
00670 else
00671 {
00672 if (zMin<pVoxelLimit.GetMinZExtent())
00673 {
00674 zMin=pVoxelLimit.GetMinZExtent();
00675 }
00676 if (zMax>pVoxelLimit.GetMaxZExtent())
00677 {
00678 zMax=pVoxelLimit.GetMaxZExtent();
00679 }
00680 }
00681 }
00682
00683 switch (pAxis)
00684 {
00685 case kXAxis:
00686 pMin=xMin;
00687 pMax=xMax;
00688 break;
00689 case kYAxis:
00690 pMin=yMin;
00691 pMax=yMax;
00692 break;
00693 case kZAxis:
00694 pMin=zMin;
00695 pMax=zMax;
00696 break;
00697 default:
00698 break;
00699 }
00700 pMin-=kCarTolerance;
00701 pMax+=kCarTolerance;
00702
00703 return true;
00704 }
00705 else
00706 {
00707
00708
00709 G4bool existsAfterClip=false;
00710 G4ThreeVectorList *vertices;
00711
00712 pMin=+kInfinity;
00713 pMax=-kInfinity;
00714
00715
00716
00717 vertices=CreateRotatedVertices(pTransform);
00718 ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
00719 ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax);
00720 ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
00721
00722 if ( (pMin!=kInfinity) || (pMax!=-kInfinity) )
00723 {
00724 existsAfterClip=true;
00725
00726
00727
00728 pMin-=kCarTolerance;
00729 pMax+=kCarTolerance;
00730 }
00731 else
00732 {
00733
00734
00735
00736
00737
00738 G4ThreeVector clipCentre(
00739 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
00740 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
00741 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
00742
00743 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
00744 {
00745 existsAfterClip=true;
00746 pMin=pVoxelLimit.GetMinExtent(pAxis);
00747 pMax=pVoxelLimit.GetMaxExtent(pAxis);
00748 }
00749 }
00750 delete vertices;
00751 return existsAfterClip;
00752 }
00753 }
00754
00755 G4ThreeVectorList*
00756 G4BREPSolid::CreateRotatedVertices(const G4AffineTransform& pTransform) const
00757 {
00758 G4Point3D Min = bbox->GetBoxMin();
00759 G4Point3D Max = bbox->GetBoxMax();
00760
00761 G4ThreeVectorList *vertices;
00762 vertices=new G4ThreeVectorList();
00763
00764 if (vertices)
00765 {
00766 vertices->reserve(8);
00767 G4ThreeVector vertex0(Min.x(),Min.y(),Min.z());
00768 G4ThreeVector vertex1(Max.x(),Min.y(),Min.z());
00769 G4ThreeVector vertex2(Max.x(),Max.y(),Min.z());
00770 G4ThreeVector vertex3(Min.x(),Max.y(),Min.z());
00771 G4ThreeVector vertex4(Min.x(),Min.y(),Max.z());
00772 G4ThreeVector vertex5(Max.x(),Min.y(),Max.z());
00773 G4ThreeVector vertex6(Max.x(),Max.y(),Max.z());
00774 G4ThreeVector vertex7(Min.x(),Max.y(),Max.z());
00775
00776 vertices->push_back(pTransform.TransformPoint(vertex0));
00777 vertices->push_back(pTransform.TransformPoint(vertex1));
00778 vertices->push_back(pTransform.TransformPoint(vertex2));
00779 vertices->push_back(pTransform.TransformPoint(vertex3));
00780 vertices->push_back(pTransform.TransformPoint(vertex4));
00781 vertices->push_back(pTransform.TransformPoint(vertex5));
00782 vertices->push_back(pTransform.TransformPoint(vertex6));
00783 vertices->push_back(pTransform.TransformPoint(vertex7));
00784 }
00785 else
00786 {
00787 G4Exception("G4BREPSolid::CreateRotatedVertices()", "GeomSolids0003",
00788 FatalException, "Out of memory - Cannot allocate vertices!");
00789 }
00790 return vertices;
00791 }
00792
00793 EInside G4BREPSolid::Inside(register const G4ThreeVector& Pt) const
00794 {
00795
00796
00797
00798 const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;
00799
00800 G4Vector3D v(1, 0, 0.01);
00801 G4Vector3D Pttmp(Pt);
00802 G4Vector3D Vtmp(v);
00803 G4Ray r(Pttmp, Vtmp);
00804
00805
00806
00807 if( !GetBBox()->Inside(Pttmp) )
00808 return kOutside;
00809
00810
00811
00812 Reset();
00813
00814
00815
00816
00817 TestSurfaceBBoxes(r);
00818
00819 G4int hits=0, samehit=0;
00820
00821 for(G4int a=0; a < nb_of_surfaces; a++)
00822 {
00823 if(SurfaceVec[a]->IsActive())
00824 {
00825
00826
00827
00828
00829
00830 if( (SurfaceVec[a]->Intersect(r)) & 1 )
00831 {
00832
00833
00834 if(SurfaceVec[a]->GetDistance() < sqrHalfTolerance)
00835 return kSurface;
00836
00837
00838
00839 for(G4int i=0; i<a; i++)
00840 if(SurfaceVec[a]->GetDistance() == SurfaceVec[i]->GetDistance())
00841 {
00842 samehit++;
00843 break;
00844 }
00845
00846
00847
00848 if(!samehit)
00849 hits++;
00850 }
00851 }
00852 }
00853
00854
00855
00856
00857 if(hits&1)
00858 return kInside;
00859 else
00860 return kOutside;
00861 }
00862
00863 G4ThreeVector G4BREPSolid::SurfaceNormal(const G4ThreeVector& Pt) const
00864 {
00865
00866
00867
00868
00869 const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;
00870 G4int iplane;
00871
00872
00873
00874 for(iplane = 0; iplane < nb_of_surfaces; iplane++)
00875 {
00876 if(SurfaceVec[iplane]->HowNear(Pt) < sqrHalfTolerance)
00877
00878 break;
00879 }
00880
00881
00882
00883 G4ThreeVector norm = SurfaceVec[iplane]->SurfaceNormal(Pt);
00884
00885 return norm.unit();
00886 }
00887
00888 G4double G4BREPSolid::DistanceToIn(const G4ThreeVector& Pt) const
00889 {
00890
00891
00892
00893
00894 G4double *dists = new G4double[nb_of_surfaces];
00895 G4int a;
00896
00897
00898
00899 Reset();
00900
00901
00902
00903
00904 for(a=0; a< nb_of_surfaces; a++)
00905 dists[a] = SurfaceVec[a]->HowNear(Pt);
00906
00907 G4double Dist = kInfinity;
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917 for(a = 0; a < nb_of_surfaces; a++)
00918 if( std::fabs(Dist) > std::fabs(dists[a]) )
00919
00920 Dist = dists[a];
00921
00922 delete[] dists;
00923
00924 if(Dist == kInfinity)
00925 return 0;
00926 else
00927 return std::fabs(Dist);
00928 }
00929
00930 G4double G4BREPSolid::DistanceToIn(register const G4ThreeVector& Pt,
00931 register const G4ThreeVector& V ) const
00932 {
00933
00934
00935
00936
00937
00938
00939 G4int a;
00940
00941
00942
00943 Reset();
00944
00945 const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;
00946 G4Vector3D Pttmp(Pt);
00947 G4Vector3D Vtmp(V);
00948 G4Ray r(Pttmp, Vtmp);
00949
00950
00951
00952
00953 TestSurfaceBBoxes(r);
00954
00955 ShortestDistance = kInfinity;
00956
00957 for(a=0; a< nb_of_surfaces; a++)
00958 {
00959 if( SurfaceVec[a]->IsActive() )
00960 {
00961
00962
00963 if( SurfaceVec[a]->Intersect(r) )
00964 {
00965 G4double surfDistance = SurfaceVec[a]->GetDistance();
00966
00967
00968
00969 if( surfDistance < ShortestDistance )
00970 {
00971 if( surfDistance > sqrHalfTolerance )
00972 {
00973 ShortestDistance = surfDistance;
00974 }
00975 else
00976 {
00977
00978
00979
00980 G4Vector3D Norm = SurfaceVec[a]->SurfaceNormal(Pttmp);
00981
00982 if( (Norm * Vtmp) < 0 )
00983 {
00984 ShortestDistance = surfDistance;
00985 }
00986 }
00987 }
00988 }
00989 }
00990 }
00991
00992
00993
00994
00995 if(ShortestDistance != kInfinity)
00996 return std::sqrt(ShortestDistance);
00997 else
00998 return kInfinity;
00999 }
01000
01001 G4double G4BREPSolid::DistanceToOut(register const G4ThreeVector& P,
01002 register const G4ThreeVector& D,
01003 const G4bool,
01004 G4bool *validNorm,
01005 G4ThreeVector* ) const
01006 {
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017 Reset();
01018
01019 const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;
01020 G4Vector3D Ptv = P;
01021 G4int a;
01022
01023 if(validNorm)
01024 *validNorm=false;
01025
01026 G4Vector3D Pttmp(Ptv);
01027 G4Vector3D Vtmp(D);
01028
01029 G4Ray r(Pttmp, Vtmp);
01030
01031
01032
01033
01034 TestSurfaceBBoxes(r);
01035
01036 ShortestDistance = kInfinity;
01037
01038 for(a=0; a< nb_of_surfaces; a++)
01039 {
01040 if(SurfaceVec[a]->IsActive())
01041 {
01042
01043
01044 if( (SurfaceVec[a]->Intersect(r)) )
01045 {
01046
01047
01048 G4double surfDistance = SurfaceVec[a]->GetDistance();
01049 if( surfDistance < ShortestDistance )
01050 {
01051 if( surfDistance > sqrHalfTolerance )
01052 {
01053 ShortestDistance = surfDistance;
01054 }
01055 else
01056 {
01057
01058 }
01059 }
01060 }
01061 }
01062 }
01063
01064
01065
01066
01067 if(ShortestDistance != kInfinity)
01068 return std::sqrt(ShortestDistance);
01069 else
01070 return 0.0;
01071 }
01072
01073 G4double G4BREPSolid::DistanceToOut(const G4ThreeVector& Pt)const
01074 {
01075
01076
01077
01078
01079 G4double *dists = new G4double[nb_of_surfaces];
01080 G4int a;
01081
01082
01083
01084 Reset();
01085
01086
01087
01088
01089 for(a=0; a< nb_of_surfaces; a++)
01090 dists[a] = SurfaceVec[a]->HowNear(Pt);
01091
01092 G4double Dist = kInfinity;
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102 for(a = 0; a < nb_of_surfaces; a++)
01103 if( std::fabs(Dist) > std::fabs(dists[a]) )
01104
01105 Dist = dists[a];
01106
01107 delete[] dists;
01108
01109 if(Dist == kInfinity)
01110 return 0;
01111 else
01112 return std::fabs(Dist);
01113 }
01114
01115 void G4BREPSolid::DescribeYourselfTo (G4VGraphicsScene& scene) const
01116 {
01117 scene.AddSolid (*this);
01118 }
01119
01120 G4Polyhedron* G4BREPSolid::CreatePolyhedron () const
01121 {
01122
01123
01124 G4Point3D Min = bbox->GetBoxMin();
01125 G4Point3D Max = bbox->GetBoxMax();
01126
01127 return new G4PolyhedronBox (Max.x(), Max.y(), Max.z());
01128 }
01129
01130 G4NURBS* G4BREPSolid::CreateNURBS () const
01131 {
01132
01133
01134 G4Point3D Min = bbox->GetBoxMin();
01135 G4Point3D Max = bbox->GetBoxMax();
01136
01137 return new G4NURBSbox (Max.x(), Max.y(), Max.z());
01138 }
01139
01140 void G4BREPSolid::CalcBBoxes()
01141 {
01142
01143
01144
01145 G4Surface* srf;
01146 G4Point3D min, max;
01147
01148 if(active)
01149 {
01150 min = PINFINITY;
01151 max = -PINFINITY;
01152
01153 for(G4int a = 0;a < nb_of_surfaces;a++)
01154 {
01155
01156
01157 srf = SurfaceVec[a];
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172 {
01173 srf->CalcBBox();
01174 G4Point3D box_min = srf->GetBBox()->GetBoxMin();
01175 G4Point3D box_max = srf->GetBBox()->GetBoxMax();
01176
01177
01178
01179
01180
01181
01182 if(max.x() < box_max.x()) max.setX(box_max.x());
01183 if(max.y() < box_max.y()) max.setY(box_max.y());
01184 if(max.z() < box_max.z()) max.setZ(box_max.z());
01185
01186
01187
01188 if(min.x() > box_min.x()) min.setX(box_min.x());
01189 if(min.y() > box_min.y()) min.setY(box_min.y());
01190 if(min.z() > box_min.z()) min.setZ(box_min.z());
01191 }
01192 }
01193 bbox = new G4BoundingBox3D(min, max);
01194 return;
01195 }
01196 G4Exception("G4BREPSolid::CalcBBoxes()", "GeomSolids1002",
01197 JustWarning, "No bbox calculated for solid.");
01198 }
01199
01200 void G4BREPSolid::RemoveHiddenFaces(register const G4Ray& rayref,
01201 G4int In ) const
01202 {
01203
01204
01205
01206
01207
01208
01209 register G4Surface* srf;
01210 register const G4Vector3D& RayDir = rayref.GetDir();
01211 register G4double Result;
01212 G4int a;
01213
01214
01215
01216 if(!In)
01217 for(a=0; a<nb_of_surfaces; a++)
01218 {
01219
01220
01221 srf = SurfaceVec[a];
01222 if(srf->MyType()==1)
01223 {
01224 const G4Vector3D& Normal = (srf->Norm())->GetDir();
01225 Result = (RayDir * Normal);
01226
01227 if( Result >= 0 )
01228 srf->Deactivate();
01229 }
01230 }
01231 else
01232 for(a=0; a<nb_of_surfaces; a++)
01233 {
01234
01235
01236
01237
01238
01239 srf = SurfaceVec[a];
01240 if(srf->MyType()==1)
01241 {
01242 const G4Vector3D& Normal = (srf->Norm())->GetDir();
01243 Result = (RayDir * Normal);
01244
01245 if( Result < 0 )
01246 srf->Deactivate();
01247 }
01248 }
01249 }
01250
01251 void G4BREPSolid::TestSurfaceBBoxes(register const G4Ray& rayref) const
01252 {
01253 register G4Surface* srf;
01254 G4int active_srfs = nb_of_surfaces;
01255
01256
01257
01258
01259 G4int intersection=0;
01260
01261 for(G4int a=0;a<nb_of_surfaces;a++)
01262 {
01263
01264
01265 srf = SurfaceVec[a];
01266
01267 if(srf->IsActive())
01268 {
01269
01270
01271 if(srf->MyType() != 1)
01272 {
01273 if(srf->GetBBox()->Test(rayref))
01274 srf->SetDistance(bbox->GetDistance());
01275 else
01276 {
01277
01278
01279 srf->Deactivate();
01280 active_srfs--;
01281 }
01282 }
01283 else
01284 {
01285
01286 intersection = srf->Intersect(rayref);
01287
01288 if(!intersection)
01289 active_srfs--;
01290 }
01291 }
01292 else
01293 active_srfs--;
01294 }
01295
01296 if(!active_srfs) Active(0);
01297 }
01298
01299
01300 G4int G4BREPSolid::Intersect(register const G4Ray& rayref) const
01301 {
01302
01303
01304
01305 const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;
01306
01307 register G4Surface* srf;
01308 G4double HitDistance = -1;
01309 const G4Point3D& RayStart = rayref.GetStart();
01310 const G4Point3D& RayDir = rayref.GetDir();
01311
01312 G4int result=1;
01313
01314
01315
01316
01317 QuickSort(SurfaceVec, 0, nb_of_surfaces-1);
01318 G4int Number=0;
01319
01320
01321
01322 for(register G4int a=0;a<nb_of_surfaces;a++)
01323 {
01324 srf = SurfaceVec[a];
01325
01326 if(srf->IsActive())
01327 {
01328 result = srf->Intersect(rayref);
01329 if(result)
01330 {
01331
01332
01333 const G4Point3D& closest_point = srf->GetClosestHit();
01334
01335
01336
01337
01338 if( !( (srf->GetDistance() < sqrHalfTolerance) ||
01339 (RayDir.dot(srf->SurfaceNormal(closest_point)) > 0) ) )
01340 {
01341
01342 if(srf->MyType()==1)
01343 HitDistance = srf->GetDistance();
01344 else
01345 {
01346
01347
01348
01349 HitDistance = RayStart.distance2(closest_point);
01350 }
01351 }
01352 }
01353 else
01354 {
01355 srf->Deactivate();
01356 }
01357 }
01358 Number++;
01359 }
01360
01361 if(HitDistance < 0)
01362 return 0;
01363
01364 QuickSort(SurfaceVec, 0, nb_of_surfaces-1);
01365
01366 if(!(SurfaceVec[0]->IsActive()))
01367 return 0;
01368
01369 ((G4BREPSolid*)this)->intersection_point = SurfaceVec[0]->GetClosestHit();
01370 bbox->SetDistance(HitDistance);
01371
01372 return 1;
01373 }
01374
01375 G4int G4BREPSolid::FinalEvaluation(register const G4Ray& rayref,
01376 G4int ToIn ) const
01377 {
01378 const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;
01379 register G4Surface* srf;
01380 G4double Dist=0;
01381
01382 ((G4BREPSolid*)this)->intersectionDistance = kInfinity;
01383
01384 for(register G4int a=0;a<nb_of_surfaces;a++)
01385 {
01386 srf = SurfaceVec[a];
01387
01388 if(srf->IsActive())
01389 {
01390 const G4Point3D& srf_intersection = srf->Evaluation(rayref);
01391
01392
01393
01394 if(srf->MyType() != 1)
01395 {
01396 G4Point3D start = rayref.GetStart();
01397 Dist = srf_intersection.distance2(start);
01398 }
01399 else
01400 Dist = srf->GetDistance();
01401
01402
01403
01404
01405 if(Dist < sqrHalfTolerance)
01406 {
01407 if(ToIn)
01408 {
01409 const G4Vector3D& Dir = rayref.GetDir();
01410 const G4Point3D& Hit = srf->GetClosestHit();
01411 const G4Vector3D& Norm = srf->SurfaceNormal(Hit);
01412
01413 if(( Dir * Norm ) >= 0)
01414 {
01415 Dist = kInfinity;
01416 srf->Deactivate();
01417 }
01418
01419
01420 }
01421 else
01422 {
01423 Dist = kInfinity;
01424 srf->Deactivate();
01425 }
01426 }
01427
01428
01429
01430
01431 if(Dist < intersectionDistance)
01432 {
01433
01434
01435
01436 const G4Point3D& Pt = rayref.GetStart();
01437 const G4Vector3D& Dir = rayref.GetDir();
01438
01439 G4Point3D TestPoint = (0.00001*Dir) + Pt;
01440 G4double TestDistance = srf_intersection.distance2(TestPoint);
01441
01442 if(TestDistance > Dist)
01443 {
01444
01445
01446 Dist = kInfinity;
01447 srf->Deactivate();
01448 }
01449 else
01450 {
01451 ((G4BREPSolid*)this)->intersectionDistance = Dist;
01452 ((G4BREPSolid*)this)->intersection_point = srf_intersection;
01453 }
01454
01455
01456
01457
01458 if(srf->IsActive())
01459 {
01460 if(a+1<nb_of_surfaces)
01461 {
01462 const G4Point3D& Hit = srf->GetClosestHit();
01463 const G4Vector3D& Norm = srf->SurfaceNormal(Hit);
01464
01465
01466
01467 if(( Dir * Norm ) < 0)
01468 {
01469 Dist = kInfinity;
01470 srf->Deactivate();
01471 }
01472
01473
01474
01475 ShortestDistance = Dist;
01476 }
01477 else
01478 {
01479 ShortestDistance = Dist;
01480 return 1;
01481 }
01482 }
01483 }
01484 }
01485 else
01486 {
01487
01488
01489
01490 }
01491 }
01492 if(intersectionDistance < kInfinity)
01493 return 1;
01494
01495 return 0;
01496 }
01497
01498 G4Point3D G4BREPSolid::Scope() const
01499 {
01500 G4Point3D scope;
01501 G4Point3D Max = bbox->GetBoxMax();
01502 G4Point3D Min = bbox->GetBoxMin();
01503
01504 scope.setX(std::fabs(Max.x()) - std::fabs(Min.x()));
01505 scope.setY(std::fabs(Max.y()) - std::fabs(Min.y()));
01506 scope.setZ(std::fabs(Max.z()) - std::fabs(Min.z()));
01507
01508 return scope;
01509 }
01510
01511 std::ostream& G4BREPSolid::StreamInfo(std::ostream& os) const
01512 {
01513 os << "-----------------------------------------------------------\n"
01514 << " *** Dump for solid - " << GetName() << " ***\n"
01515 << " ===================================================\n"
01516 << " Solid type: " << GetEntityType() << "\n"
01517 << " Parameters: \n"
01518 << " Number of solids: " << NumberOfSolids << "\n"
01519 << "-----------------------------------------------------------\n";
01520
01521 return os;
01522 }
01523
01524 G4Polyhedron* G4BREPSolid::GetPolyhedron () const
01525 {
01526 if (!fpPolyhedron ||
01527 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
01528 fpPolyhedron->GetNumberOfRotationSteps())
01529 {
01530 delete fpPolyhedron;
01531 fpPolyhedron = CreatePolyhedron();
01532 }
01533 return fpPolyhedron;
01534 }