| 1 | c-----------------------------------------------------------------------
|
|---|
| 2 | subroutine arcsrf(xml,yml,zml,nxl,nyl,nzl,ie,isid)
|
|---|
| 3 | include 'SIZE'
|
|---|
| 4 | include 'GEOM'
|
|---|
| 5 | include 'INPUT'
|
|---|
| 6 | include 'TOPOL'
|
|---|
| 7 | include 'WZ'
|
|---|
| 8 | C
|
|---|
| 9 | C ....note..... CTMP1 is used in this format in several subsequent routines
|
|---|
| 10 | C
|
|---|
| 11 | COMMON /CTMP1/ H(LX1,3,2),XCRVED(LX1),YCRVED(LY1),ZCRVED(LZ1)
|
|---|
| 12 | $ , ZGML(LX1,3),WORK(3,LX1,LZ1)
|
|---|
| 13 | DIMENSION XML(NXL,NYL,NZL,1),YML(NXL,NYL,NZL,1),ZML(NXL,NYL,NZL,1)
|
|---|
| 14 | LOGICAL IFGLJ
|
|---|
| 15 | C
|
|---|
| 16 | IFGLJ = .FALSE.
|
|---|
| 17 | IF (IFAXIS .AND. IFRZER(IE) .AND. (ISID.EQ.2 .OR. ISID.EQ.4))
|
|---|
| 18 | $IFGLJ = .TRUE.
|
|---|
| 19 | C
|
|---|
| 20 | PT1X = XC(ISID,IE)
|
|---|
| 21 | PT1Y = YC(ISID,IE)
|
|---|
| 22 | IF(ISID.EQ.4) THEN
|
|---|
| 23 | PT2X = XC(1,IE)
|
|---|
| 24 | PT2Y = YC(1,IE)
|
|---|
| 25 | ELSE IF(ISID.EQ.8) THEN
|
|---|
| 26 | PT2X = XC(5,IE)
|
|---|
| 27 | PT2Y = YC(5,IE)
|
|---|
| 28 | ELSE
|
|---|
| 29 | PT2X = XC(ISID+1,IE)
|
|---|
| 30 | PT2Y = YC(ISID+1,IE)
|
|---|
| 31 | ENDIF
|
|---|
| 32 | C
|
|---|
| 33 | C Find slope of perpendicular
|
|---|
| 34 | RADIUS=CURVE(1,ISID,IE)
|
|---|
| 35 | GAP=SQRT( (PT1X-PT2X)**2 + (PT1Y-PT2Y)**2 )
|
|---|
| 36 | IF (ABS(2.0*RADIUS).LE.GAP*1.00001) THEN
|
|---|
| 37 | write(6,10) RADIUS,ISID,IE,GAP
|
|---|
| 38 | 10 FORMAT(//,2X,'ERROR: Too small a radius (',G11.3
|
|---|
| 39 | $ ,') specified for side',I2,' of element',I4,': '
|
|---|
| 40 | $ ,G11.3,/,2X,'ABORTING during mesh generation.')
|
|---|
| 41 | call exitt
|
|---|
| 42 | ENDIF
|
|---|
| 43 | XS = PT2Y-PT1Y
|
|---|
| 44 | YS = PT1X-PT2X
|
|---|
| 45 | C Make length Radius
|
|---|
| 46 | XYS=SQRT(XS**2+YS**2)
|
|---|
| 47 | C Find Center
|
|---|
| 48 | DTHETA = ABS(ASIN(0.5*GAP/RADIUS))
|
|---|
| 49 | PT12X = (PT1X + PT2X)/2.0
|
|---|
| 50 | PT12Y = (PT1Y + PT2Y)/2.0
|
|---|
| 51 | XCENN = PT12X - XS/XYS * RADIUS*COS(DTHETA)
|
|---|
| 52 | YCENN = PT12Y - YS/XYS * RADIUS*COS(DTHETA)
|
|---|
| 53 | THETA0 = ATAN2((PT12Y-YCENN),(PT12X-XCENN))
|
|---|
| 54 | IF (IFGLJ) THEN
|
|---|
| 55 | FAC = SIGN(1.0,RADIUS)
|
|---|
| 56 | THETA1 = THETA0 - FAC*DTHETA
|
|---|
| 57 | THETA2 = THETA0 + FAC*DTHETA
|
|---|
| 58 | ENDIF
|
|---|
| 59 | C Compute perturbation of geometry
|
|---|
| 60 | ISID1 = MOD1(ISID,4)
|
|---|
| 61 | IF (IFGLJ) THEN
|
|---|
| 62 | I1 = ISID/2
|
|---|
| 63 | I2 = 2 - ISID/4
|
|---|
| 64 | DO 15 IY=1,NYL
|
|---|
| 65 | ANG = H(IY,2,I1)*THETA1 + H(IY,2,I2)*THETA2
|
|---|
| 66 | XCRVED(IY)=XCENN + ABS(RADIUS)*COS(ANG)
|
|---|
| 67 | $ - (H(IY,2,I1)*PT1X + H(IY,2,I2)*PT2X)
|
|---|
| 68 | YCRVED(IY)=YCENN + ABS(RADIUS) * SIN(ANG)
|
|---|
| 69 | $ - (H(IY,2,I1)*PT1Y + H(IY,2,I2)*PT2Y)
|
|---|
| 70 | 15 CONTINUE
|
|---|
| 71 | ELSE
|
|---|
| 72 | DO 20 IX=1,NXL
|
|---|
| 73 | IXT=IX
|
|---|
| 74 | IF (ISID1.GT.2) IXT=NXL+1-IX
|
|---|
| 75 | R=ZGML(IX,1)
|
|---|
| 76 | IF (RADIUS.LT.0.0) R=-R
|
|---|
| 77 | XCRVED(IXT) = XCENN + ABS(RADIUS) * COS(THETA0 + R*DTHETA)
|
|---|
| 78 | $ - ( H(IX,1,1)*PT1X + H(IX,1,2)*PT2X )
|
|---|
| 79 | YCRVED(IXT) = YCENN + ABS(RADIUS) * SIN(THETA0 + R*DTHETA)
|
|---|
| 80 | $ - ( H(IX,1,1)*PT1Y + H(IX,1,2)*PT2Y )
|
|---|
| 81 | 20 CONTINUE
|
|---|
| 82 | ENDIF
|
|---|
| 83 | C Points all set, add perturbation to current mesh.
|
|---|
| 84 | ISID1 = MOD1(ISID,4)
|
|---|
| 85 | ISID1 = EFACE1(ISID1)
|
|---|
| 86 | IZT = (ISID-1)/4+1
|
|---|
| 87 | IYT = ISID1-2
|
|---|
| 88 | IXT = ISID1
|
|---|
| 89 | IF (ISID1.LE.2) THEN
|
|---|
| 90 | CALL ADDTNSR(XML(1,1,1,IE),H(1,1,IXT),XCRVED,H(1,3,IZT)
|
|---|
| 91 | $ ,NXL,NYL,NZL)
|
|---|
| 92 | CALL ADDTNSR(YML(1,1,1,IE),H(1,1,IXT),YCRVED,H(1,3,IZT)
|
|---|
| 93 | $ ,NXL,NYL,NZL)
|
|---|
| 94 | ELSE
|
|---|
| 95 | CALL ADDTNSR(XML(1,1,1,IE),XCRVED,H(1,2,IYT),H(1,3,IZT)
|
|---|
| 96 | $ ,NXL,NYL,NZL)
|
|---|
| 97 | CALL ADDTNSR(YML(1,1,1,IE),YCRVED,H(1,2,IYT),H(1,3,IZT)
|
|---|
| 98 | $ ,NXL,NYL,NZL)
|
|---|
| 99 | ENDIF
|
|---|
| 100 | return
|
|---|
| 101 | end
|
|---|
| 102 | c-----------------------------------------------------------------------
|
|---|
| 103 | subroutine defsrf(xml,yml,zml,nxl,nyl,nzl,ie,iface1,ccv)
|
|---|
| 104 | include 'SIZE'
|
|---|
| 105 | include 'TOPOL'
|
|---|
| 106 | include 'GEOM'
|
|---|
| 107 | include 'WZ'
|
|---|
| 108 | COMMON /CTMP1/ H(LX1,3,2),XCRVED(LX1),YCRVED(LY1),ZCRVED(LZ1)
|
|---|
| 109 | $ , ZGML(LX1,3),WORK(3,LX1,LZ1)
|
|---|
| 110 | C
|
|---|
| 111 | DIMENSION XML(NXL,NYL,NZL,1),YML(NXL,NYL,NZL,1),ZML(NXL,NYL,NZL,1)
|
|---|
| 112 | DIMENSION X1(3),X2(3),X3(3),DX(3)
|
|---|
| 113 | DIMENSION IOPP(3),NXX(3)
|
|---|
| 114 | CHARACTER*1 CCV
|
|---|
| 115 | C
|
|---|
| 116 | CALL DSSET(NXL,NYL,NZL)
|
|---|
| 117 | IFACE = EFACE1(IFACE1)
|
|---|
| 118 | JS1 = SKPDAT(1,IFACE)
|
|---|
| 119 | JF1 = SKPDAT(2,IFACE)
|
|---|
| 120 | JSKIP1 = SKPDAT(3,IFACE)
|
|---|
| 121 | JS2 = SKPDAT(4,IFACE)
|
|---|
| 122 | JF2 = SKPDAT(5,IFACE)
|
|---|
| 123 | JSKIP2 = SKPDAT(6,IFACE)
|
|---|
| 124 | C
|
|---|
| 125 | IOPP(1) = NXL-1
|
|---|
| 126 | IOPP(2) = NXL*(NYL-1)
|
|---|
| 127 | IOPP(3) = NXL*NYL*(NZL-1)
|
|---|
| 128 | NXX(1) = NXL
|
|---|
| 129 | NXX(2) = NYL
|
|---|
| 130 | NXX(3) = NZL
|
|---|
| 131 | IDIR = 2*MOD(IFACE,2) - 1
|
|---|
| 132 | IFC2 = (IFACE+1)/2
|
|---|
| 133 | DELT = 0.0
|
|---|
| 134 | C DELT = SIDE(4,IFACE,IE)
|
|---|
| 135 | C
|
|---|
| 136 | C Compute surface deflection and perturbation due to face IFACE
|
|---|
| 137 | C
|
|---|
| 138 | DO 200 J2=JS2,JF2,JSKIP2
|
|---|
| 139 | DO 200 J1=JS1,JF1,JSKIP1
|
|---|
| 140 | JOPP = J1 + IOPP(IFC2)*IDIR
|
|---|
| 141 | X2(1) = XML(J1,J2,1,IE)
|
|---|
| 142 | X2(2) = YML(J1,J2,1,IE)
|
|---|
| 143 | X2(3) = ZML(J1,J2,1,IE)
|
|---|
| 144 | X1(1) = XML(JOPP,J2,1,IE)
|
|---|
| 145 | X1(2) = YML(JOPP,J2,1,IE)
|
|---|
| 146 | X1(3) = ZML(JOPP,J2,1,IE)
|
|---|
| 147 | CALL INTRSC(X3,X2,X1,DELT,IE,IFACE1)
|
|---|
| 148 | C
|
|---|
| 149 | DX(1) = X3(1)-X2(1)
|
|---|
| 150 | DX(2) = X3(2)-X2(2)
|
|---|
| 151 | DX(3) = X3(3)-X2(3)
|
|---|
| 152 | C
|
|---|
| 153 | NXS = NXX(IFC2)
|
|---|
| 154 | JOFF = (J1-JOPP)/(NXS-1)
|
|---|
| 155 | DO 100 IX = 2,NXS
|
|---|
| 156 | J = JOPP + JOFF*(IX-1)
|
|---|
| 157 | ZETA = 0.5*(ZGML(IX,IFC2) + 1.0)
|
|---|
| 158 | XML(J,J2,1,IE) = XML(J,J2,1,IE)+DX(1)*ZETA
|
|---|
| 159 | YML(J,J2,1,IE) = YML(J,J2,1,IE)+DX(2)*ZETA
|
|---|
| 160 | ZML(J,J2,1,IE) = ZML(J,J2,1,IE)+DX(3)*ZETA
|
|---|
| 161 | 100 CONTINUE
|
|---|
| 162 | 200 CONTINUE
|
|---|
| 163 | C
|
|---|
| 164 | return
|
|---|
| 165 | end
|
|---|
| 166 | c-----------------------------------------------------------------------
|
|---|
| 167 | subroutine intrsc(x3,x2,x1,delt,ie,iface)
|
|---|
| 168 | C
|
|---|
| 169 | DIMENSION X1(3),X2(3),X3(3)
|
|---|
| 170 | COMMON /SRFCEI/ IEL,IFCE
|
|---|
| 171 | COMMON /SRFCER/ X0(3),DX(3)
|
|---|
| 172 | COMMON /SRFCEL/ SUCCES
|
|---|
| 173 | LOGICAL SUCCES
|
|---|
| 174 | COMMON /TOLRNC/ TOL
|
|---|
| 175 | C
|
|---|
| 176 | C Load parameters for surface function FNC
|
|---|
| 177 | C
|
|---|
| 178 | IEL = IE
|
|---|
| 179 | IFCE = IFACE
|
|---|
| 180 | X0(1) = X1(1)
|
|---|
| 181 | X0(2) = X1(2)
|
|---|
| 182 | X0(3) = X1(3)
|
|---|
| 183 | DX(1) = X2(1)-X1(1)
|
|---|
| 184 | DX(2) = X2(2)-X1(2)
|
|---|
| 185 | DX(3) = X2(3)-X1(3)
|
|---|
| 186 | DIST = SQRT ( DX(1)**2 + DX(2)**2 + DX(3)**2 )
|
|---|
| 187 | C
|
|---|
| 188 | C Initial guess for bracket is given by size of element face (DELT).
|
|---|
| 189 | C
|
|---|
| 190 | ETA2 = 1.0
|
|---|
| 191 | ETA1 = ETA2 - DELT/DIST
|
|---|
| 192 | ETA1 = MAX(ETA1,0.0)
|
|---|
| 193 | CALL ZBRAC(ETA1,ETA2,SUCCES)
|
|---|
| 194 | C
|
|---|
| 195 | TOL = 1.0E-5
|
|---|
| 196 | TOLSRF = TOL*(ETA2-ETA1)
|
|---|
| 197 | IF (SUCCES) ETA3 = ZBRENT(ETA1,ETA2,TOLSRF)
|
|---|
| 198 | C
|
|---|
| 199 | X3(1) = X0(1) + DX(1)*ETA3
|
|---|
| 200 | X3(2) = X0(2) + DX(2)*ETA3
|
|---|
| 201 | X3(3) = X0(3) + DX(3)*ETA3
|
|---|
| 202 | C
|
|---|
| 203 | return
|
|---|
| 204 | end
|
|---|
| 205 | c-----------------------------------------------------------------------
|
|---|
| 206 | subroutine zbrac(x1,x2,succes)
|
|---|
| 207 | C
|
|---|
| 208 | C Given a function FNC and an initial guess (X1,X2), the routine
|
|---|
| 209 | C expands the range geometrically until a root is bracketed by the
|
|---|
| 210 | C returned range (X1,X2) [SUCCES=.TRUE.], or until the range becomes
|
|---|
| 211 | C unacceptably large [SUCCES=.FALSE.].
|
|---|
| 212 | C ( Numerical Recipes, p. 245; pff 9 Aug 1989 09:00:20 )
|
|---|
| 213 | C
|
|---|
| 214 | PARAMETER (FACTOR=1.08,NTRY=50)
|
|---|
| 215 | LOGICAL SUCCES
|
|---|
| 216 | C
|
|---|
| 217 | SUCCES = .TRUE.
|
|---|
| 218 | C
|
|---|
| 219 | IF (X1.EQ.X2) X1 = .99*X1
|
|---|
| 220 | IF (X1.EQ.0.0) X1 = 1.0E-04
|
|---|
| 221 | C
|
|---|
| 222 | F1 = FNC(X1)
|
|---|
| 223 | F2 = FNC(X2)
|
|---|
| 224 | DO 100 J=1,NTRY
|
|---|
| 225 | IF (F1*F2.LT.0.0) return
|
|---|
| 226 | IF (ABS(F1).LT.ABS(F2)) THEN
|
|---|
| 227 | X1 = X1 + FACTOR*(X1-X2)
|
|---|
| 228 | F1 = FNC(X1)
|
|---|
| 229 | ELSE
|
|---|
| 230 | X2 = X2 + FACTOR*(X2-X1)
|
|---|
| 231 | F2 = FNC(X2)
|
|---|
| 232 | ENDIF
|
|---|
| 233 | 100 CONTINUE
|
|---|
| 234 | SUCCES = .FALSE.
|
|---|
| 235 | return
|
|---|
| 236 | end
|
|---|
| 237 | FUNCTION ZBRENT(X1,X2,TOL)
|
|---|
| 238 | C
|
|---|
| 239 | C Using the Van Wijngaarden-Dekker-Brent Method, find the root
|
|---|
| 240 | C of a function FNC known to lie between X1 and X2. The root,
|
|---|
| 241 | C returned as ZBRENT, will be refined until its accuracy is TOL.
|
|---|
| 242 | C
|
|---|
| 243 | PARAMETER (ITMAX=100,EPS=3.0E-8)
|
|---|
| 244 | C
|
|---|
| 245 | A = X1
|
|---|
| 246 | B = X2
|
|---|
| 247 | FA = FNC(A)
|
|---|
| 248 | FB = FNC(B)
|
|---|
| 249 | IF (FB*FA.GT.0.0) GOTO 9000
|
|---|
| 250 | FC=FB
|
|---|
| 251 | C
|
|---|
| 252 | DO 1000 ITER=1,ITMAX
|
|---|
| 253 | IF (FB*FC.GT.0.0) THEN
|
|---|
| 254 | C = A
|
|---|
| 255 | FC = FA
|
|---|
| 256 | D = B-A
|
|---|
| 257 | E = D
|
|---|
| 258 | ENDIF
|
|---|
| 259 | IF (ABS(FC).LT.ABS(FB)) THEN
|
|---|
| 260 | A = B
|
|---|
| 261 | B = C
|
|---|
| 262 | C = A
|
|---|
| 263 | FA = FB
|
|---|
| 264 | FB = FC
|
|---|
| 265 | FC = FA
|
|---|
| 266 | ENDIF
|
|---|
| 267 | TOL1 = 2.0*EPS*ABS(B)+0.5*TOL
|
|---|
| 268 | XM = 0.5*(C-B)
|
|---|
| 269 | IF (ABS(XM).LE.TOL1.OR.FB.EQ.0.0) THEN
|
|---|
| 270 | ZBRENT = B
|
|---|
| 271 | return
|
|---|
| 272 | ENDIF
|
|---|
| 273 | IF (ABS(E).GT.TOL1.AND. ABS(FA).GT.ABS(FB)) THEN
|
|---|
| 274 | C
|
|---|
| 275 | C Attempt inverse quadratic interpolation
|
|---|
| 276 | C
|
|---|
| 277 | S=FB/FA
|
|---|
| 278 | IF (A.EQ.C) THEN
|
|---|
| 279 | P=2.0*XM*S
|
|---|
| 280 | Q=1.0-S
|
|---|
| 281 | ELSE
|
|---|
| 282 | Q=FA/FC
|
|---|
| 283 | R=FB/FC
|
|---|
| 284 | P=S*( 2.0*XM*Q*(Q-R) - (B-A)*(R-1.0) )
|
|---|
| 285 | Q=(Q-1.0)*(R-1.0)*(S-1.0)
|
|---|
| 286 | ENDIF
|
|---|
| 287 | C
|
|---|
| 288 | C Check whether in bounds...
|
|---|
| 289 | C
|
|---|
| 290 | IF (P.GT.0.0) Q = -Q
|
|---|
| 291 | P = ABS(P)
|
|---|
| 292 | IF (2.0*P.LT.MIN(3.0*XM*Q-ABS(TOL1*Q),ABS(E*Q))) THEN
|
|---|
| 293 | C Accept interpolation.
|
|---|
| 294 | E=D
|
|---|
| 295 | D=P/Q
|
|---|
| 296 | ELSE
|
|---|
| 297 | C Interpolation failed, use bisection.
|
|---|
| 298 | D=XM
|
|---|
| 299 | E=D
|
|---|
| 300 | ENDIF
|
|---|
| 301 | ELSE
|
|---|
| 302 | C Bounds decreasing too slowly, use bisection.
|
|---|
| 303 | D=XM
|
|---|
| 304 | E=D
|
|---|
| 305 | ENDIF
|
|---|
| 306 | A=B
|
|---|
| 307 | FA=FB
|
|---|
| 308 | IF (ABS(D).GT.TOL1) THEN
|
|---|
| 309 | C Evaluate new trial root
|
|---|
| 310 | B=B+D
|
|---|
| 311 | ELSE
|
|---|
| 312 | B=B+SIGN(TOL1,XM)
|
|---|
| 313 | ENDIF
|
|---|
| 314 | FB=FNC(B)
|
|---|
| 315 | 1000 CONTINUE
|
|---|
| 316 | C
|
|---|
| 317 | C
|
|---|
| 318 | 9000 CONTINUE
|
|---|
| 319 | write(6 ,*) 'Exceeding maximum number of iterations.'
|
|---|
| 320 | C WRITE(21,*) 'Exceeding maximum number of iterations.'
|
|---|
| 321 | ZBRENT=B
|
|---|
| 322 | return
|
|---|
| 323 | end
|
|---|
| 324 | FUNCTION FNC(ETA)
|
|---|
| 325 | include 'SIZE'
|
|---|
| 326 | include 'INPUT'
|
|---|
| 327 | COMMON /SRFCEI/ IEL,IFCE
|
|---|
| 328 | COMMON /SRFCER/ X0(3),DX(3)
|
|---|
| 329 | COMMON /SRFCEL/ SUCCES
|
|---|
| 330 | DIMENSION X3(3)
|
|---|
| 331 | integer icalld
|
|---|
| 332 | save icalld
|
|---|
| 333 | data icalld /0/
|
|---|
| 334 | C
|
|---|
| 335 | IF (CCURVE(IFCE,IEL).EQ.'s') THEN
|
|---|
| 336 | xctr = CURVE(1,IFCE,IEL)
|
|---|
| 337 | yctr = CURVE(2,IFCE,IEL)
|
|---|
| 338 | zctr = CURVE(3,IFCE,IEL)
|
|---|
| 339 | RADIUS = CURVE(4,IFCE,IEL)
|
|---|
| 340 | X = X0(1) + DX(1)*ETA - XCTR
|
|---|
| 341 | Y = X0(2) + DX(2)*ETA - YCTR
|
|---|
| 342 | Z = X0(3) + DX(3)*ETA - ZCTR
|
|---|
| 343 | C
|
|---|
| 344 | FNC = (RADIUS**2 - (X**2+Y**2+Z**2))/(RADIUS**2)
|
|---|
| 345 | C
|
|---|
| 346 | ENDIF
|
|---|
| 347 | C
|
|---|
| 348 | return
|
|---|
| 349 | end
|
|---|
| 350 | c-----------------------------------------------------------------------
|
|---|
| 351 | subroutine setdef
|
|---|
| 352 | C-------------------------------------------------------------------
|
|---|
| 353 | C
|
|---|
| 354 | C Set up deformed element logical switches
|
|---|
| 355 | C
|
|---|
| 356 | C-------------------------------------------------------------------
|
|---|
| 357 | include 'SIZE'
|
|---|
| 358 | include 'INPUT'
|
|---|
| 359 | DIMENSION XCC(8),YCC(8),ZCC(8)
|
|---|
| 360 | DIMENSION INDX(8)
|
|---|
| 361 | REAL VEC(3,12)
|
|---|
| 362 | LOGICAL IFVCHK
|
|---|
| 363 | C
|
|---|
| 364 | COMMON /FASTMD/ IFDFRM(LELT), IFFAST(LELT), IFH2, IFSOLV
|
|---|
| 365 | LOGICAL IFDFRM, IFFAST, IFH2, IFSOLV
|
|---|
| 366 | C
|
|---|
| 367 | C Corner notation:
|
|---|
| 368 | C
|
|---|
| 369 | C 4+-----+3 ^ Y
|
|---|
| 370 | C / /| |
|
|---|
| 371 | C / / | |
|
|---|
| 372 | C 8+-----+7 +2 +----> X
|
|---|
| 373 | C | | / /
|
|---|
| 374 | C | |/ /
|
|---|
| 375 | C 5+-----+6 Z
|
|---|
| 376 | C
|
|---|
| 377 | C
|
|---|
| 378 | DO 10 IE=1,NELT
|
|---|
| 379 | IFDFRM(IE)=.TRUE.
|
|---|
| 380 | 10 CONTINUE
|
|---|
| 381 | C
|
|---|
| 382 | IF (IFMVBD) return
|
|---|
| 383 | c
|
|---|
| 384 | c Force IFDFRM=.true. for all elements (for timing purposes only)
|
|---|
| 385 | c
|
|---|
| 386 | IF (param(59).ne.0.and.nio.eq.0)
|
|---|
| 387 | $ write(6,*) 'NOTE: All elements deformed , param(59) ^=0'
|
|---|
| 388 | IF (param(59).ne.0) return
|
|---|
| 389 | C
|
|---|
| 390 | C Check against cases which won't allow for savings in HMHOLTZ
|
|---|
| 391 | C
|
|---|
| 392 | INDX(1)=1
|
|---|
| 393 | INDX(2)=2
|
|---|
| 394 | INDX(3)=4
|
|---|
| 395 | INDX(4)=3
|
|---|
| 396 | INDX(5)=5
|
|---|
| 397 | INDX(6)=6
|
|---|
| 398 | INDX(7)=8
|
|---|
| 399 | INDX(8)=7
|
|---|
| 400 | C
|
|---|
| 401 | C Check for deformation (rotation is acceptable).
|
|---|
| 402 | C
|
|---|
| 403 | DO 500 IE=1,NELT
|
|---|
| 404 | C
|
|---|
| 405 | call rzero(vec,36)
|
|---|
| 406 | IF (IF3D) THEN
|
|---|
| 407 | DO 100 IEDG=1,8
|
|---|
| 408 | IF(CCURVE(IEDG,IE).NE.' ') THEN
|
|---|
| 409 | IFDFRM(IE)=.TRUE.
|
|---|
| 410 | GOTO 500
|
|---|
| 411 | ENDIF
|
|---|
| 412 | 100 CONTINUE
|
|---|
| 413 | C
|
|---|
| 414 | DO 105 I=1,8
|
|---|
| 415 | XCC(I)=XC(INDX(I),IE)
|
|---|
| 416 | YCC(I)=YC(INDX(I),IE)
|
|---|
| 417 | ZCC(I)=ZC(INDX(I),IE)
|
|---|
| 418 | 105 CONTINUE
|
|---|
| 419 | C
|
|---|
| 420 | DO 110 I=1,4
|
|---|
| 421 | VEC(1,I)=XCC(2*I)-XCC(2*I-1)
|
|---|
| 422 | VEC(2,I)=YCC(2*I)-YCC(2*I-1)
|
|---|
| 423 | VEC(3,I)=ZCC(2*I)-ZCC(2*I-1)
|
|---|
| 424 | 110 CONTINUE
|
|---|
| 425 | C
|
|---|
| 426 | I1=4
|
|---|
| 427 | DO 120 I=0,1
|
|---|
| 428 | DO 120 J=0,1
|
|---|
| 429 | I1=I1+1
|
|---|
| 430 | I2=4*I+J+3
|
|---|
| 431 | VEC(1,I1)=XCC(I2)-XCC(I2-2)
|
|---|
| 432 | VEC(2,I1)=YCC(I2)-YCC(I2-2)
|
|---|
| 433 | VEC(3,I1)=ZCC(I2)-ZCC(I2-2)
|
|---|
| 434 | 120 CONTINUE
|
|---|
| 435 | C
|
|---|
| 436 | I1=8
|
|---|
| 437 | DO 130 I=5,8
|
|---|
| 438 | I1=I1+1
|
|---|
| 439 | VEC(1,I1)=XCC(I)-XCC(I-4)
|
|---|
| 440 | VEC(2,I1)=YCC(I)-YCC(I-4)
|
|---|
| 441 | VEC(3,I1)=ZCC(I)-ZCC(I-4)
|
|---|
| 442 | 130 CONTINUE
|
|---|
| 443 | C
|
|---|
| 444 | DO 140 I=1,12
|
|---|
| 445 | VECLEN = VEC(1,I)**2 + VEC(2,I)**2 + VEC(3,I)**2
|
|---|
| 446 | VECLEN = SQRT(VECLEN)
|
|---|
| 447 | VEC(1,I)=VEC(1,I)/VECLEN
|
|---|
| 448 | VEC(2,I)=VEC(2,I)/VECLEN
|
|---|
| 449 | VEC(3,I)=VEC(3,I)/VECLEN
|
|---|
| 450 | 140 CONTINUE
|
|---|
| 451 | C
|
|---|
| 452 | C Check the dot product of the adjacent edges to see that it is zero.
|
|---|
| 453 | C
|
|---|
| 454 | IFDFRM(IE)=.FALSE.
|
|---|
| 455 | IF ( IFVCHK(VEC,1,5, 9) ) IFDFRM(IE)=.TRUE.
|
|---|
| 456 | IF ( IFVCHK(VEC,1,6,10) ) IFDFRM(IE)=.TRUE.
|
|---|
| 457 | IF ( IFVCHK(VEC,2,5,11) ) IFDFRM(IE)=.TRUE.
|
|---|
| 458 | IF ( IFVCHK(VEC,2,6,12) ) IFDFRM(IE)=.TRUE.
|
|---|
| 459 | IF ( IFVCHK(VEC,3,7, 9) ) IFDFRM(IE)=.TRUE.
|
|---|
| 460 | IF ( IFVCHK(VEC,3,8,10) ) IFDFRM(IE)=.TRUE.
|
|---|
| 461 | IF ( IFVCHK(VEC,4,7,11) ) IFDFRM(IE)=.TRUE.
|
|---|
| 462 | IF ( IFVCHK(VEC,4,8,12) ) IFDFRM(IE)=.TRUE.
|
|---|
| 463 | C
|
|---|
| 464 | C Check the 2D case....
|
|---|
| 465 | C
|
|---|
| 466 | ELSE
|
|---|
| 467 | C
|
|---|
| 468 | DO 200 IEDG=1,4
|
|---|
| 469 | IF(CCURVE(IEDG,IE).NE.' ') THEN
|
|---|
| 470 | IFDFRM(IE)=.TRUE.
|
|---|
| 471 | GOTO 500
|
|---|
| 472 | ENDIF
|
|---|
| 473 | 200 CONTINUE
|
|---|
| 474 | C
|
|---|
| 475 | DO 205 I=1,4
|
|---|
| 476 | XCC(I)=XC(INDX(I),IE)
|
|---|
| 477 | YCC(I)=YC(INDX(I),IE)
|
|---|
| 478 | 205 CONTINUE
|
|---|
| 479 | C
|
|---|
| 480 | VEC(1,1)=XCC(2)-XCC(1)
|
|---|
| 481 | VEC(1,2)=XCC(4)-XCC(3)
|
|---|
| 482 | VEC(1,3)=XCC(3)-XCC(1)
|
|---|
| 483 | VEC(1,4)=XCC(4)-XCC(2)
|
|---|
| 484 | VEC(1,5)=0.0
|
|---|
| 485 | VEC(2,1)=YCC(2)-YCC(1)
|
|---|
| 486 | VEC(2,2)=YCC(4)-YCC(3)
|
|---|
| 487 | VEC(2,3)=YCC(3)-YCC(1)
|
|---|
| 488 | VEC(2,4)=YCC(4)-YCC(2)
|
|---|
| 489 | VEC(2,5)=0.0
|
|---|
| 490 | C
|
|---|
| 491 | DO 220 I=1,4
|
|---|
| 492 | VECLEN = VEC(1,I)**2 + VEC(2,I)**2
|
|---|
| 493 | VECLEN = SQRT(VECLEN)
|
|---|
| 494 | VEC(1,I)=VEC(1,I)/VECLEN
|
|---|
| 495 | VEC(2,I)=VEC(2,I)/VECLEN
|
|---|
| 496 | 220 CONTINUE
|
|---|
| 497 | C
|
|---|
| 498 | C Check the dot product of the adjacent edges to see that it is zero.
|
|---|
| 499 | C
|
|---|
| 500 | IFDFRM(IE)=.FALSE.
|
|---|
| 501 | IF ( IFVCHK(VEC,1,3,5) ) IFDFRM(IE)=.TRUE.
|
|---|
| 502 | IF ( IFVCHK(VEC,1,4,5) ) IFDFRM(IE)=.TRUE.
|
|---|
| 503 | IF ( IFVCHK(VEC,2,3,5) ) IFDFRM(IE)=.TRUE.
|
|---|
| 504 | IF ( IFVCHK(VEC,2,4,5) ) IFDFRM(IE)=.TRUE.
|
|---|
| 505 | ENDIF
|
|---|
| 506 | 500 CONTINUE
|
|---|
| 507 | return
|
|---|
| 508 | end
|
|---|
| 509 | LOGICAL FUNCTION IFVCHK(VEC,I1,I2,I3)
|
|---|
| 510 | C
|
|---|
| 511 | C Take the dot product of the three components of VEC to see if it's zero.
|
|---|
| 512 | C
|
|---|
| 513 | DIMENSION VEC(3,12)
|
|---|
| 514 | LOGICAL IFTMP
|
|---|
| 515 | C
|
|---|
| 516 | IFTMP=.FALSE.
|
|---|
| 517 | EPSM=1.0E-06
|
|---|
| 518 | C
|
|---|
| 519 | DOT1=VEC(1,I1)*VEC(1,I2)+VEC(2,I1)*VEC(2,I2)+VEC(3,I1)*VEC(3,I2)
|
|---|
| 520 | DOT2=VEC(1,I2)*VEC(1,I3)+VEC(2,I2)*VEC(2,I3)+VEC(3,I2)*VEC(3,I3)
|
|---|
| 521 | DOT3=VEC(1,I1)*VEC(1,I3)+VEC(2,I1)*VEC(2,I3)+VEC(3,I1)*VEC(3,I3)
|
|---|
| 522 | C
|
|---|
| 523 | DOT1=ABS(DOT1)
|
|---|
| 524 | DOT2=ABS(DOT2)
|
|---|
| 525 | DOT3=ABS(DOT3)
|
|---|
| 526 | DOT=DOT1+DOT2+DOT3
|
|---|
| 527 | IF (DOT.GT.EPSM) IFTMP=.TRUE.
|
|---|
| 528 | C
|
|---|
| 529 | IFVCHK=IFTMP
|
|---|
| 530 | return
|
|---|
| 531 | end
|
|---|
| 532 | c-----------------------------------------------------------------------
|
|---|
| 533 | subroutine gencoor (xm3,ym3,zm3)
|
|---|
| 534 | C-----------------------------------------------------------------------
|
|---|
| 535 | C
|
|---|
| 536 | C Generate xyz coordinates for all elements.
|
|---|
| 537 | C Velocity formulation : mesh 3 is used
|
|---|
| 538 | C Stress formulation : mesh 1 is used
|
|---|
| 539 | C
|
|---|
| 540 | C-----------------------------------------------------------------------
|
|---|
| 541 | include 'SIZE'
|
|---|
| 542 | include 'GEOM'
|
|---|
| 543 | include 'INPUT'
|
|---|
| 544 | DIMENSION XM3(LX3,LY3,LZ3,1),YM3(LX3,LY3,LZ3,1),ZM3(LX3,LY3,LZ3,1)
|
|---|
| 545 | C
|
|---|
| 546 | C Select appropriate mesh
|
|---|
| 547 | C
|
|---|
| 548 | IF ( IFGMSH3 ) THEN
|
|---|
| 549 | CALL GENXYZ (XM3,YM3,ZM3,lx3,ly3,lz3)
|
|---|
| 550 | ELSE
|
|---|
| 551 | CALL GENXYZ (XM1,YM1,ZM1,lx1,ly1,lz1)
|
|---|
| 552 | ENDIF
|
|---|
| 553 | C
|
|---|
| 554 | return
|
|---|
| 555 | end
|
|---|
| 556 | c-----------------------------------------------------------------------
|
|---|
| 557 | subroutine genxyz (xml,yml,zml,nxl,nyl,nzl)
|
|---|
| 558 | C
|
|---|
| 559 | include 'SIZE'
|
|---|
| 560 | include 'WZ'
|
|---|
| 561 | include 'GEOM'
|
|---|
| 562 | include 'TOPOL'
|
|---|
| 563 | include 'INPUT'
|
|---|
| 564 | include 'PARALLEL'
|
|---|
| 565 |
|
|---|
| 566 | real xml(nxl,nyl,nzl,1),yml(nxl,nyl,nzl,1),zml(nxl,nyl,nzl,1)
|
|---|
| 567 |
|
|---|
| 568 | C Note : CTMP1 is used in this format in several subsequent routines
|
|---|
| 569 | common /ctmp1/ h(lx1,3,2),xcrved(lx1),ycrved(ly1),zcrved(lz1)
|
|---|
| 570 | $ , zgml(lx1,3),work(3,lx1,lz1)
|
|---|
| 571 |
|
|---|
| 572 | parameter (ldw=2*lx1*ly1*lz1)
|
|---|
| 573 | common /ctmp0/ w(ldw)
|
|---|
| 574 |
|
|---|
| 575 | character*1 ccv
|
|---|
| 576 |
|
|---|
| 577 | c Initialize geometry arrays with bi- triquadratic deformations
|
|---|
| 578 | call linquad(xml,yml,zml,nxl,nyl,nzl)
|
|---|
| 579 |
|
|---|
| 580 |
|
|---|
| 581 | do ie=1,nelt
|
|---|
| 582 |
|
|---|
| 583 | call setzgml (zgml,ie,nxl,nyl,nzl,ifaxis)
|
|---|
| 584 | call sethmat (h,zgml,nxl,nyl,nzl)
|
|---|
| 585 |
|
|---|
| 586 | c Deform surfaces - general 3D deformations
|
|---|
| 587 | c - extruded geometry deformations
|
|---|
| 588 | nfaces = 2*ldim
|
|---|
| 589 | do iface=1,nfaces
|
|---|
| 590 | ccv = ccurve(iface,ie)
|
|---|
| 591 | if (ccv.eq.'s')
|
|---|
| 592 | $ call sphsrf(xml,yml,zml,iface,ie,nxl,nyl,nzl,work)
|
|---|
| 593 | if (ccv.eq.'e')
|
|---|
| 594 | $ call gensrf(xml,yml,zml,iface,ie,nxl,nyl,nzl,zgml)
|
|---|
| 595 | enddo
|
|---|
| 596 |
|
|---|
| 597 | do isid=1,8
|
|---|
| 598 | ccv = ccurve(isid,ie)
|
|---|
| 599 | if (ccv.eq.'C') call arcsrf(xml,yml,zml,nxl,nyl,nzl,ie,isid)
|
|---|
| 600 | enddo
|
|---|
| 601 |
|
|---|
| 602 | enddo
|
|---|
| 603 |
|
|---|
| 604 | c call user_srf(xml,yml,zml,nxl,nyl,nzl)
|
|---|
| 605 | c call opcopy(xm1,ym1,zm1,xml,yml,zml)
|
|---|
| 606 | c call outpost(xml,yml,zml,xml,yml,' ')
|
|---|
| 607 | c call exitt
|
|---|
| 608 | C
|
|---|
| 609 | return
|
|---|
| 610 | end
|
|---|
| 611 | c-----------------------------------------------------------------------
|
|---|
| 612 | subroutine sethmat(h,zgml,nxl,nyl,nzl)
|
|---|
| 613 |
|
|---|
| 614 | include 'SIZE'
|
|---|
| 615 | include 'INPUT' ! if3d
|
|---|
| 616 |
|
|---|
| 617 | real h(lx1,3,2),zgml(lx1,3)
|
|---|
| 618 |
|
|---|
| 619 | do 10 ix=1,nxl
|
|---|
| 620 | h(ix,1,1)=(1.0-zgml(ix,1))*0.5
|
|---|
| 621 | h(ix,1,2)=(1.0+zgml(ix,1))*0.5
|
|---|
| 622 | 10 continue
|
|---|
| 623 | do 20 iy=1,nyl
|
|---|
| 624 | h(iy,2,1)=(1.0-zgml(iy,2))*0.5
|
|---|
| 625 | h(iy,2,2)=(1.0+zgml(iy,2))*0.5
|
|---|
| 626 | 20 continue
|
|---|
| 627 | if (if3d) then
|
|---|
| 628 | do 30 iz=1,nzl
|
|---|
| 629 | h(iz,3,1)=(1.0-zgml(iz,3))*0.5
|
|---|
| 630 | h(iz,3,2)=(1.0+zgml(iz,3))*0.5
|
|---|
| 631 | 30 continue
|
|---|
| 632 | else
|
|---|
| 633 | call rone(h(1,3,1),nzl)
|
|---|
| 634 | call rone(h(1,3,2),nzl)
|
|---|
| 635 | endif
|
|---|
| 636 |
|
|---|
| 637 | return
|
|---|
| 638 | end
|
|---|
| 639 | c-----------------------------------------------------------------------
|
|---|
| 640 | subroutine setzgml (zgml,e,nxl,nyl,nzl,ifaxl)
|
|---|
| 641 |
|
|---|
| 642 | include 'SIZE'
|
|---|
| 643 | include 'WZ'
|
|---|
| 644 | include 'GEOM'
|
|---|
| 645 |
|
|---|
| 646 | real zgml(lx1,3)
|
|---|
| 647 | integer e
|
|---|
| 648 | logical ifaxl
|
|---|
| 649 |
|
|---|
| 650 | call rzero (zgml,3*lx1)
|
|---|
| 651 |
|
|---|
| 652 |
|
|---|
| 653 | if (nxl.eq.3 .and. .not. ifaxl) then
|
|---|
| 654 | do k=1,3
|
|---|
| 655 | zgml(1,k) = -1
|
|---|
| 656 | zgml(2,k) = 0
|
|---|
| 657 | zgml(3,k) = 1
|
|---|
| 658 | enddo
|
|---|
| 659 | elseif (ifgmsh3.and.nxl.eq.lx3) then
|
|---|
| 660 | call copy(zgml(1,1),zgm3(1,1),lx3)
|
|---|
| 661 | call copy(zgml(1,2),zgm3(1,2),ly3)
|
|---|
| 662 | call copy(zgml(1,3),zgm3(1,3),lz3)
|
|---|
| 663 | if (ifaxl .and. ifrzer(e)) call copy(zgml(1,2),zam3,ly3)
|
|---|
| 664 | elseif (nxl.eq.lx1) then
|
|---|
| 665 | call copy(zgml(1,1),zgm1(1,1),lx1)
|
|---|
| 666 | call copy(zgml(1,2),zgm1(1,2),ly1)
|
|---|
| 667 | call copy(zgml(1,3),zgm1(1,3),lz1)
|
|---|
| 668 | if (ifaxl .and. ifrzer(e)) call copy(zgml(1,2),zam1,ly1)
|
|---|
| 669 | else
|
|---|
| 670 | call exitti('ABORT setzgml! $',nxl)
|
|---|
| 671 | endif
|
|---|
| 672 |
|
|---|
| 673 | return
|
|---|
| 674 | end
|
|---|
| 675 | c-----------------------------------------------------------------------
|
|---|
| 676 | subroutine sphsrf(xml,yml,zml,ifce,ie,nx,ny,nz,xysrf)
|
|---|
| 677 | C
|
|---|
| 678 | C 5 Aug 1988 19:29:52
|
|---|
| 679 | C
|
|---|
| 680 | C Program to generate spherical shell elements for NEKTON
|
|---|
| 681 | C input. Paul F. Fischer
|
|---|
| 682 | C
|
|---|
| 683 | include 'SIZE'
|
|---|
| 684 | include 'INPUT'
|
|---|
| 685 | include 'WZ'
|
|---|
| 686 | include 'TOPOL'
|
|---|
| 687 | DIMENSION XML(NX,NY,NZ,1),YML(NX,NY,NZ,1),ZML(NX,NY,NZ,1)
|
|---|
| 688 | DIMENSION XYSRF(3,NX,NZ)
|
|---|
| 689 | C
|
|---|
| 690 | COMMON /CTMP1/ H(LX1,3,2),XCRVED(LX1),YCRVED(LY1),ZCRVED(LZ1)
|
|---|
| 691 | $ , ZGML(LX1,3),WORK(3,LX1,LZ1)
|
|---|
| 692 | COMMON /CTMP0/ XCV(3,2,2),VN1(3),VN2(3)
|
|---|
| 693 | $ ,X1(3),X2(3),X3(3),DX(3)
|
|---|
| 694 | DIMENSION IOPP(3),NXX(3)
|
|---|
| 695 | c
|
|---|
| 696 | c
|
|---|
| 697 | c These are representative nodes on a given face, and their opposites
|
|---|
| 698 | c
|
|---|
| 699 | integer cface(2,6)
|
|---|
| 700 | save cface
|
|---|
| 701 | data cface / 1,4 , 2,1 , 3,2 , 4,3 , 1,5 , 5,1 /
|
|---|
| 702 | real vout(3),vsph(3)
|
|---|
| 703 | logical ifconcv
|
|---|
| 704 | c
|
|---|
| 705 | C
|
|---|
| 706 | C Determine geometric parameters
|
|---|
| 707 | C
|
|---|
| 708 | NXM1 = NX-1
|
|---|
| 709 | NYM1 = NY-1
|
|---|
| 710 | NXY = NX*NZ
|
|---|
| 711 | NXY3 = 3*NX*NZ
|
|---|
| 712 | XCTR = CURVE(1,IFCE,IE)
|
|---|
| 713 | YCTR = CURVE(2,IFCE,IE)
|
|---|
| 714 | ZCTR = CURVE(3,IFCE,IE)
|
|---|
| 715 | RADIUS = CURVE(4,IFCE,IE)
|
|---|
| 716 | IFACE = EFACE1(IFCE)
|
|---|
| 717 | C
|
|---|
| 718 | C Generate (normalized) corner vectors XCV(1,i,j):
|
|---|
| 719 | C
|
|---|
| 720 | CALL CRN3D(XCV,XC(1,IE),YC(1,IE),ZC(1,IE),CURVE(1,IFCE,IE),IFACE)
|
|---|
| 721 | C
|
|---|
| 722 | C Generate edge vectors on the sphere RR=1.0,
|
|---|
| 723 | C for (r,s) = (-1,*),(1,*),(*,-1),(*,1)
|
|---|
| 724 | C
|
|---|
| 725 | CALL EDG3D(XYSRF,XCV(1,1,1),XCV(1,1,2), 1, 1, 1,NY,NX,NY)
|
|---|
| 726 | CALL EDG3D(XYSRF,XCV(1,2,1),XCV(1,2,2),NX,NX, 1,NY,NX,NY)
|
|---|
| 727 | CALL EDG3D(XYSRF,XCV(1,1,1),XCV(1,2,1), 1,NX, 1, 1,NX,NY)
|
|---|
| 728 | CALL EDG3D(XYSRF,XCV(1,1,2),XCV(1,2,2), 1,NX,NY,NY,NX,NY)
|
|---|
| 729 | C
|
|---|
| 730 | C Generate intersection vectors for (i,j)
|
|---|
| 731 | C
|
|---|
| 732 | C quick check on sign of curvature: (pff , 12/08/00)
|
|---|
| 733 | c
|
|---|
| 734 | c
|
|---|
| 735 | ivtx = cface(1,ifce)
|
|---|
| 736 | ivto = cface(2,ifce)
|
|---|
| 737 | vout(1) = xc(ivtx,ie)-xc(ivto,ie)
|
|---|
| 738 | vout(2) = yc(ivtx,ie)-yc(ivto,ie)
|
|---|
| 739 | vout(3) = zc(ivtx,ie)-zc(ivto,ie)
|
|---|
| 740 | c
|
|---|
| 741 | vsph(1) = xc(ivtx,ie)-xctr
|
|---|
| 742 | vsph(2) = yc(ivtx,ie)-yctr
|
|---|
| 743 | vsph(3) = zc(ivtx,ie)-zctr
|
|---|
| 744 | ifconcv = .true.
|
|---|
| 745 | sign = DOT(vsph,vout,3)
|
|---|
| 746 | if (sign.gt.0) ifconcv = .false.
|
|---|
| 747 | c write(6,*) 'THIS IS SIGN:',sign
|
|---|
| 748 | c
|
|---|
| 749 | DO 200 J=2,NYM1
|
|---|
| 750 | CALL CROSS(VN1,XYSRF(1,1,J),XYSRF(1,NX,J))
|
|---|
| 751 | DO 200 I=2,NXM1
|
|---|
| 752 | CALL CROSS(VN2,XYSRF(1,I,1),XYSRF(1,I,NY))
|
|---|
| 753 | if (ifconcv) then
|
|---|
| 754 | c IF (IFACE.EQ.1.OR.IFACE.EQ.4.OR.IFACE.EQ.5) THEN
|
|---|
| 755 | CALL CROSS(XYSRF(1,I,J),VN2,VN1)
|
|---|
| 756 | ELSE
|
|---|
| 757 | CALL CROSS(XYSRF(1,I,J),VN1,VN2)
|
|---|
| 758 | ENDIF
|
|---|
| 759 | 200 CONTINUE
|
|---|
| 760 | C
|
|---|
| 761 | C Normalize all vectors to the unit sphere.
|
|---|
| 762 | C
|
|---|
| 763 | DO 300 I=1,NXY
|
|---|
| 764 | CALL NORM3D(XYSRF(1,I,1))
|
|---|
| 765 | 300 CONTINUE
|
|---|
| 766 | C
|
|---|
| 767 | C Scale by actual radius
|
|---|
| 768 | C
|
|---|
| 769 | CALL CMULT(XYSRF,RADIUS,NXY3)
|
|---|
| 770 | C
|
|---|
| 771 | C Add back the sphere center offset
|
|---|
| 772 | C
|
|---|
| 773 | DO 400 I=1,NXY
|
|---|
| 774 | XYSRF(1,I,1)=XYSRF(1,I,1)+XCTR
|
|---|
| 775 | XYSRF(2,I,1)=XYSRF(2,I,1)+YCTR
|
|---|
| 776 | XYSRF(3,I,1)=XYSRF(3,I,1)+ZCTR
|
|---|
| 777 | 400 CONTINUE
|
|---|
| 778 | C
|
|---|
| 779 | C
|
|---|
| 780 | C Transpose data, if necessary
|
|---|
| 781 | C
|
|---|
| 782 | IF (IFACE.EQ.1.OR.IFACE.EQ.4.OR.IFACE.EQ.5) THEN
|
|---|
| 783 | DO 500 J=1 ,NY
|
|---|
| 784 | DO 500 I=J+1,NX
|
|---|
| 785 | TMP=XYSRF(1,I,J)
|
|---|
| 786 | XYSRF(1,I,J)=XYSRF(1,J,I)
|
|---|
| 787 | XYSRF(1,J,I)=TMP
|
|---|
| 788 | TMP=XYSRF(2,I,J)
|
|---|
| 789 | XYSRF(2,I,J)=XYSRF(2,J,I)
|
|---|
| 790 | XYSRF(2,J,I)=TMP
|
|---|
| 791 | TMP=XYSRF(3,I,J)
|
|---|
| 792 | XYSRF(3,I,J)=XYSRF(3,J,I)
|
|---|
| 793 | XYSRF(3,J,I)=TMP
|
|---|
| 794 | 500 CONTINUE
|
|---|
| 795 | ENDIF
|
|---|
| 796 | C
|
|---|
| 797 | C Compute surface deflection and perturbation due to face IFACE
|
|---|
| 798 | C
|
|---|
| 799 | CALL DSSET(NX,NY,NZ)
|
|---|
| 800 | JS1 = SKPDAT(1,IFACE)
|
|---|
| 801 | JF1 = SKPDAT(2,IFACE)
|
|---|
| 802 | JSKIP1 = SKPDAT(3,IFACE)
|
|---|
| 803 | JS2 = SKPDAT(4,IFACE)
|
|---|
| 804 | JF2 = SKPDAT(5,IFACE)
|
|---|
| 805 | JSKIP2 = SKPDAT(6,IFACE)
|
|---|
| 806 | C
|
|---|
| 807 | IOPP(1) = NX-1
|
|---|
| 808 | IOPP(2) = NX*(NY-1)
|
|---|
| 809 | IOPP(3) = NX*NY*(NZ-1)
|
|---|
| 810 | NXX(1) = NX
|
|---|
| 811 | NXX(2) = NY
|
|---|
| 812 | NXX(3) = NZ
|
|---|
| 813 | IDIR = 2*MOD(IFACE,2) - 1
|
|---|
| 814 | IFC2 = (IFACE+1)/2
|
|---|
| 815 | DELT = 0.0
|
|---|
| 816 | I=0
|
|---|
| 817 | DO 700 J2=JS2,JF2,JSKIP2
|
|---|
| 818 | DO 700 J1=JS1,JF1,JSKIP1
|
|---|
| 819 | I=I+1
|
|---|
| 820 | JOPP = J1 + IOPP(IFC2)*IDIR
|
|---|
| 821 | X2(1) = XML(J1,J2,1,IE)
|
|---|
| 822 | X2(2) = YML(J1,J2,1,IE)
|
|---|
| 823 | X2(3) = ZML(J1,J2,1,IE)
|
|---|
| 824 | C
|
|---|
| 825 | DX(1) = XYSRF(1,I,1)-X2(1)
|
|---|
| 826 | DX(2) = XYSRF(2,I,1)-X2(2)
|
|---|
| 827 | DX(3) = XYSRF(3,I,1)-X2(3)
|
|---|
| 828 | C
|
|---|
| 829 | NXS = NXX(IFC2)
|
|---|
| 830 | JOFF = (J1-JOPP)/(NXS-1)
|
|---|
| 831 | DO 600 IX = 2,NXS
|
|---|
| 832 | J = JOPP + JOFF*(IX-1)
|
|---|
| 833 | ZETA = 0.5*(ZGML(IX,IFC2) + 1.0)
|
|---|
| 834 | XML(J,J2,1,IE) = XML(J,J2,1,IE)+DX(1)*ZETA
|
|---|
| 835 | YML(J,J2,1,IE) = YML(J,J2,1,IE)+DX(2)*ZETA
|
|---|
| 836 | ZML(J,J2,1,IE) = ZML(J,J2,1,IE)+DX(3)*ZETA
|
|---|
| 837 | 600 CONTINUE
|
|---|
| 838 | 700 CONTINUE
|
|---|
| 839 | C
|
|---|
| 840 | return
|
|---|
| 841 | end
|
|---|
| 842 | c-----------------------------------------------------------------------
|
|---|
| 843 | subroutine edg3d(xysrf,x1,x2,i1,i2,j1,j2,nx,ny)
|
|---|
| 844 | C
|
|---|
| 845 | C Generate XYZ vector along an edge of a surface.
|
|---|
| 846 | C
|
|---|
| 847 | include 'SIZE'
|
|---|
| 848 | COMMON /CTMP1/ H(LX1,3,2),XCRVED(LX1),YCRVED(LY1),ZCRVED(LZ1)
|
|---|
| 849 | $ , ZGML(LX1,3),WORK(3,LX1,LZ1)
|
|---|
| 850 | C
|
|---|
| 851 | DIMENSION XYSRF(3,NX,NY)
|
|---|
| 852 | DIMENSION X1(3),X2(3)
|
|---|
| 853 | REAL U1(3),U2(3),VN(3),B(3)
|
|---|
| 854 | C
|
|---|
| 855 | C Normalize incoming vectors
|
|---|
| 856 | C
|
|---|
| 857 | CALL COPY (U1,X1,3)
|
|---|
| 858 | CALL COPY (U2,X2,3)
|
|---|
| 859 | CALL NORM3D (U1)
|
|---|
| 860 | CALL NORM3D (U2)
|
|---|
| 861 | C
|
|---|
| 862 | C Find normal to the plane and tangent to the curve.
|
|---|
| 863 | C
|
|---|
| 864 | CALL CROSS(VN,X1,X2)
|
|---|
| 865 | CALL CROSS( B,VN,X1)
|
|---|
| 866 | CALL NORM3D (VN)
|
|---|
| 867 | CALL NORM3D (B)
|
|---|
| 868 | C
|
|---|
| 869 | CTHETA = DOT(U1,U2,3)
|
|---|
| 870 | THETA = ACOS(CTHETA)
|
|---|
| 871 | C
|
|---|
| 872 | IJ = 0
|
|---|
| 873 | DO 200 J=J1,J2
|
|---|
| 874 | DO 200 I=I1,I2
|
|---|
| 875 | IJ = IJ + 1
|
|---|
| 876 | THETAP = 0.5*THETA*(ZGML(IJ,1)+1.0)
|
|---|
| 877 | CTP = COS(THETAP)
|
|---|
| 878 | STP = SIN(THETAP)
|
|---|
| 879 | DO 200 IV = 1,3
|
|---|
| 880 | XYSRF(IV,I,J) = CTP*U1(IV) + STP*B(IV)
|
|---|
| 881 | 200 CONTINUE
|
|---|
| 882 | return
|
|---|
| 883 | end
|
|---|
| 884 | REAL FUNCTION DOT(V1,V2,N)
|
|---|
| 885 | C
|
|---|
| 886 | C Compute Cartesian vector dot product.
|
|---|
| 887 | C
|
|---|
| 888 | DIMENSION V1(N),V2(N)
|
|---|
| 889 | C
|
|---|
| 890 | SUM = 0
|
|---|
| 891 | DO 100 I=1,N
|
|---|
| 892 | SUM = SUM + V1(I)*V2(I)
|
|---|
| 893 | 100 CONTINUE
|
|---|
| 894 | DOT = SUM
|
|---|
| 895 | return
|
|---|
| 896 | end
|
|---|
| 897 | c-----------------------------------------------------------------------
|
|---|
| 898 | subroutine cross(v1,v2,v3)
|
|---|
| 899 | C
|
|---|
| 900 | C Compute Cartesian vector dot product.
|
|---|
| 901 | C
|
|---|
| 902 | DIMENSION V1(3),V2(3),V3(3)
|
|---|
| 903 | C
|
|---|
| 904 | V1(1) = V2(2)*V3(3) - V2(3)*V3(2)
|
|---|
| 905 | V1(2) = V2(3)*V3(1) - V2(1)*V3(3)
|
|---|
| 906 | V1(3) = V2(1)*V3(2) - V2(2)*V3(1)
|
|---|
| 907 | C
|
|---|
| 908 | return
|
|---|
| 909 | end
|
|---|
| 910 | c-----------------------------------------------------------------------
|
|---|
| 911 | subroutine norm3d(v1)
|
|---|
| 912 | C
|
|---|
| 913 | C Compute Cartesian vector dot product.
|
|---|
| 914 | C
|
|---|
| 915 | DIMENSION V1(3)
|
|---|
| 916 | C
|
|---|
| 917 | VLNGTH = DOT(V1,V1,3)
|
|---|
| 918 | VLNGTH = SQRT(VLNGTH)
|
|---|
| 919 | if (vlngth.gt.0) then
|
|---|
| 920 | V1(1) = V1(1) / VLNGTH
|
|---|
| 921 | V1(2) = V1(2) / VLNGTH
|
|---|
| 922 | V1(3) = V1(3) / VLNGTH
|
|---|
| 923 | endif
|
|---|
| 924 | C
|
|---|
| 925 | return
|
|---|
| 926 | end
|
|---|
| 927 | c-----------------------------------------------------------------------
|
|---|
| 928 | subroutine crn3d(xcv,xc,yc,zc,curve,iface)
|
|---|
| 929 | include 'SIZE'
|
|---|
| 930 | include 'TOPOL'
|
|---|
| 931 | DIMENSION XCV(3,2,2),XC(8),YC(8),ZC(8),CURVE(4)
|
|---|
| 932 | DIMENSION INDVTX(4,6)
|
|---|
| 933 | SAVE INDVTX
|
|---|
| 934 | DATA INDVTX / 1,5,3,7 , 2,4,6,8 , 1,2,5,6
|
|---|
| 935 | $ , 3,7,4,8 , 1,3,2,4 , 5,6,7,8 /
|
|---|
| 936 | C
|
|---|
| 937 | EPS = 1.0E-4
|
|---|
| 938 | XCTR = CURVE(1)
|
|---|
| 939 | YCTR = CURVE(2)
|
|---|
| 940 | ZCTR = CURVE(3)
|
|---|
| 941 | RADIUS = CURVE(4)
|
|---|
| 942 | C
|
|---|
| 943 | DO 10 I=1,4
|
|---|
| 944 | J=INDVTX(I,IFACE)
|
|---|
| 945 | K=INDX(J)
|
|---|
| 946 | XCV(1,I,1)=XC(K)-XCTR
|
|---|
| 947 | XCV(2,I,1)=YC(K)-YCTR
|
|---|
| 948 | XCV(3,I,1)=ZC(K)-ZCTR
|
|---|
| 949 | 10 CONTINUE
|
|---|
| 950 | C
|
|---|
| 951 | C Check to ensure that these points are indeed on the sphere.
|
|---|
| 952 | C
|
|---|
| 953 | IF (RADIUS.LE.0.0) THEN
|
|---|
| 954 | write(6,20) NID,XCTR,YCTR,ZCTR,IFACE
|
|---|
| 955 | 20 FORMAT(I5,'ERROR: Sphere of radius zero requested.'
|
|---|
| 956 | $ ,/,5X,'EXITING in CRN3D',3E12.4,I3)
|
|---|
| 957 | call exitt
|
|---|
| 958 | ELSE
|
|---|
| 959 | DO 40 I=1,4
|
|---|
| 960 | RADT=XCV(1,I,1)**2+XCV(2,I,1)**2+XCV(3,I,1)**2
|
|---|
| 961 | RADT=SQRT(RADT)
|
|---|
| 962 | TEST=ABS(RADT-RADIUS)/RADIUS
|
|---|
| 963 | IF (TEST.GT.EPS) THEN
|
|---|
| 964 | write(6,30) NID
|
|---|
| 965 | $ ,RADT,RADIUS,XCV(1,I,1),XCV(2,I,1),XCV(3,I,1)
|
|---|
| 966 | 30 FORMAT(I5,'ERROR: Element vertex not on requested sphere.'
|
|---|
| 967 | $ ,/,5X,'EXITING in CRN3D',5E12.4)
|
|---|
| 968 | call exitt
|
|---|
| 969 | ENDIF
|
|---|
| 970 | 40 CONTINUE
|
|---|
| 971 | ENDIF
|
|---|
| 972 | C
|
|---|
| 973 | return
|
|---|
| 974 | end
|
|---|
| 975 | c-----------------------------------------------------------------------
|
|---|
| 976 | subroutine rotxyz
|
|---|
| 977 | C-----------------------------------------------------------------------
|
|---|
| 978 | C
|
|---|
| 979 | C Establish rotation of undeformed elements.
|
|---|
| 980 | C Used for fast evaluation of D*x and DT*x.
|
|---|
| 981 | C Currently used for 2-d problems.
|
|---|
| 982 | C
|
|---|
| 983 | C-----------------------------------------------------------------------
|
|---|
| 984 | include 'SIZE'
|
|---|
| 985 | include 'INPUT'
|
|---|
| 986 | include 'GEOM'
|
|---|
| 987 | include 'ESOLV'
|
|---|
| 988 | COMMON /FASTMD/ IFDFRM(LELT), IFFAST(LELT), IFH2, IFSOLV
|
|---|
| 989 | LOGICAL IFDFRM, IFFAST, IFH2, IFSOLV
|
|---|
| 990 | C
|
|---|
| 991 | IF (IFMVBD) return
|
|---|
| 992 | IF (ldim.EQ.3) return
|
|---|
| 993 | C
|
|---|
| 994 | EPS = 1.E-6
|
|---|
| 995 | EPSINV = 1./EPS
|
|---|
| 996 | NXYZ1 = lx1*ly1*lz1
|
|---|
| 997 | DO 100 IEL=1,NELV
|
|---|
| 998 | IF (.NOT.IFDFRM(IEL)) THEN
|
|---|
| 999 | RXAVER = VLSUM(RXM1(1,1,1,IEL),NXYZ1)/(NXYZ1)
|
|---|
| 1000 | SXAVER = VLSUM(SXM1(1,1,1,IEL),NXYZ1)/(NXYZ1)
|
|---|
| 1001 | IF (SXAVER.NE.0.) THEN
|
|---|
| 1002 | ROTXY = ABS(SXAVER/RXAVER)
|
|---|
| 1003 | IF (ROTXY.LT.EPS) THEN
|
|---|
| 1004 | IFALGN(IEL) = .TRUE.
|
|---|
| 1005 | IFRSXY(IEL) = .TRUE.
|
|---|
| 1006 | ELSEIF (ROTXY.GT.EPSINV) THEN
|
|---|
| 1007 | IFALGN(IEL) = .TRUE.
|
|---|
| 1008 | IFRSXY(IEL) = .FALSE.
|
|---|
| 1009 | ELSE
|
|---|
| 1010 | IFALGN(IEL) = .FALSE.
|
|---|
| 1011 | IFRSXY(IEL) = .FALSE.
|
|---|
| 1012 | ENDIF
|
|---|
| 1013 | ELSE
|
|---|
| 1014 | IFALGN(IEL) = .TRUE.
|
|---|
| 1015 | IFRSXY(IEL) = .TRUE.
|
|---|
| 1016 | ENDIF
|
|---|
| 1017 | ENDIF
|
|---|
| 1018 | 100 CONTINUE
|
|---|
| 1019 | return
|
|---|
| 1020 | end
|
|---|
| 1021 | c-----------------------------------------------------------------------
|
|---|
| 1022 | subroutine gensrf(XML,YML,ZML,IFCE,IE,MX,MY,MZ,zgml)
|
|---|
| 1023 | C
|
|---|
| 1024 | C 9 Mar 1994
|
|---|
| 1025 | C
|
|---|
| 1026 | C Program to generate surface deformations for NEKTON
|
|---|
| 1027 | C input. Paul F. Fischer
|
|---|
| 1028 | C
|
|---|
| 1029 | c include 'basics.inc'
|
|---|
| 1030 | include 'SIZE'
|
|---|
| 1031 | include 'INPUT'
|
|---|
| 1032 | include 'WZ'
|
|---|
| 1033 | include 'TOPOL'
|
|---|
| 1034 | C
|
|---|
| 1035 | DIMENSION XML(MX,MY,MZ,1),YML(MX,MY,MZ,1),ZML(MX,MY,MZ,1)
|
|---|
| 1036 | $ ,ZGML(MX,3)
|
|---|
| 1037 | C
|
|---|
| 1038 | real IOPP(3),MXX(3),X0(3),DX(3)
|
|---|
| 1039 | C
|
|---|
| 1040 | C
|
|---|
| 1041 | C Algorithm: .Project original point onto surface S
|
|---|
| 1042 | C .Apply Gordon Hall to vector of points between x_s and
|
|---|
| 1043 | C opposite face
|
|---|
| 1044 | C
|
|---|
| 1045 | C
|
|---|
| 1046 | CALL DSSET(MX,MY,MZ)
|
|---|
| 1047 | C
|
|---|
| 1048 | IFACE = EFACE1(IFCE)
|
|---|
| 1049 | c
|
|---|
| 1050 | c Beware!! SKPDAT different from preprocessor/postprocessor!
|
|---|
| 1051 | c
|
|---|
| 1052 | JS1 = SKPDAT(1,IFACE)
|
|---|
| 1053 | JF1 = SKPDAT(2,IFACE)
|
|---|
| 1054 | JSKIP1 = SKPDAT(3,IFACE)
|
|---|
| 1055 | JS2 = SKPDAT(4,IFACE)
|
|---|
| 1056 | JF2 = SKPDAT(5,IFACE)
|
|---|
| 1057 | JSKIP2 = SKPDAT(6,IFACE)
|
|---|
| 1058 | c
|
|---|
| 1059 | IOPP(1) = MX-1
|
|---|
| 1060 | IOPP(2) = MX*(MY-1)
|
|---|
| 1061 | IOPP(3) = MX*MY*(MZ-1)
|
|---|
| 1062 | MXX(1) = MX
|
|---|
| 1063 | MXX(2) = MY
|
|---|
| 1064 | MXX(3) = MZ
|
|---|
| 1065 | IDIR = 2*MOD(IFACE,2) - 1
|
|---|
| 1066 | IFC2 = (IFACE+1)/2
|
|---|
| 1067 | I=0
|
|---|
| 1068 | C
|
|---|
| 1069 | C Find a characteristic length scale for initializing secant method
|
|---|
| 1070 | C
|
|---|
| 1071 | x0(1) = xml(js1,js2,1,ie)
|
|---|
| 1072 | x0(2) = yml(js1,js2,1,ie)
|
|---|
| 1073 | x0(3) = zml(js1,js2,1,ie)
|
|---|
| 1074 | rmin = 1.0e16
|
|---|
| 1075 | c
|
|---|
| 1076 | c
|
|---|
| 1077 | c
|
|---|
| 1078 | DO 100 J2=JS2,JF2,JSKIP2
|
|---|
| 1079 | DO 100 J1=JS1,JF1,JSKIP1
|
|---|
| 1080 | if (j1.ne.js1.or.j2.ne.js2) then
|
|---|
| 1081 | r2 = (x0(1) - xml(j1,j2,1,ie))**2
|
|---|
| 1082 | $ + (x0(2) - yml(j1,j2,1,ie))**2
|
|---|
| 1083 | $ + (x0(3) - zml(j1,j2,1,ie))**2
|
|---|
| 1084 | rmin = min(r2,rmin)
|
|---|
| 1085 | endif
|
|---|
| 1086 | 100 CONTINUE
|
|---|
| 1087 | dxc = 0.05*sqrt(rmin)
|
|---|
| 1088 | C
|
|---|
| 1089 | C Project each point on this surface onto curved surface
|
|---|
| 1090 | C
|
|---|
| 1091 | DO 300 J2=JS2,JF2,JSKIP2
|
|---|
| 1092 | DO 300 J1=JS1,JF1,JSKIP1
|
|---|
| 1093 | I=I+1
|
|---|
| 1094 | JOPP = J1 + IOPP(IFC2)*IDIR
|
|---|
| 1095 | X0(1) = XML(J1,J2,1,IE)
|
|---|
| 1096 | X0(2) = YML(J1,J2,1,IE)
|
|---|
| 1097 | X0(3) = ZML(J1,J2,1,IE)
|
|---|
| 1098 | C
|
|---|
| 1099 | call prjects(x0,dxc,curve(1,ifce,ie),ccurve(ifce,ie))
|
|---|
| 1100 | DX(1) = X0(1)-xml(j1,j2,1,ie)
|
|---|
| 1101 | DX(2) = X0(2)-yml(j1,j2,1,ie)
|
|---|
| 1102 | DX(3) = X0(3)-zml(j1,j2,1,ie)
|
|---|
| 1103 | MXS = MXX(IFC2)
|
|---|
| 1104 | JOFF = (J1-JOPP)/(MXS-1)
|
|---|
| 1105 | DO 200 IX = 2,MXS
|
|---|
| 1106 | J = JOPP + JOFF*(IX-1)
|
|---|
| 1107 | ZETA = 0.5*(ZGML(IX,1) + 1.0)
|
|---|
| 1108 | XML(J,J2,1,IE) = XML(J,J2,1,IE)+DX(1)*ZETA
|
|---|
| 1109 | YML(J,J2,1,IE) = YML(J,J2,1,IE)+DX(2)*ZETA
|
|---|
| 1110 | ZML(J,J2,1,IE) = ZML(J,J2,1,IE)+DX(3)*ZETA
|
|---|
| 1111 | 200 CONTINUE
|
|---|
| 1112 | 300 CONTINUE
|
|---|
| 1113 | C
|
|---|
| 1114 | return
|
|---|
| 1115 | end
|
|---|
| 1116 | c-----------------------------------------------------------------------
|
|---|
| 1117 | subroutine prjects(x0,dxc,c,cc)
|
|---|
| 1118 | c
|
|---|
| 1119 | c Project the point x0 onto surface described by characteristics
|
|---|
| 1120 | c given in the array c and cc.
|
|---|
| 1121 | c
|
|---|
| 1122 | c dxc - characteristic length scale used to estimate gradient.
|
|---|
| 1123 | c
|
|---|
| 1124 | real x0(3)
|
|---|
| 1125 | real c(5)
|
|---|
| 1126 | character*1 cc
|
|---|
| 1127 | real x1(3)
|
|---|
| 1128 | logical if3d
|
|---|
| 1129 | c
|
|---|
| 1130 | if3d = .true.
|
|---|
| 1131 | if (dxc.le.0.0) then
|
|---|
| 1132 | write(6,*) 'invalid dxc',dxc,x0
|
|---|
| 1133 | write(6,*) 'Abandoning prjects'
|
|---|
| 1134 | return
|
|---|
| 1135 | endif
|
|---|
| 1136 | c
|
|---|
| 1137 | call copy(x1,x0,3)
|
|---|
| 1138 | R0 = ressrf(x0,c,cc)
|
|---|
| 1139 | if (r0.eq.0) return
|
|---|
| 1140 | c
|
|---|
| 1141 | c Must at least use ctr differencing to capture symmetry!
|
|---|
| 1142 | c
|
|---|
| 1143 | x1(1) = x0(1) - dxc
|
|---|
| 1144 | R1 = ressrf(x1,c,cc)
|
|---|
| 1145 | x1(1) = x0(1) + dxc
|
|---|
| 1146 | R2 = ressrf(x1,c,cc)
|
|---|
| 1147 | x1(1) = x0(1)
|
|---|
| 1148 | Rx = 0.5*(R2-R1)/dxc
|
|---|
| 1149 | c
|
|---|
| 1150 | x1(2) = x0(2) - dxc
|
|---|
| 1151 | R1 = ressrf(x1,c,cc)/dxc
|
|---|
| 1152 | x1(2) = x0(2) + dxc
|
|---|
| 1153 | R2 = ressrf(x1,c,cc)/dxc
|
|---|
| 1154 | x1(2) = x0(2)
|
|---|
| 1155 | Ry = 0.5*(R2-R1)/dxc
|
|---|
| 1156 | c
|
|---|
| 1157 | if (if3d) then
|
|---|
| 1158 | x1(3) = x0(3) - dxc
|
|---|
| 1159 | R1 = ressrf(x1,c,cc)/dxc
|
|---|
| 1160 | x1(3) = x0(3) + dxc
|
|---|
| 1161 | R2 = ressrf(x1,c,cc)/dxc
|
|---|
| 1162 | Rz = 0.5*(R2-R1)/dxc
|
|---|
| 1163 | endif
|
|---|
| 1164 | Rnorm2 = Rx**2 + Ry**2 + Rz**2
|
|---|
| 1165 | alpha = - R0/Rnorm2
|
|---|
| 1166 | c
|
|---|
| 1167 | c Apply secant method: Use an initial segment twice expected length
|
|---|
| 1168 | c
|
|---|
| 1169 | x1(1) = x0(1) + 2.0*Rx * alpha
|
|---|
| 1170 | x1(2) = x0(2) + 2.0*Ry * alpha
|
|---|
| 1171 | x1(3) = x0(3) + 2.0*Rz * alpha
|
|---|
| 1172 | call srfind(x1,x0,c,cc)
|
|---|
| 1173 | c
|
|---|
| 1174 | c write(6,6) cc,c(2),c(3),x0,x1
|
|---|
| 1175 | c 6 format(1x,a1,1x,2f5.2,3f9.4,3x,3f9.4)
|
|---|
| 1176 | c
|
|---|
| 1177 | call copy(x0,x1,3)
|
|---|
| 1178 | c
|
|---|
| 1179 | return
|
|---|
| 1180 | end
|
|---|
| 1181 | c-----------------------------------------------------------------------
|
|---|
| 1182 | subroutine srfind(x1,x0,c,cc)
|
|---|
| 1183 | real x1(3),x0(3)
|
|---|
| 1184 | real c(5)
|
|---|
| 1185 | character*1 cc
|
|---|
| 1186 | c
|
|---|
| 1187 | c Find point on line segment that intersects the ellipsoid
|
|---|
| 1188 | c specified by:
|
|---|
| 1189 | c (x/a)**2 + (y/b)**2 + (z/b)**2 = 1
|
|---|
| 1190 | c
|
|---|
| 1191 | c
|
|---|
| 1192 | c Algorithm: 4 rounds of secant x_k+1 = x_k - f/f'
|
|---|
| 1193 | c
|
|---|
| 1194 | a0 = 0.0
|
|---|
| 1195 | a1 = 1.0
|
|---|
| 1196 | r0 = ressrf(x0,c,cc)
|
|---|
| 1197 | dx = x1(1) - x0(1)
|
|---|
| 1198 | dy = x1(2) - x0(2)
|
|---|
| 1199 | dz = x1(3) - x0(3)
|
|---|
| 1200 | c write(6,*) 'dxyz',dx,dy,dz
|
|---|
| 1201 | c write(6,*) 'cc ',x0,cc,c(2),c(3)
|
|---|
| 1202 | do 10 i=1,9
|
|---|
| 1203 | r1 = ressrf(x1,c,cc)
|
|---|
| 1204 | if (r1.ne.r0) then
|
|---|
| 1205 | da = r1*(a1-a0)/(r1-r0)
|
|---|
| 1206 | r0 = r1
|
|---|
| 1207 | a0 = a1
|
|---|
| 1208 | a1 = a1 - da
|
|---|
| 1209 | endif
|
|---|
| 1210 | x1(1) = x0(1) + a1*dx
|
|---|
| 1211 | x1(2) = x0(2) + a1*dy
|
|---|
| 1212 | x1(3) = x0(3) + a1*dz
|
|---|
| 1213 | 10 continue
|
|---|
| 1214 | c write(6,*) ' r1',r1,r0,a1
|
|---|
| 1215 | return
|
|---|
| 1216 | end
|
|---|
| 1217 | c-----------------------------------------------------------------------
|
|---|
| 1218 | function ressrf(x,c,cc)
|
|---|
| 1219 | real x(3)
|
|---|
| 1220 | real c(5)
|
|---|
| 1221 | character*1 cc
|
|---|
| 1222 | c
|
|---|
| 1223 | ressrf = 0.0
|
|---|
| 1224 | if (cc.eq.'e') then
|
|---|
| 1225 | a = c(2)
|
|---|
| 1226 | b = c(3)
|
|---|
| 1227 | ressrf = 1.0 - (x(1)/a)**2 - (x(2)/b)**2 - (x(3)/b)**2
|
|---|
| 1228 | return
|
|---|
| 1229 | endif
|
|---|
| 1230 | c
|
|---|
| 1231 | return
|
|---|
| 1232 | end
|
|---|
| 1233 | c-----------------------------------------------------------------------
|
|---|
| 1234 | subroutine linquad(xl,yl,zl,nxl,nyl,nzl)
|
|---|
| 1235 |
|
|---|
| 1236 | include 'SIZE'
|
|---|
| 1237 | include 'WZ'
|
|---|
| 1238 | include 'GEOM'
|
|---|
| 1239 | include 'TOPOL'
|
|---|
| 1240 | include 'INPUT'
|
|---|
| 1241 | include 'PARALLEL'
|
|---|
| 1242 |
|
|---|
| 1243 | real xl(nxl*nyl*nzl,1),yl(nxl*nyl*nzl,1),zl(nxl*nyl*nzl,1)
|
|---|
| 1244 |
|
|---|
| 1245 | integer e
|
|---|
| 1246 | logical ifmid
|
|---|
| 1247 |
|
|---|
| 1248 | nedge = 4 + 8*(ldim-2)
|
|---|
| 1249 |
|
|---|
| 1250 | do e=1,nelt ! Loop over all elements
|
|---|
| 1251 |
|
|---|
| 1252 | ifmid = .false.
|
|---|
| 1253 | do k=1,nedge
|
|---|
| 1254 | if (ccurve(k,e).eq.'m') ifmid = .true.
|
|---|
| 1255 | enddo
|
|---|
| 1256 |
|
|---|
| 1257 | if (lx1.eq.2) ifmid = .false.
|
|---|
| 1258 | if (ifmid) then
|
|---|
| 1259 | call xyzquad(xl(1,e),yl(1,e),zl(1,e),nxl,nyl,nzl,e)
|
|---|
| 1260 | else
|
|---|
| 1261 | call xyzlin (xl(1,e),yl(1,e),zl(1,e),nxl,nyl,nzl,e,ifaxis)
|
|---|
| 1262 | endif
|
|---|
| 1263 | enddo
|
|---|
| 1264 |
|
|---|
| 1265 | return
|
|---|
| 1266 | end
|
|---|
| 1267 | c-----------------------------------------------------------------------
|
|---|
| 1268 | subroutine xyzlin(xl,yl,zl,nxl,nyl,nzl,e,ifaxl)
|
|---|
| 1269 | c Generate bi- or trilinear mesh
|
|---|
| 1270 |
|
|---|
| 1271 | include 'SIZE'
|
|---|
| 1272 | include 'INPUT'
|
|---|
| 1273 |
|
|---|
| 1274 | real xl(nxl,nyl,nzl),yl(nxl,nyl,nzl),zl(nxl,nyl,nzl)
|
|---|
| 1275 | integer e
|
|---|
| 1276 | logical ifaxl ! local ifaxis specification
|
|---|
| 1277 |
|
|---|
| 1278 | c Preprocessor Corner notation: Symmetric Corner notation:
|
|---|
| 1279 | c
|
|---|
| 1280 | c 4+-----+3 ^ s 3+-----+4 ^ s
|
|---|
| 1281 | c / /| | / /| |
|
|---|
| 1282 | c / / | | / / | |
|
|---|
| 1283 | c 8+-----+7 +2 +----> r 7+-----+8 +2 +----> r
|
|---|
| 1284 | c | | / / | | / /
|
|---|
| 1285 | c | |/ / | |/ /
|
|---|
| 1286 | c 5+-----+6 t 5+-----+6 t
|
|---|
| 1287 |
|
|---|
| 1288 | integer indx(8)
|
|---|
| 1289 | save indx
|
|---|
| 1290 | data indx / 1,2,4,3,5,6,8,7 /
|
|---|
| 1291 |
|
|---|
| 1292 | parameter (ldw=4*lx1*ly1*lz1)
|
|---|
| 1293 | common /ctmp0/ xcb(2,2,2),ycb(2,2,2),zcb(2,2,2),w(ldw)
|
|---|
| 1294 |
|
|---|
| 1295 | common /cxyzl/ zgml(lx1,3),jx (lx1*2),jy (lx1*2),jz (lx1*2)
|
|---|
| 1296 | $ ,jxt(lx1*2),jyt(lx1*2),jzt(lx1*2)
|
|---|
| 1297 | $ ,zlin(2)
|
|---|
| 1298 | real jx,jy,jz,jxt,jyt,jzt
|
|---|
| 1299 |
|
|---|
| 1300 | call setzgml (zgml,e,nxl,nyl,nzl,ifaxl)
|
|---|
| 1301 |
|
|---|
| 1302 | zlin(1) = -1
|
|---|
| 1303 | zlin(2) = 1
|
|---|
| 1304 |
|
|---|
| 1305 | k = 1
|
|---|
| 1306 | do i=1,nxl
|
|---|
| 1307 | call fd_weights_full(zgml(i,1),zlin,1,0,jxt(k))
|
|---|
| 1308 | call fd_weights_full(zgml(i,2),zlin,1,0,jyt(k))
|
|---|
| 1309 | call fd_weights_full(zgml(i,3),zlin,1,0,jzt(k))
|
|---|
| 1310 | k=k+2
|
|---|
| 1311 | enddo
|
|---|
| 1312 | call transpose(jx,nxl,jxt,2)
|
|---|
| 1313 |
|
|---|
| 1314 | ldim2 = 2**ldim
|
|---|
| 1315 | do ix=1,ldim2 ! Convert prex notation to lexicographical
|
|---|
| 1316 | i=indx(ix)
|
|---|
| 1317 | xcb(ix,1,1)=xc(i,e)
|
|---|
| 1318 | ycb(ix,1,1)=yc(i,e)
|
|---|
| 1319 | zcb(ix,1,1)=zc(i,e)
|
|---|
| 1320 | enddo
|
|---|
| 1321 |
|
|---|
| 1322 | c Map R-S-T space into physical X-Y-Z space.
|
|---|
| 1323 |
|
|---|
| 1324 | ! NOTE: Assumes nxl=nyl=nzl !
|
|---|
| 1325 |
|
|---|
| 1326 | call tensr3(xl,nxl,xcb,2,jx,jyt,jzt,w)
|
|---|
| 1327 | call tensr3(yl,nxl,ycb,2,jx,jyt,jzt,w)
|
|---|
| 1328 | call tensr3(zl,nxl,zcb,2,jx,jyt,jzt,w)
|
|---|
| 1329 |
|
|---|
| 1330 | return
|
|---|
| 1331 | end
|
|---|
| 1332 | c-----------------------------------------------------------------------
|
|---|
| 1333 | subroutine xyzquad(xl,yl,zl,nxl,nyl,nzl,e)
|
|---|
| 1334 | c Generate bi- or trilinear mesh
|
|---|
| 1335 |
|
|---|
| 1336 | include 'SIZE'
|
|---|
| 1337 | include 'INPUT'
|
|---|
| 1338 |
|
|---|
| 1339 | real xl(nxl,nyl,nzl),yl(nxl,nyl,nzl),zl(nxl,nyl,nzl)
|
|---|
| 1340 | real xq(27),yq(27),zq(27)
|
|---|
| 1341 | integer e
|
|---|
| 1342 |
|
|---|
| 1343 | parameter (ldw=4*lx1*ly1*lz1)
|
|---|
| 1344 | common /ctmp0/ w(ldw,2),zg(3)
|
|---|
| 1345 |
|
|---|
| 1346 | c Note : CTMP1 is used in this format in several subsequent routines
|
|---|
| 1347 |
|
|---|
| 1348 | integer eindx(12) ! index of 12 edges into 3x3x3 tensor
|
|---|
| 1349 | save eindx ! Follows preprocessor notation..
|
|---|
| 1350 | data eindx / 2 , 6 , 8 , 4
|
|---|
| 1351 | $ , 20 , 24 , 26 , 22
|
|---|
| 1352 | $ , 10 , 12 , 18 , 16 / ! preproc. vtx notation
|
|---|
| 1353 |
|
|---|
| 1354 | common /cxyzl/ zgml(lx1,3),jx (lx1*3),jy (lx1*3),jz (lx1*3)
|
|---|
| 1355 | $ ,jxt(lx1*3),jyt(lx1*3),jzt(lx1*3)
|
|---|
| 1356 | $ ,zquad(3)
|
|---|
| 1357 | real jx,jy,jz,jxt,jyt,jzt
|
|---|
| 1358 |
|
|---|
| 1359 | call xyzlin(xq,yq,zq,3,3,3,e,.false.) ! map bilin to 3x3x3
|
|---|
| 1360 |
|
|---|
| 1361 | nedge = 4 + 8*(ldim-2)
|
|---|
| 1362 |
|
|---|
| 1363 | do k=1,nedge
|
|---|
| 1364 | if (ccurve(k,e).eq.'m') then
|
|---|
| 1365 | j = eindx(k)
|
|---|
| 1366 | xq(j) = curve(1,k,e)
|
|---|
| 1367 | yq(j) = curve(2,k,e)
|
|---|
| 1368 | zq(j) = curve(3,k,e)
|
|---|
| 1369 | endif
|
|---|
| 1370 | enddo
|
|---|
| 1371 |
|
|---|
| 1372 | zg(1) = -1
|
|---|
| 1373 | zg(2) = 0
|
|---|
| 1374 | zg(3) = 1
|
|---|
| 1375 |
|
|---|
| 1376 | if (if3d) then
|
|---|
| 1377 | call gh_face_extend(xq,zg,3,2,w(1,1),w(1,2)) ! 2 --> edge extend
|
|---|
| 1378 | call gh_face_extend(yq,zg,3,2,w(1,1),w(1,2))
|
|---|
| 1379 | call gh_face_extend(zq,zg,3,2,w(1,1),w(1,2))
|
|---|
| 1380 | else
|
|---|
| 1381 | call gh_face_extend_2d(xq,zg,3,2,w(1,1),w(1,2)) ! 2 --> edge extend
|
|---|
| 1382 | call gh_face_extend_2d(yq,zg,3,2,w(1,1),w(1,2))
|
|---|
| 1383 | endif
|
|---|
| 1384 | call clean_xyzq(xq,yq,zq,if3d) ! verify that midside node is in "middle"
|
|---|
| 1385 |
|
|---|
| 1386 |
|
|---|
| 1387 | c Map R-S-T space into physical X-Y-Z space.
|
|---|
| 1388 | ! NOTE: Assumes nxl=nyl=nzl !
|
|---|
| 1389 |
|
|---|
| 1390 | zquad(1) = -1
|
|---|
| 1391 | zquad(2) = 0
|
|---|
| 1392 | zquad(3) = 1
|
|---|
| 1393 |
|
|---|
| 1394 | call setzgml (zgml,e,nxl,nyl,nzl,ifaxis) ! Here we address axisymm.
|
|---|
| 1395 |
|
|---|
| 1396 | k = 1
|
|---|
| 1397 | do i=1,nxl
|
|---|
| 1398 | call fd_weights_full(zgml(i,1),zquad,2,0,jxt(k))
|
|---|
| 1399 | call fd_weights_full(zgml(i,2),zquad,2,0,jyt(k))
|
|---|
| 1400 | call fd_weights_full(zgml(i,3),zquad,2,0,jzt(k))
|
|---|
| 1401 | k=k+3
|
|---|
| 1402 | enddo
|
|---|
| 1403 | call transpose(jx,nxl,jxt,3)
|
|---|
| 1404 |
|
|---|
| 1405 | call tensr3(xl,nxl,xq,3,jx,jyt,jzt,w)
|
|---|
| 1406 | call tensr3(yl,nxl,yq,3,jx,jyt,jzt,w)
|
|---|
| 1407 | call tensr3(zl,nxl,zq,3,jx,jyt,jzt,w)
|
|---|
| 1408 |
|
|---|
| 1409 | return
|
|---|
| 1410 | end
|
|---|
| 1411 | c-----------------------------------------------------------------------
|
|---|
| 1412 | subroutine clean_xyzq(x,y,z,if3d) ! verify that midside node is in "middle"
|
|---|
| 1413 |
|
|---|
| 1414 | real x(3,3,3),y(3,3,3),z(3,3,3)
|
|---|
| 1415 | logical if3d
|
|---|
| 1416 |
|
|---|
| 1417 | integer eindx(12) ! index of 12 edges into 3x3x3 tensor
|
|---|
| 1418 | save eindx ! Follows preprocessor notation..
|
|---|
| 1419 | data eindx / 2 , 6 , 8 , 4
|
|---|
| 1420 | $ , 20 , 24 , 26 , 22
|
|---|
| 1421 | $ , 10 , 12 , 18 , 16 / ! preproc. vtx notation
|
|---|
| 1422 |
|
|---|
| 1423 |
|
|---|
| 1424 | c if (if3d) then ! 12 edges
|
|---|
| 1425 | c else
|
|---|
| 1426 | c endif
|
|---|
| 1427 |
|
|---|
| 1428 | c Here - see routine "fix_m_curve" in prenek for indexing strategy
|
|---|
| 1429 | c
|
|---|
| 1430 | c Note that "fix_m_curve" does not yet perform the actual fix, and
|
|---|
| 1431 | c it too should be updated.
|
|---|
| 1432 |
|
|---|
| 1433 | return
|
|---|
| 1434 | end
|
|---|
| 1435 | c-----------------------------------------------------------------------
|
|---|
| 1436 | subroutine mesh_metrics
|
|---|
| 1437 | include 'SIZE'
|
|---|
| 1438 | include 'TOTAL'
|
|---|
| 1439 |
|
|---|
| 1440 | parameter(nedge = 4 + 8*(ldim-2))
|
|---|
| 1441 | real ledg(nedge)
|
|---|
| 1442 |
|
|---|
| 1443 | nxyz = nx1*ny1*nz1
|
|---|
| 1444 | ntot = nxyz*nelt
|
|---|
| 1445 |
|
|---|
| 1446 | ! aspect ratio
|
|---|
| 1447 | ddmin = 1e20
|
|---|
| 1448 | ddmax = -1
|
|---|
| 1449 | ddavg = 0
|
|---|
| 1450 | do ie = 1,nelt
|
|---|
| 1451 | ledg(1) = dist_xyzc(1,2,ie)
|
|---|
| 1452 | ledg(2) = dist_xyzc(1,4,ie)
|
|---|
| 1453 | ledg(3) = dist_xyzc(2,3,ie)
|
|---|
| 1454 | ledg(4) = dist_xyzc(4,3,ie)
|
|---|
| 1455 | if (ndim.eq.3) then
|
|---|
| 1456 | ledg(5) = dist_xyzc(1,5,ie)
|
|---|
| 1457 | ledg(6) = dist_xyzc(2,6,ie)
|
|---|
| 1458 | ledg(7) = dist_xyzc(4,8,ie)
|
|---|
| 1459 | ledg(8) = dist_xyzc(3,7,ie)
|
|---|
| 1460 |
|
|---|
| 1461 | ledg(9) = dist_xyzc(5,6,ie)
|
|---|
| 1462 | ledg(10) = dist_xyzc(5,8,ie)
|
|---|
| 1463 | ledg(11) = dist_xyzc(8,7,ie)
|
|---|
| 1464 | ledg(12) = dist_xyzc(6,7,ie)
|
|---|
| 1465 | endif
|
|---|
| 1466 |
|
|---|
| 1467 | dratio = vlmax(ledg,nedge)/vlmin(ledg,nedge)
|
|---|
| 1468 | ddmin = min(ddmin,dratio)
|
|---|
| 1469 | ddmax = max(ddmax,dratio)
|
|---|
| 1470 | ddavg = ddavg + dratio
|
|---|
| 1471 | enddo
|
|---|
| 1472 | darmin = glmin(ddmin,1)
|
|---|
| 1473 | darmax = glmax(ddmax,1)
|
|---|
| 1474 | daravg = glsum(ddavg,1)/nelgt
|
|---|
| 1475 |
|
|---|
| 1476 | ! scaled Jac
|
|---|
| 1477 | ddmin = 1e20
|
|---|
| 1478 | ddmax = -1
|
|---|
| 1479 | ddavg = 0
|
|---|
| 1480 | do ie = 1,nelt
|
|---|
| 1481 | dratio = vlmin(JACM1(1,1,1,ie),nxyz)/
|
|---|
| 1482 | $ vlmax(JACM1(1,1,1,ie),nxyz)
|
|---|
| 1483 | ddmin = min(ddmin,dratio)
|
|---|
| 1484 | ddmax = max(ddmax,dratio)
|
|---|
| 1485 | ddavg = ddavg + dratio
|
|---|
| 1486 | enddo
|
|---|
| 1487 | dsjmin = glmin(ddmin,1)
|
|---|
| 1488 | dsjmax = glmax(ddmax,1)
|
|---|
| 1489 | dsjavg = glsum(ddavg,1)/nelgt
|
|---|
| 1490 |
|
|---|
| 1491 | ! dx
|
|---|
| 1492 | ddmin = 1e20
|
|---|
| 1493 | ddmax = -1
|
|---|
| 1494 | ddavg = 0
|
|---|
| 1495 | do ie = 1,nelt
|
|---|
| 1496 | ddmin = min(ddmin,dxmin_e(ie))
|
|---|
| 1497 | ddmax = max(ddmax,dxmax_e(ie))
|
|---|
| 1498 | dtmp = vlsum(JACM1(1,1,1,ie),nxyz)**1./3
|
|---|
| 1499 | ddavg = ddavg + dtmp/(lx1-1)
|
|---|
| 1500 | enddo
|
|---|
| 1501 | dxmin = glmin(ddmin,1)
|
|---|
| 1502 | dxmax = glmax(ddmax,1)
|
|---|
| 1503 | dxavg = glsum(ddavg,1)/nelgt
|
|---|
| 1504 |
|
|---|
| 1505 | if (nid.eq.0) then
|
|---|
| 1506 | write(6,*) 'mesh metrics:'
|
|---|
| 1507 | write(6,'(A,1p2E9.2)') ' GLL grid spacing min/max :',
|
|---|
| 1508 | $ dxmin,dxmax
|
|---|
| 1509 | write(6,'(A,1p3E9.2)') ' scaled Jacobian min/max/avg:',
|
|---|
| 1510 | $ dsjmin,dsjmax,dsjavg
|
|---|
| 1511 | write(6,'(A,1p3E9.2)') ' aspect ratio min/max/avg:',
|
|---|
| 1512 | $ darmin,darmax,daravg
|
|---|
| 1513 | write(6,*)
|
|---|
| 1514 | endif
|
|---|
| 1515 |
|
|---|
| 1516 | return
|
|---|
| 1517 | end
|
|---|
| 1518 | c-----------------------------------------------------------------------
|
|---|
| 1519 | real function dist_xyzc(i,j,ie)
|
|---|
| 1520 | c
|
|---|
| 1521 | c distance between two element corner points
|
|---|
| 1522 | c
|
|---|
| 1523 | include 'SIZE'
|
|---|
| 1524 | include 'INPUT'
|
|---|
| 1525 |
|
|---|
| 1526 | dist_xyzc = (xc(i,ie) - xc(j,ie))**2
|
|---|
| 1527 | dist_xyzc = dist_xyzc + (yc(i,ie) - yc(j,ie))**2
|
|---|
| 1528 | if(ndim.eq.3) dist_xyzc = dist_xyzc + (zc(i,ie) - zc(j,ie))**2
|
|---|
| 1529 | dist_xyzc = sqrt(dist_xyzc)
|
|---|
| 1530 |
|
|---|
| 1531 | return
|
|---|
| 1532 | end
|
|---|
| 1533 | c-----------------------------------------------------------------------
|
|---|
| 1534 | function dxmin_e(e)
|
|---|
| 1535 |
|
|---|
| 1536 | include 'SIZE'
|
|---|
| 1537 | include 'GEOM'
|
|---|
| 1538 | include 'INPUT'
|
|---|
| 1539 |
|
|---|
| 1540 | integer e
|
|---|
| 1541 |
|
|---|
| 1542 | d2m = 1.e20
|
|---|
| 1543 |
|
|---|
| 1544 | if (ldim.eq.3) then
|
|---|
| 1545 | do k=1,nz1-1
|
|---|
| 1546 | do j=1,ny1-1
|
|---|
| 1547 | do i=1,nx1-1
|
|---|
| 1548 | dx = xm1(i+1,j,k,e) - xm1(i,j,k,e)
|
|---|
| 1549 | dy = ym1(i+1,j,k,e) - ym1(i,j,k,e)
|
|---|
| 1550 | dz = zm1(i+1,j,k,e) - zm1(i,j,k,e)
|
|---|
| 1551 | d2 = dx*dx + dy*dy + dz*dz
|
|---|
| 1552 | d2m = min(d2m,d2)
|
|---|
| 1553 |
|
|---|
| 1554 | dx = xm1(i,j+1,k,e) - xm1(i,j,k,e)
|
|---|
| 1555 | dy = ym1(i,j+1,k,e) - ym1(i,j,k,e)
|
|---|
| 1556 | dz = zm1(i,j+1,k,e) - zm1(i,j,k,e)
|
|---|
| 1557 | d2 = dx*dx + dy*dy + dz*dz
|
|---|
| 1558 | d2m = min(d2m,d2)
|
|---|
| 1559 |
|
|---|
| 1560 | dx = xm1(i,j,k+1,e) - xm1(i,j,k,e)
|
|---|
| 1561 | dy = ym1(i,j,k+1,e) - ym1(i,j,k,e)
|
|---|
| 1562 | dz = zm1(i,j,k+1,e) - zm1(i,j,k,e)
|
|---|
| 1563 | d2 = dx*dx + dy*dy + dz*dz
|
|---|
| 1564 | d2m = min(d2m,d2)
|
|---|
| 1565 | enddo
|
|---|
| 1566 | enddo
|
|---|
| 1567 | enddo
|
|---|
| 1568 | else ! 2D
|
|---|
| 1569 | do j=1,ny1-1
|
|---|
| 1570 | do i=1,nx1-1
|
|---|
| 1571 | dx = xm1(i+1,j,1,e) - xm1(i,j,1,e)
|
|---|
| 1572 | dy = ym1(i+1,j,1,e) - ym1(i,j,1,e)
|
|---|
| 1573 | d2 = dx*dx + dy*dy
|
|---|
| 1574 | d2m = min(d2m,d2)
|
|---|
| 1575 |
|
|---|
| 1576 | dx = xm1(i,j+1,1,e) - xm1(i,j,1,e)
|
|---|
| 1577 | dy = ym1(i,j+1,1,e) - ym1(i,j,1,e)
|
|---|
| 1578 | d2 = dx*dx + dy*dy
|
|---|
| 1579 | d2m = min(d2m,d2)
|
|---|
| 1580 | enddo
|
|---|
| 1581 | enddo
|
|---|
| 1582 | endif
|
|---|
| 1583 |
|
|---|
| 1584 | dxmin_e = sqrt(d2m)
|
|---|
| 1585 |
|
|---|
| 1586 | return
|
|---|
| 1587 | end
|
|---|
| 1588 | c-----------------------------------------------------------------------
|
|---|
| 1589 | function dxmax_e(e)
|
|---|
| 1590 |
|
|---|
| 1591 | include 'SIZE'
|
|---|
| 1592 | include 'GEOM'
|
|---|
| 1593 | include 'INPUT'
|
|---|
| 1594 |
|
|---|
| 1595 | integer e
|
|---|
| 1596 |
|
|---|
| 1597 | d2m = -1.e20
|
|---|
| 1598 |
|
|---|
| 1599 | if (ldim.eq.3) then
|
|---|
| 1600 | do k=1,nz1-1
|
|---|
| 1601 | do j=1,ny1-1
|
|---|
| 1602 | do i=1,nx1-1
|
|---|
| 1603 | dx = xm1(i+1,j,k,e) - xm1(i,j,k,e)
|
|---|
| 1604 | dy = ym1(i+1,j,k,e) - ym1(i,j,k,e)
|
|---|
| 1605 | dz = zm1(i+1,j,k,e) - zm1(i,j,k,e)
|
|---|
| 1606 | d2 = dx*dx + dy*dy + dz*dz
|
|---|
| 1607 | d2m = max(d2m,d2)
|
|---|
| 1608 |
|
|---|
| 1609 | dx = xm1(i,j+1,k,e) - xm1(i,j,k,e)
|
|---|
| 1610 | dy = ym1(i,j+1,k,e) - ym1(i,j,k,e)
|
|---|
| 1611 | dz = zm1(i,j+1,k,e) - zm1(i,j,k,e)
|
|---|
| 1612 | d2 = dx*dx + dy*dy + dz*dz
|
|---|
| 1613 | d2m = max(d2m,d2)
|
|---|
| 1614 |
|
|---|
| 1615 | dx = xm1(i,j,k+1,e) - xm1(i,j,k,e)
|
|---|
| 1616 | dy = ym1(i,j,k+1,e) - ym1(i,j,k,e)
|
|---|
| 1617 | dz = zm1(i,j,k+1,e) - zm1(i,j,k,e)
|
|---|
| 1618 | d2 = dx*dx + dy*dy + dz*dz
|
|---|
| 1619 | d2m = max(d2m,d2)
|
|---|
| 1620 | enddo
|
|---|
| 1621 | enddo
|
|---|
| 1622 | enddo
|
|---|
| 1623 | else ! 2D
|
|---|
| 1624 | do j=1,ny1-1
|
|---|
| 1625 | do i=1,nx1-1
|
|---|
| 1626 | dx = xm1(i+1,j,1,e) - xm1(i,j,1,e)
|
|---|
| 1627 | dy = ym1(i+1,j,1,e) - ym1(i,j,1,e)
|
|---|
| 1628 | d2 = dx*dx + dy*dy
|
|---|
| 1629 | d2m = max(d2m,d2)
|
|---|
| 1630 |
|
|---|
| 1631 | dx = xm1(i,j+1,1,e) - xm1(i,j,1,e)
|
|---|
| 1632 | dy = ym1(i,j+1,1,e) - ym1(i,j,1,e)
|
|---|
| 1633 | d2 = dx*dx + dy*dy
|
|---|
| 1634 | d2m = max(d2m,d2)
|
|---|
| 1635 | enddo
|
|---|
| 1636 | enddo
|
|---|
| 1637 | endif
|
|---|
| 1638 |
|
|---|
| 1639 | dxmax_e = sqrt(d2m)
|
|---|
| 1640 |
|
|---|
| 1641 | return
|
|---|
| 1642 | end
|
|---|
| 1643 | c-----------------------------------------------------------------------
|
|---|
