Sources Esope | Cast3M

ddsirank
C DDSIRANK  SOURCE    CB215821  16/04/21    21:16:19     8920
      SUBROUTINE ddsigRank(sigma,d2gt,i1,i2,i3,i4,i5,i6,lcp,
     .                     parahot3,idimpara3)
 
c     This subroutine calculates the second derivative of the Rankine
c     plastic function with respect to sigma
c     gt = ft = S1 = largest principal stress
 
      IMPLICIT REAL*8 (A-B,D-H,O-Z)
      implicit integer (I-K,M,N)
      implicit logical (L)
      implicit character*10 (C)
 
****  dimension sigma(i6),dgt(i6),d2gt(i6,i6)
      dimension sigma(i6),dgt(6),d2gt(i6,i6)
****  dimension rcossigmapr(i3,i3),P1(3,3),pi1(3)
      dimension rcossigmapr(3,3),P1(3,3),pi1(3)
      dimension sigmacp(3),vloc(3),vloc2(1,3),vloc1(1)
      dimension vloc3(3,3),d2gtcp(3,3),parahot3(idimpara3)
      dimension sigmaim(6),sigmail(6),rhi(6),dgtm(6),dgtl(6)
 
 
 
*****
*     MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS)
      IF(I3.GT.3) PRINT *, ' DDSIGRANK - ERREUR  I2 = ', I2, ' > 2 '
      IF(I6.GT.6) PRINT *, ' DDSIGRANK - ERREUR  I6 = ', I6, ' > 6 '
*****
 
 
      if (lcp) then
c       In plane stress, the second derivative can be calculated analytically
        do jloc=1,6
          do iloc=1,6
            d2gt(iloc,jloc) = 0.d0
          enddo
        enddo
        P1(1,1) = 0.5d0
        P1(2,2) = 0.5d0
        P1(1,2) = -0.5d0
        P1(2,1) = -0.5d0
        do iloc=1,2
          P1(iloc,3) = 0.d0
          P1(3,iloc) = 0.d0
        enddo
        P1(3,3) = 2.d0
        pi1(1) = 1.d0
        pi1(2) = 1.d0
        pi1(3) = 0.d0
        sigmacp(1) = sigma(1)
        sigmacp(2) = sigma(2)
        sigmacp(3) = sigma(4)
        call mulAB(P1,sigmacp,vloc,3,3,3,1)
        call mulATB(sigmacp,vloc,vloc1,3,1,3,1)
        rloc = 0.5d0 * vloc1(1)
        psi1 = SQRT(rloc)
        call mulATB(sigmacp,P1,vloc2,3,1,3,3)
        call mulAB(vloc,vloc2,vloc3,3,1,1,3)
 
        test = 0.5d0*(sigmacp(1)+sigmacp(2))
        psi1lim=max(1.E-7*test,1.d-3)
        if (ABS(psi1).le.psi1lim) then
          continue
c         d2gt = 0
        else
          do jloc=1,3
            do iloc=1,3
              d2gtcp(iloc,jloc)=(P1(iloc,jloc)/(2.0d0*psi1))
     .                 - (vloc3(iloc,jloc)/(4.0d0*psi1*psi1*psi1))
            enddo
          enddo
          do iloc=1,2
            d2gt(1,iloc) = d2gtcp(1,iloc)
            d2gt(2,iloc) = d2gtcp(2,iloc)
            d2gt(4,iloc) = d2gtcp(3,iloc)
            d2gt(iloc,4) = d2gtcp(iloc,3)
          end do
          d2gt(4,4) = d2gtcp(3,3)
        endif
      else
c       In full 3D, numerical differentiation
        do jloc=1,6
          do iloc=1,6
            d2gt(iloc,jloc) = 0.d0
          enddo
        enddo
        rh = 1.0d-5
        do iloc=1,6
c         ith component of sigma => ith column of d2gt
c         calculation of   sigma + hi * ei   and   sigma - hi * ei
          do jloc=1,6
            rhi(jloc) = 0.d0
          enddo
          rhi(iloc) = 1.0d0 * rh
          do jloc=1,6
            sigmaim(jloc) = sigma(jloc) + rhi(jloc)
            sigmail(jloc) = sigma(jloc) - rhi(jloc)
          enddo
c         Calculation of q(x+hi*ei) = dft(x+hi*ei)/dsig
          call dsigRank(sigmaim,dgtm,3,6,parahot3,idimpara3,
     .                  rkappait,lcp)
          call dsigRank(sigmail,dgtl,3,6,parahot3,idimpara3,
     .                  rkappait,lcp)
          do jloc=1,6
            d2gt(jloc,iloc) = (dgtm(jloc) - dgtl(jloc))/(2.0d0*rh)
          enddo
        enddo
c       Calculation of f(x+hi*ei)
      endif
 
      RETURN
      END