37
38
39
40 USE group_param_mod
41
42
43
44#include "implicit_f.inc"
45
46
47
48#include "mvsiz_p.inc"
49
50
51
52#include "param_c.inc"
53#include "remesh_c.inc"
54#include "vect01_c.inc"
55
56
57
58 INTEGER JFT, JLT, ILEV,ISUBSTACK,IMAT,IPROP
59 INTEGER IXTG(NIXTG,*), SH3TREE(KSH3TREE,*),(NPROPMI,*),
60 . IGEO(NPROPGI,*),PM_STACK(20,*)
62 . pm(npropm,*), geo(npropg,*), px2(*),px3(*),py2(*),py3(*),
63 . stifn(*),stifr(*),thk(*),aldt(*),uparam(*),strtg(*),
64 . x1g(mvsiz), x2g(mvsiz), x3g(mvsiz),
65 . y1g(mvsiz), y2g(mvsiz), y3g(mvsiz),
66 . z1g(mvsiz), z2g(mvsiz), z3g(mvsiz),
67 . e1x(mvsiz), e1y(mvsiz), e1z(mvsiz),
68 . e2x(mvsiz), e2y(mvsiz), e2z(mvsiz),
69 . e3x(mvsiz), e3y(mvsiz), e3z(mvsiz)
70 TYPE (GROUP_PARAM_) :: GROUP_PARAM
71
72
73
74 INTEGER I, NG, N,IADB,I1,I3,IPTHK,IPPOS,IGTYP,I2,
75 . MATLY,IGMAT,IPGMAT,
77 . x21g(mvsiz), y21g(mvsiz), z21g(mvsiz),
78 . x31g(mvsiz), y31g(mvsiz), z31g(mvsiz),
79 . x32g(mvsiz), y32g(mvsiz), z32g(mvsiz),
80 . x2(mvsiz), x3(mvsiz), y2(mvsiz),y3(mvsiz),
81 . dt(mvsiz),
area(mvsiz),fac, almin,
82 . viscmx, a11, g, sti,stir,shf,viscdef,gmax,
83 . al1, al2, al3, almax, ssp,young,nu,rho,areai,
84 . c1,iz,thickt,thkly,posly,a1thk,c1thk,
85 . gthk,a11r,a12,e,ethk,nuthk,a12thk,rhog
86
87 DO i=jft,jlt
88 x21g(i) = x2g(i)-x1g(i)
89 y21g(i) = y2g(i)-y1g(i)
90 z21g(i) = z2g(i)-z1g(i)
91 x31g(i) = x3g(i)-x1g(i)
92 y31g(i) = y3g(i)-y1g(i)
93 z31g(i) = z3g(i)-z1g(i)
94 x32g(i) = x3g(i)-x2g(i)
95 y32g(i) = y3g(i)-y2g(i)
96 z32g(i) = z3g(i)-z2g(i)
97 ENDDO
98
99 DO i=jft,jlt
100 x2(i)=e1x(i)*x21g(i)+e1y(i)*y21g(i)+e1z(i)*z21g(i)
101 y2(i)=e2x(i)*x21g(i)+e2y(i)*y21g(i)+e2z(i)*z21g(i)
102 y3(i)=e2x(i)*x31g(i)+e2y(i)*y31g(i)+e2z(i)*z31g(i)
103 x3(i)=e1x(i)*x31g(i)+e1y(i)*y31g(i)+e1z(i)*z31g(i)
104 ENDDO
105
106
107
108 igtyp = igeo(11,iprop)
109 igmat = igeo(98,iprop)
110 ipgmat = 700
111
112 IF(mtn == 19)THEN
113 viscdef=fourth
114 ELSEIF(mtn == 25.OR.mtn == 27)THEN
115 viscdef=fiveem2
116 ELSE
117 viscdef=zero
118 ENDIF
119
120 DO 40 i=jft,jlt
121 al1 = x2(i) * x2(i) + y2(i) * y2(i)
122 al2 = (x3(i)-x2(i)) * (x3(i)-x2(i)) +
123 . (y3(i)-y2(i)) * (y3(i)-y2(i))
124 al3 = x3(i) * x3(i) + y3(i) * y3(i)
125 almax =
max(al1,al2,al3)
126 nu =pm(21,imat)
127 almin =
min(al1,al2,al3)
128 fac =one+zep6*(1+nu)*thk(i)*thk(i)/almin
129 almax = almax*fac
130 IF(igtyp == 11 .AND. igmat > 0) THEN
131 ssp = geo(ipgmat +9 ,iprop)
132 ELSEIF(igtyp == 52 .OR.
133 . ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0)) THEN
134 ssp = pm_stack(9 ,isubstack)
135 ELSE
136 IF(mtn<=28)THEN
137 ssp=pm(27,imat)
138 ELSEIF (mtn == 42) THEN
139 rho = pm(1 ,imat)
140 nu = pm(21,imat)
141 gmax = pm(22,imat)
142 a11 = gmax*(one + nu)/(one - nu**2)
143 ssp =
max(ssp, sqrt(a11/rho))
144 ELSEIF (mtn == 69) THEN
145 iadb = ipm(7,imat)-1
146 nu = uparam(iadb
147 gmax = uparam(iadb+1)*uparam(iadb+6)
148 . + uparam(iadb+2)*uparam(iadb+7)
149 . + uparam(iadb+3)*uparam(iadb+8)
150 . + uparam(iadb+4)*uparam(iadb+9)
151 . + uparam(iadb+5)*uparam(iadb+10)
152 rho = pm(1,imat)
153 a11 = gmax*(one + nu)/(one - nu**2)
154 ssp =
max(ssp, sqrt(a11/rho))
155 ELSEIF (mtn == 65) THEN
156 rho =pm(1,imat)
157 young=pm(20,imat)
158 ssp=sqrt(young/rho)
159 ELSE
160 rho =pm(1,imat)
161 young=pm(20,imat)
162 nu =pm(21,imat)
163 ssp=sqrt(young/(one-nu*nu)/rho)
164 ENDIF
165 ENDIF
166 viscmx = group_param%VISC_DM
167 IF (viscmx == zero) viscmx = viscdef
168 IF(mtn == 1.OR.mtn == 2.OR.mtn == 3.OR.
169 . mtn == 22.OR.mtn == 23)viscmx=zero
170 viscmx=sqrt(1.+viscmx*viscmx)-viscmx
171 aldt(i)= two*
area(i)*viscmx / sqrt(almax)
172 dt(i) = aldt(i) / ssp
173 40 CONTINUE
174
175
176
177 ipgmat = 700
178 IF(nadmesh==0)THEN
179 IF(igtyp == 11 .AND. igmat > 0) THEN
180 DO i=jft,jlt
181 a11 = geo(ipgmat + 5 ,iprop)
182 a11r = geo(ipgmat + 7 ,iprop)
183 g = geo(ipgmat + 4 ,iprop)
184 fac =
area(i)* thk(i) / (aldt(i))**2
185 sti = fac * a11
186 stir =one_over_12*fac* a11r*thk(i)**2
187 stifn(ixtg(2,i))=stifn(ixtg(2,i))+sti
188 stifn(ixtg(3,i))=stifn(ixtg(3,i))+sti
189 stifn(ixtg(4,i))=stifn(ixtg(4,i))+sti
190 stifr(ixtg(2,i))=stifr(ixtg(2,i))+stir
191 stifr(ixtg(3,i))=stifr(ixtg(3,i))+stir
192 stifr(ixtg(4,i))=stifr(ixtg(4,i))+stir
193 strtg(i) = stir
194 END DO
195 ELSEIF(igtyp == 52 .OR.
196 . ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0)) THEN
197 DO i=jft,jlt
198 a11 = pm_stack(5 ,isubstack)
199 a11r = pm_stack(7 ,isubstack)
200 g = pm_stack(4 ,isubstack)
201 fac =
area(i)* thk(i) / (aldt(i))**2
202 sti = fac * a11
203 stir =one_over_12*fac* a11r*thk(i)**2
204 stifn(ixtg(2,i))=stifn(ixtg(2,i))+sti
205 stifn(ixtg(3,i))=stifn(ixtg(3,i))+sti
206 stifn(ixtg(4,i))=stifn(ixtg(4,i))+sti
207 stifr(ixtg(2,i))=stifr(ixtg(2,i))+stir
208 stifr(ixtg(3,i))=stifr(ixtg(3,i))+stir
209 stifr(ixtg(4,i))=stifr(ixtg(4,i))+stir
210 strtg(i) = stir
211 END DO
212 ELSE
213 DO i=jft,jlt
214 a11 =geo(ipgmat +5 ,iprop)
215 a11r =geo(ipgmat +7 ,iprop)
216 g =geo(ipgmat +4 ,iprop)
217 fac =
area(i)* thk(i) / (aldt(i))**2
218 sti =fac* a11
219 stir =one_over_12*fac* a11r*thk(i)**2
220 stifn(ixtg(2,i))=stifn(ixtg(2,i))+sti
221 stifn(ixtg(3,i))=stifn(ixtg(3,i))+sti
222 stifn(ixtg(4,i))=stifn(ixtg(4,i))+sti
223 stifr(ixtg(2,i))=stifr(ixtg(2,i))+stir
224 stifr(ixtg(3,i))=stifr(ixtg(3,i))+stir
225 stifr(ixtg(4,i))=stifr(ixtg(4,i))+stir
226 strtg(i) = stir
227 END DO
228 ENDIF
229 ELSE
230 IF(igtyp == 11 .AND. igmat > 0 )THEN
231 DO i=jft,jlt
232 n=nft+i
233 IF(sh3tree(3,n) >= 0)THEN
234 a11 =geo(ipgmat +5 ,iprop)
235 a11r =geo(ipgmat +7 ,iprop)
236 g =geo(ipgmat +4 ,iprop)
237
238!
239 fac =
area(i)* thk(i) / (aldt(i))**2
240 sti =fac* a11
241 stir =one_over_12*fac* a11r*thk(i)**2
242 stifn(ixtg(2,i))=stifn(ixtg(2,i))+sti
243 stifn(ixtg(3,i))=stifn(ixtg(3,i))+sti
244 stifn(ixtg(4,i))=stifn(ixtg(4,i))+sti
245 stifr(ixtg(2,i))=stifr(ixtg(2,i))+stir
246 stifr(ixtg(3,i))=stifr(ixtg(3,i))+stir
247 stifr(ixtg(4,i))=stifr(ixtg(4,i))+stir
248 strtg(i) = stir
249 END IF
250 END DO
251 ELSEIF(igtyp == 52 .OR.
252 . ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0)) THEN
253 DO i=jft,jlt
254 n=nft+i
255 IF(sh3tree(3,n) >= 0)THEN
256 a11 = pm_stack(5 ,isubstack)
257 a11r = pm_stack(7 ,isubstack)
258 g = pm_stack(4 ,isubstack)
259 fac =
area(i)* thk(i) / (aldt(i))**2
260 sti =fac* a11
261 stir =one_over_12*fac* a11r*thk(i)**2
262 stifn(ixtg(2,i))=stifn(ixtg(2,i))+sti
263 stifn(ixtg(3,i))=stifn(ixtg(3,i))+sti
264 stifn(ixtg(4,i))=stifn(ixtg(4,i))+sti
265 stifr(ixtg(2,i))=stifr(ixtg(2,i))+stir
266 stifr(ixtg(3,i))=stifr(ixtg(3,i))+stir
267 stifr(ixtg(4,i))=stifr(ixtg(4,i))+stir
268 strtg(i) = stir
269 END IF
270 END DO
271 ELSE
272 DO i=jft,jlt
273 n=nft+i
274 IF(sh3tree(3,n) >= 0)THEN
275 a11 =pm(24,imat)
276 g =pm(22,imat)
277 sti =
area(i) * thk(i) * a11 / (aldt(i))**2
278 stir = sti * thk(i) * thk(i) / 12.
279 stifn(ixtg(2,i))=stifn(ixtg(2,i))+sti
280 stifn(ixtg(3,i))=stifn(ixtg(3,i))+sti
281 stifn(ixtg(4,i))=stifn(ixtg(4,i))+sti
282 stifr(ixtg(2,i))=stifr(ixtg(2,i))+stir
283 stifr(ixtg(3,i))=stifr(ixtg(3,i))+stir
284 stifr(ixtg(4,i))=stifr(ixtg(4,i))+stir
285 strtg(i) = stir
286 END IF
287 END DO
288 ENDIF
289 END IF
290
291
292 IF(ismstr/=3)THEN
293 DO 50 i=jft,jlt
294 px2(i) = zero
295 py2(i) = zero
296 px3(i) = zero
297 py3(i) = zero
298 50 CONTINUE
299 ELSE
300
301
302 DO i=jft,jlt
304 px2(i)=y3(i)*areai
305 py2(i)=-x3(i)*areai
306 px3(i)=-y2(i)*areai
307 py3(i)=x2(i)*areai
308 ENDDO
309
310 DO 80 i=jft,jlt
311 ng=iprop
312 IF (geo(5,ng) == zero) GOTO 80
313 geo(5,ng)=
min(geo(5,ng),dt(i))
314 80 CONTINUE
315 ENDIF
316
317
318 RETURN
319
subroutine area(d1, x, x2, y, y2, eint, stif0)