36
37
38
39 USE group_param_mod
40
41
42
43#include "implicit_f.inc"
44
45
46
47#include "mvsiz_p.inc"
48
49
50
51#include "param_c.inc"
52#include "remesh_c.inc"
53#include "vect01_c.inc"
54
55
56
57 INTEGER IMAT,IPROP
58 INTEGER IXC(NIXC,*), IGEO(NPROPGI,*), SH4TREE(KSH4TREE,*),
59 . IPM(NPROPMI,*),NLAY,ISUBSTACK,IGEO_STACK(4*NPT_STACK+2,*)
61 . pm(npropm,*), geo(npropg,*),stifn(*),stifr(*),uparam(*),
62 . px1(*),px2(*),py1(*),py2(*),thk(*),dt(*),aldt(*),pm_stack(20,*),
63 .
area(mvsiz), strc(*),
64 . x2l(mvsiz),x3l(mvsiz),x4l(mvsiz),y2l(mvsiz),y3l(mvsiz),y4l(mvsiz)
65 TYPE (GROUP_PARAM_) :: GROUP_PARAM
66
67
68
69 INTEGER I,N, IMT, IPMAT, IGTYP,IPPID,IADB,
70 . I1,I3,IPTHK,IPPOS,I2,MATLY,IGMAT,IPGMAT,IPOS,NIP,MLAWLY
72 . ssp(mvsiz), al1(mvsiz),al(mvsiz), almin(mvsiz),
73 . al2(mvsiz), al3(mvsiz), al4(mvsiz), al5(mvsiz), al6(mvsiz)
75 . viscmx,a11,a11r,a12,b1,b2,c1,vv,sti,stir,viscdef,dtdyn,rho,
76 . young,nu,gmax,thkly,posly,fac,z0
77 my_real,
DIMENSION(MVSIZ) :: zoffset
78
79 igtyp = nint(geo(12,iprop))
80 igmat = igeo(98,iprop)
81 ipgmat = 700
82 ssp(lft:llt) = zero
83 z0 = geo(199,iprop)
84 zoffset(lft:llt) = zero
85 SELECT CASE(igtyp)
86 CASE (1,9,10,11,16)
87 DO i=lft,llt
88 zoffset(i) = z0
89 ENDDO
90 CASE (17,51,52)
91 ipos = igeo(99,iprop)
92 IF(ipos == 2) THEN
93 DO i=lft,llt
94 zoffset(i) = z0 - half*thk(i)
95 ENDDO
96 ELSEIF (ipos== 3 .OR. ipos == 4) THEN
97 DO i=lft,llt
98 z0= half*thk(i)
99 zoffset(i) = z0
100 ENDDO
101 ENDIF
102 CASE DEFAULT
103 zoffset(lft:llt) = zero
104 END SELECT
105
106 IF ((igtyp == 11 .AND. igmat < 0) .OR. igtyp == 16) THEN
107 ipmat = 100
108 IF (mtn <= 28) THEN
109 DO i=lft,llt
110 DO n=1,npt
111 imt = igeo(ipmat+n,iprop)
112 ssp(i)=
max(ssp(i),pm(27,imt))
113 ENDDO
114 ENDDO
115 ELSEIF (mtn == 42) THEN
116 DO i=lft,llt
117 DO n=1,npt
118 imt = igeo(ipmat+n,iprop)
119 rho = pm(1,imt)
120 nu = pm(21,imt)
121 gmax = pm(22,imt)
122 a11 = gmax*(one + nu)/(one - nu**2)
123 ssp(i)=
max(ssp(i), sqrt(a11/rho))
124 ENDDO
125 ENDDO
126 ELSEIF (mtn == 69) THEN
127 DO i=lft,llt
128 DO n=1,npt
129 imt = igeo(ipmat+n,iprop)
130 iadb = ipm(7,imt)-1
131 nu = uparam(iadb+14)
132 gmax = uparam(iadb+1)*uparam(iadb+6)
133 . + uparam(iadb+2)*uparam(iadb+7)
134 . + uparam(iadb+3)*uparam(iadb+8)
135 . + uparam(iadb+4)*uparam(iadb+9)
136 . + uparam(iadb+5)*uparam(iadb+10)
137 rho = pm(1,imt)
138 a11 = gmax*(one + nu)/(one - nu**2)
139 ssp(i)=
max(ssp(i), sqrt(a11/rho))
140 ENDDO
141 ENDDO
142 ELSEIF (mtn == 65) THEN
143 DO i=lft,llt
144 DO n=1,npt
145 imt = igeo(ipmat+n,iprop)
146 rho =pm(1,imt)
147 young=pm(20,imt)
148 ssp(i)=
max(ssp(i), sqrt(young/rho))
149 ENDDO
150 ENDDO
151 ELSE
152 DO i=lft,llt
153 DO n=1,npt
154 imt = igeo(ipmat+n,iprop)
155 rho =pm(1,imt)
156 young=pm(20,imt)
157 nu =pm(21,imt)
158 ssp(i)=
max(ssp(i), sqrt(young/(one-nu*nu)/rho))
159 ENDDO
160 ENDDO
161 ENDIF
162 ELSEIF(igtyp == 11 .AND. igmat > 0) THEN
163 DO i=lft,llt
164 ssp(i) = geo(ipgmat +9 ,iprop)
165 ENDDO
166 ELSEIF(igtyp == 52 .OR.
167 . ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0)) THEN
168 DO i=lft,llt
169 ssp(i) = pm_stack(9 ,isubstack)
170 ENDDO
171 ELSEIF(igtyp == 17 .AND. igmat < 0) THEN
172 ippid = 100
173 nip = npt
174 ipmat = 2 + nip
175 IF(mtn<=28)THEN
176 DO i=lft,llt
177 DO n=1,nip
178 imt = igeo_stack(ipmat + n,isubstack)
179 ssp(i)=
max(ssp(i),pm(27,imt))
180 ENDDO
181 ENDDO
182 ELSEIF (mtn == 42) THEN
183 DO i=lft,llt
184 DO n=1,nip
185 imt = igeo_stack(ipmat + n,isubstack)
186 rho = pm(1,imt)
187 nu = pm(21,imt)
188 gmax = pm(22,imt)
189 a11 = gmax*(one + nu)/(one - nu**2)
190 ssp(i)=
max(ssp(i), sqrt(a11/rho))
191 ENDDO
192 ENDDO
193 ELSEIF (mtn == 69) THEN
194 DO i=lft,llt
195 DO n=1,nip
196 imt = igeo_stack(ipmat + n,isubstack)
197 iadb = ipm(7,imt)-1
198 nu = uparam(iadb+14)
199 gmax = uparam(iadb+1)*uparam(iadb+6)
200 . + uparam(iadb+2)*uparam(iadb+7)
201 . + uparam(iadb+3)*uparam(iadb+8)
202 . + uparam(iadb+4)*uparam(iadb+9)
203 . + uparam(iadb+5)*uparam(iadb+10)
204 rho = pm(1,imt)
205 a11 = gmax*(one + nu)/(one - nu**2)
206 ssp(i)=
max(ssp(i), sqrt(a11/rho))
207 ENDDO
208 ENDDO
209 ELSEIF (mtn == 69) THEN
210 DO i=lft,llt
211 DO n=1,nip
212 imt = igeo_stack(ipmat + n,isubstack)
213 iadb = ipm(7,imt)-1
214 nu = uparam(iadb+14)
215 gmax = uparam(iadb+1)*uparam(iadb+6)
216 . + uparam(iadb+2)*uparam(iadb+7)
217 . + uparam(iadb+3)*uparam(iadb+8)
218 . + uparam(iadb+4)*uparam(iadb+9)
219 . + uparam(iadb+5)*uparam(iadb+10)
220 rho = pm(1,imt)
221 a11 = gmax*(one + nu)/(one - nu**2)
222 ssp(i)=
max(ssp(i), sqrt(a11/rho))
223 ENDDO
224 ENDDO
225 ELSEIF (mtn == 65) THEN
226 DO i=lft,llt
227 DO n=1,nip
228 imt = igeo_stack(ipmat + n,isubstack)
229 rho =pm(1,imt)
230 young=pm(20,imt)
231 ssp(i)=
max(ssp(i), sqrt(young/rho))
232 ENDDO
233 ENDDO
234 ELSE
235 DO i=lft,llt
236 DO n=1,nip
237 imt = igeo_stack(ipmat + n,isubstack)
238 rho =pm(1,imt)
239 young=pm(20,imt)
240 nu =pm(21,imt)
241 ssp(i)=
max(ssp(i), sqrt(young/(one-nu*nu)/rho))
242 ENDDO
243 ENDDO
244 ENDIF
245 ELSEIF(igtyp == 51 .AND. igmat < 0) THEN
246 nip = nlay
247 ipmat = 2 + nlay
248 DO i=lft,llt
249 DO n=1,nip
250 imt = igeo_stack(ipmat + n,isubstack)
251 mlawly = nint(pm(19,imt))
252 IF (mlawly <= 28) THEN
253 ssp(i)=
max(ssp(i),pm(27,imt))
254 ELSEIF (mlawly == 42) THEN
255 rho = pm(1,imt)
256 nu = pm(21,imt)
257 gmax = pm(22,imt)
258 a11 = gmax*(one + nu)/(one - nu**2)
259 ssp(i)=
max(ssp(i), sqrt(a11/rho))
260 ELSEIF (mlawly == 69) THEN
261 iadb = ipm(7,imt)-1
262 nu = uparam(iadb+14)
263 gmax = uparam(iadb+1)*uparam(iadb+6)
264 . + uparam(iadb+2)*uparam(iadb+7)
265 . + uparam(iadb+3)*uparam(iadb+8)
266 . + uparam(iadb+4)*uparam(iadb+9)
267 . + uparam(iadb+5)*uparam(iadb+10)
268 rho = pm(1,imt)
269 a11 = gmax*(one + nu)/(one - nu**2)
270 ssp(i)=
max(ssp(i), sqrt(a11/rho))
271 ELSEIF (mlawly == 65) THEN
272 rho =pm(1,imt)
273 young=pm(20,imt)
274 ssp(i)=
max(ssp(i), sqrt(young/rho))
275 ELSE
276 rho =pm(1,imt)
277 young=pm(20,imt)
278 nu =pm(21,imt)
279 ssp(i)=
max(ssp(i), sqrt(young/(one-nu*nu)/rho))
280 ENDIF
281 ENDDO
282 ENDDO
283
284 ELSEIF (mtn<=28)THEN
285 DO i=lft,llt
286 ssp(i)=pm(27,imat)
287 ENDDO
288 ELSEIF (mtn == 42) THEN
289 DO i=lft,llt
290 rho = pm(1,imat)
291 nu = pm(21,imat)
292 gmax = pm(22,imat)
293 a11 = gmax*(one + nu)/(one - nu**2)
294 ssp(i)=
max(ssp(i), sqrt(a11/rho))
295 ENDDO
296 ELSEIF (mtn == 69) THEN
297 DO i=lft,llt
298 iadb = ipm(7,imat)-1
299 nu = uparam(iadb+14)
300 gmax = uparam(iadb+1)*uparam(iadb+6)
301 . + uparam(iadb+2)*uparam(iadb+7)
302 . + uparam(iadb+3)*uparam(iadb+8)
303 . + uparam(iadb+4)*uparam(iadb+9)
304 . + uparam(iadb+5)*uparam(iadb+10)
305 rho = pm(1,imat)
306 a11 = gmax*(one + nu)/(one - nu**2)
307 ssp(i)=
max(ssp(i), sqrt(a11/rho))
308 ENDDO
309 ELSEIF (mtn == 65) THEN
310 DO i=lft,llt
311 rho =pm(1,imat)
312 young =pm(20,imat)
313 ssp(i)=sqrt(young/rho)
314 ENDDO
315 ELSE
316 DO i=lft,llt
317 rho =pm(1,imat)
318 young=pm(20,imat)
319 nu =pm(21,imat)
320 ssp(i)=sqrt(young/(one-nu*nu)/rho)
321 ENDDO
322 ENDIF
323
324 DO 20 i=lft,llt
325 al1(i)= x2l(i) * x2l(i) + y2l(i) * y2l(i)
326 al2(i)=(x3l(i)-x2l(i))*(x3l(i)-x2l(i))+(y3l(i)-y2l(i))*(y3l(i)-y2l(i))
327 al3(i)=(x4l(i)-x3l(i))*(x4l(i)-x3l(i))+(y4l(i)-y3l(i))*(y4l(i)-y3l(i))
328 al4(i)= x4l(i) * x4l(i) + y4l(i) * y4l(i)
329 al5(i)=(x4l(i)-x2l(i))*(x4l(i)-x2l(i))+(y4l(i)-y2l(i))*(y4l
330 al6(i)= x3l(i) * x3l(i) + y3l(i) * y3l(i)
331 20 CONTINUE
332
333 DO 30 i=lft,llt
334 al(i)=
min(al1(i),al2(i),al3(i),al4(i),al5(i),al6(i))
335 IF(al3(i) == zero) al(i)=
min(al1(i),al2(i),al4(i))
336 almin(i)=sqrt(al(i))
337 30 CONTINUE
338
339 IF(mtn == 19)THEN
340 viscdef=fourth
341 ELSEIF(mtn == 25.OR.mtn == 27)THEN
342 viscdef=fiveem2
343 ELSE
344 viscdef=zero
345 ENDIF
346
347 viscmx = group_param%VISC_DM
348 IF (viscmx == zero) viscmx = viscdef
349 IF (mtn == 1 .OR.mtn == 2.OR.mtn == 3.OR.
350 . mtn == 22.OR.mtn == 23) viscmx=zero
351 viscmx = sqrt(one + viscmx*viscmx) - viscmx
352 DO i=lft,llt
353 dtdyn =
area(i)/sqrt(
max(al5(i),al6(i)))
354 aldt(i) =
max(dtdyn,almin(i))
355 dt(i) = aldt(i)*viscmx/ssp(i)
356 ENDDO
357
358
359
360 ipgmat = 700
361 IF(nadmesh == 0)THEN
362 IF (igtyp == 11 .AND. igmat > 0) THEN
363 DO i=lft,llt
364 a11 =geo(ipgmat + 5,iprop)
365 a11r =geo(ipgmat + 7,iprop)
366 b1 = px1(i)*px1(i)+py1(i)*py1(i)
367 b2 = px2(i)*px2(i)+py2(i)*py2(i)
368 vv = viscmx * viscmx
369 fac =
max(b1,b2) / (
area(i) * vv)
370 sti = fac * thk(i) * a11
371 stir = fac*a11r * thk(i)*(thk(i)**2 +
area(i))*one_over_12
372 stifn(ixc(2,i))=stifn(ixc(2,i))+sti
373 stifn(ixc(3,i))=stifn(ixc(3,i))+sti
374 stifn(ixc(4,i))=stifn(ixc(4,i))+sti
375 stifn(ixc(5,i))=stifn(ixc(5,i))+sti
376 stifr(ixc(2,i))=stifr(ixc(2,i))+stir
377 stifr(ixc(3,i))=stifr(ixc(3,i))+stir
378 stifr(ixc(4,i))=stifr(ixc(4,i))+stir
379 stifr(ixc(5,i))=stifr(ixc(5,i))+stir
380 strc(i) = stir
381 ENDDO
382 ELSEIF(igtyp == 52 .OR.
383 . ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0 )) THEN
384 DO i=lft,llt
385 a11 = pm_stack(5 ,isubstack)
386 a11r = pm_stack(7 ,isubstack)
387 b1 = px1(i)*px1(i)+py1(i)*py1(i)
388 b2 = px2(i)*px2(i)+py2(i)*py2(i)
389 vv = viscmx * viscmx
390 fac =
max(b1,b2) / (
area(i) * vv)
391 sti = fac * thk(i) * a11
392 stir = fac*a11r * thk(i)*((thk(i)**2 +
area(i))*one_over_12 +
393 . zoffset(i)*zoffset(i) )
394 stifn(ixc(2,i))=stifn(ixc(2,i))+sti
395 stifn(ixc(3,i))=stifn(ixc(3,i))+sti
396 stifn(ixc(4,i))=stifn(ixc(4,i))+sti
397 stifn(ixc(5,i))=stifn(ixc(5,i))+sti
398 stifr(ixc(2,i))=stifr(ixc(2,i))+stir
399 stifr(ixc(3,i))=stifr(ixc(3,i))+stir
400 stifr(ixc(4,i))=stifr(ixc(4,i))+stir
401 stifr(ixc(5,i))=stifr(ixc(5,i))+stir
402 strc(i) = stir
403 ENDDO
404 ELSE
405 DO i=lft,llt
406 a11 =pm(24,imat)
407 b1 = px1(i)*px1(i)+py1(i)*py1(i)
408 b2 = px2(i)*px2(i)+py2(i)*py2(i)
409 vv = viscmx * viscmx
411 . * thk(i) * a11 / (
area(i) * vv)
412 stir = sti * (thk(i)*thk(i) +
area(i)) / 12.
413 stifn(ixc(2,i))=stifn(ixc(2,i))+sti
414 stifn(ixc(3,i))=stifn(ixc(3,i))+sti
415 stifn(ixc(4,i))=stifn(ixc(4,i))+sti
416 stifn(ixc(5,i))=stifn(ixc(5,i))+sti
417 stifr(ixc(2,i))=stifr(ixc(2,i))+stir
418 stifr(ixc(3,i))=stifr(ixc(3,i))+stir
419 stifr(ixc(4,i))=stifr(ixc(4,i))+stir
420 stifr(ixc(5,i))=stifr(ixc(5,i))+stir
421 strc(i) = stir
422 ENDDO
423 ENDIF
424 ELSE
425 IF(igtyp == 11 .AND. igmat > 0) THEN
426 DO i=lft,llt
427 n=nft+i
428 IF(sh4tree(3,n) >= 0)THEN
429 a11 =geo(ipgmat + 5,iprop)
430 a11r =geo(ipgmat + 7,iprop)
431 b1 = px1(i)*px1(i)+py1(i)*py1(i)
432 b2 = px2(i)*px2(i)+py2(i)*py2(i)
433 vv = viscmx * viscmx
434 fac =
max(b1,b2) / (
area(i) * vv)
435 sti = fac * thk(i) * a11
436 stir = fac * a11r * thk(i)*(thk(i)**2 +
area(i))*one_over_12
437 stifn(ixc(2,i))=stifn(ixc(2,i))+sti
438 stifn(ixc(3,i))=stifn(ixc(3,i))+sti
439 stifn(ixc(4,i))=stifn(ixc(4,i))+sti
440 stifn(ixc(5,i))=stifn(ixc(5,i))+sti
441 stifr(ixc(2,i))=stifr(ixc(2,i))+stir
442 stifr(ixc(3,i))=stifr(ixc(3,i))+stir
443 stifr(ixc(4,i))=stifr(ixc(4,i))+stir
444 stifr(ixc(5,i))=stifr(ixc(5,i))+stir
445 strc(i) = stir
446 END IF
447 END DO
448 ELSEIF(igtyp == 52 .OR.
449 . ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0 )) THEN
450 DO i=lft,llt
451 n=nft+i
452 IF(sh4tree(3,n) >= 0)THEN
453 a11 = pm_stack(5 ,isubstack)
454 a11r = pm_stack(7 ,isubstack)
455 b1 = px1(i)*px1(i)+py1(i)*py1(i)
456 b2 = px2(i)*px2(i)+py2(i)*py2(i)
457 vv = viscmx * viscmx
458 fac =
max(b1,b2) / (
area(i) * vv)
459 sti = fac * thk(i) * a11
460 stir = fac * a11r * thk(i)*((thk(i)**2 +
area(i))*one_over_12 +
461 . zoffset(i)*zoffset(i) )
462 stifn(ixc(2,i))=stifn(ixc(2,i))+sti
463 stifn(ixc(3,i))=stifn(ixc(3,i))+sti
464 stifn(ixc(4,i))=stifn(ixc(4,i))+sti
465 stifn(ixc(5,i))=stifn(ixc(5,i))+sti
466 stifr(ixc(2,i))=stifr(ixc(2,i))+stir
467 stifr(ixc(3,i))=stifr(ixc(3,i))+stir
468 stifr(ixc(4,i))=stifr(ixc(4,i))+stir
469 stifr(ixc(5,i))=stifr(ixc(5,i))+stir
470 strc(i) = stir
471 END IF
472 END DO
473
474 ELSE
475 DO i=lft,llt
476 n=nft+i
477 IF(sh4tree(3,n) >= 0)THEN
478 a11 =pm(24,imat)
479 b1 = px1(i)*px1(i)+py1(i)*py1(i)
480 b2 = px2(i)*px2(i)+py2(i)*py2(i)
481 vv = viscmx * viscmx
483 . * thk(i) * a11 / (
area(i) * vv)
484 stir = sti * (thk(i)*thk(i) +
area(i)) / 12.
485 stifn(ixc(2,i))=stifn(ixc(2,i))+sti
486 stifn(ixc(3,i))=stifn(ixc(3,i))+sti
487 stifn(ixc(4,i))=stifn(ixc(4,i))+sti
488 stifn(ixc(5,i))=stifn(ixc(5,i))+sti
489 stifr(ixc(2,i))=stifr(ixc
490 stifr(ixc(3,i))=stifr(ixc(3,i))+stir
491 stifr(ixc(4,i))=stifr(ixc(4,i))+stir
492 stifr(ixc(5,i))=stifr(ixc(5,i))+stir
493 strc(i) = stir
494 END IF
495 END DO
496 ENDIF
497 END IF
498
499 IF(ismstr == 3)THEN
500 DO i=lft,llt
501 IF(geo(5,iprop)/=zero
502 ENDDO
503 ELSE
504 DO i=lft,llt
505 px1(i)= zero
506 px2(i)= zero
507 py1(i)= zero
508 py2(i)= zero
509 ENDDO
510 ENDIF
511
512 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)
1.15.0