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 "G4TwistTubsFlatSide.hh"
00045 #include "G4GeometryTolerance.hh"
00046
00047
00048
00049
00050 G4TwistTubsFlatSide::G4TwistTubsFlatSide(const G4String &name,
00051 const G4RotationMatrix &rot,
00052 const G4ThreeVector &tlate,
00053 const G4ThreeVector &n,
00054 const EAxis axis0 ,
00055 const EAxis axis1 ,
00056 G4double axis0min,
00057 G4double axis1min,
00058 G4double axis0max,
00059 G4double axis1max )
00060 : G4VTwistSurface(name, rot, tlate, 0, axis0, axis1,
00061 axis0min, axis1min, axis0max, axis1max)
00062 {
00063 if (axis0 == kPhi && axis1 == kRho) {
00064 G4Exception("G4TwistTubsFlatSide::G4TwistTubsFlatSide()",
00065 "GeomSolids0002", FatalErrorInArgument,
00066 "Should swap axis0 and axis1!");
00067 }
00068
00069 G4ThreeVector normal = rot.inverse()*n;
00070 fCurrentNormal.normal = normal.unit();
00071 fIsValidNorm = true;
00072
00073 SetCorners();
00074 SetBoundaries();
00075
00076 fSurfaceArea = 1 ;
00077
00078 }
00079
00080
00081
00082 G4TwistTubsFlatSide::G4TwistTubsFlatSide( const G4String &name,
00083 G4double EndInnerRadius[2],
00084 G4double EndOuterRadius[2],
00085 G4double DPhi,
00086 G4double EndPhi[2],
00087 G4double EndZ[2],
00088 G4int handedness )
00089 : G4VTwistSurface(name)
00090 {
00091 fHandedness = handedness;
00092 fAxis[0] = kRho;
00093 fAxis[1] = kPhi;
00094 G4int i = (handedness < 0 ? 0 : 1);
00095 fAxisMin[0] = EndInnerRadius[i];
00096 fAxisMax[0] = EndOuterRadius[i];
00097 fAxisMin[1] = -0.5*DPhi;
00098 fAxisMax[1] = -fAxisMin[1];
00099 fCurrentNormal.normal.set(0, 0, (fHandedness < 0 ? -1 : 1));
00100
00101 fRot.rotateZ(EndPhi[i]);
00102 fTrans.set(0, 0, EndZ[i]);
00103 fIsValidNorm = true;
00104
00105 SetCorners();
00106 SetBoundaries();
00107
00108 fSurfaceArea = 0.5*DPhi * (EndOuterRadius[i]*EndOuterRadius[i]
00109 - EndInnerRadius[i]*EndInnerRadius[i] ) ;
00110
00111 }
00112
00113
00114
00115
00116
00117 G4TwistTubsFlatSide::G4TwistTubsFlatSide( __void__& a )
00118 : G4VTwistSurface(a), fSurfaceArea(0.)
00119 {
00120 }
00121
00122
00123
00124
00125
00126 G4TwistTubsFlatSide::~G4TwistTubsFlatSide()
00127 {
00128 }
00129
00130
00131
00132
00133 G4ThreeVector G4TwistTubsFlatSide::GetNormal(const G4ThreeVector & ,
00134 G4bool isGlobal)
00135 {
00136 if (isGlobal) {
00137 return ComputeGlobalDirection(fCurrentNormal.normal);
00138 } else {
00139 return fCurrentNormal.normal;
00140 }
00141 }
00142
00143
00144
00145
00146 G4int G4TwistTubsFlatSide::DistanceToSurface(const G4ThreeVector &gp,
00147 const G4ThreeVector &gv,
00148 G4ThreeVector gxx[],
00149 G4double distance[],
00150 G4int areacode[],
00151 G4bool isvalid[],
00152 EValidate validate)
00153 {
00154 fCurStatWithV.ResetfDone(validate, &gp, &gv);
00155
00156 if (fCurStatWithV.IsDone()) {
00157 G4int i;
00158 for (i=0; i<fCurStatWithV.GetNXX(); i++) {
00159 gxx[i] = fCurStatWithV.GetXX(i);
00160 distance[i] = fCurStatWithV.GetDistance(i);
00161 areacode[i] = fCurStatWithV.GetAreacode(i);
00162 isvalid[i] = fCurStatWithV.IsValid(i);
00163 }
00164 return fCurStatWithV.GetNXX();
00165 } else {
00166
00167 G4int i;
00168 for (i=0; i<2; i++) {
00169 distance[i] = kInfinity;
00170 areacode[i] = sOutside;
00171 isvalid[i] = false;
00172 gxx[i].set(kInfinity, kInfinity, kInfinity);
00173 }
00174 }
00175
00176 G4ThreeVector p = ComputeLocalPoint(gp);
00177 G4ThreeVector v = ComputeLocalDirection(gv);
00178
00179
00180
00181
00182
00183
00184 if (std::fabs(p.z()) == 0.) {
00185 distance[0] = 0;
00186 G4ThreeVector xx = p;
00187 gxx[0] = ComputeGlobalPoint(xx);
00188
00189 if (validate == kValidateWithTol) {
00190 areacode[0] = GetAreaCode(xx);
00191 if (!IsOutside(areacode[0])) {
00192 isvalid[0] = true;
00193 }
00194 } else if (validate == kValidateWithoutTol) {
00195 areacode[0] = GetAreaCode(xx, false);
00196 if (IsInside(areacode[0])) {
00197 isvalid[0] = true;
00198 }
00199 } else {
00200 areacode[0] = sInside;
00201 isvalid[0] = true;
00202 }
00203
00204 return 1;
00205 }
00206
00207
00208
00209
00210 if (v.z() == 0) {
00211
00212 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
00213 isvalid[0], 0, validate, &gp, &gv);
00214 return 0;
00215 }
00216
00217 distance[0] = - (p.z() / v.z());
00218
00219 G4ThreeVector xx = p + distance[0]*v;
00220 gxx[0] = ComputeGlobalPoint(xx);
00221
00222 if (validate == kValidateWithTol) {
00223 areacode[0] = GetAreaCode(xx);
00224 if (!IsOutside(areacode[0])) {
00225 if (distance[0] >= 0) isvalid[0] = true;
00226 }
00227 } else if (validate == kValidateWithoutTol) {
00228 areacode[0] = GetAreaCode(xx, false);
00229 if (IsInside(areacode[0])) {
00230 if (distance[0] >= 0) isvalid[0] = true;
00231 }
00232 } else {
00233 areacode[0] = sInside;
00234 if (distance[0] >= 0) isvalid[0] = true;
00235 }
00236
00237 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
00238 isvalid[0], 1, validate, &gp, &gv);
00239
00240 #ifdef G4TWISTDEBUG
00241 G4cerr << "ERROR - G4TwistTubsFlatSide::DistanceToSurface(p,v)" << G4endl;
00242 G4cerr << " Name : " << GetName() << G4endl;
00243 G4cerr << " xx : " << xx << G4endl;
00244 G4cerr << " gxx[0] : " << gxx[0] << G4endl;
00245 G4cerr << " dist[0] : " << distance[0] << G4endl;
00246 G4cerr << " areacode[0] : " << areacode[0] << G4endl;
00247 G4cerr << " isvalid[0] : " << isvalid[0] << G4endl;
00248 #endif
00249 return 1;
00250 }
00251
00252
00253
00254
00255 G4int G4TwistTubsFlatSide::DistanceToSurface(const G4ThreeVector &gp,
00256 G4ThreeVector gxx[],
00257 G4double distance[],
00258 G4int areacode[])
00259 {
00260
00261
00262
00263
00264 fCurStat.ResetfDone(kDontValidate, &gp);
00265
00266 if (fCurStat.IsDone()) {
00267 G4int i;
00268 for (i=0; i<fCurStat.GetNXX(); i++) {
00269 gxx[i] = fCurStat.GetXX(i);
00270 distance[i] = fCurStat.GetDistance(i);
00271 areacode[i] = fCurStat.GetAreacode(i);
00272 }
00273 return fCurStat.GetNXX();
00274 } else {
00275
00276 G4int i;
00277 for (i=0; i<2; i++) {
00278 distance[i] = kInfinity;
00279 areacode[i] = sOutside;
00280 gxx[i].set(kInfinity, kInfinity, kInfinity);
00281 }
00282 }
00283
00284 G4ThreeVector p = ComputeLocalPoint(gp);
00285 G4ThreeVector xx;
00286
00287
00288
00289 if (std::fabs(p.z()) <= 0.5 * kCarTolerance) {
00290 distance[0] = 0;
00291 xx = p;
00292 } else {
00293 distance[0] = std::fabs(p.z());
00294 xx.set(p.x(), p.y(), 0);
00295 }
00296
00297 gxx[0] = ComputeGlobalPoint(xx);
00298 areacode[0] = sInside;
00299 G4bool isvalid = true;
00300 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
00301 isvalid, 1, kDontValidate, &gp);
00302 return 1;
00303
00304 }
00305
00306
00307
00308
00309 G4int G4TwistTubsFlatSide::GetAreaCode(const G4ThreeVector &xx,
00310 G4bool withTol)
00311 {
00312
00313 static const G4double rtol
00314 = 0.5*G4GeometryTolerance::GetInstance()->GetRadialTolerance();
00315
00316 G4int areacode = sInside;
00317
00318 if (fAxis[0] == kRho && fAxis[1] == kPhi) {
00319 G4int rhoaxis = 0;
00320
00321
00322 G4ThreeVector dphimin;
00323 G4ThreeVector dphimax;
00324 dphimin = GetCorner(sC0Max1Min);
00325 dphimax = GetCorner(sC0Max1Max);
00326
00327 if (withTol) {
00328
00329 G4bool isoutside = false;
00330
00331
00332
00333 if (xx.getRho() <= fAxisMin[rhoaxis] + rtol) {
00334
00335 areacode |= (sAxis0 & (sAxisRho | sAxisMin)) | sBoundary;
00336 if (xx.getRho() < fAxisMin[rhoaxis] - rtol) isoutside = true;
00337
00338 } else if (xx.getRho() >= fAxisMax[rhoaxis] - rtol) {
00339
00340 areacode |= (sAxis0 & (sAxisRho | sAxisMax)) | sBoundary;
00341 if (xx.getRho() > fAxisMax[rhoaxis] + rtol) isoutside = true;
00342
00343 }
00344
00345
00346
00347 if (AmIOnLeftSide(xx, dphimin) >= 0) {
00348
00349 areacode |= (sAxis1 & (sAxisPhi | sAxisMin));
00350 if (areacode & sBoundary) areacode |= sCorner;
00351 else areacode |= sBoundary;
00352
00353 if (AmIOnLeftSide(xx, dphimin) > 0) isoutside = true;
00354
00355 } else if (AmIOnLeftSide(xx, dphimax) <= 0) {
00356
00357 areacode |= (sAxis1 & (sAxisPhi | sAxisMax));
00358 if (areacode & sBoundary) areacode |= sCorner;
00359 else areacode |= sBoundary;
00360
00361 if (AmIOnLeftSide(xx, dphimax) < 0) isoutside = true;
00362
00363 }
00364
00365
00366
00367
00368 if (isoutside) {
00369 G4int tmpareacode = areacode & (~sInside);
00370 areacode = tmpareacode;
00371 } else if ((areacode & sBoundary) != sBoundary) {
00372 areacode |= (sAxis0 & sAxisRho) | (sAxis1 & sAxisPhi);
00373 }
00374
00375 } else {
00376
00377
00378
00379 if (xx.getRho() < fAxisMin[rhoaxis]) {
00380 areacode |= (sAxis0 & (sAxisRho | sAxisMin)) | sBoundary;
00381 } else if (xx.getRho() > fAxisMax[rhoaxis]) {
00382 areacode |= (sAxis0 & (sAxisRho | sAxisMax)) | sBoundary;
00383 }
00384
00385
00386
00387 if (AmIOnLeftSide(xx, dphimin, false) >= 0) {
00388 areacode |= (sAxis1 & (sAxisPhi | sAxisMin)) ;
00389 if (areacode & sBoundary) areacode |= sCorner;
00390 else areacode |= sBoundary;
00391
00392 } else if (AmIOnLeftSide(xx, dphimax, false) <= 0) {
00393 areacode |= (sAxis1 & (sAxisPhi | sAxisMax)) ;
00394 if (areacode & sBoundary) areacode |= sCorner;
00395 else areacode |= sBoundary;
00396
00397 }
00398
00399 if ((areacode & sBoundary) != sBoundary) {
00400 areacode |= (sAxis0 & sAxisRho) | (sAxis1 & sAxisPhi);
00401 }
00402
00403 }
00404 return areacode;
00405 } else {
00406
00407 std::ostringstream message;
00408 message << "Feature NOT implemented !" << G4endl
00409 << " fAxis[0] = " << fAxis[0] << G4endl
00410 << " fAxis[1] = " << fAxis[1];
00411 G4Exception("G4TwistTubsFlatSide::GetAreaCode()", "GeomSolids0001",
00412 FatalException, message);
00413 }
00414 return areacode;
00415 }
00416
00417
00418
00419
00420
00421 void G4TwistTubsFlatSide::SetCorners()
00422 {
00423
00424
00425 if (fAxis[0] == kRho && fAxis[1] == kPhi) {
00426
00427 G4int rhoaxis = 0;
00428 G4int phiaxis = 1;
00429
00430 G4double x, y, z;
00431
00432 x = fAxisMin[rhoaxis]*std::cos(fAxisMin[phiaxis]);
00433 y = fAxisMin[rhoaxis]*std::sin(fAxisMin[phiaxis]);
00434 z = 0;
00435 SetCorner(sC0Min1Min, x, y, z);
00436
00437 x = fAxisMax[rhoaxis]*std::cos(fAxisMin[phiaxis]);
00438 y = fAxisMax[rhoaxis]*std::sin(fAxisMin[phiaxis]);
00439 z = 0;
00440 SetCorner(sC0Max1Min, x, y, z);
00441
00442 x = fAxisMax[rhoaxis]*std::cos(fAxisMax[phiaxis]);
00443 y = fAxisMax[rhoaxis]*std::sin(fAxisMax[phiaxis]);
00444 z = 0;
00445 SetCorner(sC0Max1Max, x, y, z);
00446
00447 x = fAxisMin[rhoaxis]*std::cos(fAxisMax[phiaxis]);
00448 y = fAxisMin[rhoaxis]*std::sin(fAxisMax[phiaxis]);
00449 z = 0;
00450 SetCorner(sC0Min1Max, x, y, z);
00451
00452 } else {
00453 std::ostringstream message;
00454 message << "Feature NOT implemented !" << G4endl
00455 << " fAxis[0] = " << fAxis[0] << G4endl
00456 << " fAxis[1] = " << fAxis[1];
00457 G4Exception("G4TwistTubsFlatSide::SetCorners()", "GeomSolids0001",
00458 FatalException, message);
00459 }
00460 }
00461
00462
00463
00464
00465 void G4TwistTubsFlatSide::SetBoundaries()
00466 {
00467
00468
00469
00470 if (fAxis[0] == kRho && fAxis[1] == kPhi) {
00471
00472 G4ThreeVector direction;
00473
00474 direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
00475 direction = direction.unit();
00476 SetBoundary(sAxis0 & (sAxisPhi | sAxisMin), direction,
00477 GetCorner(sC0Min1Min), sAxisPhi);
00478
00479
00480 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
00481 direction = direction.unit();
00482 SetBoundary(sAxis0 & (sAxisPhi | sAxisMax), direction,
00483 GetCorner(sC0Max1Min), sAxisPhi);
00484
00485
00486 direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
00487 direction = direction.unit();
00488 SetBoundary(sAxis1 & (sAxisRho | sAxisMin), direction,
00489 GetCorner(sC0Min1Min), sAxisRho);
00490
00491
00492 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
00493 direction = direction.unit();
00494 SetBoundary(sAxis1 & (sAxisRho | sAxisMax), direction,
00495 GetCorner(sC0Min1Max), sAxisPhi);
00496 } else {
00497 std::ostringstream message;
00498 message << "Feature NOT implemented !" << G4endl
00499 << " fAxis[0] = " << fAxis[0] << G4endl
00500 << " fAxis[1] = " << fAxis[1];
00501 G4Exception("G4TwistTubsFlatSide::SetBoundaries()", "GeomSolids0001",
00502 FatalException, message);
00503 }
00504 }
00505
00506
00507
00508
00509 void G4TwistTubsFlatSide::GetFacets( G4int k, G4int n, G4double xyz[][3],
00510 G4int faces[][4], G4int iside )
00511 {
00512
00513 G4ThreeVector p ;
00514
00515 G4double rmin = fAxisMin[0] ;
00516 G4double rmax = fAxisMax[0] ;
00517 G4double phimin, phimax ;
00518
00519 G4double r,phi ;
00520
00521 G4int i,j ;
00522
00523 G4int nnode,nface ;
00524
00525 for ( i = 0 ; i<n ; i++ ) {
00526
00527 r = rmin + i*(rmax-rmin)/(n-1) ;
00528
00529 phimin = GetBoundaryMin(r) ;
00530 phimax = GetBoundaryMax(r) ;
00531
00532 for ( j = 0 ; j<k ; j++ )
00533 {
00534 phi = phimin + j*(phimax-phimin)/(k-1) ;
00535
00536 nnode = GetNode(i,j,k,n,iside) ;
00537 p = SurfacePoint(phi,r,true) ;
00538
00539 xyz[nnode][0] = p.x() ;
00540 xyz[nnode][1] = p.y() ;
00541 xyz[nnode][2] = p.z() ;
00542
00543 if ( i<n-1 && j<k-1 ) {
00544
00545 nface = GetFace(i,j,k,n,iside) ;
00546
00547 if (fHandedness < 0) {
00548 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1) * ( GetNode(i ,j ,k,n,iside)+1) ;
00549 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1) * ( GetNode(i ,j+1,k,n,iside)+1) ;
00550 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1) * ( GetNode(i+1,j+1,k,n,iside)+1) ;
00551 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1) * ( GetNode(i+1,j ,k,n,iside)+1) ;
00552 } else {
00553 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1) * ( GetNode(i ,j ,k,n,iside)+1) ;
00554 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1) * ( GetNode(i+1,j ,k,n,iside)+1) ;
00555 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1) * ( GetNode(i+1,j+1,k,n,iside)+1) ;
00556 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1) * ( GetNode(i ,j+1,k,n,iside)+1) ;
00557
00558 }
00559
00560
00561
00562 }
00563 }
00564 }
00565 }