OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
schour3_1.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "param_c.inc"
#include "com04_c.inc"
#include "com08_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine schour3_1 (pm, rho, off, vx1, vx2, vx3, vx4, vx5, vx6, vx7, vx8, vy1, vy2, vy3, vy4, vy5, vy6, vy7, vy8, vz1, vz2, vz3, vz4, vz5, vz6, vz7, vz8, f11, f21, f31, f12, f22, f32, f13, f23, f33, f14, f24, f34, f15, f25, f35, f16, f26, f36, f17, f27, f37, f18, f28, f38, px1h1, px1h2, px1h3, px1h4, px2h1, px2h2, px2h3, px2h4, px3h1, px3h2, px3h3, px3h4, px4h1, px4h2, px4h3, px4h4, hgx1, hgy2, hgz1, hgz2, vol, mat, cxx, pid, geo, fhour, rx0, ry0, sx0, sy0, aj5, eint, eintm, vol0, sigy, sig0, mm, nu, defp, icp, nel, mtn, nlay)

Function/Subroutine Documentation

◆ schour3_1()

subroutine schour3_1 ( pm,
rho,
off,
vx1,
vx2,
vx3,
vx4,
vx5,
vx6,
vx7,
vx8,
vy1,
vy2,
vy3,
vy4,
vy5,
vy6,
vy7,
vy8,
vz1,
vz2,
vz3,
vz4,
vz5,
vz6,
vz7,
vz8,
f11,
f21,
f31,
f12,
f22,
f32,
f13,
f23,
f33,
f14,
f24,
f34,
f15,
f25,
f35,
f16,
f26,
f36,
f17,
f27,
f37,
f18,
f28,
f38,
px1h1,
px1h2,
px1h3,
px1h4,
px2h1,
px2h2,
px2h3,
px2h4,
px3h1,
px3h2,
px3h3,
px3h4,
px4h1,
px4h2,
px4h3,
px4h4,
hgx1,
hgy2,
hgz1,
hgz2,
vol,
integer, dimension(*) mat,
cxx,
integer, dimension(*) pid,
geo,
fhour,
rx0,
ry0,
sx0,
sy0,
aj5,
eint,
eintm,
vol0,
sigy,
sig0,
mm,
nu,
defp,
integer icp,
integer nel,
integer, intent(in) mtn,
integer, intent(in) nlay )

Definition at line 28 of file schour3_1.F.

