Sources Esope | Cast3M

sigmapr3
C SIGMAPR3  SOURCE    CB215821  16/04/21    21:18:26     8920
        subroutine sigmapr3ETC(sigma,S1,S2,S3,rcossigmapr,lcp)
 
c       ******************************************************************
c       *                                                                *
c       *     Principal stresses ( Batoz, Vol. 1., 1.2.5. p 51)          *
c       *     =================================================          *
c       *                                                                *
c       ******************************************************************
 
c       *     sigma   : vector of the stresses in the global axis
c                     : Sxx, Syy, Szz, Sxy, Sxz, Syz
c       * S1          : largest principal      stress
c       * S2          : intermediate principal stress
c       * S3          : smallest principal     stress
c       * rcossigmapr :   ( lx      ly      lz )
c                         ( mx      my      mz )
c                         ( nx      ny      nz )
c               for         S1      S2      S3
 
      IMPLICIT REAL*8 (A-B,D-H,O-Z)
      implicit integer   (I-K,M,N)
      implicit logical (L)
      implicit character*10 (C)
 
      character*20 cloc1,cloc2
 
      dimension sigma(6)
      dimension rcossigmapr(3,3)
 
      i1 = 1
      i2 = 2
      i3 = 3
      i4 = 4
      i5 = 5
      i6 = 6
      i7 = 7
      i8 = 8
      i9 = 9
      i10=10
      i11=11
      i12=12
      i13=13
      i14=14
      i15=15
      i16=16
      i17=17
      i18=18
      i19=19
      i20=20
      i21=21
      i22=22
      i23=23
      i24=24
      i25=25
      i26=26
      i27=27
      i28=28
      i29=29
      i30=30
      i31=31
      i32=32
      i33=33
      i34=34
      i35=35
      i36=36
      i37=37
      i38=38
      i39=39
 
      r0 = 0.
      r1 = 1.
      r2 = 2.
      r3 = 3.
      r6 = 6.
 
      precision = 1.0d-3
      precision2 = precision*precision
      precision3 = precision2*precision
      pi = ACOS(-r1)
      racine2 = SQRT(r2)
      S1 = r0
      S2 = r0
      S3 = r0
 
      do jloc=i1,i3
        do iloc=i1,i3
          RCOSSIGMAPR(iloc,jloc) = r0
        end do
      end do
 
      Sxx = sigma(i1)
      Syy = sigma(i2)
      Szz = sigma(i3)
      Sxy = sigma(i4)
      Sxz = sigma(i5)
      Syz = sigma(i6)
 
      if (lcp) then
c       PLANE STRESS
        sxx = sigma(1)
        syy = sigma(2)
        sxy = sigma(4)
        diff = sxx-syy
        rloc = diff*diff+4.*sxy*sxy
        cloc1 = 'rloc'
        cloc2 = 'sigmapr3ETCtg'
        call check_sqrt(rloc,cloc1,cloc2,lerror)
c!       **** **********
        if (lerror) return
        S1 = (SQRT(rloc) + sxx+syy) / r2
        S2 = (-SQRT(rloc) + sxx+syy) / r2
        S3 = 0.d0
        if (S1.lt.S2) then
          rloc = S1
          S1 = S2
          S2 = rloc
        endif
 
        lfirststress = .true.
        Smax = S1
    1   continue
        if (.not.lfirststress) Smax = S2
 
c       Main direction of Smax
c       -----------------
c       condition for the case Smax = Sxx = Syy (then deltax=deltay=deltaz=0)
        if ( (ABS(Sxx-Smax).le.r1).and.(ABS(Syy-Smax).le.r1).and.
     .       (ABS(Szz-Smax).gt.r1) ) then
          if(lfirststress) then
            rnx = (SQRT(r2))/r2
          else
            rnx = - (SQRT(r2))/r2
          endif
          rny = (SQRT(r2))/r2
          rnz = r0
c       condition for the case Smax = Syy = Szz (then deltax=deltay=deltaz=0)
        else if ( (ABS(Syy-Smax).le.r1).and.(ABS(Szz-Smax).le.r1)
     .       .and.(ABS(Sxx-Smax).gt.r1) ) then
          if(lfirststress) then
            rny = (SQRT(r2))/r2
          else
            rny = - (SQRT(r2))/r2
          endif
          rnz = (SQRT(r2))/r2
          rnx = r0
