1 subroutine e3bvar (yl, ycl, BCB, shpb, shglb, 2 & xlb, lnode, g1yi, g2yi, 3 & g3yi, WdetJb, bnorm, pres, T, 4 & u1, u2, u3, rho, ei, 5 & cp, rk, 6 & rou, p, tau1n, tau2n, tau3n, 7 & heat) 8c 9c---------------------------------------------------------------------- 10c 11c This routine computes the variables at integration points for 12c the boundary element routine. 13c 14c input: 15c yl (npro,nshl,nflow) : primitive variables (perturbed, no scalars) 16c ycl (npro,nshl,ndof) : primitive variables 17c BCB (npro,nshlb,ndBCB) : Boundary Condition values 18c shpb (npro,nshl) : boundary element shape-functions 19c shglb (npro,nsd,nshl) : boundary element grad-shape-functions 20c xlb (npro,nenl,nsd) : nodal coordinates at current step 21c lnode (nenb) : local nodes on the boundary 22c 23c output: 24c g1yi (npro,nflow) : grad-v in direction 1 25c g2yi (npro,nflow) : grad-v in direction 2 26c g3yi (npro,nflow) : grad-v in direction 3 27c WdetJb (npro) : weighted Jacobian 28c bnorm (npro,nsd) : outward normal 29c pres (npro) : pressure 30c T (npro) : temperature 31c u1 (npro) : x1-velocity component 32c u2 (npro) : x2-velocity component 33c u3 (npro) : x3-velocity component 34c rho (npro) : density 35c ei (npro) : internal energy 36c cp (npro) : specific energy at constant pressure 37c rk (npro) : kinetic energy 38c rou (npro) : BC mass flux 39c p (npro) : BC pressure 40c tau1n (npro) : BC viscous flux 1 41c tau2n (npro) : BC viscous flux 2 42c tau3n (npro) : BC viscous flux 3 43c heat (npro) : BC heat flux 44c 45c 46c Zdenek Johan, Summer 1990. (Modified from e2bvar.f) 47c Zdenek Johan, Winter 1991. (Fortran 90) 48c---------------------------------------------------------------------- 49c 50 include "common.h" 51c 52 dimension yl(npro,nshl,nflow), BCB(npro,nshlb,ndBCB), 53 & ycl(npro,nshl,ndof), 54 & shpb(npro,nshl), 55 & shglb(npro,nsd,nshl), 56 & xlb(npro,nenl,nsd), 57 & lnode(27), g1yi(npro,nflow), 58 & g2yi(npro,nflow), g3yi(npro,nflow), 59 & WdetJb(npro), bnorm(npro,nsd), 60 & pres(npro), T(npro), 61 & u1(npro), u2(npro), 62 & u3(npro), rho(npro), 63 & ei(npro), cp(npro), 64 & rk(npro), 65 & rou(npro), p(npro), 66 & tau1n(npro), tau2n(npro), 67 & tau3n(npro), heat(npro) 68c 69 dimension gl1yi(npro,nflow), gl2yi(npro,nflow), 70 & gl3yi(npro,nflow), dxdxib(npro,nsd,nsd), 71 & dxidxb(npro,nsd,nsd), temp(npro), 72 & temp1(npro), temp2(npro), 73 & temp3(npro) 74c 75 dimension h(npro), cv(npro), 76 & alfap(npro), betaT(npro), 77 & gamb(npro), c(npro), 78 & tmp(npro), 79 & v1(npro,nsd), v2(npro,nsd) 80c 81c.... -------------------> integration variables <-------------------- 82c 83c.... compute the primitive variables at the integration point 84c 85 pres = zero 86 u1 = zero 87 u2 = zero 88 u3 = zero 89 T = zero 90c 91 do n = 1, nshlb 92 nodlcl = lnode(n) 93c 94 pres = pres + shpb(:,nodlcl) * yl(:,nodlcl,1) 95 u1 = u1 + shpb(:,nodlcl) * yl(:,nodlcl,2) 96 u2 = u2 + shpb(:,nodlcl) * yl(:,nodlcl,3) 97 u3 = u3 + shpb(:,nodlcl) * yl(:,nodlcl,4) 98 T = T + shpb(:,nodlcl) * yl(:,nodlcl,5) 99 enddo 100c 101c.... calculate the specific kinetic energy 102c 103 rk = pt5 * ( u1**2 + u2**2 + u3**2 ) 104c 105c.... get the thermodynamic properties 106c 107 if (iLSet .ne. 0)then 108 temp = zero 109 isc=abs(iRANS)+6 110 do n = 1, nshlb 111 temp = temp + shpb(:,n) * ycl(:,n,isc) 112 enddo 113 endif 114 115 ithm = 6 116 if (Navier .eq. 1) ithm = 7 117 call getthm (pres, T, temp, 118 & rk, rho, ei, 119 & h, tmp, cv, 120 & cp, alfap, betaT, 121 & gamb, c) 122c 123c.... ----------------------> Element Metrics <----------------------- 124c 125c.... compute the deformation gradient 126c 127 dxdxib = zero 128c 129 do n = 1, nenl 130 dxdxib(:,1,1) = dxdxib(:,1,1) + xlb(:,n,1) * shglb(:,1,n) 131 dxdxib(:,1,2) = dxdxib(:,1,2) + xlb(:,n,1) * shglb(:,2,n) 132 dxdxib(:,1,3) = dxdxib(:,1,3) + xlb(:,n,1) * shglb(:,3,n) 133 dxdxib(:,2,1) = dxdxib(:,2,1) + xlb(:,n,2) * shglb(:,1,n) 134 dxdxib(:,2,2) = dxdxib(:,2,2) + xlb(:,n,2) * shglb(:,2,n) 135 dxdxib(:,2,3) = dxdxib(:,2,3) + xlb(:,n,2) * shglb(:,3,n) 136 dxdxib(:,3,1) = dxdxib(:,3,1) + xlb(:,n,3) * shglb(:,1,n) 137 dxdxib(:,3,2) = dxdxib(:,3,2) + xlb(:,n,3) * shglb(:,2,n) 138 dxdxib(:,3,3) = dxdxib(:,3,3) + xlb(:,n,3) * shglb(:,3,n) 139 enddo 140c 141c.... compute the normal to the boundary 142c 143c.... compute the normal to the boundary. This is achieved by taking 144c the cross product of two vectors in the plane of the 2-d 145c boundary face. 146 v1 = xlb(:,2,:) - xlb(:,1,:) 147 v2 = xlb(:,3,:) - xlb(:,1,:) 148c 149c.....The following are done in order to correct temp1..3 150c based on the results from compressible code. This is done only 151c for wedges, depending on the bounary face.(tri or quad) 152c 153 if (lcsyst .eq. 4) then 154 temp1 = dxdxib(:,2,1) * dxdxib(:,3,3) - 155 & dxdxib(:,2,3) * dxdxib(:,3,1) 156 temp2 = dxdxib(:,3,1) * dxdxib(:,1,3) - 157 & dxdxib(:,3,3) * dxdxib(:,1,1) 158 temp3 = dxdxib(:,1,1) * dxdxib(:,2,3) - 159 & dxdxib(:,1,3) * dxdxib(:,2,1) 160 161 elseif (lcsyst .eq. 1) then 162 temp1 = v1(:,2) * v2(:,3) - v2(:,2) * v1(:,3) 163 temp2 = v2(:,1) * v1(:,3) - v1(:,1) * v2(:,3) 164 temp3 = v1(:,1) * v2(:,2) - v2(:,1) * v1(:,2) 165 else 166 temp1 = - v1(:,2) * v2(:,3) + v2(:,2) * v1(:,3) 167 temp2 = - v2(:,1) * v1(:,3) + v1(:,1) * v2(:,3) 168 temp3 = - v1(:,1) * v2(:,2) + v2(:,1) * v1(:,2) 169 endif 170c 171 temp = one / sqrt ( temp1**2 + temp2**2 + temp3**2 ) 172 bnorm(:,1) = temp1 * temp 173 bnorm(:,2) = temp2 * temp 174 bnorm(:,3) = temp3 * temp 175c 176 177 if (lcsyst .eq. 3) then 178 WdetJb = (1 - Qwtb(lcsyst,intp)) / (four*temp) 179 elseif (lcsyst .eq. 4) then 180c 181c quad face wedges have a conflict in lnode ordering that makes the 182c normal negative 183c 184c bnorm=-bnorm 185c 186 WdetJb = Qwtb(lcsyst,intp) / temp 187 else 188 WdetJb = Qwtb(lcsyst,intp) / (four*temp) 189 endif 190c 191c.... --------------------------> Grad-V <---------------------------- 192c 193c.... compute grad-v for Navier-Stokes terms 194c 195 if (Navier .eq. 1) then 196c 197c.... compute the inverse of deformation gradient 198c 199 dxidxb(:,1,1) = dxdxib(:,2,2) * dxdxib(:,3,3) 200 & - dxdxib(:,3,2) * dxdxib(:,2,3) 201 dxidxb(:,1,2) = dxdxib(:,3,2) * dxdxib(:,1,3) 202 & - dxdxib(:,1,2) * dxdxib(:,3,3) 203 dxidxb(:,1,3) = dxdxib(:,1,2) * dxdxib(:,2,3) 204 & - dxdxib(:,1,3) * dxdxib(:,2,2) 205 temp = one / ( dxidxb(:,1,1) * dxdxib(:,1,1) 206 & + dxidxb(:,1,2) * dxdxib(:,2,1) 207 & + dxidxb(:,1,3) * dxdxib(:,3,1) ) 208 dxidxb(:,1,1) = dxidxb(:,1,1) * temp 209 dxidxb(:,1,2) = dxidxb(:,1,2) * temp 210 dxidxb(:,1,3) = dxidxb(:,1,3) * temp 211 dxidxb(:,2,1) = (dxdxib(:,2,3) * dxdxib(:,3,1) 212 & - dxdxib(:,2,1) * dxdxib(:,3,3)) * temp 213 dxidxb(:,2,2) = (dxdxib(:,1,1) * dxdxib(:,3,3) 214 & - dxdxib(:,3,1) * dxdxib(:,1,3)) * temp 215 dxidxb(:,2,3) = (dxdxib(:,2,1) * dxdxib(:,1,3) 216 & - dxdxib(:,1,1) * dxdxib(:,2,3)) * temp 217 dxidxb(:,3,1) = (dxdxib(:,2,1) * dxdxib(:,3,2) 218 & - dxdxib(:,2,2) * dxdxib(:,3,1)) * temp 219 dxidxb(:,3,2) = (dxdxib(:,3,1) * dxdxib(:,1,2) 220 & - dxdxib(:,1,1) * dxdxib(:,3,2)) * temp 221 dxidxb(:,3,3) = (dxdxib(:,1,1) * dxdxib(:,2,2) 222 & - dxdxib(:,1,2) * dxdxib(:,2,1)) * temp 223c 224c.... compute local-grad-Y 225c 226 gl1yi = zero 227 gl2yi = zero 228 gl3yi = zero 229c 230 do n = 1, nshl 231 gl1yi(:,1) = gl1yi(:,1) + shglb(:,1,n) * yl(:,n,1) 232 gl1yi(:,2) = gl1yi(:,2) + shglb(:,1,n) * yl(:,n,2) 233 gl1yi(:,3) = gl1yi(:,3) + shglb(:,1,n) * yl(:,n,3) 234 gl1yi(:,4) = gl1yi(:,4) + shglb(:,1,n) * yl(:,n,4) 235 gl1yi(:,5) = gl1yi(:,5) + shglb(:,1,n) * yl(:,n,5) 236c 237 gl2yi(:,1) = gl2yi(:,1) + shglb(:,2,n) * yl(:,n,1) 238 gl2yi(:,2) = gl2yi(:,2) + shglb(:,2,n) * yl(:,n,2) 239 gl2yi(:,3) = gl2yi(:,3) + shglb(:,2,n) * yl(:,n,3) 240 gl2yi(:,4) = gl2yi(:,4) + shglb(:,2,n) * yl(:,n,4) 241 gl2yi(:,5) = gl2yi(:,5) + shglb(:,2,n) * yl(:,n,5) 242c 243 gl3yi(:,1) = gl3yi(:,1) + shglb(:,3,n) * yl(:,n,1) 244 gl3yi(:,2) = gl3yi(:,2) + shglb(:,3,n) * yl(:,n,2) 245 gl3yi(:,3) = gl3yi(:,3) + shglb(:,3,n) * yl(:,n,3) 246 gl3yi(:,4) = gl3yi(:,4) + shglb(:,3,n) * yl(:,n,4) 247 gl3yi(:,5) = gl3yi(:,5) + shglb(:,3,n) * yl(:,n,5) 248 enddo 249c 250c.... convert local-grads to global-grads 251c 252 g1yi(:,2) = dxidxb(:,1,1) * gl1yi(:,2) + 253 & dxidxb(:,2,1) * gl2yi(:,2) + 254 & dxidxb(:,3,1) * gl3yi(:,2) 255 g2yi(:,2) = dxidxb(:,1,2) * gl1yi(:,2) + 256 & dxidxb(:,2,2) * gl2yi(:,2) + 257 & dxidxb(:,3,2) * gl3yi(:,2) 258 g3yi(:,2) = dxidxb(:,1,3) * gl1yi(:,2) + 259 & dxidxb(:,2,3) * gl2yi(:,2) + 260 & dxidxb(:,3,3) * gl3yi(:,2) 261c 262 g1yi(:,3) = dxidxb(:,1,1) * gl1yi(:,3) + 263 & dxidxb(:,2,1) * gl2yi(:,3) + 264 & dxidxb(:,3,1) * gl3yi(:,3) 265 g2yi(:,3) = dxidxb(:,1,2) * gl1yi(:,3) + 266 & dxidxb(:,2,2) * gl2yi(:,3) + 267 & dxidxb(:,3,2) * gl3yi(:,3) 268 g3yi(:,3) = dxidxb(:,1,3) * gl1yi(:,3) + 269 & dxidxb(:,2,3) * gl2yi(:,3) + 270 & dxidxb(:,3,3) * gl3yi(:,3) 271c 272 g1yi(:,4) = dxidxb(:,1,1) * gl1yi(:,4) + 273 & dxidxb(:,2,1) * gl2yi(:,4) + 274 & dxidxb(:,3,1) * gl3yi(:,4) 275 g2yi(:,4) = dxidxb(:,1,2) * gl1yi(:,4) + 276 & dxidxb(:,2,2) * gl2yi(:,4) + 277 & dxidxb(:,3,2) * gl3yi(:,4) 278 g3yi(:,4) = dxidxb(:,1,3) * gl1yi(:,4) + 279 & dxidxb(:,2,3) * gl2yi(:,4) + 280 & dxidxb(:,3,3) * gl3yi(:,4) 281c 282 g1yi(:,5) = dxidxb(:,1,1) * gl1yi(:,5) + 283 & dxidxb(:,2,1) * gl2yi(:,5) + 284 & dxidxb(:,3,1) * gl3yi(:,5) 285 g2yi(:,5) = dxidxb(:,1,2) * gl1yi(:,5) + 286 & dxidxb(:,2,2) * gl2yi(:,5) + 287 & dxidxb(:,3,2) * gl3yi(:,5) 288 g3yi(:,5) = dxidxb(:,1,3) * gl1yi(:,5) + 289 & dxidxb(:,2,3) * gl2yi(:,5) + 290 & dxidxb(:,3,3) * gl3yi(:,5) 291c 292c.... end grad-v 293c 294 endif 295c 296c.... --------------------> Boundary Conditions <-------------------- 297c 298c.... compute the Euler boundary conditions 299c 300 rou = zero 301 p = zero 302c 303 do n = 1, nshlb 304 nodlcl = lnode(n) 305c 306 rou = rou + shpb(:,nodlcl) * BCB(:,n,1) 307 p = p + shpb(:,nodlcl) * BCB(:,n,2) 308 enddo 309c 310c.... compute the Navier-Stokes boundary conditions 311c 312 if (Navier .eq. 1) then 313c 314 tau1n = zero 315 tau2n = zero 316 tau3n = zero 317 heat = zero 318c 319 do n = 1, nshlb 320 nodlcl = lnode(n) 321c 322 tau1n = tau1n + shpb(:,nodlcl) * BCB(:,n,3) 323 tau2n = tau2n + shpb(:,nodlcl) * BCB(:,n,4) 324 tau3n = tau3n + shpb(:,nodlcl) * BCB(:,n,5) 325 heat = heat + shpb(:,nodlcl) * BCB(:,n,6) 326 enddo 327c 328c.... flop count 329c 330 flops = flops + (184+30*nshl+8*nshlb)*npro 331c 332 endif 333c 334c.... flop count 335c 336 flops = flops + (27+18*nshl+14*nshlb)*npro 337c 338c.... return 339c 340 return 341 end 342c 343c 344c 345 subroutine e3bvarSclr(ycl, BCB, shpb, shglb, 346 & xlb, lnode, 347 & u1, u2, u3, 348 & g1yti, g2yti, g3yti, WdetJb, 349 & bnorm, T, rho, cp, rou, 350 & Sclr, SclrF) 351c 352c---------------------------------------------------------------------- 353c 354c This routine computes the variables at integration points for 355c the boundary element routine. 356c 357c input: 358c ycl (npro,nshl,ndof) : Y variables 359c BCB (npro,nenbl,ndBCB) : Boundary Condition values 360c shpb (npro,nen) : boundary element shape-functions 361c shglb (nsd,nen) : boundary element grad-shape-functions 362c xlb (npro,nshl,nsd) : nodal coordinates at current step 363c lnode (nenb) : local nodes on the boundary 364c 365c output: 366c g1yti (npro) 367c g2yti (npro) 368c g3yti (npro) 369c WdetJb (npro) : weighted Jacobian 370c bnorm (npro,nsd) : outward normal 371c T (npro) : temperature 372c rho (npro) : density 373c cp (npro) : specific energy at constant pressure 374c rou (npro) : BC mass flux 375c SclrF (npro) : BC Scalar flux 376c 377c Zdenek Johan, Summer 1990. (Modified from e2bvar.f) 378c Zdenek Johan, Winter 1991. (Fortran 90) 379c---------------------------------------------------------------------- 380c 381 include "common.h" 382c 383 dimension ycl(npro,nshl,ndof), BCB(npro,nshlb,ndBCB), 384 & shpb(npro,nshl), shglb(npro,nsd,nshl), 385 & xlb(npro,nshl,nsd), 386 & lnode(27), 387 & g1yti(npro), g2yti(npro), 388 & g3yti(npro), 389 & WdetJb(npro), bnorm(npro,nsd), 390 & pres(npro), T(npro), 391 & u1(npro), u2(npro), 392 & u3(npro), rho(npro), 393 & ei(npro), cp(npro), 394 & rk(npro), Sclr(npro), 395 & rou(npro), 396 & SclrF(npro) 397c 398 dimension dxdxib(npro,nsd,nsd), 399 & dxidxb(npro,nsd,nsd), temp(npro), 400 & temp1(npro), temp2(npro), 401 & temp3(npro), gl1yti(npro), 402 & gl2yti(npro), gl3yti(npro) 403c 404 dimension h(npro), cv(npro), 405 & alfap(npro), betaT(npro), 406 & gamb(npro), c(npro), 407 & tmp(npro), v1(npro,nsd), 408 & v2(npro,nsd) 409c 410c.... -------------------> integration variables <-------------------- 411c 412c.... compute the primitive variables at the integration point 413c 414 pres = zero 415 u1 = zero 416 u2 = zero 417 u3 = zero 418 T = zero 419 Sclr = zero 420 421 id = isclr+5 422 ibb = isclr+6 423c 424 do n = 1, nshlb 425 nodlcl = lnode(n) 426c 427 pres = pres + shpb(:,nodlcl) * ycl(:,nodlcl,1) 428 u1 = u1 + shpb(:,nodlcl) * ycl(:,nodlcl,2) 429 u2 = u2 + shpb(:,nodlcl) * ycl(:,nodlcl,3) 430 u3 = u3 + shpb(:,nodlcl) * ycl(:,nodlcl,4) 431 T = T + shpb(:,nodlcl) * ycl(:,nodlcl,5) 432 Sclr = Sclr + shpb(:,nodlcl) * ycl(:,nodlcl,id) 433 enddo 434c 435c.... calculate the specific kinetic energy 436c 437 rk = pt5 * ( u1**2 + u2**2 + u3**2 ) 438c 439c.... get the thermodynamic properties 440c 441 ithm = 6 442 if (Navier .eq. 1) ithm = 7 443 call getthm (pres, T, Sclr, 444 & rk, rho, ei, 445 & h, tmp, cv, 446 & cp, alfap, betaT, 447 & gamb, c) 448c 449 if (iconvsclr.eq.2) rho=one 450c 451c.... ----------------------> Element Metrics <----------------------- 452c 453c.... compute the deformation gradient 454c 455 dxdxib = zero 456c 457 do n = 1, nshl 458 dxdxib(:,1,1) = dxdxib(:,1,1) + xlb(:,n,1) * shglb(:,1,n) 459 dxdxib(:,1,2) = dxdxib(:,1,2) + xlb(:,n,1) * shglb(:,2,n) 460 dxdxib(:,1,3) = dxdxib(:,1,3) + xlb(:,n,1) * shglb(:,3,n) 461 dxdxib(:,2,1) = dxdxib(:,2,1) + xlb(:,n,2) * shglb(:,1,n) 462 dxdxib(:,2,2) = dxdxib(:,2,2) + xlb(:,n,2) * shglb(:,2,n) 463 dxdxib(:,2,3) = dxdxib(:,2,3) + xlb(:,n,2) * shglb(:,3,n) 464 dxdxib(:,3,1) = dxdxib(:,3,1) + xlb(:,n,3) * shglb(:,1,n) 465 dxdxib(:,3,2) = dxdxib(:,3,2) + xlb(:,n,3) * shglb(:,2,n) 466 dxdxib(:,3,3) = dxdxib(:,3,3) + xlb(:,n,3) * shglb(:,3,n) 467 enddo 468c 469c.... compute the normal to the boundary 470c 471c 472 v1 = xlb(:,2,:) - xlb(:,1,:) 473 v2 = xlb(:,3,:) - xlb(:,1,:) 474c 475c.... compute the normal to the boundary. This is achieved by taking 476c the cross product of two vectors in the plane of the 2-d 477c boundary face. 478c 479c.....The following are done in order to correct temp1..3 480c based on the results from compressible code. This is done only 481c for wedges, depending on the bounary face.(tri or quad) 482c 483 if (lcsyst .eq. 4) then 484 temp1 = dxdxib(:,2,1) * dxdxib(:,3,3) - 485 & dxdxib(:,2,3) * dxdxib(:,3,1) 486 temp2 = dxdxib(:,3,1) * dxdxib(:,1,3) - 487 & dxdxib(:,3,3) * dxdxib(:,1,1) 488 temp3 = dxdxib(:,1,1) * dxdxib(:,2,3) - 489 & dxdxib(:,1,3) * dxdxib(:,2,1) 490 491 elseif (lcsyst .eq. 1) then 492 temp1 = v1(:,2) * v2(:,3) - v2(:,2) * v1(:,3) 493 temp2 = v2(:,1) * v1(:,3) - v1(:,1) * v2(:,3) 494 temp3 = v1(:,1) * v2(:,2) - v2(:,1) * v1(:,2) 495 else 496 temp1 = - v1(:,2) * v2(:,3) + v2(:,2) * v1(:,3) 497 temp2 = - v2(:,1) * v1(:,3) + v1(:,1) * v2(:,3) 498 temp3 = - v1(:,1) * v2(:,2) + v2(:,1) * v1(:,2) 499 endif 500c 501 temp = one / sqrt ( temp1**2 + temp2**2 + temp3**2 ) 502 bnorm(:,1) = temp1 * temp 503 bnorm(:,2) = temp2 * temp 504 bnorm(:,3) = temp3 * temp 505c 506 507 if (lcsyst .eq. 3) then 508 WdetJb = (1 - Qwtb(lcsyst,intp)) / (four*temp) 509 elseif (lcsyst .eq. 4) then 510 WdetJb = Qwtb(lcsyst,intp)/ temp 511 else 512 WdetJb =Qwtb(lcsyst,intp) / (four*temp) 513 endif 514c 515c.... --------------------------> Grad-V <---------------------------- 516c 517c.... compute grad-v for Navier-Stokes terms 518c 519 if (Navier .eq. 1) then 520c 521c.... compute the inverse of deformation gradient 522c 523 dxidxb(:,1,1) = dxdxib(:,2,2) * dxdxib(:,3,3) 524 & - dxdxib(:,3,2) * dxdxib(:,2,3) 525 dxidxb(:,1,2) = dxdxib(:,3,2) * dxdxib(:,1,3) 526 & - dxdxib(:,1,2) * dxdxib(:,3,3) 527 dxidxb(:,1,3) = dxdxib(:,1,2) * dxdxib(:,2,3) 528 & - dxdxib(:,1,3) * dxdxib(:,2,2) 529 temp = one / ( dxidxb(:,1,1) * dxdxib(:,1,1) 530 & + dxidxb(:,1,2) * dxdxib(:,2,1) 531 & + dxidxb(:,1,3) * dxdxib(:,3,1) ) 532 dxidxb(:,1,1) = dxidxb(:,1,1) * temp 533 dxidxb(:,1,2) = dxidxb(:,1,2) * temp 534 dxidxb(:,1,3) = dxidxb(:,1,3) * temp 535 dxidxb(:,2,1) = (dxdxib(:,2,3) * dxdxib(:,3,1) 536 & - dxdxib(:,2,1) * dxdxib(:,3,3)) * temp 537 dxidxb(:,2,2) = (dxdxib(:,1,1) * dxdxib(:,3,3) 538 & - dxdxib(:,3,1) * dxdxib(:,1,3)) * temp 539 dxidxb(:,2,3) = (dxdxib(:,2,1) * dxdxib(:,1,3) 540 & - dxdxib(:,1,1) * dxdxib(:,2,3)) * temp 541 dxidxb(:,3,1) = (dxdxib(:,2,1) * dxdxib(:,3,2) 542 & - dxdxib(:,2,2) * dxdxib(:,3,1)) * temp 543 dxidxb(:,3,2) = (dxdxib(:,3,1) * dxdxib(:,1,2) 544 & - dxdxib(:,1,1) * dxdxib(:,3,2)) * temp 545 dxidxb(:,3,3) = (dxdxib(:,1,1) * dxdxib(:,2,2) 546 & - dxdxib(:,1,2) * dxdxib(:,2,1)) * temp 547 548c 549c.... compute local-grad-Y 550c 551 gl1yti = zero 552 gl2yti = zero 553 gl3yti = zero 554c 555 do n = 1, nshl 556 gl1yti(:) = gl1yti(:) + shglb(:,1,n) * ycl(:,n,id) 557 gl2yti(:) = gl2yti(:) + shglb(:,2,n) * ycl(:,n,id) 558 gl3yti(:) = gl3yti(:) + shglb(:,3,n) * ycl(:,n,id) 559 enddo 560c 561c.... convert local-grads to global-grads 562c 563 g1yti(:) = dxidxb(:,1,1) * gl1yti(:) + 564 & dxidxb(:,2,1) * gl2yti(:) + 565 & dxidxb(:,3,1) * gl3yti(:) 566 g2yti(:) = dxidxb(:,1,2) * gl1yti(:) + 567 & dxidxb(:,2,2) * gl2yti(:) + 568 & dxidxb(:,3,2) * gl3yti(:) 569 g3yti(:) = dxidxb(:,1,3) * gl1yti(:) + 570 & dxidxb(:,2,3) * gl2yti(:) + 571 & dxidxb(:,3,3) * gl3yti(:) 572 573c 574c.... end grad-Sclr 575 endif 576c 577c.... --------------------> Boundary Conditions <-------------------- 578c 579c.... compute the Euler boundary conditions 580c 581 rou = zero 582 do n = 1, nshlb 583 nodlcl = lnode(n) 584 rou = rou + shpb(:,nodlcl) * BCB(:,n,1) 585 enddo 586c 587c.... impose scalar flux boundary conditions 588 SclrF = zero 589 do n=1,nshlb 590 nodlcl = lnode(n) 591 SclrF = SclrF + shpb(:,nodlcl) * BCB(:,n,ibb) 592 enddo 593 594c 595c.... flop count 596c 597 flops = flops + (27+18*nshl+14*nenbl)*npro 598c 599c.... return 600c 601 return 602 end 603 604