1 subroutine hessian ( y, x, 2 & shp, shgl, iBC, 3 & shpb, shglb, iper, 4 & ilwork, uhess, gradu ) 5 use pointer_data ! brings in the pointers for the blocked arrays 6 7 include "common.h" 8c 9 dimension y(nshg,ndof), 10 & x(numnp,nsd), iBC(nshg), 11 & iper(nshg) 12c 13 dimension shp(MAXTOP,maxsh,MAXQPT), 14 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 15 & shpb(MAXTOP,maxsh,MAXQPT), 16 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 17c 18 dimension gradu(nshg,9), rmass(nshg), 19 & uhess(nshg,27) 20c 21 dimension ilwork(nlwork) 22 23c 24 gradu = zero 25 rmass = zero 26 27 do iblk = 1, nelblk 28 iel = lcblk(1,iblk) 29 lelCat = lcblk(2,iblk) 30 lcsyst = lcblk(3,iblk) 31 iorder = lcblk(4,iblk) 32 nenl = lcblk(5,iblk) ! no. of vertices per element 33 nshl = lcblk(10,iblk) 34 mattyp = lcblk(7,iblk) 35 ndofl = lcblk(8,iblk) 36 nsymdl = lcblk(9,iblk) 37 npro = lcblk(1,iblk+1) - iel 38 ngauss = nint(lcsyst) 39c 40 call velocity_gradient ( y, 41 & x, 42 & shp(lcsyst,1:nshl,:), 43 & shgl(lcsyst,:,1:nshl,:), 44 & mien(iblk)%p, 45 & gradu, 46 & rmass ) 47 48 end do 49 50c 51 call reconstruct( rmass, gradu, iBC, iper, ilwork, 9 ) 52 53 uhess = zero 54 rmass = zero 55 56 do iblk = 1, nelblk 57 iel = lcblk(1,iblk) 58 lelCat = lcblk(2,iblk) 59 lcsyst = lcblk(3,iblk) 60 iorder = lcblk(4,iblk) 61 nenl = lcblk(5,iblk) ! no. of vertices per element 62 nshl = lcblk(10,iblk) 63 mattyp = lcblk(7,iblk) 64 ndofl = lcblk(8,iblk) 65 nsymdl = lcblk(9,iblk) 66 npro = lcblk(1,iblk+1) - iel 67 ngauss = nint(lcsyst) 68c 69 70 call velocity_hessian ( gradu, 71 & x, 72 & shp(lcsyst,1:nshl,:), 73 & shgl(lcsyst,:,1:nshl,:), 74 & mien(iblk)%p, 75 & uhess, 76 & rmass ) 77 end do 78 79 80 call reconstruct( rmass, uhess, iBC, iper, ilwork, 27 ) 81c 82 return 83 end 84 85c----------------------------------------------------------------------------- 86 87 subroutine velocity_gradient ( y, x, shp, shgl, 88 & ien, gradu, rmass ) 89 90 include "common.h" 91c 92 dimension y(nshg,ndof), x(numnp,nsd), 93 & shp(nshl,ngauss), shgl(nsd,nshl,ngauss), 94 & ien(npro,nshl), gradu(nshg,9), 95 & shdrv(npro,nsd,nshl), shape( npro, nshl ), 96 & gradul(npro,9) , rmass( nshg ) 97c 98 dimension yl(npro,nshl,ndof), xl(npro,nenl,nsd), 99 & ql(npro,nshl,9), dxidx(npro,nsd,nsd), 100 & WdetJ(npro), rmassl(npro,nshl) 101c 102c 103 dimension sgn(npro,nshl) 104c 105c.... create the matrix of mode signs for the hierarchic basis 106c functions. 107c 108 do i=1,nshl 109 where ( ien(:,i) < 0 ) 110 sgn(:,i) = -one 111 elsewhere 112 sgn(:,i) = one 113 endwhere 114 enddo 115 116c 117c.... gather the variables 118c 119 120 call localy (y, yl, ien, ndof, 'gather ') 121 call localx (x, xl, ien, nsd, 'gather ') 122c 123c.... get the element residuals 124c 125 ql = zero 126 rmassl = zero 127 128 do intp = 1, ngauss 129 130 if ( Qwt( lcsyst, intp ) .eq. zero ) cycle 131 132 gradul = zero 133 call getshp( shp, shgl, sgn, shape, shdrv ) 134 call local_gradient( yl(:,:,2:4), 3, shdrv, xl, 135 & gradul , dxidx, WdetJ ) 136 137c.... assemble contribution of gradu to each element node 138c 139 do i=1,nshl 140 do j = 1, 9 141 ql(:,i,j) = ql(:,i,j)+shape(:,i)*WdetJ*gradul(:,j) 142 end do 143 144 rmassl(:,i) = rmassl(:,i) + shape(:,i)*WdetJ 145 146 end do 147 148 end do 149c 150c 151 call local (gradu, ql, ien, 9, 'scatter ') 152 call local (rmass, rmassl, ien, 1, 'scatter ') 153c 154c.... end 155c 156 return 157 end 158 159 160c----------------------------------------------------------------------------- 161 162 subroutine velocity_hessian ( gradu, x, shp, shgl, 163 & ien, uhess, rmass ) 164 165 include "common.h" 166c 167 dimension gradu(nshg,9), x(numnp,nsd), 168 & shp(nshl,ngauss), shgl(nsd,nshl,ngauss), 169 & ien(npro,nshl), uhess(nshg,27), 170 & shdrv(npro,nsd,nshl), shape( npro, nshl ), 171 & uhessl(npro,27), rmass( nshg ) 172c 173 dimension gradul(npro,nshl,9), xl(npro,nenl,nsd), 174 & ql(npro,nshl,27), dxidx(npro,nsd,nsd), 175 & WdetJ(npro), rmassl(npro, nshl) 176c 177c 178 dimension sgn(npro,nshl) 179c 180c.... create the matrix of mode signs for the hierarchic basis 181c functions. 182c 183 do i=1,nshl 184 where ( ien(:,i) < 0 ) 185 sgn(:,i) = -one 186 elsewhere 187 sgn(:,i) = one 188 endwhere 189 enddo 190 191c 192c.... gather the variables 193c 194 195 call local (gradu, gradul, ien, 9 , 'gather ') 196 call localx (x, xl, ien, nsd, 'gather ') 197c 198c.... get the element residuals 199c 200 ql = zero 201 rmassl = zero 202 203 do intp = 1, ngauss 204 205 if ( Qwt( lcsyst, intp ) .eq. zero ) cycle 206 207 uhessl = zero 208 call getshp( shp, shgl, sgn, shape, shdrv ) 209 call local_gradient( gradul, 9, shdrv, xl, 210 & uhessl , dxidx, WdetJ ) 211 212c.... assemble contribution of gradu ., 213c 214 do i=1,nshl 215 do j = 1,27 216 ql(:,i,j)=ql(:,i,j)+shape(:,i)*WdetJ*uhessl(:,j ) 217 end do 218 219 rmassl(:,i) = rmassl(:,i) + shape(:,i)*WdetJ 220 end do 221 222 end do 223c 224c 225 call local (uhess, ql, ien, 27, 'scatter ') 226 call local (rmass, rmassl, ien, 1, 'scatter ') 227c 228c.... end 229c 230 return 231 end 232 233 234c-------------------------------------------------------------------- 235 subroutine reconstruct( rmass, qres, iBC, iper, ilwork, vsize ) 236 237 include "common.h" 238 239 integer vsize 240 dimension rmass(nshg), qres( nshg, vsize), 241 & iBC(nshg), iper(nshg), ilwork(nlwork) 242c 243c 244c.... compute qi for node A, i.e., qres <-- qres/rmass 245c 246 if (numpe > 1) then 247 call commu ( qres , ilwork, vsize , 'in ') 248 call commu ( rmass , ilwork, 1 , 'in ') 249 endif 250c 251c take care of periodic boundary conditions 252c 253 do j= 1,nshg 254 if ((btest(iBC(j),10))) then 255 i = iper(j) 256 rmass(i) = rmass(i) + rmass(j) 257 qres(i,:) = qres(i,:) + qres(j,:) 258 endif 259 enddo 260 261 do j= 1,nshg 262 if ((btest(iBC(j),10))) then 263 i = iper(j) 264 rmass(j) = rmass(i) 265 qres(j,:) = qres(i,:) 266 endif 267 enddo 268c 269c.... invert the diagonal mass matrix and find q 270c 271 rmass = one/rmass 272 273 do i=1,vsize 274 qres(:,i) = rmass*qres(:,i) 275 enddo 276 277 if(numpe > 1) then 278 call commu (qres, ilwork, vsize, 'out') 279 endif 280 281c.... return 282c 283 return 284 end 285 286c------------------------------------------------------------------------- 287 288 subroutine local_gradient ( vector, vsize, shgl, xl, 289 & gradient, dxidx, WdetJ ) 290c 291c 292 include "common.h" 293c 294c passed arrays 295 296 integer vsize 297c 298 dimension vector(npro,nshl,vsize), 299 & shgl(npro,nsd,nshl), xl(npro,nenl,nsd), 300 & gradient(npro,vsize*3), shg(npro,nshl,nsd), 301 & dxidx(npro,nsd,nsd), WdetJ(npro) 302c 303c local arrays 304c 305 dimension tmp(npro), dxdxi(npro,nsd,nsd) 306 307c 308c.... compute the deformation gradient 309c 310 dxdxi = zero 311c 312 do n = 1, nenl 313 dxdxi(:,1,1) = dxdxi(:,1,1) + xl(:,n,1) * shgl(:,1,n) 314 dxdxi(:,1,2) = dxdxi(:,1,2) + xl(:,n,1) * shgl(:,2,n) 315 dxdxi(:,1,3) = dxdxi(:,1,3) + xl(:,n,1) * shgl(:,3,n) 316 dxdxi(:,2,1) = dxdxi(:,2,1) + xl(:,n,2) * shgl(:,1,n) 317 dxdxi(:,2,2) = dxdxi(:,2,2) + xl(:,n,2) * shgl(:,2,n) 318 dxdxi(:,2,3) = dxdxi(:,2,3) + xl(:,n,2) * shgl(:,3,n) 319 dxdxi(:,3,1) = dxdxi(:,3,1) + xl(:,n,3) * shgl(:,1,n) 320 dxdxi(:,3,2) = dxdxi(:,3,2) + xl(:,n,3) * shgl(:,2,n) 321 dxdxi(:,3,3) = dxdxi(:,3,3) + xl(:,n,3) * shgl(:,3,n) 322 enddo 323c 324c.... compute the inverse of deformation gradient 325c 326 dxidx(:,1,1) = dxdxi(:,2,2) * dxdxi(:,3,3) 327 & - dxdxi(:,3,2) * dxdxi(:,2,3) 328 dxidx(:,1,2) = dxdxi(:,3,2) * dxdxi(:,1,3) 329 & - dxdxi(:,1,2) * dxdxi(:,3,3) 330 dxidx(:,1,3) = dxdxi(:,1,2) * dxdxi(:,2,3) 331 & - dxdxi(:,1,3) * dxdxi(:,2,2) 332 tmp = one / ( dxidx(:,1,1) * dxdxi(:,1,1) 333 & + dxidx(:,1,2) * dxdxi(:,2,1) 334 & + dxidx(:,1,3) * dxdxi(:,3,1) ) 335 dxidx(:,1,1) = dxidx(:,1,1) * tmp 336 dxidx(:,1,2) = dxidx(:,1,2) * tmp 337 dxidx(:,1,3) = dxidx(:,1,3) * tmp 338 dxidx(:,2,1) = (dxdxi(:,2,3) * dxdxi(:,3,1) 339 & - dxdxi(:,2,1) * dxdxi(:,3,3)) * tmp 340 dxidx(:,2,2) = (dxdxi(:,1,1) * dxdxi(:,3,3) 341 & - dxdxi(:,3,1) * dxdxi(:,1,3)) * tmp 342 dxidx(:,2,3) = (dxdxi(:,2,1) * dxdxi(:,1,3) 343 & - dxdxi(:,1,1) * dxdxi(:,2,3)) * tmp 344 dxidx(:,3,1) = (dxdxi(:,2,1) * dxdxi(:,3,2) 345 & - dxdxi(:,2,2) * dxdxi(:,3,1)) * tmp 346 dxidx(:,3,2) = (dxdxi(:,3,1) * dxdxi(:,1,2) 347 & - dxdxi(:,1,1) * dxdxi(:,3,2)) * tmp 348 dxidx(:,3,3) = (dxdxi(:,1,1) * dxdxi(:,2,2) 349 & - dxdxi(:,1,2) * dxdxi(:,2,1)) * tmp 350c 351 WdetJ = Qwt(lcsyst,intp)/ tmp 352 353c 354c.... ---------------------> Global Gradients <----------------------- 355c 356 gradient = zero 357c 358c 359 do n = 1, nshl 360c 361c.... compute the global gradient of shape-function 362c 363c ! N_{a,x_i}= N_{a,xi_i} xi_{i,x_j} 364c 365 shg(:,n,1) = shgl(:,1,n) * dxidx(:,1,1) + 366 & shgl(:,2,n) * dxidx(:,2,1) + 367 & shgl(:,3,n) * dxidx(:,3,1) 368 shg(:,n,2) = shgl(:,1,n) * dxidx(:,1,2) + 369 & shgl(:,2,n) * dxidx(:,2,2) + 370 & shgl(:,3,n) * dxidx(:,3,2) 371 shg(:,n,3) = shgl(:,1,n) * dxidx(:,1,3) + 372 & shgl(:,2,n) * dxidx(:,2,3) + 373 & shgl(:,3,n) * dxidx(:,3,3) 374c 375c 376c Y_{,x_i}=SUM_{a=1}^nenl (N_{a,x_i}(int) Ya) 377c 378 do i = 1, 3 379 do j = 1, vsize 380 k = (i-1)*vsize+j 381 gradient(:,k) = gradient(:,k) + shg(:,n,i)*vector(:,n,j) 382 end do 383 end do 384 385 end do 386 387c 388c.... return 389c 390 return 391 end 392