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