Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // G4TwistTrapParallelSide implementation 26 // G4TwistTrapParallelSide implementation 27 // 27 // 28 // Author: Oliver Link (Oliver.Link@cern.ch) 28 // Author: Oliver Link (Oliver.Link@cern.ch) 29 // ------------------------------------------- 29 // -------------------------------------------------------------------- 30 30 31 #include <cmath> 31 #include <cmath> 32 32 33 #include "G4TwistTrapParallelSide.hh" 33 #include "G4TwistTrapParallelSide.hh" 34 #include "G4PhysicalConstants.hh" 34 #include "G4PhysicalConstants.hh" 35 #include "G4JTPolynomialSolver.hh" 35 #include "G4JTPolynomialSolver.hh" 36 36 37 //============================================ 37 //===================================================================== 38 //* constructors ----------------------------- 38 //* constructors ------------------------------------------------------ 39 39 40 G4TwistTrapParallelSide::G4TwistTrapParallelSi 40 G4TwistTrapParallelSide::G4TwistTrapParallelSide(const G4String& name, 41 G4double PhiTwist, 41 G4double PhiTwist, // twist angle 42 G4double pDz, 42 G4double pDz, // half z lenght 43 G4double pTheta, 43 G4double pTheta, // direction between end planes 44 G4double pPhi, 44 G4double pPhi, // by polar and azimutal angles 45 G4double pDy1, 45 G4double pDy1, // half y length at -pDz 46 G4double pDx1, 46 G4double pDx1, // half x length at -pDz,-pDy 47 G4double pDx2, 47 G4double pDx2, // half x length at -pDz,+pDy 48 G4double pDy2, 48 G4double pDy2, // half y length at +pDz 49 G4double pDx3, 49 G4double pDx3, // half x length at +pDz,-pDy 50 G4double pDx4, 50 G4double pDx4, // half x length at +pDz,+pDy 51 G4double pAlph, 51 G4double pAlph, // tilt angle at +pDz 52 G4double AngleSide 52 G4double AngleSide // parity 53 53 ) 54 : G4VTwistSurface(name) 54 : G4VTwistSurface(name) 55 { 55 { 56 56 57 fAxis[0] = kXAxis; // in local coordinate 57 fAxis[0] = kXAxis; // in local coordinate system 58 fAxis[1] = kZAxis; 58 fAxis[1] = kZAxis; 59 fAxisMin[0] = -kInfinity ; // X Axis bounda 59 fAxisMin[0] = -kInfinity ; // X Axis boundary 60 fAxisMax[0] = kInfinity ; // depends on 60 fAxisMax[0] = kInfinity ; // depends on z !! 61 fAxisMin[1] = -pDz ; // Z Axis bounda 61 fAxisMin[1] = -pDz ; // Z Axis boundary 62 fAxisMax[1] = pDz ; 62 fAxisMax[1] = pDz ; 63 63 64 fDx1 = pDx1 ; 64 fDx1 = pDx1 ; 65 fDx2 = pDx2 ; 65 fDx2 = pDx2 ; 66 fDx3 = pDx3 ; 66 fDx3 = pDx3 ; 67 fDx4 = pDx4 ; 67 fDx4 = pDx4 ; 68 68 69 fDy1 = pDy1 ; 69 fDy1 = pDy1 ; 70 fDy2 = pDy2 ; 70 fDy2 = pDy2 ; 71 71 72 fDz = pDz ; 72 fDz = pDz ; 73 73 74 fAlph = pAlph ; 74 fAlph = pAlph ; 75 fTAlph = std::tan(fAlph) ; 75 fTAlph = std::tan(fAlph) ; 76 76 77 fTheta = pTheta ; 77 fTheta = pTheta ; 78 fPhi = pPhi ; 78 fPhi = pPhi ; 79 79 80 // precalculate frequently used parameters 80 // precalculate frequently used parameters 81 // 81 // 82 fDx4plus2 = fDx4 + fDx2 ; 82 fDx4plus2 = fDx4 + fDx2 ; 83 fDx4minus2 = fDx4 - fDx2 ; 83 fDx4minus2 = fDx4 - fDx2 ; 84 fDx3plus1 = fDx3 + fDx1 ; 84 fDx3plus1 = fDx3 + fDx1 ; 85 fDx3minus1 = fDx3 - fDx1 ; 85 fDx3minus1 = fDx3 - fDx1 ; 86 fDy2plus1 = fDy2 + fDy1 ; 86 fDy2plus1 = fDy2 + fDy1 ; 87 fDy2minus1 = fDy2 - fDy1 ; 87 fDy2minus1 = fDy2 - fDy1 ; 88 88 89 fa1md1 = 2*fDx2 - 2*fDx1 ; 89 fa1md1 = 2*fDx2 - 2*fDx1 ; 90 fa2md2 = 2*fDx4 - 2*fDx3 ; 90 fa2md2 = 2*fDx4 - 2*fDx3 ; 91 91 92 fPhiTwist = PhiTwist ; // dphi 92 fPhiTwist = PhiTwist ; // dphi 93 fAngleSide = AngleSide ; // 0,90,180,270 de 93 fAngleSide = AngleSide ; // 0,90,180,270 deg 94 94 95 fdeltaX = 2*fDz*std::tan(fTheta)*std::cos(fP 95 fdeltaX = 2*fDz*std::tan(fTheta)*std::cos(fPhi); // dx in surface equation 96 fdeltaY = 2*fDz*std::tan(fTheta)*std::sin(fP 96 fdeltaY = 2*fDz*std::tan(fTheta)*std::sin(fPhi); // dy in surface equation 97 97 98 fRot.rotateZ( AngleSide ) ; 98 fRot.rotateZ( AngleSide ) ; 99 99 100 fTrans.set(0, 0, 0); // No Translation 100 fTrans.set(0, 0, 0); // No Translation 101 fIsValidNorm = false; 101 fIsValidNorm = false; 102 102 103 SetCorners() ; 103 SetCorners() ; 104 SetBoundaries() ; 104 SetBoundaries() ; 105 } 105 } 106 106 107 //============================================ 107 //===================================================================== 108 //* Fake default constructor ----------------- 108 //* Fake default constructor ------------------------------------------ 109 109 110 G4TwistTrapParallelSide::G4TwistTrapParallelSi 110 G4TwistTrapParallelSide::G4TwistTrapParallelSide( __void__& a ) 111 : G4VTwistSurface(a) << 111 : G4VTwistSurface(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.), >> 112 fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fAlph(0.), fTAlph(0.), fPhiTwist(0.), >> 113 fAngleSide(0.), fdeltaX(0.), fdeltaY(0.), fDx4plus2(0.), fDx4minus2(0.), >> 114 fDx3plus1(0.), fDx3minus1(0.), fDy2plus1(0.), fDy2minus1(0.), fa1md1(0.), >> 115 fa2md2(0.) 112 { 116 { 113 } 117 } 114 118 115 //============================================ 119 //===================================================================== 116 //* destructor ------------------------------- 120 //* destructor -------------------------------------------------------- 117 121 118 G4TwistTrapParallelSide::~G4TwistTrapParallelS << 122 G4TwistTrapParallelSide::~G4TwistTrapParallelSide() >> 123 { >> 124 } 119 125 120 //============================================ 126 //===================================================================== 121 //* GetNormal -------------------------------- 127 //* GetNormal --------------------------------------------------------- 122 128 123 G4ThreeVector G4TwistTrapParallelSide::GetNorm 129 G4ThreeVector G4TwistTrapParallelSide::GetNormal(const G4ThreeVector& tmpxx, 124 130 G4bool isGlobal) 125 { 131 { 126 // GetNormal returns a normal vector at a s 132 // GetNormal returns a normal vector at a surface (or very close 127 // to surface) point at tmpxx. 133 // to surface) point at tmpxx. 128 // If isGlobal=true, it returns the normal 134 // If isGlobal=true, it returns the normal in global coordinate. 129 // 135 // 130 136 131 G4ThreeVector xx; 137 G4ThreeVector xx; 132 if (isGlobal) 138 if (isGlobal) 133 { 139 { 134 xx = ComputeLocalPoint(tmpxx); 140 xx = ComputeLocalPoint(tmpxx); 135 if ((xx - fCurrentNormal.p).mag() < 0.5 141 if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance) 136 { 142 { 137 return ComputeGlobalDirection(fCurren 143 return ComputeGlobalDirection(fCurrentNormal.normal); 138 } 144 } 139 } 145 } 140 else 146 else 141 { 147 { 142 xx = tmpxx; 148 xx = tmpxx; 143 if (xx == fCurrentNormal.p) 149 if (xx == fCurrentNormal.p) 144 { 150 { 145 return fCurrentNormal.normal; 151 return fCurrentNormal.normal; 146 } 152 } 147 } 153 } 148 154 149 G4double phi ; 155 G4double phi ; 150 G4double u ; 156 G4double u ; 151 157 152 GetPhiUAtX(xx,phi,u) ; // phi,u for point 158 GetPhiUAtX(xx,phi,u) ; // phi,u for point xx close to surface 153 159 154 G4ThreeVector normal = NormAng(phi,u) ; / 160 G4ThreeVector normal = NormAng(phi,u) ; // the normal vector at phi,u 155 161 156 #ifdef G4TWISTDEBUG 162 #ifdef G4TWISTDEBUG 157 G4cout << "normal vector = " << normal << 163 G4cout << "normal vector = " << normal << G4endl ; 158 G4cout << "phi = " << phi << " , u = " << u 164 G4cout << "phi = " << phi << " , u = " << u << G4endl ; 159 #endif 165 #endif 160 166 161 // normal = normal/normal.mag() ; 167 // normal = normal/normal.mag() ; 162 168 163 if (isGlobal) 169 if (isGlobal) 164 { 170 { 165 fCurrentNormal.normal = ComputeGlobalDir 171 fCurrentNormal.normal = ComputeGlobalDirection(normal.unit()); 166 } 172 } 167 else 173 else 168 { 174 { 169 fCurrentNormal.normal = normal.unit(); 175 fCurrentNormal.normal = normal.unit(); 170 } 176 } 171 return fCurrentNormal.normal; 177 return fCurrentNormal.normal; 172 } 178 } 173 179 174 //============================================ 180 //===================================================================== 175 //* DistanceToSurface ------------------------ 181 //* DistanceToSurface ------------------------------------------------- 176 182 177 G4int G4TwistTrapParallelSide::DistanceToSurfa 183 G4int G4TwistTrapParallelSide::DistanceToSurface(const G4ThreeVector& gp, 178 184 const G4ThreeVector& gv, 179 185 G4ThreeVector gxx[], 180 186 G4double distance[], 181 187 G4int areacode[], 182 188 G4bool isvalid[], 183 189 EValidate validate) 184 { 190 { 185 static const G4double pihalf = pi/2 ; 191 static const G4double pihalf = pi/2 ; 186 const G4double ctol = 0.5 * kCarTolerance; 192 const G4double ctol = 0.5 * kCarTolerance; 187 193 188 G4bool IsParallel = false ; 194 G4bool IsParallel = false ; 189 G4bool IsConverged = false ; 195 G4bool IsConverged = false ; 190 196 191 G4int nxx = 0 ; // number of physical solut 197 G4int nxx = 0 ; // number of physical solutions 192 198 193 fCurStatWithV.ResetfDone(validate, &gp, &gv) 199 fCurStatWithV.ResetfDone(validate, &gp, &gv); 194 200 195 if (fCurStatWithV.IsDone()) 201 if (fCurStatWithV.IsDone()) 196 { 202 { 197 for (G4int i=0; i<fCurStatWithV.GetNXX(); 203 for (G4int i=0; i<fCurStatWithV.GetNXX(); ++i) 198 { 204 { 199 gxx[i] = fCurStatWithV.GetXX(i); 205 gxx[i] = fCurStatWithV.GetXX(i); 200 distance[i] = fCurStatWithV.GetDistance( 206 distance[i] = fCurStatWithV.GetDistance(i); 201 areacode[i] = fCurStatWithV.GetAreacode( 207 areacode[i] = fCurStatWithV.GetAreacode(i); 202 isvalid[i] = fCurStatWithV.IsValid(i); 208 isvalid[i] = fCurStatWithV.IsValid(i); 203 } 209 } 204 return fCurStatWithV.GetNXX(); 210 return fCurStatWithV.GetNXX(); 205 } 211 } 206 else // initialize 212 else // initialize 207 { 213 { 208 for (G4int i=0; i<G4VSURFACENXX ; ++i) 214 for (G4int i=0; i<G4VSURFACENXX ; ++i) 209 { 215 { 210 distance[i] = kInfinity; 216 distance[i] = kInfinity; 211 areacode[i] = sOutside; 217 areacode[i] = sOutside; 212 isvalid[i] = false; 218 isvalid[i] = false; 213 gxx[i].set(kInfinity, kInfinity, kInfini 219 gxx[i].set(kInfinity, kInfinity, kInfinity); 214 } 220 } 215 } 221 } 216 222 217 G4ThreeVector p = ComputeLocalPoint(gp); 223 G4ThreeVector p = ComputeLocalPoint(gp); 218 G4ThreeVector v = ComputeLocalDirection(gv); 224 G4ThreeVector v = ComputeLocalDirection(gv); 219 225 220 #ifdef G4TWISTDEBUG 226 #ifdef G4TWISTDEBUG 221 G4cout << "Local point p = " << p << G4endl 227 G4cout << "Local point p = " << p << G4endl ; 222 G4cout << "Local direction v = " << v << G4e 228 G4cout << "Local direction v = " << v << G4endl ; 223 #endif 229 #endif 224 230 225 G4double phi,u ; // parameters 231 G4double phi,u ; // parameters 226 232 227 // temporary variables 233 // temporary variables 228 234 229 G4double tmpdist = kInfinity ; 235 G4double tmpdist = kInfinity ; 230 G4ThreeVector tmpxx; 236 G4ThreeVector tmpxx; 231 G4int tmpareacode = sOutside ; 237 G4int tmpareacode = sOutside ; 232 G4bool tmpisvalid = false ; 238 G4bool tmpisvalid = false ; 233 239 234 std::vector<Intersection> xbuf ; 240 std::vector<Intersection> xbuf ; 235 Intersection xbuftmp ; 241 Intersection xbuftmp ; 236 242 237 // prepare some variables for the intersecti 243 // prepare some variables for the intersection finder 238 244 239 G4double L = 2*fDz ; 245 G4double L = 2*fDz ; 240 246 241 G4double phixz = fPhiTwist * ( p.x() * v.z() 247 G4double phixz = fPhiTwist * ( p.x() * v.z() - p.z() * v.x() ) ; 242 G4double phiyz = fPhiTwist * ( p.y() * v.z() 248 G4double phiyz = fPhiTwist * ( p.y() * v.z() - p.z() * v.y() ) ; 243 249 244 // special case vz = 0 250 // special case vz = 0 245 251 246 if ( v.z() == 0. ) 252 if ( v.z() == 0. ) 247 { 253 { 248 if ( std::fabs(p.z()) <= L ) // inte 254 if ( std::fabs(p.z()) <= L ) // intersection possible in z 249 { 255 { 250 phi = p.z() * fPhiTwist / L ; // phi is 256 phi = p.z() * fPhiTwist / L ; // phi is determined by the z-position 251 257 252 u = (2*(fdeltaY*phi*v.x() - fPhiTwist*p. 258 u = (2*(fdeltaY*phi*v.x() - fPhiTwist*p.y()*v.x() - fdeltaX*phi*v.y() 253 + fPhiTwist*p.x()*v.y()) + (fDy2 259 + fPhiTwist*p.x()*v.y()) + (fDy2plus1*fPhiTwist 254 + 2*fDy2minus1*phi)*(v.x()*std:: 260 + 2*fDy2minus1*phi)*(v.x()*std::cos(phi) + v.y()*std::sin(phi))) 255 / (2.* fPhiTwist*(v.y()*std::cos(phi) 261 / (2.* fPhiTwist*(v.y()*std::cos(phi) - v.x()*std::sin(phi))); 256 262 257 xbuftmp.phi = phi ; 263 xbuftmp.phi = phi ; 258 xbuftmp.u = u ; 264 xbuftmp.u = u ; 259 xbuftmp.areacode = sOutside ; 265 xbuftmp.areacode = sOutside ; 260 xbuftmp.distance = kInfinity ; 266 xbuftmp.distance = kInfinity ; 261 xbuftmp.isvalid = false ; 267 xbuftmp.isvalid = false ; 262 268 263 xbuf.push_back(xbuftmp) ; // store it to 269 xbuf.push_back(xbuftmp) ; // store it to xbuf 264 } 270 } 265 else // no in 271 else // no intersection possible 266 { 272 { 267 distance[0] = kInfinity; 273 distance[0] = kInfinity; 268 gxx[0].set(kInfinity,kInfinity,kInfinity 274 gxx[0].set(kInfinity,kInfinity,kInfinity); 269 isvalid[0] = false ; 275 isvalid[0] = false ; 270 areacode[0] = sOutside ; 276 areacode[0] = sOutside ; 271 fCurStatWithV.SetCurrentStatus(0, gxx[0] 277 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], 272 areacode[ 278 areacode[0], isvalid[0], 273 0, valida 279 0, validate, &gp, &gv); 274 280 275 return 0; 281 return 0; 276 } // end std::fabs(p.z() <= L 282 } // end std::fabs(p.z() <= L 277 } // end v.z() == 0 283 } // end v.z() == 0 278 else // general solution for non-zero vz 284 else // general solution for non-zero vz 279 { 285 { 280 G4double c[9],srd[8],si[8] ; 286 G4double c[9],srd[8],si[8] ; 281 287 282 c[8] = -3600*(-2*phiyz + fDy2plus1*fPhiTwi 288 c[8] = -3600*(-2*phiyz + fDy2plus1*fPhiTwist*v.z()) ; 283 c[7] = -7200*(phixz - 2*fDz*v.y() + (fdelt 289 c[7] = -7200*(phixz - 2*fDz*v.y() + (fdeltaY + fDy2minus1)*v.z()) ; 284 c[6] = 120*(-52*phiyz - 120*fDz*v.x() + 60 290 c[6] = 120*(-52*phiyz - 120*fDz*v.x() + 60*fdeltaX*v.z() 285 + 11*fDy2plus1*fPhiTwist*v.z()) ; 291 + 11*fDy2plus1*fPhiTwist*v.z()) ; 286 c[5] = 240*(16*phixz - 52*fDz*v.y() + 26*f 292 c[5] = 240*(16*phixz - 52*fDz*v.y() + 26*fdeltaY*v.z() 287 + 11*fDy2minus1*v.z()) ; 293 + 11*fDy2minus1*v.z()) ; 288 c[4] = 12*(127*phiyz + 640*fDz*v.x() - 320 294 c[4] = 12*(127*phiyz + 640*fDz*v.x() - 320*fdeltaX*v.z() 289 + 4*fDy2plus1*fPhiTwist*v.z()) ; 295 + 4*fDy2plus1*fPhiTwist*v.z()) ; 290 c[3] = -404*phixz + 3048*fDz*v.y() - 1524* 296 c[3] = -404*phixz + 3048*fDz*v.y() - 1524*fdeltaY*v.z() 291 + 96*fDy2minus1*v.z() ; 297 + 96*fDy2minus1*v.z() ; 292 c[2] = -72*phiyz + 404*(-2*fDz*v.x() + fde 298 c[2] = -72*phiyz + 404*(-2*fDz*v.x() + fdeltaX*v.z()) ; 293 c[1] = 12*(phixz - 12*fDz*v.y() + 6*fdelta 299 c[1] = 12*(phixz - 12*fDz*v.y() + 6*fdeltaY*v.z()) ; 294 c[0] = 24*fDz*v.x() - 12*fdeltaX*v.z() ; 300 c[0] = 24*fDz*v.x() - 12*fdeltaX*v.z() ; 295 301 296 302 297 #ifdef G4TWISTDEBUG 303 #ifdef G4TWISTDEBUG 298 G4cout << "coef = " << c[0] << " " 304 G4cout << "coef = " << c[0] << " " 299 << c[1] << " " 305 << c[1] << " " 300 << c[2] << " " 306 << c[2] << " " 301 << c[3] << " " 307 << c[3] << " " 302 << c[4] << " " 308 << c[4] << " " 303 << c[5] << " " 309 << c[5] << " " 304 << c[6] << " " 310 << c[6] << " " 305 << c[7] << " " 311 << c[7] << " " 306 << c[8] << G4endl ; 312 << c[8] << G4endl ; 307 #endif 313 #endif 308 314 309 G4JTPolynomialSolver trapEq ; 315 G4JTPolynomialSolver trapEq ; 310 G4int num = trapEq.FindRoots(c,8,srd,si); 316 G4int num = trapEq.FindRoots(c,8,srd,si); 311 317 312 for (G4int i = 0 ; i<num ; ++i ) // loop 318 for (G4int i = 0 ; i<num ; ++i ) // loop over all mathematical solutions 313 { 319 { 314 if ( si[i]==0.0 ) // only real solutio 320 if ( si[i]==0.0 ) // only real solutions 315 { 321 { 316 #ifdef G4TWISTDEBUG 322 #ifdef G4TWISTDEBUG 317 G4cout << "Solution " << i << " : " << 323 G4cout << "Solution " << i << " : " << srd[i] << G4endl ; 318 #endif 324 #endif 319 phi = std::fmod(srd[i] , pihalf) ; 325 phi = std::fmod(srd[i] , pihalf) ; 320 u = (1/std::cos(phi)*(2*phixz + 4*fDz* 326 u = (1/std::cos(phi)*(2*phixz + 4*fDz*phi*v.x() 321 - 2*fdeltaX*phi*v.z() + (fDy2plus1*f 327 - 2*fdeltaX*phi*v.z() + (fDy2plus1*fPhiTwist 322 + 2*fDy2minus1*phi)*v.z()* std::sin( 328 + 2*fDy2minus1*phi)*v.z()* std::sin(phi)))/(2.*fPhiTwist*v.z()) ; 323 329 324 xbuftmp.phi = phi ; 330 xbuftmp.phi = phi ; 325 xbuftmp.u = u ; 331 xbuftmp.u = u ; 326 xbuftmp.areacode = sOutside ; 332 xbuftmp.areacode = sOutside ; 327 xbuftmp.distance = kInfinity ; 333 xbuftmp.distance = kInfinity ; 328 xbuftmp.isvalid = false ; 334 xbuftmp.isvalid = false ; 329 335 330 xbuf.push_back(xbuftmp) ; // store it 336 xbuf.push_back(xbuftmp) ; // store it to xbuf 331 337 332 #ifdef G4TWISTDEBUG 338 #ifdef G4TWISTDEBUG 333 G4cout << "solution " << i << " = " << 339 G4cout << "solution " << i << " = " << phi << " , " << u << G4endl ; 334 #endif 340 #endif 335 341 336 } // end if real solution 342 } // end if real solution 337 } // end loop i 343 } // end loop i 338 } // end general case 344 } // end general case 339 345 340 nxx = (G4int)xbuf.size() ; // save the numb << 346 nxx = xbuf.size() ; // save the number of solutions 341 347 342 G4ThreeVector xxonsurface ; // point 348 G4ThreeVector xxonsurface ; // point on surface 343 G4ThreeVector surfacenormal ; // normal 349 G4ThreeVector surfacenormal ; // normal vector 344 G4double deltaX ; // distance between inter 350 G4double deltaX ; // distance between intersection point and point on surface 345 G4double theta ; // angle 351 G4double theta ; // angle between track and surfacenormal 346 G4double factor ; // a scal 352 G4double factor ; // a scaling factor 347 G4int maxint = 30 ; // number 353 G4int maxint = 30 ; // number of iterations 348 354 349 355 350 for (auto & k : xbuf) << 356 for ( size_t k = 0 ; k<xbuf.size() ; ++k ) 351 { 357 { 352 #ifdef G4TWISTDEBUG 358 #ifdef G4TWISTDEBUG 353 G4cout << "Solution " << k << " : " 359 G4cout << "Solution " << k << " : " 354 << "reconstructed phiR = " << xbuf[ 360 << "reconstructed phiR = " << xbuf[k].phi 355 << ", uR = " << xbuf[k].u << G4endl 361 << ", uR = " << xbuf[k].u << G4endl ; 356 #endif 362 #endif 357 363 358 phi = k.phi ; // get the stored values fo << 364 phi = xbuf[k].phi ; // get the stored values for phi and u 359 u = k.u ; << 365 u = xbuf[k].u ; 360 366 361 IsConverged = false ; // no convergence 367 IsConverged = false ; // no convergence at the beginning 362 368 363 for ( G4int i = 1 ; i<maxint ; ++i ) 369 for ( G4int i = 1 ; i<maxint ; ++i ) 364 { 370 { 365 xxonsurface = SurfacePoint(phi,u) ; 371 xxonsurface = SurfacePoint(phi,u) ; 366 surfacenormal = NormAng(phi,u) ; 372 surfacenormal = NormAng(phi,u) ; 367 tmpdist = DistanceToPlaneWithV(p, v, xxo 373 tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx); 368 deltaX = ( tmpxx - xxonsurface ).mag() ; 374 deltaX = ( tmpxx - xxonsurface ).mag() ; 369 theta = std::fabs(std::acos(v*surfacenor 375 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ; 370 if ( theta < 0.001 ) 376 if ( theta < 0.001 ) 371 { 377 { 372 factor = 50 ; 378 factor = 50 ; 373 IsParallel = true ; 379 IsParallel = true ; 374 } 380 } 375 else 381 else 376 { 382 { 377 factor = 1 ; 383 factor = 1 ; 378 } 384 } 379 385 380 #ifdef G4TWISTDEBUG 386 #ifdef G4TWISTDEBUG 381 G4cout << "Step i = " << i << ", distanc 387 G4cout << "Step i = " << i << ", distance = " 382 << tmpdist << ", " << deltaX << G 388 << tmpdist << ", " << deltaX << G4endl ; 383 G4cout << "X = " << tmpxx << G4endl ; 389 G4cout << "X = " << tmpxx << G4endl ; 384 #endif 390 #endif 385 391 386 GetPhiUAtX(tmpxx, phi, u); // new point 392 GetPhiUAtX(tmpxx, phi, u); // new point xx is accepted and phi/u replaced 387 393 388 #ifdef G4TWISTDEBUG 394 #ifdef G4TWISTDEBUG 389 G4cout << "approximated phi = " << phi < 395 G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ; 390 #endif 396 #endif 391 397 392 if ( deltaX <= factor*ctol ) { IsConverg 398 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; } 393 } // end iterative loop (i) 399 } // end iterative loop (i) 394 400 395 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 401 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ; 396 402 397 #ifdef G4TWISTDEBUG 403 #ifdef G4TWISTDEBUG 398 G4cout << "refined solution " << phi << " 404 G4cout << "refined solution " << phi << " , " << u << G4endl ; 399 G4cout << "distance = " << tmpdist << G4en 405 G4cout << "distance = " << tmpdist << G4endl ; 400 G4cout << "local X = " << tmpxx << G4endl 406 G4cout << "local X = " << tmpxx << G4endl ; 401 #endif 407 #endif 402 408 403 tmpisvalid = false ; // init 409 tmpisvalid = false ; // init 404 410 405 if ( IsConverged ) 411 if ( IsConverged ) 406 { 412 { 407 if (validate == kValidateWithTol) 413 if (validate == kValidateWithTol) 408 { 414 { 409 tmpareacode = GetAreaCode(tmpxx); 415 tmpareacode = GetAreaCode(tmpxx); 410 if (!IsOutside(tmpareacode)) 416 if (!IsOutside(tmpareacode)) 411 { 417 { 412 if (tmpdist >= 0) tmpisvalid = true; 418 if (tmpdist >= 0) tmpisvalid = true; 413 } 419 } 414 } 420 } 415 else if (validate == kValidateWithoutTol 421 else if (validate == kValidateWithoutTol) 416 { 422 { 417 tmpareacode = GetAreaCode(tmpxx, false 423 tmpareacode = GetAreaCode(tmpxx, false); 418 if (IsInside(tmpareacode)) 424 if (IsInside(tmpareacode)) 419 { 425 { 420 if (tmpdist >= 0) tmpisvalid = true; 426 if (tmpdist >= 0) tmpisvalid = true; 421 } 427 } 422 } 428 } 423 else // kDontValidate 429 else // kDontValidate 424 { 430 { 425 G4Exception("G4TwistTrapParallelSide:: 431 G4Exception("G4TwistTrapParallelSide::DistanceToSurface()", 426 "GeomSolids0001", FatalExc 432 "GeomSolids0001", FatalException, 427 "Feature NOT implemented ! 433 "Feature NOT implemented !"); 428 } 434 } 429 } 435 } 430 else 436 else 431 { 437 { 432 tmpdist = kInfinity; // no convergen 438 tmpdist = kInfinity; // no convergence after 10 steps 433 tmpisvalid = false ; // solution is 439 tmpisvalid = false ; // solution is not vaild 434 } 440 } 435 441 436 // store the found values 442 // store the found values 437 k.xx = tmpxx ; << 443 xbuf[k].xx = tmpxx ; 438 k.distance = tmpdist ; << 444 xbuf[k].distance = tmpdist ; 439 k.areacode = tmpareacode ; << 445 xbuf[k].areacode = tmpareacode ; 440 k.isvalid = tmpisvalid ; << 446 xbuf[k].isvalid = tmpisvalid ; 441 } // end loop over physical solutions (vari 447 } // end loop over physical solutions (variable k) 442 448 443 std::sort(xbuf.begin() , xbuf.end(), Distanc 449 std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ; // sorting 444 450 445 #ifdef G4TWISTDEBUG 451 #ifdef G4TWISTDEBUG 446 G4cout << G4endl << "list xbuf after sorting 452 G4cout << G4endl << "list xbuf after sorting : " << G4endl ; 447 G4cout << G4endl << G4endl ; 453 G4cout << G4endl << G4endl ; 448 #endif 454 #endif 449 455 450 // erase identical intersection (within kCar 456 // erase identical intersection (within kCarTolerance) 451 xbuf.erase(std::unique(xbuf.begin(),xbuf.end 457 xbuf.erase(std::unique(xbuf.begin(),xbuf.end(),EqualIntersection),xbuf.end()); 452 458 453 459 454 // add guesses 460 // add guesses 455 461 456 auto nxxtmp = (G4int)xbuf.size() ; << 462 G4int nxxtmp = xbuf.size() ; 457 463 458 if ( nxxtmp<2 || IsParallel ) 464 if ( nxxtmp<2 || IsParallel ) 459 { 465 { 460 // positive end 466 // positive end 461 #ifdef G4TWISTDEBUG 467 #ifdef G4TWISTDEBUG 462 G4cout << "add guess at +z/2 .. " << G4end 468 G4cout << "add guess at +z/2 .. " << G4endl ; 463 #endif 469 #endif 464 470 465 phi = fPhiTwist/2 ; 471 phi = fPhiTwist/2 ; 466 u = 0 ; 472 u = 0 ; 467 473 468 xbuftmp.phi = phi ; 474 xbuftmp.phi = phi ; 469 xbuftmp.u = u ; 475 xbuftmp.u = u ; 470 xbuftmp.areacode = sOutside ; 476 xbuftmp.areacode = sOutside ; 471 xbuftmp.distance = kInfinity ; 477 xbuftmp.distance = kInfinity ; 472 xbuftmp.isvalid = false ; 478 xbuftmp.isvalid = false ; 473 479 474 xbuf.push_back(xbuftmp) ; // store it to 480 xbuf.push_back(xbuftmp) ; // store it to xbuf 475 481 476 #ifdef G4TWISTDEBUG 482 #ifdef G4TWISTDEBUG 477 G4cout << "add guess at -z/2 .. " << G4end 483 G4cout << "add guess at -z/2 .. " << G4endl ; 478 #endif 484 #endif 479 485 480 phi = -fPhiTwist/2 ; 486 phi = -fPhiTwist/2 ; 481 u = 0 ; 487 u = 0 ; 482 488 483 xbuftmp.phi = phi ; 489 xbuftmp.phi = phi ; 484 xbuftmp.u = u ; 490 xbuftmp.u = u ; 485 xbuftmp.areacode = sOutside ; 491 xbuftmp.areacode = sOutside ; 486 xbuftmp.distance = kInfinity ; 492 xbuftmp.distance = kInfinity ; 487 xbuftmp.isvalid = false ; 493 xbuftmp.isvalid = false ; 488 494 489 xbuf.push_back(xbuftmp) ; // store it to 495 xbuf.push_back(xbuftmp) ; // store it to xbuf 490 496 491 for ( std::size_t k = nxxtmp ; k<xbuf.size << 497 for ( size_t k = nxxtmp ; k<xbuf.size() ; ++k ) 492 { 498 { 493 #ifdef G4TWISTDEBUG 499 #ifdef G4TWISTDEBUG 494 G4cout << "Solution " << k << " : " 500 G4cout << "Solution " << k << " : " 495 << "reconstructed phiR = " << xbu 501 << "reconstructed phiR = " << xbuf[k].phi 496 << ", uR = " << xbuf[k].u << G4en 502 << ", uR = " << xbuf[k].u << G4endl ; 497 #endif 503 #endif 498 504 499 phi = xbuf[k].phi ; // get the stored v 505 phi = xbuf[k].phi ; // get the stored values for phi and u 500 u = xbuf[k].u ; 506 u = xbuf[k].u ; 501 507 502 IsConverged = false ; // no convergenc 508 IsConverged = false ; // no convergence at the beginning 503 509 504 for ( G4int i = 1 ; i<maxint ; ++i ) 510 for ( G4int i = 1 ; i<maxint ; ++i ) 505 { 511 { 506 xxonsurface = SurfacePoint(phi,u) ; 512 xxonsurface = SurfacePoint(phi,u) ; 507 surfacenormal = NormAng(phi,u) ; 513 surfacenormal = NormAng(phi,u) ; 508 tmpdist = DistanceToPlaneWithV(p, v, x 514 tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx); 509 deltaX = ( tmpxx - xxonsurface ).mag() 515 deltaX = ( tmpxx - xxonsurface ).mag() ; 510 theta = std::fabs(std::acos(v*surfacen 516 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ; 511 if ( theta < 0.001 ) 517 if ( theta < 0.001 ) 512 { 518 { 513 factor = 50 ; 519 factor = 50 ; 514 } 520 } 515 else 521 else 516 { 522 { 517 factor = 1 ; 523 factor = 1 ; 518 } 524 } 519 525 520 #ifdef G4TWISTDEBUG 526 #ifdef G4TWISTDEBUG 521 G4cout << "Step i = " << i << ", dista 527 G4cout << "Step i = " << i << ", distance = " 522 << tmpdist << ", " << deltaX << 528 << tmpdist << ", " << deltaX << G4endl ; 523 G4cout << "X = " << tmpxx << G4endl ; 529 G4cout << "X = " << tmpxx << G4endl ; 524 #endif 530 #endif 525 531 526 GetPhiUAtX(tmpxx, phi, u) ; // new poi 532 GetPhiUAtX(tmpxx, phi, u) ; // new point xx accepted and phi/u replaced 527 533 528 #ifdef G4TWISTDEBUG 534 #ifdef G4TWISTDEBUG 529 G4cout << "approximated phi = " << phi 535 G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ; 530 #endif 536 #endif 531 537 532 if ( deltaX <= factor*ctol ) { IsConve 538 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; } 533 } // end iterative loop (i) 539 } // end iterative loop (i) 534 540 535 if ( std::fabs(tmpdist)<ctol ) tmpdist = 541 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ; 536 542 537 #ifdef G4TWISTDEBUG 543 #ifdef G4TWISTDEBUG 538 G4cout << "refined solution " << phi << 544 G4cout << "refined solution " << phi << " , " << u << G4endl ; 539 G4cout << "distance = " << tmpdist << G4 545 G4cout << "distance = " << tmpdist << G4endl ; 540 G4cout << "local X = " << tmpxx << G4end 546 G4cout << "local X = " << tmpxx << G4endl ; 541 #endif 547 #endif 542 548 543 tmpisvalid = false ; // init 549 tmpisvalid = false ; // init 544 550 545 if ( IsConverged ) 551 if ( IsConverged ) 546 { 552 { 547 if (validate == kValidateWithTol) 553 if (validate == kValidateWithTol) 548 { 554 { 549 tmpareacode = GetAreaCode(tmpxx); 555 tmpareacode = GetAreaCode(tmpxx); 550 if (!IsOutside(tmpareacode)) 556 if (!IsOutside(tmpareacode)) 551 { 557 { 552 if (tmpdist >= 0) tmpisvalid = tru 558 if (tmpdist >= 0) tmpisvalid = true; 553 } 559 } 554 } 560 } 555 else if (validate == kValidateWithoutT 561 else if (validate == kValidateWithoutTol) 556 { 562 { 557 tmpareacode = GetAreaCode(tmpxx, fal 563 tmpareacode = GetAreaCode(tmpxx, false); 558 if (IsInside(tmpareacode)) 564 if (IsInside(tmpareacode)) 559 { 565 { 560 if (tmpdist >= 0) tmpisvalid = tru 566 if (tmpdist >= 0) tmpisvalid = true; 561 } 567 } 562 } 568 } 563 else // kDontValidate 569 else // kDontValidate 564 { 570 { 565 G4Exception("G4TwistedBoxSide::Dista 571 G4Exception("G4TwistedBoxSide::DistanceToSurface()", 566 "GeomSolids0001", FatalE 572 "GeomSolids0001", FatalException, 567 "Feature NOT implemented 573 "Feature NOT implemented !"); 568 } 574 } 569 575 570 } 576 } 571 else 577 else 572 { 578 { 573 tmpdist = kInfinity; // no converg 579 tmpdist = kInfinity; // no convergence after 10 steps 574 tmpisvalid = false ; // solution i 580 tmpisvalid = false ; // solution is not vaild 575 } 581 } 576 582 577 // store the found values 583 // store the found values 578 xbuf[k].xx = tmpxx ; 584 xbuf[k].xx = tmpxx ; 579 xbuf[k].distance = tmpdist ; 585 xbuf[k].distance = tmpdist ; 580 xbuf[k].areacode = tmpareacode ; 586 xbuf[k].areacode = tmpareacode ; 581 xbuf[k].isvalid = tmpisvalid ; 587 xbuf[k].isvalid = tmpisvalid ; 582 588 583 } // end loop over physical solutions 589 } // end loop over physical solutions 584 } // end less than 2 solutions 590 } // end less than 2 solutions 585 591 586 // sort again 592 // sort again 587 std::sort(xbuf.begin() , xbuf.end(), Distanc 593 std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ; // sorting 588 594 589 // erase identical intersection (within kCar 595 // erase identical intersection (within kCarTolerance) 590 xbuf.erase(std::unique(xbuf.begin(),xbuf.end 596 xbuf.erase(std::unique(xbuf.begin(),xbuf.end(),EqualIntersection),xbuf.end()); 591 597 592 #ifdef G4TWISTDEBUG 598 #ifdef G4TWISTDEBUG 593 G4cout << G4endl << "list xbuf after sorting 599 G4cout << G4endl << "list xbuf after sorting : " << G4endl ; 594 G4cout << G4endl << G4endl ; 600 G4cout << G4endl << G4endl ; 595 #endif 601 #endif 596 602 597 nxx = (G4int)xbuf.size() ; // determine nu << 603 nxx = xbuf.size() ; // determine number of solutions again. 598 604 599 for ( G4int i = 0 ; i<(G4int)xbuf.size() ; + << 605 for ( size_t i = 0 ; i<xbuf.size() ; ++i ) 600 { 606 { 601 distance[i] = xbuf[i].distance; 607 distance[i] = xbuf[i].distance; 602 gxx[i] = ComputeGlobalPoint(xbuf[i].x 608 gxx[i] = ComputeGlobalPoint(xbuf[i].xx); 603 areacode[i] = xbuf[i].areacode ; 609 areacode[i] = xbuf[i].areacode ; 604 isvalid[i] = xbuf[i].isvalid ; 610 isvalid[i] = xbuf[i].isvalid ; 605 611 606 fCurStatWithV.SetCurrentStatus(i, gxx[i], 612 fCurStatWithV.SetCurrentStatus(i, gxx[i], distance[i], areacode[i], 607 isvalid[i 613 isvalid[i], nxx, validate, &gp, &gv); 608 #ifdef G4TWISTDEBUG 614 #ifdef G4TWISTDEBUG 609 G4cout << "element Nr. " << i 615 G4cout << "element Nr. " << i 610 << ", local Intersection = " << xbu 616 << ", local Intersection = " << xbuf[i].xx 611 << ", distance = " << xbuf[i].dista 617 << ", distance = " << xbuf[i].distance 612 << ", u = " << xbuf[i].u 618 << ", u = " << xbuf[i].u 613 << ", phi = " << xbuf[i].phi 619 << ", phi = " << xbuf[i].phi 614 << ", isvalid = " << xbuf[i].isvali 620 << ", isvalid = " << xbuf[i].isvalid 615 << G4endl ; 621 << G4endl ; 616 #endif 622 #endif 617 } // end for( i ) loop 623 } // end for( i ) loop 618 624 619 #ifdef G4TWISTDEBUG 625 #ifdef G4TWISTDEBUG 620 G4cout << "G4TwistTrapParallelSide finished 626 G4cout << "G4TwistTrapParallelSide finished " << G4endl ; 621 G4cout << nxx << " possible physical solutio 627 G4cout << nxx << " possible physical solutions found" << G4endl ; 622 for ( G4int k= 0 ; k< nxx ; ++k ) 628 for ( G4int k= 0 ; k< nxx ; ++k ) 623 { 629 { 624 G4cout << "global intersection Point found 630 G4cout << "global intersection Point found: " << gxx[k] << G4endl ; 625 G4cout << "distance = " << distance[k] << 631 G4cout << "distance = " << distance[k] << G4endl ; 626 G4cout << "isvalid = " << isvalid[k] << G4 632 G4cout << "isvalid = " << isvalid[k] << G4endl ; 627 } 633 } 628 #endif 634 #endif 629 635 630 return nxx ; 636 return nxx ; 631 } 637 } 632 638 633 //============================================ 639 //===================================================================== 634 //* DistanceToSurface ------------------------ 640 //* DistanceToSurface ------------------------------------------------- 635 641 636 G4int G4TwistTrapParallelSide::DistanceToSurfa 642 G4int G4TwistTrapParallelSide::DistanceToSurface(const G4ThreeVector& gp, 637 643 G4ThreeVector gxx[], 638 644 G4double distance[], 639 645 G4int areacode[]) 640 { 646 { 641 const G4double ctol = 0.5 * kCarTolerance; 647 const G4double ctol = 0.5 * kCarTolerance; 642 648 643 fCurStat.ResetfDone(kDontValidate, &gp); 649 fCurStat.ResetfDone(kDontValidate, &gp); 644 650 645 if (fCurStat.IsDone()) 651 if (fCurStat.IsDone()) 646 { 652 { 647 for (G4int i=0; i<fCurStat.GetNXX(); ++i 653 for (G4int i=0; i<fCurStat.GetNXX(); ++i) 648 { 654 { 649 gxx[i] = fCurStat.GetXX(i); 655 gxx[i] = fCurStat.GetXX(i); 650 distance[i] = fCurStat.GetDistance(i) 656 distance[i] = fCurStat.GetDistance(i); 651 areacode[i] = fCurStat.GetAreacode(i) 657 areacode[i] = fCurStat.GetAreacode(i); 652 } 658 } 653 return fCurStat.GetNXX(); 659 return fCurStat.GetNXX(); 654 } 660 } 655 else // initialize 661 else // initialize 656 { 662 { 657 for (G4int i=0; i<G4VSURFACENXX; ++i) 663 for (G4int i=0; i<G4VSURFACENXX; ++i) 658 { 664 { 659 distance[i] = kInfinity; 665 distance[i] = kInfinity; 660 areacode[i] = sOutside; 666 areacode[i] = sOutside; 661 gxx[i].set(kInfinity, kInfinity, kInf 667 gxx[i].set(kInfinity, kInfinity, kInfinity); 662 } 668 } 663 } 669 } 664 670 665 G4ThreeVector p = ComputeLocalPoint(gp); 671 G4ThreeVector p = ComputeLocalPoint(gp); 666 G4ThreeVector xx; // intersection point 672 G4ThreeVector xx; // intersection point 667 G4ThreeVector xxonsurface ; // interpolated 673 G4ThreeVector xxonsurface ; // interpolated intersection point 668 674 669 // the surfacenormal at that surface point 675 // the surfacenormal at that surface point 670 G4double phiR = 0 ; // 676 G4double phiR = 0 ; // 671 G4double uR = 0 ; 677 G4double uR = 0 ; 672 678 673 G4ThreeVector surfacenormal ; 679 G4ThreeVector surfacenormal ; 674 G4double deltaX ; 680 G4double deltaX ; 675 681 676 G4int maxint = 20 ; 682 G4int maxint = 20 ; 677 683 678 for ( G4int i = 1 ; i<maxint ; ++i ) 684 for ( G4int i = 1 ; i<maxint ; ++i ) 679 { 685 { 680 xxonsurface = SurfacePoint(phiR,uR) ; 686 xxonsurface = SurfacePoint(phiR,uR) ; 681 surfacenormal = NormAng(phiR,uR) ; 687 surfacenormal = NormAng(phiR,uR) ; 682 distance[0] = DistanceToPlane(p, xxonsurf 688 distance[0] = DistanceToPlane(p, xxonsurface, surfacenormal, xx); // new XX 683 deltaX = ( xx - xxonsurface ).mag() ; 689 deltaX = ( xx - xxonsurface ).mag() ; 684 690 685 #ifdef G4TWISTDEBUG 691 #ifdef G4TWISTDEBUG 686 G4cout << "i = " << i << ", distance = " 692 G4cout << "i = " << i << ", distance = " 687 << distance[0] << ", " << deltaX < 693 << distance[0] << ", " << deltaX << G4endl ; 688 G4cout << "X = " << xx << G4endl ; 694 G4cout << "X = " << xx << G4endl ; 689 #endif 695 #endif 690 696 691 // the new point xx is accepted and phi/p 697 // the new point xx is accepted and phi/psi replaced 692 GetPhiUAtX(xx, phiR, uR) ; 698 GetPhiUAtX(xx, phiR, uR) ; 693 699 694 if ( deltaX <= ctol ) { break ; } 700 if ( deltaX <= ctol ) { break ; } 695 } 701 } 696 702 697 // check validity of solution ( valid phi,p 703 // check validity of solution ( valid phi,psi ) 698 704 699 G4double halfphi = 0.5*fPhiTwist ; 705 G4double halfphi = 0.5*fPhiTwist ; 700 G4double uMax = GetBoundaryMax(phiR) ; 706 G4double uMax = GetBoundaryMax(phiR) ; 701 G4double uMin = GetBoundaryMin(phiR) ; 707 G4double uMin = GetBoundaryMin(phiR) ; 702 708 703 if ( phiR > halfphi ) phiR = halfphi ; 709 if ( phiR > halfphi ) phiR = halfphi ; 704 if ( phiR < -halfphi ) phiR = -halfphi ; 710 if ( phiR < -halfphi ) phiR = -halfphi ; 705 if ( uR > uMax ) uR = uMax ; 711 if ( uR > uMax ) uR = uMax ; 706 if ( uR < uMin ) uR = uMin ; 712 if ( uR < uMin ) uR = uMin ; 707 713 708 xxonsurface = SurfacePoint(phiR,uR) ; 714 xxonsurface = SurfacePoint(phiR,uR) ; 709 distance[0] = ( p - xx ).mag() ; 715 distance[0] = ( p - xx ).mag() ; 710 if ( distance[0] <= ctol ) { distance[0] = 716 if ( distance[0] <= ctol ) { distance[0] = 0 ; } 711 717 712 // end of validity 718 // end of validity 713 719 714 #ifdef G4TWISTDEBUG 720 #ifdef G4TWISTDEBUG 715 G4cout << "refined solution " << phiR << " 721 G4cout << "refined solution " << phiR << " , " << uR << " , " << G4endl ; 716 G4cout << "distance = " << distance[0] << G 722 G4cout << "distance = " << distance[0] << G4endl ; 717 G4cout << "X = " << xx << G4endl ; 723 G4cout << "X = " << xx << G4endl ; 718 #endif 724 #endif 719 725 720 G4bool isvalid = true; 726 G4bool isvalid = true; 721 gxx[0] = ComputeGlobalPoint(xx); 727 gxx[0] = ComputeGlobalPoint(xx); 722 728 723 #ifdef G4TWISTDEBUG 729 #ifdef G4TWISTDEBUG 724 G4cout << "intersection Point found: " << g 730 G4cout << "intersection Point found: " << gxx[0] << G4endl ; 725 G4cout << "distance = " << distance[0] << G 731 G4cout << "distance = " << distance[0] << G4endl ; 726 #endif 732 #endif 727 733 728 fCurStat.SetCurrentStatus(0, gxx[0], distan 734 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0], 729 isvalid, 1, kDontV 735 isvalid, 1, kDontValidate, &gp); 730 return 1; 736 return 1; 731 } 737 } 732 738 733 //============================================ 739 //===================================================================== 734 //* GetAreaCode ------------------------------ 740 //* GetAreaCode ------------------------------------------------------- 735 741 736 G4int G4TwistTrapParallelSide::GetAreaCode(con 742 G4int G4TwistTrapParallelSide::GetAreaCode(const G4ThreeVector& xx, 737 G4b 743 G4bool withTol) 738 { 744 { 739 // We must use the function in local coordi 745 // We must use the function in local coordinate system. 740 // See the description of DistanceToSurface 746 // See the description of DistanceToSurface(p,v). 741 747 742 const G4double ctol = 0.5 * kCarTolerance; 748 const G4double ctol = 0.5 * kCarTolerance; 743 749 744 G4double phi ; 750 G4double phi ; 745 G4double yprime ; 751 G4double yprime ; 746 GetPhiUAtX(xx, phi,yprime ) ; 752 GetPhiUAtX(xx, phi,yprime ) ; 747 753 748 G4double fXAxisMax = GetBoundaryMax(phi) ; 754 G4double fXAxisMax = GetBoundaryMax(phi) ; 749 G4double fXAxisMin = GetBoundaryMin(phi) ; 755 G4double fXAxisMin = GetBoundaryMin(phi) ; 750 756 751 #ifdef G4TWISTDEBUG 757 #ifdef G4TWISTDEBUG 752 G4cout << "GetAreaCode: phi = " << phi << G 758 G4cout << "GetAreaCode: phi = " << phi << G4endl ; 753 G4cout << "GetAreaCode: yprime = " << yprim 759 G4cout << "GetAreaCode: yprime = " << yprime << G4endl ; 754 G4cout << "Intervall is " << fXAxisMin << " 760 G4cout << "Intervall is " << fXAxisMin << " to " << fXAxisMax << G4endl ; 755 #endif 761 #endif 756 762 757 G4int areacode = sInside; 763 G4int areacode = sInside; 758 764 759 if (fAxis[0] == kXAxis && fAxis[1] == kZAxi 765 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) 760 { 766 { 761 G4int zaxis = 1; 767 G4int zaxis = 1; 762 768 763 if (withTol) 769 if (withTol) 764 { 770 { 765 G4bool isoutside = false; 771 G4bool isoutside = false; 766 772 767 // test boundary of xaxis 773 // test boundary of xaxis 768 774 769 if (yprime < fXAxisMin + ctol) 775 if (yprime < fXAxisMin + ctol) 770 { 776 { 771 areacode |= (sAxis0 & (sAxisX | sA 777 areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary; 772 if (yprime <= fXAxisMin - ctol) is 778 if (yprime <= fXAxisMin - ctol) isoutside = true; 773 779 774 } 780 } 775 else if (yprime > fXAxisMax - ctol) 781 else if (yprime > fXAxisMax - ctol) 776 { 782 { 777 areacode |= (sAxis0 & (sAxisX | sA 783 areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary; 778 if (yprime >= fXAxisMax + ctol) i 784 if (yprime >= fXAxisMax + ctol) isoutside = true; 779 } 785 } 780 786 781 // test boundary of z-axis 787 // test boundary of z-axis 782 788 783 if (xx.z() < fAxisMin[zaxis] + ctol) 789 if (xx.z() < fAxisMin[zaxis] + ctol) 784 { 790 { 785 areacode |= (sAxis1 & (sAxisZ | sA 791 areacode |= (sAxis1 & (sAxisZ | sAxisMin)); 786 792 787 if ((areacode & sBoundary) != 0) << 793 if (areacode & sBoundary) areacode |= sCorner; // xx is on corner. 788 else areaco 794 else areacode |= sBoundary; 789 if (xx.z() <= fAxisMin[zaxis] - ct 795 if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true; 790 796 791 } 797 } 792 else if (xx.z() > fAxisMax[zaxis] - c 798 else if (xx.z() > fAxisMax[zaxis] - ctol) 793 { 799 { 794 areacode |= (sAxis1 & (sAxisZ | sA 800 areacode |= (sAxis1 & (sAxisZ | sAxisMax)); 795 801 796 if ((areacode & sBoundary) != 0) << 802 if (areacode & sBoundary) areacode |= sCorner; // xx is on corner. 797 else areaco 803 else areacode |= sBoundary; 798 if (xx.z() >= fAxisMax[zaxis] + ct 804 if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true; 799 } 805 } 800 806 801 // if isoutside = true, clear inside 807 // if isoutside = true, clear inside bit. 802 // if not on boundary, add axis infor 808 // if not on boundary, add axis information. 803 809 804 if (isoutside) 810 if (isoutside) 805 { 811 { 806 G4int tmpareacode = areacode & (~s 812 G4int tmpareacode = areacode & (~sInside); 807 areacode = tmpareacode; 813 areacode = tmpareacode; 808 } 814 } 809 else if ((areacode & sBoundary) != sB 815 else if ((areacode & sBoundary) != sBoundary) 810 { 816 { 811 areacode |= (sAxis0 & sAxisX) | (s 817 areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ); 812 } 818 } 813 819 814 } 820 } 815 else 821 else 816 { 822 { 817 // boundary of y-axis 823 // boundary of y-axis 818 824 819 if (yprime < fXAxisMin ) 825 if (yprime < fXAxisMin ) 820 { 826 { 821 areacode |= (sAxis0 & (sAxisX | sA 827 areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary; 822 } 828 } 823 else if (yprime > fXAxisMax) 829 else if (yprime > fXAxisMax) 824 { 830 { 825 areacode |= (sAxis0 & (sAxisX | sA 831 areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary; 826 } 832 } 827 833 828 // boundary of z-axis 834 // boundary of z-axis 829 835 830 if (xx.z() < fAxisMin[zaxis]) 836 if (xx.z() < fAxisMin[zaxis]) 831 { 837 { 832 areacode |= (sAxis1 & (sAxisZ | sA 838 areacode |= (sAxis1 & (sAxisZ | sAxisMin)); 833 if ((areacode & sBoundary) != 0) << 839 if (areacode & sBoundary) areacode |= sCorner; // xx is on corner. 834 else areaco 840 else areacode |= sBoundary; 835 841 836 } 842 } 837 else if (xx.z() > fAxisMax[zaxis]) 843 else if (xx.z() > fAxisMax[zaxis]) 838 { 844 { 839 areacode |= (sAxis1 & (sAxisZ | sA 845 areacode |= (sAxis1 & (sAxisZ | sAxisMax)) ; 840 if ((areacode & sBoundary) != 0) << 846 if (areacode & sBoundary) areacode |= sCorner; // xx is on corner. 841 else areaco 847 else areacode |= sBoundary; 842 } 848 } 843 849 844 if ((areacode & sBoundary) != sBounda 850 if ((areacode & sBoundary) != sBoundary) 845 { 851 { 846 areacode |= (sAxis0 & sAxisX) | (s 852 areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ); 847 } 853 } 848 } 854 } 849 return areacode; 855 return areacode; 850 } 856 } 851 else 857 else 852 { 858 { 853 G4Exception("G4TwistTrapParallelSide::Ge 859 G4Exception("G4TwistTrapParallelSide::GetAreaCode()", 854 "GeomSolids0001", FatalExcep 860 "GeomSolids0001", FatalException, 855 "Feature NOT implemented !") 861 "Feature NOT implemented !"); 856 } 862 } 857 return areacode; 863 return areacode; 858 } 864 } 859 865 860 //============================================ 866 //===================================================================== 861 //* SetCorners() ----------------------------- 867 //* SetCorners() ------------------------------------------------------ 862 868 863 void G4TwistTrapParallelSide::SetCorners() 869 void G4TwistTrapParallelSide::SetCorners() 864 { 870 { 865 871 866 // Set Corner points in local coodinate. 872 // Set Corner points in local coodinate. 867 873 868 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis 874 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) 869 { 875 { 870 G4double x, y, z; 876 G4double x, y, z; 871 877 872 // corner of Axis0min and Axis1min 878 // corner of Axis0min and Axis1min 873 879 874 x = -fdeltaX/2. + (-fDx2 + fDy1*fTAlph)*st 880 x = -fdeltaX/2. + (-fDx2 + fDy1*fTAlph)*std::cos(fPhiTwist/2.) 875 + fDy1*std::sin(fPhiTwist/2.) ; 881 + fDy1*std::sin(fPhiTwist/2.) ; 876 y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/ 882 y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/2.) 877 + (fDx2 - fDy1*fTAlph)*std::sin(fPhiTwis 883 + (fDx2 - fDy1*fTAlph)*std::sin(fPhiTwist/2.) ; 878 z = -fDz ; 884 z = -fDz ; 879 885 880 SetCorner(sC0Min1Min, x, y, z); 886 SetCorner(sC0Min1Min, x, y, z); 881 887 882 // corner of Axis0max and Axis1min 888 // corner of Axis0max and Axis1min 883 889 884 x = -fdeltaX/2. + (fDx2 + fDy1*fTAlph)*std 890 x = -fdeltaX/2. + (fDx2 + fDy1*fTAlph)*std::cos(fPhiTwist/2.) 885 + fDy1*std::sin(fPhiTwist/2.) ; 891 + fDy1*std::sin(fPhiTwist/2.) ; 886 y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/ 892 y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/2.) 887 - (fDx2 + fDy1*fTAlph)*std::sin(fPhiTwis 893 - (fDx2 + fDy1*fTAlph)*std::sin(fPhiTwist/2.) ; 888 z = -fDz; 894 z = -fDz; 889 895 890 SetCorner(sC0Max1Min, x, y, z); 896 SetCorner(sC0Max1Min, x, y, z); 891 897 892 // corner of Axis0max and Axis1max 898 // corner of Axis0max and Axis1max 893 x = fdeltaX/2. + (fDx4 + fDy2*fTAlph)*std: 899 x = fdeltaX/2. + (fDx4 + fDy2*fTAlph)*std::cos(fPhiTwist/2.) 894 - fDy2*std::sin(fPhiTwist/2.) ; 900 - fDy2*std::sin(fPhiTwist/2.) ; 895 y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2 901 y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2.) 896 + (fDx4 + fDy2*fTAlph)*std::sin(fPhiTwis 902 + (fDx4 + fDy2*fTAlph)*std::sin(fPhiTwist/2.) ; 897 z = fDz ; 903 z = fDz ; 898 904 899 SetCorner(sC0Max1Max, x, y, z); 905 SetCorner(sC0Max1Max, x, y, z); 900 906 901 // corner of Axis0min and Axis1max 907 // corner of Axis0min and Axis1max 902 x = fdeltaX/2. + (-fDx4 + fDy2*fTAlph)*std 908 x = fdeltaX/2. + (-fDx4 + fDy2*fTAlph)*std::cos(fPhiTwist/2.) 903 - fDy2*std::sin(fPhiTwist/2.) ; 909 - fDy2*std::sin(fPhiTwist/2.) ; 904 y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2 910 y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2.) 905 + (-fDx4 + fDy2*fTAlph)*std::sin(fPhiTwi 911 + (-fDx4 + fDy2*fTAlph)*std::sin(fPhiTwist/2.) ; 906 z = fDz ; 912 z = fDz ; 907 913 908 SetCorner(sC0Min1Max, x, y, z); 914 SetCorner(sC0Min1Max, x, y, z); 909 } 915 } 910 else 916 else 911 { 917 { 912 G4Exception("G4TwistTrapParallelSide::SetC 918 G4Exception("G4TwistTrapParallelSide::SetCorners()", 913 "GeomSolids0001", FatalExcepti 919 "GeomSolids0001", FatalException, 914 "Method NOT implemented !"); 920 "Method NOT implemented !"); 915 } 921 } 916 } 922 } 917 923 918 //============================================ 924 //===================================================================== 919 //* SetBoundaries() -------------------------- 925 //* SetBoundaries() --------------------------------------------------- 920 926 921 void G4TwistTrapParallelSide::SetBoundaries() 927 void G4TwistTrapParallelSide::SetBoundaries() 922 { 928 { 923 // Set direction-unit vector of boundary-li 929 // Set direction-unit vector of boundary-lines in local coodinate. 924 // 930 // 925 931 926 G4ThreeVector direction; 932 G4ThreeVector direction; 927 933 928 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis 934 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) 929 { 935 { 930 // sAxis0 & sAxisMin 936 // sAxis0 & sAxisMin 931 direction = GetCorner(sC0Min1Max) - GetCor 937 direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min); 932 direction = direction.unit(); 938 direction = direction.unit(); 933 SetBoundary(sAxis0 & (sAxisX | sAxisMin), 939 SetBoundary(sAxis0 & (sAxisX | sAxisMin), direction, 934 GetCorner(sC0Min1Min), sAxisZ) 940 GetCorner(sC0Min1Min), sAxisZ) ; 935 941 936 // sAxis0 & sAxisMax 942 // sAxis0 & sAxisMax 937 direction = GetCorner(sC0Max1Max) - GetCor 943 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min); 938 direction = direction.unit(); 944 direction = direction.unit(); 939 SetBoundary(sAxis0 & (sAxisX | sAxisMax), 945 SetBoundary(sAxis0 & (sAxisX | sAxisMax), direction, 940 GetCorner(sC0Max1Min), sAxisZ) 946 GetCorner(sC0Max1Min), sAxisZ); 941 947 942 // sAxis1 & sAxisMin 948 // sAxis1 & sAxisMin 943 direction = GetCorner(sC0Max1Min) - GetCor 949 direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min); 944 direction = direction.unit(); 950 direction = direction.unit(); 945 SetBoundary(sAxis1 & (sAxisZ | sAxisMin), 951 SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction, 946 GetCorner(sC0Min1Min), sAxisX) 952 GetCorner(sC0Min1Min), sAxisX); 947 953 948 // sAxis1 & sAxisMax 954 // sAxis1 & sAxisMax 949 direction = GetCorner(sC0Max1Max) - GetCor 955 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max); 950 direction = direction.unit(); 956 direction = direction.unit(); 951 SetBoundary(sAxis1 & (sAxisZ | sAxisMax), 957 SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction, 952 GetCorner(sC0Min1Max), sAxisX) 958 GetCorner(sC0Min1Max), sAxisX); 953 } 959 } 954 else 960 else 955 { 961 { 956 G4Exception("G4TwistTrapParallelSide::SetC 962 G4Exception("G4TwistTrapParallelSide::SetCorners()", 957 "GeomSolids0001", FatalExcepti 963 "GeomSolids0001", FatalException, 958 "Feature NOT implemented !"); 964 "Feature NOT implemented !"); 959 } 965 } 960 } 966 } 961 967 962 //============================================ 968 //===================================================================== 963 //* GetPhiUAtX() ----------------------------- 969 //* GetPhiUAtX() ------------------------------------------------------ 964 970 965 void 971 void 966 G4TwistTrapParallelSide::GetPhiUAtX( const G4T << 972 G4TwistTrapParallelSide::GetPhiUAtX( G4ThreeVector p, 967 G4double& 973 G4double& phi, G4double& u ) 968 { 974 { 969 // find closest point XX on surface for a gi 975 // find closest point XX on surface for a given point p 970 // X0 is a point on the surface, d is the d 976 // X0 is a point on the surface, d is the direction 971 // ( both for a fixed z = pz) 977 // ( both for a fixed z = pz) 972 978 973 // phi is given by the z coordinate of p 979 // phi is given by the z coordinate of p 974 980 975 phi = p.z()/(2*fDz)*fPhiTwist ; 981 phi = p.z()/(2*fDz)*fPhiTwist ; 976 982 977 u = ((-(fdeltaX*phi) + fPhiTwist*p.x())* std 983 u = ((-(fdeltaX*phi) + fPhiTwist*p.x())* std::cos(phi) 978 + (-(fdeltaY*phi) + fPhiTwist*p.y())*std:: 984 + (-(fdeltaY*phi) + fPhiTwist*p.y())*std::sin(phi))/fPhiTwist ; 979 } 985 } 980 986 981 //============================================ 987 //===================================================================== 982 //* ProjectPoint() --------------------------- 988 //* ProjectPoint() ---------------------------------------------------- 983 989 984 G4ThreeVector G4TwistTrapParallelSide::Project 990 G4ThreeVector G4TwistTrapParallelSide::ProjectPoint(const G4ThreeVector& p, 985 991 G4bool isglobal) 986 { 992 { 987 // Get Rho at p.z() on Hyperbolic Surface. 993 // Get Rho at p.z() on Hyperbolic Surface. 988 G4ThreeVector tmpp; 994 G4ThreeVector tmpp; 989 if (isglobal) 995 if (isglobal) 990 { 996 { 991 tmpp = fRot.inverse()*p - fTrans; 997 tmpp = fRot.inverse()*p - fTrans; 992 } 998 } 993 else 999 else 994 { 1000 { 995 tmpp = p; 1001 tmpp = p; 996 } 1002 } 997 1003 998 G4double phi ; 1004 G4double phi ; 999 G4double u ; 1005 G4double u ; 1000 1006 1001 GetPhiUAtX( tmpp, phi, u ) ; // calculate 1007 GetPhiUAtX( tmpp, phi, u ) ; // calculate (phi, u) for p close to surface 1002 1008 1003 G4ThreeVector xx = SurfacePoint(phi,u) ; // 1009 G4ThreeVector xx = SurfacePoint(phi,u) ; // transform back to Cartesian coords 1004 1010 1005 if (isglobal) 1011 if (isglobal) 1006 { 1012 { 1007 return (fRot * xx + fTrans); 1013 return (fRot * xx + fTrans); 1008 } 1014 } 1009 else 1015 else 1010 { 1016 { 1011 return xx; 1017 return xx; 1012 } 1018 } 1013 } 1019 } 1014 1020 1015 //=========================================== 1021 //===================================================================== 1016 //* GetFacets() ----------------------------- 1022 //* GetFacets() ------------------------------------------------------- 1017 1023 1018 void G4TwistTrapParallelSide::GetFacets( G4in 1024 void G4TwistTrapParallelSide::GetFacets( G4int k, G4int n, G4double xyz[][3], 1019 G4in 1025 G4int faces[][4], G4int iside ) 1020 { 1026 { 1021 G4double phi ; 1027 G4double phi ; 1022 G4double z, u ; // the two parameters f 1028 G4double z, u ; // the two parameters for the surface equation 1023 G4ThreeVector p ; // a point on the surfac 1029 G4ThreeVector p ; // a point on the surface, given by (z,u) 1024 1030 1025 G4int nnode ; 1031 G4int nnode ; 1026 G4int nface ; 1032 G4int nface ; 1027 1033 1028 G4double umin, umax ; 1034 G4double umin, umax ; 1029 1035 1030 // calculate the (n-1)*(k-1) vertices 1036 // calculate the (n-1)*(k-1) vertices 1031 1037 1032 for ( G4int i = 0 ; i<n ; ++i ) 1038 for ( G4int i = 0 ; i<n ; ++i ) 1033 { 1039 { 1034 z = -fDz+i*(2.*fDz)/(n-1) ; 1040 z = -fDz+i*(2.*fDz)/(n-1) ; 1035 phi = z*fPhiTwist/(2*fDz) ; 1041 phi = z*fPhiTwist/(2*fDz) ; 1036 umin = GetBoundaryMin(phi) ; 1042 umin = GetBoundaryMin(phi) ; 1037 umax = GetBoundaryMax(phi) ; 1043 umax = GetBoundaryMax(phi) ; 1038 1044 1039 for ( G4int j = 0 ; j<k ; ++j ) 1045 for ( G4int j = 0 ; j<k ; ++j ) 1040 { 1046 { 1041 nnode = GetNode(i,j,k,n,iside) ; 1047 nnode = GetNode(i,j,k,n,iside) ; 1042 u = umax - j*(umax-umin)/(k-1) ; 1048 u = umax - j*(umax-umin)/(k-1) ; 1043 p = SurfacePoint(phi,u,true) ; // surf 1049 p = SurfacePoint(phi,u,true) ; // surface point in global coords 1044 1050 1045 xyz[nnode][0] = p.x() ; 1051 xyz[nnode][0] = p.x() ; 1046 xyz[nnode][1] = p.y() ; 1052 xyz[nnode][1] = p.y() ; 1047 xyz[nnode][2] = p.z() ; 1053 xyz[nnode][2] = p.z() ; 1048 1054 1049 if ( i<n-1 && j<k-1 ) // conterclock 1055 if ( i<n-1 && j<k-1 ) // conterclock wise filling 1050 { 1056 { 1051 nface = GetFace(i,j,k,n,iside) ; 1057 nface = GetFace(i,j,k,n,iside) ; 1052 faces[nface][0] = GetEdgeVisibility(i 1058 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1) 1053 * (GetNode(i ,j ,k, 1059 * (GetNode(i ,j ,k,n,iside)+1) ; // fortran numbering 1054 faces[nface][1] = GetEdgeVisibility(i 1060 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1) 1055 * (GetNode(i ,j+1,k, 1061 * (GetNode(i ,j+1,k,n,iside)+1) ; 1056 faces[nface][2] = GetEdgeVisibility(i 1062 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1) 1057 * (GetNode(i+1,j+1,k, 1063 * (GetNode(i+1,j+1,k,n,iside)+1) ; 1058 faces[nface][3] = GetEdgeVisibility(i 1064 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1) 1059 * (GetNode(i+1,j ,k, 1065 * (GetNode(i+1,j ,k,n,iside)+1) ; 1060 } 1066 } 1061 } 1067 } 1062 } 1068 } 1063 } 1069 } 1064 1070