sigeps60.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    sigeps60    ../engine/source/materials/mat/mat060/sigeps60.F
!||--- called by ------------------------------------------------------
!||    mulaw       ../engine/source/materials/mat_share/mulaw.F90
!||    mulaw8      ../engine/source/materials/mat_share/mulaw8.f90
!||--- calls      -----------------------------------------------------
!||    finter      ../engine/source/tools/curve/finter.F
!||    finter2     ../engine/source/tools/curve/vinter.F
!||    inter_rat   ../engine/source/materials/mat/mat060/sigeps60c.F
!||    vinter      ../engine/source/tools/curve/vinter.F
!||    vinter2     ../engine/source/tools/curve/vinter.F
!||====================================================================
      SUBROUTINE sigeps60(
     1     NEL    ,NUPARAM,NUVAR   ,MFUNC   ,KFUNC   ,NPF    ,
     2     TF     ,TIME   ,TIMESTEP,UPARAM  ,RHO0    ,RHO    ,
     3     VOLUME ,EINT   ,
     4     EPSPXX ,EPSPYY ,EPSPZZ  ,EPSPXY  ,EPSPYZ  ,EPSPZX ,
     5     DEPSXX ,DEPSYY ,DEPSZZ  ,DEPSXY  ,DEPSYZ  ,DEPSZX ,
     6     EPSXX  ,EPSYY  ,EPSZZ   ,EPSXY   ,EPSYZ   ,EPSZX  ,
     7     SIGOXX ,SIGOYY ,SIGOZZ  ,SIGOXY  ,SIGOYZ  ,SIGOZX ,
     8     SIGNXX ,SIGNYY ,SIGNZZ  ,SIGNXY  ,SIGNYZ  ,SIGNZX ,
     9     SIGVXX ,SIGVYY ,SIGVZZ  ,SIGVXY  ,SIGVYZ  ,SIGVZX ,
     A     SOUNDSP,VISCMAX,UVAR    ,OFF     ,NGL     ,IPT    ,
     B     IPM     ,MAT    ,EPSP    ,IPLA    ,YLD     ,PLA   ,
     C     DPLA    ,AMU    )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "comlock.inc"
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "mvsiz_p.inc"
C---------+---------+---+---+--------------------------------------------
C VAR     | SIZE    |TYP| RW| DEFINITION
C---------+---------+---+---+--------------------------------------------
C NEL     |  1      | I | R | SIZE OF THE ELEMENT GROUP NEL 
C NUPARAM |  1      | I | R | SIZE OF THE USER PARAMETER ARRAY
C NUVAR   |  1      | I | R | NUMBER OF USER ELEMENT VARIABLES
C---------+---------+---+---+--------------------------------------------
C MFUNC   |  1      | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
C KFUNC   | NFUNC   | I | R | FUNCTION INDEX not used
C NPF     |  *      | I | R | FUNCTION ARRAY   
C TF      |  *      | F | R | FUNCTION ARRAY 
C---------+---------+---+---+--------------------------------------------
C TIME    |  1      | F | R | CURRENT TIME
C TIMESTEP|  1      | F | R | CURRENT TIME STEP
C UPARAM  | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
C RHO0    | NEL     | F | R | INITIAL DENSITY
C RHO     | NEL     | F | R | DENSITY
C VOLUME  | NEL     | F | R | VOLUME
C EINT    | NEL     | F | R | TOTAL INTERNAL ENERGY
C EPSPXX  | NEL     | F | R | STRAIN RATE XX
C EPSPYY  | NEL     | F | R | STRAIN RATE YY
C ...     |         |   |   |
C DEPSXX  | NEL     | F | R | STRAIN INCREMENT XX
C DEPSYY  | NEL     | F | R | STRAIN INCREMENT YY
C ...     |         |   |   |
C EPSXX   | NEL     | F | R | STRAIN XX
C EPSYY   | NEL     | F | R | STRAIN YY
C ...     |         |   |   |
C SIGOXX  | NEL     | F | R | OLD ELASTO PLASTIC STRESS XX 
C SIGOYY  | NEL     | F | R | OLD ELASTO PLASTIC STRESS YY
C ...     |         |   |   |    
C---------+---------+---+---+--------------------------------------------
C SIGNXX  | NEL     | F | W | NEW ELASTO PLASTIC STRESS XX
C SIGNYY  | NEL     | F | W | NEW ELASTO PLASTIC STRESS YY
C ...     |         |   |   |
C SIGVXX  | NEL     | F | W | VISCOUS STRESS XX
C SIGVYY  | NEL     | F | W | VISCOUS STRESS YY
C ...     |         |   |   |
C SOUNDSP | NEL     | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
C VISCMAX | NEL     | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
C---------+---------+---+---+--------------------------------------------
C UVAR    |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
C OFF     | NEL     | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
C---------+---------+---+---+--------------------------------------------
#include      "param_c.inc"
#include      "parit_c.inc"
#include      "scr17_c.inc"
#include      "com01_c.inc"
#include      "com08_c.inc"
#include      "units_c.inc"
C
      INTEGER NEL, NUPARAM, NUVAR,IPT,
     .   NGL(NEL),MAT(NEL),IPLA, IPM(NPROPMI,*)
      my_real
     .   TIME,TIMESTEP,UPARAM(*),
     .   RHO(NEL),RHO0(NEL),VOLUME(NEL),EINT(NEL),
     .   EPSPXX(NEL),EPSPYY(NEL),EPSPZZ(NEL),
     .   EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
     .   DEPSXX(NEL),DEPSYY(NEL),DEPSZZ(NEL),
     .   DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
     .   EPSXX(NEL) ,EPSYY(NEL) ,EPSZZ(NEL) ,
     .   EPSXY(NEL) ,EPSYZ(NEL) ,EPSZX(NEL) ,
     .   SIGOXX(NEL),SIGOYY(NEL),SIGOZZ(NEL),
     .   sigoxy(nel),sigoyz(nel),sigozx(nel),
     .   epsp(nel)  ,amu(nel)
