33 1 NEL , NUPARAM, NUVAR , NFUNC , IFUNC , NPF ,
34 2 TF , TIME , TIMESTEP, UPARAM, RHO0 , RHO ,
35 3 VOLUME , EINT , NGL ,
36 4 EPSPXX , EPSPYY , EPSPZZ , EPSPXY, EPSPYZ, EPSPZX,
37 5 DEPSXX , DEPSYY , DEPSZZ , DEPSXY, DEPSYZ, DEPSZX,
38 6 EPSXX , EPSYY , EPSZZ , EPSXY , EPSYZ , EPSZX ,
39 7 SIGOXX , SIGOYY , SIGOZZ , SIGOXY, SIGOYZ, SIGOZX,
40 8 SIGNXX , SIGNYY , SIGNZZ , SIGNXY, SIGNYZ, SIGNZX,
41 9 SIGVXX , SIGVYY , SIGVZZ , SIGVXY, SIGVYZ, SIGVZX,
42 A SOUNDSP, VISCMAX, UVAR , OFF , ISMSTR, ISRATE ,
43 B ASRATE , ET ,IHET ,OFFG , EPSTH3, IEXPAN ,
48#include "implicit_f.inc"
57 INTEGER NEL, NUPARAM, NUVAR,ISMSTR,NGL(*),IHET,
60 . TIME , TIMESTEP , UPARAM(NUPARAM),
61 . RHO (NEL), RHO0 (), VOLUME(NEL), EINT(NEL),
62 . EPSPXX(NEL), EPSPYY(NEL), EPSPZZ(NEL),
63 . EPSPXY(NEL), EPSPYZ(NEL), EPSPZX(NEL),
64 . DEPSXX(NEL), DEPSYY(NEL), DEPSZZ(NEL),
65 . DEPSXY(NEL), DEPSYZ(NEL), DEPSZX(NEL),
66 . EPSXX (NEL), EPSYY (NEL), EPSZZ (NEL),
67 . EPSXY (NEL), EPSYZ (NEL), EPSZX (NEL),
68 . SIGOXX(NEL), SIGOYY(NEL), SIGOZZ(NEL),
69 . sigoxy(nel), sigoyz(nel), sigozx(nel),
70 . asrate,offg(nel),epsth3(nel),epsd(nel)
75 . signxx(nel), signyy(nel), signzz(nel),
76 . signxy(nel), signyz(nel), signzx(nel),
77 . sigvxx(nel), sigvyy(nel), sigvzz(nel),
78 . sigvxy(nel), sigvyz(nel), sigvzx(nel),
79 . soundsp(nel), viscmax(nel), et(nel)
84 . uvar(nel,nuvar), off(nel)
88 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
89 my_real FINTER,FINTTE,TF(*),FINT2V
90 EXTERNAL FINTER,FINTTE
94 INTEGER I,J,ITENS,NLOAD, NUNLOAD,NE_LOAD,NE_UNLOAD,IUNLOAD,
95 . jj,j1,j2,iunl_for,kk,nn,k,icase,idir,
96 . indx_l(nel),indx_unl(nel),itype,ifd(3),itag(nel)
98 . av(6,nel),ev(3,nel),ee(3,nel),edot(3,nel),
99 . evv(3,nel),dirprv(3,3,nel),erate(nel),rv(nel),
100 . yfac(nfunc),rate(nfunc),t1(nel),t2(nel),t3(nel),
101 . ee0(6,nel),ee1(6,nel),fg1(6),gmax(nel),fqstat(nel,3),
103 . deint,p,invr,df,fac,d,yfac_unl,f01,f02,f03,eps1,
104 . rbulk,nu,gs,hys,shape,f(3),f1,f2,dtinv,
105 . pvisc,dd,beta,dd1,dd2,dd3,aa,
106 . de1,de2,de3,de4,de5,de6, ert11,ert12,ert13,ert21,ert22,ert23,
107 . ert31,ert32,ert33,dep1,dep2,dep3,df1,df2,df3, xint,yint, lam_tens,
108 . t_tens, scale(3),es(3),evs(3),fstar,e2,f3,f_int,xinc,yinc,
109 . lam_comp,t_comp,e_tens,e_comp,yintl,bb,fc,ft,ec_start,et_start,
110 . ecurent(nel),epsp(3,nel),epe,dmax,sum,ft1,ft2,fc1,fc2,fcd,ftd,fcs,fts,
111 . epsdold(nel),yint0,yinc0,fts_unl,fcs_unl,dmoy,eintv,epseq(nel),
112 . t_tens_unl, t_comp_unl,scale_unl,scale_l,dam_f(nel),lam_xinc,
119 nload = int(uparam(4))
120 iunl_for = int(uparam(5))
123 icase = nint(uparam(9))
128 rate(j) = uparam(9 + 2*j - 1 )
129 yfac(j) = uparam(9 + 2*j )
131 yfac_unl = uparam(9 + 2*nload + 2 )
154 av(4,i)=half*epsxy(i)
155 av(5,i)=half*epsyz(i)
156 av(6,i)=half*epszx(i)
164 CALL valpvec(av,evv,dirprv,nel)
170 IF(ismstr==0.OR.ismstr==2.OR.ismstr==4)
THEN
173 ev(1,i)=exp(evv(1,i))
174 ev(2,i)=exp(evv(2,i))
175 ev(3,i)=exp(evv(3,i))
177 ELSEIF(ismstr==10.OR.ismstr==12)
THEN
179 IF(offg(i)<=one)
THEN
180 ev(1,i)=sqrt(evv(1,i)+ one)
181 ev(2,i)=sqrt(evv(2,i)+ one)
182 ev(3,i)=sqrt(evv(3,i)+ one)
184 ev(1,i)=evv(1,i)+ one
185 ev(2,i)=evv(2,i)+ one
186 ev(3,i)=evv(3,i)+ one
192 ev(1,i)=evv(1,i) + one
193 ev(2,i)=evv(2,i) + one
194 ev(3,i)=evv(3,i) + one
199 ee(1,i) = ev(1,i) - one
200 ee(2,i) = ev(2,i) - one
201 ee(3,i) = ev(3,i) - one
203 ee1(1,i) = dirprv(1,1,i)*dirprv(1,1,i)*ee(1,i)
205 . + dirprv(1,3,i)*dirprv(1,3,i)*ee(3,i)
207 ee1(2,i) = dirprv(2,2,i)*dirprv(2,2,i)*ee(2,i)
208 . + dirprv(2,3,i)*dirprv(2,3,i)*ee(3,i)
209 . + dirprv(2,1,i)*dirprv(2,1,i)*ee(1,i)
211 ee1(3,i) = dirprv(3,3,i)*dirprv(3,3,i)*ee(3,i)
212 . + dirprv(3,1,i)*dirprv(3,1,i)*ee(1,i)
213 . + dirprv(3,2,i)*dirprv(3,2,i)*ee(2,i)
215 ee1(4,i) = dirprv(1,1,i)*dirprv(2,1,i)*ee(1,i)
216 . + dirprv(1,2,i)*dirprv(2,2,i)*ee(2,i)
217 . + dirprv(1,3,i)*dirprv(2,3,i)*ee(3,i)
219 ee1(5,i) = dirprv(2,2,i)*dirprv(3,2,i)*ee(2,i)
220 . + dirprv(2,3,i)*dirprv(3,3,i)*ee(3,i)
221 . + dirprv(2,1,i)*dirprv(3,1,i)*ee(1,i)
223 ee1(6,i) = dirprv(3,3,i)*dirprv(1,3,i)*ee(3,i)
224 . + dirprv(3,1,i)*dirprv(1,1,i)*ee(1,i)
225 . + dirprv(3,2,i)*dirprv(1,2,i)*ee(2,i)
229 dtinv =timestep/
max(timestep**2,em20)
240 de1 = ee1(1,i) - ee0(1,i)
241 de2 = ee1(2,i) - ee0(2,i)
242 de3 = ee1(3,i) - ee0(3,i)
243 de4 = ee1(4,i) - ee0(4,i)
244 de5 = ee1(5,i) - ee0(5,i)
245 de6 = ee1(6,i) - ee0(6,i)
247 ert11 =dirprv(1,1,i)*de1 + dirprv(2,1,i)*de4 + dirprv(3,1,i)*de6
248 ert12 =dirprv(1,2,i)*de1 + dirprv(2,2,i)*de4 + dirprv(3,2,i)*de6
249 ert13 =dirprv(1,3,i)*de1 + dirprv(2,3,i)*de4 + dirprv(3,3,i)*de6
251 ert21 =dirprv(1,1,i)*de4 + dirprv(2,1,i)*de2 + dirprv(3,1,i)*de5
252 ert22 =dirprv(1,2,i)*de4 + dirprv(2,2,i)*de2 + dirprv(3,2,i)*de5
253 ert23 =dirprv(1,3,i)*de4 + dirprv(2,3,i)*de2 + dirprv(3,3,i)*de5
255 ert31 =dirprv(1,1,i)*de6 + dirprv(2,1,i)*de5 + dirprv(3,1,i)*de3
256 ert32 =dirprv(1,2,i)*de6 + dirprv(2,2,i)*de5 + dirprv(3,2,i)*de3
257 ert33 =dirprv(1,3,i)*de6 + dirprv(2,3,i)*de5 + dirprv(3,3,i)*de3
259 dep1 = dirprv(1,1,i)*ert11 + dirprv(2,1,i)*ert21
260 . + dirprv(3,1,i)*ert31
261 dep2 = dirprv(1,2,i)*ert12 + dirprv(2,2,i)*ert22
262 . + dirprv(3,2,i)*ert32
263 dep3 = dirprv(1,3,i)*ert13 + dirprv(2,3,i)*ert23
264 . + dirprv(3,3,i)*ert33
266 edot(1,i) = abs(dep1)*dtinv
267 edot(2,i) = abs(dep2)*dtinv
268 edot(3,i) = abs(dep3)*dtinv
270 erate(i) =
max(edot(1,i),edot(2,i),edot(3,i))
271 rv(i) = ev(1,i)*ev(2,i)*ev(3,i)
274 IF(iexpan > 0.AND.(ismstr==10.OR.ismstr==11.OR.ismstr==12))
THEN
276 rv(i) = rv(i) - epsth3(i)
282 edot(1,i) = asrate*edot(1,i) + (one - asrate)*uvar(i,7)
283 edot(2,i) = asrate*edot(2,i) + (one - asrate)*uvar(i,8)
284 edot(3,i) = asrate*edot(3,i) + (one - asrate)*uvar(i,9)
285 uvar(i,7) = edot(1,i)
286 uvar(i,8) = edot(2,i)
287 uvar(i,9) = edot(3,i)
289 epsd(i) =
max(edot(1,i),edot(2,i),edot(3,i))
297 uvar(i,7) = edot(1,i)
298 uvar(i,8) = edot(2,i)
299 uvar(i,9) = edot(3,i)
301 epsd(i) =
max(edot(1,i),edot(2,i),edot(3,i))
316 f(1) = ev(1,i)*yfac(1)*finter(ifunc(1),ee(1,i),npf,tf,df)
317 f(2) = ev(2,i)*yfac(1)*finter(ifunc(1),ee(2,i),npf,tf,df)
318 f(3) = ev(3,i)*yfac(1)*finter(ifunc(1),ee(3,i),npf,tf,df)
320 epseq(i) = sqrt(ee(1,i)**2 + ee(2,i)**2 + ee(3,i)**2)
322 invr = one/
max(em20,rv(i))
323 p = invr*rbulk*(rv(i) - one)
324 t1(i) =invr*(two_third*f(1) - third*(f(2) + f(3))) + p
325 t2(i) =invr*(two_third*f(2) - third*(f(1) + f(3))) + p
326 t3(i) =invr*(two_third*f(3) - third*(f(1) + f(2))) + p
327 ecurent(i)= half*(ee(1,i)*t1(i) + ee(2,i)*t2(i) + ee(3,i)*t3(i))
328 ecurent(i) =
max(em20,ecurent(i))
329 IF(ecurent(i) >= uvar(i,17)) uvar(i,20) = zero
330 uvar(i,17)=
max(uvar(i,17), ecurent(i))
331 deint = ecurent(i) - uvar(i,16)
332 uvar(i,16) = ecurent(i)
337 IF(uvar(i,19) == -one)
THEN
338 IF(deint/uvar(i,16) >= 1e-7)
THEN
339 ne_load = ne_load + 1
344 IF(iunl_for > 1)
THEN
353 ne_unload = ne_unload + 1
354 indx_unl(ne_unload) = i
358 IF(deint/
max(em20,uvar(i,17)) >= zero)
THEN
359 ne_load = ne_load + 1
363 ne_unload = ne_unload + 1
364 indx_unl(ne_unload) = i
368 lam_tens =
max(one,ev(1,i),ev(2,i),ev(3,i))
369 e_tens =
max(zero,ee(1,i),ee(2,i),ee(3,i))
370 t_tens=
max(zero,f(1),f(2),f(3))
371 lam_comp =
min(one,ev(1,i),ev(2,i),ev(3,i))
372 e_comp =
min(zero,ee(1,i),ee(2,i),ee(3,i))
373 t_comp =
min(zero,f(1),f(2),f(3))
374 IF(iunl_for == 1)
THEN
377 IF(uvar(i,20) == zero )
THEN
381 IF( epsdold(i) >= rate(k) )j1 = k
384 aa = (epsdold(i) - rate(j1))/(rate
387 ft1 =lam_tens*yfac(j1)*finter(ifunc(j1),e_tens,npf,tf,df1)
388 ft2 =lam_tens*yfac(j2)*finter(ifunc(j2),e_tens,npf,tf,df2)
389 fts =lam_tens*yfac_unl*finter(ifunc(kk),e_tens,npf,tf,df)
391 fc1 =lam_comp*yfac(j1)*finter(ifunc(j1),e_comp,npf,tf,df1)
392 fc2 =lam_comp*yfac(j2)*finter(ifunc(j2),e_comp,npf,tf,df2)
393 fcs =lam_comp*yfac_unl*finter(ifunc(kk),e_comp,npf,tf,df1)
394 ftd = ft1 + aa*(ft2 - ft1)
395 fcd = fc1 + aa*(fc2 - fc1)
398 uvar(i,22)= uvar(i,22)/e_tens
399 uvar(i,25)= uvar(i,25)/abs(e_comp)
403 IF( epsdold(i) >= rate(k) )j1 = k
406 aa = (epsdold(i) - rate(j1))/(rate(j2) - rate(j1))
409 ft1 =lam_tens*yfac(j1)*finter(ifunc(j1),e_tens,npf,tf,df1)
410 ft2 =lam_tens*yfac(j2)*finter(ifunc(j2),e_tens,npf,tf,df2)
411 fts =lam_tens*yfac(1)*finter(ifunc(1),e_tens,npf,tf,df)
413 fc1 =lam_comp*yfac(j1)*finter(ifunc(j1),e_comp,npf,tf
414 fc2 =lam_comp*yfac(j2)*finter(ifunc(j2),e_comp,npf,tf,df2)
415 fcs =lam_comp*yfac(1)*finter(ifunc(1),e_comp,npf,tf,df1)
419 dd1 = lam_tens**fac - one
420 ft1 =ft1 + (dd1 + one)*yfac(j1)*finter(ifunc(j1),dd1,npf,tf,df)
421 ft2 =ft2 + (dd1 + one)*yfac(j2)*finter(ifunc(j2),dd1,npf,tf,df)
422 fts =fts + (dd1 + one)*yfac(1)*finter(ifunc(1),dd1,npf,tf,df)
423 dd1 = lam_comp**fac - one
424 fc1 =fc1 + (dd1 + one)*yfac(j1)*finter(ifunc(j1),dd1,npf,tf,df)
425 fc2 =fc2 + (dd1 + one)*yfac(j2)*finter(ifunc(j2),dd1,npf,tf,df)
426 fcs =fcs + (dd1 + one)*yfac(1)*finter(ifunc(1),dd1,npf,tf,df)
428 ftd = ft1 + aa*(ft2 - ft1)
429 fcd = fc1 + aa*(fc2 - fc1)
430 uvar(i,22) = (ftd/fts - one)/e_tens
431 uvar(i,25) = (fcd/fcs - one)/e_comp
440 IF(ne_load > 0 )
THEN
442 IF(iunl_for == 1)
THEN
446 invr = one/
max(em20,rv(i))
447 dam(i) = uvar(i,16)/uvar(i,17)
449 f(1) = ev(1,i)*yfac(1)*finter(ifunc(1),ee(1,i),npf,tf,df1)
450 f(2) = ev(2,i)*yfac(1)*finter(ifunc(1),ee
451 f(3) = ev(3,i)*yfac(1)*finter(ifunc(1),ee(3,i),npf,tf,df3)
452 gmax(i) =
max(gmax(i),yfac(1)*df1,yfac(1)*df2,
454 IF(dam(i) /= one)
THEN
455 funl(1)=ev(1,i)*yfac_unl*finter(ifunc(kk),ee(1,i),npf,tf,df)
456 funl(2)=ev(2,i)*yfac_unl*finter(ifunc(kk),ee(2,i),npf,tf,df)
457 funl(3)=ev(3,i)*yfac_unl*finter(ifunc(kk),ee(3,i),npf,tf,df)
459 IF(f(1) /= zero) dam_f(i) =
max(dam_f(i),funl(1)/f(1))
460 IF(f(2) /= zero) dam_f(i) =
max(dam_f(i),funl(2)/f(2))
461 IF(f(3) /= zero) dam_f(i) =
max(dam_f(i),funl(3)/f(3))
462 IF(dam_f(i) <= one)
THEN
463 dam(i) =
max(dam(i),dam_f(i))
468 dam(i) =
max(dam(i),uvar(i,18))
469 dam(i) =
min(one, dam(i))
473 dd1 = ev(1,i)**fac - one
474 dd2 = ev(2,i)**fac - one
475 dd3 = ev(3,i)**fac - one
477 . (dd1 + one)*yfac(1)*finter(ifunc(1),dd1,npf,tf,df)
479 . (dd2 + one)*yfac(1)*finter(ifunc(1),dd2,npf,tf,df)
481 . (dd3 + one)*yfac(1)*finter(ifunc(1),dd3,npf,tf,df)
484 p = invr*rbulk*(rv(i) - one)
485 t1(i) = dam(i)*(invr*(two_third*f(1) - third*(f(2) + f(3))) )+ p
486 t2(i) = dam(i)*(invr*(two_third*f(2) - third*(f(1) + f(3))) )+ p
487 t3(i) = dam(i)*(invr*(two_third*f(3) - third*(f(1) + f(2))) )+ p
493 invr = one/
max(em20,rv(i))
494 f(1) = ev(1,i)*yfac(1)*finter(ifunc(1),ee(1,i),npf,tf,df1)
495 f(2) = ev(2,i)*yfac(1)*finter(ifunc(1),ee(2,i),npf,tf,df2)
496 f(3) = ev(3,i)*yfac(1)*finter(ifunc(1),ee(3,i),npf,tf,df3)
497 gmax(i) =
max(gmax(i),yfac(1)*df1,yfac(1)*df2,
502 dd1 = ev(1,i)**fac - one
503 dd2 = ev(2,i)**fac - one
504 dd3 = ev(3,i)**fac - one
506 . (dd1 + one)*yfac(1)*finter(ifunc(1),dd1,npf,tf,df)
508 . (dd2 + one)*yfac(1)*finter(ifunc(1),dd2,npf,tf,df)
510 . (dd3 + one)*yfac(1)*finter(ifunc(1),dd3,npf,tf,df)
512 p = invr*rbulk*(rv(i) - one)
513 t1(i) =uvar(i,22)*invr*(two_third*f(1) - third*(f(2) + f(3))) + p
514 t2(i) =uvar(i,22)*invr*(two_third*f(2) - third*(f(1) + f(3))) + p
515 t3(i) =uvar(i,22)*invr*(two_third*f(3) - third*(f(1) + f(2))) + p
520 IF(iunl_for == 1)
THEN
524 IF(uvar(i,20) == one .AND. ecurent(i) <= uvar(i,17))
THEN
525 f(1)=ev(1,i)*yfac_unl*finter(ifunc(kk),ee(1,i),npf,tf,df)
526 f(2)=ev(2,i)*yfac_unl*finter(ifunc(kk),ee(2,i),npf,tf,df)
527 f(3)=ev(3,i)*yfac_unl*finter(ifunc(kk),ee(3,i),npf,tf,df)
528 e_tens =
max(zero,ee(1,i),ee(2,i),ee(3,i))
529 dam(i)= uvar(i,22)*e_tens
530 IF(itype == -1 .OR. (itype == 2 .AND. rv(i) < one))
THEN
531 e_comp =
min(zero,ee(1,i),ee(2,i),ee(3,i))
532 dam(i)= uvar(i,25)*abs(e_comp)
537 dd1 = ev(1,i)**fac - one
538 dd2 = ev(2,i)**fac - one
539 dd3 = ev(3,i)**fac - one
541 . (dd1 + one)*yfac_unl*finter(ifunc(kk),dd1,npf,tf,df)
543 . (dd2 + one)*yfac_unl*finter(ifunc(kk),dd2,npf,tf,df)
545 . (dd3 + one)*yfac_unl*finter(ifunc(kk),dd3,npf,tf,df)
547 invr = one/
max(em20,rv(i))
548 p = invr*rbulk*(rv(i) - one)
549 t1(i) = invr*dam(i)*(two_third*f(1) - third*(f(2) + f(3))) + p
550 t2(i) = invr*dam(i)*(two_third*f(2) - third*(f(1) + f(3))) + p
551 t3(i) = invr*dam(i)*(two_third*f(3) - third*(f(1) + f(2))) + p
556 IF( edot(j,i) >= rate(k) )j1 = k
559 aa = (edot(j,i) - rate(j1))/(rate(j2) - rate(j1))
562 f1 =ev(j,i)*yfac(j1)*finter(ifunc(j1),ee(j,i),npf,tf,df1)
563 f2 =ev(j,i)*yfac(j2)*finter(ifunc(j2),ee(j,i),npf,tf,df2)
566 df =
max(df, df1 + aa*(df2 - df1))
567 gmax(i) =
max(gmax(i),yfac(j1)*df1,yfac(j2)*df2 )
571 dd1 = ev(j,i)**fac - one
572 f1 =f1+(dd1 + one)*yfac(j1)*finter(ifunc(j1),dd1,npf,tf,df)
573 f2 =f2+(dd1 + one)*yfac(j2)*finter(ifunc(j2),dd1,npf,tf,df)
575 f(j) = f1 + aa*(f2 - f1)
577 invr = one/
max(em20,rv(i))
578 p = invr*rbulk*(rv(i) - one)
579 t1(i) =invr*(two_third*f(1) - third*(f(2) + f(3))) + p
580 t2(i) =invr*(two_third*f(2) - third*(f(1) + f(3))) + p
581 t3(i) =invr*(two_third*f(3) - third*(f(1) + f(2))) + p
590 IF( edot(j,i) >= rate(k) )j1 = k
593 aa = (edot(j,i) - rate(j1))/(rate(j2) - rate(j1))
596 f1 =ev(j,i)*yfac(j1)*finter(ifunc(j1),ee(j,i),npf,tf,df1)
597 f2 =ev(j,i)*yfac(j2)*finter(ifunc(j2),ee(j,i),npf,tf,df2)
600 df =
max(df, df1 + aa*(df2 - df1))
601 gmax(i) =
max(gmax(i),yfac(j1)*df1,yfac
605 dd1 = ev(j,i)**fac - one
606 f1 =f1+(dd1 + one)*yfac
607 f2 =f2+(dd1 + one)*yfac(j2)*finter(ifunc(j2),dd1,npf,tf,df)
609 f(j) = f1 + aa*(f2 - f1)
611 invr = one/
max(em20,rv(i))
612 p = invr*rbulk*(rv(i) - one)
613 t1(i) =invr*(two_third*f(1) - third*(f(2) + f(3))) + p
614 t2(i) =invr*(two_third*f(2) - third*(f(1) + f(3))) + p
615 t3(i) =invr*(two_third*f(3) - third*(f(1) + f(2))) + p
623 IF(ne_unload > 0 )
THEN
625 SELECT CASE (iunl_for)
630 f(1) = ev(1,i)*yfac(1)*finter(ifunc(1),ee(1,i),npf,tf,df1)
631 f(2) = ev(2,i)*yfac(1)*finter(ifunc(1),ee(2,i),npf,tf,df2)
632 f(3) = ev(3,i)*yfac(1)*finter(ifunc(1),ee(3,i),npf,tf,df3)
633 gmax(i) =
max(gmax(i),yfac(1)*df1,yfac
638 dd1 = ev(1,i)**fac - one
639 dd2 = ev(2,i)**fac - one
640 dd3 = ev(3,i)**fac - one
642 . (dd1 + one)*yfac(1)*finter(ifunc
644 . (dd2 + one)*yfac(1)*finter(ifunc(1),dd2,npf
646 . (dd3 + one)*yfac(1)*finter(ifunc(1)
648 invr = one/
max(em20,rv(i))
649 p = invr*rbulk*(rv(i) - one)
650 t1(i) = invr*(two_third*f(1) - third*(f(2) + f(3))) + p
651 t2(i) = invr*(two_third*f(2) - third*(f(1) + f(3))) + p
652 t3(i) = invr*(two_third*f(3) - third*(f(1) + f(2))) + p
662 dam(i) = uvar(i,16)/uvar(i,17)
663 f(1) = ev(1,i)*yfac(1)*finter(ifunc(1),ee(1,i),npf,tf,df1)
664 f(2) = ev(2,i)*yfac(1)*finter(ifunc(1),ee(2,i),npf,tf,df2)
665 f(3) = ev(3,i)*yfac(1)*finter(ifunc(1),ee(3,i),npf,tf,df3)
666 gmax(i) =
max(gmax(i),yfac(1)*df1,yfac(1)*df2,
668 funl(1)=ev(1,i)*yfac_unl*finter(ifunc(kk),ee(1,i),npf,tf,df)
669 funl(2)=ev(2,i)*yfac_unl*finter(ifunc(kk),ee(2,i),npf,tf,df)
670 funl(3)=ev(3,i)*yfac_unl*finter(ifunc(kk),ee(3,i),npf,tf,df)
672 IF(f(1)/=zero) dam_f(i) =
max(dam_f(i),funl(1)/f(1))
673 IF(f(2)/=zero) dam_f(i) =
max(dam_f(i),funl(2)/f(2))
674 IF(f(3)/=zero) dam_f(i) =
max(dam_f(i),funl(3)/f(3))
675 IF(dam_f(i) <= one)
THEN
676 dam(i) =
max(dam(i),dam_f(i))
680 dam(i) =
min(one,dam(i))
685 dd1 = ev(1,i)**fac - one
686 dd2 = ev(2,i)**fac - one
687 dd3 = ev(3,i)**fac - one
689 . (dd1 + one)*yfac(1)*finter(ifunc(
691 . (dd2 + one)*yfac(1)*finter(ifunc(1),dd2,npf,tf,df)
695 invr = one/
max(em20,rv(i))
696 p = invr*rbulk*(rv(i) - one)
697 t1(i) = dam(i)*(invr*(two_third*f(1) - third*(f(2) + f(3))) )+ p
698 t2(i) = dam(i)*(invr*(two_third*f(2) - third*(f(1) + f(3))) )+ p
699 t3(i) = dam(i)*(invr*(two_third*f(3) - third*(f(1) + f(2))) )+ p
706 funl(1)=ev(1,i)*yfac_unl*finter(ifunc(kk),ee(1,i),npf,tf,df)
707 funl(2)=ev(2,i)*yfac_unl*finter(ifunc(kk),ee(2,i),npf,tf,df)
708 funl(3)=ev(3,i)*yfac_unl*finter(ifunc(kk),ee(3,i),npf,tf,df)
710 f(1)=ev(1,i)*yfac(1)*finter(ifunc(1),ee(1,i),npf,tf,df)
711 f(2)=ev(2,i)*yfac(1)*finter(ifunc(1),ee(2,i),npf,tf,df)
712 f(3)=ev(3,i)*yfac(1)*finter(ifunc(1),ee(3,i),npf,tf,df)
714 IF(f(1)/=zero) dam(i) =
715 IF(f(2)/=zero) dam(i) =
max(dam(i),funl(2)/f(2))
716 IF(f(3)/=zero) dam(i) =
max(dam(i),funl(3)/f(3))
717 IF(itype == -1 .OR. (itype == 2 .AND. rv(i) < one))
THEN
718 e_comp =
min(zero,ee(1,i),ee(2,i),ee(3,i))
719 dam(i)= uvar(i,25)*dam(i)*abs(e_comp)
721 e_tens =
max(zero,ee(1,i),ee(2,i),ee(3,i))
722 dam(i)= uvar(i,22)*dam(i)*e_tens
728 dd1 = ev(1,i)**fac - one
729 dd2 = ev(2,i)**fac - one
730 dd3 = ev(3,i)**fac - one
732 . (dd1 + one)*yfac(1)*finter(ifunc(1),dd1,npf,tf,df)
734 . (dd2 + one)*yfac(1)*finter(ifunc(1),dd2,npf,tf,df)
736 . (dd3 + one)*yfac(1)*finter(ifunc(1),dd3,npf,tf,df)
738 invr = one/
max(em20,rv(i))
739 p = invr*rbulk*(rv(i) - one)
740 t1(i) = invr*dam(i)*(two_third*f(1) - third*(f(2) + f(3))) + p
741 t2(i) = invr*dam(i)*(two_third*f(2) - third*(f(1) + f(3))) + p
742 t3(i) = invr*dam(i)*(two_third*f(3) - third*(f(1) + f(2))) + p
753 dam(i) = one - (uvar(i,16)/
max(em20,uvar(i,17)))**shape
754 dam(i) = one - (one - hys)*dam(i)
760 f(1)=ev(1,i)*yfac(1)*finter(ifunc(1),ee(1,i),npf,tf,df)
761 f(2)=ev(2,i)*yfac(1)*finter(ifunc(1),ee(2,i),npf,tf,df)
762 f(3)=ev(3,i)*yfac(1)*finter(ifunc(1),ee(3,i),npf,tf,df)
766 dd1 = ev(1,i)**fac - one
767 dd2 = ev(2,i)**fac - one
768 dd3 = ev(3,i)**fac - one
770 . (dd1 + one)*yfac(1)*finter(ifunc(1),dd1,npf,tf,df)
772 . (dd2 + one)*yfac(1)*finter(ifunc(1),dd2,npf,tf,df)
774 . (dd3 + one)*yfac(1)*finter(ifunc(1),dd3,npf,tf,df)
776 invr = one/
max(em20,rv(i))
777 p = invr*rbulk*(rv(i) - one)
778 t1(i) = invr*dam(i)*(two_third*f(1) - third*(f(2) + f(3))) + p
779 t2(i) = invr*dam(i)*(two_third*f(2) - third*(f(1) + f(3))) + p
780 t3(i) = invr*dam(i)*(two_third*f(3) - third*(f(1) + f(2))) + p
790 dam(i) = one - (uvar(i,16)/
max(em20,uvar(i,17)))**shape
791 dam(i) = one - (one - hys)*dam(i)
798 scale(1) = (uvar(i,22)* ee(1,i) + one)*yfac(1)
800 scale(1) = (uvar(i,25)* ee(1,i) + one)*yfac(1)
802 IF(ee(2,i)>= zero)
THEN
803 scale(2) = (uvar(i,22)* ee(2,i) + one)*yfac(1)
805 scale(2) = (uvar(i,25)* ee(2,i) + one)*yfac(1)
807 IF(ee(3,i) >= zero)
THEN
808 scale(3) = (uvar(i,22)* ee(3,i) + one)*yfac(1)
810 scale(3) = (uvar(i,25)* ee(3,i) + one)*yfac(1)
813 f(1) = ev(1,i)*scale(1)*finter(ifunc(1),ee(1,i),npf,tf,df1)
814 f(2) = ev(2,i)*scale(2)*finter(ifunc(1),ee(2,i),npf,tf,df2)
815 f(3) = ev(3,i)*scale(3)*finter(ifunc(1),ee(3,i),npf,tf,df3)
817 gmax(i) =
max(gmax(i),scale(1)*df1,scale(2)*df2,
822 dd1 = ev(1,i)**fac - one
823 dd2 = ev(2,i)**fac - one
824 dd3 = ev(3,i)**fac - one
826 . (dd1 + one)*scale(1)*finter(ifunc(1),dd1,npf,tf,df)
828 . (dd2 + one)*scale(2)*finter(ifunc(1),dd2,npf,tf,df)
830 . (dd3 + one)*scale(3)*finter(ifunc(1),dd3,npf,tf,df)
832 invr = one/
max(em20,rv(i))
833 p = invr*rbulk*(rv(i) - one)
834 t1(i) = invr*dam(i)*(two_third*f(1) - third*(f(2) + f(3))) + p
835 t2(i) = invr*dam(i)*(two_third*f(2) - third*(f(1) + f(3))) + p
836 t3(i) = invr*dam(i)*(two_third*f(3) - third*(f(1) + f(2))) + p
845 signxx(i) = dirprv(1,1,i)*dirprv(1,1,i)*t1(i)
846 . + dirprv(1,2,i)*dirprv(1,2,i)*t2(i)
847 . + dirprv(1,3,i)*dirprv(1,3,i)*t3(i)
849 signyy(i) = dirprv(2,2,i)*dirprv(2,2,i)*t2
850 . + dirprv(2,3,i)*dirprv(2,3,i)*t3(i)
851 . + dirprv(2,1,i)*dirprv(2,1,i)*t1(i)
853 signzz(i) = dirprv(3,3,i)*dirprv(3,3,i)*t3(i)
854 . + dirprv(3,1,i)*dirprv(3,1,i)*t1(i)
855 . + dirprv(3,2,i)*dirprv(3,2,i)*t2(i)
857 signxy(i) = dirprv(1,1,i)*dirprv(2,1,i)*t1(i)
858 . + dirprv(1,2,i)*dirprv(2,2,i)*t2(i)
859 . + dirprv(1,3,i)*dirprv(2,3,i)*t3(i)
861 signyz(i) = dirprv(2,2,i)*dirprv(3,2,i)*t2(i)
862 . + dirprv(2,3,i)*dirprv(3,3,i)*t3(i)
863 . + dirprv(2,1,i)*dirprv(3,1,i)*t1(i)
865 signzx(i) = dirprv(3,3,i)*dirprv(1,3,i)*t3(i)
866 . + dirprv(3,1,i)*dirprv(1,1,i)*t1(i)
867 . + dirprv(3,2,i)*dirprv(1,2,i)*t2(i)
870 soundsp(i)=sqrt((four_over_3*gs + rbulk)/rho(i))
874 IF (impl_s > 0 .OR. ihet > 1)
THEN
876 et(i)=
max(one,gmax(i)/gs)