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