C-----------------------------------------------
C   O U T P U T   A r g u m e n t s
C-----------------------------------------------
      my_real
     .    signxx(nel),signyy(nel),signzz(nel),
     .    signxy(nel),signyz(nel),signzx(nel),
     .    sigvxx(nel),sigvyy(nel),sigvzz(nel),
     .    sigvxy(nel),sigvyz(nel),sigvzx(nel),
     .    soundsp(nel),viscmax(nel),dpla(nel)
C-----------------------------------------------
C   I N P U T   O U T P U T   A r g u m e n t s 
C-----------------------------------------------
      my_real 
     .        uvar(nel,nuvar), off(nel),  yld(nel),pla(nel)
C-----------------------------------------------
C   VARIABLES FOR FUNCTION INTERPOLATION 
C-----------------------------------------------
      INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
      my_real
     .   FINTER2, TF(*),FINTER
      EXTERNAL FINTER2,FINTER
C     EXTERNAL FINTER, FINTER2
C        Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
C        Y       : y = f(x)
C        X       : x
C        DYDX    : f'(x) = dy/dx
C        IFUNC(J): FUNCTION INDEX
C              J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
C        NPF,TF  : FUNCTION PARAMETER
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,J,IADBUFV,JJ(MVSIZ),NFUNC,
     .        NRATE,IPOS1(MVSIZ),IPOS2(MVSIZ),IAD1(MVSIZ),
     .        ILEN1(MVSIZ),IAD2(MVSIZ),ILEN2(MVSIZ),
     .        IFUNC(100),NFUNCV,PFUN,
     .        IPOSP(MVSIZ),IADP(MVSIZ),ILENP(MVSIZ),NFUNCM,NRATEM,
     .        IPFLAG,IPARAM,NPAR,IPOS3(MVSIZ),IPOS4(MVSIZ),IAD3(MVSIZ),
     .        ILEN3(MVSIZ),IAD4(MVSIZ),ILEN4(MVSIZ),J1, J2, J3, J4,
     .        NINDX,INDX(MVSIZ), NN,OPTE,IFUNCE,MX
      my_real
     .        E(MVSIZ),NU,P,DAV,VM,R,FAC,EPST,EP1,EP2,EP3,EP4,EP5,EP6,
     .        E1,E2,E3,E4,E5,E6,C,CC,D,Y,YP,E42,E52,E62,EPST2,
     .        C1(MVSIZ),G(MVSIZ),G2(MVSIZ),G3(MVSIZ),
     .        epsmax,rate(mvsiz,4),yfac(mvsiz,4),
     .        y1(mvsiz),y2(mvsiz),h(mvsiz),dydx1(mvsiz),
     .        dydx2(mvsiz),fail(mvsiz),epsr1,
     .        epsr2,p0(mvsiz),pfac(mvsiz),escale(mvsiz),
     .        dfdp(mvsiz),einf,ce,dydxe(mvsiz),
     .        y11, y22, y33, y44,x,
     .        x1, x2, x3, x4,y3(mvsiz),y4(mvsiz),
     .        dydx3(mvsiz),dydx4(mvsiz),pscale(mvsiz)