c       condition for the case Smax = Sxx = Szz (then deltax=deltay=deltaz=0)
        else if ( (ABS(Sxx-Smax).le.r1).and.(ABS(Szz-Smax).le.r1)
     .       .and.(ABS(Syy-Smax).gt.r1) ) then
          if(lfirststress) then
            rnx = (SQRT(r2))/r2
          else
            rnx = - (SQRT(r2))/r2
          endif
          rnz = (SQRT(r2))/r2
          rny = r0
c       condition for the case Smax = Sxx = Syy = Szz (then deltax=deltay=deltaz=0)
        else if ( (ABS(Sxx-Smax).le.r1).and.(ABS(Syy-Smax).le.r1)
     .       .and.(ABS(Szz-Smax).le.r1).and.(ABS(Sxy).le.r1).and.
     .        (ABS(Sxz).le.r1).and.(ABS(Syz).le.r1) ) then
          rnx = (SQRT(r3))/r3
          rny = (SQRT(r3))/r3
          rnz = (SQRT(r3))/r3
        else
          deltax = (Szz-Smax)*(Syy-Smax)-Syz*Syz
          deltay = (Sxx-Smax)*(Szz-Smax)-Sxz*Sxz
          deltaz = (Sxx-Smax)*(Syy-Smax)-Sxy*Sxy
          deltamax = ABS(deltax)
          ldeltay = .false.
          ldeltaz = .false.
          if (ABS(deltay).gt.deltamax) then
            deltamax = ABS(deltay)
            ldeltay = .true.
          endif
          if (ABS(deltaz).gt.deltamax) then
            deltamax = ABS(deltaz)
            ldeltay = .false.
            ldeltaz = .true.
          endif
 
          if (ldeltaz) then
c           Main orientation OZ
            aa = ( Sxy*Syz - Sxz*(Syy-Smax) ) / deltaz
            bb = ( Sxz*Sxy - Syz*(Sxx-Smax) ) / deltaz
            rloc = SQRT( aa*aa + bb*bb + r1 )
            rnz = r1/rloc
            rnx = aa/rloc
            rny = bb/rloc
          else if (ldeltay) then
c           Main orientation OY
            aa = ( Sxz*Syz - Sxy*(Szz-Smax) ) / deltay
            bb = ( Sxz*Sxy - Syz*(Sxx-Smax) ) / deltay
            rloc = SQRT( aa*aa + bb*bb + r1 )
            rny = r1/rloc
            rnx = aa/rloc
            rnz = bb/rloc
          else
c           Main orientation OX
            aa = ( Sxz*Syz - Sxy*(Szz-Smax) ) / deltax
            bb = ( Syz*Sxy - Sxz*(Syy-Smax) ) / deltax
            rloc = SQRT( aa*aa + bb*bb + r1 )
            rnx = r1/rloc
            rny = aa/rloc
            rnz = bb/rloc
          endif
        endif
 
c       The 3 following lines to try to solve the problem of unsymetry
c       in routine STRESS
        if (ABS(rnx).lt.precision3) rnx = r0
        if (ABS(rny).lt.precision3) rny = r0
        if (ABS(rnz).lt.precision3) rnz = r0
 
        if(lfirststress) then
c         go and calculate for S2
c         -----------------------
          lfirststress = .false.
          rcossigmapr(i1,i1) = rnx
          rcossigmapr(i2,i1) = rny
          rcossigmapr(i3,i1) = rnz
          go to 1
        else
c         S2 has been calculated
c         ----------------------
          rcossigmapr(i1,i2) = rnx
          rcossigmapr(i2,i2) = rny
          rcossigmapr(i3,i2) = rnz
        endif
c       S1 and S2 have been calculated
c       ------------------------------
 
c       calculate S3.
c       -------------
        rcossigmapr(i1,i3) = rcossigmapr(i2,i1)*rcossigmapr(i3,i2)
     .                   - rcossigmapr(i3,i1)*rcossigmapr(i2,i2)
        rcossigmapr(i2,i3) = rcossigmapr(i3,i1)*rcossigmapr(i1,i2)
     .                   - rcossigmapr(i1,i1)*rcossigmapr(i3,i2)
        rcossigmapr(i3,i3) = rcossigmapr(i1,i1)*rcossigmapr(i2,i2)
     .                   - rcossigmapr(i2,i1)*rcossigmapr(i1,i2)
 
      else
