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 "scr14_c.inc"
63#include "scr16_c.inc"
64#include "scr17_c.inc"
65#include "sms_c.inc"
66
67
68
69 INTEGER IXC(NIXC,*),IPARTC(*), JFT, JLT,
70 . IHBE ,NFT ,ISMSTR,IGTYP,IGMAT,NEL
71
73 . thk(*), hour(nel,5), off(*),eani(*),partsav(npsav,*),
74 . px1(*), px2(*), py1(*), py2(*), sigy(*),dt1c(*),stir(*),
75 . ssp(mvsiz), rho(mvsiz),sti(mvsiz),z2(mvsiz), shf(mvsiz),
76 . thk0(mvsiz),viscmx(mvsiz), alpe(mvsiz), a11r(*)
77
79 . b11(mvsiz), b12(mvsiz), b13(mvsiz), b14(mvsiz),
80 . b21(mvsiz), b22(mvsiz), b23(mvsiz), b24(mvsiz),
81 . h11(mvsiz), h12(mvsiz), h13(mvsiz), h14(mvsiz),
82 . h21(mvsiz), h22(mvsiz), h23(mvsiz), h24(mvsiz),
83 . h31(mvsiz), h32(mvsiz), h33(mvsiz), h34(mvsiz),
84 . rx1(mvsiz), rx2(mvsiz), rx3(mvsiz), rx4(mvsiz),
85 . ry1(mvsiz), ry2(mvsiz), ry3(mvsiz), ry4(mvsiz),
86 . vx1(mvsiz), vx2(mvsiz), vx3(mvsiz), vx4(mvsiz),
87 . vy1(mvsiz), vy2(mvsiz), vy3(mvsiz), vy4(mvsiz),
88 . vz1(mvsiz), vz2(mvsiz), vz3(mvsiz), vz4(mvsiz),
89 .
area(mvsiz),vhx(mvsiz),
90 . h1(mvsiz), h2(mvsiz), h3(mvsiz), vhy(mvsiz),
91 . ym(mvsiz), nu(mvsiz), thk02(mvsiz),srh3(*)
92
93
94
95 INTEGER I, MX
96
98 . gama1(mvsiz), gama2(mvsiz), gama3(mvsiz), gama4(mvsiz),
99 . h3l(mvsiz), h3q(mvsiz),
100 . hg1(mvsiz), hg2(mvsiz) ,hp1(mvsiz), hp2(mvsiz) ,
101 . g(mvsiz) , b1(mvsiz), b2(mvsiz),
102 . a11(mvsiz) ,ehou(mvsiz),
103 . px1v, px2v, py1v, py2v,vv,
104 . shfpr3, fac, hour1(mvsiz), hour2(mvsiz), hour3(mvsiz),
105 . ehourt, fmax(mvsiz),
106 .
107 . scale(mvsiz)
108
109
110 ehourt = zero
111
112 IF(ismstr/=3.AND.ihbe/=0)THEN
113 DO i=jft,jlt
114 px1v = px1(i)*vhx(i)
115 px2v = px2(i)*vhx(i)
116 py1v = py1(i)*vhy(i)
117 py2v = py2(i)*vhy(i)
118 gama1(i)= off(i)*( one - px1v-py1v)
119 gama3(i)= off(i)*( one + px1v+py1v)
120 gama2(i)= off(i)*(-one - px2v-py2v)
121 gama4(i)= off(i)*(-one + px2v+py2v)
122 ENDDO
123 ELSE
124 DO i=jft,jlt
125 gama1(i)= off(i)
126 gama3(i)= off(i)
127 gama2(i)= -off(i)
128 gama4(i)= -off(i)
129 ENDDO
130 ENDIF
131
132 DO i=jft,jlt
133 hp1(i) =h1(i)
134 hp2(i) =h2(i)
135 h3(i) =srh3(i)
136 shfpr3 = shf(i)/(three*(one+nu(i)))
137 fac = fourth*rho(i)*thk02(i)
138 h3l(i) = h3(i)*zep072169*fac*
area(i)
139 h3q(i) = h3(i)*h3l(i)*hundred
140 h3l(i) = h3l(i)*ssp(i)
141
142 fac = ym(i)*thk0(i)*dt1c(i)*one_over_8
143 h1(i) = hp1(i) * fac
144 b1(i) = px1(i)*px1(i)+py1(i)*py1(i)
145 b2(i) = px2(i)*px2(i)+py2(i)*py2(i)
146 h2(i) = hp2(i) * fac*thk02(i)*shfpr3/(b1(i)+b2(i))
147 ENDDO
148
149
150
151 DO i=jft,jlt
152 IF(ixc(4,i)==ixc(5,i))THEN
153 h3q(i)=zero
154 h3l(i)=zero
155 h1(i)=zero
156 h2(i)=zero
157 ENDIF
158 ENDDO
159
160
161
162 DO i=jft,jlt
163 scale(i) = zero
164 ENDDO
165 IF(nodadt/=0.OR.idt1sh==1.OR.idtmins==2)THEN
166 DO i=jft,jlt
167 scale(i) = +
max(gama1(i)*gama1(i),gama2(i)*gama2(i),
168 . gama3(i)*gama3(i),gama4(i)*gama4(i)) *
169 . dt1c(i)*
max(h1(i),h2(i),h3l(i)) /
170 .
max(dt1c(i)*dt1c(i),em20)
171 sti(i)=sti(i) + scale(i)
172 ENDDO
173
174
175 IF(igtyp == 52 .OR.
176 . ((igtyp == 11 .OR. igtyp == 17 .OR. igtyp == 51)
177 . .AND. igmat > 0)) THEN
178 IF(nadmesh==0)THEN
179 DO i=jft,jlt
180 IF (off(i)==zero) THEN
181 sti(i) = zero
182 stir(i) = zero
183 ELSE
184 vv = viscmx(i) * viscmx(i) * alpe(i)
185 fac =
max(b1(i),b2(i))/(
area(i) * vv)
186 sti(i) = sti(i) + fac * thk0(i) * a11(i)
187 stir(i) = fac * a11r(i)*thk0(i)*thk02(i)*one_over_12 +
188 . fac * a11(i) *thk0(i)*
area(i) *one_over_9 +
189 . scale(i)*(thk02(i)* one_over_12 +
area(i)*one_over_9)
190 ENDIF
191 ENDDO
192 ELSE
193 DO i=jft,jlt
194 IF (off(i)==zero) THEN
195 sti(i) = zero
196 stir(i) = zero
197 ELSE
198 vv = viscmx(i) * viscmx(i) * alpe(i)
199 fac =
max(b1(i),b2(i))/(
area(i) * vv)
200 sti(i) = sti(i) + fac * thk0(i) * a11(i)
201 stir(i) = fac * a11r(i)*thk0(i)*thk02(i)* one_over_12 +
202 . scale(i)*thk02(i)* one_over_12
203 ENDIF
204 ENDDO
205 END IF
206 ELSE
207 IF(nadmesh==0)THEN
208 DO i=jft,jlt
209 IF (off(i)==zero) THEN
210 sti(i) = zero
211 stir(i) = zero
212 ELSE
213 vv = viscmx(i) * viscmx(i) * alpe(i)
214 sti(i) = sti(i) +
max(b1(i),b2(i))
215 . * thk0(i) * a11(i) / (
area(i) * vv)
216 stir(i) = sti(i) * (thk02(i)* one_over_12 +
area(i) * one_over_9)
217
218 ENDIF
219 ENDDO
220 ELSE
221 DO i=jft,jlt
222 IF (off(i)==zero) THEN
223 sti(i) = zero
224 stir(i) = zero
225 ELSE
226 vv = viscmx(i) * viscmx(i) * alpe(i)
227 sti(i) = sti(i) +
max(b1(i),b2(i))
228 . * thk0(i) * a11(i) / (
area(i) * vv)
229 stir(i) = sti(i) * thk02(i)* one_over_12
230 ENDIF
231 ENDDO
232 END IF
233 endif
234 ENDIF
235
236
237
238 DO i=jft,jlt
239 hg1(i)=vx1(i)*gama1(i)+vx2(i)*gama2(i)
240 . +vx3(i)*gama3(i)+vx4(i)*gama4(i)
241 hg2(i)=vy1(i)*gama1(i)+vy2(i)*gama2(i)
242 . +vy3(i)*gama3(i)+vy4(i)*gama4(i)
243 ENDDO
244
245 DO i=jft,jlt
246 hour(i,1)=hour(i,1)+hg1(i)*h1(i)
247 hour(i,2)=hour(i,2)+hg2(i)*h1(i)
248 fmax(i) = half*sigy(i)*thk0(i)*sqrt(
area(i))*hp1(i)
249 hour(i,1)=sign(
min( abs(hour(i,1)),fmax(i) ) , hour(i,1))
250 hour(i,2)=sign(
min( abs(hour(i,2)),fmax(i) ) , hour(i,2))
251 hour1(i)=hour(i,1)
252 h11(i)=hour1(i)*gama1(i)
253 h12(i)=hour1(i)*gama2(i)
254 h13(i)=hour1(i)*gama3(i)
255 h14(i)=hour1(i)*gama4(i)
256 hour2(i)=hour(i,2)
257 h21(i)=hour2(i)*gama1(i)
258 h22(i)=hour2(i)*gama2(i)
259 h23(i)=hour2(i)*gama3(i)
260 h24(i)=hour2(i)*gama4(i)
261 ehou(i) = hour1(i)*hg1(i) + hour2(i)*hg2(i)
262 ENDDO
263
264
265
266
267 DO i=jft,jlt
268 hg1(i)=vz1(i)*gama1(i)+vz2(i)*gama2(i)
269 . +vz3(i)*gama3(i)+vz4(i)*gama4(i)
270 ENDDO
271 DO i=jft,jlt
272 hour(i,3)=hour(i,3)+hg1(i)*h2(i)
273 fmax(i) = fourth*sigy(i)*thk0(i)*thk0(i)*hp2(i)
274 hour(i,3)=sign(
min( abs(hour(i,3)),fmax(i) ) , hour(i,3))
275 hour3(i) =hour(i,3)
276 h31(i)=hour3(i)*gama1(i)
277 h32(i)=hour3(i)*gama2(i)
278 h33(i)=hour3(i)*gama3(i)
279 h34(i)=hour3(i)*gama4(i)
280 ehou(i) = ehou(i) + hour3(i)*hg1(i)
281 ENDDO
282
283
284
285
286 DO i=jft,jlt
287 hg1(i)=rx1(i)-rx2(i)+rx3(i)-rx4(i)
288 hg2(i)=ry1(i)-ry2(i)+ry3(i)-ry4(i)
289 ENDDO
290
291 DO i=jft,jlt
292 hour(i,4)=hg1(i)*(h3l(i)+h3q(i)*abs(hg1(i)))
293 hour(i,5)=hg2(i)*(h3l(i)+h3q(i)*abs(hg2(i)))
294 ehou(i) = ehou(i) + hour(i,4)*hg1(i) + hour(i,5)*hg2(i)
295 ehou(i) = dt1c(i) * ehou(i)*off(i)
296 ehourt = ehourt + ehou(i)
297 b11(i)= hour(i,4)*off(i)
298 b12(i)=-hour(i,4)*off(i)
299 b13(i)= hour(i,4)*off(i)
300 b14(i)=-hour(i,4)*off(i)
301 b21(i)= hour(i,5)*off(i)
302 b22(i)=-hour(i,5)*off(i)
303 b23(i)= hour(i,5)*off(i)
304 b24(i)=-hour(i,5)*off(i)
305 ENDDO
306
307 DO i=jft,jlt
308 mx = ipartc(i)
309 partsav(8,mx)=partsav(8,mx) + ehou(i)
310 ENDDO
311
312
313 ehour = ehour + ehourt
314
315 DO i=jft,jlt
316 eani(nft+numels+i) = eani(nft+numels+i)+ehou(i)
317 ENDDO
318
319 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)