C-----------------------------------------------
C     USER VARIABLES INITIALIZATION
C-----------------------------------------------
      IF (ivector/=1) THEN
       mx = mat(1)
       nfunc  = ipm(10,mx)
       DO j=1,nfunc
        ifunc(j)=ipm(10+j,mx)        
       ENDDO
       pfun= ifunc(nfunc-1)  !! -1
C
      ELSE
C ivector = 1 only
       mx = mat(1)
       nfuncm = 0
       nfuncv = ipm(10,mx)
       nfuncm = max(nfuncm,nfuncv)
       DO j=1,nfuncm
            IF(nfuncv>=j) THEN
              ifunc(j)=ipm(10+j,mx)
            ENDIF
       ENDDO
       pfun= ifunc(nfuncv-1) !  -1
      ENDIF
C
      nratem = 0
C
C
      mx = mat(1)
      iadbufv = ipm(7,mx)-1
      nu  = uparam(iadbufv+6)
C
      nrate = nint(uparam(iadbufv+1))
      nratem = max(nratem,nrate)
     
      epsmax=uparam(iadbufv+6+2*nrate+1)
      IF(epsmax==zero)THEN
       IF(tf(npf(ifunc(1)+1)-1)==zero)THEN
         epsmax=tf(npf(ifunc(1)+1)-2)
       ELSE
         epsmax= ep30
       ENDIF
      ENDIF
c ---
      opte    = uparam(iadbufv+2*nrate+18)
      einf    = uparam(iadbufv+2*nrate+19)
      ce      = uparam(iadbufv+2*nrate+20)
      epsr1 =uparam(iadbufv+6+2*nrate+2)
      IF(epsr1==zero)epsr1= ep30
      epsr2 =uparam(iadbufv+6+2*nrate+3)
      IF(epsr2==zero)epsr2= twoep30
      DO i=1,nel
        e(i)   = uparam(iadbufv+2)
        g(i)   = uparam(iadbufv+5)
        g2(i)  = two*g(i)
        g3(i)  = three*g(i)
        c1(i)  = e(i)/three/(one - two*nu)
        pscale(i)= uparam(iadbufv+16+2*nrate)
