44 use element_mod , only : nixc
45
46
47
48#include "implicit_f.inc"
49#include "comlock.inc"
50
51
52
53#include "mvsiz_p.inc"
54
55
56
57#include "com04_c.inc"
58#include "com06_c.inc"
59#include "param_c.inc"
60#include "remesh_c.inc"
61#include "scr02_c.inc"
62#include "scr06_c.inc"
63#include "scr14_c.inc"
64#include "scr16_c.inc"
65#include "scr17_c.inc"
66#include "sms_c.inc"
67
68
69
70 INTEGER IXC(NIXC,*),IPARTC(*), JFT, JLT,
71 . IHBE ,NFT ,ISMSTR,,IGTYP,
72 . IGMAT,NEL
73
75 . thk(*), hour(nel,5), off(*),partsav(npsav,*),
76 . px1(*), px2(*), py1(*), py2(*),dt1c(*),eani(*),
77 . ssp(mvsiz), rho(mvsiz),sti(mvsiz),stir(mvsiz),
78 . vx1(mvsiz), vx2(mvsiz), vx3(mvsiz), vx4(mvsiz),
79 . vy1(mvsiz), vy2(mvsiz), vy3(mvsiz), vy4(mvsiz),
80 . vz1(mvsiz), vz2(mvsiz), vz3(mvsiz), vz4(mvsiz),
81 .
area(mvsiz),thk0(mvsiz),vhx(mvsiz), vhy(mvsiz),
82 . shf(mvsiz),z2(mvsiz),viscmx(mvsiz),g(mvsiz),
83 . h1(mvsiz), h2(mvsiz), h3(mvsiz), a11(mvsiz),
84 . b11(mvsiz), b12(mvsiz), b13(mvsiz), b14(mvsiz),
85 . b21(mvsiz), b22(mvsiz), b23(mvsiz), b24(mvsiz),
86 . h11(mvsiz), h12(mvsiz), h13(mvsiz),
87 . h21(mvsiz), h22(mvsiz), h23(mvsiz),
88 . h31(mvsiz), h32(mvsiz), h33(mvsiz),
89 . rx1(mvsiz), rx2(mvsiz), rx3(mvsiz), rx4(mvsiz),
90 . ry1(mvsiz), ry2(mvsiz), ry3(mvsiz), ry4(mvsiz),
91 . thk02(mvsiz),ym(mvsiz), nu(mvsiz), alpe(mvsiz),
92 . srh1(*) ,srh2(*) ,srh3(*),a11r(*)
93
94
95
96 INTEGER I,MX,J, II, IC, JST(MVSIZ+1)
97
99 . h1l(mvsiz), h2l(mvsiz), h3l(mvsiz), h1q(mvsiz), h2q(mvsiz),
100 . h3q(mvsiz), hg1(mvsiz), hg2(mvsiz),
101 . hh1(mvsiz), hh2(mvsiz), hh3(mvsiz),
102 . a1(mvsiz), a2(mvsiz), a3(mvsiz), a4(mvsiz),
103 . a5(mvsiz), a6(mvsiz), a7(mvsiz), a8(mvsiz),
104 . b1(mvsiz), b2(mvsiz), ehou(mvsiz),
105 .
106 . gama1(mvsiz), gama2(mvsiz), gama3(mvsiz), gama4(mvsiz)
108 . px1v, px2v, py1v, py2v, hvish1, hvish2,
109 . shfpr3, fac, ehourt,
110 . r0,r1,t2a,tsa,td,vv,
111 . hour1a,hour2a,hour3a,
112 .
113 . sr2d2,srshfpr3, inv9, inv12,scale(mvsiz)
114
115 ehourt = zero
116 sr2d2 = sqrt(two)* half
117 inv9 = one_over_9
118 inv12 = one_over_12
119
120
121
122 IF(ismstr/=1.AND.ismstr/=11.AND.ihbe>=1)THEN
123 DO i=jft,jlt
124 px1v = px1(i)*vhx(i)
125 px2v = px2(i)*vhx(i)
126 py1v = py1(i)*vhy(i)
127 py2v = py2(i)*vhy(i)
128 gama1(i)= off(i)*( one - px1v-py1v)
129 gama3(i)= off(i)*( one + px1v+py1v)
130 gama2(i)= off(i)*(-one - px2v-py2v)
131 gama4(i)= off(i)*(-one + px2v+py2v)
132 ENDDO
133 ELSE
134 DO i=jft,jlt
135 gama1(i)= off(i)
136 gama3(i)= off(i)
137 gama2(i)= -off(i)
138 gama4(i)= -off(i)
139 ENDDO
140 ENDIF
141
142 DO i=jft,jlt
143 shfpr3 = shf(i)/(three*(one + nu(i)))
144 hvish1 = hvisc*h1(i)
145 hvish2 = hvisc*h2(i)
146 r0 = fourth*rho(i)
147 r1 = r0*hundred
148 r0 = r0*hvlin
149 a1(i) = r1*hvish1
150 a2(i) = r0*sr2d2*srh1(i)
151 srshfpr3 = sqrt(shfpr3)
152 a3(i) = r1*hvish2*srshfpr3
153 a4(i) = r0*sr2d2*srh2(i)*srshfpr3
154 hh3(i) = helas*h3(i)
155 a5(i) = hh3(i)*r1*zep072169
156 hh3(i) = sr2d2*srh3(i)
157 a6(i) = hh3(i)*r0*zep072169
158 r0 = fourth*ym(i)*helas
159 a7(i) = h1(i)*r0
160 a8(i) = h2(i)*r0*shfpr3
161 ENDDO
162
163 DO i=jft,jlt
164 t2a = thk02(i)*
area(i)
165 tsa = sqrt(t2a)
166 h1q(i) = a1(i)*tsa
167 h1l(i) = a2(i)*ssp(i)*tsa
168 h2q(i) = a3(i)*thk02(i)
169 h2l(i) = a4(i)*ssp(i)*thk02(i)
170 h3q(i) = a5(i)*t2a
171 h3l(i) = a6(i)*ssp(i)*t2a
172 td = thk0(i)*dt1c(i)
173 hh1(i) = a7(i)*td
174 b1(i) = px1(i)*px1(i)+py1(i)*py1(i)
175 b2(i) = px2(i)*px2(i)+py2(i)*py2(i)
176 hh2(i) = a8(i)*thk02(i)*td/(b1(i)+b2(i))
177 ENDDO
178
179
180
181 DO i=jft,jlt
182 IF(ixc(4,i)==ixc(5,i))THEN
183 h1q(i)=zero
184 h1l(i)=zero
185 h2q(i)=zero
186 h2l(i)=zero
187 h3q(i)=zero
188 h3l(i)=zero
189 hh1(i)=zero
190 hh2(i)=zero
191 END IF
192 END DO
193
194
195
196 DO i=jft,jlt
197 scale(i) = zero
198 ENDDO
199 IF(nodadt/=0.OR.idt1sh==1.OR.idtmins==2)THEN
200 DO i=jft,jlt
201 scale(i) =
max(gama1(i)*gama1(i),gama2(i)*gama2(i),
202 . gama3(i)*gama3(i),gama4(i)*gama4(i)) *
203 . dt1c(i)*
max(hh1(i)+h1l(i),hh2(i)+h2l(i),h3l(i))
204 . /
max(dt1c(i)*dt1c(i),em20)
205 sti(i) = sti(i) + scale(i)
206 ENDDO
207
208
209 IF(igtyp == 52 .OR.
210 . ((igtyp == 11 .OR. igtyp == 17 .OR. igtyp == 51)
211 . .AND. igmat > 0 )) THEN
212 IF(nadmesh==0)THEN
213 DO i=jft,jlt
214 IF (off(i)==zero) THEN
215 sti(i) = zero
216 stir(i) = zero
217 ELSE
218 vv = viscmx(i) * viscmx(i) * alpe(i)
219 fac =
max(b1(i),b2(i))/(
area(i) * vv)
220 sti(i) = sti(i) + fac * thk0(i) * a11(i)
221 stir(i) = fac * a11r(i)*one_over_12*thk0(i)**3 +
222 . fac * a11(i)*thk0(i)*
area(i)*inv9 +
223 . fac*scale(i)*(one_over_12*thk0(i)**2 +
area(i)*inv9)
224 ENDIF
225 ENDDO
226 ELSE
227 DO i=jft,jlt
228 IF (off(i)==zero) THEN
229 sti(i) = zero
230 stir(i) = zero
231 ELSE
232 vv = viscmx(i) * viscmx(i) * alpe(i)
233 fac =
max(b1(i),b2(i))/(
area(i) * vv)
234 sti(i) = sti(i) + fac * thk0(i) * a11(i)
235 stir(i) = fac * a11r(i)*one_over_12*thk0(i)**3 +
236 . fac*scale(i)*one_over_12*thk0(i)**2
237 ENDIF
238 ENDDO
239 END IF
240 ELSE
241 IF(nadmesh==0)THEN
242 DO i=jft,jlt
243 IF (off(i)==zero) THEN
244 sti(i) = zero
245 stir(i) = zero
246 ELSE
247 vv = viscmx(i) * viscmx(i) * alpe(i)
248 sti(i) = sti(i) +
max(b1(i),b2(i))
249 . * thk0(i) * a11(i) / (
area(i) * vv)
250 stir(i) = sti(i) * (thk02(i)*inv12 +
area(i)*inv9)
251
252 ENDIF
253 ENDDO
254 ELSE
255 DO i=jft,jlt
256 IF (off(i)==zero) THEN
257 sti(i) = zero
258 stir(i) = zero
259 ELSE
260 vv = viscmx(i) * viscmx(i) * alpe(i)
261 sti(i) = sti(i) +
max(b1(i),b2(i))
262 . * thk0(i) * a11(i) / (
area(i) * vv)
263 stir(i) = sti(i) * thk02(i)*inv12
264 ENDIF
265 ENDDO
266 END IF
267 endif
268 ENDIF
269
270
271
272 IF(ismstr==1.OR.ismstr==11.OR.ihbe<1)THEN
273 DO i=jft,jlt
274 hg1(i)=(vx1(i)-vx2(i)+vx3(i)-vx4(i))*off(i)
275 hg2(i)=(vy1(i)-vy2(i)+vy3(i)-vy4(i))*off(i)
276 ENDDO
277 ELSE
278 DO i=jft,jlt
279 hg1(i)=vx1(i)*gama1(i)+vx2(i)*gama2(i)
280 . +vx3(i)*gama3(i)+vx4(i)*gama4(i)
281 hg2(i)=vy1(i)*gama1(i)+vy2(i)*gama2(i)
282 . +vy3(i)*gama3(i)+vy4(i)*gama4(i)
283 ENDDO
284 ENDIF
285 DO i=jft,jlt
286 hour(i,1)=hour(i,1)+hg1(i)*hh1(i)
287 hour(i,2)=hour(i,2)+hg2(i)*hh1(i)
288 hour1a =hour(i,1)+hg1(i)*(h1l(i)+h1q(i)*abs(hg1(i)))
289 h11(i)=hour1a*gama1(i)
290 h12(i)=hour1a*gama2(i)
291 h13(i)=hour1a*gama3(i)
292
293 hour2a =hour(i,2)+hg2(i)*(h1l(i)+h1q(i)*abs(hg2(i)))
294 h21(i)=hour2a*gama1(i)
295 h22(i)=hour2a*gama2(i)
296 h23(i)=hour2a*gama3(i)
297
298 ehou(i) = hour1a*hg1(i) + hour2a*hg2(i)
299 ENDDO
300
301
302
303
304 IF(ismstr==1.OR.ismstr==11.OR.ihbe<1)THEN
305 DO i=jft,jlt
306 hg1(i)=(vz1(i)-vz2(i)+vz3(i)-vz4(i))*off(i)
307 ENDDO
308 ELSE
309 DO i=jft,jlt
310 hg1(i)=vz1(i)*gama1(i)+vz2(i)*gama2(i)
311 . +vz3(i)*gama3(i)+vz4(i)*gama4(i)
312 ENDDO
313 ENDIF
314 DO i=jft,jlt
315 hour(i,3)=hour(i,3)+hg1(i)*hh2(i)
316 hour3a =hour(i,3)+hg1(i)*(h2l(i)+h2q(i)*abs(hg1(i)))
317 h31(i)=hour3a*gama1(i)
318 h32(i)=hour3a*gama2(i)
319 h33(i)=hour3a*gama3(i)
320
321 ehou(i) = ehou(i) + hour3a*hg1(i)
322 ENDDO
323
324
325
326
327 DO i=jft,jlt
328 hg1(i)=rx1(i)-rx2(i)+rx3(i)-rx4(i)
329 hg2(i)=ry1(i)-ry2(i)+ry3(i)-ry4(i)
330 ENDDO
331
332 DO i=jft,jlt
333 hour(i,4)=hg1(i)*(h3l(i)+h3q(i)*abs(hg1(i)))
334 hour(i,5)=hg2(i)*(h3l(i)+h3q(i)*abs(hg2(i)))
335 ehou(i) = ehou(i) +
336 . hour(i,4)*hg1(i) + hour(i,5)*hg2(i)
337 ehou(i) = dt1c(i) * ehou(i) * off(i)
338
339 b11(i)= hour(i,4)*off(i)
340 b12(i)=-hour(i,4)*off(i)
341 b13(i)= hour(i,4)*off(i)
342 b14(i)=-hour(i,4)*off(i)
343 b21(i)= hour(i,5)*off(i)
344 b22(i)=-hour(i,5)*off(i)
345 b23(i)= hour(i,5)*off(i)
346 b24(i)=-hour(i,5)*off(i)
347 ENDDO
348
349 DO i=jft,jlt
350 ehourt = ehourt + ehou(i)
351 ENDDO
352
353 ic=1
354 jst(ic)=jft
355 DO j=jft+1,jlt
356 IF (ipartc(j)/=ipartc(j-1)) THEN
357 ic=ic+1
358 jst(ic)=j
359 ENDIF
360 ENDDO
361
362 jst(ic+1)=jlt+1
363 IF(ic==1) THEN
364 mx = ipartc(jft)
365 DO i=jft,jlt
366 partsav(8,mx)=partsav(8,mx) + ehou(i)
367 ENDDO
368
369 ELSEIF(ic==2.AND.kfts>0)THEN
370 mx = ipartc(jft)
371 DO i=jft, kfts-1
372 partsav(8,mx)=partsav(8,mx) + ehou(i)
373 ENDDO
374 mx = ipartc(jlt)
375 DO i=kfts,jlt
376 partsav(8,mx)=partsav(8,mx) + ehou(i)
377 ENDDO
378
379 ELSE
380
381 DO ii=1,ic
382 mx=ipartc(jst(ii))
383 IF (jst(ii+1)-jst(ii)>15) THEN
384 DO j=jst(ii),jst(ii+1)-1
385 partsav(8,mx)=partsav(8,mx)+ehou(j)
386 ENDDO
387 ELSE
388 DO j=jst(ii),jst(ii+1)-1
389 partsav(8,mx)=partsav(8,mx)+ehou(j)
390 ENDDO
391 ENDIF
392 ENDDO
393 ENDIF
394
395
396 ehour = ehour + ehourt
397
398 DO i=jft,jlt
399 eani(nft+numels+i) = eani(nft+numels+i) + ehou(i)
400 ENDDO
401
402 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)