c       FULL 3D
c       Hydrostatic stress
c       ------------------
        rI1 = Sxx+Syy+Szz
        Sh = rI1/r3
        SumS = rI1+Sxy+Sxz+Syz
 
c       2nd invariant of the deviator J
c       -------------------------------
        rJ2 = (   (Sxx-Syy) * (Sxx-Syy)
     .          + (Sxx-Szz) * (Sxx-Szz)
     .          + (Syy-Szz) * (Syy-Szz)  ) / r6
     .        + Sxy*Sxy + Sxz*Sxz + Syz*Syz
 
c       octahedral shear stress
c       -----------------------
        toc = r2*rJ2/r3
        cloc1 = 'toc'
        cloc2 = 'Sigmapr3'
        toc = SQRT(toc)
 
        smallS = SumS*precision2
 
        if (ABS(toc).gt.ABS(smallS)) then
 
c         3th invariant of the deviator J
c         -------------------------------
          rI2 =   Sxx*Syy + Sxx*Szz + Syy*Szz
     .          - Sxy*Sxy - Sxz*Sxz - Syz*Syz
          rI3 =   Sxx*Syy*Szz + r2*Sxy*Syz*Sxz
     .          - Sxx*Syz*Syz - Syy*Sxz*Sxz - Szz*Sxy*Sxy
          rJ3 = rI3 - rI1*rI2/r3 + r2*rI1*rI1*rI1/27d0
 
c         angle omega
c         -----------
          rloc = ( racine2*rJ3 / (toc*toc*toc) )
          if (rloc.ge.1d0) then
            omega = r0
          else if (rloc.le.(-1d0)) then
            omega = ( ACOS( (-1d0) ) )/r3
          else
            omega = ( ACOS( racine2*rJ3 / (toc*toc*toc) ) )/r3
            if ((omega.lt.r0).or.(omega.gt.pi/r3)) then
              print *,'Subroutine sigmapr3 - ERROR in omega'
              print *,'omega = ',omega
              stop
            endif
          endif
 
c         principal stresses
c         ------------------
          S1 = racine2*toc*COS(omega-r2*pi/r3)+Sh
          S2 = racine2*toc*COS(omega+r2*pi/r3)+Sh
          S3 = racine2*toc*COS(omega)         +Sh
 
c         sort the principal stresses  : S1 > S2 > S3
          if (S2.lt.S3) then
            rloc = S2
            S2 = S3
            S3 = rloc
          endif
          if (S1.lt.S2) then
            rloc = S1
            S1 = S2
            S2 = rloc
          endif
          if (S2.lt.S3) then
            rloc = S2
            S2 = S3
            S3 = rloc
          endif
 
          lfirststress = .true.
          Smax = S1
    2     continue
          if (.not.lfirststress) Smax = S2
 
c         Main direction of Smax
c         -----------------
c         condition for the case Smax = Sxx = Syy (then deltax=deltay=deltaz=0)
          if ( (ABS(Sxx-Smax).le.r1).and.(ABS(Syy-Smax).le.r1).and.
     .         (ABS(Szz-Smax).gt.r1) ) then
            if(lfirststress) then
              rnx = (SQRT(r2))/r2
            else
              rnx = - (SQRT(r2))/r2
            endif
            rny = (SQRT(r2))/r2
            rnz = r0
c         condition for the case Smax = Syy = Szz (then deltax=deltay=deltaz=0)
          else if ( (ABS(Syy-Smax).le.r1).and.(ABS(Szz-Smax).le.r1)
     .         .and.(ABS(Sxx-Smax).gt.r1) ) then
            if(lfirststress) then
              rny = (SQRT(r2))/r2
            else
              rny = - (SQRT(r2))/r2
            endif
            rnz = (SQRT(r2))/r2
            rnx = r0
c         condition for the case Smax = Sxx = Szz (then deltax=deltay=deltaz=0)
          else if ( (ABS(Sxx-Smax).le.r1).and.(ABS(Szz-Smax).le.r1)
     .         .and.(ABS(Syy-Smax).gt.r1) ) then
            if(lfirststress) then
              rnx = (SQRT(r2))/r2
            else
              rnx = - (SQRT(r2))/r2
            endif
            rnz = (SQRT(r2))/r2
            rny = r0
