1 module lhsGkeep 2 3 real*8, allocatable :: lhsG(:) 4 5 end module 6 7c---------------------------------------------------------------------------- 8 9 subroutine keeplhsG 10 11 use lhsGkeep 12 13 include "common.h" 14 15 allocate ( lhsG(nnz*nshg) ) 16 17 return 18 19 end 20 subroutine cmass (shp, shgl, xl, em) 21c 22c---------------------------------------------------------------------- 23c 24c This subroutine computes the consistent mass matrices 25c 26c Ken Jansen, Spring 2000 27c---------------------------------------------------------------------- 28c 29c 30 include "common.h" 31c 32 integer ne, na, nb, nodlcla, nodlclb, iel 33 dimension 34 & shp(nshl,MAXQPT), shgl(nsd,nshl,MAXQPT), 35 & em(npro,nshl,nshl), 36 & xl(npro,nenl,nsd) 37c 38 dimension shape(npro,nshl), shdrv(npro,nsd,nshl), 39 & sgn(npro,nshl), dxidx(npro,nsd,nsd), 40 & shg(npro,nshl,nsd), 41 & WdetJ(npro) 42c 43 em = zero 44c 45c.... loop through the integration points 46c 47 do intp = 1, ngauss ! (these are in common.h) 48c 49c.... get the hierarchic shape functions at this int point 50c 51 call getshp(shp, shgl, sgn, 52 & shape, shdrv, intp) 53c 54c.... calculate the determinant of the jacobian and weight it 55c 56 call e3metric( xl, shdrv,dxidx,shg,WdetJ) 57c 58 do iel = 1, npro 59 do na = 1, nshl 60 do nb = 1, nshl 61 shp2 = shape(iel,na) * shape(iel,nb) 62 em(iel,na,nb) = em(iel,na,nb) + shp2*WdetJ(iel) 63 enddo 64 enddo 65 enddo 66 enddo 67c 68c.... return 69c 70 return 71 end 72 73 subroutine cmassl (shp, shgl, shpf, shglf, xl, em, Qwtf) 74c 75c---------------------------------------------------------------------- 76c 77c This subroutine computes the consistent mass matrices 78c 79c Ken Jansen, Spring 2000 80c---------------------------------------------------------------------- 81c 82c 83 include "common.h" 84c 85 integer ne, na, nb, nodlcla, nodlclb, iel 86 dimension 87 & shp(nshl,MAXQPT), shgl(nsd,nshl,MAXQPT), 88 & shpf(nshl,MAXQPT), shglf(nsd,nshl,MAXQPT), 89 & em(npro,nshl,nshl), eml(npro,nshl), 90 & xl(npro,nenl,nsd) 91c 92 dimension shape(npro,nshl), shdrv(npro,nsd,nshl), 93 & sgn(npro,nshl), dxidx(npro,nsd,nsd), 94 & shg(npro,nshl,nsd), Qwtf(ngaussf), 95 & WdetJ(npro) 96c 97 em = zero 98 eml= zero 99 100 if (ifproj.eq.1)then 101 nods = nshl 102 else 103 nods = nenl 104 endif 105 106c----------------> Get the lumped mass matrix <----------------------- 107 108c 109c.... loop through the integration points 110c 111 do intp = 1, ngaussf ! (these are in common.h) 112c 113c.... get the hierarchic shape functions at this int point 114c 115 call getshp(shpf, shglf, sgn, 116 & shape, shdrv, intp) 117c 118c.... calculate the determinant of the jacobian and weight it 119c 120 call e3metricf( xl, shdrv,dxidx,shg,WdetJ,Qwtf) 121c 122 do i=1,nods !nenl !nshl 123 eml(:,i) = eml(:,i) + shape(:,i)*WdetJ(:) 124 enddo 125 126 enddo ! End loop over quad points. 127 128 129c--------------> Get the consistent mass matrix <------------------------ 130 131 132 shape = zero 133 shdrv = zero 134 dxidx = zero 135 WdetJ = zero 136 shg = zero 137 138c 139c.... loop through the integration points 140c 141 do intp = 1, ngauss ! (these are in common.h) 142 143c.... get the hierarchic shape functions at this int point 144c 145c 146c.... for the mass matrix to be consistent shp and shgl must be 147c.... evaluated with at least higher quadrature than one-pt. quad. 148 149 call getshp(shp, shgl, sgn, 150 & shape, shdrv, intp) 151 152c 153c.... calculate the determinant of the jacobian and weight it 154c 155 call e3metric( xl, shdrv,dxidx,shg,WdetJ) 156c 157 158 do iel = 1, npro 159 do na = 1, nods !nenl !nshl 160 do nb = 1, nods !nenl !nshl 161 shp2 = shape(iel,na) * shape(iel,nb) 162 em(iel,na,nb) = em(iel,na,nb) + shp2*WdetJ(iel) 163 enddo 164 enddo 165 enddo 166 167 enddo ! End loop over quadrature points 168 169 170 171c----------> Obtain a mixed (lumped/consistent) mass matrix <------------ 172 173c... Different combinations of the lump and mass matrices yield 174c... filters of varying widths. In the limiting case were 175c... the entire matrix is lumped, we obtain the same filter as 176c... in getdmc.f. Note that in these equivalent ways of 177c... filtering one-point quadrature is used for shpf and shgl.. 178 179 180 em = (one-flump)*em 181 182 do iel = 1, npro 183 do na = 1, nods !nenl !nshl 184 em(iel,na,na) = em(iel,na,na)+flump*eml(iel,na) 185 enddo 186 enddo 187 188c 189c.... return 190c 191 return 192 end 193 194 195 196 subroutine cmasstl (shp, shgl, shpf, shglf, xl, em, Qwtf) 197c 198c---------------------------------------------------------------------- 199c 200c This subroutine computes the consistent mass matrices 201c 202c Ken Jansen, Spring 2000 203c---------------------------------------------------------------------- 204c 205c 206 include "common.h" 207c 208 integer ne, na, nb, nodlcla, nodlclb, iel 209 dimension 210 & shp(nshl,MAXQPT), shgl(nsd,nshl,MAXQPT), 211 & shpf(nshl,MAXQPT), shglf(nsd,nshl,MAXQPT), 212 & em(npro,nshl,nshl), eml(npro,nshl), 213 & xl(npro,nenl,nsd) 214c 215 dimension shape(npro,nshl), shdrv(npro,nsd,nshl), 216 & sgn(npro,nshl), dxidx(npro,nsd,nsd), 217 & shg(npro,nshl,nsd), Qwtf(ngaussf), 218 & WdetJ(npro) 219c 220 em = zero 221 eml= zero 222 223c----------------> Get the lumped mass matrix <----------------------- 224 225c 226c.... loop through the integration points 227c 228 do intp = 1, ngaussf ! (these are in common.h) 229c 230c.... get the hierarchic shape functions at this int point 231c 232 call getshp(shpf, shglf, sgn, 233 & shape, shdrv, intp) 234c 235c.... calculate the determinant of the jacobian and weight it 236c 237 call e3metricf( xl, shdrv,dxidx,shg,WdetJ,Qwtf) 238c 239 do i=1,nshl 240 eml(:,i) = eml(:,i) + shape(:,i)*WdetJ(:) 241 enddo 242 243 enddo ! End loop over quad points. 244 245 246c--------------> Get the consistent mass matrix <------------------------ 247 248 249 shape = zero 250 shdrv = zero 251 dxidx = zero 252 WdetJ = zero 253 shg = zero 254 255c 256c.... loop through the integration points 257c 258 do intp = 1, ngauss ! (these are in common.h) 259 260c.... get the hierarchic shape functions at this int point 261c 262c 263c.... for the mass matrix to be consistent shp and shgl must be 264c.... evaluated with at least higher quadrature than one-pt. quad. 265 266 call getshp(shp, shgl, sgn, 267 & shape, shdrv, intp) 268 269c 270c.... calculate the determinant of the jacobian and weight it 271c 272 call e3metric( xl, shdrv,dxidx,shg,WdetJ) 273c 274 275 do iel = 1, npro 276 do na = 1, nshl 277 do nb = 1, nshl 278 shp2 = shape(iel,na) * shape(iel,nb) 279 em(iel,na,nb) = em(iel,na,nb) + shp2*WdetJ(iel) 280 enddo 281 enddo 282 enddo 283 284 enddo ! End loop over quadrature points 285 286 287 288c----------> Obtain a mixed (lumped/consistent) mass matrix <------------ 289 290c... Different combinations of the lump and mass matrices yield 291c... filters of varying widths. In the limiting case were 292c... the entire matrix is lumped, we obtain the same filter as 293c... in getdmc.f. Note that in these equivalent ways of 294c... filtering one-point quadrature is used for shpf and shgl.. 295 296 297 do iel = 1, npro 298 do na = 1, nshl 299 em(iel,na,na) = eml(iel,na) 300 enddo 301 enddo 302 303c 304c.... return 305c 306 return 307 end 308 309c----------------------------------------------------------------------- 310c 311c compute the metrics of the mapping from global to local 312c coordinates and the jacobian of the mapping (weighted by 313c the quadrature weight 314c 315c----------------------------------------------------------------------- 316 subroutine e3metricf( xl, shgl, dxidx, 317 & shg, WdetJ, Qwtf) 318 319 include "common.h" 320 321 real*8 xl(npro,nenl,nsd), shgl(npro,nsd,nshl), 322 & dxidx(npro,nsd,nsd), shg(npro,nshl,nsd), 323 & WdetJ(npro), Qwtf(ngaussf) 324 325 real*8 dxdxi(npro,nsd,nsd), tmp(npro) 326 327c 328c.... compute the deformation gradient 329c 330 dxdxi = zero 331c 332 do n = 1, nenl 333 dxdxi(:,1,1) = dxdxi(:,1,1) + xl(:,n,1) * shgl(:,1,n) 334 dxdxi(:,1,2) = dxdxi(:,1,2) + xl(:,n,1) * shgl(:,2,n) 335 dxdxi(:,1,3) = dxdxi(:,1,3) + xl(:,n,1) * shgl(:,3,n) 336 dxdxi(:,2,1) = dxdxi(:,2,1) + xl(:,n,2) * shgl(:,1,n) 337 dxdxi(:,2,2) = dxdxi(:,2,2) + xl(:,n,2) * shgl(:,2,n) 338 dxdxi(:,2,3) = dxdxi(:,2,3) + xl(:,n,2) * shgl(:,3,n) 339 dxdxi(:,3,1) = dxdxi(:,3,1) + xl(:,n,3) * shgl(:,1,n) 340 dxdxi(:,3,2) = dxdxi(:,3,2) + xl(:,n,3) * shgl(:,2,n) 341 dxdxi(:,3,3) = dxdxi(:,3,3) + xl(:,n,3) * shgl(:,3,n) 342 enddo 343c 344c.... compute the inverse of deformation gradient 345c 346 dxidx(:,1,1) = dxdxi(:,2,2) * dxdxi(:,3,3) 347 & - dxdxi(:,3,2) * dxdxi(:,2,3) 348 dxidx(:,1,2) = dxdxi(:,3,2) * dxdxi(:,1,3) 349 & - dxdxi(:,1,2) * dxdxi(:,3,3) 350 dxidx(:,1,3) = dxdxi(:,1,2) * dxdxi(:,2,3) 351 & - dxdxi(:,1,3) * dxdxi(:,2,2) 352 tmp = one / ( dxidx(:,1,1) * dxdxi(:,1,1) 353 & + dxidx(:,1,2) * dxdxi(:,2,1) 354 & + dxidx(:,1,3) * dxdxi(:,3,1) ) 355 dxidx(:,1,1) = dxidx(:,1,1) * tmp 356 dxidx(:,1,2) = dxidx(:,1,2) * tmp 357 dxidx(:,1,3) = dxidx(:,1,3) * tmp 358 dxidx(:,2,1) = (dxdxi(:,2,3) * dxdxi(:,3,1) 359 & - dxdxi(:,2,1) * dxdxi(:,3,3)) * tmp 360 dxidx(:,2,2) = (dxdxi(:,1,1) * dxdxi(:,3,3) 361 & - dxdxi(:,3,1) * dxdxi(:,1,3)) * tmp 362 dxidx(:,2,3) = (dxdxi(:,2,1) * dxdxi(:,1,3) 363 & - dxdxi(:,1,1) * dxdxi(:,2,3)) * tmp 364 dxidx(:,3,1) = (dxdxi(:,2,1) * dxdxi(:,3,2) 365 & - dxdxi(:,2,2) * dxdxi(:,3,1)) * tmp 366 dxidx(:,3,2) = (dxdxi(:,3,1) * dxdxi(:,1,2) 367 & - dxdxi(:,1,1) * dxdxi(:,3,2)) * tmp 368 dxidx(:,3,3) = (dxdxi(:,1,1) * dxdxi(:,2,2) 369 & - dxdxi(:,1,2) * dxdxi(:,2,1)) * tmp 370c 371c WdetJ = Qwt(lcsyst,intp) / tmp 372 373 WdetJ = Qwtf(intp) / tmp 374c 375c.... compute the global gradient of shape-functions 376c 377 do n = 1, nshl 378 shg(:,n,1) = shgl(:,1,n) * dxidx(:,1,1) + 379 & shgl(:,2,n) * dxidx(:,2,1) + 380 & shgl(:,3,n) * dxidx(:,3,1) 381 shg(:,n,2) = shgl(:,1,n) * dxidx(:,1,2) + 382 & shgl(:,2,n) * dxidx(:,2,2) + 383 & shgl(:,3,n) * dxidx(:,3,2) 384 shg(:,n,3) = shgl(:,1,n) * dxidx(:,1,3) + 385 & shgl(:,2,n) * dxidx(:,2,3) + 386 & shgl(:,3,n) * dxidx(:,3,3) 387 enddo 388 389 return 390 end 391 392 393 394