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