c ---
        pla(i) = uvar(i,1)
       ENDDO        
       IF (opte == 1)THEN
        DO i=1,nel      
          IF(pla(i) > zero)THEN    
             ifunce = uparam(iadbufv+2*nrate+17)                                 
             escale(i) = finter(kfunc(ifunce),pla(i),npf,tf,dydxe(i))          
             e(i) =  escale(i)* e(i)                                       
             g(i) =  half*e(i)/(one+nu) 
             g2(i) = two*g(i)                                  
             g3(i) = three*g(i) 
             c1(i) =e(i)/three/(one - two*nu)                                              
             soundsp(i) = sqrt((c1(i) + four_over_3*g(i))/rho0(i))                
           ENDIF
         ENDDO           
       ELSEIF ( ce /= zero) THEN
        DO i=1,nel       
           IF(pla(i) > zero)THEN                                                        
             e(i) = e(i)-(e(i)-einf)*(one-exp(-ce*pla(i)))                      
             g(i) =  half*e(i)/(one+nu)   
             g2(i) = two*g(i)                                  
             g3(i) = three*g(i)                                               
             c1(i) =e(i)/three/(one - two*nu)                                              
             soundsp(i) = sqrt((c1(i) + four_over_3*g(i))/rho0(i))                               
           ENDIF    
        ENDDO           
      ENDIF
C +++
      ipflag = 0
      DO i=1,nel
        pfac(i) = one
          IF (pfun>0) ipflag = ipflag + 1
      ENDDO
C ---
C
C
      IF (isigi==0) THEN
      IF(time==0.0)THEN
        DO i=1,nel
         uvar(i,1)=zero
         uvar(i,2)=zero
         DO j=1,nrate
           uvar(i,j+2)=0
         ENDDO
         IF (pfun>0) uvar(i,nrate+5)=0
        ENDDO
      ENDIF
      ENDIF
C-----------------------------------------------
C
      DO i=1,nel
C
        pla(i) = uvar(i,1)
        p0(i) = -(sigoxx(i)+sigoyy(i)+sigozz(i))*third
        dav = (depsxx(i)+depsyy(i)+depszz(i))*third
        signxx(i)=sigoxx(i)+p0(i)+g2(i)*(depsxx(i)-dav)
        signyy(i)=sigoyy(i)+p0(i)+g2(i)*(depsyy(i)-dav)
        signzz(i)=sigozz(i)+p0(i)+g2(i)*(depszz(i)-dav)
        signxy(i)=sigoxy(i)+g(i) *depsxy(i)
        signyz(i)=sigoyz(i)+g(i) *depsyz(i)
        signzx(i)=sigozx(i)+g(i) *depszx(i)
        pscale(i) = pscale(i)*p0(i)
        soundsp(i) = sqrt((c1(i)+four_over_3*g(i))/rho0(i))
        viscmax(i) = zero
      ENDDO
C-------------------
C     STRAIN RATE
C-------------------
      DO i=1,nel
C-------------------
C     STRAIN principal 1, 4 newton iterations 
C-------------------
        dav = (epsxx(i)+epsyy(i)+epszz(i))*third
        e1 = epsxx(i) - dav
        e2 = epsyy(i) - dav
        e3 = epszz(i) - dav
        e4 = 0.5*epsxy(i)
        e5 = 0.5*epsyz(i)
        e6 = 0.5*epszx(i)
C        -y = (e1-x)(e2-x)(e3-x)
C           - e5^2(e1-x) - e6^2(e2-x) - e4^2(e3-x)
C           + 2e4 e5 e6
C         e1 + e2 + e3 = 0 => terme en x^2 = 0
C         y = x^3 + c x + d
c         yp= 3 x^2 + c
        e42 = e4*e4
        e52 = e5*e5
        e62 = e6*e6
        c = - e1*e1 - e2*e2 - e3*e3 - e42 - e52 - e62
        d = - e1*e2*e3 + e1*e52 + e2*e62 + e3*e42
     &      - two*e4*e5*e6 
        cc = c*third
        epst = sqrt(-cc)
        epst2 = epst * epst
        y = (epst2 + c)* epst + d
        IF(abs(y)>em08)THEN
          epst = onep75 * epst
          epst2 = epst * epst
          y = (epst2 + c)* epst + d
          yp = three*epst2 + c
          IF(yp/=zero)epst = epst - y/yp
          epst2 = epst * epst
          y = (epst2 + c)* epst + d
          yp = three*epst2 + c
          IF(yp/=zero)epst = epst - y/yp
          epst2 = epst * epst
          y = (epst2 + c)* epst + d
          yp = three*epst2 + c
          IF(yp/=zero)epst = epst - y/yp
          epst2 = epst * epst
          y = (epst2 + c)* epst + d
          yp = three*epst2 + c
          IF(yp/=zero)epst = epst - y/yp
        ENDIF
          epst = epst + dav
