Esope source files

   1 : C DSIGGC    SOURCE    AM        12/02/28    00:04:15     7289           
   2 :       SUBROUTINE dsiggc(sigma,dgc,ag,P2,pi2,i1,i6,fc,fb,lcp)
   3 : 
   4 : c     This subroutine calculates the derivative of the plastic potential in compression
   5 : c     (necessary for the calculation of the plastic strain in compression)                                               
   6 : c                                           
   7 : c     Plastic potential: gc = dsqrt (3 J2) + ag I1  with  J2 : 2nd invariant of the deviatoric stress                                                 
   8 : c                                                   and   I1 : 1st invariant of the stress
   9 : 
  10 :       IMPLICIT REAL*8 (A-B,D-H,O-Z)
  11 :       implicit integer (I-K,M,N)
  12 :       implicit logical (L)
  13 :       implicit character*10 (C)      
  14 : 
  15 : ****  dimension vloc(i6),vloc1(i1),pi2(i6),vloc2(i6)
  16 :       dimension vloc(6),vloc1(1),pi2(I6),vloc2(6)
  17 :       dimension dgc(i6),p2(i6,i6),sigma(i6)
  18 : 
  19 : 
  20 : 
  21 : *****
  22 : *     MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS)
  23 :       IF(I1.GT.1) PRINT *, ' DSIGGC - ERREUR  I1 = ', I1, ' > 1 '
  24 :       IF(I6.GT.6) PRINT *, ' DSIGGC - ERREUR  I6 = ', I6, ' > 6 '
  25 : *****
  26 : 
  27 : 
  28 : 
  29 : 
  30 :       r1 = 1.
  31 :       r2 = 2.
  32 :       alpha  = (fb-fc)/(2*fb-fc)
  33 :       
  33 :       call mulAB(P2,sigma,vloc,i6,i6,i6,i1)             
  34 : c          [6;1]=[6;6]x[6;1]
  35 :       do iloc=1,6
  36 :         vloc2(iloc) = vloc(iloc)
  37 :       end do 
  37 :       call mulATB(sigma,vloc,vloc1,i6,i1,i6,i1)         
  38 : c          [1]=[1;6]x[6;1]
  39 :       rloc = (r1/r2)* vloc1(i1)
  40 : 
  41 :       if (rloc.le.0.d0) then
  42 :         psi2 = 0.d0
  43 :       else
  44 :         psi2 = dsqrt(rloc)
  45 : c       psi2 = dsqrt (3 * J2)
  46 :       endif
  47 :      
  47 :       call mulATB(pi2,sigma,vloc1,i6,i1,i6,i1)          
  48 : c          [1]=[1;6]x[6;1]
  49 :       test = alpha*vloc1(i1)
  50 : c     test = alpha * I1
  51 : 
  52 :       psi2lim=max(dabs(1.E-7*test),1.d-3)
  53 :       f_compr = DPfunc(sigma,alpha)
  54 :       if (dabs(psi2).le.psi2lim) then
  55 :         do iloc=i1,i6
  56 :           dgc(iloc) = ((vloc2(iloc)/(2.*psi2lim)) + ag * pi2(iloc))
  57 :         end do
  58 :       else if (dabs(f_compr).le.100.) then      
  59 : c       APEX of the Drucker-Prager surface
  60 :         do iloc=i1,i6
  61 :           dgc(iloc) = ((vloc2(iloc)/(2.*psi2)) + ag * pi2(iloc))
  62 :         end do
  63 :       else
  64 :         do iloc=i1,i6
  65 :           dgc(iloc) = ((vloc2(iloc)/(2.*psi2)) + ag * pi2(iloc))
  66 :         end do
  67 :       endif
  68 : 
  69 :       RETURN
  70 :       END
  71 : 
  72 :