OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
redef3_law113.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23C----s---1----+----2----+----3----+----4----+----5----+----6----+----7--
24!||====================================================================
25!|| redef3_law113 ../engine/source/elements/spring/redef3_law113.F
26!||--- calls -----------------------------------------------------
27!|| vinter2 ../engine/source/tools/curve/vinter.F
28!|| vinter2dp ../engine/source/tools/curve/vinter.F
29!|| vinter_mixed ../engine/source/tools/curve/vinter_mixed.F90
30!||--- uses -----------------------------------------------------
31!|| python_funct_mod ../common_source/modules/python_mod.F90
32!|| vinter_mixed_mod ../engine/source/tools/curve/vinter_mixed.F90
33!||====================================================================
34 SUBROUTINE redef3_law113(PYTHON,NEL ,NFT,
35 . FX ,XK ,DX ,FXEP ,DXOLD ,
36 . DPX ,TF ,NPF ,XC ,OFF ,
37 . E ,DPX2 ,ANIM ,IANI ,POS ,
38 . XL0 ,DMN ,DMX ,DVX ,
39 . FF ,LSCALE ,EE ,GF3 ,IFUNC3,
40 . YIELD ,AK ,B ,D ,
41 . IECROU ,IFUNC ,IFV ,IFUNC2 ,IFUNC4 ,
42 . EPLA ,XX_OLD )
43C-----------------------------------------------
44 USE python_funct_mod
45 USE vinter_mixed_mod
46C-----------------------------------------------
47C I m p l i c i t T y p e s
48C-----------------------------------------------
49#include "implicit_f.inc"
50C-----------------------------------------------
51C G l o b a l P a r a m e t e r s
52C-----------------------------------------------
53#include "mvsiz_p.inc"
54C-----------------------------------------------
55C C o m m o n B l o c k s
56C-----------------------------------------------
57#include "com08_c.inc"
58#include "scr05_c.inc"
59#include "impl1_c.inc"
60#include "com04_c.inc"
61#include "tabsiz_c.inc"
62C-----------------------------------------------
63C D u m m y A r g u m e n t
64C-----------------------------------------------
65 type(python_), intent(inout) :: PYTHON
66 INTEGER, INTENT(IN) :: NPF(SNPC),IANI,NEL, NFT
67 INTEGER, DIMENSION(MVSIZ),INTENT(INOUT) :: IECROU,IFUNC,IFV,
68 . IFUNC2,IFUNC3,IFUNC4
69C REAL
70 my_real, DIMENSION(MVSIZ), INTENT(IN):: AK,B,XL0,XK,FF,LSCALE,DMN,
71 . DMX,OFF,XC,EE,D
72 my_real, DIMENSION(MVSIZ), INTENT(INOUT) :: EPLA,DXOLD
73 my_real, DIMENSION(MVSIZ), INTENT(OUT) :: dvx
74 my_real, DIMENSION(NEL) , INTENT(INOUT) :: fx,dx,dpx,fxep,
75 . e,dpx2, yield,xx_old
76 my_real, INTENT(INOUT) :: pos(6,nel), anim(numelr*iani)
77 my_real, INTENT(IN) :: tf(stf)
78C-----------------------------------------------
79C L o c a l V a r i a b l e s
80C-----------------------------------------------
81 INTEGER, DIMENSION(MVSIZ) :: JPOS , JLEN,JAD,JPOS2,
82 . jlen2,jpos3, jlen3,jad2,
83 . ic1, jlen4,ic2,ic0,
84 . j2pos,j2len,
85 . j2ad,jpos4, j3pos, j3len,j3ad
86
87 INTEGER
88 . JFUNC,JFUNC2,JDMP,JECROU(-1:11),J2DMP,K1,NP1,NP2,
89 . i, j, ii, interp, k, func,fund,j1,j2func,j3func,j2k
90c REAL ou REAL*8
91 my_real, DIMENSION(MVSIZ) :: ddx,fold, gx, dxela, dyd,
92 . xx, xx2, xx3, yy, yy2, yy3,dydx,
93 . dydx2, dydx3, dydxv, dperm,fmax,
94 . dvxs,gf3,dydxv2,fmin,gx2,fy0,
95 . xfy0,ycut,xcut,xn3fy0,
96 . ul,fyield,dpl,yy22,xfyn1,an3y0,ddpx,
97 . xfyn33,xfyn11,dpx3,dpx22,xx22,
98 . edecal2,edecal3,ecrou,kx2,dkdx2,
99 . tmp1, tmp2
100 my_real
101 . dvv, dfac, dt11, damp, damm,b1,xi1,xi2,yi1,yi2 ,
102 . s1,s2,t1c,t2c,x1,x2,y1,y2
103 INTEGER :: PYID
104 LOGICAL :: ANY_PYTHON_FUNCTION
105C=======================================================================
106 ANY_PYTHON_FUNCTION = .false.
107 dt11 = dt1
108 IF(dt11==zero)dt11 = ep30
109 DO i=1,nel
110 dx(i)=dx(i)/xl0(i)
111 dxold(i)=dxold(i)/xl0(i)
112 dpx(i)=dpx(i)/xl0(i)
113 dpx2(i)=dpx2(i)/xl0(i)
114 e(i)=e(i)/xl0(i)
115 ENDDO
116C
117 DO 80 i=1,nel
118 fold(i)=fx(i)
119 ddx(i)= (dx(i)-dxold(i))
120 dvx(i)= ddx(i)/ dt11
121 dvxs(i)= dvx(i)*ff(i)
122 80 CONTINUE
123C
124C
125 IF(iani/=0)THEN
126 DO i=1,nel
127 ii=i + nft
128 damp=dx(i)/max(dmx(i),em15)
129 damm=dx(i)/min(dmn(i),-em15)
130 anim(ii)=max(anim(ii),damp,damm)
131 anim(ii)=min(anim(ii),one)
132 ENDDO
133 ENDIF
134C-------------------------------------
135C VECTOR INTERPOLATION (ADRESS)
136C-------------------------------------
137 jecrou(-1) = 0
138 jecrou(0) = 0
139 jecrou(1) = 0
140 jecrou(2) = 0
141 jecrou(3) = 0
142 jecrou(4) = 0
143 jecrou(5) = 0
144 jecrou(6) = 0
145 jecrou(7) = 0
146 jecrou(8) = 0
147 jecrou(9) = 0
148 jecrou(10) = 0
149 jecrou(11) = 0
150 interp = 0
151 jdmp = 0
152 j2dmp = 0
153 j2k = 0
154C
155 DO i=1,nel
156 IF(ifunc(i)==0)THEN ! ifunc =IGEO(101)-FCT_id1
157 jecrou(-1) = jecrou(-1) + 1
158c modif pour vectorisation
159 ELSEIF(iecrou(i)==0)THEN
160 jecrou(0) = jecrou(0) + 1
161 interp = 1
162 ELSEIF(iecrou(i)==1)THEN
163 jecrou(1) = jecrou(1) + 1
164 interp = 1
165 ELSEIF(iecrou(i)==2)THEN
166 jecrou(2) = jecrou(2) + 1
167 interp = 1
168 ELSEIF(iecrou(i)==3)THEN
169 jecrou(3) = jecrou(3) + 1
170 interp = 1
171 ELSEIF(iecrou(i)==4)THEN
172 jecrou(4) = jecrou(4) + 1
173 interp = 1
174 ELSEIF(iecrou(i)==5)THEN
175 jecrou(5) = jecrou(5) + 1
176 interp = 1
177 ELSEIF(iecrou(i)==6)THEN
178 jecrou(6) = jecrou(6) + 1
179 interp = 1
180 ELSEIF(iecrou(i)==7)THEN
181 jecrou(7) = jecrou(7) + 1
182 interp = 1
183 ENDIF
184 IF(ifv(i)/=0) jdmp = jdmp + 1
185 IF(ifunc3(i)/=0) j2dmp = j2dmp + 1
186 IF(ifunc4(i)/=0) j2k = j2k + 1
187 ENDDO
188C
189 IF(interp > 0)THEN
190 DO i=1,nel
191 jpos(i) = nint(pos(1,i))
192 jpos2(i) = nint(pos(2,i))
193 jpos3(i) = nint(pos(3,i))
194 jfunc =max(1,ifunc(i))
195 pyid = python_funct_id(nfunct,jfunc,npf)
196 IF(pyid <=0) THEN ! not python function
197 jfunc =max(1,ifunc(i))
198 jad(i) = npf(jfunc) / 2 + 1
199 jlen(i) = npf(jfunc+1) / 2 - jad(i) - jpos(i)
200 jlen3(i) = npf(jfunc+1) / 2 - jad(i) - jpos3(i)
201 ELSE
202 any_python_function = .true.
203 jad(i) = -pyid
204 jlen(i) = -pyid
205 jlen3(i) = -pyid
206 ENDIF
207 jfunc2=max(1,ifunc2(i))
208 pyid = python_funct_id(nfunct, jfunc2,npf)
209 IF(pyid <=0) THEN ! not python function
210 jfunc2=max(1,ifunc2(i))
211 jad2(i) = npf(jfunc2) / 2 + 1
212 jlen2(i) = npf(jfunc2+1) / 2 - jad2(i) - jpos2(i)
213 ELSE
214 jad2(i) = -pyid
215 jlen2(i) = -pyid
216 any_python_function = .true.
217 ENDIF
218 xx(i) =zero
219 xx2(i)=zero
220 xx3(i)=zero
221 ENDDO
222 ENDIF
223C-------------------------------------
224C NON LINEAR ELASTIC
225C NL ELASTO PLASTIC (UNCOUPLED TENSION COMPRESSION )
226C-------------------------------------
227 IF(jecrou(0)+jecrou(2)== nel)THEN
228 DO i=1,nel
229 xx(i)=dx(i)
230 ENDDO
231 ELSEIF(jecrou(0)+jecrou(2)>0)THEN
232 DO i=1,nel
233 IF(ifunc(i)/=0.AND.(iecrou(i)==0.OR.iecrou(i)==2))THEN
234 xx(i)=dx(i)
235 ENDIF
236 ENDDO
237 ENDIF
238C-------------------------------------
239C ELASTO-PLASTIC (ISOTROP)
240C ELASTO PLASTIC (6DOF COUPLED TENSION COMPRESSION )
241C-------------------------------------
242 IF(jecrou(1)+jecrou(3)== nel )THEN
243 DO i=1,nel
244 fx(i)=fxep(i)+xk(i)*ddx(i)
245 IF(fx(i)>=0.)THEN
246 xx(i)=dpx(i)+fx(i)/xk(i)
247 ELSE
248 xx(i)=-dpx(i)+fx(i)/xk(i)
249 ENDIF
250 ENDDO
251 ELSEIF(jecrou(1)+jecrou(3)>0)THEN
252 DO i=1,nel
253 IF(ifunc(i)/=0.AND.(iecrou(i)==1.OR.iecrou(i)==3))THEN
254 fx(i)=fxep(i)+xk(i)*ddx(i)
255 IF(fx(i)>=zero)THEN
256 xx(i)=dpx(i)+fx(i)/xk(i)
257 ELSE
258 xx(i)=-dpx(i)+fx(i)/xk(i)
259 ENDIF
260 ENDIF
261 ENDDO
262 ENDIF
263C-------------------------------------
264C ELASTO PLASTIC (KINEMATIC HARDENING )
265C-------------------------------------
266 IF(jecrou(4)== nel )THEN
267 DO i=1,nel
268 interp = max(2,interp)
269 xx(i) =dx(i)
270 xx2(i)=dx(i)
271 ENDDO
272 ELSEIF(jecrou(4)>0)THEN
273 DO i=1,nel
274 IF(ifunc(i)/=0.AND.iecrou(i)==4)THEN
275 interp = max(2,interp)
276 xx(i) =dx(i)
277 xx2(i)=dx(i)
278 ENDIF
279 ENDDO
280 ENDIF
281C-------------------------------------
282C ELASTO PLASTIC (UNCOUPLED TENSION COMPRESSION )
283C (D/R) NON LINEAR RELOADING
284C DPX = MAXIMAL DISPLACEMENT
285C-------------------------------------
286 IF(jecrou(5)== nel )THEN
287 DO i=1,nel
288 xx(i)=dx(i)
289 IF(dx(i)>zero)THEN
290 interp = max(3,interp)
291 xx2(i)=dpx(i)
292 xx3(i)=dpx(i)
293 ELSEIF(dx(i)<zero)THEN
294 interp = max(3,interp)
295 xx2(i)=dpx2(i)
296 xx3(i)=dpx2(i)
297 ELSE
298 interp = max(1,interp)
299 ENDIF
300 ENDDO
301 ELSEIF(jecrou(5)>0)THEN
302 DO i=1,nel
303 IF(ifunc(i)/=0.AND.iecrou(i)==5)THEN
304 xx(i)=dx(i)
305 IF(dx(i)>zero)THEN
306 interp = max(3,interp)
307 xx2(i)=dpx(i)
308 xx3(i)=dpx(i)
309 ELSEIF(dx(i)<zero)THEN
310 interp = max(3,interp)
311 xx2(i)=dpx2(i)
312 xx3(i)=dpx2(i)
313 ELSE
314 interp = max(1,interp)
315 ENDIF
316 ENDIF
317 ENDDO
318 ENDIF
319C-------------------------------------
320C ELASTO PLASTIC (ISTROPIC HARDENING )
321C-------------------------------------
322 IF(any_python_function .AND. jecrou(6) > 0) THEN
323 DO i=1,nel
324 IF(ifunc2(i) < 0 ) THEN
325 fund = -jlen2(i) ! Curve ID
326 ! Obtain the derivative directly at FXEP(I)
327 CALL python_deriv_funct1d(python, fund, fxep(i), an3y0(i))
328 ! Calculate FX(I) and XX(I) based on the derivative
329 fx(i) = fxep(i) + an3y0(i) * ddx(i)
330 xx(i) = sign(abs(xx_old(i)), fx(i))
331 xx(i) = xx(i) + ddx(i)
332 interp = 0
333 ENDIF
334 ENDDO
335 ENDIF
336
337 IF(jecrou(6)== nel)THEN
338 DO i=1,nel
339 fund = ifunc2(i) ! courbe N3 de unload
340 IF(jlen2(i) < 0) cycle !if JLEN2(I) < 0, FUNC2 is a python function
341 np2 = (npf(fund+1)-npf(fund))/2
342 an3y0(i)= zero
343 DO k=2,np2
344 k1=2*(k-2)
345 x1=tf(npf(fund)+k1)
346 x2=tf(npf(fund)+k1+2)
347 y1=tf(npf(fund)+k1+1)
348 y2=tf(npf(fund)+k1+3)
349 IF((fxep(i)< y2.AND.fxep(i)>=y1))THEN
350 an3y0(i)=(y2-y1)/ (x2-x1)
351 xn3fy0(i)=(fxep(i)-y1)/an3y0(i) + x1 !ABS DE N3
352 EXIT
353 ENDIF
354 ENDDO
355 IF (an3y0(i)== zero)THEN ! extrapolation (exterieur aux points de l input)
356 x1=tf(npf(fund)+(np2-2)*2)
357 x2=tf(npf(fund)+(np2-2)*2+2)
358 y1=tf(npf(fund)+(np2-2)*2+1)
359 y2=tf(npf(fund)+(np2-2)*2+3)
360C
361 xi1=tf(npf(fund))
362 xi2=tf(npf(fund)+2)
363 yi1=tf(npf(fund)+1)
364 yi2=tf(npf(fund)+3)
365 IF(fxep(i)>y2)an3y0(i)=(y2-y1)/ (x2-x1)
366 IF(fxep(i)<yi1)an3y0(i)=(yi2-yi1)/ (xi2-xi1)
367 ENDIF
368 fx(i)=fxep(i)+an3y0(i)*ddx(i)
369 xx(i)=sign(abs(xx_old(i)),fx(i))
370 xx(i)=xx(i)+ddx(i)
371 interp=0
372 ENDDO
373 ELSEIF(jecrou(6)>0)THEN
374 DO i=1,nel
375 IF(jlen(i) < 0) cycle ! python function
376 IF(ifunc(i) /= 0.AND.iecrou(i)== 6)THEN
377 fund = ifunc2(i)
378 np2 = (npf(fund+1)-npf(fund))/2
379 an3y0(i)= zero
380 DO k=2,np2
381 k1=2*(k-2)
382 x1=tf(npf(fund)+k1)
383 x2=tf(npf(fund)+k1+2)
384 y1=tf(npf(fund)+k1+1)
385 y2=tf(npf(fund)+k1+3)
386 IF((fxep(i)< y2.AND.fxep(i)>=y1))THEN
387 an3y0(i)=(y2-y1)/ (x2-x1)
388 xn3fy0(i)=(fxep(i)-y1)/an3y0(i) + x1 !ABS DE N3
389 EXIT
390 ENDIF
391 ENDDO
392 IF (an3y0(i)== zero)THEN
393 x1=tf(npf(fund)+(np2-2)*2)
394 x2=tf(npf(fund)+(np2-2)*2+2)
395 y1=tf(npf(fund)+(np2-2)*2+1)
396 y2=tf(npf(fund)+(np2-2)*2+3)
397C
398 xi1=tf(npf(fund))
399 xi2=tf(npf(fund)+2)
400 yi1=tf(npf(fund)+1)
401 yi2=tf(npf(fund)+3)
402 IF(fxep(i)>y2)an3y0(i)=(y2-y1)/ (x2-x1)
403 IF(fxep(i)<yi1)an3y0(i)=(yi2-yi1)/ (xi2-xi1)
404 ENDIF
405 fx(i)=fxep(i)+an3y0(i)*ddx(i)
406 xx(i)=sign(abs(xx_old(i)),fx(i))
407 xx(i)=xx(i)+ddx(i)
408 ENDIF
409 ENDDO
410 ENDIF
411c-------------------------------------
412c ELASTO PLASTIC TWO CURVES FOR LOAD AND UNLOAD
413c-------------------------------------
414 IF(jecrou(7)== nel)THEN
415 DO i=1,nel
416 interp = max(2,interp)
417 xx(i) =dx(i)
418 xx2(i)=dx(i)
419 ENDDO
420 ELSEIF(jecrou(7)>0)THEN
421 DO i=1,nel
422 IF(ifunc(i)/=0.AND.iecrou(i)==7)THEN
423 interp = max(2,interp)
424 xx(i) =dx(i)
425 xx2(i)=dx(i)
426 ENDIF
427 ENDDO
428 ENDIF
429C-------------------------------------
430c VECTOR INTERPOLATION
431C-------------------------------------
432 DO i=1,nel
433 xx(i) = xx(i) *lscale(i)
434 xx2(i) = xx2(i)*lscale(i)
435 xx3(i) = xx3(i)*lscale(i)
436 ENDDO
437C----s---1----+----2----+----3----+----4----+----5----+----6----+----7--
438 IF(any_python_function) THEN
439 IF(interp>=1) CALL vinter_mixed(python, tf,jad,jpos,jlen,nel,xx,dydx,yy)
440 IF(interp>=2) CALL vinter_mixed(python, tf,jad2,jpos2,jlen2,nel,xx2,dydx2,yy2)
441 IF(interp>=3) CALL vinter_mixed(python, tf,jad ,jpos3,jlen3,nel,xx3,dydx3,yy3)
442 ELSEIF (iresp==1) THEN
443 IF(interp>=1)
444 . CALL vinter2dp(tf,jad ,jpos ,jlen ,nel,xx ,dydx ,yy)
445 IF(interp>=2)
446 . CALL vinter2dp(tf,jad2,jpos2,jlen2,nel,xx2,dydx2,yy2)
447 IF(interp>=3)
448 . CALL vinter2dp(tf,jad ,jpos3,jlen3,nel,xx3,dydx3,yy3)
449 ELSE
450 IF(interp>=1)
451 . CALL vinter2(tf,jad ,jpos ,jlen ,nel,xx ,dydx ,yy )
452 IF(interp>=2)
453 . CALL vinter2(tf,jad2,jpos2,jlen2,nel,xx2,dydx2,yy2)
454 IF(interp>=3)
455 . CALL vinter2(tf,jad ,jpos3,jlen3,nel,xx3,dydx3,yy3)
456 ENDIF
457 IF(interp>0)THEN
458 DO i=1,nel
459 pos(1,i) = jpos(i)
460 pos(2,i) = jpos2(i)
461 pos(3,i) = jpos3(i)
462 ENDDO
463 ENDIF
464C-------------------------------------
465C LINEAR ELASTIC
466C-------------------------------------
467 IF(jecrou(-1) == nel)THEN
468 DO i=1,nel
469 fx(i)=xk(i)*dx(i)
470 ENDDO
471 ELSEIF(jecrou(-1)>0)THEN
472 DO i=1,nel
473 IF(ifunc(i)==0)THEN
474 fx(i)=xk(i)*dx(i)
475 ENDIF
476 ENDDO
477 ENDIF
478C-------------------------------------
479C ELASTIC F = f(total length)
480C-------------------------------------
481 IF (jecrou(8) == nel )THEN
482 DO i=1,nel
483 fx(i)=yy(i)
484 ENDDO
485 ELSEIF (jecrou(8) > 0)THEN
486 DO i=1,nel
487 IF (iecrou(i) == 8) fx(i)=yy(i)
488 ENDDO
489 ENDIF
490C-------------------------------------
491C NON LINEAR ELASTIC
492C-------------------------------------
493 IF(jecrou(0)== nel)THEN
494 DO i=1,nel
495 fx(i)=yy(i)
496 ENDDO
497 ELSEIF(jecrou(0)>0)THEN
498 DO i=1,nel
499 IF(ifunc(i)/=0.AND.iecrou(i)==0) fx(i)=yy(i)
500 ENDDO
501 ENDIF
502C-------------------------------------
503C ELASTO PLASTIC (ISOTROP)
504C-------------------------------------
505 IF(jecrou(1)== nel )THEN
506 DO i=1,nel
507 IF(fx(i)>=zero.AND.fx(i)>yy(i))THEN
508 dpx(i)=dpx(i)+(fx(i)-yy(i))/xk(i)
509 fx(i)=yy(i)
510 ELSEIF(fx(i)<zero.AND.fx(i)<yy(i))THEN
511 dpx(i)=dpx(i)+(yy(i)-fx(i))/xk(i)
512 fx(i)=yy(i)
513 ENDIF
514 fxep(i)=fx(i)
515 ENDDO
516 ELSEIF(jecrou(1)>0)THEN
517 DO i=1,nel
518 IF(ifunc(i)/=0.AND.iecrou(i)==1)THEN
519 IF(fx(i)>=zero.AND.fx(i)>yy(i))THEN
520 dpx(i)=dpx(i)+(fx(i)-yy(i))/xk(i)
521 fx(i)=yy(i)
522 ELSEIF(fx(i)<zero.AND.fx(i)<yy(i))THEN
523 dpx(i)=dpx(i)+(yy(i)-fx(i))/xk(i)
524 fx(i)=yy(i)
525 ENDIF
526 fxep(i)=fx(i)
527 ENDIF
528 ENDDO
529 ENDIF
530C-------------------------------------
531C ELASTO PLASTIC (UNCOUPLED TENSION COMPRESSION )
532C-------------------------------------
533 IF(jecrou(2)== nel )THEN
534 DO i=1,nel
535 IF(dx(i)>dpx(i))THEN
536 fx(i) = xk(i) * (dx(i)-dpx(i))
537 fxep(i)= yy(i)
538 fx(i) = min(fx(i),fxep(i))
539 dpx(i) = dx(i) - fx(i) / xk(i)
540 ELSEIF(dx(i)<dpx2(i))THEN
541 fx(i) = xk(i) * (dx(i)-dpx2(i))
542 fxep(i) = yy(i)
543 fx(i) = max(fx(i),fxep(i))
544 dpx2(i) = dx(i) - fx(i) / xk(i)
545 ELSE
546 fx(i) = zero
547 ENDIF
548 ENDDO
549 ELSEIF(jecrou(2)>0)THEN
550 DO i=1,nel
551 IF(ifunc(i)/=0.AND.iecrou(i)==2)THEN
552 IF(dx(i)>dpx(i))THEN
553 fx(i) = xk(i) * (dx(i)-dpx(i))
554 fxep(i)= yy(i)
555 fx(i) = min(fx(i),fxep(i))
556 dpx(i) = dx(i) - fx(i) / xk(i)
557 ELSEIF(dx(i)<dpx2(i))THEN
558 fx(i) = xk(i) * (dx(i)-dpx2(i))
559 fxep(i) = yy(i)
560 fx(i) = max(fx(i),fxep(i))
561 dpx2(i) = dx(i) - fx(i) / xk(i)
562 ELSE
563 fx(i) = zero
564 ENDIF
565 ENDIF
566 ENDDO
567 ENDIF
568C-------------------------------------
569C ELASTO PLASTIC (6DOF COUPLED TENSION COMPRESSION )
570C-------------------------------------
571 IF(jecrou(3)== nel )THEN
572 DO i=1,nel
573 IF(fx(i)>=zero.AND.fx(i)>yy(i))THEN
574 epla(i)=epla(i)+abs(yy(i)*(fx(i)-yy(i))/xk(i))
575 fx(i)=yy(i)
576 ELSEIF(fx(i)<zero.AND.fx(i)<yy(i))THEN
577 epla(i)=epla(i)+abs(yy(i)*(yy(i)-fx(i))/xk(i))
578 fx(i)=yy(i)
579 ENDIF
580 fxep(i)=fx(i)
581 ENDDO
582 ELSEIF(jecrou(3)>0)THEN
583 DO i=1,nel
584 IF(ifunc(i)/=0.AND.iecrou(i)==3)THEN
585 IF(fx(i)>=zero.AND.fx(i)>yy(i))THEN
586 epla(i)=epla(i)+abs(yy(i)*(fx(i)-yy(i))/xk(i))
587 fx(i)=yy(i)
588 ELSEIF(fx(i)<zero.AND.fx(i)<yy(i))THEN
589 epla(i)=epla(i)+abs(yy(i)*(yy(i)-fx(i))/xk(i))
590 fx(i)=yy(i)
591 ENDIF
592 fxep(i)=fx(i)
593 ENDIF
594 ENDDO
595 ENDIF
596C-------------------------------------
597C ELASTO-PLASTIC (KINEMATIC HARDENING)
598C-------------------------------------
599 IF(jecrou(4)== nel )THEN
600 DO i=1,nel
601 fx(i) = fxep(i) + xk(i)*ddx(i)
602 ic2(i)= 0
603 IF(fx(i)>yy(i))THEN
604 ic2(i)=1
605 fx(i) = yy(i)
606 ENDIF
607 IF(fx(i)<yy2(i))THEN
608 ic2(i)=2
609 fx(i) = yy2(i)
610 ENDIF
611 ENDDO
612 DO i=1,nel
613 fxep(i)=fx(i)
614 dpx(i) = dx(i) - fx(i) / xk(i)
615 ENDDO
616 ELSEIF(jecrou(4)>0)THEN
617 DO i=1,nel
618 IF(ifunc(i)/=0.AND.iecrou(i)==4)THEN
619 fx(i) = fxep(i) + xk(i)*ddx(i)
620 ic2(i)= 0
621 IF(fx(i)>yy(i))THEN
622 ic2(i)=1
623 fx(i) = yy(i)
624 ENDIF
625 IF(fx(i)<yy2(i))THEN
626 ic2(i)=2
627 fx(i) = yy2(i)
628 ENDIF
629 ENDIF
630 ENDDO
631 DO i=1,nel
632 IF(ifunc(i)/=0.AND.iecrou(i)==4)THEN
633 fxep(i)=fx(i)
634 dpx(i) = dx(i) - fx(i) / xk(i)
635 ENDIF
636 ENDDO
637 ENDIF
638C-------------------------------------
639C ELASTO-PLASTIC (TENSION COMPRESSION UNCOUPLED)
640C NONLINEAR LOADING/UNLOADING
641C DPX = MAXIMAL DISPLACEMENT (ELASTIC)
642C-------------------------------------
643 IF(jecrou(5)== nel )THEN
644 DO i=1,nel
645 IF(dx(i)>dpx(i))THEN
646 fx(i)=yy(i)
647 dpx(i) = dx(i)
648 ELSEIF(dx(i)>zero)THEN
649 dperm(i)=max(yy2(i),zero)
650 tmp1(i) = dperm(i)
651 IF(dx(i)>dperm(i).AND.yy3(i)/=zero)THEN
652 fmax(i)=yy3(i)/lscale(i)
653 tmp2(i) = fmax(i)
654 dperm(i)=min(dperm(i),dpx(i)- fmax(i) / xk(i))
655 b1 = (dpx(i)-dperm(i))*xk(i)/fmax(i)
656 fmin(i)=fmax(i)*
657 . ( (dx(i)-dperm(i))/(dpx(i)-dperm(i)) )**b1
658 fmax(i) = fmax(i)*(dx(i)-dperm(i))/(dpx(i)-dperm(i))
659 fx(i)=fxep(i)+xk(i)*ddx(i)
660 fx(i)=max(fx(i),fmin(i),zero)
661 fx(i)=min(fx(i),fmax(i),yy(i))
662 ELSE
663 fx(i) = zero
664 ENDIF
665 ELSEIF(dx(i)<dpx2(i))THEN
666 fx(i)=yy(i)
667 dpx2(i) = dx(i)
668 ELSEIF(dx(i)<zero)THEN
669 dperm(i)=min(yy2(i),zero)
670 tmp1(i) = dperm(i)
671 IF(dx(i)<dperm(i).AND.yy3(i)/=zero)THEN
672 fmax(i)=yy3(i)/lscale(i)
673 tmp2(i) = fmax(i)
674 dperm(i)=max(dperm(i),dpx2(i)- fmax(i) / xk(i))
675 b1 = (dpx2(i)-dperm(i))*xk(i)/fmax(i)
676 fmin(i)= fmax(i)*
677 . ( (-dx(i)+dperm(i)) / (-dpx2(i)+dperm(i)) )**b1
678 fmax(i) = fmax(i)*(dx(i)-dperm(i))/(dpx2(i)-dperm(i))
679 fx(i)=fxep(i)+xk(i)*ddx(i)
680 fx(i)=min(fx(i),fmin(i),zero)
681 fx(i)=max(fx(i),fmax(i),yy(i))
682 ELSE
683 fx(i) = zero
684 ENDIF
685 ENDIF
686 fxep(i)=fx(i)
687 ENDDO
688 ELSEIF(jecrou(5)>0)THEN
689 DO i=1,nel
690 IF(ifunc(i)/=0.AND.iecrou(i)==5)THEN
691 IF(dx(i)>dpx(i))THEN
692 fx(i)=yy(i)
693 dpx(i) = dx(i)
694 ELSEIF(dx(i)>zero)THEN
695 dperm(i)=max(yy2(i),zero)
696 IF(dx(i)>dperm(i).AND.yy3(i)/=zero)THEN
697 fmax(i)=yy3(i)/lscale(i)
698 dperm(i)=min(dperm(i),dpx(i)- fmax(i) / xk(i))
699C y = a (x-x1)^b
700 b1 = (dpx(i)-dperm(i))*xk(i)/fmax(i)
701 fmin(i) = fmax(i) *
702 . ( (dx(i)-dperm(i))/(dpx(i)-dperm(i)) )**b1
703 fmax(i) = fmax(i)*(dx(i)-dperm(i))/(dpx(i)-dperm(i))
704 fx(i)=fxep(i)+xk(i)*ddx(i)
705 fx(i)=max(fx(i),fmin(i),zero)
706 fx(i)=min(fx(i),fmax(i),yy(i))
707 ELSE
708 fx(i) = zero
709 ENDIF
710 ELSEIF(dx(i)<dpx2(i))THEN
711 fx(i)=yy(i)
712 dpx2(i) = dx(i)
713 ELSEIF(dx(i)<zero)THEN
714 dperm(i)=yy2(i)
715 dperm(i)=min(dperm(i),zero)
716 IF(dx(i)<dperm(i).AND.yy3(i)/=zero)THEN
717 fmax(i)=yy3(i)/lscale(i)
718 dperm(i)=max(dperm(i),dpx2(i)- fmax(i) / xk(i))
719C y = a (x-x1)^b
720 b1 = (dpx2(i)-dperm(i))*xk(i)/fmax(i)
721 fmin(i) = fmax(i)*
722 . ( (-dx(i)+dperm(i))/(-dpx2(i)+dperm(i)) )**b1
723 fmax(i) = fmax(i)*(dx(i)-dperm(i))/(dpx2(i)-dperm(i))
724 fx(i)=fxep(i)+xk(i)*ddx(i)
725 fx(i)=min(fx(i),fmin(i),zero)
726 fx(i)=max(fx(i),fmax(i),yy(i))
727 ELSE
728 fx(i) = zero
729 ENDIF
730 ENDIF
731 fxep(i)=fx(i)
732 ENDIF
733 ENDDO
734 ENDIF
735C-------------------------------------
736C ELASTO PLASTIC (ISOTROPIC HARDENING )
737C-------------------------------------
738 IF(jecrou(6) == nel )THEN
739 if(any_python_function) then
740 CALL vinter_mixed(python, tf,jad,jpos,jlen,nel,xx,dydx,yy)
741 ELSE
742 CALL vinter2(tf,jad ,jpos ,jlen ,nel ,xx ,dydx ,yy )
743 ENDIF
744 DO i=1,nel
745 IF(fx(i)>= zero.AND.fx(i)>yield(i))THEN
746 pos(1,i) = jpos(i)
747C-- COMPUTE PLASTIC AND ELASTIC DEFORMATION (TOTAL)
748 dpx(i)=dpx(i)+(fx(i)-yy(i))/an3y0(i)
749 dxela(i)=dx(i)-dpx(i)
750 fx(i)=yy(i) !H1
751 yield(i)=fx(i)
752C-- ECR variable for hardening/softening - always incremented with positive value
753 xx_old(i) = xx_old(i) + abs(ddx(i))
754C---FX< O
755 ELSEIF(fx(i)< zero.AND.fx(i)< -yield(i))THEN
756 pos(1,i) = jpos(i)
757C-- COMPUTE PLASTIC AND ELASTIC DEFORMATION (TOTAL)
758 dpx(i)=dpx(i)+(yy(i)-fx(i))/an3y0(i)
759 dxela(i)=dx(i)-dpx(i)
760 fx(i)=yy(i)
761 yield(i)=-fx(i)
762C-- ECR variable for hardening/softening - always incremented with positive value
763 xx_old(i) = xx_old(i) + abs(ddx(i))
764 ENDIF
765 fxep(i)=fx(i)
766 ENDDO
767 ELSEIF(jecrou(6)>0)THEN
768 if(any_python_function) then
769 CALL vinter_mixed(python, tf,jad,jpos,jlen,nel,xx,dydx,yy)
770 ELSE
771 CALL vinter2(tf,jad ,jpos ,jlen ,nel ,xx ,dydx ,yy )
772 ENDIF
773 DO i=1,nel
774 IF(ifunc(i)/= 0.AND.iecrou(i)== 6)THEN
775 IF(fx(i)>= zero.AND.fx(i)>yield(i))THEN
776 pos(1,i) = jpos(i)
777C-- COMPUTE PLASTIC AND ELASTIC DEFORMATION (TOTAL)
778 dpx(i)=dpx(i)+(fx(i)-yy(i))/an3y0(i)
779 dxela(i)=dx(i)-dpx(i)
780 fx(i)=yy(i) !H1
781 yield(i)=fx(i)
782C-- ECR variable for hardening/softening - always incremented with positive value
783 xx_old(i) = xx_old(i) + abs(ddx(i))
784C---FX< O
785 ELSEIF(fx(i)< zero.AND.fx(i)< -yield(i))THEN
786 pos(1,i) = jpos(i)
787C-- COMPUTE PLASTIC AND ELASTIC DEFORMATION (TOTAL)
788 dpx(i)=dpx(i)+(yy(i)-fx(i))/an3y0(i)
789 dxela(i)=dx(i)-dpx(i)
790 fx(i)=yy(i)
791 yield(i)=-fx(i)
792C-- ECR variable for hardening/softening - always incremented with positive value
793 xx_old(i) = xx_old(i) + abs(ddx(i))
794 ENDIF
795 fxep(i)=fx(i)
796 ENDIF
797 ENDDO
798 ENDIF
799C
800C-------------------------------------
801C ELSTO-PLASTIC
802C-------------------------------------
803 IF (jecrou(7)== nel ) THEN
804 DO i=1,nel
805 fx(i) = fxep(i) + xk(i)*ddx(i)
806 IF (dx(i)>= dxold(i).AND.dx(i)>=0)THEN
807 IF (fx(i)>yy(i)) fx(i) = yy(i)
808 ELSEIF(dx(i)< dxold(i).AND.dx(i)>=0)THEN
809 IF (fx(i)<yy2(i))fx(i) = yy2(i)
810 ELSEIF(dx(i)>= dxold(i).AND.dx(i)<0)THEN
811 IF(fx(i)> yy2(i))fx(i) = yy2(i)
812 ELSEIF(dx(i)< dxold(i).AND.dx(i)<0)THEN
813 IF(fx(i)< yy(i))fx(i) = yy(i)
814 ENDIF
815 ENDDO
816 DO i=1,nel
817 fxep(i)= fx(i)
818 dpx(i) = dx(i) - fx(i) / xk(i)
819 ENDDO
820 ELSEIF (jecrou(7) > 0) THEN
821 DO i=1,nel
822 IF (ifunc(i)/=0 .AND. iecrou(i)==7) THEN
823 fx(i) = fxep(i) + xk(i)*ddx(i)
824 IF (dx(i)>= dxold(i) .AND. dx(i)>=0) THEN
825 IF (fx(i)>yy(i)) fx(i) = yy(i)
826 ELSEIF (dx(i)< dxold(i) .AND. dx(i)>= 0) THEN
827 IF (fx(i) < yy2(i)) fx(i) = yy2(i)
828 ELSEIF (dx(i)>= dxold(i) .AND. dx(i)<0) THEN
829 IF (fx(i)> yy2(i)) fx(i) = yy2(i)
830 ELSEIF (dx(i)< dxold(i) .AND. dx(i)<0) THEN
831 IF (fx(i)< yy(i)) fx(i) = yy(i)
832 ENDIF
833 fxep(i) = fx(i)
834 dpx(i) = dx(i) - fx(i) / xk(i)
835 ENDIF
836 ENDDO
837 ENDIF
838C--------------------------------------------------------------------
839C NON LINEAR DAMPING
840C--------------------------------------------------------------------
841 IF(impl_s==0.OR.idyna>0) THEN
842 IF(jdmp>0)THEN
843 DO i=1,nel
844 jfunc=max(ifv(i),1)
845 pyid = python_funct_id(nfunct,jfunc,npf)
846 IF(pyid > 0) THEN
847 jpos(i) = -pyid
848 jad(i) = -pyid
849 any_python_function = .true.
850 ELSE
851 jpos(i) = nint(pos(4,i))
852 jfunc=max(ifv(i),1)
853 jad(i) = npf(jfunc) / 2 + 1
854 jlen(i) = npf(jfunc+1) / 2 - jad(i) - jpos(i)
855 ENDIF
856 ENDDO
857C
858 IF(any_python_function) THEN
859 CALL vinter_mixed(python, tf,jad,jpos,jlen,nel,dvxs,dydxv,gx)
860 ELSE
861 CALL vinter2(tf,jad,jpos,jlen,nel ,dvxs,dydxv,gx)
862 ENDIF
863C
864 DO i=1,nel
865 pos(4,i) = jpos(i)
866 ENDDO
867 ENDIF
868c---------------------------G * funct_id_4
869 IF(j2dmp>0)THEN
870 DO i=1,nel
871 j2func=max(ifunc3(i),1)
872 pyid = python_funct_id(nfunct,j2func,npf)
873 IF(pyid > 0) then ! python function
874 j2pos(i) = -pyid
875 j3ad(i) = -pyid
876 any_python_function = .true.
877 ELSE ! tabulated function
878 j2pos(i) = nint(pos(5,i))
879 j2func=max(ifunc3(i),1)
880 j2ad(i) = npf(j2func) / 2 + 1
881 j2len(i) = npf(j2func+1) / 2 - j2ad(i) - j2pos(i)
882 ENDIF
883 ENDDO
884 IF(any_python_function) THEN
885 CALL vinter_mixed(python, tf,j2ad,j2pos,j2len,nel,dvxs,dydxv2,gx2)
886 ELSE
887 CALL vinter2(tf,j2ad,j2pos,j2len,nel,dvxs,dydxv2,gx2)
888 ENDIF
889 DO i=1,nel
890 pos(5,i) = j2pos(i)
891 ENDDO
892 ENDIF
893c---------------------------K * funct_id_5
894 IF(j2k > 0)THEN
895 DO i=1,nel
896 j3func = max(ifunc4(i),1)
897 pyid = python_funct_id(nfunct,j3func,npf)
898 IF(pyid > 0) then ! python function
899 j3pos(i) = -pyid
900 j3ad(i) = -pyid
901 any_python_function = .true.
902 ELSE ! TABULATED FUNCTION
903 j3pos(i) = nint(pos(6,i))
904 j3func = max(ifunc4(i),1)
905 j3ad(i) = npf(j3func) / 2 + 1
906 j3len(i) = npf(j3func+1) / 2 - j3ad(i) - j3pos(i)
907 ENDIF
908 ENDDO
909 IF(any_python_function) THEN
910 CALL vinter_mixed(python, tf,j3ad,j3pos,j3len,nel,dvxs,dkdx2,kx2)
911 ELSE
912 CALL vinter2(tf,j3ad,j3pos,j3len,nel,xx,dkdx2,kx2)
913 ENDIF
914 DO i=1,nel
915 pos(6,i) = j3pos(i)
916 ENDDO
917 ENDIF
918c-------------------------
919 IF(j2k /= nel)THEN
920 DO i=1,nel
921 IF(ifunc4(i) == 0) kx2(i) = one
922 ENDDO
923 ENDIF
924 IF(jdmp/= nel)THEN
925 DO i=1,nel
926 IF(ifv(i)==0) gx(i)=zero
927 ENDDO
928 ENDIF
929 IF(j2dmp/= nel)THEN
930 DO i=1,nel
931 IF(ifunc3(i)==0) gx2(i)=zero
932 gx2(i)=gx2(i)*kx2(i)
933 ENDDO
934 ELSE
935 DO i=1,nel
936 gx2(i)=gx2(i)*kx2(i)
937 ENDDO
938 ENDIF
939
940 DO i=1,nel
941 dvv = max(one,abs(dvx(i)/d(i)))
942 dfac = ak(i) + b(i) * log(dvv) + ee(i)*gx(i)
943 fx(i)= ( dfac*fx(i) + xc(i)*dvx(i) + gf3(i)*gx2(i) ) *off(i)
944 e(i) = e(i) + (dx(i)-dxold(i)) * (fx(i)+fold(i)) * half
945 ENDDO
946 ELSE
947 DO i=1,nel
948 fx(i)= fx(i) *ak(i)* off(i)
949 e(i) = e(i) + (dx(i)-dxold(i)) * (fx(i)+fold(i)) * half
950 ENDDO
951 ENDIF
952 DO i=1,nel
953 dx(i)=dx(i)*xl0(i)
954 dxold(i)=dxold(i)*xl0(i)
955 dpx(i)=dpx(i)*xl0(i)
956 dpx2(i)=dpx2(i)*xl0(i)
957 e(i)=e(i)*xl0(i)
958 ENDDO
959C
960C----
961 RETURN
962 END
interp
subroutine interp(tf, tt, npoint, f, tg)
Definition interp.F:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
redef3_law113
subroutine redef3_law113(python, nel, nft, fx, xk, dx, fxep, dxold, dpx, tf, npf, xc, off, e, dpx2, anim, iani, pos, xl0, dmn, dmx, dvx, ff, lscale, ee, gf3, ifunc3, yield, ak, b, d, iecrou, ifunc, ifv, ifunc2, ifunc4, epla, xx_old)
Definition redef3_law113.F:43
vinter2
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter.F:143
vinter2dp
subroutine vinter2dp(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter.F:214