44
45
46
47 USE elbufdef_mod
48
49
50
51#include "implicit_f.inc"
52
53
54
55#include "scr17_c.inc"
56#include "com08_c.inc"
57#include "units_c.inc"
58#include "comlock.inc"
59
60
61
62 INTEGER NEL,NUPARAM,NUVAR,NFUNC,INLOC,
63 . IPG,NPG,NPF(*),NGL(NEL),
64 . IFUNC(NFUNC)
66 . time,timestep,tf(*),uparam(nuparam)
67 my_real,
DIMENSION(NEL),
INTENT(IN) ::
68 . rho,dt,
69 . depsxx,depsyy,depszz,depsxy,depsyz,depszx,
70 . epspxx,epspyy,epspzz,epspxy,epspyz,epspzx,
71 . sigoxx,sigoyy,sigozz,sigoxy,sigoyz,sigozx
72 my_real ,
DIMENSION(NEL),
INTENT(OUT) ::
73 . soundsp,signxx,signyy,signzz,signxy,signyz,signzx
75 . sigy,et
76 my_real ,
DIMENSION(NEL),
INTENT(INOUT) ::
77 . pla,dpla,epsd,varnl,loff,off
78 my_real ,
DIMENSION(NEL,NUVAR),
INTENT(INOUT) ::
79 . uvar
80 TYPE (ELBUF_STRUCT_), TARGET :: ELBUF_TAB
81
82
83
84 INTEGER I,J,IRES,Ifail,NINDX,NINDX2,INDX(NEL),INDX2(NEL),
85 . IPOS(NEL),IAD(NEL),ILEN(NEL),IR,IS,IT
86 my_real dtmin,xscale_fail,yscale_fail,s1,s2,q,s11,
87 . s22,r_inter,dfdepsd(nel),fail(nel),seq(nel),
88 . i1,i2,i3,r,psi,s33
89
90
91 ires = nint(uparam(11))
92
93
94 ifail = nint(uparam(13))
95
96
97
98
99
100 dtmin = uparam(15)
101 xscale_fail = uparam(22)
102 yscale_fail = uparam(23)
103
104 SELECT CASE (ires)
105
106 CASE(1)
107
109 1 nel ,ngl ,nuparam ,nuvar ,nfunc ,ifunc ,npf ,
110 2 tf ,timestep,time ,uparam ,uvar ,rho ,pla ,
111 3 dpla ,soundsp ,epsd ,off ,
112 4 depsxx ,depsyy ,depszz ,depsxy ,depsyz ,depszx ,
113 5 epspxx ,epspyy ,epspzz ,epspxy ,epspyz ,epspzx ,
114 6 sigoxx ,sigoyy ,sigozz ,sigoxy ,sigoyz ,sigozx ,
115 7 signxx ,signyy ,signzz ,signxy ,signyz ,signzx ,
116 8 sigy ,et ,seq )
117
118 CASE(2)
119
121 1 nel ,ngl ,nuparam ,nuvar ,nfunc ,ifunc ,npf ,
122 2 tf ,timestep,time ,uparam ,uvar ,rho ,pla ,
123 3 dpla ,soundsp ,epsd ,off ,
124 4 depsxx ,depsyy ,depszz ,depsxy ,depsyz ,depszx ,
125 5 epspxx ,epspyy ,epspzz ,epspxy ,epspyz ,epspzx ,
126 6 sigoxx ,sigoyy ,sigozz ,sigoxy ,sigoyz ,sigozx ,
127 7 signxx ,signyy ,signzz ,signxy ,signyz ,signzx ,
128 8 sigy ,et ,seq )
129
130 END SELECT
131
132
133
134
135
136
137 IF (ifunc(4) > 0) THEN
138 ipos(1:nel) = 1
139 iad(1:nel) = npf(ifunc(4)) / 2 + 1
140 ilen(1:nel) = npf(ifunc(4)+1) / 2 - iad(1:nel) - ipos(1:nel)
141 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_fail,dfdepsd,fail)
142 fail(1:nel) = yscale_fail*fail(1:nel)
143 ENDIF
144
145
146 nindx = 0
147 nindx2 = 0
148 IF (dtmin > zero .OR. ifunc(4) > 0) THEN
149 DO i = 1,nel
150
151
152 IF ((dt(i) > em20).AND.(dt(i) < dtmin).AND.(off(i) == one)) THEN
153 off(i) = zero
154 nindx = nindx+1
155 indx(nindx) = i
156 ENDIF
157
158
159 IF ((ifunc(4) > 0).AND.(off(i) == one)) THEN
160
161 IF (loff(i) < em01) loff(i) = zero
162 IF (loff(i) < one) loff(i) = loff(i)*four_over_5
163
164
165 IF (npg > 1) THEN
166
167 IF (ipg == npg) THEN
168
169 off(i) = zero
170
171 DO ir = 1, elbuf_tab%NPTR
172 DO is = 1, elbuf_tab%NPTS
173 DO it = 1, elbuf_tab%NPTT
174
175 IF (elbuf_tab%BUFLY(1)%LBUF(ir,is,it)%OFF(i)>zero) off(i) = one
176 ENDDO
177 ENDDO
178 ENDDO
179 ENDIF
180
181 ELSE
182
183 IF (ipg == 1) THEN
184 off(i) = zero
185 ENDIF
186
187 IF (loff(i)>zero) off(i) = one
188 ENDIF
189
190 IF (loff(i) == one) THEN
191
192 IF (ifail == 0) THEN
193 IF (seq(i) >= fail(i)) loff(i) = four_over_5
194
195 ELSEIF (ifail == 1) THEN
196 IF (pla(i) >= fail(i)) loff(i) = four_over_5
197
198 ELSEIF (ifail == 2) THEN
199
200 i1 = signxx(i)+signyy(i)+signzz(i)
201 i2 = signxx(i)*signyy(i)+signyy(i)*signzz(i)+signzz(i)*signxx(i)-
202 . signxy(i)*signxy(i)-signzx(i)*signzx(i)-signyz(i)*signyz(i)
203 i3 = signxx(i)*signyy(i)*signzz(i)-signxx(i)*signyz(i)*signyz(i)-
204 . signyy(i)*signzx(i)*signzx(i)-signzz(i)*signxy
205 . two*signxy(i)*signzx(i)*signyz(i)
206 q = (three*i2 - i1*i1)/nine
207 r = (two*i1*i1*i1-nine*i1*i2+twenty7*i3)/cinquante4
208 r_inter =
min(r/sqrt(
max(em20,(-q**3))),one)
209 psi = acos(
max(r_inter,-one))
210 s11 = two*sqrt(-q)*cos(psi/three)+third*i1
211 s22 = two*sqrt(-q)*cos((psi+two*pi)/three)+third*i1
212 s33 = two*sqrt(-q)*cos((psi+four*pi)/three)+third*i1
213
214 IF (s11 < s22) THEN
215 r_inter = s11
216 s11 = s22
217 s22 = r_inter
218 ENDIF
219 IF (s22 < s33)THEN
220 r_inter = s22
221 s22 = s33
222 s33 = r_inter
223 ENDIF
224 IF (s11 < s22)THEN
225 r_inter = s11
226 s11 = s22
227 s22 = r_inter
228 ENDIF
229 IF ((s11>=fail(i)).OR.(abs(s33)>=fail(i))) loff(i) = four_over_5
230
231 ELSEIF (ifail == 3) THEN
232
233 i1 = signxx(i)+signyy(i)+signzz(i)
234 i2 = signxx(i)*signyy(i)+signyy(i)*signzz(i)+signzz(i)*signxx(i)-
235 . signxy(i)*signxy(i)-signzx(i)*signzx(i)-signyz(i)*signyz(i)
236 i3 = signxx(i)*signyy(i)*signzz(i)-signxx(i)*signyz(i)*signyz(i)-
237 . signyy(i)*signzx(i)*signzx(i)-signzz(i)*signxy(i)*signxy(i)+
238 . two*signxy(i)*signzx(i)*signyz(i)
239 q = (three*i2 - i1*i1)/nine
240 r = (two*i1*i1*i1-nine*i1*i2+twenty7*i3)/cinquante4
241 r_inter =
min(r/sqrt(
max(em20,(-q**3))),one)
242 psi = acos(
max(r_inter,-one))
243 s11 = two*sqrt(-q)*cos(psi/three)+third*i1
244 s22 = two*sqrt(-q)*cos((psi+two*pi)/three)+third*i1
245 s33 = two*sqrt(-q)*cos((psi+four*pi)/three)+third*i1
246
247 IF (s11 < s22) THEN
248 r_inter = s11
249 s11 = s22
250 s22 = r_inter
251 ENDIF
252 IF (s22 < s33)THEN
253 r_inter = s22
254 s22 = s33
255 s33 = r_inter
256 ENDIF
257 IF (s11 < s22)THEN
258 r_inter = s11
259 s11 = s22
260 s22 = r_inter
261 ENDIF
262 IF (s11>=fail(i)) loff(i) = four_over_5
263 ENDIF
264
265 IF (loff(i) == four_over_5) THEN
266 nindx2 = nindx2+1
267 indx2(nindx2) = i
268 ENDIF
269 ENDIF
270 ENDIF
271 ENDDO
272 ENDIF
273
274
275 IF (nindx>0) THEN
276 DO j=1,nindx
277#include "lockon.inc"
278 WRITE(iout, 1000) ngl(indx(j))
279 WRITE(istdo,1100) ngl(indx(j)),tt
280#include "lockoff.inc"
281 ENDDO
282 ENDIF
283
284
285 IF (nindx2>0) THEN
286 DO j=1,nindx2
287#include "lockon.inc"
288 WRITE(iout, 2000) ngl(indx2(j)),ipg
289 WRITE(istdo,2100) ngl(indx2(j)),ipg,tt
290#include "lockoff.inc"
291 ENDDO
292 ENDIF
293
294 1000 FORMAT(1x,'MINIMUM TIMESTEP (PLAS_RATE) REACHED, DELETED SOLID ELEMENT ',i10)
295 1100 FORMAT(1x,'MINIMUM TIMESTEP (PLAS_RATE) REACHED, DELETED SOLID ELEMENT ',i10,1x,'AT TIME :',1pe12.4)
296 2000 FORMAT(1x,'FAILURE (PLAS_RATE) IN SOLID ELEMENT ',i10,1x,',GAUSS PT',i2)
297 2100 FORMAT(1x,'FAILURE (PLAS_RATE) IN SOLID ELEMENT ',i10,1x,',GAUSS PT',i2,1x,'AT TIME :',1pe12.4)
298
299
300 RETURN
subroutine mat121_newton(nel, ngl, nuparam, nuvar, nfunc, ifunc, npf, tf, timestep, time, uparam, uvar, rho, pla, dpla, soundsp, epsd, off, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigy, et, seq)
subroutine mat121_nice(nel, ngl, nuparam, nuvar, nfunc, ifunc, npf, tf, timestep, time, uparam, uvar, rho, pla, dpla, soundsp, epsd, off, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigy, et, seq)
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)