1 subroutine e3bvar (yl, acl, ul, 2 & shpb, shglb, 3 & xlb, lnode, 4 & WdetJb, bnorm, pres, 5 & u1, u2, u3, rmu, 6 & unm, tau1n, tau2n, tau3n, 7 & vdot, rlKwall, 8 & xKebe, rKwall_glob) 9c 10c---------------------------------------------------------------------- 11c 12c This routine computes the variables at integration points for 13c the boundary element routine. 14c 15c input: 16c yl (npro,nshl,ndof) : primitive variables (local) 17c ndof: 5[p,v1,v2,v3,T]+number of scalars solved 18c acl (npro,nshl,ndof) : acceleration (local) 19c ul (npro,nshlb,nsd) : displacement (local) 20c shpb (nen) : boundary element shape-functions 21c shglb (nsd,nen) : boundary element grad-shape-functions 22c xlb (npro,nenl,nsd) : nodal coordinates at current step 23c lnode (nenb) : local nodes on the boundary 24c 25c output: 26c g1yi (npro,ndof) : grad-v in direction 1 27c g2yi (npro,ndof) : grad-v in direction 2 28c g3yi (npro,ndof) : grad-v in direction 3 29c WdetJb (npro) : weighted Jacobian 30c bnorm (npro,nsd) : outward normal 31c pres (npro) : pressure 32c u1 (npro) : x1-velocity component 33c u2 (npro) : x2-velocity component 34c u3 (npro) : x3-velocity component 35c unm (npro) : BC u dot n 36c p (npro) : BC pressure 37c tau1n (npro) : BC viscous flux 1 38c tau2n (npro) : BC viscous flux 2 39c tau3n (npro) : BC viscous flux 3 40c vdot (npro,nsd) : acceleration at quadrature points 41c rlKwall(npro,nshlb,nsd) : wall stiffness contribution to the local residual 42c 43c Zdenek Johan, Summer 1990. (Modified from e2bvar.f) 44c Zdenek Johan, Winter 1991. (Fortran 90) 45c Alberto Figueroa, Winter 2004. CMM-FSI 46c---------------------------------------------------------------------- 47c 48 use turbsa 49 include "common.h" 50c 51 dimension yl(npro,nshl,ndof), rmu(npro), 52 & shpb(npro,nshl), shglb(npro,nsd,nshl), 53 & xlb(npro,nenl,nsd), 54 & lnode(27), g1yi(npro,ndof), 55 & g2yi(npro,ndof), g3yi(npro,ndof), 56 & WdetJb(npro), bnorm(npro,nsd), 57 & pres(npro), 58 & u1(npro), u2(npro), 59 & u3(npro), 60 & unm(npro), 61 & tau1n(npro), tau2n(npro), 62 & tau3n(npro), 63 & acl(npro,nshl,ndof), ul(npro,nshl,nsd), 64 & vdot(npro,nsd), rlKwall(npro,nshlb,nsd) 65c 66 dimension gl1yi(npro,ndof), gl2yi(npro,ndof), 67 & gl3yi(npro,ndof), dxdxib(npro,nsd,nsd), 68 & dxidxb(npro,nsd,nsd), temp(npro), 69 & temp1(npro), temp2(npro), 70 & temp3(npro), 71 & v1(npro,nsd), v2(npro,nsd), 72 & v3(npro,nsd), 73 & rotnodallocal(npro,nsd,nsd), 74 & x1rot(npro,nsd), x2rot(npro,nsd), 75 & x3rot(npro,nsd), detJacrot(npro), 76 & B1(npro,5,3), B2(npro,5,3), 77 & B3(npro,5,3), Dmatrix(npro,5,5), 78 & DtimesB1(npro,5,3), DtimesB2(npro,5,3), 79 & DtimesB3(npro,5,3), 80 & rKwall_local11(npro,nsd,nsd), 81 & rKwall_local12(npro,nsd,nsd), 82 & rKwall_local13(npro,nsd,nsd), 83 & rKwall_local21(npro,nsd,nsd), 84 & rKwall_local22(npro,nsd,nsd), 85 & rKwall_local23(npro,nsd,nsd), 86 & rKwall_local31(npro,nsd,nsd), 87 & rKwall_local32(npro,nsd,nsd), 88 & rKwall_local33(npro,nsd,nsd), 89 & rKwall_glob11(npro,nsd,nsd), 90 & rKwall_glob12(npro,nsd,nsd), 91 & rKwall_glob13(npro,nsd,nsd), 92 & rKwall_glob21(npro,nsd,nsd), 93 & rKwall_glob22(npro,nsd,nsd), 94 & rKwall_glob23(npro,nsd,nsd), 95 & rKwall_glob31(npro,nsd,nsd), 96 & rKwall_glob32(npro,nsd,nsd), 97 & rKwall_glob33(npro,nsd,nsd) 98c 99 dimension rKwall_glob(npro,9,nshl,nshl), 100 & xKebe(npro,9,nshl,nshl) 101c 102 real*8 lhmFctvw, tsFctvw(npro) 103 104 dimension tmp1(npro) 105c 106 real*8 Turb(npro), xki, 107 & xki3, fv1 108c 109 integer e, i, j 110c 111 integer aa, b 112 113 114c 115c.... -------------------> integration variables <-------------------- 116c 117c.... compute the primitive variables at the integration point 118c 119 pres = zero 120 u1 = zero 121 u2 = zero 122 u3 = zero 123c 124 125 do n = 1, nshlb 126 nodlcl = lnode(n) 127c 128 pres = pres + shpb(:,nodlcl) * yl(:,nodlcl,1) 129 u1 = u1 + shpb(:,nodlcl) * yl(:,nodlcl,2) 130 u2 = u2 + shpb(:,nodlcl) * yl(:,nodlcl,3) 131 u3 = u3 + shpb(:,nodlcl) * yl(:,nodlcl,4) 132 133 enddo 134c 135c.... ----------------------> Element Metrics <----------------------- 136c 137c.... compute the deformation gradient 138c 139 dxdxib = zero 140c 141 do n = 1, nenl 142 dxdxib(:,1,1) = dxdxib(:,1,1) + xlb(:,n,1) * shglb(:,1,n) 143 dxdxib(:,1,2) = dxdxib(:,1,2) + xlb(:,n,1) * shglb(:,2,n) 144 dxdxib(:,1,3) = dxdxib(:,1,3) + xlb(:,n,1) * shglb(:,3,n) 145 dxdxib(:,2,1) = dxdxib(:,2,1) + xlb(:,n,2) * shglb(:,1,n) 146 dxdxib(:,2,2) = dxdxib(:,2,2) + xlb(:,n,2) * shglb(:,2,n) 147 dxdxib(:,2,3) = dxdxib(:,2,3) + xlb(:,n,2) * shglb(:,3,n) 148 dxdxib(:,3,1) = dxdxib(:,3,1) + xlb(:,n,3) * shglb(:,1,n) 149 dxdxib(:,3,2) = dxdxib(:,3,2) + xlb(:,n,3) * shglb(:,2,n) 150 dxdxib(:,3,3) = dxdxib(:,3,3) + xlb(:,n,3) * shglb(:,3,n) 151 enddo 152c 153c.... compute the normal to the boundary 154c 155c$$$ if (lcsyst .eq. 4) then ! wedge-quad 156c$$$ temp1 = dxdxib(:,2,1) * dxdxib(:,3,3) - 157c$$$ & dxdxib(:,2,3) * dxdxib(:,3,1) 158c$$$ temp2 = dxdxib(:,3,1) * dxdxib(:,1,3) - 159c$$$ & dxdxib(:,3,3) * dxdxib(:,1,1) 160c$$$ temp3 = dxdxib(:,1,1) * dxdxib(:,2,3) - 161c$$$ & dxdxib(:,1,3) * dxdxib(:,2,1) 162c$$$ elseif( lcyst .eq. 6) then ! pyr-tri face 163c$$$ temp1 = dxdxib(:,2,1) * dxdxib(:,3,3) - 164c$$$ & dxdxib(:,2,3) * dxdxib(:,3,1) 165c$$$ temp2 = dxdxib(:,3,1) * dxdxib(:,1,3) - 166c$$$ & dxdxib(:,3,3) * dxdxib(:,1,1) 167c$$$ temp3 = dxdxib(:,1,1) * dxdxib(:,2,3) - 168c$$$ & dxdxib(:,1,3) * dxdxib(:,2,1) 169c$$$ elseif( lcyst .eq. 1) then !usual wrong way tets 170c$$$ temp1 = -dxdxib(:,2,2) * dxdxib(:,3,1) + 171c$$$ & dxdxib(:,2,1) * dxdxib(:,3,2) 172c$$$ temp2 = -dxdxib(:,3,2) * dxdxib(:,1,1) + 173c$$$ & dxdxib(:,3,1) * dxdxib(:,1,2) 174c$$$ temp3 = -dxdxib(:,1,2) * dxdxib(:,2,1) + 175c$$$ & dxdxib(:,1,1) * dxdxib(:,2,2) 176c$$$ else 177c$$$ temp1 = dxdxib(:,2,2) * dxdxib(:,3,1) - 178c$$$ & dxdxib(:,2,1) * dxdxib(:,3,2) 179c$$$ temp2 = dxdxib(:,3,2) * dxdxib(:,1,1) - 180c$$$ & dxdxib(:,3,1) * dxdxib(:,1,2) 181c$$$ temp3 = dxdxib(:,1,2) * dxdxib(:,2,1) - 182c$$$ & dxdxib(:,1,1) * dxdxib(:,2,2) 183c$$$ endif 184c$$$c 185c$$$ temp = one / sqrt ( temp1**2 + temp2**2 + temp3**2 ) 186c$$$ bnorm(:,1) = temp1 * temp 187c$$$ bnorm(:,2) = temp2 * temp 188c$$$ bnorm(:,3) = temp3 * temp 189c$$$c 190c$$$ WdetJb = Qwtb(lcsyst,intp) / temp 191c$$$ if(lcsyst .eq. 3) WdetJb=WdetJb*two 192c 193c.... compute the normal to the boundary. This is achieved by taking 194c the cross product of two vectors in the plane of the 2-d 195c boundary face. 196c 197 if(lcsyst.eq.1) then ! set to curl into element all others out 198 ipt2=2 199 ipt3=3 200 elseif(lcsyst.eq.2) then 201 ipt2=4 202 ipt3=2 203 elseif(lcsyst.eq.3) then 204 ipt2=3 205 ipt3=2 206 elseif(lcsyst.eq.4) then 207 ipt2=2 208 ipt3=4 209 elseif(lcsyst.eq.5) then 210 ipt2=4 211 ipt3=2 212 elseif(lcsyst.eq.6) then 213 ipt2=2 214 ipt3=5 215 endif 216 v1 = xlb(:,ipt2,:) - xlb(:,1,:) 217 v2 = xlb(:,ipt3,:) - xlb(:,1,:) 218c 219c compute cross product 220c 221 temp1 = v1(:,2) * v2(:,3) - v2(:,2) * v1(:,3) 222 temp2 = v2(:,1) * v1(:,3) - v1(:,1) * v2(:,3) 223 temp3 = v1(:,1) * v2(:,2) - v2(:,1) * v1(:,2) 224c 225c mag is area for quads, twice area for tris 226c 227 temp = one / sqrt ( temp1**2 + temp2**2 + temp3**2 ) 228 bnorm(:,1) = temp1 * temp 229 bnorm(:,2) = temp2 * temp 230 bnorm(:,3) = temp3 * temp 231c 232 233 if (lcsyst .eq. 1) then 234 WdetJb = Qwtb(lcsyst,intp) / (four*temp) 235 elseif (lcsyst .eq. 2) then 236 WdetJb = Qwtb(lcsyst,intp) / (four*temp) 237 elseif (lcsyst .eq. 3) then 238 WdetJb = Qwtb(lcsyst,intp) / (two*temp) 239 elseif (lcsyst .eq. 4) then 240 WdetJb = Qwtb(lcsyst,intp) / (four*temp) 241 elseif (lcsyst .eq. 5) then 242 WdetJb = Qwtb(lcsyst,intp) / (four*temp) 243 elseif (lcsyst .eq. 6) then 244 WdetJb = Qwtb(lcsyst,intp) / (two*temp) 245 endif 246c 247c.... --------------------------> Grad-V <---------------------------- 248c 249c.... compute grad-v for Navier-Stokes terms 250c 251 if (Navier .eq. 1) then 252c 253c.... compute the inverse of deformation gradient 254c 255 dxidxb(:,1,1) = dxdxib(:,2,2) * dxdxib(:,3,3) 256 & - dxdxib(:,3,2) * dxdxib(:,2,3) 257 dxidxb(:,1,2) = dxdxib(:,3,2) * dxdxib(:,1,3) 258 & - dxdxib(:,1,2) * dxdxib(:,3,3) 259 dxidxb(:,1,3) = dxdxib(:,1,2) * dxdxib(:,2,3) 260 & - dxdxib(:,1,3) * dxdxib(:,2,2) 261 temp = one / ( dxidxb(:,1,1) * dxdxib(:,1,1) 262 & + dxidxb(:,1,2) * dxdxib(:,2,1) 263 & + dxidxb(:,1,3) * dxdxib(:,3,1) ) 264 dxidxb(:,1,1) = dxidxb(:,1,1) * temp 265 dxidxb(:,1,2) = dxidxb(:,1,2) * temp 266 dxidxb(:,1,3) = dxidxb(:,1,3) * temp 267 dxidxb(:,2,1) = (dxdxib(:,2,3) * dxdxib(:,3,1) 268 & - dxdxib(:,2,1) * dxdxib(:,3,3)) * temp 269 dxidxb(:,2,2) = (dxdxib(:,1,1) * dxdxib(:,3,3) 270 & - dxdxib(:,3,1) * dxdxib(:,1,3)) * temp 271 dxidxb(:,2,3) = (dxdxib(:,2,1) * dxdxib(:,1,3) 272 & - dxdxib(:,1,1) * dxdxib(:,2,3)) * temp 273 dxidxb(:,3,1) = (dxdxib(:,2,1) * dxdxib(:,3,2) 274 & - dxdxib(:,2,2) * dxdxib(:,3,1)) * temp 275 dxidxb(:,3,2) = (dxdxib(:,3,1) * dxdxib(:,1,2) 276 & - dxdxib(:,1,1) * dxdxib(:,3,2)) * temp 277 dxidxb(:,3,3) = (dxdxib(:,1,1) * dxdxib(:,2,2) 278 & - dxdxib(:,1,2) * dxdxib(:,2,1)) * temp 279c 280c.... compute local-grad-Y 281c 282 gl1yi = zero 283 gl2yi = zero 284 gl3yi = zero 285c 286 do n = 1, nshl 287 gl1yi(:,1) = gl1yi(:,1) + shglb(:,1,n) * yl(:,n,1) 288 gl1yi(:,2) = gl1yi(:,2) + shglb(:,1,n) * yl(:,n,2) 289 gl1yi(:,3) = gl1yi(:,3) + shglb(:,1,n) * yl(:,n,3) 290 gl1yi(:,4) = gl1yi(:,4) + shglb(:,1,n) * yl(:,n,4) 291c 292 gl2yi(:,1) = gl2yi(:,1) + shglb(:,2,n) * yl(:,n,1) 293 gl2yi(:,2) = gl2yi(:,2) + shglb(:,2,n) * yl(:,n,2) 294 gl2yi(:,3) = gl2yi(:,3) + shglb(:,2,n) * yl(:,n,3) 295 gl2yi(:,4) = gl2yi(:,4) + shglb(:,2,n) * yl(:,n,4) 296c 297 gl3yi(:,1) = gl3yi(:,1) + shglb(:,3,n) * yl(:,n,1) 298 gl3yi(:,2) = gl3yi(:,2) + shglb(:,3,n) * yl(:,n,2) 299 gl3yi(:,3) = gl3yi(:,3) + shglb(:,3,n) * yl(:,n,3) 300 gl3yi(:,4) = gl3yi(:,4) + shglb(:,3,n) * yl(:,n,4) 301 enddo 302c 303c.... convert local-grads to global-grads 304c 305 g1yi(:,2) = dxidxb(:,1,1) * gl1yi(:,2) + 306 & dxidxb(:,2,1) * gl2yi(:,2) + 307 & dxidxb(:,3,1) * gl3yi(:,2) 308 g2yi(:,2) = dxidxb(:,1,2) * gl1yi(:,2) + 309 & dxidxb(:,2,2) * gl2yi(:,2) + 310 & dxidxb(:,3,2) * gl3yi(:,2) 311 g3yi(:,2) = dxidxb(:,1,3) * gl1yi(:,2) + 312 & dxidxb(:,2,3) * gl2yi(:,2) + 313 & dxidxb(:,3,3) * gl3yi(:,2) 314c 315 g1yi(:,3) = dxidxb(:,1,1) * gl1yi(:,3) + 316 & dxidxb(:,2,1) * gl2yi(:,3) + 317 & dxidxb(:,3,1) * gl3yi(:,3) 318 g2yi(:,3) = dxidxb(:,1,2) * gl1yi(:,3) + 319 & dxidxb(:,2,2) * gl2yi(:,3) + 320 & dxidxb(:,3,2) * gl3yi(:,3) 321 g3yi(:,3) = dxidxb(:,1,3) * gl1yi(:,3) + 322 & dxidxb(:,2,3) * gl2yi(:,3) + 323 & dxidxb(:,3,3) * gl3yi(:,3) 324c 325 g1yi(:,4) = dxidxb(:,1,1) * gl1yi(:,4) + 326 & dxidxb(:,2,1) * gl2yi(:,4) + 327 & dxidxb(:,3,1) * gl3yi(:,4) 328 g2yi(:,4) = dxidxb(:,1,2) * gl1yi(:,4) + 329 & dxidxb(:,2,2) * gl2yi(:,4) + 330 & dxidxb(:,3,2) * gl3yi(:,4) 331 g3yi(:,4) = dxidxb(:,1,3) * gl1yi(:,4) + 332 & dxidxb(:,2,3) * gl2yi(:,4) + 333 & dxidxb(:,3,3) * gl3yi(:,4) 334c 335c.... end grad-v 336c 337 endif 338 339c 340c.... mass flux 341c 342 unm = bnorm(:,1) * u1 +bnorm(:,2) * u2 +bnorm(:,3) * u3 343! no rho in continuity eq. 344 345 346c 347c.... viscous flux 348c 349 tau1n = bnorm(:,1) * two * rmu * g1yi(:,2) 350 & + bnorm(:,2) * (rmu * (g2yi(:,2) + g1yi(:,3))) 351 & + bnorm(:,3) * (rmu * (g3yi(:,2) + g1yi(:,4))) 352 tau2n = bnorm(:,1) * (rmu * (g2yi(:,2) + g1yi(:,3))) 353 & + bnorm(:,2) * two * rmu * g2yi(:,3) 354 & + bnorm(:,3) * (rmu * (g3yi(:,3) + g2yi(:,4))) 355 tau3n = bnorm(:,1) * (rmu * (g3yi(:,2) + g1yi(:,4))) 356 & + bnorm(:,2) * (rmu * (g3yi(:,3) + g2yi(:,4))) 357 & + bnorm(:,3) * two * rmu * g3yi(:,4) 358c 359 temp1 = bnorm(:,1) * tau1n 360 & + bnorm(:,2) * tau2n 361 & + bnorm(:,3) * tau3n 362 363 pres = pres - temp1 364 365 tau1n = tau1n - bnorm(:,1) * temp1 366 tau2n = tau2n - bnorm(:,2) * temp1 367 tau3n = tau3n - bnorm(:,3) * temp1 368 369c 370c.... viscous flux control 371c 372c if iviscflux = 1, we consider the viscous flux on the RHS 373c otherwise, if iviscflux = 0, we eliminate this term: stability in 374c pressure-flow coupled boundaries for cardiovascular applications in 375c situations of flow reversal 376 tau1n = tau1n * iviscflux 377 tau2n = tau2n * iviscflux 378 tau3n = tau3n * iviscflux 379 380 vdot = zero 381 rlKwall = zero 382 if (intp.eq.ngaussb) then ! do this only for the last gauss point 383 rKwall_glob = zero 384 endif 385 386 if(ideformwall.eq.1) then 387 do n = 1, nshlb 388 nodlcl = lnode(n) 389c 390 vdot(:,1) = vdot(:,1) + shpb(:,nodlcl) * acl(:,nodlcl,2) 391 vdot(:,2) = vdot(:,2) + shpb(:,nodlcl) * acl(:,nodlcl,3) 392 vdot(:,3) = vdot(:,3) + shpb(:,nodlcl) * acl(:,nodlcl,4) 393 394 enddo 395 vdot = vdot * thicknessvw * rhovw 396c 397c.... ---------------------> Stiffness matrix & residual <----------------- 398c 399c.... B^t * D * B formulation for plane stress enhanced membrane 400c 401c 402c.... rotation matrix 403c 404 v1 = xlb(:,ipt2,:) - xlb(:,1,:) 405 temp = one / sqrt ( v1(:,1)**2 + v1(:,2)**2 + v1(:,3)**2 ) 406 v1(:,1) = v1(:,1) * temp 407 v1(:,2) = v1(:,2) * temp 408 v1(:,3) = v1(:,3) * temp 409 410 v2 = xlb(:,ipt3,:) - xlb(:,1,:) 411 412c compute cross product 413 temp1 = v1(:,2) * v2(:,3) - v2(:,2) * v1(:,3) 414 temp2 = v2(:,1) * v1(:,3) - v1(:,1) * v2(:,3) 415 temp3 = v1(:,1) * v2(:,2) - v2(:,1) * v1(:,2) 416 417 temp = one / sqrt ( temp1**2 + temp2**2 + temp3**2 ) 418 v3(:,1) = temp1 * temp 419 v3(:,2) = temp2 * temp 420 v3(:,3) = temp3 * temp 421 422c cross product again for v2 423 temp1 = v3(:,2) * v1(:,3) - v1(:,2) * v3(:,3) 424 temp2 = v1(:,1) * v3(:,3) - v3(:,1) * v1(:,3) 425 temp3 = v3(:,1) * v1(:,2) - v1(:,1) * v3(:,2) 426 427 temp = one / sqrt ( temp1**2 + temp2**2 + temp3**2 ) 428 v2(:,1) = temp1 * temp 429 v2(:,2) = temp2 * temp 430 v2(:,3) = temp3 * temp 431 432 do j = 1, nsd 433 rotnodallocal(:,1,j) = v1(:,j) 434 rotnodallocal(:,2,j) = v2(:,j) 435 rotnodallocal(:,3,j) = v3(:,j) 436 enddo 437 438c 439c.... rotated coordinates 440c 441 x1rot = zero 442 x2rot = zero 443 x3rot = zero 444 445 do i = 1, nsd 446 do j = 1, nsd 447 x1rot(:,i) = x1rot(:,i)+rotnodallocal(:,i,j)*xlb(:,1,j) 448 x2rot(:,i) = x2rot(:,i)+rotnodallocal(:,i,j)*xlb(:,ipt2,j) 449 x3rot(:,i) = x3rot(:,i)+rotnodallocal(:,i,j)*xlb(:,ipt3,j) 450 enddo 451 enddo 452 453c 454c.... B matrices 455c 456 B1 = zero 457 B2 = zero 458 B3 = zero 459 detJacrot = (x2rot(:,1)-x1rot(:,1)) * (x3rot(:,2)-x1rot(:,2)) - 460 & (x3rot(:,1)-x1rot(:,1)) * (x2rot(:,2)-x1rot(:,2)) 461 462 B1(:,1,1) = (x2rot(:,2)-x3rot(:,2))/detJacrot(:) 463 B1(:,2,2) = (x3rot(:,1)-x2rot(:,1))/detJacrot(:) 464 B1(:,3,1) = (x3rot(:,1)-x2rot(:,1))/detJacrot(:) 465 B1(:,3,2) = (x2rot(:,2)-x3rot(:,2))/detJacrot(:) 466 B1(:,4,3) = (x2rot(:,2)-x3rot(:,2))/detJacrot(:) 467 B1(:,5,3) = (x3rot(:,1)-x2rot(:,1))/detJacrot(:) 468 469 B2(:,1,1) = (x3rot(:,2)-x1rot(:,2))/detJacrot(:) 470 B2(:,2,2) = (x1rot(:,1)-x3rot(:,1))/detJacrot(:) 471 B2(:,3,1) = (x1rot(:,1)-x3rot(:,1))/detJacrot(:) 472 B2(:,3,2) = (x3rot(:,2)-x1rot(:,2))/detJacrot(:) 473 B2(:,4,3) = (x3rot(:,2)-x1rot(:,2))/detJacrot(:) 474 B2(:,5,3) = (x1rot(:,1)-x3rot(:,1))/detJacrot(:) 475 476 B3(:,1,1) = (x1rot(:,2)-x2rot(:,2))/detJacrot(:) 477 B3(:,2,2) = (x2rot(:,1)-x1rot(:,1))/detJacrot(:) 478 B3(:,3,1) = (x2rot(:,1)-x1rot(:,1))/detJacrot(:) 479 B3(:,3,2) = (x1rot(:,2)-x2rot(:,2))/detJacrot(:) 480 B3(:,4,3) = (x1rot(:,2)-x2rot(:,2))/detJacrot(:) 481 B3(:,5,3) = (x2rot(:,1)-x1rot(:,1))/detJacrot(:) 482 483C B1 = B1 / detJacrot 484C B2 = B2 / detJacrot 485C B3 = B3 / detJacrot 486 487c 488c.... D matrix 489c 490 Dmatrix = zero 491 temp1 = evw / (1.0d0 - rnuvw*rnuvw) 492 temp2 = rnuvw * temp1 493 temp3 = pt5 * (1.0d0 - rnuvw) * temp1 494 Dmatrix(:,1,1) = temp1 495 Dmatrix(:,1,2) = temp2 496 Dmatrix(:,2,1) = temp2 497 Dmatrix(:,2,2) = temp1 498 Dmatrix(:,3,3) = temp3 499 Dmatrix(:,4,4) = temp3*rshearconstantvw 500 Dmatrix(:,5,5) = temp3*rshearconstantvw 501c 502c.... D * [B1|B2|B3] 503c 504 DtimesB1 = zero 505 DtimesB2 = zero 506 DtimesB3 = zero 507 do i = 1, 5 508 do j = 1, 3 509 do k = 1, 5 510 DtimesB1(:,i,j) = DtimesB1(:,i,j) 511 & + Dmatrix(:,i,k) * B1(:,k,j) 512 DtimesB2(:,i,j) = DtimesB2(:,i,j) 513 & + Dmatrix(:,i,k) * B2(:,k,j) 514 DtimesB3(:,i,j) = DtimesB3(:,i,j) 515 & + Dmatrix(:,i,k) * B3(:,k,j) 516 enddo 517 enddo 518 enddo 519c 520c.... [B1|B2|B3]^T * D * [B1|B2|B3] 521c 522 rKwall_local11 = zero 523 rKwall_local12 = zero 524 rKwall_local13 = zero 525 rKwall_local21 = zero 526 rKwall_local22 = zero 527 rKwall_local23 = zero 528 rKwall_local31 = zero 529 rKwall_local32 = zero 530 rKwall_local33 = zero 531 532 do i = 1, 3 ! i is a node index: i=1, nenbl=3 533 do j = 1, 3 ! same is true for j 534 do k = 1, 5 535 rKwall_local11(:,i,j) = rKwall_local11(:,i,j) 536 & + B1(:,k,i) * DtimesB1(:,k,j) 537 rKwall_local12(:,i,j) = rKwall_local12(:,i,j) 538 & + B1(:,k,i) * DtimesB2(:,k,j) 539 rKwall_local13(:,i,j) = rKwall_local13(:,i,j) 540 & + B1(:,k,i) * DtimesB3(:,k,j) 541 rKwall_local21(:,i,j) = rKwall_local21(:,i,j) 542 & + B2(:,k,i) * DtimesB1(:,k,j) 543 rKwall_local22(:,i,j) = rKwall_local22(:,i,j) 544 & + B2(:,k,i) * DtimesB2(:,k,j) 545 rKwall_local23(:,i,j) = rKwall_local23(:,i,j) 546 & + B2(:,k,i) * DtimesB3(:,k,j) 547 rKwall_local31(:,i,j) = rKwall_local31(:,i,j) 548 & + B3(:,k,i) * DtimesB1(:,k,j) 549 rKwall_local32(:,i,j) = rKwall_local32(:,i,j) 550 & + B3(:,k,i) * DtimesB2(:,k,j) 551 rKwall_local33(:,i,j) = rKwall_local33(:,i,j) 552 & + B3(:,k,i) * DtimesB3(:,k,j) 553 enddo 554 enddo 555 enddo 556 557c 558c.... Now we need to rotate each of these submatrices to the global frame 559c 560 call rotatestiff(rKwall_local11, rotnodallocal, rKwall_glob11) 561 call rotatestiff(rKwall_local12, rotnodallocal, rKwall_glob12) 562 call rotatestiff(rKwall_local13, rotnodallocal, rKwall_glob13) 563 call rotatestiff(rKwall_local21, rotnodallocal, rKwall_glob21) 564 call rotatestiff(rKwall_local22, rotnodallocal, rKwall_glob22) 565 call rotatestiff(rKwall_local23, rotnodallocal, rKwall_glob23) 566 call rotatestiff(rKwall_local31, rotnodallocal, rKwall_glob31) 567 call rotatestiff(rKwall_local32, rotnodallocal, rKwall_glob32) 568 call rotatestiff(rKwall_local33, rotnodallocal, rKwall_glob33) 569 570c multiply the nodal matrices by the area and the thickness 571 do i =1, nsd 572 do j = 1, nsd 573 rKwall_glob11(:,i,j) = rKwall_glob11(:,i,j) * detJacrot(:) 574 & * pt5 * thicknessvw 575 rKwall_glob12(:,i,j) = rKwall_glob12(:,i,j) * detJacrot(:) 576 & * pt5 * thicknessvw 577 rKwall_glob13(:,i,j) = rKwall_glob13(:,i,j) * detJacrot(:) 578 & * pt5 * thicknessvw 579 rKwall_glob21(:,i,j) = rKwall_glob21(:,i,j) * detJacrot(:) 580 & * pt5 * thicknessvw 581 rKwall_glob22(:,i,j) = rKwall_glob22(:,i,j) * detJacrot(:) 582 & * pt5 * thicknessvw 583 rKwall_glob23(:,i,j) = rKwall_glob23(:,i,j) * detJacrot(:) 584 & * pt5 * thicknessvw 585 rKwall_glob31(:,i,j) = rKwall_glob31(:,i,j) * detJacrot(:) 586 & * pt5 * thicknessvw 587 rKwall_glob32(:,i,j) = rKwall_glob32(:,i,j) * detJacrot(:) 588 & * pt5 * thicknessvw 589 rKwall_glob33(:,i,j) = rKwall_glob33(:,i,j) * detJacrot(:) 590 & * pt5 * thicknessvw 591 enddo 592 enddo 593 594c 595c.... Final K * u product (in global coordinates) to get the residual 596c 597 do i = 1, 3 ! now i is a spatial index: i=1, nsd=3 598 rlKwall(:,1,1) = rlKwall(:,1,1) 599 & + rKwall_glob11(:,1,i) * ul(:,1,i) 600 & + rKwall_glob12(:,1,i) * ul(:,2,i) 601 & + rKwall_glob13(:,1,i) * ul(:,3,i) 602 rlKwall(:,1,2) = rlKwall(:,1,2) 603 & + rKwall_glob11(:,2,i) * ul(:,1,i) 604 & + rKwall_glob12(:,2,i) * ul(:,2,i) 605 & + rKwall_glob13(:,2,i) * ul(:,3,i) 606 rlKwall(:,1,3) = rlKwall(:,1,3) 607 & + rKwall_glob11(:,3,i) * ul(:,1,i) 608 & + rKwall_glob12(:,3,i) * ul(:,2,i) 609 & + rKwall_glob13(:,3,i) * ul(:,3,i) 610 rlKwall(:,2,1) = rlKwall(:,2,1) 611 & + rKwall_glob21(:,1,i) * ul(:,1,i) 612 & + rKwall_glob22(:,1,i) * ul(:,2,i) 613 & + rKwall_glob23(:,1,i) * ul(:,3,i) 614 rlKwall(:,2,2) = rlKwall(:,2,2) 615 & + rKwall_glob21(:,2,i) * ul(:,1,i) 616 & + rKwall_glob22(:,2,i) * ul(:,2,i) 617 & + rKwall_glob23(:,2,i) * ul(:,3,i) 618 rlKwall(:,2,3) = rlKwall(:,2,3) 619 & + rKwall_glob21(:,3,i) * ul(:,1,i) 620 & + rKwall_glob22(:,3,i) * ul(:,2,i) 621 & + rKwall_glob23(:,3,i) * ul(:,3,i) 622 rlKwall(:,3,1) = rlKwall(:,3,1) 623 & + rKwall_glob31(:,1,i) * ul(:,1,i) 624 & + rKwall_glob32(:,1,i) * ul(:,2,i) 625 & + rKwall_glob33(:,1,i) * ul(:,3,i) 626 rlKwall(:,3,2) = rlKwall(:,3,2) 627 & + rKwall_glob31(:,2,i) * ul(:,1,i) 628 & + rKwall_glob32(:,2,i) * ul(:,2,i) 629 & + rKwall_glob33(:,2,i) * ul(:,3,i) 630 rlKwall(:,3,3) = rlKwall(:,3,3) 631 & + rKwall_glob31(:,3,i) * ul(:,1,i) 632 & + rKwall_glob32(:,3,i) * ul(:,2,i) 633 & + rKwall_glob33(:,3,i) * ul(:,3,i) 634 enddo 635c 636c.... --------------> End of Stiffness matrix & residual <----------------- 637c 638 639c 640c.... -----> Wall Stiffness and Mass matrices for implicit LHS <----------- 641c 642 643c.... Here we just add the mass matrix contribution. The stiffness contribution 644c.... is added in e3b 645 646c lhmFct = almi * (one - flmpl) Maybe we have to define flmplW: lumped 647 ! mass parameter for the wall 648 lhmFctvw = almi * (one - flmpl) 649c 650c.... scale variables for efficiency 651c 652 tsFctvw = lhmFctvw * WdetJb * rhovw * thicknessvw 653c 654c.... compute mass and convection terms 655c 656c.... NOTE: the wall mass contributions should only have 3 nodal components 657c.... since the fourth node is an interior node... therefore, the loops should 658c.... be done from 1 to nshlb=3... 659 660 do b = 1, nshlb 661 do aa = 1, nshlb 662 tmp1 = tsFctvw * shpb(:,aa) * shpb(:,b) 663c 664c tmp1=alpha_m*(1-lmp)*WdetJ*N^aN^b*rho*thickness the time term 665c 666 xKebe(:,1,aa,b) = xKebe(:,1,aa,b) + tmp1 667 xKebe(:,5,aa,b) = xKebe(:,5,aa,b) + tmp1 668 xKebe(:,9,aa,b) = xKebe(:,9,aa,b) + tmp1 669 enddo 670 enddo 671 672c 673c.... assemble the nodal stiffness into the element stiffness matrix rKwall_glob 674c 675c.... We have passed the integer intp to make this operation only once: we are 676c.... not using the gauss points structure to compute the stiffness of the wall 677c.... elements, so we don't want to be redundant and calculate ngaussb times the 678c.... stiffness matrix which is constant for linear triangles... 679 680c.... This is ugly, but I will fix it later... 681 682 if (intp.eq.ngaussb) then ! do this only for the last gauss point 683 rKwall_glob(:,1,1,1) = rKwall_glob11(:,1,1) 684 rKwall_glob(:,2,1,1) = rKwall_glob11(:,1,2) 685 rKwall_glob(:,3,1,1) = rKwall_glob11(:,1,3) 686 rKwall_glob(:,4,1,1) = rKwall_glob11(:,2,1) 687 rKwall_glob(:,5,1,1) = rKwall_glob11(:,2,2) 688 rKwall_glob(:,6,1,1) = rKwall_glob11(:,2,3) 689 rKwall_glob(:,7,1,1) = rKwall_glob11(:,3,1) 690 rKwall_glob(:,8,1,1) = rKwall_glob11(:,3,2) 691 rKwall_glob(:,9,1,1) = rKwall_glob11(:,3,3) 692 693 rKwall_glob(:,1,1,2) = rKwall_glob12(:,1,1) 694 rKwall_glob(:,2,1,2) = rKwall_glob12(:,1,2) 695 rKwall_glob(:,3,1,2) = rKwall_glob12(:,1,3) 696 rKwall_glob(:,4,1,2) = rKwall_glob12(:,2,1) 697 rKwall_glob(:,5,1,2) = rKwall_glob12(:,2,2) 698 rKwall_glob(:,6,1,2) = rKwall_glob12(:,2,3) 699 rKwall_glob(:,7,1,2) = rKwall_glob12(:,3,1) 700 rKwall_glob(:,8,1,2) = rKwall_glob12(:,3,2) 701 rKwall_glob(:,9,1,2) = rKwall_glob12(:,3,3) 702 703 rKwall_glob(:,1,1,3) = rKwall_glob13(:,1,1) 704 rKwall_glob(:,2,1,3) = rKwall_glob13(:,1,2) 705 rKwall_glob(:,3,1,3) = rKwall_glob13(:,1,3) 706 rKwall_glob(:,4,1,3) = rKwall_glob13(:,2,1) 707 rKwall_glob(:,5,1,3) = rKwall_glob13(:,2,2) 708 rKwall_glob(:,6,1,3) = rKwall_glob13(:,2,3) 709 rKwall_glob(:,7,1,3) = rKwall_glob13(:,3,1) 710 rKwall_glob(:,8,1,3) = rKwall_glob13(:,3,2) 711 rKwall_glob(:,9,1,3) = rKwall_glob13(:,3,3) 712 713 rKwall_glob(:,1,2,1) = rKwall_glob21(:,1,1) 714 rKwall_glob(:,2,2,1) = rKwall_glob21(:,1,2) 715 rKwall_glob(:,3,2,1) = rKwall_glob21(:,1,3) 716 rKwall_glob(:,4,2,1) = rKwall_glob21(:,2,1) 717 rKwall_glob(:,5,2,1) = rKwall_glob21(:,2,2) 718 rKwall_glob(:,6,2,1) = rKwall_glob21(:,2,3) 719 rKwall_glob(:,7,2,1) = rKwall_glob21(:,3,1) 720 rKwall_glob(:,8,2,1) = rKwall_glob21(:,3,2) 721 rKwall_glob(:,9,2,1) = rKwall_glob21(:,3,3) 722 723 rKwall_glob(:,1,2,2) = rKwall_glob22(:,1,1) 724 rKwall_glob(:,2,2,2) = rKwall_glob22(:,1,2) 725 rKwall_glob(:,3,2,2) = rKwall_glob22(:,1,3) 726 rKwall_glob(:,4,2,2) = rKwall_glob22(:,2,1) 727 rKwall_glob(:,5,2,2) = rKwall_glob22(:,2,2) 728 rKwall_glob(:,6,2,2) = rKwall_glob22(:,2,3) 729 rKwall_glob(:,7,2,2) = rKwall_glob22(:,3,1) 730 rKwall_glob(:,8,2,2) = rKwall_glob22(:,3,2) 731 rKwall_glob(:,9,2,2) = rKwall_glob22(:,3,3) 732 733 rKwall_glob(:,1,2,3) = rKwall_glob23(:,1,1) 734 rKwall_glob(:,2,2,3) = rKwall_glob23(:,1,2) 735 rKwall_glob(:,3,2,3) = rKwall_glob23(:,1,3) 736 rKwall_glob(:,4,2,3) = rKwall_glob23(:,2,1) 737 rKwall_glob(:,5,2,3) = rKwall_glob23(:,2,2) 738 rKwall_glob(:,6,2,3) = rKwall_glob23(:,2,3) 739 rKwall_glob(:,7,2,3) = rKwall_glob23(:,3,1) 740 rKwall_glob(:,8,2,3) = rKwall_glob23(:,3,2) 741 rKwall_glob(:,9,2,3) = rKwall_glob23(:,3,3) 742 743 rKwall_glob(:,1,3,1) = rKwall_glob31(:,1,1) 744 rKwall_glob(:,2,3,1) = rKwall_glob31(:,1,2) 745 rKwall_glob(:,3,3,1) = rKwall_glob31(:,1,3) 746 rKwall_glob(:,4,3,1) = rKwall_glob31(:,2,1) 747 rKwall_glob(:,5,3,1) = rKwall_glob31(:,2,2) 748 rKwall_glob(:,6,3,1) = rKwall_glob31(:,2,3) 749 rKwall_glob(:,7,3,1) = rKwall_glob31(:,3,1) 750 rKwall_glob(:,8,3,1) = rKwall_glob31(:,3,2) 751 rKwall_glob(:,9,3,1) = rKwall_glob31(:,3,3) 752 753 rKwall_glob(:,1,3,2) = rKwall_glob32(:,1,1) 754 rKwall_glob(:,2,3,2) = rKwall_glob32(:,1,2) 755 rKwall_glob(:,3,3,2) = rKwall_glob32(:,1,3) 756 rKwall_glob(:,4,3,2) = rKwall_glob32(:,2,1) 757 rKwall_glob(:,5,3,2) = rKwall_glob32(:,2,2) 758 rKwall_glob(:,6,3,2) = rKwall_glob32(:,2,3) 759 rKwall_glob(:,7,3,2) = rKwall_glob32(:,3,1) 760 rKwall_glob(:,8,3,2) = rKwall_glob32(:,3,2) 761 rKwall_glob(:,9,3,2) = rKwall_glob32(:,3,3) 762 763 rKwall_glob(:,1,3,3) = rKwall_glob33(:,1,1) 764 rKwall_glob(:,2,3,3) = rKwall_glob33(:,1,2) 765 rKwall_glob(:,3,3,3) = rKwall_glob33(:,1,3) 766 rKwall_glob(:,4,3,3) = rKwall_glob33(:,2,1) 767 rKwall_glob(:,5,3,3) = rKwall_glob33(:,2,2) 768 rKwall_glob(:,6,3,3) = rKwall_glob33(:,2,3) 769 rKwall_glob(:,7,3,3) = rKwall_glob33(:,3,1) 770 rKwall_glob(:,8,3,3) = rKwall_glob33(:,3,2) 771 rKwall_glob(:,9,3,3) = rKwall_glob33(:,3,3) 772 773 rKwall_glob = rKwall_glob*betai*Delt(itseq)*Delt(itseq)*alfi 774 775 else 776c.... nothing happens 777 goto 123 778 endif 779 780123 continue 781 782 endif 783c 784c.... return 785c 786 return 787 end 788 789c--------------------------------------------------------------------- 790c 791c variables for boundary elements 792c 793c--------------------------------------------------------------------- 794 subroutine e3bvarSclr (yl, shdrv, xlb, 795 & shape, WdetJb, bnorm, 796 & flux, dwl ) 797 798 include "common.h" 799c 800 dimension yl(npro,nshl,ndof), shdrv(npro,nsd,nshl), 801 & xlb(npro,nenl,nsd), shape(npro,nshl), 802 & WdetJb(npro), bnorm(npro,nsd), 803 & flux(npro) 804c 805 dimension dxdxib(npro,nsd,nsd), 806 & dxidxb(npro,nsd,nsd), temp(npro), 807 & temp1(npro), temp2(npro), 808 & temp3(npro), 809 & v1(npro,nsd), v2(npro,nsd), 810 & gradSl(npro,nsd), gradS(npro,nsd) 811 812 real*8 diffus(npro), dwl(npro,nshl) 813 814 call getdiffsclr(shape,dwl,yl,diffus) 815c 816c.... ----------------------> Element Metrics <----------------------- 817c 818c.... compute the deformation gradient 819c 820 dxdxib = zero 821c 822 do n = 1, nenl 823 dxdxib(:,1,1) = dxdxib(:,1,1) + xlb(:,n,1) * shdrv(:,1,n) 824 dxdxib(:,1,2) = dxdxib(:,1,2) + xlb(:,n,1) * shdrv(:,2,n) 825 dxdxib(:,1,3) = dxdxib(:,1,3) + xlb(:,n,1) * shdrv(:,3,n) 826 dxdxib(:,2,1) = dxdxib(:,2,1) + xlb(:,n,2) * shdrv(:,1,n) 827 dxdxib(:,2,2) = dxdxib(:,2,2) + xlb(:,n,2) * shdrv(:,2,n) 828 dxdxib(:,2,3) = dxdxib(:,2,3) + xlb(:,n,2) * shdrv(:,3,n) 829 dxdxib(:,3,1) = dxdxib(:,3,1) + xlb(:,n,3) * shdrv(:,1,n) 830 dxdxib(:,3,2) = dxdxib(:,3,2) + xlb(:,n,3) * shdrv(:,2,n) 831 dxdxib(:,3,3) = dxdxib(:,3,3) + xlb(:,n,3) * shdrv(:,3,n) 832 enddo 833c 834c.... compute the normal to the boundary. This is achieved by taking 835c the cross product of two vectors in the plane of the 2-d 836c boundary face. 837c 838 v1 = xlb(:,2,:) - xlb(:,1,:) 839 v2 = xlb(:,3,:) - xlb(:,1,:) 840 841c 842c.....The following are done in order to correct temp1..3 843c based on the results from compressible code. This is done only 844c for wedges, depending on the bounary face.(tri or quad) 845c 846 if (lcsyst .eq. 4) then 847 temp1 = dxdxib(:,2,1) * dxdxib(:,3,3) - 848 & dxdxib(:,2,3) * dxdxib(:,3,1) 849 temp2 = dxdxib(:,3,1) * dxdxib(:,1,3) - 850 & dxdxib(:,3,3) * dxdxib(:,1,1) 851 temp3 = dxdxib(:,1,1) * dxdxib(:,2,3) - 852 & dxdxib(:,1,3) * dxdxib(:,2,1) 853 854 elseif (lcsyst .eq. 1) then 855 temp1 = v1(:,2) * v2(:,3) - v2(:,2) * v1(:,3) 856 temp2 = v2(:,1) * v1(:,3) - v1(:,1) * v2(:,3) 857 temp3 = v1(:,1) * v2(:,2) - v2(:,1) * v1(:,2) 858 else 859 temp1 = - v1(:,2) * v2(:,3) + v2(:,2) * v1(:,3) 860 temp2 = - v2(:,1) * v1(:,3) + v1(:,1) * v2(:,3) 861 temp3 = - v1(:,1) * v2(:,2) + v2(:,1) * v1(:,2) 862 endif 863c 864 temp = one / sqrt ( temp1**2 + temp2**2 + temp3**2 ) 865 bnorm(:,1) = temp1 * temp 866 bnorm(:,2) = temp2 * temp 867 bnorm(:,3) = temp3 * temp 868c 869 870 if (lcsyst .eq. 3) then 871 WdetJb = (1 - Qwtb(lcsyst,intp)) / (four*temp) 872 elseif (lcsyst .eq. 4) then 873 WdetJb = Qwtb(lcsyst,intp) / temp 874 else 875 WdetJb = Qwtb(lcsyst,intp) / (four*temp) 876 endif 877c 878c.... --------------------------> Grad-V <---------------------------- 879c 880c.... compute grad-v for Navier-Stokes terms 881c 882 if (Navier .eq. 1) then 883c 884c.... compute the inverse of deformation gradient 885c 886 dxidxb(:,1,1) = dxdxib(:,2,2) * dxdxib(:,3,3) 887 & - dxdxib(:,3,2) * dxdxib(:,2,3) 888 dxidxb(:,1,2) = dxdxib(:,3,2) * dxdxib(:,1,3) 889 & - dxdxib(:,1,2) * dxdxib(:,3,3) 890 dxidxb(:,1,3) = dxdxib(:,1,2) * dxdxib(:,2,3) 891 & - dxdxib(:,1,3) * dxdxib(:,2,2) 892 temp = one / ( dxidxb(:,1,1) * dxdxib(:,1,1) 893 & + dxidxb(:,1,2) * dxdxib(:,2,1) 894 & + dxidxb(:,1,3) * dxdxib(:,3,1) ) 895 dxidxb(:,1,1) = dxidxb(:,1,1) * temp 896 dxidxb(:,1,2) = dxidxb(:,1,2) * temp 897 dxidxb(:,1,3) = dxidxb(:,1,3) * temp 898 dxidxb(:,2,1) = (dxdxib(:,2,3) * dxdxib(:,3,1) 899 & - dxdxib(:,2,1) * dxdxib(:,3,3)) * temp 900 dxidxb(:,2,2) = (dxdxib(:,1,1) * dxdxib(:,3,3) 901 & - dxdxib(:,3,1) * dxdxib(:,1,3)) * temp 902 dxidxb(:,2,3) = (dxdxib(:,2,1) * dxdxib(:,1,3) 903 & - dxdxib(:,1,1) * dxdxib(:,2,3)) * temp 904 dxidxb(:,3,1) = (dxdxib(:,2,1) * dxdxib(:,3,2) 905 & - dxdxib(:,2,2) * dxdxib(:,3,1)) * temp 906 dxidxb(:,3,2) = (dxdxib(:,3,1) * dxdxib(:,1,2) 907 & - dxdxib(:,1,1) * dxdxib(:,3,2)) * temp 908 dxidxb(:,3,3) = (dxdxib(:,1,1) * dxdxib(:,2,2) 909 & - dxdxib(:,1,2) * dxdxib(:,2,1)) * temp 910c 911c.... compute local-grad-Y 912c 913c 914 gradSl = zero 915 isc=5+isclr 916 do n = 1, nshl 917 gradSl(:,1) = gradSl(:,1) + shdrv(:,1,n) * yl(:,n,isc) 918 gradSl(:,2) = gradSl(:,2) + shdrv(:,2,n) * yl(:,n,isc) 919 gradSl(:,3) = gradSl(:,3) + shdrv(:,3,n) * yl(:,n,isc) 920 enddo 921c 922c.... convert local-grads to global-grads 923c 924 gradS(:,1) = dxidxb(:,1,1) * gradSl(:,1) + 925 & dxidxb(:,2,1) * gradSl(:,2) + 926 & dxidxb(:,3,1) * gradSl(:,3) 927 928c 929 gradS(:,2) = dxidxb(:,1,2) * gradSl(:,1) + 930 & dxidxb(:,2,2) * gradSl(:,2) + 931 & dxidxb(:,3,2) * gradSl(:,3) 932 933 gradS(:,3) = dxidxb(:,1,3) * gradSl(:,1) + 934 & dxidxb(:,2,3) * gradSl(:,2) + 935 & dxidxb(:,3,3) * gradSl(:,3) 936c 937c.... end grad-T 938c 939 endif 940 941 flux = diffus * ( gradS(:,1) * bnorm(:,1) 942 & + gradS(:,2) * bnorm(:,2) 943 & + gradS(:,3) * bnorm(:,3) ) 944c 945c.... return 946c 947 return 948 end 949 950 951c--------------------------------------------------------------------- 952c 953c rotates the local nodal stiffnesses to the goblal frame 954c 955c--------------------------------------------------------------------- 956 subroutine rotatestiff(rKlocal, rotation, 957 & rKglobal) 958 959 include "common.h" 960 961 dimension rKlocal(npro,nsd,nsd), rotation(npro,nsd,nsd), 962 & rKglobal(npro,nsd,nsd) 963 964 dimension tempm(npro,nsd,nsd) 965 966 tempm = zero 967 do i = 1, 3 968 do j = 1, 3 969 do k = 1, 3 970 tempm(:,i,j) = tempm(:,i,j) 971 & + rKlocal(:,i,k) * rotation(:,k,j) 972 enddo 973 enddo 974 enddo 975 976 rKglobal = zero 977 do i = 1, 3 978 do j = 1, 3 979 do k = 1, 3 980 rKglobal(:,i,j) = rKglobal(:,i,j) 981 & + rotation(:,k,i) * tempm(:,k,j) 982 enddo 983 enddo 984 enddo 985 986 return 987 end 988