C-------------------
C     tension failure 
C-------------------
        fail(i) = max(zero,min(one,
     .            (epsr2-epst)/(epsr2-epsr1) ))
      ENDDO
C-------------------
C     CRITERE
C-------------------
      DO i=1,nel
            jj(i) = 1
      ENDDO
      IF (ivector/=1) THEN
       DO i=1,nel
        DO j=2,nrate-1
          IF(epsp(i)>=uparam(iadbufv+6+j)) jj(i) = j
        ENDDO
       ENDDO
C
      ELSE
        DO j=2,nratem-1
          DO i=1,nel
            IF(nrate-1>=j) THEN
              IF(epsp(i)>=uparam(iadbufv+6+j)) jj(i) = j
            ENDIF
          ENDDO
        ENDDO
      ENDIF
C
      DO i=1,nel   
         IF(jj(i)==1)THEN 
           j1 = jj(i)
           j2 = j1+1
           j3 = j2+1
           j4 = j3+1 
          ELSEIF(jj(i)==(nrate-1))THEN
           j1 = jj(i) - 2
           j2 = j1+1 
           j3 = j2+1 
           j4 = j3+1 
          ELSE
           j1 = jj(i) - 1
           j2 = j1+1
           j3 = j2+1
           j4 = j3+1 
          ENDIF        
        rate(i,1)=uparam(iadbufv + 6 +  j1 )
        rate(i,2)=uparam(iadbufv + 6 +  j2 )
        rate(i,3)=uparam(iadbufv + 6 +  j3 )
        rate(i,4)=uparam(iadbufv + 6 +  j4 )
        yfac(i,1)=uparam(iadbufv+6+nrate + j1 )     
        yfac(i,2)=uparam(iadbufv+6+nrate + j2 )
        yfac(i,3)=uparam(iadbufv+6+nrate + j3 )
        yfac(i,4)=uparam(iadbufv+6+nrate + j4 )         
      ENDDO
C
C
 
      DO i=1,nel 
         IF(jj(i)==1)THEN 
           j1 = jj(i)
           j2 = j1+1
           j3 = j2+1
           j4 = j3+1 
          ELSEIF(jj(i)==(nrate-1))THEN
           j1 = jj(i) - 2
           j2 = j1+1 
           j3 = j2+1 
           j4 = j3+1 
          ELSE
           j1 = jj(i) - 1
           j2 = j1+1
           j3 = j2+1
           j4 = j3+1 
          ENDIF
        ipos1(i) = nint(uvar(i,j1+2))
        iad1(i)  = npf(ifunc(j1)) / 2 + 1
        ilen1(i) = npf(ifunc(j1)+1) / 2 - iad1(i) - ipos1(i)
        ipos2(i) = nint(uvar(i,j2+2))
        iad2(i)  = npf(ifunc(j2)) / 2 + 1
        ilen2(i) = npf(ifunc(j2)+1) / 2 - iad2(i) - ipos2(i)
        ipos3(i) = nint(uvar(i,j3+4))
        iad3(i)  = npf(ifunc(j3)) / 2 + 1
        ilen3(i) = npf(ifunc(j3)+1) / 2 - iad3(i) - ipos3(i)
        ipos4(i) = nint(uvar(i,j4+4))
        iad4(i)  = npf(ifunc(j4)) / 2 + 1
        ilen4(i) = npf(ifunc(j4)+1) / 2 - iad4(i) - ipos4(i)   
      ENDDO
