1 subroutine e3wmlt (shp, shg, WdetJ, 2 & ri, rmi, rl, 3 & rml, stiff, EGmass ) 4c 5c---------------------------------------------------------------------- 6c 7c Up to now most of the terms have not been multiplied by the 8c shape function from the weight space. Now that we have collected 9c all the terms that the shape function (and its derivatives for the 10c terms that were integrated by parts) multiplies, in this routine 11c we carry out that multiplication. After these operations we will 12c have this quadrature points contribution to the integral for the 13c residual (i.e. g^e_b(xi^l) D(xi^l) QW(xi^l) where e is the element, 14c b is the local node number on the element, l is the quadrature point 15c D is the determinate of the Jacobian of the mapping, and QW is this 16c quadrature points weight....WHEW). When we add up all of the 17c integration points we get G^e_b which we will assemble in the 18c subroutine LOCAL to form G_B. 19c 20c This routine also forms the LHS contribution from the LS term 21c and the diffusion term which has been partially constructed and 22c placed in stiff. Those familiar with elasticity might recognize 23c this naming convention since this is like a stiffness matrix that 24c if you had a linear problem would be calculated once and saved for 25c all time. 26c 27c ri: LS, Diffusion, Convective and mass, 28c rmi: LS, Diffusion and mass, 29c stiff: LS, Diffusion. 30c 31c input: 32c shp (nshl) : element shape-functions 33c shg (npro,nshl,nsd) : element grad-shape-functions 34c WdetJ (npro) : weighted Jacobian 35c ri (npro,nflow*(nsd+1)) : partial residual 36c rmi (npro,nflow*(nsd+1)) : partial modified residual 37c stiff (npro,nsd*nflow,nsd*nflow) : stiffness matrix 38c ( K_ij + A_i tau A_j ) 39c 40c output: 41c rl (npro,nshl,nflow) : residual 42c rml (npro,nshl,nflow) : modified residual 43c EGmass (npro,nedof,nedof) : element LHS tangent mass matrix 44c 45c 46c Zdenek Johan, Summer 1990. (Modified from e2assm.f) 47c Zdenek Johan, Winter 1991. (Fortran 90) 48c Kenneth Jansen, Winter 1997 Prim. variables 49c---------------------------------------------------------------------- 50c 51 include "common.h" 52c 53 dimension shp(npro,nshl), 54 & shg(npro,nshl,nsd), 55 & WdetJ(npro), ri(npro,nflow*(nsd+1)), 56 & rmi(npro,nflow*(nsd+1)), stiff(npro,3*nflow,3*nflow), 57 & rl(npro,nshl,nflow), rml(npro,nshl,nflow), 58 & EGmass(npro,nedof,nedof) 59c 60 dimension shg1(npro), shg2(npro), 61 & shg3(npro), stif1(npro,nflow,nflow), 62 & stif2(npro,nflow,nflow), stif3(npro,nflow,nflow) 63 64 ttim(29) = ttim(29) - secs(0.0) 65c 66c.... ----------------------------> RHS <---------------------------- 67c 68c.... add spatial contribution to rl and rml 69c 70c.... ires = 1 or 3 71c 72 if ((ires .eq. 1) .or. (ires .eq. 3)) then 73c 74 do i = 1, nshl 75 rl(:,i,1) = rl(:,i,1) + WdetJ * ( 76 & shg(:,i,1) * ri(:, 1) + shg(:,i,2) * ri(:, 6) 77 & + shg(:,i,3) * ri(:,11) ) 78 rl(:,i,2) = rl(:,i,2) + WdetJ * ( 79 & shg(:,i,1) * ri(:, 2) + shg(:,i,2) * ri(:, 7) 80 & + shg(:,i,3) * ri(:,12) ) 81 rl(:,i,3) = rl(:,i,3) + WdetJ * ( 82 & shg(:,i,1) * ri(:, 3) + shg(:,i,2) * ri(:, 8) 83 & + shg(:,i,3) * ri(:,13) ) 84 rl(:,i,4) = rl(:,i,4) + WdetJ * ( 85 & shg(:,i,1) * ri(:, 4) + shg(:,i,2) * ri(:, 9) 86 & + shg(:,i,3) * ri(:,14) ) 87 rl(:,i,5) = rl(:,i,5) + WdetJ * ( 88 & shg(:,i,1) * ri(:, 5) + shg(:,i,2) * ri(:,10) 89 & + shg(:,i,3) * ri(:,15) ) 90 enddo 91c 92! flops = flops + 36*nshl*npro 93 endif 94c 95c.... ires = 2 or 3 96c 97 if ((ires .eq. 2) .or. (ires .eq. 3)) then 98 do i = 1, nshl 99 rml(:,i,1) = rml(:,i,1) + WdetJ * ( 100 & shg(:,i,1) * rmi(:, 1) + shg(:,i,2) * rmi(:, 6) 101 & + shg(:,i,3) * rmi(:,11) 102 & + shp(:,i) * rmi(:,16) ) 103 rml(:,i,2) = rml(:,i,2) + WdetJ * ( 104 & shg(:,i,1) * rmi(:, 2) + shg(:,i,2) * rmi(:, 7) 105 & + shg(:,i,3) * rmi(:,12) 106 & + shp(:,i) * rmi(:,17) ) 107 rml(:,i,3) = rml(:,i,3) + WdetJ * ( 108 & shg(:,i,1) * rmi(:, 3) + shg(:,i,2) * rmi(:, 8) 109 & + shg(:,i,3) * rmi(:,13) 110 & + shp(:,i) * rmi(:,18) ) 111 rml(:,i,4) = rml(:,i,4) + WdetJ * ( 112 & shg(:,i,1) * rmi(:, 4) + shg(:,i,2) * rmi(:, 9) 113 & + shg(:,i,3) * rmi(:,14) 114 & + shp(:,i) * rmi(:,19) ) 115 rml(:,i,5) = rml(:,i,5) + WdetJ * ( 116 & shg(:,i,1) * rmi(:, 5) + shg(:,i,2) * rmi(:,10) 117 & + shg(:,i,3) * rmi(:,15) 118 & + shp(:,i) * rmi(:,20) ) 119 enddo 120c 121! flops = flops + (15+45*nshl)*npro 122 endif 123c 124c.... add temporal contribution to rl 125c 126 if (ngauss .eq. 1 .and. nshl.eq.4) then !already ex. integ mass 127 if(matflg(5,1).ge.1) then 128 do i = 1, nshl 129 rl(:,i,2) = rl(:,i,2) + shp(:,i) * WdetJ * ri(:,17) 130 rl(:,i,3) = rl(:,i,3) + shp(:,i) * WdetJ * ri(:,18) 131 rl(:,i,4) = rl(:,i,4) + shp(:,i) * WdetJ * ri(:,19) 132 rl(:,i,5) = rl(:,i,5) + shp(:,i) * WdetJ * ri(:,20) 133 enddo 134 endif 135 else !check for a body force 136 do i = 1, nshl 137 rl(:,i,1) = rl(:,i,1) + shp(:,i) * WdetJ * ri(:,16) 138 rl(:,i,2) = rl(:,i,2) + shp(:,i) * WdetJ * ri(:,17) 139 rl(:,i,3) = rl(:,i,3) + shp(:,i) * WdetJ * ri(:,18) 140 rl(:,i,4) = rl(:,i,4) + shp(:,i) * WdetJ * ri(:,19) 141 rl(:,i,5) = rl(:,i,5) + shp(:,i) * WdetJ * ri(:,20) 142 enddo 143 144! flops = flops + 11*nshl*npro 145 endif 146 147c 148c.... ----------------------------> LHS <---------------------------- 149c 150 if (lhs .eq. 1) then 151c 152c.... loop through columns (nodes j) 153c 154 do j = 1, nshl 155 j0 = nflow * (j - 1) 156c 157c.... set up factors 158c 159 shg1 = WdetJ * shg(:,j,1) 160 shg2 = WdetJ * shg(:,j,2) 161 shg3 = WdetJ * shg(:,j,3) 162c 163c.... loop through d.o.f.'s 164c 165 do jdof = 1, nflow 166 do idof = 1, nflow 167 idof2 = idof + nflow 168 jdof2 = jdof + nflow 169 170 idof3 = idof2 + nflow 171 jdof3 = jdof2 + nflow 172c 173c.... calculate weighted stiffness matrix (first part) 174c 175 stif1(:,idof,jdof) = shg1 * stiff(:,idof,jdof) 176 & + shg2 * stiff(:,idof,jdof2) 177 & + shg3 * stiff(:,idof,jdof3) 178 stif2(:,idof,jdof) = shg1 * stiff(:,idof2,jdof) 179 & + shg2 * stiff(:,idof2,jdof2) 180 & + shg3 * stiff(:,idof2,jdof3) 181 stif3(:,idof,jdof) = shg1 * stiff(:,idof3,jdof) 182 & + shg2 * stiff(:,idof3,jdof2) 183 & + shg3 * stiff(:,idof3,jdof3) 184 enddo 185 enddo 186c 187c.... loop through rows (nodes i) 188c 189 do i = 1, nshl 190 i0 = nflow * (i - 1) 191c 192c.... add contribution of stiffness to EGmass 193c 194 do jdof = 1, nflow 195 EGmass(:,i0+1,j0+jdof) = EGmass(:,i0+1,j0+jdof) 196 & + shg(:,i,1) * stif1(:,1,jdof) 197 & + shg(:,i,2) * stif2(:,1,jdof) 198 & + shg(:,i,3) * stif3(:,1,jdof) 199 EGmass(:,i0+2,j0+jdof) = EGmass(:,i0+2,j0+jdof) 200 & + shg(:,i,1) * stif1(:,2,jdof) 201 & + shg(:,i,2) * stif2(:,2,jdof) 202 & + shg(:,i,3) * stif3(:,2,jdof) 203 EGmass(:,i0+3,j0+jdof) = EGmass(:,i0+3,j0+jdof) 204 & + shg(:,i,1) * stif1(:,3,jdof) 205 & + shg(:,i,2) * stif2(:,3,jdof) 206 & + shg(:,i,3) * stif3(:,3,jdof) 207 EGmass(:,i0+4,j0+jdof) = EGmass(:,i0+4,j0+jdof) 208 & + shg(:,i,1) * stif1(:,4,jdof) 209 & + shg(:,i,2) * stif2(:,4,jdof) 210 & + shg(:,i,3) * stif3(:,4,jdof) 211 EGmass(:,i0+5,j0+jdof) = EGmass(:,i0+5,j0+jdof) 212 & + shg(:,i,1) * stif1(:,5,jdof) 213 & + shg(:,i,2) * stif2(:,5,jdof) 214 & + shg(:,i,3) * stif3(:,5,jdof) 215 enddo 216c 217c.... end loop on rows 218c 219 enddo 220c 221c.... end loop on columns 222c 223 enddo 224c 225c.... end left hand side assembly 226c 227 endif 228 229 ttim(29) = ttim(29) + secs(0.0) 230c 231c.... return 232c 233 return 234 end 235c 236c 237c 238 subroutine e3wmltSclr(shp, shg, WdetJ, 239 & rti, rtl, 240 & stifft, EGmasst ) 241c 242c---------------------------------------------------------------------- 243c 244c Up to now most of the terms have not been multiplied by the 245c shape function from the weight space. Now that we have collected 246c all the terms that the shape function (and its derivatives for the 247c terms that were integrated by parts) multiplies, in this routine 248c we carry out that multiplication. After these operations we will 249c have this quadrature points contribution to the integral for the 250c residual (i.e. g^e_b(xi^l) D(xi^l) QW(xi^l) where e is the element, 251c b is the local node number on the element, l is the quadrature point 252c D is the determinate of the Jacobian of the mapping, and QW is this 253c quadrature points weight....WHEW). When we add up all of the 254c integration points we get G^e_b which we will assemble in the 255c subroutine LOCAL to form G_B. 256c 257c This routine also forms the LHS contribution from the LS term 258c and the diffusion term which has been partially constructed and 259c placed in stiff. Those familiar with elasticity might recognize 260c this naming convention since this is like a stiffness matrix that 261c if you had a linear problem would be calculated once and saved for 262c all time. 263c 264c ri: LS, Diffusion, Convective and mass, 265c stiff: LS, Diffusion. 266c 267c input: 268c shp (npro,nshl) : element shape-functions 269c shg (npro,nshl,nsd) : element grad-shape-functions 270c WdetJ (npro) : weighted Jacobian 271c rti (npro,nsd+1) : partial residual 272c stifft (npro,nsd,nsd) : stiffness matrix 273c ( K_ij + A_i tau A_j ) 274c 275c output: 276c rtl (npro,nshl) : residual (will end up being G^e_b) 277c EGmasst (npro,nshape,nshape) : element LHS tangent mass matrix 278c (will end up being dG^e_b/dY_a) 279c 280c 281c Zdenek Johan, Summer 1990. (Modified from e2assm.f) 282c Zdenek Johan, Winter 1991. (Fortran 90) 283c Kenneth Jansen, Winter 1997 Prim. variables 284c---------------------------------------------------------------------- 285c 286 include "common.h" 287c 288 dimension shp(npro,nshl), shg(npro,nshl,nsd), 289 & WdetJ(npro), rti(npro,nsd+1), 290 & rmti(npro,nsd+1), stifft(npro,nsd,nsd), 291 & rtl(npro,nshl), rmtl(npro,nshl), 292 & EGmasst(npro,nshape,nshape) 293c 294 dimension shg1(npro), shg2(npro), 295 & shg3(npro), stift1(npro), 296 & stift2(npro), stift3(npro) 297 298 ttim(29) = ttim(29) - tmr(0.0) 299c 300c.... ----------------------------> RHS <---------------------------- 301c 302c.... add spatial contribution to rl and rml 303c 304c.... ires = 1 or 3 305c 306 if ((ires .eq. 1) .or. (ires .eq. 3)) then 307c 308 do i = 1, nshl 309 rtl(:,i) = rtl(:,i) + WdetJ * ( 310 & shg(:,i,1) * rti(:,1) + shg(:,i,2) * rti(:,2) 311 & + shg(:,i,3) * rti(:,3) ) 312 313 enddo 314! flops = flops + 36*nshl*npro 315 endif 316c 317c.... ires = 2 or 3 318c 319c if ((ires .eq. 2) .or. (ires .eq. 3)) then 320c do i = 1, nshl 321c rml(:,i,1) = rml(:,i,1) + WdetJ * ( 322c & shg(:,i,1) * rmi(:, 1) + shg(:,i,2) * rmi(:, 6) 323c & + shg(:,i,3) * rmi(:,11) + shp(:,i) * rmi(:,16) ) 324c rml(:,i,2) = rml(:,i,2) + WdetJ * ( 325c & shg(:,i,1) * rmi(:, 2) + shg(:,i,2) * rmi(:, 7) 326c & + shg(:,i,3) * rmi(:,12) + shp(:,i) * rmi(:,17) ) 327c rml(:,i,3) = rml(:,i,3) + WdetJ * ( 328c & shg(:,i,1) * rmi(:, 3) + shg(:,i,2) * rmi(:, 8) 329c & + shg(:,i,3) * rmi(:,13) + shp(:,i) * rmi(:,18) ) 330c rml(:,i,4) = rml(:,i,4) + WdetJ * ( 331c & shg(:,i,1) * rmi(:, 4) + shg(:,i,2) * rmi(:, 9) 332c & + shg(:,i,3) * rmi(:,14) + shp(:,i) * rmi(:,19) ) 333c rml(:,i,5) = rml(:,i,5) + WdetJ * ( 334c & shg(:,i,1) * rmi(:, 5) + shg(:,i,2) * rmi(:,10) 335c & + shg(:,i,3) * rmi(:,15) + shp(:,i) * rmi(:,20) ) 336c enddo 337c 338c ! flops = flops + (15+45*nshl)*npro 339c endif 340 341c 342c.... add temporal contribution to rl 343c 344 345 do i = 1, nshl 346 rtl(:,i) = rtl(:,i) + shp(:,i) * WdetJ * rti(:,4) 347 enddo 348 349! flops = flops + 11*nshl*npro 350 351c 352c.... ----------------------------> LHS <---------------------------- 353c 354 if (lhs .eq. 1) then 355c 356c.... loop through columns (nodes j) 357c 358 do j = 1, nshl 359c 360c.... set up factors 361c 362 shg1 = WdetJ * shg(:,j,1) 363 shg2 = WdetJ * shg(:,j,2) 364 shg3 = WdetJ * shg(:,j,3) 365c 366c.... calculate weighted stiffness matrix (first part) 367c 368 stift1(:) = shg1 * stifft(:,1,1) 369 & + shg2 * stifft(:,1,2) 370 & + shg3 * stifft(:,1,3) 371 stift2(:) = shg1 * stifft(:,2,1) 372 & + shg2 * stifft(:,2,2) 373 & + shg3 * stifft(:,2,3) 374 stift3(:) = shg1 * stifft(:,3,1) 375 & + shg2 * stifft(:,3,2) 376 & + shg3 * stifft(:,3,3) 377c 378c.... loop through rows (nodes i) 379c 380 do i = 1, nshl 381c 382c.... add contribution of stiffness to EGmass 383c 384 EGmasst(:,i,j) = EGmasst(:,i,j) 385 & + shg(:,i,1) * stift1(:) 386 & + shg(:,i,2) * stift2(:) 387 & + shg(:,i,3) * stift3(:) 388 389c 390c.... end loop on rows 391c 392 enddo 393c 394c.... end loop on columns 395c 396 enddo 397c 398c.... end left hand side assembly 399c 400 endif 401 402 ttim(29) = ttim(29) + tmr(0.0) 403c 404c.... return 405c 406 return 407 end 408