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