53C-----------------------------------------------
54C I m p l i c i t T y p e s
55C-----------------------------------------------
56#include "implicit_f.inc"
57C-----------------------------------------------
58C G l o b a l P a r a m e t e r s
59C-----------------------------------------------
60#include "mvsiz_p.inc"
61C-----------------------------------------------
62C C o m m o n B l o c k s
63C-----------------------------------------------
64#include "param_c.inc"
65#include "com04_c.inc"
66#include "com08_c.inc"
67C-----------------------------------------------
68C D u m m y A r g u m e n t s
69C-----------------------------------------------
70 INTEGER, INTENT(IN) :: MTN,NLAY
71 INTEGER MAT(*),PID(*),ICP,NEL
72 my_real
73 . pm(npropm,*),geo(npropg,*), rho(*),off(*),
74 . vx1(*),vx2(*),vx3(*),vx4(*),vx5(*),vx6(*),vx7(*),vx8(*),
75 . vy1(*),vy2(*),vy3(*),vy4(*),vy5(*),vy6(*),vy7(*),vy8(*),
76 . vz1(*),vz2(*),vz3(*),vz4(*),vz5(*),vz6(*),vz7(*),vz8(*),
77 . f11(*),f21(*),f31(*),f12(*),f22(*),f32(*),
78 . f13(*),f23(*),f33(*),f14(*),f24(*),f34(*),
79 . f15(*),f25(*),f35(*),f16(*),f26(*),f36(*),
80 . f17(*),f27(*),f37(*),f18(*),f28(*),f38(*),
81 . px1h1(*), px1h2(*), px1h3(*),px1h4(*),
82 . px2h1(*), px2h2(*), px2h3(*),px2h4(*),
83 . px3h1(*), px3h2(*), px3h3(*),px3h4(*),
84 . px4h1(*), px4h2(*), px4h3(*),px4h4(*),
85 . hgx1(*), hgy2(*),hgz1(*), hgz2(*),
86 . sig0(nel,6),mm(mvsiz,2),
87 . vol(*),cxx(*),
88 . fhour(nel,12),rx0(*),ry0(*),sx0(*),sy0(*),aj5(*),
89 . eint(*),eintm(*),sigy(*) ,vol0(*),defp(*),nu(*)
90C-----------------------------------------------
91C L o c a l V a r i a b l e s
92C-----------------------------------------------
93 INTEGER I, MX, J,K
94 my_real
95 . caq(mvsiz), fcl(mvsiz), fcq(mvsiz),
96 . hgx3(mvsiz), hgx4(mvsiz),hgy3(mvsiz), hgy4(mvsiz),
97 . hgz3(mvsiz), hgz4(mvsiz),thk(mvsiz),thk_1(mvsiz),
98 . g11(mvsiz),g21(mvsiz),g31(mvsiz),g41(mvsiz),
99 . g51(mvsiz),g61(mvsiz),g71(mvsiz),g81(mvsiz),
100 . g12(mvsiz),g22(mvsiz),g32(mvsiz),g42(mvsiz),
101 . g52(mvsiz),g62(mvsiz),g72(mvsiz),g82(mvsiz),
102 . g13(mvsiz),g23(mvsiz),g33(mvsiz),g43(mvsiz),
103 . g53(mvsiz),g63(mvsiz),g73(mvsiz),g83(mvsiz),
104 . g14(mvsiz),g24(mvsiz),g34(mvsiz),g44(mvsiz),
105 . g54(mvsiz),g64(mvsiz),g74(mvsiz),g84(mvsiz),
106 . nfhx3(mvsiz),nfhy3(mvsiz),nfhz3(mvsiz),
107 . nfhx4(mvsiz),nfhy4(mvsiz),nfhz4(mvsiz),
108 . nfhx1(mvsiz),nfhz1(mvsiz),nfhy2(mvsiz),nfhz2(mvsiz),
109 . hhx3(mvsiz),hhy3(mvsiz),hxy3(mvsiz),e0(mvsiz),
110 . g_3dt(mvsiz),nu1,nu2,dh_xe1,dh_ze1 ,dh_yk2, dh_zk2 ,
111 . dh_xe3,dh_yk3,dh_zk3 ,dh_ze3,dh_xe4,dh_yk4,dh_z4 ,
112 . c11,c22,e,cc,sigy2,s1,s2,s3,svm0,sr1,sr2,sr3,sr4,
113 . sr5,rx1,rxy1,sy1,syx1,vol_3,vmm0,svmm,svm1,svm2,
114 . svmr,svm,svmrst,fac,facm,facb,f_z,coef,coefh,se,sk,dp
115C coef= (3*0.8)^2 , FZ=0.5*0.95^2
116 my_real
117 . a1(mvsiz), a2(mvsiz), a3(mvsiz),facp,f_max,
118 . ss1,ss2,ss3,faczm
119 DATA f_z/.45/,coef/5.76/,coefh/0.5/,f_max/0.998/
120C-----------------------------------------------
121C OPEN(UNIT=12,FILE='DEBHC.TMP',STATUS='UNKNOWN',FORM='FORMATTED')
122C ----: r->ksi; s->eta; t->zeta------------
123! IF(MTN==11 .OR. MTN==17 .OR. MTN==6
124! . .OR.MTN==46.OR.MTN==47.OR.MTN==37.OR.MTN==51) RETURN !TS: a tester au Starter.
125C
126 IF(invstr>=35)THEN
127 DO i=1,nel
128 caq(i)=fourth*off(i)*geo(13,pid(i))
129 ENDDO
130 ELSE
131 DO i=1,nel
132 caq(i)=fourth*off(i)*pm(4,mat(i))
133 ENDDO
134 ENDIF
135 DO i=1,nel
136 thk(i) =one_over_8*aj5(i)
137 thk_1(i) =one/thk(i)
138 e0(i)=three*(one-two*nu(i))*pm(32,mat(i))
139 g_3dt(i)=one_over_6*e0(i)*off(i)*dt1/(one+nu(i))
140 ENDDO
141C
142 DO i=1,nel
143 fcl(i)=caq(i)*rho(i)*vol(i)**four_over_3
144 fcl(i)=zep01666666667*fcl(i)*cxx(i)
145 ENDDO
146C
147 DO i=1,nel
148C alpha =1
149C 1 1 -1 -1 -1 -1 1 1
150 g11(i)= one_over_8-px1h1(i)
151 g21(i)= one_over_8-px2h1(i)
152 g31(i)=-one_over_8-px3h1(i)
153 g41(i)=-one_over_8-px4h1(i)
154 g51(i)=-one_over_8+px3h1(i)
155 g61(i)=-one_over_8+px4h1(i)
156 g71(i)= one_over_8+px1h1(i)
157 g81(i)= one_over_8+px2h1(i)
158 ENDDO
159 DO i=1,nel
160C alpha =2
161C 1 -1 -1 1 -1 1 1 -1
162 g12(i)= one_over_8-px1h2(i)
163 g22(i)=-one_over_8-px2h2(i)
164 g32(i)=-one_over_8-px3h2(i)
165 g42(i)= one_over_8-px4h2(i)
166 g52(i)=-one_over_8+px3h2(i)
167 g62(i)= one_over_8+px4h2(i)
168 g72(i)= one_over_8+px1h2(i)
169 g82(i)=-one_over_8+px2h2(i)
170 ENDDO
171 DO i=1,nel
172C alpha =3
173C 1 -1 1 -1 1 -1 1 -1
174 g13(i)= one_over_8 -px1h3(i)
175 g23(i)=-one_over_8 -px2h3(i)
176 g33(i)= one_over_8 -px3h3(i)
177 g43(i)=-one_over_8 -px4h3(i)
178 g53(i)= one_over_8 +px3h3(i)
179 g63(i)=-one_over_8 +px4h3(i)
180 g73(i)= one_over_8 +px1h3(i)
181 g83(i)=-one_over_8 +px2h3(i)
182 hgx3(i)=
183 & g13(i)*vx1(i)+g23(i)*vx2(i)+g33(i)*vx3(i)+g43(i)*vx4(i)
184 & +g53(i)*vx5(i)+g63(i)*vx6(i)+g73(i)*vx7(i)+g83(i)*vx8(i)
185 hgy3(i)=
186 & g13(i)*vy1(i)+g23(i)*vy2(i)+g33(i)*vy3(i)+g43(i)*vy4(i)
187 & +g53(i)*vy5(i)+g63(i)*vy6(i)+g73(i)*vy7(i)+g83(i)*vy8(i)
188 hgz3(i)=
189 & g13(i)*vz1(i)+g23(i)*vz2(i)+g33(i)*vz3(i)+g43(i)*vz4(i)
190 & +g53(i)*vz5(i)+g63(i)*vz6(i)+g73(i)*vz7(i)+g83(i)*vz8(i)
191 ENDDO
192C
193C 1 -1 1 -1 -1 1 -1 1
194C alpha =4
195C -1 1 -1 1 1 -1 1 -1
196 DO i=1,nel
197 g14(i)=-one_over_8-px1h4(i)
198 g24(i)= one_over_8-px2h4(i)
199 g34(i)=-one_over_8-px3h4(i)
200 g44(i)= one_over_8-px4h4(i)
201 g54(i)= one_over_8+px3h4(i)
202 g64(i)=-one_over_8+px4h4(i)
203 g74(i)= one_over_8+px1h4(i)
204 g84(i)=-one_over_8+px2h4(i)
205 hgx4(i)=
206 & g14(i)*vx1(i)+g24(i)*vx2(i)+g34(i)*vx3(i)+g44(i)*vx4(i)
207 & +g54(i)*vx5(i)+g64(i)*vx6(i)+g74(i)*vx7(i)+g84(i)*vx8(i)
208 hgy4(i)=
209 & g14(i)*vy1(i)+g24(i)*vy2(i)+g34(i)*vy3(i)+g44(i)*vy4(i)
210 & +g54(i)*vy5(i)+g64(i)*vy6(i)+g74(i)*vy7(i)+g84(i)*vy8(i)
211 hgz4(i)=
212 & g14(i)*vz1(i)+g24(i)*vz2(i)+g34(i)*vz3(i)+g44(i)*vz4(i)
213 & +g54(i)*vz5(i)+g64(i)*vz6(i)+g74(i)*vz7(i)+g84(i)*vz8(i)
214 ENDDO
215c IF (ICP==1) THEN
216c DO I=1,NEL
217c A1(I) =HALF*(ONE-NU(I))
218c A2(I) =HALF*(NU(I)-ONE)
219c A3(I) = ZERO
220c ENDDO
221c ELSEIF (ICP==2) THEN
222c DO I=1,NEL
223c FACP=ONE-DEFP(I)/(SIGY(I)/E0(I)+DEFP(I))
224c CC = FACP*(ONE+NU(I))
225c A1(I) =HALF*(ONE-NU(I)+CC)
226c A2(I) =HALF*(NU(I)-ONE+CC)
227c A3(I) = ZERO
228c ENDDO
229c ELSE
230c DO I=1,NEL
231c A1(I) =ONE
232c A2(I) =NU(I)
233c A3(I) = ONE
234c ENDDO
235c ENDIF
236 DO i=1,nel
237 a1(i) =one
238 a2(i) =nu(i)
239 a3(i) = one
240 ENDDO
241 !-------elstic increament, attention il y a un fac=1/3----
242 DO i=1,nel
243 c11 = two*g_3dt(i)/(one-nu(i))
244 cc = g_3dt(i)*thk_1(i)
245 dh_xe1 = cc*hgx1(i)
246 dh_yk2 = cc*hgy2(i)
247 !-------reduit 0.1->0.2---------
248 e = zep4*g_3dt(i)*(one +nu(i))
249 cc =e*thk_1(i)
250 dh_ze1 = cc*hgz1(i)
251 dh_zk2 = cc*hgz2(i)
252 !-------reduit 0.01----------
253 dh_z4 = 0.33*cc*hgz4(i)
254c DH_Z4 = EM01*CC*HGZ4(I)
255 cc = c11*sy0(i)
256 dh_xe3 = cc*hgx3(i)
257 dh_xe4 = cc*hgx4(i)
258 cc = c11*rx0(i)
259 dh_yk3 = cc*hgy3(i)
260 dh_yk4 = cc*hgy4(i)
261 dh_zk3 = g_3dt(i)*rx0(i)*hgz3(i)
262 dh_ze3 = g_3dt(i)*sy0(i)*hgz3(i)
263C
264C -----------FHOUR(1-11,I):XE1,ZE1,YK2,ZK2,XE3,YK3,ZK3,ZE3,XE4,YK4,Z4-------
265 fhour(i,1) = off(i)*fhour(i,1) + dh_xe1
266 fhour(i,2) = off(i)*fhour(i,2) + dh_ze1
267 fhour(i,3) = off(i)*fhour(i,3) + dh_yk2
268 fhour(i,4) = off(i)*fhour(i,4) + dh_zk2
269 fhour(i,5) = off(i)*fhour(i,5) + dh_xe3
270 fhour(i,6) = off(i)*fhour(i,6) + dh_yk3
271 fhour(i,7) = off(i)*fhour(i,7) + dh_zk3
272 fhour(i,8) = off(i)*fhour(i,8) + dh_ze3
273 IF (nlay>1) THEN
274 fhour(i,9) = off(i)*fhour(i,9) + dh_xe4
275 fhour(i,10) = off(i)*fhour(i,10) + dh_yk4
276 END IF
277 fhour(i,11) = off(i)*fhour(i,11) + dh_z4
278C------ SIGY(I) : now min(sigy(ip)) and =sigy(e_plas=0) at beginning in law2
279 IF (sigy(i)<zep9ep30) THEN
280 se =(one-nu(i))*fhour(i,5)
281 sk =(one-nu(i))*fhour(i,6)
282 sr1 =se*se+sk*sk
283 svm =-se*sk
284 se =nu(i)*fhour(i,5)
285c SE =NU(I)*FHOUR(I,5)-FHOUR(I,2)
286 sk =fhour(i,6)
287c SK =FHOUR(I,6)-FHOUR(I,4)
288 sr2 =se*se+sk*sk
289 svm = svm+se*sk
290 se =fhour(i,5)
291c SE =FHOUR(I,5)-FHOUR(I,2)
292 sk =nu(i)*fhour(i,6)
293c SK =NU(I)*FHOUR(I,6)-FHOUR(I,4)
294 sr3 =se*se+sk*sk
295 svm = svm+se*sk
296 sr4 =(fhour(i,8)+fhour(i,1))*a3(i)
297 sr5 =(fhour(i,7)+fhour(i,3))*a3(i)
298 svmr = half*(sr1+sr2+sr3)+3*(sr4*sr4+sr5*sr5)
299 svm0 =mm(i,2)
300 svmr = svmr + max(svm,-svm)
301 svm2 = sqrt(svm0+coef*svmr)
302c SVM2 = SQRT(MM(I,2)+COEF*(SVMR+MIN(SVM,-SVM)))
303!-----------memb :ELASTIC-PLASTIC yield criterion, FACM=(1-r_m)------------
304 IF (svm2>sigy(i)) THEN
305 IF (fhour(i,12)==zero) fhour(i,12) = sigy(i)
306 IF (svm2>fhour(i,12)) THEN
307 fac = (svm2-fhour(i,12))/svm2
308 ELSE
309 fac = (svm2-sigy(i))/svm2
310 END IF
311 facm = min(f_max,one-fac)
312 fhour(i,12) = fhour(i,12)+facm*(svm2-fhour(i,12))
313 ELSE
314 facm = zero
315 fhour(i,12) = svm2
316 ENDIF
317C------bending behaviours :
318 IF (nlay>1) THEN
319 se =fhour(i,9)*fhour(i,9)
320 sk =fhour(i,10)*fhour(i,10)
321 sr1 = (one +nu(i)*nu(i))*(se+sk)
322 sr2 = sr1 - two*nu(i)*(se+sk)
323 sr3 = abs((six*nu(i)-two)*fhour(i,9)*fhour(i,10))
324 svmm = f_z*(sr1+sr2+sr3)
325 faczm = svmr/max(em20,(svmm+svmr))
326 IF (mm(i,1)>sigy(i)*sigy(i)) THEN
327 IF (svmm>em20) THEN
328 fac = svmr/(svmm+svmr)
329 facb =one-min(sqrt(fac),one-facm)
330 ELSE
331 facb = facm
332 END IF
333 facb = min(facb,f_max)
334 ELSE
335 facb = facm
336 ENDIF
337 fhour(i,9) = fhour(i,9) - facb*dh_xe4
338 fhour(i,10) = fhour(i,10) - facb*dh_yk4
339 END IF !(NLAY>1) THEN
340c DP = THIRD*(ONE+NU(I))*(DH_XE3+DH_YK3)
341 dp = zero
342 IF (facm>zero) THEN
343 !-----------Membrane terms-----------
344 fhour(i,1) = fhour(i,1) - facm*dh_xe1*a3(i)
345 fhour(i,3) = fhour(i,3) - facm*dh_yk2*a3(i)
346 fhour(i,5) = fhour(i,5) - facm*(dh_xe3-dp)
347 fhour(i,6) = fhour(i,6) - facm*(dh_yk3-dp)
348 fhour(i,7) = fhour(i,7) - facm*dh_zk3*a3(i)
349 fhour(i,8) = fhour(i,8) - facm*dh_ze3*a3(i)
350C------only elsto-plastic in memb dominate case
351 IF (faczm >0.99) THEN
352 fhour(i,2) = fhour(i,2) - faczm*dh_ze1
353 fhour(i,4) = fhour(i,4) - faczm*dh_zk2
354 END IF
355 ENDIF
356 ENDIF
357 ENDDO
358C
359 DO i=1,nel
360 cc = thk_1(i)*vol(i)
361 nfhx1(i) = cc*(fhour(i,1)+fhour(i,8))
362 nfhz1(i) = cc*fhour(i,2)
363 nfhy2(i) = cc*(fhour(i,3)+fhour(i,7))
364 nfhz2(i) = cc*fhour(i,4)
365 rx1 = sx0(i)*sx0(i)/rx0(i)+rx0(i)
366 rxy1 = sx0(i)*sy0(i)/rx0(i)+ry0(i)
367 sy1 = ry0(i)*ry0(i)/sy0(i)+sy0(i)
368 hhx3(i) = rx1*rx0(i)
369 hhy3(i) = sy1*sy0(i)
370 hxy3(i) = rxy1*rx0(i)*sqrt(nu(i))
371 nfhz3(i) = vol(i)*(sy1*fhour(i,8)+rx1*fhour(i,7)
372 . +sy0(i)*fhour(i,1)+rx0(i)*fhour(i,3))
373 sy1 = sy1*a1(i)
374 rx1 = rx1*a1(i)
375 rxy1 = rxy1*a2(i)
376 syx1 = rxy1*rx0(i)/sy0(i)
377 vol_3 = third*vol(i)
378 nfhx3(i) = vol(i)*(sy1*fhour(i,5)+rxy1*fhour(i,6))
379 nfhy3(i) = vol(i)*(rx1*fhour(i,6)+syx1*fhour(i,5))
380 nfhx4(i) = vol_3*(sy1*fhour(i,9)+rxy1*fhour(i,10))
381 nfhy4(i) = vol_3*(rx1*fhour(i,10)+syx1*fhour(i,9))
382 nfhz4(i) =thk_1(i)*vol_3*fhour(i,11)
383 ENDDO
384 !--------visco-termes-------
385 DO i=1,nel
386 cc = fcl(i)*thk_1(i)/sqrt(1+nu(i))
387 nfhx1(i) = nfhx1(i)+cc*(thk_1(i)*hgx1(i)+sy0(i)*hgz3(i))
388 nfhy2(i) = nfhy2(i)+cc*(thk_1(i)*hgy2(i)+rx0(i)*hgz3(i))
389 nfhz3(i) = nfhz3(i)+cc*(thk(i)*(hhx3(i)+hhy3(i))*hgz3(i)
390 . +sy0(i)*hgx1(i)+rx0(i)*hgy2(i))
391 cc = fifteen*fcl(i)*thk_1(i)*thk_1(i)
392 nfhz1(i) = nfhz1(i)+cc*hgz1(i)
393 nfhz2(i) = nfhz2(i)+cc*hgz2(i)
394 nfhz4(i) = nfhz4(i)+cc*hgz4(i)
395 cc = fcl(i)/sqrt(1-nu(i)*nu(i))
396 nfhx3(i) = nfhx3(i)+cc*(hhx3(i)*hgx3(i)+hxy3(i)*hgy3(i))
397 nfhy3(i) = nfhy3(i)+cc*(hhy3(i)*hgy3(i)+hxy3(i)*hgx3(i))
398 cc = third*cc
399 nfhx4(i) = nfhx4(i)+cc*(hhx3(i)*hgx4(i)+hxy3(i)*hgy4(i))
400 nfhy4(i) = nfhy4(i)+cc*(hhy3(i)*hgy4(i)+hxy3(i)*hgx4(i))
401 ENDDO
402C
403 DO i=1,nel
404 f11(i) =f11(i)-g11(i)*nfhx1(i)
405 . -g13(i)*nfhx3(i)-g14(i)*nfhx4(i)
406 f12(i) =f12(i)-g21(i)*nfhx1(i)
407 . -g23(i)*nfhx3(i)-g24(i)*nfhx4(i)
408 f13(i) =f13(i)-g31(i)*nfhx1(i)
409 . -g33(i)*nfhx3(i)-g34(i)*nfhx4(i)
410 f14(i) =f14(i)-g41(i)*nfhx1(i)
411 . -g43(i)*nfhx3(i)-g44(i)*nfhx4(i)
412 f15(i) =f15(i)-g51(i)*nfhx1(i)
413 . -g53(i)*nfhx3(i)-g54(i)*nfhx4(i)
414 f16(i) =f16(i)-g61(i)*nfhx1(i)
415 . -g63(i)*nfhx3(i)-g64(i)*nfhx4(i)
416 f17(i) =f17(i)-g71(i)*nfhx1(i)
417 . -g73(i)*nfhx3(i)-g74(i)*nfhx4(i)
418 f18(i) =f18(i)-g81(i)*nfhx1(i)
419 . -g83(i)*nfhx3(i)-g84(i)*nfhx4(i)
420C
421 f21(i) =f21(i)-g12(i)*nfhy2(i)
422 . -g13(i)*nfhy3(i)-g14(i)*nfhy4(i)
423 f22(i) =f22(i)-g22(i)*nfhy2(i)
424 . -g23(i)*nfhy3(i)-g24(i)*nfhy4(i)
425 f23(i) =f23(i)-g32(i)*nfhy2(i)
426 . -g33(i)*nfhy3(i)-g34(i)*nfhy4(i)
427 f24(i) =f24(i)-g42(i)*nfhy2(i)
428 . -g43(i)*nfhy3(i)-g44(i)*nfhy4(i)
429 f25(i) =f25(i)-g52(i)*nfhy2(i)
430 . -g53(i)*nfhy3(i)-g54(i)*nfhy4(i)
431 f26(i) =f26(i)-g62(i)*nfhy2(i)
432 . -g63(i)*nfhy3(i)-g64(i)*nfhy4(i)
433 f27(i) =f27(i)-g72(i)*nfhy2(i)
434 . -g73(i)*nfhy3(i)-g74(i)*nfhy4(i)
435 f28(i) =f28(i)-g82(i)*nfhy2(i)
436 . -g83(i)*nfhy3(i)-g84(i)*nfhy4(i)
437C
438 f31(i) =f31(i)-g11(i)*nfhz1(i)-g12(i)*nfhz2(i)
439 . -g13(i)*nfhz3(i)-g14(i)*nfhz4(i)
440 f32(i) =f32(i)-g21(i)*nfhz1(i)-g22(i)*nfhz2(i)
441 . -g23(i)*nfhz3(i)-g24(i)*nfhz4(i)
442 f33(i) =f33(i)-g31(i)*nfhz1(i)-g32(i)*nfhz2(i)
443 . -g33(i)*nfhz3(i)-g34(i)*nfhz4(i)
444 f34(i) =f34(i)-g41(i)*nfhz1(i)-g42(i)*nfhz2(i)
445 . -g43(i)*nfhz3(i)-g44(i)*nfhz4(i)
446 f35(i) =f35(i)-g51(i)*nfhz1(i)-g52(i)*nfhz2(i)
447 . -g53(i)*nfhz3(i)-g54(i)*nfhz4(i)
448 f36(i) =f36(i)-g61(i)*nfhz1(i)-g62(i)*nfhz2(i)
449 . -g63(i)*nfhz3(i)-g64(i)*nfhz4(i)
450 f37(i) =f37(i)-g71(i)*nfhz1(i)-g72(i)*nfhz2(i)
451 . -g73(i)*nfhz3(i)-g74(i)*nfhz4(i)
452 f38(i) =f38(i)-g81(i)*nfhz1(i)-g82(i)*nfhz2(i)
453 . -g83(i)*nfhz3(i)-g84(i)*nfhz4(i)
454 ENDDO
455 !----hourglass energy is included in internal energy------
456 DO i=1,nel
457 eint(i)= eint(i)+dt1*(
458 . nfhx1(i)*hgx1(i) + nfhz1(i)*hgz1(i) +
459 . nfhy2(i)*hgy2(i) + nfhz2(i)*hgz2(i) +
460 . nfhz3(i)*hgz3(i) + nfhz4(i)*hgz4(i) +
461 . nfhx3(i)*hgx3(i) + nfhx4(i)*hgx4(i) +
462 . nfhy3(i)*hgy3(i) + nfhy4(i)*hgy4(i) )
463 . /max(em20,vol0(i)) +eintm(i)
464 ENDDO
465C
466 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21