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