OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps60c.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/.
23!||====================================================================
24!|| sigeps60c ../engine/source/materials/mat/mat060/sigeps60c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!|| inter_rat_v ../engine/source/materials/mat/mat060/sigeps60c.F
30!|| vinter ../engine/source/tools/curve/vinter.F
31!||====================================================================
32 SUBROUTINE sigeps60c(
33 1 NEL0 ,NUPARAM,NUVAR ,MFUNC ,KFUNC ,
34 2 NPF ,NPT0 ,IPT ,IFLAG ,
35 2 TF ,TIME ,TIMESTEP,UPARAM ,RHO0 ,
36 3 AREA ,EINT ,THKLY ,ISRATE ,ASRATE ,
37 4 EPSPXX ,EPSPYY ,EPSPXY ,EPSPYZ ,EPSPZX ,
38 5 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
39 6 EPSXX ,EPSYY ,EPSXY ,EPSYZ ,EPSZX ,
40 7 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
41 8 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
42 9 SIGVXX ,SIGVYY ,SIGVXY ,SIGVYZ ,SIGVZX ,
43 A SOUNDSP,VISCMAX,THK ,PLA ,UVAR ,
44 B OFF ,NGL ,IPM ,MAT ,ETSE ,
45 C GS ,YLD ,EPSD_PG ,EPSP ,DPLA_I ,
46 D SHF ,INLOC ,DPLANL ,LOFF )
47C-----------------------------------------------
48C I m p l i c i t T y p e s
49C-----------------------------------------------
50#include "implicit_f.inc"
51C-----------------------------------------------
52C G l o b a l P a r a m e t e r s
53C-----------------------------------------------
54#include "mvsiz_p.inc"
55C---------+---------+---+---+--------------------------------------------
56C VAR | SIZE |TYP| RW| DEFINITION
57C---------+---------+---+---+--------------------------------------------
58C NEL0 | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL0
59C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
60C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
61C---------+---------+---+---+--------------------------------------------
62C NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
63C IFUNC | NFUNC | I | R | FUNCTION INDEX
64C NPF | * | I | R | FUNCTION ARRAY
65C NPT0 | 1 | I | R | NUMBER OF LAYERS OR INTEGRATION POINTS
66C IPT | 1 | I | R | LAYER OR INTEGRATION POINT NUMBER
67C IFLAG | * | I | R | GEOMETRICAL FLAGS
68C TF | * | F | R | FUNCTION ARRAY
69C---------+---------+---+---+--------------------------------------------
70C TIME | 1 | F | R | CURRENT TIME
71C TIMESTEP| 1 | F | R | CURRENT TIME STEP
72C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
73C RHO0 | NEL0 | F | R | INITIAL DENSITY
74C AREA | NEL0 | F | R | AREA
75C EINT | 2*NEL0 | F | R | INTERNAL ENERGY(MEMBRANE,BENDING)
76C THKLY | NEL0 | F | R | LAYER THICKNESS
77C EPSPXX | NEL0 | F | R | STRAIN RATE XX
78C EPSPYY | NEL0 | F | R | STRAIN RATE YY
79C ... | | | |
80C DEPSXX | NEL0 | F | R | STRAIN INCREMENT XX
81C DEPSYY | NEL0 | F | R | STRAIN INCREMENT YY
82C ... | | | |
83C EPSXX | NEL0 | F | R | STRAIN XX
84C EPSYY | NEL0 | F | R | STRAIN YY
85C ... | | | |
86C SIGOXX | NEL0 | F | R | OLD ELASTO PLASTIC STRESS XX
87C SIGOYY | NEL0 | F | R | OLD ELASTO PLASTIC STRESS YY
88C ... | | | |
89C---------+---------+---+---+--------------------------------------------
90C SIGNXX | NEL0 | F | W | NEW ELASTO PLASTIC STRESS XX
91C SIGNYY | NEL0 | F | W | NEW ELASTO PLASTIC STRESS YY
92C ... | | | |
93C SIGVXX | NEL0 | F | W | VISCOUS STRESS XX
94C SIGVYY | NEL0 | F | W | VISCOUS STRESS YY
95C ... | | | |
96C SOUNDSP | NEL0 | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
97C VISCMAX | NEL0 | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
98C---------+---------+---+---+--------------------------------------------
99C THK | NEL0 | F |R/W| THICKNESS
100C PLA | NEL0 | F |R/W| PLASTIC STRAIN
101C UVAR |NEL0*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
102C OFF | NEL0 | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
103C---------+---------+---+---+--------------------------------------------
104C C o m m o n B l o c k s
105C-----------------------------------------------
106#include "param_c.inc"
107#include "com01_c.inc"
108C-----------------------------------------------
109C I N P U T A r g u m e n t s
110C-----------------------------------------------
111C
112 INTEGER NEL0, NUPARAM, NUVAR, NPT0, IPT,IFLAG(*),
113 . NGL(NEL0),MAT(NEL0),ISRATE,IPM(NPROPMI,*),INLOC
114 my_real ,INTENT(IN) :: ASRATE
115 my_real ,DIMENSION(NEL0) ,INTENT(IN) :: EPSD_PG !< global element strain rate in Gauss pt
116 my_real ,DIMENSION(NEL0) ,INTENT(INOUT) :: EPSP !< local strain rate used in equations
117 my_real TIME,TIMESTEP,UPARAM(*),
118 . AREA(NEL0),RHO0(NEL0),EINT(NEL0,2),
119 . THKLY(NEL0),PLA(NEL0),
120 . EPSPXX(NEL0),EPSPYY(NEL0),
121 . EPSPXY(NEL0),EPSPYZ(NEL0),EPSPZX(NEL0),
122 . DEPSXX(NEL0),DEPSYY(NEL0),
123 . DEPSXY(NEL0),DEPSYZ(NEL0),DEPSZX(NEL0),
124 . EPSXX(NEL0) ,EPSYY(NEL0) ,
125 . EPSXY(NEL0) ,EPSYZ(NEL0) ,EPSZX(NEL0) ,
126 . sigoxx(nel0),sigoyy(nel0),
127 . sigoxy(nel0),sigoyz(nel0),sigozx(nel0),
128 . gs(*) ,shf(nel0), dplanl(nel0)
129 my_real, DIMENSION(NEL0), INTENT(IN) :: loff
130C-----------------------------------------------
131C O U T P U T A r g u m e n t s
132C-----------------------------------------------
133 my_real
134 . signxx(nel0),signyy(nel0),
135 . signxy(nel0),signyz(nel0),signzx(nel0),
136 . sigvxx(nel0),sigvyy(nel0),
137 . sigvxy(nel0),sigvyz(nel0),sigvzx(nel0),
138 . soundsp(nel0),viscmax(nel0),etse(nel0),
139 . dpla_i(nel0)
140C-----------------------------------------------
141C I N P U T O U T P U T A r g u m e n t s
142C-----------------------------------------------
143 my_real uvar(nel0,nuvar), off(nel0),thk(nel0),yld(nel0)
144C-----------------------------------------------
145C VARIABLES FOR FUNCTION INTERPOLATION
146C-----------------------------------------------
147 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
148 my_real FINTER ,TF(*)
149 EXTERNAL FINTER
150C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
151C Y : y = f(x)
152C X : x
153C DYDX : f'(x) = dy/dx
154C IFUNC(J): FUNCTION INDEX
155C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
156C NPF,TF : FUNCTION PARAMETER
157C-----------------------------------------------
158C L o c a l V a r i a b l e s
159C-----------------------------------------------
160 INTEGER I,J,NRATE(MVSIZ),J1(MVSIZ),J2(MVSIZ),N,NINDX,NMAX,IADBUF,NFUNC,
161 . iad1(mvsiz),ipos1(mvsiz),ilen1(mvsiz),
162 . iad2(mvsiz),ipos2(mvsiz),ilen2(mvsiz),
163 . iad3(mvsiz),ipos3(mvsiz),ilen3(mvsiz),
164 . iad4(mvsiz),ipos4(mvsiz),ilen4(mvsiz),
165 . jj(mvsiz),index(mvsiz),ifunc(mvsiz,100),nratem,
166 . nrate1,ifunc1(100),iadbufv(mvsiz),nfuncv(mvsiz),
167 . nfuncm, j3(mvsiz), j4(mvsiz), nn, index1,mopte,meinf,mce,ifunce
168 my_real
169 . e(mvsiz),a1(mvsiz),a2(mvsiz),g(mvsiz),g3(mvsiz),
170 . dydx1(mvsiz),dydx2(mvsiz),rate(mvsiz,4),svm(mvsiz),
171 . y1(mvsiz),y2(mvsiz),dr(mvsiz),
172 . yfac(mvsiz,4),nnu1(mvsiz),nu1(mvsiz),
173 . nu2(mvsiz),nu3(mvsiz),nu4(mvsiz),nu5(mvsiz),nu6(mvsiz),
174 . aa(mvsiz),bb(mvsiz),dpla_j(mvsiz),
175 . pp(mvsiz),qq(mvsiz),fail(mvsiz),svmo(mvsiz),h(mvsiz),
176 . epsmax(mvsiz),epsr1(mvsiz),epsr2(mvsiz),fisokin(mvsiz),
177 . sigexx(mvsiz),sigeyy(mvsiz),sigexy(mvsiz),
178 . hk(mvsiz),nu(mvsiz),y3(mvsiz),y4(mvsiz),
179 . dydx3(mvsiz),dydx4(mvsiz),einf(mvsiz),
180 . ce(mvsiz),escale(mvsiz),e0(mvsiz),dydxe(mvsiz)
181 my_real
182 . r,umr,a,b,c,amu,s11,s22,s12,p,p2,fac,dezz,
183 . sigz,s1,s2,s3,dpla,vm2,epst,nnu2,ce1, einf1,
184 . err,f,df,pla_i,q2,yld_i,sigpxx,sigpyy,sigpxy,
185 . alpha,
186 . e1,a11,a21,g1,g31,nnu11,nu11,nu21,nu31,nu41,nu51,nu61,
187 . epsmax1,epsr11,epsr21,fisokin1,
188 . dsxx,dsyy,dsxy,dexx,deyy,dexy,nux, x(mvsiz)
189 my_real
190 . me, ma1, ma2, mg, mnu,
191 . mepsmax,mepsr1,mepsr2,mfisokin,
192 . y11(mvsiz), y22(mvsiz), y33(mvsiz), y44(mvsiz), y(mvsiz),yp(mvsiz),
193 . x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz)
194 INTEGER JST(MVSIZ+1), IC, MNRATE, OPTE1,OPTE(MVSIZ)
195C
196 DATA nmax/3/
197C-----------------------------------------------
198C USER VARIABLES INITIALIZATION
199C-----------------------------------------------
200 nfunc = ipm(10,mat(1))
201 DO j=1,nfunc
202 ifunc1(j)=ipm(10+j,mat(1))
203 ENDDO
204C
205 iadbuf = ipm(7,mat(1)) - 1
206 e1 = uparam(iadbuf+2)
207 a11 = uparam(iadbuf+3)
208 a21 = uparam(iadbuf+4)
209 g1 = uparam(iadbuf+5)
210 g31 = 3.*g1
211 nux = uparam(iadbuf+6)
212 nrate1 = nint(uparam(iadbuf+1))
213 epsmax1=uparam(iadbuf+6+2*nrate1+1)
214c-----------------------------
215 opte1 = uparam(iadbuf+2*nrate1 + 18)
216 einf1 = uparam(iadbuf+2*nrate1 + 19)
217 ce1 = uparam(iadbuf+2*nrate1 + 20)
218c-----------------------------
219 DO i=1,nel0
220 e(i) = e1
221 a1(i) = a11
222 a2(i) = a21
223 g(i) = g1
224 g3(i) = g31
225 ENDDO
226 IF (opte1 == 1)THEN
227 ifunce = uparam(iadbuf+2*nrate1 + 17)
228 DO i=1,nel0
229 IF(pla(i) > zero)THEN
230 escale(i) = finter(kfunc(ifunce),pla(i),npf,tf,dydxe(i))
231 ENDIF
232 ENDDO
233 DO i=1,nel0
234 IF(pla(i) > zero)THEN
235 e(i) = escale(i)* e1
236 a1(i) = e(i)/(one - nux*nux)
237 a2(i) = nux*a1(i)
238 g(i) = half*e(i)/(one+nux)
239 g3(i) = three*g(i)
240 gs(i) = g(i) * shf(i)
241 ENDIF
242 ENDDO
243 ELSEIF ( ce1 /= zero) THEN
244 DO i=1,nel0
245 IF(pla(i) > zero)THEN
246 e(i) = e1-(e1-einf1)*(one-exp(-ce1*pla(i)))
247 a1(i) = e(i)/(one - nux*nux)
248 a2(i) = nux*a1(i)
249 g(i) = half*e(i)/(one+nux)
250 g3(i) = three*g(i)
251 gs(i) = g(i) * shf(i)
252 ENDIF
253 ENDDO
254 ENDIF
255 IF(epsmax1==0.)THEN
256 IF(tf(npf(ifunc1(1)+1)-1)==zero)THEN
257 epsmax1=tf(npf(ifunc1(1)+1)-2)
258 ELSE
259 epsmax1= ep30
260 ENDIF
261 ENDIF
262C
263 nnu11 = nux / (one - nux)
264 nnu2 = nnu11*nnu11
265 nu11 = one/(one - nux)
266 nu21 = one/(one + nux)
267 nu31 = one - nnu11
268 nu41 = one + nnu2 + nnu11
269 nu51 = one + nnu2 - two*nnu11
270 nu61 = half - nnu2 + half*nnu11
271C
272 epsr11 =uparam(iadbuf+6+2*nrate1+2)
273 epsr21 =uparam(iadbuf+6+2*nrate1+3)
274 fisokin1=uparam(iadbuf+6+2*nrate1+8)
275C
276 IF (isigi==0) THEN
277 IF(time==zero)THEN
278 DO i=1,nel0
279 uvar(i,1)=zero
280 uvar(i,2)=zero
281 uvar(i,3)=zero
282 uvar(i,4)=zero
283 ENDDO
284 DO j=1,nrate1
285 DO i=1,nel0
286 uvar(i,j+4)=0
287 ENDDO
288 ENDDO
289 ENDIF
290 ENDIF
291C------------------------------------------
292C ECROUISSAGE CINE
293C------------------------------------------
294 DO i=1,nel0
295! sigoxx(i) = sigoxx(i) - uvar(i,2)
296! SIGOYY(I) = SIGOYY(I) - UVAR(I,3)
297! SIGOXY(I) = SIGOXY(I) - UVAR(I,4)
298 ENDDO
299C-----------------------------------------------
300C
301 DO i=1,nel0
302C
303 signxx(i)=sigoxx(i) - uvar(i,2) +a1(i)*depsxx(i)+a2(i)*depsyy(i)
304 signyy(i)=sigoyy(i) - uvar(i,3) +a2(i)*depsxx(i)+a1(i)*depsyy(i)
305 signxy(i)=sigoxy(i) - uvar(i,4) +g(i) *depsxy(i)
306 signyz(i)=sigoyz(i)+gs(i) *depsyz(i)
307 signzx(i)=sigozx(i)+gs(i) *depszx(i)
308 sigexx(i) = signxx(i)
309 sigeyy(i) = signyy(i)
310 sigexy(i) = signxy(i)
311C
312 soundsp(i) = sqrt(a1(i)/rho0(i))
313 viscmax(i) = zero
314 etse(i) = one
315C-------------------
316C STRAIN RATE
317C-------------------
318C
319 IF (israte==0) THEN
320 epsp(i) = half*( abs(epspxx(i)+epspyy(i))
321 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
322 . + epspxy(i)*epspxy(i) ) )
323 ELSE
324 epsp(i) = asrate*epsd_pg(i) + (one-asrate)*epsp(i)
325 ENDIF
326C-------------------
327C STRAIN
328C-------------------
329C
330 epst = half*( epsxx(i)+epsyy(i)
331 . + sqrt( (epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
332 . + epsxy(i)*epsxy(i) ) )
333 fail(i) = max(zero,min(one,(epsr21-epst)/(epsr21-epsr11)))
334C
335 ENDDO
336C-------------------
337C CRITERE
338C-------------------
339 DO i=1,nel0
340 jj(i) = 1
341 ENDDO
342 iadbuf = ipm(7,mat(1)) - 1
343C inversion boucles
344 DO j=2,nrate1-1
345 DO i=1,nel0
346 IF(epsp(i)>=uparam(iadbuf+6+j)) jj(i) = j
347 ENDDO
348 ENDDO
349#include "vectorize.inc"
350 DO i=1,nel0
351 IF(jj(i)==1)THEN
352 j1(i) = jj(i)
353 j2(i) = j1(i)+1
354 j3(i) = j2(i)+1
355 j4(i) = j3(i)+1
356 ELSEIF(jj(i)==(nrate1-1))THEN
357 j1(i) = jj(i) - 2
358 j2(i) = j1(i)+1
359 j3(i) = j2(i)+1
360 j4(i) = j3(i)+1
361 ELSE
362 j1(i) = jj(i) - 1
363 j2(i) = j1(i)+1
364 j3(i) = j2(i)+1
365 j4(i) = j3(i)+1
366 ENDIF
367 ENDDO
368#include "vectorize.inc"
369 DO i=1,nel0
370 rate(i,1)=uparam(iadbuf + 6 + j1(i) )
371 rate(i,2)=uparam(iadbuf + 6 + j2(i) )
372 rate(i,3)=uparam(iadbuf + 6 + j3(i) )
373 rate(i,4)=uparam(iadbuf + 6 + j4(i) )
374 yfac(i,1)=uparam(iadbuf + 6 + nrate1 + j1(i) )
375 yfac(i,2)=uparam(iadbuf + 6 + nrate1 + j2(i) )
376 yfac(i,3)=uparam(iadbuf + 6 + nrate1 + j3(i) )
377 yfac(i,4)=uparam(iadbuf + 6 + nrate1 + j4(i) )
378 ENDDO
379C
380#include "vectorize.inc"
381 DO i=1,nel0
382 IF(jj(i)==1)THEN
383 j1(i) = jj(i)
384 j2(i) = j1(i)+1
385 j3(i) = j2(i)+1
386 j4(i) = j3(i)+1
387 ELSEIF(jj(i)==(nrate1-1))THEN
388 j1(i) = jj(i) - 2
389 j2(i) = j1(i)+1
390 j3(i) = j2(i)+1
391 j4(i) = j3(i)+1
392 ELSE
393 j1(i) = jj(i) - 1
394 j2(i) = j1(i)+1
395 j3(i) = j2(i)+1
396 j4(i) = j3(i)+1
397 ENDIF
398 ENDDO
399#include "vectorize.inc"
400 DO i=1,nel0
401 ipos1(i) = nint(uvar(i,j1(i)+4))
402 iad1(i) = npf(ifunc1(j1(i))) / 2 + 1
403 ilen1(i) = npf(ifunc1(j1(i))+1) / 2 - iad1(i) - ipos1(i)
404 ipos2(i) = nint(uvar(i,j2(i)+4))
405 iad2(i) = npf(ifunc1(j2(i))) / 2 + 1
406 ilen2(i) = npf(ifunc1(j2(i))+1) / 2 - iad2(i) - ipos2(i)
407 ipos3(i) = nint(uvar(i,j3(i)+4))
408 iad3(i) = npf(ifunc1(j3(i))) / 2 + 1
409 ilen3(i) = npf(ifunc1(j3(i))+1) / 2 - iad3(i) - ipos3(i)
410 ipos4(i) = nint(uvar(i,j4(i)+4))
411 iad4(i) = npf(ifunc1(j4(i))) / 2 + 1
412 ilen4(i) = npf(ifunc1(j4(i))+1) / 2 - iad4(i) - ipos4(i)
413 ENDDO
414C
415 CALL vinter(tf,iad1,ipos1,ilen1,nel0,pla,dydx1,y1)
416 CALL vinter(tf,iad2,ipos2,ilen2,nel0,pla,dydx2,y2)
417 CALL vinter(tf,iad3,ipos3,ilen3,nel0,pla,dydx3,y3)
418 CALL vinter(tf,iad4,ipos4,ilen4,nel0,pla,dydx4,y4)
419C
420C
421 IF (fisokin1==0.) THEN
422! ------------------------
423#include "vectorize.inc"
424 DO i=1,nel0
425 IF(jj(i)==1)THEN
426 j1(i) = jj(i)
427 j2(i) = j1(i)+1
428 j3(i) = j2(i)+1
429 j4(i) = j3(i)+1
430 dydx1(i)= dydx1(i)*yfac(i,1)
431 dydx2(i) = dydx2(i)*yfac(i,2)
432 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
433 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
434 ELSEIF(jj(i)==(nrate1-1))THEN
435 j1(i) = jj(i) - 2
436 j2(i) = j1(i)+1
437 j3(i) = j2(i)+1
438 j4(i) = j3(i)+1
439 dydx3(i) = dydx3(i)*yfac(i,3)
440 dydx4(i) = dydx4(i)*yfac(i,4)
441 fac = (epsp(i) - rate(i,3))/(rate(i,4) - rate(i,3))
442 h(i) = fail(i)*(dydx3(i) + fac*(dydx4(i)-dydx3(i)))
443 ELSE
444 j1(i) = jj(i) - 1
445 j2(i) = j1(i)+1
446 j3(i) = j2(i)+1
447 j4(i) = j3(i)+1
448 dydx2(i) = dydx2(i)*yfac(i,2)
449 dydx3(i) = dydx3(i)*yfac(i,3)
450 fac = (epsp(i) - rate(i,2))/(rate(i,3) - rate(i,2))
451 h(i) = fail(i)*(dydx2(i) + fac*(dydx3(i)-dydx2(i)))
452 ENDIF
453 ENDDO
454 DO i=1,nel0
455 uvar(i,j1(i)+4) = ipos1(i)
456 uvar(i,j2(i)+4) = ipos2(i)
457 uvar(i,j3(i)+4) = ipos3(i)
458 uvar(i,j4(i)+4) = ipos4(i)
459 ENDDO
460 x1(1:nel0) = rate(1:nel0,1)
461 x2(1:nel0) = rate(1:nel0,2)
462 x3(1:nel0) = rate(1:nel0,3)
463 x4(1:nel0) = rate(1:nel0,4)
464 x(1:nel0) = epsp(1:nel0)
465 y11(1:nel0) = y1(1:nel0)*yfac(1:nel0,1)
466 y22(1:nel0) = y2(1:nel0)*yfac(1:nel0,2)
467 y33(1:nel0) = y3(1:nel0)*yfac(1:nel0,3)
468 y44(1:nel0) = y4(1:nel0)*yfac(1:nel0,4)
469 CALL inter_rat_v(nel0,x1,x2,x3,x4,y11,y22,y33,y44,
470 . x,y,yp,jj, nrate1)
471 yld(1:nel0) = fail(1:nel0)*y(1:nel0)
472 yld(1:nel0) = max(yld(1:nel0),em20)
473! ------------------------
474 ELSEIF (fisokin1==one) THEN
475#include "vectorize.inc"
476 DO i=1,nel0
477 IF(jj(i)==1)THEN
478 j1(i) = jj(i)
479 j2(i) = j1(i)+1
480 j3(i) = j2(i)+1
481 j4(i) = j3(i)+1
482 dydx1(i)= dydx1(i)*yfac(i,1)
483 dydx2(i) = dydx2(i)*yfac(i,2)
484 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
485 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
486 ELSEIF(jj(i)==(nrate1-1))THEN
487 j1(i) = jj(i) - 2
488 j2(i) = j1(i)+1
489 j3(i) = j2(i)+1
490 j4(i) = j3(i)+1
491 dydx3(i) = dydx3(i)*yfac(i,3)
492 dydx4(i) = dydx4(i)*yfac(i,4)
493 fac = (epsp(i) - rate(i,3))/(rate(i,4) - rate(i,3))
494 h(i) = fail(i)*(dydx3(i) + fac*(dydx4(i)-dydx3(i)))
495 ELSE
496 j1(i) = jj(i) - 1
497 j2(i) = j1(i)+1
498 j3(i) = j2(i)+1
499 j4(i) = j3(i)+1
500 dydx2(i) = dydx2(i)*yfac(i,2)
501 dydx3(i) = dydx3(i)*yfac(i,3)
502 fac = (epsp(i) - rate(i,2))/(rate(i,3) - rate(i,2))
503 h(i) = fail(i)*(dydx2(i) + fac*(dydx3(i)-dydx2(i)))
504 ENDIF
505 ENDDO
506 DO i=1,nel0
507 uvar(i,j1(i)+4) = ipos1(i)
508 uvar(i,j2(i)+4) = ipos2(i)
509 uvar(i,j3(i)+4) = ipos3(i)
510 uvar(i,j4(i)+4) = ipos4(i)
511 ENDDO
512C ECROUISSAGE CINEMATIQUE
513 y11(1:nel0)=tf(npf(ifunc1(j1(1:nel0)))+1)*yfac(1:nel0,1)
514 y22(1:nel0)=tf(npf(ifunc1(j2(1:nel0)))+1)*yfac(1:nel0,2)
515 y33(1:nel0)=tf(npf(ifunc1(j3(1:nel0)))+1)*yfac(1:nel0,3)
516 y44(1:nel0)=tf(npf(ifunc1(j4(1:nel0)))+1)*yfac(1:nel0,4)
517 x1(1:nel0) = rate(1:nel0,1)
518 x2(1:nel0) = rate(1:nel0,2)
519 x3(1:nel0) = rate(1:nel0,3)
520 x4(1:nel0) = rate(1:nel0,4)
521 x(1:nel0) = epsp(1:nel0)
522 CALL inter_rat_v(nel0,x1,x2,x3,x4,y11,y22,y33,y44,
523 . x,y,yp,jj,nrate1)
524 yld(1:nel0) = fail(1:nel0)*y(1:nel0)
525 ELSE
526#include "vectorize.inc"
527 DO i=1,nel0
528 IF(jj(i)==1)THEN
529 j1(i) = jj(i)
530 j2(i) = j1(i)+1
531 j3(i) = j2(i)+1
532 j4(i) = j3(i)+1
533 dydx1(i)= dydx1(i)*yfac(i,1)
534 dydx2(i) = dydx2(i)*yfac(i,2)
535 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
536 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
537 ELSEIF(jj(i)==(nrate1-1))THEN
538 j1(i) = jj(i) - 2
539 j2(i) = j1(i)+1
540 j3(i) = j2(i)+1
541 j4(i) = j3(i)+1
542 dydx3(i) = dydx3(i)*yfac(i,3)
543 dydx4(i) = dydx4(i)*yfac(i,4)
544 fac = (epsp(i) - rate(i,3))/(rate(i,4) - rate(i,3))
545 h(i) = fail(i)*(dydx3(i) + fac*(dydx4(i)-dydx3(i)))
546 ELSE
547 j1(i) = jj(i) - 1
548 j2(i) = j1(i)+1
549 j3(i) = j2(i)+1
550 j4(i) = j3(i)+1
551 dydx2(i) = dydx2(i)*yfac(i,2)
552 dydx3(i) = dydx3(i)*yfac(i,3)
553 fac = (epsp(i) - rate(i,2))/(rate(i,3) - rate(i,2))
554 h(i) = fail(i)*(dydx2(i) + fac*(dydx3(i)-dydx2(i)))
555 ENDIF
556 ENDDO
557 x1(1:nel0) = rate(1:nel0,1)
558 x2(1:nel0) = rate(1:nel0,2)
559 x3(1:nel0) = rate(1:nel0,3)
560 x4(1:nel0) = rate(1:nel0,4)
561 x(1:nel0) = epsp(1:nel0)
562 y11(1:nel0) = y1(1:nel0)*yfac(1:nel0,1)
563 y22(1:nel0) = y2(1:nel0)*yfac(1:nel0,2)
564 y33(1:nel0) = y3(1:nel0)*yfac(1:nel0,3)
565 y44(1:nel0) = y4(1:nel0)*yfac(1:nel0,4)
566 CALL inter_rat_v(nel0,x1,x2,x3,x4,y11,y22,y33,y44,
567 . x,y,yp,jj,nrate1)
568 yld(1:nel0) = fail(1:nel0)*y(1:nel0)
569 yld(1:nel0) = max(yld(1:nel0),em20)
570
571 DO i=1,nel0
572 uvar(i,j1(i)+4) = ipos1(i)
573 uvar(i,j2(i)+4) = ipos2(i)
574 uvar(i,j3(i)+4) = ipos3(i)
575 uvar(i,j4(i)+4) = ipos4(i)
576 ENDDO
577C ECROUISSAGE CINEMATIQUE
578 y11(1:nel0)=tf(npf(ifunc1(j1(1:nel0)))+1)*yfac(1:nel0,1)
579 y22(1:nel0)=tf(npf(ifunc1(j2(1:nel0)))+1)*yfac(1:nel0,2)
580 y33(1:nel0)=tf(npf(ifunc1(j3(1:nel0)))+1)*yfac(1:nel0,3)
581 y44(1:nel0)=tf(npf(ifunc1(j4(1:nel0)))+1)*yfac(1:nel0,4)
582 CALL inter_rat_v(nel0,x1,x2,x3,x4,y11,y22,y33,y44,
583 . x,y,yp,jj,nrate1)
584 yld(1:nel0) = (one -fisokin1) * yld(1:nel0) + fisokin1 * (fail(1:nel0)*y(1:nel0))
585 ENDIF
586
587C-------------------
588C PROJECTION
589C-------------------
590 IF(iflag(1)==0)THEN
591C projection radiale
592 DO i=1,nel0
593 svm(i)=sqrt(signxx(i)*signxx(i)
594 . +signyy(i)*signyy(i)
595 . -signxx(i)*signyy(i)
596 . + three*signxy(i)*signxy(i))
597 r = min(one,yld(i)/max(em20,svm(i)))
598 signxx(i)=signxx(i)*r
599 signyy(i)=signyy(i)*r
600 signxy(i)=signxy(i)*r
601 umr = one - r
602 dpla_i(i) = off(i)*svm(i)*umr/(g3(i)+h(i))
603 pla(i) = pla(i) + dpla_i(i)
604 s1=half*(signxx(i)+signyy(i))
605 IF (inloc == 0) THEN
606 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) / yld(i)
607 dezz=-(depsxx(i)+depsyy(i))*nnu11-nu31*dezz
608 thk(i) = thk(i) + dezz*thkly(i)*off(i)
609 ENDIF
610 IF(r<one) etse(i)= h(i)/(h(i)+e(i))
611 ENDDO
612C
613 ELSEIF(iflag(1)==1)THEN
614C-------------------------
615C CRITERE DE VON MISES
616C-------------------------
617 DO i=1,nel0
618 h(i) = max(zero,h(i))
619 s1=signxx(i)+signyy(i)
620 s2=signxx(i)-signyy(i)
621 s3=signxy(i)
622 aa(i)=fourth*s1*s1
623 bb(i)=three_over_4*s2*s2+three*s3*s3
624 svm(i)=sqrt(aa(i)+bb(i))
625 IF (inloc == 0) THEN
626 dezz = -(depsxx(i)+depsyy(i))*nnu11
627 thk(i) = thk(i) + dezz*thkly(i)*off(i)
628 ENDIF
629 ENDDO
630C-------------------------
631C GATHER PLASTIC FLOW
632C-------------------------
633 nindx=0
634 DO i=1,nel0
635 IF(svm(i)>yld(i).AND.off(i)==one) THEN
636 nindx=nindx+1
637 index(nindx)=i
638 ENDIF
639 ENDDO
640C
641 IF(nindx/=0) THEN
642C---------------------------
643C DEP EN CONTRAINTE PLANE
644C---------------------------
645#include "vectorize.inc"
646 DO j=1,nindx
647 i=index(j)
648 dpla_j(i)=(svm(i)-yld(i))/(g3(i)+h(i))
649 etse(i)= h(i)/(h(i)+e(i))
650 hk(i) = h(i)*(one-fisokin1)
651 ENDDO
652C
653 DO n=1,nmax
654#include "vectorize.inc"
655 DO j=1,nindx
656 i=index(j)
657 dpla_i(i)=dpla_j(i)
658 yld_i =yld(i)+hk(i)*dpla_i(i)
659 dr(i) =half*e(i)*dpla_i(i)/yld_i
660 pp(i) =one/(one + dr(i)*nu11)
661 qq(i) =one/(one + three*dr(i)*nu21)
662 p2 =pp(i)*pp(i)
663 q2 =qq(i)*qq(i)
664 f =aa(i)*p2+bb(i)*q2-yld_i*yld_i
665 df =-(aa(i)*nu11*p2*pp(i)+ three*bb(i)*nu21*q2*qq(i))
666 . *(e(i)- two*dr(i)*hk(i))/yld_i
667 . - two*hk(i)*yld_i
668 df = sign(max(abs(df),em20),df)
669 IF(dpla_i(i)>zero) THEN
670 dpla_j(i)=max(zero,dpla_i(i)-f/df)
671 ELSE
672 dpla_j(i)=zero
673 ENDIF
674 ENDDO
675 ENDDO
676C------------------------------------------
677C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
678C------------------------------------------
679#include "vectorize.inc"
680 DO j=1,nindx
681 i=index(j)
682 pla(i) = pla(i) + dpla_i(i)
683 s1=(signxx(i)+signyy(i))*pp(i)
684 s2=(signxx(i)-signyy(i))*qq(i)
685 signxx(i)=half*(s1+s2)
686 signyy(i)=half*(s1-s2)
687 signxy(i)=signxy(i)*qq(i)
688 IF (inloc == 0) THEN
689 dezz = - nu31*dr(i)*s1/e(i)
690 thk(i) = thk(i) + dezz*thkly(i)*off(i)
691 ENDIF
692 ENDDO
693 ENDIF
694C-------------------------------------------
695 ELSEIF(iflag(1)==2)THEN
696C projection radial sur le deviateur sur un critere reduit
697C projection elastique en z => sig33 = 0
698C le coef. de reduction du critere est tel que
699C l'on se trouve sur le critere apres les 2 projections
700 DO i=1,nel0
701 p = -(signxx(i)+signyy(i))*third
702 s11 = signxx(i)+p
703 s22 = signyy(i)+p
704C s33 = p = -(S11 + S22)
705 s12 = signxy(i)
706C
707 p2 = p*p
708 vm2= three*(s12*s12 - s11*s22)
709 a = p2*nu41 + vm2
710 vm2= three*p2 + vm2
711 b = p2*nu61
712 c = p2*nu51 - yld(i)*yld(i)
713 r = min(one,(-b + sqrt(max(zero,b*b-a*c)))/max(a ,em20))
714 signxx(i) = s11*r - p
715 signyy(i) = s22*r - p
716 signxy(i) = s12*r
717 umr = one - r
718 sigz = nnu11*p*umr
719 signxx(i) = signxx(i) + sigz
720 signyy(i) = signyy(i) + sigz
721 svm(i)=sqrt(vm2)
722 dpla_i(i) = off(i)*svm(i)*umr/(g3(i)+h(i))
723 pla(i) = pla(i) + dpla_i(i)
724 IF (inloc == 0) THEN
725 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) / yld(i)
726 dezz=-(depsxx(i)+depsyy(i))*nnu11-nu31*dezz
727 thk(i) = thk(i) + dezz*thkly(i)*off(i)
728 ENDIF
729 IF(r<one) etse(i)= h(i)/(h(i)+e(i))
730 ENDDO
731 ENDIF
732C
733 DO i=1,nel0
734 IF(pla(i)>epsmax1.AND.off(i)==one)off(i)=four_over_5
735 ENDDO
736C------------------------------------------
737C ECROUISSAGE CINE
738C------------------------------------------
739C test bypass
740 IF (fisokin1/=zero) THEN
741 DO i=1,nel0
742 dsxx = sigexx(i) - signxx(i)
743 dsyy = sigeyy(i) - signyy(i)
744 dsxy = sigexy(i) - signxy(i)
745 dexx = (dsxx - nux*dsyy)
746 deyy = (dsyy - nux*dsxx)
747 dexy = two*(one + nux)*dsxy
748 alpha = fisokin1*h(i)/(e(i)+h(i))/three
749 sigpxx = alpha*(four*dexx + two*deyy)
750 sigpyy = alpha*(four*deyy + two*dexx)
751 sigpxy = alpha*dexy
752 signxx(i) = signxx(i) + uvar(i,2)
753 signyy(i) = signyy(i) + uvar(i,3)
754 signxy(i) = signxy(i) + uvar(i,4)
755 uvar(i,2) = uvar(i,2) + sigpxx
756 uvar(i,3) = uvar(i,3) + sigpyy
757 uvar(i,4) = uvar(i,4) + sigpxy
758 ENDDO
759 ENDIF
760C--------------------------------
761C NON-LOCAL THICKNESS VARIATION
762C--------------------------------
763 IF (inloc > 0) THEN
764 DO i = 1,nel0
765 IF (loff(i) == one) THEN
766 svm(i) = sqrt(signxx(i)*signxx(i)
767 . + signyy(i)*signyy(i)
768 . - signxx(i)*signyy(i)
769 . + three*signxy(i)*signxy(i))
770 dezz = max(dplanl(i),zero)*half*(signxx(i)+signyy(i))/max(svm(i),em20)
771 dezz = -nux*((signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/e1) - dezz
772 thk(i) = thk(i) + dezz*thkly(i)*off(i)
773 ENDIF
774 ENDDO
775 ENDIF
776C
777 RETURN
778 END
779!||====================================================================
780!|| inter_rat ../engine/source/materials/mat/mat060/sigeps60c.F
781!||--- called by ------------------------------------------------------
782!|| sigeps60 ../engine/source/materials/mat/mat060/sigeps60.f
783!|| sigeps60g ../engine/source/materials/mat/mat060/sigeps60g.F
784!||====================================================================
785 SUBROUTINE inter_rat(X0,X1,X2,X3,Y0,Y1,Y2,Y3,X,Y,YP,I,N)
786C
787C INTERPOLATION (RATIONAL FUNCTION)
788C
789C-----------------------------------------------
790C I m p l i c i t T y p e s
791C-----------------------------------------------
792#include "implicit_f.inc"
793C-----------------------------------------------
794C D u m m y A r g u m e n t s
795C-----------------------------------------------
796 INTEGER I, N
797C REAL
798 my_real
799 . x0, x1, x2, x3, y0, y1, y2, y3, x, y, yp
800C-----------------------------------------------
801C L o c a l V a r i a b l e s
802C-----------------------------------------------
803C REAL
804 my_real
805 . q, d, r, s, sp, c2, dm, sm, c1, c6, c4, c3, c3d, c5,
806 . a1, a0, a2, t, y11
807C
808 q = x-x1
809 d = x2-x1
810 r = d-q
811 s = (y2-y1) / d
812 sp = (y3-y2) / (x3-x2)
813 c2 = (sp-s) / (x3-x1)
814 dm = x1-x0
815 dm = sign(max(em30,abs(dm)),dm)
816 sm = (y1-y0) / dm
817 c1 = (s-sm) / (d+dm)
818 c6 = 0.
819C
820 IF(i==1)THEN
821 IF(x<=x0) THEN
822 c1 = zero
823 c2 = zero
824 sm = zero
825 ELSE
826 c2 = c1
827 c1 = sm/(x1 - x0)
828 ENDIF
829 r = (x1 - x)
830 q = x - x0
831 d = x1 - x0
832 c3 = abs(c2*r)
833 c3d = c3+abs(c1*q)
834 c5 = zero
835 IF(c3d>zero) THEN
836 c3 = c3/c3d
837 c5 = c3*(c1-c2)
838 ENDIF
839 c4 = c2+c5
840 c6 = d*c5*(one - c3)
841 y = y0 + q*(sm-r*c4)
842 yp = sm+(q-r)*c4+c6
843 ELSE IF(i==n-1)THEN
844 IF(sp==zero.OR.(x>x3)) THEN
845 c1 = zero
846 c2 = zero
847 ELSE
848 c1 = (sp - s)/(x3 - x1)
849 c2 = zero
850 ENDIF
851 r = (x3 - x)
852 q = x - x2
853 d = x3 - x2
854 c3 = abs(c2*r)
855 c3d = c3+abs(c1*q)
856 c5 = zero
857 IF(c3d>zero) THEN
858 c3 = c3/c3d
859 c5 = c3*(c1-c2)
860 ENDIF
861 c4 = c2+c5
862 c6 = d*c5*(one - c3)
863 y = y2 + (x - x2)*(sp-(x3 - x)*c4)
864 yp = sp+(q-r)*c4+c6
865 ELSE
866 IF(i==2.AND.sm*(sm-dm*c1)<=zero)c1=(s-sm-sm) / d
867 c3 = abs(c2*r)
868 c3d = c3+abs(c1*q)
869 c5 = zero
870 IF(c3d>zero) THEN
871 c3 = c3/c3d
872 c5 = c3*(c1-c2)
873 ENDIF
874 c4 = c2+c5
875 c6 = d*c5*(one - c3)
876 y = y1 + q*(s-r*c4)
877 yp = s+(q-r)*c4+c6
878C
879 ENDIF
880
881 END
882
883
884!||====================================================================
885!|| inter_rat_v ../engine/source/materials/mat/mat060/sigeps60c.F
886!||--- called by ------------------------------------------------------
887!|| sigeps60c ../engine/source/materials/mat/mat060/sigeps60c.F
888!||====================================================================
889 SUBROUTINE inter_rat_v(NEL0,X0,X1,X2,X3,Y0,Y1,Y2,Y3,X,Y,YP,JJ,N)
890!
891! INTERPOLATION (RATIONAL FUNCTION)
892!
893C-----------------------------------------------
894C I m p l i c i t T y p e s
895C-----------------------------------------------
896#include "implicit_f.inc"
897C-----------------------------------------------
898C D u m m y A r g u m e n t s
899C-----------------------------------------------
900 INTEGER NEL0,JJ(NEL0),N
901C REAL
902 my_real
903 . x0(nel0), x1(nel0), x2(nel0), x3(nel0),
904 . y0(nel0), y1(nel0), y2(nel0), y3(nel0),
905 . x(nel0), y(nel0), yp(nel0)
906C-----------------------------------------------
907C L o c a l V a r i a b l e s
908C-----------------------------------------------
909 INTEGER :: I
910C REAL
911 my_real
912 . Q(NEL0), D(NEL0), R(NEL0), S(NEL0), SP(NEL0),
913 . C2(NEL0), DM(NEL0), SM(NEL0), C1(NEL0), C6(NEL0),
914 . c4(nel0), c3(nel0), c3d(nel0), c5(nel0)
915C-----------------------------------------------
916
917 q(1:nel0) = x(1:nel0)-x1(1:nel0)
918 d(1:nel0) = x2(1:nel0)-x1(1:nel0)
919 r(1:nel0) = d(1:nel0)-q(1:nel0)
920 s(1:nel0) = (y2(1:nel0)-y1(1:nel0)) / d(1:nel0)
921 sp(1:nel0) = (y3(1:nel0)-y2(1:nel0)) / (x3(1:nel0)-x2(1:nel0))
922 c2(1:nel0) = (sp(1:nel0)-s(1:nel0)) / (x3(1:nel0)-x1(1:nel0))
923 dm(1:nel0) = x1(1:nel0)-x0(1:nel0)
924 dm(1:nel0) = sign(max(em30,abs(dm(1:nel0))),dm(1:nel0))
925 sm(1:nel0) = (y1(1:nel0)-y0(1:nel0)) / dm(1:nel0)
926 c1(1:nel0) = (s(1:nel0)-sm(1:nel0)) / (d(1:nel0)+dm(1:nel0))
927 c6(1:nel0) = zero
928C
929#include "vectorize.inc"
930 DO i=1,nel0
931 ! -------------------
932 IF(jj(i)==1)THEN
933 IF(x(i)<=x0(i)) THEN
934 c1(i) = zero
935 c2(i) = zero
936 sm(i) = zero
937 ELSE
938 c2(i) = c1(i)
939 c1(i) = sm(i)/(x1(i) - x0(i))
940 ENDIF
941 r(i) = (x1(i) - x(i))
942 q(i) = x(i) - x0(i)
943 d(i) = x1(i) - x0(i)
944 c3(i) = abs(c2(i)*r(i))
945 c3d(i) = c3(i)+abs(c1(i)*q(i))
946 c5(i) = zero
947 IF(c3d(i)>zero) THEN
948 c3(i) = c3(i)/c3d(i)
949 c5(i) = c3(i)*(c1(i)-c2(i))
950 ENDIF
951 c4(i) = c2(i)+c5(i)
952 c6(i) = d(i)*c5(i)*(one - c3(i))
953 y(i) = y0(i) + q(i)*(sm(i)-r(i)*c4(i))
954 yp(i) = sm(i)+(q(i)-r(i))*c4(i)+c6(i)
955 ! -------------------
956 ELSE IF(jj(i)==n-1)THEN
957 IF(sp(i)==zero.OR.(x(i)>x3(i))) THEN
958 c1(i) = zero
959 c2(i) = zero
960 ELSE
961 c1(i) = (sp(i) - s(i))/(x3(i) - x1(i))
962 c2(i) = zero
963 ENDIF
964 r(i) = (x3(i) - x(i))
965 q(i) = x(i) - x2(i)
966 d(i) = x3(i) - x2(i)
967 c3(i) = abs(c2(i)*r(i))
968 c3d(i) = c3(i)+abs(c1(i)*q(i))
969 c5(i) = zero
970 IF(c3d(i)>zero) THEN
971 c3(i) = c3(i)/c3d(i)
972 c5(i) = c3(i)*(c1(i)-c2(i))
973 ENDIF
974 c4(i) = c2(i)+c5(i)
975 c6(i) = d(i)*c5(i)*(one - c3(i))
976 y(i) = y2(i) + (x(i) - x2(i))*(sp(i)-(x3(i) - x(i))*c4(i))
977 yp(i) = sp(i)+(q(i)-r(i))*c4(i)+c6(i)
978 ! -------------------
979 ELSE
980 IF(jj(i)==2.AND.sm(i)*(sm(i)-dm(i)*c1(i))<=zero)c1(i)=(s(i)-sm(i)-sm(i)) / d(i)
981 c3(i) = abs(c2(i)*r(i))
982 c3d(i) = c3(i)+abs(c1(i)*q(i))
983 c5(i) = zero
984 IF(c3d(i)>zero) THEN
985 c3(i) = c3(i)/c3d(i)
986 c5(i) = c3(i)*(c1(i)-c2(i))
987 ENDIF
988 c4(i) = c2(i)+c5(i)
989 c6(i) = d(i)*c5(i)*(one - c3(i))
990 y(i) = y1(i) + q(i)*(s(i)-r(i)*c4(i))
991 yp(i) = s(i)+(q(i)-r(i))*c4(i)+c6(i)
992 ENDIF
993 ! -------------------
994 ENDDO
995 RETURN
996 END SUBROUTINE inter_rat_v
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sigeps60
subroutine sigeps60(nel, nuparam, nuvar, mfunc, kfunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off, ngl, ipt, ipm, mat, epsp, ipla, yld, pla, dpla, amu)
Definition sigeps60.F:48
inter_rat_v
subroutine inter_rat_v(nel0, x0, x1, x2, x3, y0, y1, y2, y3, x, y, yp, jj, n)
Definition sigeps60c.F:890
inter_rat
subroutine inter_rat(x0, x1, x2, x3, y0, y1, y2, y3, x, y, yp, i, n)
Definition sigeps60c.F:786
sigeps60c
subroutine sigeps60c(nel0, nuparam, nuvar, mfunc, kfunc, npf, npt0, ipt, iflag, tf, time, timestep, uparam, rho0, area, eint, thkly, israte, asrate, epspxx, epspyy, epspxy, epspyz, epspzx, depsxx, depsyy, depsxy, depsyz, depszx, epsxx, epsyy, epsxy, epsyz, epszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, sigvxx, sigvyy, sigvxy, sigvyz, sigvzx, soundsp, viscmax, thk, pla, uvar, off, ngl, ipm, mat, etse, gs, yld, epsd_pg, epsp, dpla_i, shf, inloc, dplanl, loff)
Definition sigeps60c.F:47
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72