1 module rlssave 2 3 real*8, allocatable :: rls(:,:) 4 5 end module 6 7c----------------------------------------------------------------------------- 8 9 10 11c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 12c||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| 13c||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| 14 15 16 17 18 subroutine setrls 19 20 use rlssave 21 22 include "common.h" 23 24 allocate ( rls(nshg,6) ) 25 26 return 27 28 end 29 30 31 32 33c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 34c||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| 35c||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| 36 37 38 39 40 subroutine bardmc (y, shgl, shp, 41 & iper, ilwork, 42 & nsons, ifath, x) 43 44 use pointer_data 45 use rlssave ! rls(nshg,6) is retrieved and modified here. 46 47 use quadfilt ! This module gives us shglf(maxtp,nsd,maxsh,ngaussf), 48c shpf(maxtp,maxsh,ngaussf), and Qwtf(maxtp,ngaussf). 49c Shpf and shglf are the shape funciotns and their 50c gradient evaluated using the quadrature rule desired 51c for computing the dmod. Qwtf contains the weights of the 52c quad. points. 53 54 include "common.h" 55 include "mpif.h" 56 include "auxmpi.h" 57 58c 59 dimension fres(nshg,33), fwr(nshg), 60 & strnrm(nshg), cdelsq(nshg), 61 & xnum(nshg), xden(nshg), 62 & xmij(nshg,6), xlij(nshg,6), 63 & xrij(nshg,6), xhij(nshg,6), 64 & xnude(nfath,2), xnuder(nfath,2), 65 & nsons(nfath), 66 & strl(numel,maxnint), 67 & y(nshg,5), yold(nshg,5), 68 & ifath(nshg), iper(nshg), 69 & ilwork(nlwork), 70 & x(numnp,3), 71 & shgl(MAXTOP,nsd,maxsh,MAXQPT), shp(MAXTOP,maxsh,MAXQPT), 72 & xnutf(nfath), 73 & hfres(nshg,11), 74 & rho(nshg), 75 & eps(3), epst(3), 76 & tracerls(nshg) 77 78c 79c 80c setup the weights for time averaging of cdelsq (now in quadfilt module) 81c 82 denom=max((one*lstep),one) 83 if(dtavei.lt.0) then 84 wcur=one/denom 85 else 86 wcur=dtavei 87 endif 88 whist=1.0-wcur 89 90 fres = zero 91 hfres = zero 92 93 yold(:,1)=y(:,4) 94 yold(:,2:4)=y(:,1:3) 95 yold(:,5) = y(:,5) 96 97 if(matflg(1,1).eq.0) then ! compressible 98 rho(:) = yold(:,1) / (Rgas * yold(:,5)) ! get density at nodes. 99 else 100 rho(:) = one ! unit density 101 endif 102 103c 104c hack in an interesting velocity field (uncomment to test dmod) 105c 106c yold(:,5) = 1.0 ! Debugging 107c yold(:,2) = 2.0*x(:,1) - 3*x(:,2) 108c yold(:,3) = 3.0*x(:,1) + 4.0*x(:,2) 109c yold(:,4) = 4.0*x(:,1) + x(:,2) + x(:,3) 110c yold(:,1) = Rgas * yold(:,5) ! Necessary to make model suitable 111c suitable for the 112 113 do iblk = 1,nelblk 114 lcsyst = lcblk(3,iblk) 115 iel = lcblk(1,iblk) !Element number where this block begins 116 npro = lcblk(1,iblk+1) - iel 117 lelCat = lcblk(2,iblk) 118 nenl = lcblk(5,iblk) 119 nshl = lcblk(10,iblk) 120 inum = iel + npro - 1 121 122 ngauss = nint(lcsyst) 123 ngaussf = nintf(lcsyst) 124 125 call hfilter (yold, x, mien(iblk)%p, hfres, 126 & shglf(lcsyst,:,1:nshl,:), 127 & shpf(lcsyst,1:nshl,:),Qwtf(lcsyst,1:ngaussf)) 128 129 enddo 130 131 if(numpe>1) call commu (hfres, ilwork, 11, 'in ') 132c 133c... account for periodicity in filtered variables 134c 135 do j = 1,nshg ! Add on-processor slave contribution to masters 136 i = iper(j) 137 if (i .ne. j) then 138 hfres(i,:) = hfres(i,:) + hfres(j,:) 139 endif 140 enddo 141 do j = 1,nshg ! Set on-processor slaves to be the same as masters 142 i = iper(j) 143 if (i .ne. j) then 144 hfres(j,:) = hfres(i,:) 145 endif 146 enddo 147 148c... Set off-processor slaves to be the same as their masters 149 150 if(numpe>1) call commu (hfres, ilwork, 11, 'out') 151 152 153 hfres(:,5) = one / hfres(:,5) ! one/(volume of hat filter kernel) 154 155 do j = 1, 4 156 hfres(:,j) = hfres(:,j) * hfres(:,5) 157 enddo 158 do j = 6, 11 159 hfres(:,j) = hfres(:,j) * hfres(:,5) 160 enddo 161 162c Uncomment to test dmod 163 164c do j = 1, 11 165c do i = 1, numnp 166c hfres(i,j) = hfres(14,j) 167c enddo 168c enddo 169 170c End of dmod test 171 172c... Form the resolved Leonard stress (rls) at each node. 173 174 rls(:,1) = hfres(:,6) - hfres(:,1)*hfres(:,1)/rho(:) 175 rls(:,2) = hfres(:,7) - hfres(:,2)*hfres(:,2)/rho(:) 176 rls(:,3) = hfres(:,8) - hfres(:,3)*hfres(:,3)/rho(:) 177 178 tracerls=(rls(:,1)+rls(:,2)+rls(:,3))/three 179 180 rls(:,1) = rls(:,1) - tracerls 181 rls(:,2) = rls(:,2) - tracerls 182 rls(:,3) = rls(:,3) - tracerls 183 184 rls(:,4) = hfres(:,9) - hfres(:,1)*hfres(:,2)/rho(:) 185 rls(:,5) = hfres(:,10) - hfres(:,1)*hfres(:,3)/rho(:) 186 rls(:,6) = hfres(:,11) - hfres(:,2)*hfres(:,3)/rho(:) 187 188c rls(:,1) = 4.0d0 ! Debugging 189c rls(:,2) = 5.0d0 190c rls(:,3) = 1.0d0 191c rls(:,4) = 2.1d0 192c rls(:,5) = 1.5d0 193c rls(:,6) = 3.0d0 194 195c... Stored the resolved Leonard stresses in a module. 196 197c... Compute dissipation of resolved kinetic energy (u_{i} u_{i}) 198c... due to molecular stresses (epst(3)), to modeled residual stresses 199c... (epst(1)), and to resolved Leonard stresses (epst(2)). 200 201c if (lstep .gt. 0) then 202 203c eps = zero 204c do iblk = 1,nelblk 205c lcsyst = lcblk(3,iblk) 206c iel = lcblk(1,iblk) !Element number where this block begins 207c npro = lcblk(1,iblk+1) - iel 208c lelCat = lcblk(2,iblk) 209c nenl = lcblk(5,iblk) 210c nshl = lcblk(10,iblk) 211c inum = iel + npro - 1 212 213c ngauss = nint(lcsyst) 214 215c call disPtion ( yold, x, mien(iblk)%p, 216c & shgl(lcsyst,:,1:nshl,:), shp(lcsyst,1:nshl,:), 217c & eps, mxmudmi(iblk)%p ) 218 219c enddo 220 221c if(numpe.gt.1)then 222c call drvAllreduce( eps, epst, 3 ) 223c else 224c epst = eps 225c endif 226 227c if (myrank .eq. master) then 228c write(22,*) one*(lstep+1), epst(1), epst(2) 229c call flush(22) 230c endif 231 232c endif 233 234 235c... Done w/ h-filtering. Begin 2h-filtering. 236 237 do iblk = 1,nelblk 238 lcsyst = lcblk(3,iblk) 239 iel = lcblk(1,iblk) !Element number where this block begins 240 npro = lcblk(1,iblk+1) - iel 241 lelCat = lcblk(2,iblk) 242 nenl = lcblk(5,iblk) 243 nshl = lcblk(10,iblk) 244 inum = iel + npro - 1 245 246 ngauss = nint(lcsyst) 247 ngaussf = nintf(lcsyst) 248 249 call twohfilter (yold, x, strl(iel:inum,:), mien(iblk)%p, 250 & fres, hfres, shgl(lcsyst,:,1:nshl,:), 251 & shp(lcsyst,1:nshl,:),Qwtf(lcsyst,1:ngaussf)) 252 253 enddo 254c 255 256 do iblk = 1,nelblk 257 lcsyst = lcblk(3,iblk) 258 iel = lcblk(1,iblk) !Element number where this block begins 259 npro = lcblk(1,iblk+1) - iel 260 lelCat = lcblk(2,iblk) 261 nenl = lcblk(5,iblk) 262 nshl = lcblk(10,iblk) 263 inum = iel + npro - 1 264 265 ngauss = nint(lcsyst) 266 ngaussf = nintf(lcsyst) 267 if (ngauss .ne. ngaussf) then 268 call getstrl (yold, x, mien(iblk)%p, 269 & strl(iel:inum,:), shgl(lcsyst,:,1:nshl,:), 270 & shp(lcsyst,1:nshl,:)) 271 endif 272 enddo 273c 274c 275C must fix for abc and dynamic model 276c if(any(btest(iBC,12))) !are there any axisym bc's 277c & call rotabc(res, iBC, 'in ') 278c 279 if(numpe>1) call commu (fres, ilwork, 33, 'in ') 280c 281c account for periodicity in filtered variables 282c 283 do j = 1,nshg 284 i = iper(j) 285 if (i .ne. j) then 286 fres(i,:) = fres(i,:) + fres(j,:) 287 fres(j,:) = fres(i,:) 288 endif 289 enddo 290 291 if (nfath .eq. 0) then 292 fres = fres(iper,:) 293 if(numpe>1) call commu (fres, ilwork, 33, 'out') 294c again we will need a call 295c if(any(btest(iBC,12))) !are there any axisym bc's 296c & rotabc(y, iBC, 'out') 297 endif 298 299 fres(:,23) = one / fres(:,23) 300 do j = 1,22 301 fres(:,j) = fres(:,j) * fres(:,23) 302 enddo 303c fres(:,24) = fres(:,24) * fres(:,23) 304 do j = 25,33 305 fres(:,j) = fres(:,j) * fres(:,23) 306 enddo 307c 308c.....at this point fres is really all of our filtered quantities 309c at the nodes 310c 311 312 strnrm = sqrt( 313 & two * (fres(:,10)**2 + fres(:,11)**2 + fres(:,12)**2) 314 & + four * ( fres(:,13)**2 + fres(:,14)**2 + fres(:,15)**2 ) ) 315 316 fwr = fwr1 * fres(:,22) * strnrm 317 318 if(matflg(1,1).eq.0) then ! compressible 319 320 xmij(:,1) = -fwr 321 & * pt33 * (two*fres(:,10) - fres(:,11) - fres(:,12)) 322 & + pt33 * (two*fres(:,16) - fres(:,17) - fres(:,18)) 323c 324 xmij(:,2) = -fwr 325 & * pt33 * (two*fres(:,11) - fres(:,10) - fres(:,12)) 326 & + pt33 * (two*fres(:,17) - fres(:,16) - fres(:,18)) 327c 328 xmij(:,3) = -fwr 329 & * pt33 * (two*fres(:,12) - fres(:,10) - fres(:,11)) 330 & + pt33 * (two*fres(:,18) - fres(:,16) - fres(:,17)) 331 332 else 333 334 xmij(:,1) = -fwr 335 & * fres(:,10) + fres(:,16) 336 xmij(:,2) = -fwr 337 & * fres(:,11) + fres(:,17) 338 xmij(:,3) = -fwr 339 & * fres(:,12) + fres(:,18) 340 341 endif 342 343 xmij(:,4) = -fwr * fres(:,13) + fres(:,19) 344 xmij(:,5) = -fwr * fres(:,14) + fres(:,20) 345 xmij(:,6) = -fwr * fres(:,15) + fres(:,21) 346 347 fres(:,22) = one / fres(:,22) 348 349 xlij(:,1) = fres(:,4) - fres(:,1) * fres(:,1) * fres(:,22) 350 xlij(:,2) = fres(:,5) - fres(:,2) * fres(:,2) * fres(:,22) 351 xlij(:,3) = fres(:,6) - fres(:,3) * fres(:,3) * fres(:,22) 352 xlij(:,4) = fres(:,7) - fres(:,1) * fres(:,2) * fres(:,22) 353 xlij(:,5) = fres(:,8) - fres(:,1) * fres(:,3) * fres(:,22) 354 xlij(:,6) = fres(:,9) - fres(:,2) * fres(:,3) * fres(:,22) 355 356 xhij(:,1) = fres(:,28) - fres(:,25)*fres(:,25) 357 xhij(:,2) = fres(:,29) - fres(:,26)*fres(:,26) 358 xhij(:,3) = fres(:,30) - fres(:,27)*fres(:,27) 359 xhij(:,4) = fres(:,31) - fres(:,25)*fres(:,26) 360 xhij(:,5) = fres(:,32) - fres(:,25)*fres(:,27) 361 xhij(:,6) = fres(:,33) - fres(:,26)*fres(:,27) 362 363 xrij(:,1) = xlij(:,1) - xhij(:,1) 364 xrij(:,2) = xlij(:,2) - xhij(:,2) 365 xrij(:,3) = xlij(:,3) - xhij(:,3) 366 xrij(:,4) = xlij(:,4) - xhij(:,4) 367 xrij(:,5) = xlij(:,5) - xhij(:,5) 368 xrij(:,6) = xlij(:,6) - xhij(:,6) 369 370 xnum = xrij(:,1) * xmij(:,1) + xrij(:,2) * xmij(:,2) 371 & + xrij(:,3) * xmij(:,3) 372 & + two * (xrij(:,4) * xmij(:,4) + xrij(:,5) * xmij(:,5) 373 & + xrij(:,6) * xmij(:,6)) 374 xden = xmij(:,1) * xmij(:,1) + xmij(:,2) * xmij(:,2) 375 & + xmij(:,3) * xmij(:,3) 376 & + two * (xmij(:,4) * xmij(:,4) + xmij(:,5) * xmij(:,5) 377 & + xmij(:,6) * xmij(:,6)) 378 xden = two * xden 379 380c zero on processor periodic nodes so that they will not be added twice 381 do j = 1,numnp 382 i = iper(j) 383 if (i .ne. j) then 384 xnum(j) = zero 385 xden(j) = zero 386 endif 387 enddo 388 389 if (numpe.gt.1 .and. nsons(1).gt.1) then 390 391 numtask = ilwork(1) 392 itkbeg = 1 393 394c zero the nodes that are "solved" on the other processors 395 do itask = 1, numtask 396 397 iacc = ilwork (itkbeg + 2) 398 numseg = ilwork (itkbeg + 4) 399 400 if (iacc .eq. 0) then 401 do is = 1,numseg 402 isgbeg = ilwork (itkbeg + 3 + 2*is) 403 lenseg = ilwork (itkbeg + 4 + 2*is) 404 isgend = isgbeg + lenseg - 1 405 xnum(isgbeg:isgend) = zero 406 xden(isgbeg:isgend) = zero 407 enddo 408 endif 409 410 itkbeg = itkbeg + 4 + 2*numseg 411 412 enddo 413 414 endif 415c 416c Description of arrays. Each processor has an array of length equal 417c to the total number of fathers times 2 xnude(nfathers,2). One to collect 418c the numerator and one to collect the denominator. There is also an array 419c of length nshg on each processor which tells the father number of each 420c on processor node, ifath(nnshg). Finally, there is an arry of length 421c nfathers to tell the total (on all processors combined) number of sons 422c for each father. 423c 424c Now loop over nodes and accumlate the numerator and the denominator 425c to the father nodes. Only on processor addition at this point. 426c Note that serrogate fathers are collect some for the case where some 427c sons are on another processor 428c 429 if(nsons(1) .gt. 1) then 430 xnude = zero 431 do i = 1,nshg 432 xnude(ifath(i),1) = xnude(ifath(i),1) + xnum(i) 433 xnude(ifath(i),2) = xnude(ifath(i),2) + xden(i) 434 enddo 435 else ! no homogeneity assumed 436 xnude(:,1)=xnum(:) 437 xnude(:,2)=xden(:) 438 endif 439 440c 441c Now the true fathers and serrogates combine results and update 442c each other. 443c 444c ONLY DO THIS IF 1) multiprocessors AND 2) homogenous directions exist 445 446 if(numpe .gt. 1 .and. nsons(1).ne.1 )then 447 call drvAllreduce(xnude, xnuder,2*nfath) 448c 449c xnude is the sum of the sons for each father on this processor 450c 451c xnuder is the sum of the sons for each father on all processor combined 452c (the same as if we had not partitioned the mesh for each processor) 453c 454c For each father we have precomputed the number of sons (including 455c the sons off processor). 456c 457c Now divide by number of sons to get the average (not really necessary 458c for dynamic model since ratio will cancel nsons at each father) 459c 460c xnuder(:,1) = xnuder(:,1) ! / nsons(:) 461c xnuder(:,2) = xnuder(:,2) ! / nsons(:) 462c 463c the next line is c \Delta^2 464c 465 numNden(:,1) = whist*numNden(:,1)+wcur*xnuder(ifath(:),1) 466 numNden(:,2) = whist*numNden(:,2)+wcur*xnuder(ifath(:),2) 467 cdelsq(:) = numNden(:,1) / (numNden(:,2) + 1.d-09) 468 else 469c 470c the next line is c \Delta^2, not nu_T but we want to save the 471c memory 472c 473c$$$ 474c$$$ write(400+myrank,555) (numNden(j*500,1),j=1,5) 475c$$$ write(410+myrank,555) (numNden(j*500,2),j=1,5) 476c$$$ write(500+myrank,555) (numNden(j+500,1),j=1,5) 477c$$$ write(510+myrank,555) (numNden(j+500,2),j=1,5) 478c$$$ 479c$$$ write(430+myrank,555) (xnude(j*500,1),j=1,5) 480c$$$ write(440+myrank,555) (xnude(j*500,2),j=1,5) 481c$$$ write(530+myrank,555) (xnude(j+500,1),j=1,5) 482c$$$ write(540+myrank,555) (xnude(j+500,2),j=1,5) 483 484 numNden(:,1) = whist*numNden(:,1)+wcur*xnude(ifath(:),1) 485 numNden(:,2) = whist*numNden(:,2)+wcur*xnude(ifath(:),2) 486 cdelsq(:) = numNden(:,1) / (numNden(:,2) + 1.d-09) 487 488c 489c to get started we hold cdelsq fixed 490c 491c cdelsq(:) = 2.27e-4 492 493 494 495c$$$ write(450+myrank,555) (cdelsq(j*500),j=1,5) 496c$$$ write(460+myrank,555) (y(j*500,1),j=1,5) 497c$$$ write(470+myrank,555) (y(j*500,2),j=1,5) 498c$$$ write(480+myrank,555) (y(j*500,3),j=1,5) 499c$$$ write(490+myrank,555) (strnrm(j*500),j=1,5) 500c$$$ 501c$$$ write(550+myrank,555) (cdelsq(j+500),j=1,5) 502c$$$ write(560+myrank,555) (y(j+500,1),j=1,5) 503c$$$ write(570+myrank,555) (y(j+500,2),j=1,5) 504c$$$ write(580+myrank,555) (y(j+500,3),j=1,5) 505c$$$ write(590+myrank,555) (strnrm(j+500),j=1,5) 506c$$$ 507c$$$ call flush(400+myrank) 508c$$$ call flush(410+myrank) 509c$$$ call flush(430+myrank) 510c$$$ call flush(440+myrank) 511c$$$ call flush(450+myrank) 512c$$$ call flush(460+myrank) 513c$$$ call flush(470+myrank) 514c$$$ call flush(480+myrank) 515c$$$ call flush(490+myrank) 516c$$$ call flush(500+myrank) 517c$$$ call flush(510+myrank) 518c$$$ call flush(530+myrank) 519c$$$ call flush(540+myrank) 520c$$$ call flush(550+myrank) 521c$$$ call flush(560+myrank) 522c$$$ call flush(570+myrank) 523c$$$ call flush(580+myrank) 524c$$$ call flush(590+myrank) 525 endif 526 555 format(5(2x,e14.7)) 527 528c $$$$$$$$$$$$$$$$$$$$$$$$$$$ 529 tmp1 = MINVAL(cdelsq) 530 tmp2 = MAXVAL(cdelsq) 531 if(numpe>1) then 532 call MPI_REDUCE (tmp1, tmp3, 1,MPI_DOUBLE_PRECISION, 533 & MPI_MIN, master, MPI_COMM_WORLD, ierr) 534 call MPI_REDUCE (tmp2, tmp4, 1, MPI_DOUBLE_PRECISION, 535 & MPI_MAX, master, MPI_COMM_WORLD, ierr) 536 tmp1=tmp3 537 tmp2=tmp4 538 endif 539 if (myrank .EQ. master) then !print CDelta^2 range 540 write(34,*)lstep,tmp1,tmp2 541 call flush(34) 542 endif 543c $$$$$$$$$$$$$$$$$$$$$$$$$$$ 544 545 if (myrank .eq. master) then 546c write(*,*)'xnut=',sum(cdelsq)/nshg 547 write(*,*)'cdelsq=',cdelsq(1), cdelsq(2) 548 endif 549 550 cdeslq = zero ! Debugging 551 do iblk = 1,nelblk 552 lcsyst = lcblk(3,iblk) 553 iel = lcblk(1,iblk) 554 npro = lcblk(1,iblk+1) - iel 555 lelCat = lcblk(2,iblk) 556 inum = iel + npro - 1 557 558 ngauss = nint(lcsyst) 559 560 call scatnu (mien(iblk)%p, strl(iel:inum,:), 561 & mxmudmi(iblk)%p,cdelsq,shp(lcsyst,1:nshl,:)) 562 enddo 563c $$$$$$$$$$$$$$$$$$$$$$$$$$$ 564c$$$ tmp1 = MINVAL(xmudmi) 565c$$$ tmp2 = MAXVAL(xmudmi) 566c$$$ if(numpe>1) then 567c$$$ call MPI_REDUCE (tmp1, tmp3, 1, MPI_DOUBLE_PRECISION, 568c$$$ & MPI_MIN, master, MPI_COMM_WORLD, ierr) 569c$$$ call MPI_REDUCE (tmp2, tmp4, 1, MPI_DOUBLE_PRECISION, 570c$$$ & MPI_MAX, master, MPI_COMM_WORLD, ierr) 571c$$$ tmp1=tmp3 572c$$$ tmp2=tmp4 573c$$$ endif 574c$$$ if (myrank .EQ. master) then 575c$$$ write(35,*) lstep,tmp1,tmp2 576c$$$ call flush(35) 577c$$$ endif 578c $$$$$$$$$$$$$$$$$$$$$$$$$$$ 579 580c 581c if flag set, write a restart file with info (reuse xmij's memory) 582c 583 if(irs.eq.11) then 584 lstep=999 585 xmij(:,1)=xnum(:) 586 xmij(:,2)=xden(:) 587 xmij(:,3)=cdelsq(:) 588 xmij(:,5)=xlij(:,4) !leave M_{12} in 4 and put L_{12} here 589 call restar('out ',xmij,xlij) !also dump all of L_{ij} in ac 590 stop 591 endif 592c 593c local clipping moved to scatnu with the creation of mxmudmi pointers 594c 595c$$$ rmu=datmat(1,2,1) 596c$$$ xmudmi=min(xmudmi,1000.0*rmu) !don't let it get larger than 1000 mu 597c$$$ xmudmi=max(xmudmi, -rmu) ! don't let (xmudmi + mu) < 0 598c stop !uncomment to test dmod 599c 600 601 602c write out the nodal values of xnut (estimate since we don't calc strain 603c there and must use the filtered strain). 604c 605 606c if ((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) then 607c 608c collect the average strain into xnude(2) 609c 610c xnude(:,2) = zero 611c do i = 1,numnp 612c xnude(ifath(i),2) = xnude(ifath(i),2) + strnrm(i) 613c enddo 614 615c 616c get the instantaneous cdelsq of the fathers 617c 618c xnuder(:,1)=xnuder(:,1)/(xnuder(:,2)+1.0e-9) 619 620c if(numpe .gt. 1 .and. nsons(1).ne.1 )then 621c call drvAllreduce(xnude(:,2), xnuder(:,2),nfath) 622c else 623c xnuder=xnude 624c endif 625c 626c nut= cdelsq * |S| 627c 628c xnutf=xnuder(:,1)*xnuder(:,2)/nsons 629c 630c collect the x and y coords into xnude 631c 632c xnude = zero 633c do i = 1,numnp 634c xnude(ifath(i),1) = xnude(ifath(i),1) + x(i,1) 635c xnude(ifath(i),2) = xnude(ifath(i),2) + x(i,2) 636c enddo 637 638c if(numpe .gt. 1 .and. nsons(1).ne.1 )then 639c call drvAllreduce(xnude, xnuder,2*nfath) 640c xnuder(:,1)=xnuder(:,1)/nsons(:) 641c xnuder(:,2)=xnuder(:,2)/nsons(:) 642c else 643c xnuder=xnude 644c endif 645c 646c xnude is the sum of the sons for each father on this processor 647c 648c if((myrank.eq.master)) then 649c do i=1,nfath ! cdelsq * |S| 650c write(444,*) xnuder(i,1),xnuder(i,2),xnutf(i) 651c enddo 652c call flush(444) 653c endif 654c endif 655 656 return 657 end 658 659 660 661 662 663 664 665c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 666c||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| 667c||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| 668 669 670 671 672 subroutine hfilter (y, x, ien, hfres, shgl, shp, Qwtf) 673 674 include "common.h" 675 676 dimension y(nshg,5), hfres(nshg,11) 677 dimension x(numnp,3), xl(npro,nenl,3) 678 dimension ien(npro,nshl), yl(npro,nshl,nflow), 679 & fresl(npro,nshl,11), WdetJ(npro), 680 & u1(npro), u2(npro), 681 & u3(npro), rho(npro), 682 & dxdxi(npro,nsd,nsd), dxidx(npro,nsd,nsd), 683 & shgl(nsd,nshl,maxsh), shg(npro,nshl,nsd), 684 & shp(nshl,maxsh), Qwtf(ngaussf), 685 & fresli(npro,nshl,11) 686 687 dimension tmp(npro) 688 689 call local (y, yl, ien, 5, 'gather ') 690 call localx (x, xl, ien, 3, 'gather ') 691c 692c 693 if(matflg(1,1).eq.0) then ! compressible 694 yl (:,:,1) = yl(:,:,1) / (Rgas * yl(:,:,5)) !get density 695 else 696 yl (:,:,1) = one 697 endif 698c 699 fresl = zero 700c 701 do intp = 1, ngaussf 702 703c calculate the metrics 704c 705c 706c.... ---------------------> Element Metrics <----------------------- 707c 708c.... compute the deformation gradient 709c 710 dxdxi = zero 711c 712 do n = 1, nenl 713 dxdxi(:,1,1) = dxdxi(:,1,1) + xl(:,n,1) * shgl(1,n,intp) 714 dxdxi(:,1,2) = dxdxi(:,1,2) + xl(:,n,1) * shgl(2,n,intp) 715 dxdxi(:,1,3) = dxdxi(:,1,3) + xl(:,n,1) * shgl(3,n,intp) 716 dxdxi(:,2,1) = dxdxi(:,2,1) + xl(:,n,2) * shgl(1,n,intp) 717 dxdxi(:,2,2) = dxdxi(:,2,2) + xl(:,n,2) * shgl(2,n,intp) 718 dxdxi(:,2,3) = dxdxi(:,2,3) + xl(:,n,2) * shgl(3,n,intp) 719 dxdxi(:,3,1) = dxdxi(:,3,1) + xl(:,n,3) * shgl(1,n,intp) 720 dxdxi(:,3,2) = dxdxi(:,3,2) + xl(:,n,3) * shgl(2,n,intp) 721 dxdxi(:,3,3) = dxdxi(:,3,3) + xl(:,n,3) * shgl(3,n,intp) 722 enddo 723c 724c.... compute the inverse of deformation gradient 725c 726 dxidx(:,1,1) = dxdxi(:,2,2) * dxdxi(:,3,3) 727 & - dxdxi(:,3,2) * dxdxi(:,2,3) 728 dxidx(:,1,2) = dxdxi(:,3,2) * dxdxi(:,1,3) 729 & - dxdxi(:,1,2) * dxdxi(:,3,3) 730 dxidx(:,1,3) = dxdxi(:,1,2) * dxdxi(:,2,3) 731 & - dxdxi(:,1,3) * dxdxi(:,2,2) 732 tmp = one / ( dxidx(:,1,1) * dxdxi(:,1,1) 733 & + dxidx(:,1,2) * dxdxi(:,2,1) 734 & + dxidx(:,1,3) * dxdxi(:,3,1) ) 735 dxidx(:,1,1) = dxidx(:,1,1) * tmp 736 dxidx(:,1,2) = dxidx(:,1,2) * tmp 737 dxidx(:,1,3) = dxidx(:,1,3) * tmp 738 dxidx(:,2,1) = (dxdxi(:,2,3) * dxdxi(:,3,1) 739 & - dxdxi(:,2,1) * dxdxi(:,3,3)) * tmp 740 dxidx(:,2,2) = (dxdxi(:,1,1) * dxdxi(:,3,3) 741 & - dxdxi(:,3,1) * dxdxi(:,1,3)) * tmp 742 dxidx(:,2,3) = (dxdxi(:,2,1) * dxdxi(:,1,3) 743 & - dxdxi(:,1,1) * dxdxi(:,2,3)) * tmp 744 dxidx(:,3,1) = (dxdxi(:,2,1) * dxdxi(:,3,2) 745 & - dxdxi(:,2,2) * dxdxi(:,3,1)) * tmp 746 dxidx(:,3,2) = (dxdxi(:,3,1) * dxdxi(:,1,2) 747 & - dxdxi(:,1,1) * dxdxi(:,3,2)) * tmp 748 dxidx(:,3,3) = (dxdxi(:,1,1) * dxdxi(:,2,2) 749 & - dxdxi(:,1,2) * dxdxi(:,2,1)) * tmp 750c 751 wght = Qwtf(intp) 752 WdetJ(:) = wght / tmp(:) 753 754 rho = zero 755 do i=1,nshl 756 rho = rho+shp(i,intp)*yl(:,i,1) !density at qpt 757 enddo 758 759 rho = rho * WdetJ !rho * WdetJ 760 761 u1=zero 762 u2=zero 763 u3=zero 764 do i=1,nshl ! velocities at qpt 765 u1 = u1 + shp(i,intp)*yl(:,i,2) 766 u2 = u2 + shp(i,intp)*yl(:,i,3) 767 u3 = u3 + shp(i,intp)*yl(:,i,4) 768 enddo 769 770 do i = 1, nshl 771 772 773 fresli(:,i,1) = rho * u1 * shp(i,intp) !rho * u1 * WdetJ 774 fresli(:,i,2) = rho * u2 * shp(i,intp) !rho * u2 * WdetJ 775 fresli(:,i,3) = rho * u3 * shp(i,intp) !rho * u3 * WdetJ 776 fresli(:,i,4) = rho * shp(i,intp) !rho * WdetJ 777 fresli(:,i,5) = WdetJ*shp(i,intp) !Integral of filter kernel 778c over the element 779 fresli(:,i,6) = rho * u1 * u1 * shp(i,intp) !rho * u1 * u1 * WdetJ 780 fresli(:,i,7) = rho * u2 * u2 * shp(i,intp) !rho * u2 * u2 * WdetJ 781 fresli(:,i,8) = rho * u3 * u3 * shp(i,intp) !rho * u3 * u3 * WdetJ 782 fresli(:,i,9) = rho * u1 * u2 * shp(i,intp) !rho * u1 * u2 * WdetJ 783 fresli(:,i,10)= rho * u1 * u3 * shp(i,intp) !rho * u1 * u3 * WdetJ 784 fresli(:,i,11)= rho * u2 * u3 * shp(i,intp) !rho * u2 * u3 * WdetJ 785 786 787 do j = 1, 11 788 fresl(:,i,j) = fresl(:,i,j) + fresli(:,i,j) 789 enddo 790 791 enddo !end loop over element nodes 792 793 enddo !end of loop over integration points 794 795c do j = 1, n 796c do i = 1, nenl 797c do nel = 1, npro 798c hfres(ien(nel,i),j) = hfres(ien(nel,i),j) + fresl(nel,i,j) 799c enddo 800c enddo 801c enddo 802 803 call local (hfres, fresl, ien, 11, 'scatter ') 804 805 return 806 end 807 808 809 810 811 812 813 814c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 815c||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| 816c||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| 817 818 819 820 821 822 823 824 825 subroutine twohfilter (y, x, strnrm, ien, fres, 826 & hfres, shgl, shp, Qwtf) 827 828 include "common.h" 829 830 dimension y(nshg,ndof), fres(nshg,33) 831 dimension x(numnp,nsd), xl(npro,nenl,nsd) 832 dimension ien(npro,nshl), yl(npro,nshl,nflow), 833 & fresl(npro,33), WdetJ(npro), 834 & u1(npro), u2(npro), 835 & u3(npro), dxdxi(npro,nsd,nsd), 836 & strnrm(npro,ngauss), dxidx(npro,nsd,nsd), 837 & shgl(nsd,nshl,maxsh), shg(npro,nshl,nsd), 838 & shp(nshl,maxsh), 839 & fresli(npro,33), Qwtf(ngaussf), 840 & hfres(nshg,11), hfresl(npro,nshl,11) 841 842 dimension tmp(npro) 843 844 call local (y, yl, ien, 5, 'gather ') 845 call localx (x, xl, ien, 3, 'gather ') 846 call local (hfres, hfresl, ien, 11, 'gather ') 847c 848 if(matflg(1,1).eq.0) then ! compressible 849 yl (:,:,1) = yl(:,:,1) / (Rgas * yl(:,:,5)) !get density 850 else 851 yl (:,:,1) = one 852 endif 853c 854 fresl = zero 855 856 do intp = 1, ngauss 857 858 859c calculate the metrics 860c 861c 862c.... ---------------------> Element Metrics <----------------------- 863c 864c.... compute the deformation gradient 865c 866 dxdxi = zero 867c 868 do n = 1, nenl 869 dxdxi(:,1,1) = dxdxi(:,1,1) + xl(:,n,1) * shgl(1,n,intp) 870 dxdxi(:,1,2) = dxdxi(:,1,2) + xl(:,n,1) * shgl(2,n,intp) 871 dxdxi(:,1,3) = dxdxi(:,1,3) + xl(:,n,1) * shgl(3,n,intp) 872 dxdxi(:,2,1) = dxdxi(:,2,1) + xl(:,n,2) * shgl(1,n,intp) 873 dxdxi(:,2,2) = dxdxi(:,2,2) + xl(:,n,2) * shgl(2,n,intp) 874 dxdxi(:,2,3) = dxdxi(:,2,3) + xl(:,n,2) * shgl(3,n,intp) 875 dxdxi(:,3,1) = dxdxi(:,3,1) + xl(:,n,3) * shgl(1,n,intp) 876 dxdxi(:,3,2) = dxdxi(:,3,2) + xl(:,n,3) * shgl(2,n,intp) 877 dxdxi(:,3,3) = dxdxi(:,3,3) + xl(:,n,3) * shgl(3,n,intp) 878 enddo 879c 880c.... compute the inverse of deformation gradient 881c 882 dxidx(:,1,1) = dxdxi(:,2,2) * dxdxi(:,3,3) 883 & - dxdxi(:,3,2) * dxdxi(:,2,3) 884 dxidx(:,1,2) = dxdxi(:,3,2) * dxdxi(:,1,3) 885 & - dxdxi(:,1,2) * dxdxi(:,3,3) 886 dxidx(:,1,3) = dxdxi(:,1,2) * dxdxi(:,2,3) 887 & - dxdxi(:,1,3) * dxdxi(:,2,2) 888 tmp = one / ( dxidx(:,1,1) * dxdxi(:,1,1) 889 & + dxidx(:,1,2) * dxdxi(:,2,1) 890 & + dxidx(:,1,3) * dxdxi(:,3,1) ) 891 dxidx(:,1,1) = dxidx(:,1,1) * tmp 892 dxidx(:,1,2) = dxidx(:,1,2) * tmp 893 dxidx(:,1,3) = dxidx(:,1,3) * tmp 894 dxidx(:,2,1) = (dxdxi(:,2,3) * dxdxi(:,3,1) 895 & - dxdxi(:,2,1) * dxdxi(:,3,3)) * tmp 896 dxidx(:,2,2) = (dxdxi(:,1,1) * dxdxi(:,3,3) 897 & - dxdxi(:,3,1) * dxdxi(:,1,3)) * tmp 898 dxidx(:,2,3) = (dxdxi(:,2,1) * dxdxi(:,1,3) 899 & - dxdxi(:,1,1) * dxdxi(:,2,3)) * tmp 900 dxidx(:,3,1) = (dxdxi(:,2,1) * dxdxi(:,3,2) 901 & - dxdxi(:,2,2) * dxdxi(:,3,1)) * tmp 902 dxidx(:,3,2) = (dxdxi(:,3,1) * dxdxi(:,1,2) 903 & - dxdxi(:,1,1) * dxdxi(:,3,2)) * tmp 904 dxidx(:,3,3) = (dxdxi(:,1,1) * dxdxi(:,2,2) 905 & - dxdxi(:,1,2) * dxdxi(:,2,1)) * tmp 906c 907 wght=Qwt(lcsyst,intp) ! may be different now 908c wght=Qwtf(intp) 909 WdetJ = wght / tmp 910c 911 fresli=zero 912c 913 do i=1,nshl 914 fresli(:,22) = fresli(:,22)+shp(i,intp)*yl(:,i,1) ! density at qpt 915 fresli(:,25) = fresli(:,25)+shp(i,intp)*hfresl(:,i,1) !bar(rho u1) 916 fresli(:,26) = fresli(:,26)+shp(i,intp)*hfresl(:,i,2) !bar(rho u2) 917 fresli(:,27) = fresli(:,27)+shp(i,intp)*hfresl(:,i,3) !bar(rho u3) 918 enddo 919c 920c...fresli(:,28) = WdetJ * bar(rho u1) * bar(rho u1) / rho 921c...fresli(:,29) = WdetJ * bar(rho u2) * bar(rho u2) / rho 922c...fresli(:,30) = WdetJ * bar(rho u3) * bar(rho u3) / rho 923c...fresli(:,31) = WdetJ * bar(rho u1) * bar(rho u2) / rho 924c...fresli(:,32) = WdetJ * bar(rho u1) * bar(rho u3) / rho 925c...fresli(:,33) = WdetJ * bar(rho u2) * bar(rho u3) / rho 926 927 fresli(:,28) = WdetJ(:) * fresli(:,25) * 928 & fresli(:,25) / fresli(:,22) 929 fresli(:,29) = WdetJ(:) * fresli(:,26) * 930 & fresli(:,26) / fresli(:,22) 931 fresli(:,30) = WdetJ(:) * fresli(:,27) * 932 & fresli(:,27) / fresli(:,22) 933 fresli(:,31) = WdetJ(:) * fresli(:,25) * 934 & fresli(:,26) / fresli(:,22) 935 fresli(:,32) = WdetJ(:) * fresli(:,25) * 936 & fresli(:,27) / fresli(:,22) 937 fresli(:,33) = WdetJ(:) * fresli(:,26) * 938 & fresli(:,27) / fresli(:,22) 939 940 fresli(:,25) = fresli(:,25) * WdetJ(:) ! WdetJ*bar(rho u1) 941 fresli(:,26) = fresli(:,26) * WdetJ(:) ! WdetJ*bar(rho u2) 942 fresli(:,27) = fresli(:,27) * WdetJ(:) ! WdetJ*bar(rho u3) 943c 944 do n = 1,nshl 945 shg(:,n,1) = (shgl(1,n,intp) * dxidx(:,1,1) 946 & + shgl(2,n,intp) * dxidx(:,2,1) 947 & + shgl(3,n,intp) * dxidx(:,3,1)) 948 shg(:,n,2) = (shgl(1,n,intp) * dxidx(:,1,2) 949 & + shgl(2,n,intp) * dxidx(:,2,2) 950 & + shgl(3,n,intp) * dxidx(:,3,2)) 951 shg(:,n,3) = (shgl(1,n,intp) * dxidx(:,1,3) 952 & + shgl(2,n,intp) * dxidx(:,2,3) 953 & + shgl(3,n,intp) * dxidx(:,3,3)) 954 enddo 955 956 do j=10,12 ! normal strainrate u_{i,i} no sum on i 957 ig=j-9 958 iv=j-8 959 do i=1,nshl 960 fresli(:,j) = fresli(:,j)+shg(:,i,ig)*yl(:,i,iv) 961 enddo 962 enddo 963 964c shear stresses NOTE there may be faster ways to do this 965c check agains CM5 code for speed WTP 966 967 do i=1,nshl 968 fresli(:,13) = fresli(:,13)+shg(:,i,2)*yl(:,i,2) 969 & +shg(:,i,1)*yl(:,i,3) 970 fresli(:,14) = fresli(:,14)+shg(:,i,3)*yl(:,i,2) 971 & +shg(:,i,1)*yl(:,i,4) 972 fresli(:,15) = fresli(:,15)+shg(:,i,3)*yl(:,i,3) 973 & +shg(:,i,2)*yl(:,i,4) 974 enddo 975 976 fresli(:,13) = pt5 * fresli(:,13) 977 fresli(:,14) = pt5 * fresli(:,14) 978 fresli(:,15) = pt5 * fresli(:,15) 979 980 strnrm(:,intp) = fresli(:,22) * sqrt( 981 & two * (fresli(:,10)**2 + fresli(:,11)**2 + fresli(:,12)**2) 982 & + four * ( fresli(:,13)**2 + fresli(:,14)**2 + 983 & fresli(:,15)**2 ) ) 984 985c 986c S_ij 987c 988 989 fresli(:,10) = fresli(:,10) * WdetJ ! u_{1,1}*WdetJ 990 fresli(:,11) = fresli(:,11) * WdetJ ! u_{2,2}*WdetJ 991 fresli(:,12) = fresli(:,12) * WdetJ ! u_{3,3}*WdetJ 992 fresli(:,13) = fresli(:,13) * WdetJ ! (1/2)*(u_{1,2}+u_{2,1})*WdetJ 993 fresli(:,14) = fresli(:,14) * WdetJ ! (1/2)*(u_{1,3}+u_{3,1})*WdetJ 994 fresli(:,15) = fresli(:,15) * WdetJ ! (1/2)*(u_{2,3}+u_{3,2})*WdetJ 995 996 fresli(:,22) = fresli(:,22) * WdetJ !rho * WdetJ 997c fresli(:,24) = fresli(:,24) * WdetJ 998 999 u1=zero 1000 u2=zero 1001 u3=zero 1002 do i=1,nshl 1003 u1 = u1 + shp(i,intp)*yl(:,i,2) 1004 u2 = u2 + shp(i,intp)*yl(:,i,3) 1005 u3 = u3 + shp(i,intp)*yl(:,i,4) 1006 enddo 1007 1008 fresli(:,1) = fresli(:,22) * u1 !rho u1 * WdetJ 1009 fresli(:,2) = fresli(:,22) * u2 !rho u2 * WdetJ 1010 fresli(:,3) = fresli(:,22) * u3 !rho u3 * WdetJ 1011 1012 fresli(:,4) = fresli(:,1) * u1 !rho u1 u1 *WdetJ 1013 fresli(:,5) = fresli(:,2) * u2 !rho u2 u2 *WdetJ 1014 fresli(:,6) = fresli(:,3) * u3 !rho u3 u3 *WdetJ 1015 fresli(:,7) = fresli(:,1) * u2 !rho u1 u2 *WdetJ 1016 fresli(:,8) = fresli(:,1) * u3 !rho u1 u3 *WdetJ 1017 fresli(:,9) = fresli(:,2) * u3 !rho u2 u3 *WdetJ 1018 1019 fresli(:,16) = strnrm(:,intp) * fresli(:,10) ! rho *|Eps| *Eps11 *WdetJ 1020 fresli(:,17) = strnrm(:,intp) * fresli(:,11) ! rho *|Eps| *Eps22 *WdetJ 1021 fresli(:,18) = strnrm(:,intp) * fresli(:,12) ! rho *|Eps| *Eps33 *WdetJ 1022 fresli(:,19) = strnrm(:,intp) * fresli(:,13) ! rho *|Eps| *Eps12 *WdetJ 1023 fresli(:,20) = strnrm(:,intp) * fresli(:,14) ! rho *|Eps| *Eps13 *WdetJ 1024 fresli(:,21) = strnrm(:,intp) * fresli(:,15) ! rho *|Eps| *Eps23 *WdetJ 1025 1026 fresli(:,23) = WdetJ ! Integral of 1 over the element 1027c 1028 do i = 1, 33 1029 fresl(:,i) = fresl(:,i) + fresli(:,i) 1030 enddo 1031 1032 enddo !end of loop over integration points 1033c 1034 do j = 1,nshl 1035 do nel = 1,npro 1036 fres(ien(nel,j),:) = fres(ien(nel,j),:) + fresl(nel,:) 1037 enddo 1038 enddo 1039 1040 return 1041 end 1042 1043 1044 1045 1046 1047 1048 1049 1050c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 1051c||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| 1052c||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| 1053 1054 1055 1056 1057 1058 1059 1060 subroutine disPtion (y, x, ien, shgl, shp, eps, xmudmi) 1061 1062 use rlssave ! Use the resolved Leonard stresses at the nodes. 1063 1064 include "common.h" 1065 1066 dimension xmudmi(npro,ngauss), y(nshg,ndof), 1067 & x(numnp,nsd), ien(npro,nshl), 1068 & shg(npro,nshl,nsd), 1069 & shgl(nsd,nshl,maxsh), shp(nshl,maxsh), 1070 & dxdxi(npro,nsd,nsd), dxidx(npro,nsd,nsd), 1071 & WdetJ(npro), 1072 & eps(3), fresli(npro,33), 1073 & epsli(npro,3), epsl(npro,3) 1074 1075 dimension yl(npro,nshl,nflow), xl(npro,nenl,nsd), 1076 & strnrm(npro,ngauss), rlsl(npro,nshl,6) 1077 1078 dimension tmp(npro) 1079 1080 call local (y, yl, ien, 5, 'gather ') 1081 call localx (x, xl, ien, 3, 'gather ') 1082 call local (rls, rlsl, ien, 6, 'gather ') 1083 1084 epsl = zero 1085 1086 yl (:,:,1) = one ! Unit density 1087 1088 do intp = 1, ngauss 1089 1090 fresli=zero 1091 1092c calculate the metrics 1093c 1094c 1095c.... ---------------------> Element Metrics <----------------------- 1096c 1097c.... compute the deformation gradient 1098c 1099 dxdxi = zero 1100c 1101 do n = 1, nenl 1102 dxdxi(:,1,1) = dxdxi(:,1,1) + xl(:,n,1) * shgl(1,n,intp) 1103 dxdxi(:,1,2) = dxdxi(:,1,2) + xl(:,n,1) * shgl(2,n,intp) 1104 dxdxi(:,1,3) = dxdxi(:,1,3) + xl(:,n,1) * shgl(3,n,intp) 1105 dxdxi(:,2,1) = dxdxi(:,2,1) + xl(:,n,2) * shgl(1,n,intp) 1106 dxdxi(:,2,2) = dxdxi(:,2,2) + xl(:,n,2) * shgl(2,n,intp) 1107 dxdxi(:,2,3) = dxdxi(:,2,3) + xl(:,n,2) * shgl(3,n,intp) 1108 dxdxi(:,3,1) = dxdxi(:,3,1) + xl(:,n,3) * shgl(1,n,intp) 1109 dxdxi(:,3,2) = dxdxi(:,3,2) + xl(:,n,3) * shgl(2,n,intp) 1110 dxdxi(:,3,3) = dxdxi(:,3,3) + xl(:,n,3) * shgl(3,n,intp) 1111 enddo 1112c 1113c.... compute the inverse of deformation gradient 1114c 1115 dxidx(:,1,1) = dxdxi(:,2,2) * dxdxi(:,3,3) 1116 & - dxdxi(:,3,2) * dxdxi(:,2,3) 1117 dxidx(:,1,2) = dxdxi(:,3,2) * dxdxi(:,1,3) 1118 & - dxdxi(:,1,2) * dxdxi(:,3,3) 1119 dxidx(:,1,3) = dxdxi(:,1,2) * dxdxi(:,2,3) 1120 & - dxdxi(:,1,3) * dxdxi(:,2,2) 1121 tmp = one / ( dxidx(:,1,1) * dxdxi(:,1,1) 1122 & + dxidx(:,1,2) * dxdxi(:,2,1) 1123 & + dxidx(:,1,3) * dxdxi(:,3,1) ) 1124 dxidx(:,1,1) = dxidx(:,1,1) * tmp 1125 dxidx(:,1,2) = dxidx(:,1,2) * tmp 1126 dxidx(:,1,3) = dxidx(:,1,3) * tmp 1127 dxidx(:,2,1) = (dxdxi(:,2,3) * dxdxi(:,3,1) 1128 & - dxdxi(:,2,1) * dxdxi(:,3,3)) * tmp 1129 dxidx(:,2,2) = (dxdxi(:,1,1) * dxdxi(:,3,3) 1130 & - dxdxi(:,3,1) * dxdxi(:,1,3)) * tmp 1131 dxidx(:,2,3) = (dxdxi(:,2,1) * dxdxi(:,1,3) 1132 & - dxdxi(:,1,1) * dxdxi(:,2,3)) * tmp 1133 dxidx(:,3,1) = (dxdxi(:,2,1) * dxdxi(:,3,2) 1134 & - dxdxi(:,2,2) * dxdxi(:,3,1)) * tmp 1135 dxidx(:,3,2) = (dxdxi(:,3,1) * dxdxi(:,1,2) 1136 & - dxdxi(:,1,1) * dxdxi(:,3,2)) * tmp 1137 dxidx(:,3,3) = (dxdxi(:,1,1) * dxdxi(:,2,2) 1138 & - dxdxi(:,1,2) * dxdxi(:,2,1)) * tmp 1139c 1140 wght=Qwt(lcsyst,intp) 1141 WdetJ = wght / tmp 1142 1143c 1144c 1145 do n = 1,nshl 1146 shg(:,n,1) = (shgl(1,n,intp) * dxidx(:,1,1) 1147 & + shgl(2,n,intp) * dxidx(:,2,1) 1148 & + shgl(3,n,intp) * dxidx(:,3,1)) 1149 shg(:,n,2) = (shgl(1,n,intp) * dxidx(:,1,2) 1150 & + shgl(2,n,intp) * dxidx(:,2,2) 1151 & + shgl(3,n,intp) * dxidx(:,3,2)) 1152 shg(:,n,3) = (shgl(1,n,intp) * dxidx(:,1,3) 1153 & + shgl(2,n,intp) * dxidx(:,2,3) 1154 & + shgl(3,n,intp) * dxidx(:,3,3)) 1155 enddo 1156 1157 1158c 1159c 1160 do i=1,nshl 1161 fresli(:,22) = fresli(:,22)+shp(i,intp) ! unit density at qpt 1162 enddo 1163 1164c 1165c 1166 do j=10,12 ! normal strainrate u_{i,i} no sum on i 1167 ig=j-9 1168 iv=j-8 1169 do i=1,nshl 1170 fresli(:,j) = fresli(:,j)+shg(:,i,ig)*yl(:,i,iv) 1171 enddo 1172 enddo 1173 1174c shear stresses NOTE there may be faster ways to do this 1175c check agains CM5 code for speed WTP 1176 1177 do i=1,nshl 1178 fresli(:,13) = fresli(:,13)+shg(:,i,2)*yl(:,i,2) 1179 & +shg(:,i,1)*yl(:,i,3) 1180 fresli(:,14) = fresli(:,14)+shg(:,i,3)*yl(:,i,2) 1181 & +shg(:,i,1)*yl(:,i,4) 1182 fresli(:,15) = fresli(:,15)+shg(:,i,3)*yl(:,i,3) 1183 & +shg(:,i,2)*yl(:,i,4) 1184 enddo 1185 1186 fresli(:,13) = pt5 * fresli(:,13) 1187 fresli(:,14) = pt5 * fresli(:,14) 1188 fresli(:,15) = pt5 * fresli(:,15) 1189 1190 strnrm(:,intp) = fresli(:,22) * ( 1191 & two * (fresli(:,10)**2 + fresli(:,11)**2 + fresli(:,12)**2) 1192 & + four * ( fresli(:,13)**2 + fresli(:,14)**2 + 1193 & fresli(:,15)**2 ) ) 1194 1195c 1196c 1197 1198 fresli(:,10) = fresli(:,10) * WdetJ ! u_{1,1}*WdetJ 1199 fresli(:,11) = fresli(:,11) * WdetJ ! u_{2,2}*WdetJ 1200 fresli(:,12) = fresli(:,12) * WdetJ ! u_{3,3}*WdetJ 1201 fresli(:,13) = fresli(:,13) * WdetJ ! (1/2)*(u_{1,2}+u_{2,1})*WdetJ 1202 fresli(:,14) = fresli(:,14) * WdetJ ! (1/2)*(u_{1,3}+u_{3,1})*WdetJ 1203 fresli(:,15) = fresli(:,15) * WdetJ ! (1/2)*(u_{2,3}+u_{3,2})*WdetJ 1204c 1205 strnrm(:,intp) = strnrm(:,intp) * WdetJ ! ( |Eps|^2 )*WdetJ 1206 1207 epsli(:,1) = -xmudmi(:,intp)*strnrm(:,intp) 1208 1209 epsli(:,2) = rlsl(:,intp,1)*fresli(:,10) + 1210 & rlsl(:,intp,2)*fresli(:,11) + 1211 & rlsl(:,intp,3)*fresli(:,12) + 1212 & two*( rlsl(:,intp,4)*fresli(:,13)+ 1213 & rlsl(:,intp,5)*fresli(:,14)+ 1214 & rlsl(:,intp,6)*fresli(:,15) ) 1215 1216 epsl(:,1) = epsl(:,1) + epsli(:,1) 1217 epsl(:,2) = epsl(:,2) + epsli(:,2) 1218 1219 enddo ! end loop over integ. pts 1220 1221 do i = 1, npro 1222 eps(1) = eps(1) + epsl(i,1) 1223 eps(2) = eps(2) + epsl(i,2) 1224 enddo 1225 1226 return 1227 end 1228 1229 1230 1231