46
47
48
49#include "implicit_f.inc"
50
51
52
53#include "mvsiz_p.inc"
54
55
56
57 INTEGER G_PLA,G_EPSD,NEL
58
60 . sig(nel,6),
61 . px1(*), px2(*), px3(*), px4(*),
62 . py1(*), py2(*), py3(*), py4(*),
63 . pz1(*), pz2(*), pz3(*), pz4(*),
64 . px5(*), px6(*), px7(*), px8(*),
65 . py5(*), py6(*), py7(*), py8(*),
66 . pz5(*), pz6(*), pz7(*), pz8(*),
67 . f11(*),f21(*),f31(*),f12(*),f22(*),f32(*),
68 . f13(*),f23(*),f33(*),f14(*),f24(*),f34(*),
69 . f15(*),f25(*),f35(*),f16(*),f26(*),f36(*),
70 . f17(*),f27(*),f37(*),f18(*),f28(*),f38(*),
71 . vol(*),qvis(*),
72 . px1h1(*), px1h2(*), px2h1(*), px2h2(*),
73 . px3h1(*), px3h2(*), px4h1(*), px4h2(*),
74 . rx0(*), ry0(*), sx0(*), sy0(*),
75 . eint(*),rho(*),q(*),eplasm(*),eplas(*),
76 . sigm(nel,6),eintm(*),einto(*),rhom(*),qm(*),epsd(*),epsdm(*),
77 . nu(*),zi,wi,volg(*),mm(mvsiz,2),off(*),
78 + vol0(*),vol0g(*)
79 my_real,
DIMENSION(MVSIZ,6),
INTENT(INOUT) :: svis
80 INTEGER, INTENT(IN) :: G_WPLA_FLAG
81 my_real,
DIMENSION(NEL*G_WPLA_FLAG),
INTENT(INOUT) :: g_wpla
82 my_real,
DIMENSION(NEL*G_WPLA_FLAG),
INTENT(IN) :: l_wpla
83
84
85
86 INTEGER I, J
87
89 . s1(mvsiz), s2(mvsiz), s3(mvsiz),
90 . s4(mvsiz), s5(mvsiz), s6(mvsiz),
91 . rx1(mvsiz), ry1(mvsiz), sx1(mvsiz), sy1(mvsiz),
92 . ss1,ss2,ss3,svm0,fint,fint1,fint2,
93 . xh1i,xh2i,yh1i,yh2i,nu1(mvsiz),fac(mvsiz)
94
95 DO i=1,nel
96 s1(i)=(sig(i,1)+svis(i,1)-qvis(i))*vol(i)
97 s2(i)=(sig(i,2)+svis(i,2)-qvis(i))*vol(i)
98 s3(i)=(sig(i,3)+svis(i,3)-qvis(i))*vol(i)
99 s4(i)=(sig(i,4)+svis(i,4))*vol(i)
100 s5(i)=(sig(i,5)+svis(i,5))*vol(i)
101 s6(i)=(sig(i,6)+svis(i,6))*vol(i)
102 ENDDO
103 DO i=1,nel
104 fint=s1(i)*px1(i)+s4(i)*py1(i)+s6(i)*pz1(i)
105 f11(i)=f11(i)-fint
106 f17(i)=f17(i)+fint
107 fint=s2(i)*py1(i)+s4(i)*px1(i)+s5(i)*pz1(i)
108 f21(i)=f21(i)-fint
109 f27(i)=f27(i)+fint
110 fint=s3(i)*pz1(i)+s6(i)*px1(i)+s5(i)*py1(i)
111 f31(i)=f31(i)-fint
112 f37(i)=f37(i)+fint
113
114 fint=s1(i)*px2(i)+s4(i)*py2(i)+s6(i)*pz2(i)
115 f12(i)=f12(i)-fint
116 f18(i)=f18(i)+fint
117 fint=s2(i)*py2(i)+s4(i)*px2(i)+s5(i)*pz2(i)
118 f22(i)=f22(i)-fint
119 f28(i)=f28(i)+fint
120 fint=s3(i)*pz2(i)+s6(i)*px2(i)+s5(i)*py2(i)
121 f32(i)=f32(i)-fint
122 f38(i)=f38(i)+fint
123
124 fint=s1(i)*px3(i)+s4(i)*py3(i)+s6(i)*pz3(i)
125 f13(i)=f13(i)-fint
126 f15(i)=f15(i)+fint
127 fint=s2(i)*py3(i)+s4(i)*px3(i)+s5(i)*pz3(i)
128 f23(i)=f23(i)-fint
129 f25(i)=f25(i)+fint
130 fint=s3(i)*pz3(i)+s6(i)*px3(i)+s5(i)*py3(i)
131 f33(i)=f33(i)-fint
132 f35(i)=f35(i)+fint
133
134 fint=s1(i)*px4(i)+s4(i)*py4(i)+s6(i)*pz4(i)
135 f14(i)=f14(i)-fint
136 f16(i)=f16(i)+fint
137 fint=s2(i)*py4(i)+s4(i)*px4(i)+s5(i)*pz4(i)
138 f24(i)=f24(i)-fint
139 f26(i)=f26(i)+fint
140 fint=s3(i)*pz4(i)+s6(i)*px4(i)+s5(i)*py4(i)
141 f34(i)=f34(i)-fint
142 f36(i)=f36(i)+fint
143 ENDDO
144
145 DO i=1,nel
146 ry1(i)= zi*ry0(i)
147 rx1(i)= zi*rx0(i)
148 sy1(i)= zi*sy0(i)
149 sx1(i)= zi*sx0(i)
150 nu1(i) = nu(i)/(one -nu(i))
151 ENDDO
152
153 DO i=1,nel
154
155 xh1i =-ry1(i)*px1h1(i)
156 xh2i = sy1(i)*px1h2(i)
157 yh1i = rx1(i)*px1h1(i)
158 yh2i =-sx1(i)*px1h2(i)
159 fint=s1(i)*(xh1i+xh2i)-nu(i)*xh1i*s2(i)
160 fint=fint-(nu1(i)*xh2i+nu(i)*xh1i)*s3(i)+s4(i)*(yh1i+yh2i)
161 f11(i)=f11(i)+fint
162 f17(i)=f17(i)-fint
163 fint=s2(i)*(yh1i+yh2i)-nu(i)*yh2i*s1(i)
164 fint=fint-(nu(i)*yh2i+nu1(i)*yh1i)*s3(i)+s4(i)*(xh1i+xh2i)
165 f21(i)=f21(i)+fint
166 f27(i)=f27(i)-fint
167 fint= s6(i)*xh1i+s5(i)*yh2i
168 f31(i)=f31(i)+fint
169 f37(i)=f37(i)-fint
170
171 xh1i =-ry1(i)*px2h1(i)
172 xh2i = sy1(i)*px2h2(i)
173 yh1i = rx1(i)*px2h1(i)
174 yh2i =-sx1(i)*px2h2(i)
175 fint=s1(i)*(xh1i+xh2i)-nu(i)*xh1i*s2(i)
176 fint=fint-(nu(i)*xh1i+nu1(i)*xh2i)*s3(i)+s4(i)*(yh1i+yh2i)
177 f12(i)=f12(i)+fint
178 f18(i)=f18(i)-fint
179 fint=s2(i)*(yh1i+yh2i)-nu(i)*yh2i*s1(i)
180 fint=fint-(nu(i)*yh2i+nu1(i)*yh1i)*s3(i)+s4(i)*(xh1i+xh2i)
181 f22(i)=f22(i)+fint
182 f28(i)=f28(i)-fint
183 fint= s6(i)*xh1i+s5(i)*yh2i
184 f32(i)=f32(i)+fint
185 f38(i)=f38(i)-fint
186
187 xh1i =-ry1(i)*px3h1(i)
188 xh2i = sy1(i)*px3h2(i)
189 yh1i = rx1(i)*px3h1(i)
190 yh2i =-sx1(i)*px3h2(i)
191 fint=s1(i)*(xh1i+xh2i)-nu(i)*xh1i*s2(i)
192 fint=fint-(nu(i)*xh1i+nu1(i)*xh2i)*s3(i)+s4(i)*(yh1i+yh2i)
193 f13(i)=f13(i)+fint
194 f15(i)=f15(i)-fint
195 fint=s2(i)*(yh1i+yh2i)-nu(i)*yh2i*s1(i)
196 fint=fint-(nu(i)*yh2i+nu1(i)*yh1i)*s3(i)+s4(i)*(xh1i+xh2i)
197 f23(i)=f23(i)+fint
198 f25(i)=f25(i)-fint
199 fint= s6(i)*xh1i+s5(i)*yh2i
200 f33(i)=f33(i)+fint
201 f35(i)=f35(i)-fint
202
203 xh1i =-ry1(i)*px4h1(i)
204 xh2i = sy1(i)*px4h2(i)
205 yh1i = rx1(i)*px4h1(i)
206 yh2i =-sx1(i)*px4h2(i)
207 fint=s1(i)*(xh1i+xh2i)-nu(i)*xh1i*s2(i)
208 fint=fint-(nu(i)*xh1i+nu1(i)*xh2i)*s3(i)+s4(i)*(yh1i+yh2i)
209 f14(i)=f14(i)+fint
210 f16(i)=f16(i)-fint
211 fint=s2(i)*(yh1i+yh2i)-nu(i)*yh2i*s1(i)
212 fint=fint-(nu(i)*yh2i+nu1(i)*yh1i)*s3(i)+s4(i)*(xh1i+xh2i)
213 f24(i)=f24(i)+fint
214 f26(i)=f26(i)-fint
215 fint= s6(i)*xh1i+s5(i)*yh2i
216 f34(i)=f34(i)+fint
217 f36(i)=f36(i)-fint
218 ENDDO
219
220
221
222 DO i=1,nel
223 xh1i =-one_over_8*ry1(i)
224 xh2i = one_over_8*sy1(i)
225 yh1i = one_over_8*rx1(i)
226 yh2i =-one_over_8*sx1(i)
227
228 fint1=(s1(i)-nu(i)*s2(i)-nu(i)*s3(i))*xh1i+s4(i)*yh1i
229 fint2=(s1(i)-nu1(i)*s3(i))*xh2i+s4(i)*yh2i
230
231 f11(i)=f11(i)-fint1-fint2
232 f12(i)=f12(i)-fint1+fint2
233 f13(i)=f13(i)+fint1+fint2
234 f14(i)=f14(i)+fint1-fint2
235 f15(i)=f15(i)+fint1+fint2
236 f16(i)=f16(i)+fint1-fint2
237 f17(i)=f17(i)-fint1-fint2
238 f18(i)=f18(i)-fint1+fint2
239
240 fint1=(s2(i)-nu1(i)*s3(i))*yh1i+s4(i)*xh1i
241 fint2=(s2(i)-nu(i)*s1(i)-nu(i)*s3(i))*yh2i+s4(i)*xh2i
242
243 f21(i)=f21(i)-fint1-fint2
244 f22(i)=f22(i)-fint1+fint2
245 f23(i)=f23(i)+fint1+fint2
246 f24(i)=f24(i)+fint1-fint2
247 f25(i)=f25(i)+fint1+fint2
248 f26(i)=f26(i)+fint1-fint2
249 f27(i)=f27(i)-fint1-fint2
250 f28(i)=f28(i)-fint1+fint2
251
252 fint1=s6(i)*xh1i
253 fint2=s5(i)*yh2i
254
255 f31(i)=f31(i)-fint1-fint2
256 f32(i)=f32(i)-fint1+fint2
257 f33(i)=f33(i)+fint1+fint2
258 f34(i)=f34(i)+fint1-fint2
259 f35(i)=f35(i)+fint1+fint2
260 f36(i)=f36(i)+fint1-fint2
261 f37(i)=f37(i)-fint1-fint2
262 f38(i)=f38(i)-fint1+fint2
263 ENDDO
264
265
266
267 DO i=1,nel
268 ss1 =sig(i,1)-sig(i,2)
269 ss2 =sig(i,2)-sig(i,3)
270 ss3 =sig(i,1)-sig(i,3)
271 svm0 =(ss1*ss1+ss2*ss2+ss3*ss3)*half + three*(sig(i,4)*sig(i,4)
272 . + sig(i,5)*sig(i,5)+sig(i,6)*sig(i,6))
273 mm(i,1) =
max(svm0,mm(i,1))
274 mm(i,2) =
min(svm0,mm(i,2))
275 ENDDO
276
277
278 DO i=1,nel
279 fac(i) = off(i)*vol(i)/volg(i)
280 sigm(i,1) = sigm(i,1) + fac(i) * sig(i,1)
281 sigm(i,2) = sigm(i,2) + fac(i) * sig(i,2)
282 sigm(i,3) = sigm(i,3) + fac(i) * sig(i,3)
283 sigm(i,4) = sigm(i,4) + fac(i) * sig(i,4)
284 sigm(i,5) = sigm(i,5) + fac(i) * sig(i,5)
285 sigm(i,6) = sigm(i,6) + fac(i) * sig(i,6)
286 rhom(i) = rhom(i) + fac(i) * rho(i)
287 eintm(i) = eintm(i)+off(i)*vol0(i)/vol0g(i)*(eint(i)-einto(i))
288 IF (g_wpla_flag > 0) g_wpla(i) = g_wpla(i) + l_wpla(i)
289 qm(i) = qm(i) + fac(i) * q(i)
290 ENDDO
291
292 IF (g_pla > 0) THEN
293 DO i=1,nel
294 eplasm(i) = eplasm(i)
295 ENDDO
296 ENDIF
297 IF (g_epsd > 0) THEN
298 DO i=1,nel
299 epsdm(i) = epsdm(i) + fac
300 ENDDO
301 ENDIF
302
303 RETURN