c         condition for the case Smax = Sxx = Syy = Szz (then deltax=deltay=deltaz=0)
          else if ( (ABS(Sxx-Smax).le.r1).and.(ABS(Syy-Smax).le.r1)
     .         .and.(ABS(Szz-Smax).le.r1).and.(ABS(Sxy).le.r1).and.
     .          (ABS(Sxz).le.r1).and.(ABS(Syz).le.r1) ) then
            rnx = (SQRT(r3))/r3
            rny = (SQRT(r3))/r3
            rnz = (SQRT(r3))/r3
          else
            deltax = (Szz-Smax)*(Syy-Smax)-Syz*Syz
            deltay = (Sxx-Smax)*(Szz-Smax)-Sxz*Sxz
            deltaz = (Sxx-Smax)*(Syy-Smax)-Sxy*Sxy
            deltamax = ABS(deltax)
            ldeltay = .false.
            ldeltaz = .false.
            if (ABS(deltay).gt.deltamax) then
              deltamax = ABS(deltay)
              ldeltay = .true.
            endif
            if (ABS(deltaz).gt.deltamax) then
              deltamax = ABS(deltaz)
              ldeltay = .false.
              ldeltaz = .true.
            endif
 
            if (ldeltaz) then
c             Main orientation OZ
              aa = ( Sxy*Syz - Sxz*(Syy-Smax) ) / deltaz
              bb = ( Sxz*Sxy - Syz*(Sxx-Smax) ) / deltaz
              rloc = SQRT( aa*aa + bb*bb + r1 )
              rnz = r1/rloc
              rnx = aa/rloc
              rny = bb/rloc
            else if (ldeltay) then
c             Main orientation OY
              aa = ( Sxz*Syz - Sxy*(Szz-Smax) ) / deltay
              bb = ( Sxz*Sxy - Syz*(Sxx-Smax) ) / deltay
              rloc = SQRT( aa*aa + bb*bb + r1 )
              rny = r1/rloc
              rnx = aa/rloc
              rnz = bb/rloc
            else
c             Main orientation OX
              aa = ( Sxz*Syz - Sxy*(Szz-Smax) ) / deltax
              bb = ( Syz*Sxy - Sxz*(Syy-Smax) ) / deltax
              rloc = SQRT( aa*aa + bb*bb + r1 )
              rnx = r1/rloc
              rny = aa/rloc
              rnz = bb/rloc
            endif
          endif
 
c         The 3 following lines to try to solve the problem of unsymetry
c         in routine STRESS
          if (ABS(rnx).lt.precision3) rnx = r0
          if (ABS(rny).lt.precision3) rny = r0
          if (ABS(rnz).lt.precision3) rnz = r0
 
          if(lfirststress) then
c           go and calculate for S2
c           -----------------------
            lfirststress = .false.
            rcossigmapr(i1,i1) = rnx
            rcossigmapr(i2,i1) = rny
            rcossigmapr(i3,i1) = rnz
            go to 2
          else
c           S2 has been calculated
c           ----------------------
            rcossigmapr(i1,i2) = rnx
            rcossigmapr(i2,i2) = rny
            rcossigmapr(i3,i2) = rnz
          endif
c         S1 and S2 have been calculated
c         ------------------------------
 
c         calculate S3.
c         -------------
          rcossigmapr(i1,i3) = rcossigmapr(i2,i1)*rcossigmapr(i3,i2)
     .                     - rcossigmapr(i3,i1)*rcossigmapr(i2,i2)
          rcossigmapr(i2,i3) = rcossigmapr(i3,i1)*rcossigmapr(i1,i2)
     .                     - rcossigmapr(i1,i1)*rcossigmapr(i3,i2)
          rcossigmapr(i3,i3) = rcossigmapr(i1,i1)*rcossigmapr(i2,i2)
     .                    - rcossigmapr(i2,i1)*rcossigmapr(i1,i2)
 
        else
          S1 = sigma(i1)
          S2 = sigma(i2)
          S3 = sigma(i3)
c         case where S1 = S2 = S3
          rcossigmapr(1,1)=1.0d0
          rcossigmapr(2,2)=1.0d0
          rcossigmapr(3,3)=1.0d0
 
        endif
 
      endif
 
      return
      end