Sources Esope | Cast3M

dsiggc
C DSIGGC    SOURCE    CB215821  16/04/21    21:16:32     8920
      SUBROUTINE dsiggc(sigma,dgc,ag,P2,pi2,i1,i6,fc,fb,lcp)
 
c     This subroutine calculates the derivative of the plastic potential in compression
c     (necessary for the calculation of the plastic strain in compression)
c
c     Plastic potential: gc = SQRT(3 J2) + ag I1  with  J2 : 2nd invariant of the deviatoric stress
c                                                   and   I1 : 1st invariant of the stress
 
      IMPLICIT REAL*8 (A-B,D-H,O-Z)
      implicit integer (I-K,M,N)
      implicit logical (L)
      implicit character*10 (C)
 
****  dimension vloc(i6),vloc1(i1),pi2(i6),vloc2(i6)
      dimension vloc(6),vloc1(1),pi2(I6),vloc2(6)
      dimension dgc(i6),p2(i6,i6),sigma(i6)
 
 
 
*****
*     MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS)
      IF(I1.GT.1) PRINT *, ' DSIGGC - ERREUR  I1 = ', I1, ' > 1 '
      IF(I6.GT.6) PRINT *, ' DSIGGC - ERREUR  I6 = ', I6, ' > 6 '
*****
 
 
 
 
      r1 = 1.
      r2 = 2.
      alpha  = (fb-fc)/(2*fb-fc)
 
      call mulAB(P2,sigma,vloc,i6,i6,i6,i1)
c          [6;1]=[6;6]x[6;1]
      do iloc=1,6
        vloc2(iloc) = vloc(iloc)
      end do
      call mulATB(sigma,vloc,vloc1,i6,i1,i6,i1)
c          [1]=[1;6]x[6;1]
      rloc = (r1/r2)* vloc1(i1)
 
      if (rloc.le.0.d0) then
        psi2 = 0.d0
      else
        psi2 = SQRT(rloc)
c       psi2 = SQRT(3 * J2)
      endif
 
      call mulATB(pi2,sigma,vloc1,i6,i1,i6,i1)
c          [1]=[1;6]x[6;1]
      test = alpha*vloc1(i1)
c     test = alpha * I1
 
      psi2lim=max(ABS(1.E-7*test),1.d-3)
      f_compr = DPfunc(sigma,alpha)
      if (ABS(psi2).le.psi2lim) then
        do iloc=i1,i6
          dgc(iloc) = ((vloc2(iloc)/(2.*psi2lim)) + ag * pi2(iloc))
        end do
      else if (ABS(f_compr).le.100.) then
c       APEX of the Drucker-Prager surface
        do iloc=i1,i6
          dgc(iloc) = ((vloc2(iloc)/(2.*psi2)) + ag * pi2(iloc))
        end do
      else
        do iloc=i1,i6
          dgc(iloc) = ((vloc2(iloc)/(2.*psi2)) + ag * pi2(iloc))
        end do
      endif
 
      RETURN
      END