pri324.F File Reference

#include "implicit_f.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	pri324 (sig, epstot, eps, vec)

Function/Subroutine Documentation

pri324()

subroutine pri324	(		sig,
			epstot,
			eps,
			vec )

Definition at line 28 of file pri324.F.

C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      my_real
     .   sig(*), epstot(*), eps(*), vec(*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER IPERM(3), I, L, LMAX
      my_real
     .   cs(6), str(3), a(3,3), v(3,3), b(3,3), xmag(3), pr, aa, bb,
     .   cc, angp, dd, ftpi, ttpi, strmax, tol1, tol2, xmax, vmag, s11,
     .   s21, s31, s12, s22, s32, s13, s23, s33, a11, a12, a13, a21,
     .   a22, a23, a31, a32, a33
C-----------------------------------------------
      DATA ftpi,ttpi / 4.188790205, 2.094395102 /
      DATA iperm/2,3,1/
C=======================================================================
C       DEVIATEUR PRINCIPAL DE CONTRAINTE
C . . . . . . . . . . . . . . . . . . .
      DO i=1,6
        vec(i)=zero
        cs(i)=sig(i)
      ENDDO
      vec(1)=1
      vec(5)=1
C
      pr = -(cs(1)+cs(2)+cs(3)) * third                               
      cs(1)=cs(1) + pr                                                
      cs(2)=cs(2) + pr                                                
      cs(3)=cs(3) + pr                                                
C
      aa = cs(4)**2 + cs(5)**2 + cs(6)**2 - cs(1)*cs(2) - cs(2)*cs(3) 
     &   - cs(1)*cs(3)                                                
 
      IF (aa < em20) RETURN                                           
C
      bb = cs(1)*cs(5)**2 + cs(2)*cs(6)**2 + cs(3)*cs(4)**2             
     &   - cs(1)*cs(2)*cs(3) - two*cs(4)*cs(5)*cs(6)                  
C
      cc=-sqrt(twenty7/aa)*bb*half/aa                                
      cc= min(cc,one)                                                  
      cc= max(cc,-one)                                                 
      angp=acos(cc) * third                                           
      dd=two*sqrt(aa*third)                                          
      str(1)=dd*cos(angp)                                             
      str(2)=dd*cos(angp+ftpi)                                        
      str(3)=dd*cos(angp+ttpi)                                        
C . . . . . . . . . . .
C      VECTEURS PROPRES
C . . . . . . . . . . .
        strmax= max(abs(str(1)),abs(str(3)))
 
        tol1= max(em20,sixem4*strmax**2)
        tol2= twoem4*strmax
 
        a(1,1)=cs(1)-str(1)
        a(2,2)=cs(2)-str(1)
        a(3,3)=cs(3)-str(1)
        a(1,2)=cs(4)
        a(2,1)=cs(4)
        a(2,3)=cs(5)
        a(3,2)=cs(5)
        a(1,3)=cs(6)
        a(3,1)=cs(6)
C
        DO l=1,3
          b(1,l) = a(2,l)*a(3,iperm(l))-a(3,l)*a(2,iperm(l))
          b(2,l) = a(3,l)*a(1,iperm(l))-a(1,l)*a(3,iperm(l))
          b(3,l) = a(1,l)*a(2,iperm(l))-a(2,l)*a(1,iperm(l))
          xmag(l)= sqrt(b(1,l)**2+b(2,l)**2+b(3,l)**2)
        ENDDO
 
        xmax=zero
        lmax = 1
        DO l=1,3
          IF (xmag(l) > xmax) THEN
            xmax=xmag(l)
            lmax=l
          ENDIF
        ENDDO
C
        IF (xmax > tol1) THEN
           v(1,1)=b(1,lmax)/xmax
           v(2,1)=b(2,lmax)/xmax
           v(3,1)=b(3,lmax)/xmax
           a(1,1)=cs(1)-str(3)
           a(2,2)=cs(2)-str(3)
           a(3,3)=cs(3)-str(3)
           a(1,2)=cs(4)
           a(2,1)=cs(4)
           a(2,3)=cs(5)
           a(3,2)=cs(5)
           a(1,3)=cs(6)
           a(3,1)=cs(6)
C
           DO l=1,3
             b(1,l)=a(2,l)*v(3,1)-a(3,l)*v(2,1)
             b(2,l)=a(3,l)*v(1,1)-a(1,l)*v(3,1)
             b(3,l)=a(1,l)*v(2,1)-a(2,l)*v(1,1)
             xmag(l)=sqrt(b(1,l)**2+b(2,l)**2+b(3,l)**2)
           ENDDO
           xmax=zero
           DO l=1,3
             IF(xmag(l) > xmax)THEN
               xmax=xmag(l)
               lmax=l
             ENDIF
           ENDDO
C
           IF (xmax > tol2) THEN
             v(1,3)= b(1,lmax)/xmax
             v(2,3)= b(2,lmax)/xmax
             v(3,3)= b(3,lmax)/xmax
             v(1,2)= v(2,3)*v(3,1)-v(2,1)*v(3,3)
             v(2,2)= v(3,3)*v(1,1)-v(3,1)*v(1,3)
             v(3,2)= v(1,3)*v(2,1)-v(1,1)*v(2,3)
             vmag  = sqrt(v(1,2)**2 + v(2,2)**2 + v(3,2)**2)
             v(1,2)= v(1,2)/vmag
             v(2,2)= v(2,2)/vmag
             v(3,2)= v(3,2)/vmag
           ELSE
             vmag  = sqrt(v(1,1)**2 + v(2,1)**2 + v(3,1)**2)
             IF (vmag > tol2/strmax) THEN
               v(1,2)=-v(2,1)/vmag
               v(2,2)= v(1,1)/vmag
               v(3,2)= zero
             ELSE
               v(1,2)=one
               v(2,2)=zero
               v(3,2)=zero     
             ENDIF
           ENDIF
        ELSE
C . . . . . . . . . . . . .
C         SOLUTION DOUBLE
C . . . . . . . . . . . . .
          DO l=1,3
            xmag(l) = sqrt(a(1,l)**2 + a(2,l)**2)
          ENDDO
          xmax = zero              
          DO l=1,3             
            IF(xmag(l) > xmax)THEN 
              lmax=l               
              xmax=xmag(l)         
            ENDIF                  
          ENDDO
C
         IF(max(abs(a(3,1)),abs(a(3,2)),abs(a(3,3))) < tol2)THEN
           v(1,1) = zero
           v(2,1) = zero
           v(3,1) = one   
           v(1,2) = -a(2,lmax)/xmax
           v(2,2) =  a(1,lmax)/xmax
           v(3,2) =  zero
C
         ELSEIF(xmax > tol2)THEN
           v(1,1) = -a(2,lmax)/xmax
           v(2,1) =  a(1,lmax)/xmax
           v(3,1) =  zero
           v(1,2) = -a(3,lmax)*v(2,1)
           v(2,2) =  a(3,lmax)*v(1,1)
           v(3,2) =  a(1,lmax)*v(2,1)-a(2,lmax)*v(1,1)
           vmag = sqrt(v(1,2)**2 + v(2,2)**2 + v(3,2)**2)
           v(1,2)=v(1,2)/vmag
           v(2,2)=v(2,2)/vmag
           v(3,2)=v(3,2)/vmag
         ELSE
           v(1,1) = one
           v(2,1) = zero
           v(3,1) = zero
           v(1,2) = zero
           v(2,2) = one
           v(3,2) = zero     
         ENDIF
      ENDIF
C
      vec(1)=v(1,1)
      vec(2)=v(2,1)
      vec(3)=v(3,1)
      vec(4)=v(1,2)
      vec(5)=v(2,2)
      vec(6)=v(3,2)
C . . . . . . . . . . . .
C     ROTATION EPS EPSTOT
C . . . . . . . . . . . .
      s11=vec(1)
      s21=vec(2)
      s31=vec(3)
      s12=vec(4)
      s22=vec(5)
      s32=vec(6)
      s13=s21*s32-s31*s22
      s23=s31*s12-s11*s32
  5   s33=s11*s22-s21*s12
C
      a11=sig(1)*s11+sig(4)*s21+sig(6)*s31
      a12=sig(1)*s12+sig(4)*s22+sig(6)*s32
      a13=sig(1)*s13+sig(4)*s23+sig(6)*s33
      a21=sig(4)*s11+sig(2)*s21+sig(5)*s31
      a22=sig(4)*s12+sig(2)*s22+sig(5)*s32
      a23=sig(4)*s13+sig(2)*s23+sig(5)*s33
      a31=sig(6)*s11+sig(5)*s21+sig(3)*s31
      a32=sig(6)*s12+sig(5)*s22+sig(3)*s32
      a33=sig(6)*s13+sig(5)*s23+sig(3)*s33
      sig(4)=s11*a12+s21*a22+s31*a32
      sig(5)=s12*a13+s22*a23+s32*a33
      sig(6)=s11*a13+s21*a23+s31*a33
C
      eps(4)=half*eps(4)
      eps(5)=half*eps(5)
      eps(6)=half*eps(6)      
      a11=eps(1)*s11+eps(4)*s21+eps(6)*s31
      a12=eps(1)*s12+eps(4)*s22+eps(6)*s32
      a13=eps(1)*s13+eps(4)*s23+eps(6)*s33
      a21=eps(4)*s11+eps(2)*s21+eps(5)*s31
      a22=eps(4)*s12+eps(2)*s22+eps(5)*s32
      a23=eps(4)*s13+eps(2)*s23+eps(5)*s33
      a31=eps(6)*s11+eps(5)*s21+eps(3)*s31
      a32=eps(6)*s12+eps(5)*s22+eps(3)*s32
      a33=eps(6)*s13+eps(5)*s23+eps(3)*s33
      eps(1)=s11*a11+s21*a21+s31*a31
      eps(2)=s12*a12+s22*a22+s32*a32
      eps(3)=s13*a13+s23*a23+s33*a33
      epstot(4)=half*epstot(4)
      epstot(5)=half*epstot(5)
      epstot(6)=half*epstot(6)      
      a11=epstot(1)*s11+epstot(4)*s21+epstot(6)*s31
      a12=epstot(1)*s12+epstot(4)*s22+epstot(6)*s32
      a13=epstot(1)*s13+epstot(4)*s23+epstot(6)*s33
      a21=epstot(4)*s11+epstot(2)*s21+epstot(5)*s31
      a22=epstot(4)*s12+epstot(2)*s22+epstot(5)*s32
      a23=epstot(4)*s13+epstot(2)*s23+epstot(5)*s33
      a31=epstot(6)*s11+epstot(5)*s21+epstot(3)*s31
      a32=epstot(6)*s12+epstot(5)*s22+epstot(3)*s32
      a33=epstot(6)*s13+epstot(5)*s23+epstot(3)*s33
      epstot(1)=s11*a11+s21*a21+s31*a31
      epstot(2)=s12*a12+s22*a22+s32*a32
      epstot(3)=s13*a13+s23*a23+s33*a33
c-----------
      RETURN