29 SUBROUTINE cdleni(PM ,GEO ,STIFN ,STIFR,IXC,
30 . PX1 ,PX2 ,PY1 ,PY2 ,THK,
31 . IGEO,DT ,SH4TREE,ALDT,UPARAM ,
32 . IPM ,NLAY,PM_STACK,ISUBSTACK,STRC,
34 . X2L ,X3L ,X4L ,Y2L ,Y3L ,Y4L ,
35 . IGEO_STACK ,GROUP_PARAM)
40 use element_mod ,
only : nixc
44#include "implicit_f.inc"
53#include "remesh_c.inc"
54#include "vect01_c.inc"
59 INTEGER IXC(NIXC,*), IGEO(NPROPGI,*), SH4TREE(KSH4TREE,*),
60 . IPM(NPROPMI,*),NLAY,ISUBSTACK,IGEO_STACK(4*NPT_STACK+2,*)
64 .
area(mvsiz), strc(*),
65 . x2l(mvsiz),x3l(mvsiz),x4l(mvsiz),y2l(mvsiz),y3l(mvsiz),y4l(mvsiz)
66 TYPE (GROUP_PARAM_) :: GROUP_PARAM
70 INTEGER I,N, IMT, IPMAT, IGTYP,IPPID,IADB,
71 . IGMAT,IPGMAT,IPOS,,MLAWLY
73 . ssp(mvsiz), al1(mvsiz),al(mvsiz), almin(mvsiz),
74 . al2(mvsiz), al3(mvsiz), al4(mvsiz), al5(mvsiz), al6(mvsiz)
76 . viscmx,a11,a11r,b1,b2,vv,sti,stir,viscdef,dtdyn,rho,
77 . young,nu,gmax,fac,z0
78 my_real,
DIMENSION(MVSIZ) :: zoffset
80 igtyp = nint(geo(12,iprop))
81 igmat = igeo(98,iprop)
85 zoffset(lft:llt) = zero
95 zoffset(i) = z0 - half*thk(i)
97 ELSEIF (ipos== 3 .OR. ipos == 4)
THEN
104 zoffset(lft:llt) = zero
107 IF ((igtyp == 11 .AND. igmat < 0) .OR. igtyp == 16)
THEN
112 imt = igeo(ipmat+n,iprop)
113 ssp(i)=
max(ssp(i),pm(27,imt))
116 ELSEIF (mtn == 42)
THEN
119 imt = igeo(ipmat+n,iprop)
123 a11 = gmax*(one + nu)/(one - nu**2)
124 ssp(i)=
max(ssp(i), sqrt(a11/rho))
127 ELSEIF (mtn == 69)
THEN
130 imt = igeo(ipmat+n,iprop)
133 gmax = uparam(iadb+1)*uparam(iadb+6)
134 . + uparam(iadb+2)*uparam(iadb+7)
135 . + uparam(iadb+3)*uparam(iadb+8)
136 . + uparam(iadb+4)*uparam(iadb+9)
137 . + uparam(iadb+5)*uparam(iadb+10)
139 a11 = gmax*(one + nu)/(one - nu**2)
140 ssp(i)=
max(ssp(i), sqrt(a11/rho))
143 ELSEIF (mtn == 65)
THEN
146 imt = igeo(ipmat+n,iprop)
149 ssp(i)=
max(ssp(i), sqrt(young/rho))
155 imt = igeo(ipmat+n,iprop)
159 ssp(i)=
max(ssp(i), sqrt(young/(one-nu*nu)/rho))
163 ELSEIF(igtyp == 11 .AND. igmat > 0)
THEN
165 ssp(i) = geo(ipgmat +9 ,iprop)
167 ELSEIF(igtyp == 52 .OR.
168 . ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0))
THEN
170 ssp(i) = pm_stack(9 ,isubstack)
172 ELSEIF(igtyp == 17 .AND. igmat < 0)
THEN
179 imt = igeo_stack(ipmat + n,isubstack)
180 ssp(i)=
max(ssp(i),pm(27,imt))
183 ELSEIF (mtn == 42)
THEN
186 imt = igeo_stack(ipmat + n,isubstack)
190 a11 = gmax*(one + nu)/(one - nu**2)
191 ssp(i)=
max(ssp(i), sqrt(a11/rho))
194 ELSEIF (mtn == 69)
THEN
197 imt = igeo_stack(ipmat + n,isubstack)
200 gmax = uparam(iadb+1)*uparam(iadb+6)
201 . + uparam(iadb+2)*uparam(iadb+7)
202 . + uparam(iadb+3)*uparam(iadb+8)
203 . + uparam(iadb+4)*uparam(iadb+9)
204 . + uparam(iadb+5)*uparam(iadb+10)
206 a11 = gmax*(one + nu)/(one - nu**2)
207 ssp(i)=
max(ssp(i), sqrt(a11/rho))
210 ELSEIF (mtn == 69)
THEN
213 imt = igeo_stack(ipmat + n,isubstack)
216 gmax = uparam(iadb+1)*uparam(iadb+6)
217 . + uparam(iadb+2)*uparam(iadb+7)
218 . + uparam(iadb+3)*uparam
219 . + uparam(iadb+4)*uparam
220 . + uparam(iadb+5)*uparam(iadb+10)
222 a11 = gmax*(one + nu)/(one - nu**2)
226 ELSEIF (mtn == 65)
THEN
229 imt = igeo_stack(ipmat + n,isubstack)
232 ssp(i)=
max(ssp(i), sqrt(young/rho))
238 imt = igeo_stack(ipmat + n,isubstack)
242 ssp(i)=
max(ssp(i), sqrt(young/(one-nu*nu)/rho))
246 ELSEIF(igtyp == 51 .AND. igmat < 0)
THEN
251 imt = igeo_stack(ipmat + n,isubstack)
252 mlawly = nint(pm(19,imt))
253 IF (mlawly <= 28)
THEN
254 ssp(i)=
max(ssp(i),pm(27,imt))
255 ELSEIF (mlawly == 42)
THEN
259 a11 = gmax*(one + nu)/(one - nu**2)
260 ssp(i)=
max(ssp(i), sqrt(a11/rho))
261 ELSEIF (mlawly == 69)
THEN
264 gmax = uparam(iadb+1)*uparam(iadb+6)
265 . + uparam(iadb+2)*uparam(iadb+7)
266 . + uparam(iadb+3)*uparam(iadb+8)
267 . + uparam(iadb+4)*uparam(iadb+9)
268 . + uparam(iadb+5)*uparam(iadb+10)
270 a11 = gmax*(one + nu)/(one - nu**2)
271 ssp(i)=
max(ssp(i), sqrt(a11/rho))
272 ELSEIF (mlawly == 65)
THEN
275 ssp(i)=
max(ssp(i), sqrt(young/rho))
280 ssp(i)=
max(ssp(i), sqrt(young/(one-nu*nu)/rho))
289 ELSEIF (mtn == 42)
THEN
294 a11 = gmax*(one + nu)/(one - nu**2)
295 ssp(i)=
max(ssp(i), sqrt(a11/rho))
297 ELSEIF (mtn == 69)
THEN
301 gmax = uparam(iadb+1)*uparam(iadb+6)
302 . + uparam(iadb+2)*uparam(iadb+7)
303 . + uparam(iadb+3)*uparam(iadb+8)
304 . + uparam(iadb+4)*uparam(iadb+9)
305 . + uparam(iadb+5)*uparam(iadb+10)
307 a11 = gmax*(one + nu)/(one - nu**2)
308 ssp(i)=
max(ssp(i), sqrt(a11/rho))
310 ELSEIF (mtn == 65)
THEN
314 ssp(i)=sqrt(young/rho)
321 ssp(i)=sqrt(young/(one-nu*nu)/rho)
326 al1(i)= x2l(i) * x2l(i) + y2l(i) * y2l(i)
327 al2(i)=(x3l(i)-x2l(i))*(x3l(i)-x2l(i))+(y3l(i)-y2l(i))*(y3l(i)-y2l(i))
328 al3(i)=(x4l(i)-x3l(i))*(x4l(i)-x3l(i))+(y4l(i)-y3l(i))*(y4l(i)-y3l(i))
329 al4(i)= x4l(i) * x4l(i) + y4l(i) * y4l(i)
330 al5(i)=(x4l(i)-x2l(i))*(x4l(i)-x2l(i))+(y4l(i)-y2l(i))*(y4l(i)-y2l(i))
331 al6(i)= x3l(i) * x3l(i) + y3l(i) * y3l(i)
335 al(i)=
min(al1(i),al2(i),al3(i),al4(i),al5(i),al6(i))
336 IF(al3(i) == zero) al(i)=
min(al1(i),al2(i),al4(i))
342 ELSEIF(mtn == 25.OR.mtn == 27 .OR. mtn == 125 .OR. mtn == 127)
THEN
348 viscmx = group_param%VISC_DM
349 IF (viscmx == zero) viscmx = viscdef
350 IF (mtn == 1 .OR.mtn == 2.OR.mtn == 3.OR.
351 . mtn == 22.OR.mtn == 23) viscmx=zero
352 viscmx = sqrt(one + viscmx*viscmx) - viscmx
354 dtdyn =
area(i)/sqrt(
max(al5(i),al6(i)))
355 aldt(i) =
max(dtdyn,almin(i))
356 dt(i) = aldt(i)*viscmx/ssp(i)
363 IF (igtyp == 11 .AND. igmat > 0)
THEN
365 a11 =geo(ipgmat + 5,iprop)
366 a11r =geo(ipgmat + 7,iprop)
367 b1 = px1(i)*px1(i)+py1(i)*py1(i)
368 b2 = px2(i)*px2(i)+py2(i)*py2(i)
370 fac =
max(b1,b2) / (
area(i) * vv)
371 sti = fac * thk(i) * a11
372 stir = fac*a11r * thk(i)*(thk(i)**2 +
area(i))*one_over_12
373 stifn(ixc(2,i))=stifn(ixc(2,i))+sti
374 stifn(ixc(3,i))=stifn(ixc(3,i))+sti
375 stifn(ixc(4,i))=stifn(ixc(4,i))+sti
376 stifn(ixc(5,i))=stifn(ixc(5,i))+sti
377 stifr(ixc(2,i))=stifr(ixc(2,i))+stir
378 stifr(ixc(3,i))=stifr(ixc(3,i))+stir
379 stifr(ixc(4,i))=stifr(ixc(4,i))+stir
380 stifr(ixc(5,i))=stifr(ixc(5,i))+stir
383 ELSEIF(igtyp == 52 .OR.
384 . ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0 ))
THEN
386 a11 = pm_stack(5 ,isubstack)
387 a11r = pm_stack(7 ,isubstack)
388 b1 = px1(i)*px1(i)+py1(i)*py1(i)
389 b2 = px2(i)*px2(i)+py2(i)*py2(i)
391 fac =
max(b1,b2) / (
area(i) * vv)
392 sti = fac * thk(i) * a11
393 stir = fac*a11r * thk(i)*((thk(i)**2 +
area(i))*one_over_12 +
394 . zoffset(i)*zoffset(i) )
395 stifn(ixc(2,i))=stifn(ixc(2,i))+sti
396 stifn(ixc(3,i))=stifn(ixc(3,i))+sti
397 stifn(ixc(4,i))=stifn(ixc(4,i))+sti
398 stifn(ixc(5,i))=stifn(ixc(5,i))+sti
399 stifr(ixc(2,i))=stifr(ixc(2,i))+stir
400 stifr(ixc(3,i))=stifr(ixc(3,i))+stir
402 stifr(ixc(5,i))=stifr(ixc(5,i))+stir
408 b1 = px1(i)*px1(i)+py1(i)*py1(i)
409 b2 = px2(i)*px2(i)+py2(i)*py2(i)
412 . * thk(i) * a11 / (
area(i) * vv)
413 stir = sti * (thk(i)*thk(i) +
area(i)) / 12.
414 stifn(ixc(2,i))=stifn(ixc(2,i))+sti
415 stifn(ixc(3,i))=stifn(ixc(3,i))+sti
416 stifn(ixc(4,i))=stifn(ixc(4,i))+sti
417 stifn(ixc(5,i))=stifn(ixc(5,i))+sti
418 stifr(ixc(2,i))=stifr(ixc(2,i))+stir
419 stifr(ixc(3,i))=stifr(ixc(3,i))+stir
420 stifr(ixc(4,i))=stifr(ixc(4,i))+stir
421 stifr(ixc(5,i))=stifr(ixc(5,i))+stir
426 IF(igtyp == 11 .AND. igmat > 0)
THEN
429 IF(sh4tree(3,n) >= 0)
THEN
430 a11 =geo(ipgmat + 5,iprop)
431 a11r =geo(ipgmat + 7,iprop)
432 b1 = px1(i)*px1(i)+py1(i)*py1(i)
433 b2 = px2(i)*px2(i)+py2(i)*py2(i)
435 fac =
max(b1,b2) / (
area(i) * vv)
436 sti = fac * thk(i) * a11
437 stir = fac * a11r * thk(i)*(thk(i)**2 +
area(i))*one_over_12
438 stifn(ixc(2,i))=stifn(ixc(2,i))+sti
439 stifn(ixc(3,i))=stifn(ixc(3,i))+sti
440 stifn(ixc(4,i))=stifn(ixc(4,i))+sti
441 stifn(ixc(5,i))=stifn(ixc(5,i))+sti
442 stifr(ixc(2,i))=stifr(ixc(2,i))+stir
443 stifr(ixc(3,i))=stifr(ixc(3,i))+stir
444 stifr(ixc(4,i))=stifr(ixc(4,i))+stir
445 stifr(ixc(5,i))=stifr(ixc(5,i))+stir
449 ELSEIF(igtyp == 52 .OR.
450 . ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0 ))
THEN
453 IF(sh4tree(3,n) >= 0)
THEN
454 a11 = pm_stack(5 ,isubstack)
455 a11r = pm_stack(7 ,isubstack)
456 b1 = px1(i)*px1(i)+py1(i)*py1(i)
457 b2 = px2(i)*px2(i)+py2(i)*py2(i)
460 sti = fac * thk(i) * a11
461 stir = fac * a11r * thk(i)*((thk(i)**2 +
area(i))*one_over_12 +
462 . zoffset(i)*zoffset(i) )
463 stifn(ixc(2,i))=stifn(ixc(2,i))+sti
464 stifn(ixc(3,i))=stifn(ixc(3,i))+sti
465 stifn(ixc(4,i))=stifn(ixc(4,i))+sti
466 stifn(ixc(5,i))=stifn(ixc(5,i))+sti
467 stifr(ixc(2,i))=stifr(ixc(2,i))+stir
468 stifr(ixc(3,i))=stifr(ixc(3,i))+stir
469 stifr(ixc(4,i))=stifr(ixc(4,i))+stir
470 stifr(ixc(5,i))=stifr(ixc(5,i))+stir
478 IF(sh4tree(3,n) >= 0)
THEN
480 b1 = px1(i)*px1(i)+py1(i)*py1(i)
481 b2 = px2(i)*px2(i)+py2(i)*py2(i)
484 . * thk(i) * a11 / (
area(i) * vv)
485 stir = sti * (thk(i)*thk(i) +
area(i)) / 12.
487 stifn(ixc(3,i))=stifn(ixc(3,i))+sti
488 stifn(ixc(4,i))=stifn(ixc(4,i))+sti
489 stifn(ixc(5,i))=stifn(ixc(5,i))+sti
490 stifr(ixc(2,i))=stifr(ixc(2,i))+stir
491 stifr(ixc(3,i))=stifr(ixc(3,i))+stir
492 stifr(ixc(4,i))=stifr(ixc(4,i))+stir
493 stifr(ixc(5,i))=stifr(ixc(5,i))+stir
502 IF(geo(5,iprop)/=zero)geo(5,iprop)=
min(geo(5,iprop),dt(i))