52
53
54
55#include "implicit_f.inc"
56
57
58
59#include "mvsiz_p.inc"
60
61
62
63#include "com08_c.inc"
64#include "impl1_c.inc"
65
66
67
68 INTEGER, INTENT(IN) :: NEL
69 INTEGER, INTENT(IN) :: ISMSTR
70 INTEGER ICP,I_SH,IDEGE(*)
71
73 . dxx(*), dxy(*), dxz(*),
74 . dyx(*), dyy(*), dyz(*),
75 . dzx(*), dzy(*), dzz(*), d4(*), d5(*), d6(*),
76 . wxx(*), wyy(*), wzz(*), offs(*),dsv(*),
77 . volo(*),off(*),eint(*),fac(*),sdv(*),
78 . dn_x(mvsiz,8),dn_y(mvsiz,8),dn_z(mvsiz,8),
79 . dn_r(8), dn_s(8), dn_t(8), invj(9,mvsiz),
80 . ulx1(mvsiz), ulx2(mvsiz), ulx3(mvsiz), ulx4(mvsiz),
81 . ulx5(mvsiz), ulx6(mvsiz), ulx7(mvsiz), ulx8(mvsiz),
82 . uly1(mvsiz), uly2(mvsiz), uly3(mvsiz), uly4(mvsiz),
83 . uly5(mvsiz), uly6(mvsiz), uly7(mvsiz), uly8(mvsiz),
84 . ulz1(mvsiz), ulz2(mvsiz), ulz3(mvsiz), ulz4(mvsiz),
85 . ulz5(mvsiz), ulz6(mvsiz), ulz7(mvsiz), ulz8(mvsiz),
86 . x1(*),x2(*),x3(*),x4(*),x5(*),x6(*)
87 . y1(*),y2(*),y3(*),y4(*),y5(*),y6(*),y7(*),y8(*),
88 . z1(*),z2(*),z3(*),z4(*),z5(*),z6(*),z7(*),z8(*),
89 . a11(mvsiz), a12(mvsiz), a13(mvsiz),
90 . a21(mvsiz), a22(mvsiz), a23(mvsiz),
91 . a31(mvsiz), a32(mvsiz), a33(mvsiz),
92 . bb(6
93 . ksi,eta,zeta
94
95
96
97 INTEGER I,J,JJ,L,M,LM
98
100 . dt1d2,dvc(mvsiz),tol,dt1d,udt1,
101 . b(6,24),um,up,c(6,6),
102 . a1a1,a1a2,a1a3,a2a1,a2a2,a2a3,
103 . a1b1,a1b2,a1b3,a2b1,a2b2,a2b3,
104 . b1a1,b1a2,b1a3,b2a1,b2a2,b2a3,
105 . b1b1,b1b2,b1b3,b2b1,b2b2,b2b3,
106 . c1a1,c1a2,c1a3,c2a1,c2a2,c2a3,
107 . c1b1,c1b2,c1b3,c2b1,c2b2,c2b3,
108 . ca11,ca21,ca31,ca12,ca22,ca32,ca13,ca23,ca33,
109 . d11,d22,d33,d12,d13,d23
110 double precision
111 . dvp,dv,dv1
112
113 tol = one-em20
114
115
116 DO i=1,nel
117 DO j=1,8
118 jj = 3*(j-1)
119 b(1,jj+1) = a11(i)*dn_t(j)
120 b(1,jj+2) = a12(i)*dn_t(j)
121 b(1,jj+3) = a13(i)*dn_t(j)
122 b(2,jj+1) = a21(i)*dn_r(j)
123 b(2,jj+2) = a22(i)*dn_r(j)
124 b(2,jj+3) = a23(i)*dn_r(j)
125 b(3,jj+1) = a31(i)*dn_s(j)
126 b(3,jj+2) = a32(i)*dn_s(j)
127 b(3,jj+3) = a33(i)*dn_s(j)
128 ENDDO
129
130
131 a1a1 = (x2(i)+x3(i)-x1(i)-x4(i))*fourth
132 a1a2 = (y2(i)+y3(i)-y1(i)-y4(i))*fourth
133 a1a3 = (z2(i)+z3(i)-z1(i)-z4(i))*fourth
134 a2a1 = (x4(i)+x3(i)-x1(i)-x2(i))*fourth
135 a2a2 = (y4(i)+y3(i)-y1(i)-y2(i))*fourth
136 a2a3 = (z4(i)+z3(i)-z1(i)-z2(i))*fourth
137 a1b1 = (x6(i)+x7(i)-x5(i)-x8(i))*fourth
138 a1b2 = (y6(i)+y7(i)-y5(i)-y8(i))*fourth
139 a1b3 = (z6(i)+z7(i)-z5(i)-z8(i))*fourth
140 a2b1 = (x8(i)+x7(i)-x5(i)-x6(i))*fourth
141 a2b2 = (y8(i)+y7(i)-y5(i)-y6(i))*fourth
142 a2b3 = (z8(i)+z7(i)-z5(i)-z6(i))*fourth
143
144 um = (one-zeta)/eight
145 up = (one+zeta)/eight
146 b(4,1) = -um*(a1a1+a2a1)
147 b(4,2) = -um*(a1a2+a2a2)
148 b(4,3) = -um*(a1a3+a2a3)
149 b(4,4) = um*(a2a1-a1a1)
150 b(4,5) = um*(a2a2-a1a2)
151 b(4,6) = um*(a2a3-a1a3)
152 b(4,7) = um*(a2a1+a1a1)
153 b(4,8) = um*(a2a2+a1a2)
154 b(4,9) = um*(a2a3+a1a3)
155 b(4,10)= um*(a1a1-a2a1)
156 b(4,11)= um*(a1a2-a2a2)
157 b(4,12)= um*(a1a3-a2a3)
158 b(4,13) =-up*(a1b1+a2b1)
159 b(4,14) =-up*(a1b2+a2b2)
160 b(4,15) =-up*(a1b3+a2b3)
161 b(4,16) = up*(a2b1-a1b1)
162 b(4,17) = up*(a2b2-a1b2)
163 b(4,18) = up*(a2b3-a1b3)
164 b(4,19) = up*(a2b1+a1b1)
165 b(4,20) = up*(a2b2+a1b2)
166 b(4,21) = up*(a2b3+a1b3)
167 b(4,22) = up*(a1b1-a2b1)
168 b(4,23) = up*(a1b2-a2b2)
169 b(4,24) = up*(a1b3-a2b3)
170
171
172 b1a1 = (x2(i)+x6(i)-x1(i)-x5(i))*fourth
173 b1a2 = (y2(i)+y6(i)-y1(i)-y5(i))*fourth
174 b1a3 = (z2(i)+z6(i)-z1(i)-z5(i))*fourth
175 b2a1 = (x5(i)+x6(i)-x1(i)-x2(i))*fourth
176 b2a2 = (y5(i)+y6(i)-y1(i)-y2(i))*fourth
177 b2a3 = (z5(i)+z6(i)-z1(i)-z2(i))*fourth
178 b1b1 = (x3(i)+x7(i)-x4(i)-x8(i))*fourth
179 b1b2 = (y3(i)+y7(i)-y4(i)-y8(i))*fourth
180 b1b3 = (z3(i)+z7(i)-z4(i)-z8(i))*fourth
181 b2b1 = (x7(i)+x8(i)-x3(i)-x4(i))*fourth
182 b2b2 = (y7(i)+y8(i)-y3(i)-y4(i))*fourth
183 b2b3 = (z7(i)+z8(i)-z3(i)-z4(i))*fourth
184
185 um = (one-eta)/eight
186 up = (one+eta)/eight
187 b(5,1) = -um*(b1a1+b2a1)
188 b(5,2) = -um*(b1a2+b2a2)
189 b(5,3) = -um*(b1a3+b2a3)
190 b(5,4) = um*(b2a1-b1a1)
191 b(5,5) = um*(b2a2-b1a2)
192 b(5,6) = um*(b2a3-b1a3)
193 b(5,7) = up*(b2b1-b1b1)
194 b(5,8) = up*(b2b2-b1b2)
195 b(5,9) = up*(b2b3-b1b3)
196 b(5,10)= -up*(b1b1+b2b1)
197 b(5,11)= -up*(b1b2+b2b2)
198 b(5,12)= -up*(b1b3+b2b3)
199 b(5,13)= um*(b1a1-b2a1)
200 b(5,14)= um*(b1a2-b2a2)
201 b(5,15)= um*(b1a3-b2a3)
202 b(5,16)= um*(b1a1+b2a1)
203 b(5,17)= um*(b1a2+b2a2)
204 b(5,18)= um*(b1a3+b2a3)
205 b(5,19)= up*(b1b1+b2b1)
206 b(5,20)= up*(b1b2+b2b2)
207 b(5,21)= up*(b1b3+b2b3)
208 b(5,22)= up*(b1b1-b2b1)
209 b(5,23)= up*(b1b2-b2b2)
210 b(5,24)= up*(b1b3-b2b3)
211
212
213 c1a1 = (x4(i)+x8(i)-x1(i)-x5(i))*fourth
214 c1a2 = (y4(i)+y8(i)-y1(i)-y5(i))*fourth
215 c1a3 = (z4(i)+z8(i)-z1(i)-z5(i))*fourth
216 c2a1 = (x5(i)+x8(i)-x1(i)-x4(i))*fourth
217 c2a2 = (y5(i)+y8(i)-y1(i)-y4(i))*fourth
218 c2a3 = (z5(i)+z8(i)-z1(i)-z4(i))*fourth
219 c1b1 = (x3(i)+x7(i)-x2(i)-x6(i))*fourth
220 c1b2 = (y3(i)+y7(i)-y2(i)-y6(i))*fourth
221 c1b3 = (z3(i)+z7(i)-z2(i)-z6(i))*fourth
222 c2b1 = (x6(i)+x7(i)-x2(i)-x3(i))*fourth
223 c2b2 = (y6(i)+y7(i)-y2(i)-y3(i))*fourth
224 c2b3 = (z6(i)+z7(i)-z2(i)-z3(i))*fourth
225
226 um = (one-ksi)/eight
227 up = (one+ksi)/eight
228 b(6,1) = -um*(c1a1+c2a1)
229 b(6,2) = -um*(c1a2+c2a2)
230 b(6,3) = -um*(c1a3+c2a3)
231 b(6,4) = -up*(c1b1+c2b1)
232 b(6,5) = -up*(c1b2+c2b2)
233 b(6,6) = -up*(c1b3+c2b3)
234 b(6,7) = up*(c2b1-c1b1)
235 b(6,8) = up*(c2b2-c1b2)
236 b(6,9) = up*(c2b3-c1b3)
237 b(6,10)= um*(c2a1-c1a1)
238 b(6,11)= um*(c2a2-c1a2)
239 b(6,12)= um*(c2a3-c1a3)
240 b(6,13)= um*(c1a1-c2a1)
241 b(6,14)= um*(c1a2-c2a2)
242 b(6,15)= um*(c1a3-c2a3)
243 b(6,16)= up*(c1b1-c2b1)
244 b(6,17)= up*(c1b2-c2b2)
245 b(6,18)= up*(c1b3-c2b3)
246 b(6,19)= up*(c1b1+c2b1)
247 b(6,20)= up*(c1b2+c2b2)
248 b(6,21)= up*(c1b3+c2b3)
249 b(6,22)= um*(c1a1+c2a1)
250 b(6,23)= um*(c1a2+c2a2)
251 b(6,24)= um*(c1a3+c2a3)
252
253
254
255 ca11 = invj(3,i)
256 ca21 = invj(1,i)
257 ca31 = invj(2,i)
258 ca12 = invj(6,i)
259 ca22 = invj(4,i)
260 ca32 = invj(5,i)
261 ca13 = invj(9,i)
262 ca23 = invj(7,i)
263 ca33 = invj(8,i)
264
265 c(1,1) = ca11*ca11
266 c(2,1) = ca12*ca12
267 c(3,1) = ca13*ca13
268 c(4,1) = two*ca11*ca12
269 c(5,1) = two*ca11*ca13
270 c(6,1) = two*ca12*ca13
271 c(1,2) = ca21*ca21
272 c(2,2) = ca22*ca22
273 c(3,2) = ca23*ca23
274 c(4,2) = two*ca21*ca22
275 c(5,2) = two*ca21*ca23
276 c(6,2) = two*ca22*ca23
277 c(1,3) = ca31*ca31
278 c(2,3) = ca32*ca32
279 c(3,3) = ca33*ca33
280 c(4,3) = two*ca31*ca32
281 c(5,3) = two*ca31*ca33
282 c(6,3) = two*ca32*ca33
283 c(1,4) = ca11*ca21
284 c(2,4) = ca12*ca22
285 c(3,4) = ca13*ca23
286 c(4,4) = ca12*ca21+ca11*ca22
287 c(5,4) = ca13*ca21+ca11*ca23
288 c(6,4) = ca12*ca23+ca13*ca22
289 c(1,5) = ca11*ca31
290 c(2,5) = ca12*ca32
291 c(3,5) = ca13*ca33
292 c(4,5) = ca12*ca31+ca11*ca32
293 c(5,5) = ca13*ca31+ca11*ca33
294 c(6,5) = ca12*ca33+ca13*ca32
295 c(1,6) = ca21*ca31
296 c(2,6) = ca22*ca32
297 c(3,6) = ca23*ca33
298 c(4,6) = ca22*ca31+ca21*ca32
299 c(5,6) = ca23*ca31+ca21*ca33
300 c(6,6) = ca23*ca32+ca22*ca33
301
302 d11 =b(1,1)*ulx1(i)+b(1,2)*uly1(i)+b(1,3)*ulz1(i)
303 . +b(1,4)*ulx2(i)+b(1,5)*uly2(i)+b(1,6)*ulz2(i)
304 . +b(1,7)*ulx3(i)+b(1,8)*uly3(i)+b(1,9)*ulz3(i)
305 . +b(1,10)*ulx4(i)+b(1,11)*uly4(i)+b(1,12)*ulz4(i)
306 . +b(1,13)*ulx5(i)+b(1,14)*uly5(i)+b(1,15)*ulz5(i)
307 . +b(1,16)*ulx6(i)+b(1,17)*uly6(i)+b(1,18)*ulz6(i)
308 . +b(1,19)*ulx7(i)+b(1,20)*uly7(i)+b(1,21)*ulz7(i)
309 . +b(1,22)*ulx8(i)+b(1,23)*uly8(i)+b(1,24)*ulz8(i)
310 d22 =b(2,1)*ulx1(i)+b(2,2)*uly1(i)+b(2,3)*ulz1(i)
311 . +b(2,4)*ulx2(i)+b(2,5)*uly2(i)+b(2,6)*ulz2(i)
312 . +b(2,7)*ulx3(i)+b(2,8)*uly3(i)+b(2,9)*ulz3(i)
313 . +b(2,10)*ulx4(i)+b(2,11)*uly4(i)+b(2,12)*ulz4(i)
314 . +b(2,13)*ulx5(i)+b(2,14)*uly5(i)+b(2,15)*ulz5(i)
315 . +b(2,16)*ulx6(i)+b(2,17)*uly6(i)+b(2,18)*ulz6(i)
316 . +b(2,19)*ulx7(i)+b(2,20)*uly7(i)+b(2,21)*ulz7(i)
317 . +b(2,22)*ulx8(i)+b(2,23)*uly8(i)+b(2,24)*ulz8(i)
318 d33 =b(3,1)*ulx1(i)+b(3,2)*uly1(i)+b(3,3)*ulz1(i)
319 . +b(3,4)*ulx2(i)+b(3,5)*uly2(i)+b(3,6)*ulz2(i)
320 . +b(3,7)*ulx3(i)+b(3,8)*uly3(i)+b(3,9)*ulz3(i)
321 . +b(3,10)*ulx4(i)+b(3,11)*uly4(i)+b(3,12)*ulz4(i)
322 . +b(3,13)*ulx5(i)+b(3,14)*uly5(i)+b(3,15)*ulz5(i)
323 . +b(3,16)*ulx6(i)+b(3,17)*uly6(i)+b(3,18)*ulz6(i)
324 . +b(3,19)*ulx7(i)+b(3,20)*uly7(i)+b(3,21)*ulz7(i)
325 . +b(3,22)*ulx8(i)+b(3,23)*uly8(i)+b(3,24)*ulz8(i)
326 d12 =b(4,1)*ulx1(i)+b(4,2)*uly1(i)+b(4,3)*ulz1(i)
327 . +b(4,4)*ulx2(i)+b(4,5)*uly2(i)+b(4,6)*ulz2(i)
328 . +b(4,7)*ulx3(i)+b(4,8)*uly3(i)+b(4,9)*ulz3(i)
329 . +b(4,10)*ulx4(i)+b(4,11)*uly4(i)+b(4,12)*ulz4(i)
330 . +b(4,13)*ulx5(i)+b(4,14)*uly5(i)+b(4,15)*ulz5(i)
331 . +b(4,16)*ulx6(i)+b(
332 . +b(4,19)*ulx7(i)+b(4,20)*uly7(i)+b(4,21)*ulz7(i)
333 . +b(4,22)*ulx8(i)+b(4,23)*uly8(i)+b(4,24)*ulz8(i)
334 d13 =b(5,1)*ulx1(i)+b(5,2)*uly1(i)+b(5,3)*ulz1(i)
335 . +b(5,4)*ulx2(i)+b(5,5)*uly2(i)+b(5,6)*ulz2(i)
336 . +b(5,7)*ulx3(i)+b(5,8)*uly3(i)+b(5,9)*ulz3(i)
337 . +b(5,10)*ulx4(i)+b(5,11)*uly4(i)+b(5,
338 . +b(5,13)*ulx5(i)+b(5,14)*uly5(i)+b(5,15)*ulz5(i)
339 . +b(5,16)*ulx6(i)+b(5,17)*uly6(i)+b(5,18)*ulz6(i)
340 . +b(5,19)*ulx7(i)+b(5,20)*uly7(i)+b(5,21)*ulz7(i)
341 . +b(5,22)*ulx8(i)+b(5,23)*uly8(i)+b(5,24)*ulz8(i)
342 d23 =b(6,1)*ulx1(i)+b(6,2)*uly1(i)+b(6,3)*ulz1(i)
343 . +b(6,4)*ulx2(i)+b(6,5)*uly2(i)+b(6,6)*ulz2(i)
344 . +b(6,7)*ulx3(i)+b(6,8)*uly3(i)+b(6,9)*ulz3(i)
345 . +b(6,10)*ulx4(i)+b(6,11)*uly4(i)+b(6,12)*ulz4(i)
346 . +b(6,13)*ulx5(i)+b(6,14)*uly5(i)+b(6,15)*ulz5(i)
347 . +b(6,16)*ulx6(i)+b(6,17)*uly6(i)+b(6,18)*ulz6(i)
348 . +b(6,19)*ulx7(i)+b(6,20)*uly7(i)+b(6,21)*ulz7(i)
349 . +b(6,22)*ulx8(i)+b(6,23)*uly8(i)+b(6,24)*ulz8(i)
350
351 dxx(i) =c(1,1)*d11+c(1,2)*d22+c(1,3)*d33
352 . +c(1,4)*d12+c(1,5)*d13+c(1,6)*d23
353 dyy(i) =c(2,1)*d11+c(2,2)*d22+c(2,3)*d33
354 . +c(2,4)*d12+c(2,5)*d13+c(2,6)*d23
355 dzz(i) =c(3,1)*d11+c(3,2)*d22+c(3,3)*d33
356 . +c(3,4)*d12+c(3,5)*d13+c(3,6)*d23
357 d4(i) =c(4,1)*d11+c(4,2)*d22+c(4,3)*d33
358 . +c(4,4)*d12+c(4,5)*d13+c(4,6)*d23
359 d6(i) =c(5,1)*d11+c(5,2)*d22+c(5,3)*d33
360 . +c(5,4)*d12+c(5,5)*d13+c(5,6)*d23
361 d5(i) =c(6,1)*d11+c(6,2)*d22+c(6,3)*d33
362 . +c(6,4)*d12+c(6,5)*d13+c(6,6)*d23
363 DO l=1,6
364 DO m=1,24
365 bb(l,m,i) = zero
366 DO lm=1,6
367 bb(l,m,i) = bb(l,m,i)+c(l,lm)*b(lm,m)
368 ENDDO
369 ENDDO
370 ENDDO
371 dvc(i) =zero
372 ENDDO
373
374 IF (icp==2) THEN
375 DO i=1,nel
376 dvc(i)=(dsv(i)-dvc(i))*fac(i)*dt1
377 ENDDO
378 ELSEIF (icp==1) THEN
379 DO i=1,nel
380 dvc(i)=(dsv(i)-dvc(i))*dt1
381 ENDDO
382 ENDIF
383 IF ((icp>0.AND.ismstr/=10.AND.ismstr/=12).OR.i_sh>1) THEN
384 DO i=1,nel
385 dv =dvc(i)*off(i)
386 IF(
idege(i)>10) dv = zero
387 sdv(i) =dv
388 IF (dv>tol) THEN
389 dv =zero
390 ENDIF
391 IF(offs(i)==two.OR.ismstr==11) cycle
392 dv1 = one- dv
393 volo(i) = volo(i)*dv1
394 eint(i) = eint(i)/dv1
395
396 ENDDO
397 ELSE
398 DO i=1,nel
399 sdv(i) =zero
400 ENDDO
401 ENDIF
402
403 dt1d2=half*dt1
404 IF (iscau>0)dt1d2=dt1
405 dt1d=two*dt1d2
406 udt1 = one
407 IF (dt1 > zero) udt1 = one/dt1
408 DO i=1,nel
409 dxx(i)= dxx(i)*udt1
410 dyy(i)= dyy(i)*udt1
411 dzz(i)= dzz(i)*udt1
412 d4(i) = d4(i)*udt1
413 d5(i) = d5(i)*udt1
414 d6(i) = d6(i)*udt1
415 ENDDO
416
417 RETURN
subroutine idege(x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, a, amax, fac, it4, it, indx, n_indx)