51
52
53
54 USE redef3_mod
57 USE python_funct_mod
58 use element_mod , only : nixr
59
60
61
62#include "implicit_f.inc"
63#include "comlock.inc"
64
65
66
67#include "mvsiz_p.inc"
68
69
70
71#include "units_c.inc"
72#include "scr17_c.inc"
73#include "param_c.inc"
74#include "com04_c.inc"
75#include "com08_c.inc"
76#include "scr14_c.inc"
77#include "com01_c.inc"
78#include "impl1_c.inc"
79
80
81
82 type(python_), intent(inout) :: PYTHON
83 INTEGER, INTENT(IN) :: SANIN
84 INTEGER, INTENT(IN) :: STF
85 INTEGER, INTENT(IN) :: NFT
86 INTEGER, INTENT(IN) :: NEL
87 INTEGER, INTENT(IN) :: IRESP
88 INTEGER, INTENT(IN) :: SNPC
89 INTEGER NPF(SNPC),IXR(NIXR,*),IGEO(NPROPGI,*),NGL(*),MGN(*),
90 . NC1(*),NC2(*),NC3(*),NUVAR
91
93 . geo(npropg,*), f(*), al0(*), e(*), dl(*), tf(stf), off(*),
94 . dpl(*), fep(*), dpx2(*),dfs(*),v(3,*),df(*),anim(*),
95 . ipos(*),al0_err(*),yield(*),inifric(*),
96 . ex(mvsiz),ey(mvsiz),ez(mvsiz),xk(mvsiz),
97 . xm(mvsiz),xc(mvsiz),ak(mvsiz),ex2(mvsiz),ey2(mvsiz),ez2(mvsiz),
98 . uvar(nuvar,*),dl0(*)
99 DOUBLE PRECISION X1DP(3,*),X2DP(3,*),X3DP(3,*)
100 TYPE(TTABLE) TABLE(*)
101 TARGET :: uvar
102
103
104
105 INTEGER IFUNC2(MVSIZ),
106 . IECROU(MVSIZ), IFUNC(MVSIZ), IFV(MVSIZ),
107 . NINDX, INDX(MVSIZ),I,ILENG,IPID,J,
108 . IFUNC3(MVSIZ),ITABF(MVSIZ),IFRIC(MVSIZ)
109
111 . xl0(mvsiz),dmn(mvsiz),dmx(mvsiz),
112 . dlold(mvsiz),b(mvsiz), d(mvsiz), dv(mvsiz),
113 . ff(mvsiz),lscale(mvsiz),ee(mvsiz),
114 . gf3(mvsiz),fric(mvsiz),xscalf(mvsiz),yscalf(mvsiz),
115 . f_min(mvsiz),f_max(mvsiz),epla(mvsiz)
117 . beta, mu,fmax, ddx ,vx21,vy21,vz21,vl21,
118 . vx32,vy32,vz32,vl32,not_used,not_used2(2)
120 . DIMENSION(1) :: xx
121 DOUBLE PRECISION EXDP(MVSIZ) ,EYDP(MVSIZ) ,EZDP(MVSIZ),
122 . EX2DP(MVSIZ),EY2DP(MVSIZ),EZ2DP(MVSIZ),
123 . AL1DP(MVSIZ),AL2DP(MVSIZ),ALDP(MVSIZ),
124 . AL0DP(MVSIZ)
126 my_real ,
DIMENSION(:),
POINTER :: xx_old
127 TARGET :: not_used2
128
129
130 not_used = zero
131 not_used2 = zero
132
133 DO i=1,nel
134 epla(i) = zero
135 ipid = mgn(i)
136 xm(i) =geo(1,ipid)
137 xk(i) =geo(2,ipid)
138 xc(i) =geo(3,ipid)
139 iecrou(i)=nint(geo(7,ipid))
140 ak(i) = geo(10,ipid)
141 b(i) = geo(11,ipid)
142 d(i) = geo(13,ipid)
143 ee(i) = geo(40,ipid)
144 gf3(i) = geo(132,ipid)
145 ff(i) = geo(18,ipid)
146 dmn(i) = geo(15,ipid)
147 dmx(i) = geo(16,ipid)
148 fric(i) = geo(17,ipid)
149 xscalf(i)= geo(20,ipid)
150 lscale(i)= geo(39,ipid)
151 ifunc(i) = igeo(101,ipid)
152 ifv(i) = igeo(102,ipid)
153 ifunc2(i)= igeo(103,ipid)
154 ifunc3(i)= igeo(119,ipid)
155 itabf(i) = igeo(201,ipid)
156 ifric(i) = igeo(202,ipid)
157 f_min(i) = geo(138,ipid)
158 f_max(i) = geo(139,ipid)
159 yscalf(i)= geo(140,ipid)
160 max_slope(i) = geo(141,ipid)
161 ENDDO
162
163 DO i=1,nel
164 exdp(i)=x2dp(1,i)-x1dp(1,i)
165 eydp(i)=x2dp(2,i)-x1dp(2,i)
166 ezdp(i)=x2dp(3,i)-x1dp(3,i)
167 ex2dp(i)=x2dp(1,i)-x3dp(1,i)
168 ey2dp(i)=x2dp(2,i)-x3dp(2,i)
169 ez2dp(i)=x2dp(3,i)-x3dp(3,i)
170 dlold(i)=dl(i)
171 ENDDO
172
173 IF (inispri /= 0 .and. tt == zero) THEN
174 DO i=1,nel
175 dlold(i)=dl0(i)
176 ENDDO
177 ENDIF
178
179 DO i=1,nel
180 al1dp(i)=sqrt(exdp(i)*exdp(i)+eydp(i)*eydp(i)+ezdp(i)*ezdp(i))
181 al2dp(i)=sqrt(ex2dp(i)*ex2dp(i)+ey2dp(i)*ey2dp(i)+
182 . ez2dp(i)*ez2dp(i))
183 aldp(i) =al1dp(i) + al2dp(i)
184 al1dp(i)=
max(al1dp(i),em15)
185 exdp(i) = exdp(i)/al1dp(i)
186 eydp(i) = eydp(i)/al1dp(i)
187 ezdp(i) = ezdp(i)/al1dp(i)
188 al2dp(i)=
max(al2dp(i),em15)
189 ex2dp(i)= ex2dp(i)/al2dp(i)
190 ey2dp(i)= ey2dp(i)/al2dp(i)
191 ez2dp(i)= ez2dp(i)/al2dp(i)
192 ENDDO
193
194 DO i=1,nel
195 ex(i) = exdp(i)
196 ey(i) = eydp(i)
197 ez(i) = ezdp(i)
198 ex2(i)= ex2dp(i)
199 ey2(i)= ey2dp(i)
200 ez2(i)= ez2dp(i)
201 ENDDO
202
203 IF (inispri /= 0 .and. tt == zero) THEN
204 DO i=1,nel
205 xl0(i)= al0(i)
206
207 IF (xl0(i) == zero) xl0(i) = aldp(i)
208 ENDDO
209 ENDIF
210
211 IF (tt == zero) THEN
212 DO i=1,nel
213 al0(i)= aldp(i)
214 ENDDO
215 ENDIF
216
217 IF (scodver >= 101) THEN
218 IF (tt == zero) THEN
219 DO i=1,nel
220 al0_err(i)=aldp(i)-al0(i)
221 ENDDO
222 ENDIF
223 ENDIF
224
225 IF ( inispri /= 0 .and. tt == zero) THEN
226 DO i=1,nel
227 al0(i)= xl0(i)
228 ENDDO
229 ENDIF
230
231 DO i=1,nel
232 al0dp(i)= al0(i)
233 ENDDO
234
235 IF (scodver >= 101) THEN
236 DO i=1,nel
237 al0dp(i)= al0dp(i)+al0_err(i)
238 ENDDO
239 ENDIF
240
241 IF (ismdisp > 0) THEN
242 DO i=1,nel
243 vx21= v(1,nc2(i)) - v(1,nc1(i))
244 vy21= v(2,nc2(i)) - v(2,nc1(i))
245 vz21= v(3,nc2(i)) - v(3,nc1(i))
246 vl21 = vx21*ex(i)+vy21*ey(i)+vz21*ez(i)
247 vx32= v(1,nc2(i)) - v(1,nc3(i))
248 vy32= v(2,nc2(i)) - v(2,nc3(i))
249 vz32= v(3,nc2(i)) - v(3,nc3(i))
250 vl32 = vx32*ex2(i)+vy32*ey2(i)+vz32*ez2(i)
251 dl(i)= dlold(i)+(vl21+vl32)*dt1
252 ENDDO
253 ELSE
254 DO i=1,nel
255 dl(i)= (aldp(i)-al0dp(i))
256 ENDDO
257 ENDIF
258
259 DO i=1,nel
260 ileng=nint(geo(93,mgn(i)))
261 IF (ileng /= 0) THEN
262 xl0(i)=al0dp(i)
263 ELSE
264 xl0(i) = one
265 ENDIF
266 ENDDO
267
268 IF (nuvar > 0) THEN
269 xx_old => uvar(1,1:nel)
270 ELSE
271 xx_old => not_used2
272 ENDIF
273 CALL redef3(python,
274 1 f, xk, dl, fep,
275 2 dlold, dpl, tf, npf,
276 3 xc, off, e, dpx2,
277 4 anim, anim_fe(11),ipos,
278 5 xl0, dmn, dmx, dv,
279 6 ff, lscale, ee, gf3,
280 7 ifunc3, yield, aldp, ak,
281 8 b, d, iecrou, ifunc,
282 9 ifv, ifunc2, epla, xx_old,
283 a nel, nft, stf, sanin, dt1,
284 b iresp, impl_s, idyna, snpc,
285 c max_slope=max_slope)
286 nindx = 0
287 DO i=1,nel
288 IF (off(i) == one .AND. dmx(i) /= zero .AND. dmn(i) /= zero) THEN
289 IF (dl(i) > dmx(i) .OR. dl(i) < dmn(i)) THEN
290 off(i)=zero
291 nindx = nindx + 1
292 indx(nindx) = i
293 idel7nok = 1
294 ENDIF
295 ENDIF
296 ENDDO
297
298 DO j=1,nindx
299 i = indx(j)
300#include "lockon.inc"
301 WRITE(iout, 1000) ngl(i)
302 WRITE(istdo,1100) ngl(i),tt
303#include "lockoff.inc"
304 ENDDO
305
306 DO i=1,nel
307 xm(i)=xm(i)*xl0(i)
308 xk(i)=xk(i)/xl0(i)
309
310
311
312 xc(i) = (xc(i) + max_slope(i)) / xl0(i)
313 ENDDO
314
315 DO i=1,nel
316 IF (ifric(i) == 0) THEN
317 IF (itabf(i) > zero) THEN
318 ddx = exdp(i) * (v(1,nc2(i)) - v(1,nc1(i)))
319 . + eydp(i) * (v(2,nc2(i)) - v(2,nc1(i)))
320 . + ezdp(i) * (v(3,nc2(i)) - v(3,nc1(i)))
321 . + ex2dp(i)* (v(1,nc3(i)) - v(1,nc2(i)))
322 . + ey2dp(i)* (v(2,nc3(i)) - v(2,nc2(i)))
323 . + ez2dp(i)* (v(3,nc3(i)) - v(3,nc2(i)))
324 dfs(i)= dfs(i) + ddx * dt1 * xk(i)
325 beta = pi - acos(exdp(i)*ex2dp(i)+eydp(i)*ey2dp(i)+
326 . ezdp(i)*ez2dp(i))
327 xx(1) = abs(dfs(i)*xscalf(i))
329 mu = mu*yscalf(i)
330 fmax =
max(zero,f(i) * tanh(half*mu*beta))
331 fmax =
min(fmax,abs(dfs(i)))
332 dfs(i)= sign(fmax,dfs(i))
333 df(i) = dfs(i)
334 ELSEIF (fric(i) > zero) THEN
335 ddx = exdp(i) * (v(1,nc2(i)) - v(1,nc1(i)))
336 . + eydp(i) * (v(2,nc2(i)) - v(2,nc1(i)))
337 . + ezdp(i) * (v(3,nc2(i)) - v(3,nc1(i)))
338 . + ex2dp(i)* (v(1,nc3(i)) - v(1,nc2(i)))
339 . + ey2dp(i)* (v(2,nc3(i)) - v(2,nc2(i)))
340 . + ez2dp(i)* (v(3,nc3(i)) - v(3,nc2(i)))
341 dfs(i)= dfs(i) + ddx * dt1 * xk(i)
342 beta = pi - acos(exdp(i)*ex2dp(i)+eydp(i)*ey2dp(i)+
343 . ezdp(i)*ez2dp(i))
344 fmax =
max(zero,f(i) * tanh(half*fric(i)*beta))
345 fmax =
min(fmax,abs(dfs(i)))
346 dfs(i)= sign(fmax,dfs(i))
347 df(i) = dfs(i)
348
349
350 ELSE
351 df(i)= zero
352 ENDIF
353 ELSEIF (ifric(i) > 0) THEN
354 IF (itabf(i) > zero) THEN
355 ddx = exdp(i) * (v(1,nc2(i)) - v(1,nc1(i)))
356 . + eydp(i) * (v(2,nc2(i)) - v(2,nc1(i)))
357 . + ezdp(i) * (v(3,nc2(i)) - v(3,nc1(i)))
358 . + ex2dp(i)* (v(1,nc3(i)) - v(1,nc2(i)))
359 . + ey2dp(i)* (v(2,nc3(i)) - v(2,nc2(i)))
360 . + ez2dp(i)* (v(3,nc3(i)) - v(3,nc2(i)))
361 dfs(i)= dfs(i) + ddx * dt1 * xk(i)
362 beta = pi - acos(exdp(i)*ex2dp(i)+eydp(i)*ey2dp(i)+
363 . ezdp(i)*ez2dp(i))
364
365 xx(1) = dfs(i)*xscalf(i)
367 mu = mu*yscalf(i)
368
369
370
371
372
373
374
375 IF (dfs(i) <= f_min(i) .OR. dfs(i) >= f_max(i)) inifric(i) = one
376 IF (inifric(i) == one) mu = fric(i)
377
378 fmax =
max(zero,f(i) * tanh(half*mu*beta))
379 fmax =
min(fmax,abs(dfs(i)))
380 dfs(i)= sign(fmax,dfs(i))
381 df(i) = dfs(i)
382 ELSEIF (fric(i) > zero .AND. itabf(i) == 0) THEN
383 ddx = exdp(i) * (v(1,nc2(i)) - v(1,nc1(i)))
384 . + eydp(i) * (v(2,nc2(i)) - v(2,nc1(i)))
385 . + ezdp(i) * (v(3,nc2(i)) - v(3,nc1(i)))
386 . + ex2dp(i)* (v(1,nc3(i)) - v(1,nc2(i)))
387 . + ey2dp(i)* (v(2,nc3(i)) - v(2,nc2(i)))
388 . + ez2dp(i)* (v(3,nc3(i)) - v(3,nc2(i)))
389 dfs(i)= dfs(i) + ddx * dt1 * xk(i)
390 beta = pi - acos(exdp(i)*ex2dp(i)+eydp(i)*ey2dp(i)+
391 . ezdp(i)*ez2dp(i))
392 fmax =
max(zero,f(i) * tanh(half*fric(i)*beta))
393 fmax =
min(fmax,abs(dfs(i)))
394 dfs(i)= sign(fmax,dfs(i))
395 df(i) = dfs(i)
396 ELSE
397 df(i)= zero
398 ENDIF
399 ENDIF
400 ENDDO
401
402 1000 FORMAT(1x,'-- RUPTURE OF SPRING ELEMENT NUMBER ',i10)
403 1100 FORMAT(1x,'-- RUPTURE OF SPRING ELEMENT :',i10,' AT TIME :',g11.4)
404
405 RETURN