C
      CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
      CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
      CALL vinter(tf,iad3,ipos3,ilen3,nel,pla,dydx3,y3)
      CALL vinter(tf,iad4,ipos4,ilen4,nel,pla,dydx4,y4)
C
C
C
C---  PRESSURE DEPENDENT YIELD FUNCTION FACTOR
      IF (ipflag==nel.AND.(iparit==0)) THEN
C------- to optimize in SPMD & PARIT/ON case add a global flag 
C        in a case of homogeneous element groups (to do if needed)
        DO i=1,nel
          iposp(i) = nint(uvar(i,nrate+5))
          iadp(i)  = npf(pfun) / 2 + 1
          ilenp(i) = npf(pfun+1) / 2 - iadp(i) - iposp(i)
        ENDDO
        CALL vinter2(tf,iadp,iposp,ilenp,nel,pscale,dfdp,pfac)
        DO i=1,nel
          uvar(i,nrate+5) = iposp(i)
        ENDDO
      ELSEIF (ipflag>0) THEN
        DO i=1,nel
          IF (pfun>0) THEN
            iposp(i) = nint(uvar(i,nrate+5))
            iadp(i)  = npf(pfun) / 2 + 1
            ilenp(i) = npf(pfun+1) / 2 - iadp(i)-iposp(i)
            pfac(i)  = finter2(tf   ,iadp(i),iposp(i),ilenp(i),
     .                         pscale(i),dfdp(i))
            uvar(i,nrate+5) = iposp(i)
          ENDIF
        ENDDO
      ENDIF
C---
      DO i=1,nel
        IF(jj(i)==1)THEN 
         j1 = jj(i)
         j2 = j1+1
         j3 = j2+1
         j4 = j3+1
         dydx1(i)= dydx1(i)*yfac(i,1)
         dydx2(i) = dydx2(i)*yfac(i,2)          
         fac   = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1)) 
         h(i)   = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
        ELSEIF(jj(i)==(nrate-1))THEN
         j1 = jj(i) - 2
         j2 = j1+1 
         j3 = j2+1 
         j4 = j3+1
         dydx3(i) = dydx3(i)*yfac(i,3)
         dydx4(i) = dydx4(i)*yfac(i,4)
         fac   = (epsp(i) - rate(i,3))/(rate(i,4) - rate(i,3)) 
         h(i)   = fail(i)*(dydx3(i) + fac*(dydx4(i)-dydx3(i)))  
        ELSE
         j1 = jj(i) - 1
         j2 = j1+1
         j3 = j2+1
         j4 = j3+1 
         dydx2(i) = dydx2(i)*yfac(i,2)
         dydx3(i) = dydx3(i)*yfac(i,3)
         fac   = (epsp(i) - rate(i,2))/(rate(i,3) - rate(i,2))  
         h(i)   = fail(i)*(dydx2(i) + fac*(dydx3(i)-dydx2(i)))  
        ENDIF      
        nn = nrate
        x1 = rate(i,1)
        x2 = rate(i,2)
        x3 = rate(i,3)
        x4 = rate(i,4)
        x  = epsp(i)
        y11 = y1(i)*yfac(i,1)
        y22 = y2(i)*yfac(i,2)
        y33 = y3(i)*yfac(i,3)
        y44 = y4(i)*yfac(i,4)
        CALL  inter_rat(x1,x2,x3,x4,y11,y22,y33,y44,
     .                                           x,y,yp,jj(i),nn)
        yld(i) = fail(i)*y
Cc        Y11 = DYDX1(I)*YFAC(I,1)
Cc        Y22 = DYDX2(I)*YFAC(I,2)
Cc        Y33 = DYDX3(I)*YFAC(I,3)
Cc        Y44 = DYDX4(I)*YFAC(I,4)
Cc        CALL  INTER_RAT(X1,X2,X3,X4,Y11,Y22,Y33,Y44,
Cc     .                                           X,Y,YP,JJ(I),NN)
Cc        H(I)   = FAIL(I)*Y         
C
        yld(i) = yld(i) * max(zero, pfac(i))
        uvar(i,j1+2) = ipos1(i)
        uvar(i,j2+2) = ipos2(i)
        uvar(i,j3+2) = ipos3(i)
        uvar(i,j4+2) = ipos4(i)
      ENDDO
