45
46
47
49
50
51
52#include "implicit_f.inc"
53#include "comlock.inc"
54
55
56
57#include "mvsiz_p.inc"
58
59
60
61#include "impl1_c.inc"
62#include "scr17_c.inc"
63#include "scr07_c.inc"
64#include "scr18_c.inc"
65
66
67
68 INTEGER, INTENT(IN) :: NPT
69 INTEGER, INTENT(IN) :: ISMSTR
70 INTEGER, INTENT(IN) :: JLAG
71 INTEGER NGL(*), NC(MVSIZ,10), NEL
72
73 double precision
74 . xx(mvsiz,10), yy(mvsiz,10), zz(mvsiz,10),sav(nel,30),voldp(mvsiz,5)
76 . off(nel),vol(mvsiz,5),deltax(*),deltax2(*),
77 . rx(*),ry(*),rz(*), sx(*),sy(*),sz(*), tx(*),ty(*),tz(*),
78 . nx(mvsiz,10,5),voln(*),volg(mvsiz),
79 . px(mvsiz,10,5),py(mvsiz,10,5),pz(mvsiz,10,5),
80 . wip(5),alph(5),beta(5),offg(nel)
81
82
83
84 INTEGER I,IP,N,K1,K2,K3,K4,,K6,K7,K8,K9,K10,
85 . M,IPERM(10,4),ICOR,NNEGA,INDEX(MVSIZ),J,NN
86 INTEGER ITER,NITER
87 double precision
88 . xa(mvsiz,10),ya(mvsiz,10),za(mvsiz,10),
89 . xb(mvsiz,10),yb(mvsiz,10),zb(mvsiz,10),
90 . a4, b4, a4m1, b4m1
91 DATA iperm/
92 . 2, 4, 3, 1, 9,10, 6, 5, 8, 7,
93 . 4, 1, 3, 2, 8, 7,10, 9, 5, 6,
94 . 1, 4, 2, 3, 8, 9, 5, 7,10, 6,
95 . 1, 2, 3, 4, 5, 6, 7, 8, 9,10/
96
97 a4 = four * alph(1)
98 b4 = four * beta(1)
99 a4m1 = a4 - one
100 b4m1 = b4 - one
101
102 DO i=1,nel
103 rx(i) = xx(i,1) - xx(i,4)
104 ry(i) = yy(i,1) - yy(i,4)
105 rz(i) = zz(i,1) - zz(i,4)
106 sx(i) = xx(i,2) - xx(i,4)
107 sy(i) = yy(i,2) - yy(i,4)
108 sz(i) = zz(i,2) - zz(i,4)
109 tx(i) = xx(i,3) - xx(i,4)
110 ty(i) = yy(i,3) - yy(i,4)
111 tz(i) = zz(i,3) - zz(i,4)
112 volg(i) =zero
113 ENDDO
114
115 DO n=1,4
116 DO i=1,nel
117 xa(i,n) = a4m1*xx(i,n)
118 ya(i,n) = a4m1*yy(i,n)
119 za(i,n) = a4m1*zz(i,n)
120
121 xb(i,n) = b4m1*xx(i,n)
122 yb(i,n) = b4m1*yy(i,n)
123 zb(i,n) = b4m1*zz(i,n)
124 ENDDO
125 ENDDO
126
127 DO n=5,10
128 DO i=1,nel
129 xa(i,n) = a4*xx(i,n)
130 ya(i,n) = a4*yy(i,n)
131 za(i,n) = a4*zz(i,n)
132
133 xb(i,n) = b4*xx(i,n)
134 yb(i,n) = b4*yy(i,n)
135 zb(i,n) = b4*zz(i,n)
136 ENDDO
137 ENDDO
138
139 DO ip=1,4
140 k1 = iperm(1,ip)
141 k2 = iperm(2,ip)
142 k3 = iperm(3,ip)
143 k4 = iperm(4,ip)
144 k5 = iperm(5,ip)
145 k6 = iperm(6,ip)
146 k7 = iperm(7,ip)
147 k8 = iperm(8,ip)
148 k9 = iperm(9,ip)
149 k10= iperm(10,ip)
151 1 alph(ip), beta(ip), wip(ip), xb(1,k1),
152 2 xb(1,k2), xb(1,k3), xa(1,k4), xb(1,k5),
153 3 xb(1,k6), xb(1,k7), xb(1,k8), xb(1,k9),
154 4 xb(1,k10), xa(1,k8), xa(1,k9), xa(1,k10),
155 5 yb(1,k1), yb(1,k2), yb(1,k3), ya(1,k4),
156 6 yb(1,k5), yb(1,k6), yb(1,k7), yb(1,k8),
157 7 yb(1,k9), yb(1,k10), ya(1,k8), ya(1,k9),
158 8 ya(1,k10), zb(1,k1), zb(1,k2), zb(1,k3),
159 9 za(1,k4), zb(1,k5), zb(1,k6), zb(1,k7),
160 a zb(1,k8), zb(1,k9), zb(1,k10), za(1,k8),
161 b za(1,k9), za(1,k10), px(1,k1,ip), px(1,k2,ip),
162 c px(1,k3,ip), px(1,k4,ip), px(1,k5,ip), px(1,k6,ip),
163 d px(1,k7,ip), px(1,k8,ip), px(1,k9,ip), px(1,k10,ip),
164 e py(1,k1,ip), py(1,k2,ip), py(1,k3,ip), py(1,k4,ip),
165 f py(1,k5,ip), py(1,k6,ip), py(1,k7,ip), py(1,k8,ip),
166 g py(1,k9,ip), py(1,k10,ip),pz(1,k1,ip), pz(1,k2,ip),
167 h pz(1,k3,ip), pz(1,k4,ip), pz(1,k5,ip), pz(1,k6,ip),
168 i pz(1,k7,ip), pz(1,k8,ip), pz(1,k9,ip), pz(1,k10,ip),
169 j nx(1,k1,ip), nx(1,k2,ip), nx(1,k3,ip), nx(1,k4,ip),
170 k nx(1,k5,ip), nx(1,k6,ip), nx(1,k7,ip), nx(1,k8,ip),
171 l nx(1,k9,ip), nx(1,k10,ip),vol(1,ip), voldp(1,ip),
172 m nel, offg)
173
174 ENDDO
175
176
177 IF(npt==5)THEN
178 ip = 5
179
180 DO i=1,nel
181 xa(i,1) = zero
182 ENDDO
183
185 1 alph(ip), beta(ip), wip(ip), xa(1,1),
186 2 xa(1,1), xa(1,1), xa(1,1), xx(1,k5),
187 3 xx(1,k6), xx(1,k7), xx(1,k8), xx(1,k9),
188 4 xx(1,k10), xx(1,k8), xx(1,k9), xx(1,k10),
189 5 xa(1,1), xa(1,1), xa(1,1), xa(1,1),
190 6 yy(1,k5), yy(1,k6), yy(1,k7), yy(1,k8),
191 7 yy(1,k9), yy(1,k10), yy(1,k8), yy(1,k9),
192 8 yy(1,k10), xa(1,1), xa(1,1), xa(1,1),
193 9 xa(1,1), zz(1,k5), zz(1,k6), zz(1,k7),
194 a zz(1,k8), zz(1,k9), zz(1,k10), zz(1,k8),
195 b zz(1,k9), zz(1,k10), px(1,k1,ip), px(1,k2,ip),
196 c px(1,k3,ip), px(1,k4,ip), px(1,k5,ip), px(1,k6,ip),
197 d px(1,k7,ip), px(1,k8,ip), px(1,k9,ip), px(1,k10,ip),
198 e py(1,k1,ip), py(1,k2,ip), py(1,k3,ip), py(1,k4,ip),
199 f py(1,k5,ip), py(1,k6,ip), py(1,k7,ip), py(1,k8,ip),
200 g py(1,k9,ip), py(1,k10,ip),pz(1,k1,ip), pz(1,k2,ip),
201 h pz(1,k3,ip), pz(1,k4,ip), pz(1,k5,ip), pz(1,k6,ip),
202 i pz(1,k7,ip), pz(1,k8,ip), pz(1,k9,ip), pz(1,k10,ip),
203 j nx(1,k1,ip), nx(1,k2,ip), nx(1,k3,ip), nx(1,k4,ip),
204 k nx(1,k5,ip), nx(1,k6,ip), nx(1,k7,ip), nx(1,k8,ip),
205 l nx(1,k9,ip), nx(1,k10,ip),vol(1,ip), voldp(1,ip),
206 m nel, offg)
207 ENDIF
208
209
210 nnega = 0
211 IF(jlag/=0.AND.(ismstr==10.OR.(ismstr==12.AND.idtmin(1)/=3))) THEN
212 DO i=1,nel
213 IF(offg(i) > one)THEN
214 nnega=nnega+1
215 index(nnega)=i
216 END IF
217 ENDDO
218 END IF
219 icor=0
220 DO ip=1,npt
221 DO i=1,nel
222 IF(off(i)==0.)THEN
223 vol(i,ip)=one
224 ELSEIF(off(i)> one)THEN
225 ELSEIF(vol(i,ip)<=zero)THEN
226 icor=1
227 ENDIF
228 ENDDO
229 ENDDO
230
231 IF(jlag/=0)THEN
232 IF(icor/=0.AND.inconv==1)THEN
233 DO ip=1,npt
234 DO i=1,nel
235 IF(off(i) == zero.OR.offg(i) > one)THEN
236 ELSEIF(vol(i,ip)<=zero)THEN
237 nnega=nnega+1
238 index(nnega)=i
239#include "lockon.inc"
240 IF(ismstr<10) THEN
241 CALL ancmsg(msgid=260,anmode=aninfo,
242 . i1=ngl(i))
243 ELSE
244 CALL ancmsg(msgid=262,anmode=aninfo,
245 . i1=ngl(i))
246 END IF
247#include "lockoff.inc"
248 offg(i) = two
249 ENDIF
250 ENDDO
251 ENDDO
253 ENDIF
254
255 IF(nnega >0 )THEN
256 DO n=1,10
257#include "vectorize.inc"
258 DO j=1,nnega
259 i = index(j)
260 nn = nc(i,n)
261 xx(i,n)=sav(i,n)
262 yy(i,n)=sav(i,n+10)
263 zz(i,n)=sav(i,n+20)
264 ENDDO
265 END DO
266#include "vectorize.inc"
267 DO j=1,nnega
268 i = index(j)
269 rx(i) = xx(i,1) - xx(i,4)
270 ry(i) = yy(i,1) - yy(i,4)
271 rz(i) = zz(i,1) - zz(i,4)
272 sx(i) = xx(i,2) - xx(i,4)
273 sy(i) = yy(i,2) - yy(i,4)
274 sz(i) = zz(i,2) - zz(i,4)
275 tx(i) = xx(i,3) - xx(i,4)
276 ty(i) = yy(i,3) - yy(i,4)
277 tz(i) = zz(i,3) - zz(i,4)
278 ENDDO
279
280 DO n=1,4
281#include "vectorize.inc"
282 DO j=1,nnega
283 i = index(j)
284 xa(i,n) = a4m1*xx(i,n)
285 ya(i,n) = a4m1*yy(i,n)
286 za(i,n) = a4m1*zz(i,n)
287
288 xb(i,n) = b4m1*xx(i,n)
289 yb(i,n) = b4m1*yy(i,n)
290 zb(i,n) = b4m1*zz(i,n)
291 ENDDO
292 ENDDO
293
294 DO n=5,10
295#include "vectorize.inc"
296 DO j=1,nnega
297 i = index(j)
298 xa(i,n) = a4*xx(i,n)
299 ya(i,n) = a4*yy(i,n)
300 za(i,n) = a4*zz(i,n)
301
302 xb(i,n) = b4*xx(i,n)
303 yb(i,n) = b4*yy(i,n)
304 zb(i,n) = b4*zz(i,n)
305 ENDDO
306 ENDDO
307 DO ip=1,4
308 k1 = iperm(1,ip)
309 k2 = iperm(2,ip)
310 k3 = iperm(3,ip)
311 k4 = iperm(4,ip)
312 k5 = iperm(5,ip)
313 k6 = iperm(6,ip)
314 k7 = iperm(7,ip)
315 k8 = iperm(8,ip)
316 k9 = iperm(9,ip)
317 k10= iperm(10,ip)
318 CALL s10jacob1(alph(ip),beta(ip),wip(ip),
319 . xb(1,k1),xb(1,k2),xb(1,k3),xa(1,k4),xb(1,k5),xb(1,k6),xb(1,k7),
320 . xb(1,k8),xb(1,k9),xb(1,k10),xa(1,k8),xa(1,k9),xa(1,k10),
321 . yb(1,k1),yb(1,k2),yb(1,k3),ya(1,k4),yb(1,k5),yb(1,k6),yb(1,k7),
322 . yb(1,k8),yb(1,k9),yb(1,k10),ya(1,k8),ya(1,k9),ya(1,k10),
323 . zb(1,k1),zb(1,k2),zb(1,k3),za(1,k4),zb(1,k5),zb(1,k6),zb(1,k7),
324 . zb(1,k8),zb(1,k9),zb(1,k10),za(1,k8),za(1,k9),za(1,k10),
325 . px(1,k1,ip) ,px(1,k2,ip),px(1,k3,ip),px(1,k4,ip),px(1,k5,ip),
326 . px(1,k6,ip) ,px(1,k7,ip),px(1,k8,ip),px(1,k9,ip),px(1,k10,ip),
327 . py(1,k1,ip) ,py(1,k2,ip),py(1,k3,ip),py(1,k4,ip),py(1,k5,ip),
328 . py(1,k6,ip) ,py(1,k7,ip),py(1,k8,ip),py(1,k9,ip),py(1,k10,ip),
329 . pz(1,k1,ip) ,pz(1,k2,ip),pz(1,k3,ip),pz(1,k4,ip),pz(1,k5,ip),
330 . pz(1,k6,ip) ,pz(1,k7,ip),pz(1,k8,ip),pz(1,k9,ip),pz(1,k10,ip),
331 . nx(1,k1,ip) ,nx(1,k2,ip),nx(1,k3,ip),nx(1,k4,ip),nx(1,k5,ip),
332 . nx(1,k6,ip) ,nx(1,k7,ip),nx(1,k8,ip),nx(1,k9,ip),nx(1,k10,ip),
333 . vol(1,ip) ,nnega, index ,voldp(1,ip))
334 ENDDO
335
336
337 IF(npt==5)THEN
338 ip = 5
339
340#include "vectorize.inc"
341 DO j=1,nnega
342 i = index(j)
343 xa(i,1) = zero
344 ENDDO
345
346 CALL s10jacob1(alph(ip),beta(ip),wip(ip),
347 . xa(1,1) ,xa(1,1) ,xa(1,1) ,xa(1,1) ,xx(1,k5),
348 . xx(1,k6),xx(1,k7),xx(1,k8),xx(1,k9),xx(1,k10),
349 . xx(1,k8),xx(1,k9),xx(1,k10),
350 . xa(1,1) ,xa(1,1) ,xa(1,1) ,xa(1,1) ,yy(1,k5),
351 . yy(1,k6),yy(1,k7),yy(1,k8),yy(1,k9),yy(1,k10),
352 . yy(1,k8),yy(1,k9),yy(1,k10),
353 . xa(1,1) ,xa(1,1) ,xa(1,1) ,xa(1,1) ,zz(1,k5),
354 . zz(1,k6),zz(1,k7),zz(1,k8),zz(1,k9),zz(1,k10),
355 . zz(1,k8),zz(1,k9),zz(1,k10),
356 . px(1,k1,ip) ,px(1,k2,ip),px(1,k3,ip),px(1,k4,ip),px(1,k5,ip
357 . px(1,k6,ip) ,px(1,k7,ip),px(1,k8,ip),px(1,k9,ip),px(1,k10,ip),
358 . py(1,k1,ip) ,py(1,k2,ip),py(1,k3,ip),py(1,k4,ip),py(1,k5,ip),
359 . py(1,k6,ip) ,py(1,k7,ip),py(1,k8,ip),py(1,k9,ip),py(1,k10,ip),
360 . pz(1,k1,ip) ,pz(1,k2,ip),pz(1,k3,ip),pz(1,k4,ip),pz(1,k5,ip),
361 . pz(1,k6,ip) ,pz(1,k7,ip),pz(1,k8,ip),pz(1,k9,ip),pz(1,k10,ip),
362 . nx(1,k1,ip) ,nx(1,k2,ip),nx(1,k3,ip),nx(1,k4,ip),nx(1,k5,ip),
363 . nx(1,k6,ip) ,nx(1,k7,ip),nx(1,k8,ip),nx(1,k9,ip),nx(1,k10,ip),
364 . vol(1,ip) ,nnega, index ,voldp(1,ip))
365 ENDIF
366 IF(ineg_v==0)THEN
367 CALL ancmsg(msgid=280,anmode=aninfo)
368 mstop = 1
369 END IF
370 ENDIF
371
372 DO ip=1,npt
373 DO i=1,nel
374 volg(i) =volg(i) + vol(i,ip)
375 ENDDO
376 ENDDO
377
378 DO i=1,nel
379 voln(i) =volg(i)/npt
380 ENDDO
381
382 1000 FORMAT(/' ZERO OR NEGATIVE VOLUME : 10 NODES TETRAHEDRON NB ',i10,
383 . ' INTEGRATION POINT NB ',i1/)
384 1100 FORMAT(/' ZERO OR NEGATIVE VOLUME : 4 NODES TETRAHEDRON NB ',i10,
385 . ' INTEGRATION POINT NB ',i1/)
386 2000 FORMAT(/' ZERO OR NEGATIVE VOLUME : DELETE 3D-ELEMENT NB',i10/)
387 3000 FORMAT(/' ZERO OR NEGATIVE VOLUME : 3D-ELEMENT NB:',i10/,
388 + ' ELEMENT IS SWITCHED TO SMALL STRAIN OPTION'/)
389 4000 FORMAT(/' ZERO OR NEGATIVE VOLUME : 3D-ELEMENT NB:',i10/)
390
391 RETURN
if(complex_arithmetic) id
subroutine s10jacob1(alph, beta, w, x1b, x2b, x3b, x4a, x5b, x6b, x7b, x8b, x9b, x10b, x8a, x9a, x10a, y1b, y2b, y3b, y4a, y5b, y6b, y7b, y8b, y9b, y10b, y8a, y9a, y10a, z1b, z2b, z3b, z4a, z5b, z6b, z7b, z8b, z9b, z10b, z8a, z9a, z10a, px1, px2, px3, px4, px5, px6, px7, px8, px9, px10, py1, py2, py3, py4, py5, py6, py7, py8, py9, py10, pz1, pz2, pz3, pz4, pz5, pz6, pz7, pz8, pz9, pz10, nx1, nx2, nx3, nx4, nx5, nx6, nx7, nx8, nx9, nx10, vol, nnega, index, voldp)
subroutine s10jacob(alph, beta, w, x1b, x2b, x3b, x4a, x5b, x6b, x7b, x8b, x9b, x10b, x8a, x9a, x10a, y1b, y2b, y3b, y4a, y5b, y6b, y7b, y8b, y9b, y10b, y8a, y9a, y10a, z1b, z2b, z3b, z4a, z5b, z6b, z7b, z8b, z9b, z10b, z8a, z9a, z10a, px1, px2, px3, px4, px5, px6, px7, px8, px9, px10, py1, py2, py3, py4, py5, py6, py7, py8, py9, py10, pz1, pz2, pz3, pz4, pz5, pz6, pz7, pz8, pz9, pz10, nx1, nx2, nx3, nx4, nx5, nx6, nx7, nx8, nx9, nx10, vol, voldp)
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)