C-------------------
C     PROJECTION
C-------------------
      IF(ipla==0)THEN
       DO i=1,nel
        vm =half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
     1          +signxy(i)**2+signyz(i)**2+signzx(i)**2
        vm =sqrt(three*vm)
        r = min(one,yld(i)/ max(vm,em20))
c        P = C1(I) * (RHO(I)/RHO0(I) - ONE)
        p = c1(i) * amu(i)
        signxx(i)=signxx(i)*r-p
        signyy(i)=signyy(i)*r-p
        signzz(i)=signzz(i)*r-p
        signxy(i)=signxy(i)*r
        signyz(i)=signyz(i)*r
        signzx(i)=signzx(i)*r
        pla(i)=pla(i) + (one - r)*vm/max(g3(i)+h(i),em20)
        uvar(i,1)=pla(i)
        dpla(i) =   (one - r)*vm/max(g3(i)+h(i),em20)       
       ENDDO
      ELSEIF(ipla==2)THEN
       DO i=1,nel
        vm =half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
     1          +signxy(i)**2+signyz(i)**2+signzx(i)**2
        vm =sqrt(three*vm)
        r = min(one,yld(i)/ max(vm,em20))
c        P = C1(I) * (RHO(I)/RHO0(I)- ONE)
        p = c1(i) * amu(i)
        signxx(i)=signxx(i)*r-p
        signyy(i)=signyy(i)*r-p
        signzz(i)=signzz(i)*r-p
        signxy(i)=signxy(i)*r
        signyz(i)=signyz(i)*r
        signzx(i)=signzx(i)*r
        pla(i)=pla(i) + (one - r)*vm/max(g3(i),em20)
        uvar(i,1)=pla(i)
        dpla(i) = (one - r)*vm/max(g3(i),em20)        
       ENDDO
      ELSEIF(ipla==1)THEN
       DO i=1,nel
        vm =half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
     1          +signxy(i)**2+signyz(i)**2+signzx(i)**2
        vm =sqrt(three*vm)
        r = min(one,yld(i)/ max(vm,em20))
C     plastic strain increment.
       dpla(i) = (one - r)*vm/max(g3(i)+h(i),em20)
C     actual yield stress.
         yld(i)= max(yld(i)+dpla(i)*h(i),zero)
        r = min(one,yld(i)/ max(vm,em20))
C
c        P = C1(I) * (RHO(I)/RHO0(I) - ONE)
        p = c1(i) * amu(i)
        signxx(i)=signxx(i)*r-p
        signyy(i)=signyy(i)*r-p
        signzz(i)=signzz(i)*r-p
        signxy(i)=signxy(i)*r
        signyz(i)=signyz(i)*r
        signzx(i)=signzx(i)*r
       pla(i)=pla(i) + dpla(i)
        uvar(i,1)=pla(i)
       ENDDO
      ENDIF
C---
C-------------------
 
      DO i=1,nel
        IF(off(i)<em01) off(i)=zero
        IF(off(i)<one) off(i)=off(i)*four_over_5
      ENDDO
C
C
C
      nindx=0
      DO i=1,nel
        IF(pla(i)>epsmax.AND.off(i)==one) THEN
          off(i)= four_over_5
          nindx=nindx+1
          indx(nindx)=i
        ENDIF
      ENDDO
      IF(nindx>0)THEN
        DO j=1,nindx
#include "lockon.inc"
        WRITE(iout, 1000) ngl(indx(j))
        WRITE(istdo,1100) ngl(indx(j)),tt
#include "lockoff.inc"
        ENDDO
      ENDIF
 
 1000 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10)
 1100 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10,
     .          ' AT TIME :',g11.4)
C
C
      RETURN
      END
C