PROGRAM PRINCIPAL_STRESS
C=======================================================================
C COMPUTATION OF PRINCIPAL STESSES USING various METHODS
C EIJI FUKUMORI
C 30 DECEMBER 2012
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
INCLUDE 'PARAM.DAT'
CCCCC PARAMETER ( MXE=40000,MXN=53000,MXB=5000, ND=8 )
C=======================================================================
C ARRAYS
CHARACTER INPFILE*14, STRESSF*12
DIMENSION GLBTAU(7,MXN)
C=======================================================================
EPS = 1.D-14
C=======================================================================
INPFILE = 'NONE'
IF ( ND .EQ. 8 ) THEN
INPFILE = 'STATIC3D08.DAT'
STRESSF = 'STRESS08.OUT'
END IF
IF ( ND .EQ. 27 ) THEN
INPFILE = 'STATIC3D27.DAT'
STRESSF = 'STRESS27.OUT'
END IF
IF ( INPFILE .EQ. 'NONE' ) STOP 'NO INPUT FILE AT MAIN'
C=======================================================================
WRITE(*,*) 'PRINCIPAL STESSES COMPUTATION PROGRAM'
WRITE (*,*)' READING IN DATA FILES'
CALL INPUT ( INPFILE,STRESSF,MXN,NNODE,GLBTAU )
C=======================================================================
WRITE (*,*) 'FILEMK '
CALL FILEMK ( EPS,ND,MXN, NNODE, GLBTAU )
STOP' NORMAL TERMINATION'
END
C
C
SUBROUTINE FILEMK ( EPS,ND,MXN,NNODE,GLBTAU )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION GLBTAU(7,MXN)
DIMENSION TAUJ(3), VECJ(3,3), TAUP(3), VECP(3,3)
DIMENSION TAUN(3), VECN(3,3), TAUC(3), VECC(3,3)
DIMENSION TAUB(3), VECB(3,3)
CHARACTER OUTFILE*15
LOGICAL YES
C========> GLBTAU(1,I)=TAUXX = MU*(2DU/DX)+LAMBDA*DIVV
C========> GLBTAU(2,I)=TAUYY = MU*(2DV/DY)+LAMBDA*DIVV
C========> GLBTAU(3,I)=TAUZZ = MU*(2DW/DZ)+LAMBDA*DIVV
C========> GLBTAU(4,I)=TAUXY = MU*(DU/DY+DU/DX)
C========> GLBTAU(5,I)=TAUXZ = MU*(DU/DZ+DW/DX)
C========> GLBTAU(6,I)=TAUYZ = MU*(DV/DZ+DW/DY)
C========> GLBTAU(7,I)= DIVV
C--------0---------0---------0---------0---------0---------0---------0--
C
OUTFILE = 'NONE'
IF ( ND .EQ. 8 ) THEN
OUTFILE = 'PRINCIPAL08.OUT'
END IF
IF ( ND .EQ. 27 ) THEN
OUTFILE = 'PRINCIPAL27.OUT'
END IF
IF ( OUTFILE .EQ. 'NONE' ) STOP 'NO FILE ASIGNED'
C
OPEN ( 1, FILE=OUTFILE, STATUS='UNKNOWN' )
WRITE (1,*) '# TXX TYY TZZ TXY TXZ TYZ ',
* ' ITJ T1J T2J T3J ',
* ' A1XJ A1YJ A1ZJ A2XJ A2YJ A2ZJ A3XJ A3YJ A3ZJ ',
* ' ITP T1P T2P T3P ',
* ' A1XP A1YP A1ZP A2XP A2YP A2ZP A3XP A3YP A3ZP ',
* ' ITN T1N T2N T3N ',
* ' A1XN A1YN A1ZN A2XN A2YN A2ZN A3XN A3YN A3ZN ',
* ' ITC T1C T2C T3C ',
* ' A1XC A1YC A1ZC A2XC A2YC A2ZC A3XC A3YC A3ZC ',
* ' ITB T1B T2B T3B ',
* ' A1XB A1YB A1ZB A2XB A2YB A2ZB A3XB A3YB A3ZB '
C
DO I = 1 , NNODE
C----- PRINCIPAL STRESSES TAU'S ----------
CALL JACOBI (EPS,GLBTAU(1,I),GLBTAU(2,I),GLBTAU(3,I),
* GLBTAU(4,I),GLBTAU(5,I),GLBTAU(6,I), ICJ, TAUJ, VECJ )
C
CALL POWER (EPS,GLBTAU(1,I),GLBTAU(2,I),GLBTAU(3,I),
* GLBTAU(4,I),GLBTAU(5,I),GLBTAU(6,I), ICP, TAUP, VECP )
C
CALL NEWTON (EPS,GLBTAU(1,I),GLBTAU(2,I),GLBTAU(3,I),
* GLBTAU(4,I),GLBTAU(5,I),GLBTAU(6,I), ICN, TAUN, VECN )
C
CALL CARDANO (EPS,GLBTAU(1,I),GLBTAU(2,I),GLBTAU(3,I),
* GLBTAU(4,I),GLBTAU(5,I),GLBTAU(6,I), ICC, TAUC, VECC )
C
CALL BSECTION (EPS,GLBTAU(1,I),GLBTAU(2,I),GLBTAU(3,I),
* GLBTAU(4,I),GLBTAU(5,I),GLBTAU(6,I), ICB, TAUB, VECB )
C
C----- FILE--------VECX(J,MODE)
WRITE (1,*)
* I, (GLBTAU(J,I),J=1,6),
* ICJ, (TAUJ(J),J=1,3),((VECJ(J,MODE),J=1,3),MODE=1,3),
* ICP, (TAUP(J),J=1,3),((VECP(J,MODE),J=1,3),MODE=1,3),
* ICN, (TAUN(J),J=1,3),((VECN(J,MODE),J=1,3),MODE=1,3),
* ICC, (TAUC(J),J=1,3),((VECC(J,MODE),J=1,3),MODE=1,3),
* ICB, (TAUB(J),J=1,3),((VECB(J,MODE),J=1,3),MODE=1,3)
END DO
CLOSE (1)
RETURN
END
C
C
SUBROUTINE INPUT ( INPFILE,STRESSF,MXN,NNODE,GLBTAU )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION GLBTAU(7,MXN)
CHARACTER*14 INPFILE*14, STRESSF*12
LOGICAL YES
C========> FILENAME INPFILE
OPEN ( 1, FILE=INPFILE, STATUS='UNKNOWN' )
C========> PARAMETERS
READ (1,*) YOUNG, POISSON
C========> ELEMENTS
READ (1,*) NE
DO I = 1 , NE
READ (1,*)
END DO
CC========>
READ (1,*) NNODE
CLOSE (1)
C========> FILENAME STRESSF
OPEN ( 1, FILE=STRESSF, STATUS='UNKNOWN' )
READ (1,*)
DO I = 1 , NNODE
READ (1,*) NODE,(GLBTAU(J,I),J=1,7)
END DO
CLOSE (1)
RETURN
END
C
C
SUBROUTINE JACOBI (EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, IC, TAU, U )
C=======================================================================
C PRINCIPAL STRESSES USING JACOB METHOD
C=======================================================================
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER ( MAXTRTN=1000 )
DIMENSION A(3,3),U(3,3), TAU(3), UTEMP(3)
C---------------------- RESETTING EIGEN VECTORS ------------------------
IC = 0
NNODE = 3
DO I=1,NNODE
DO J=1,NNODE
U(I,J) = 0.0D0
END DO
U(I,I) = 1.0D0
END DO
C------------ SUBSTITUTION OF STRESSES INTO [A] MATRIX -----------------
A(1,1) = TXX
A(2,2) = TYY
A(3,3) = TZZ
A(1,2) = TXY
A(1,3) = TXZ
A(2,3) = TYZ
A(2,1) = A(1,2)
A(3,1) = A(1,3)
A(3,2) = A(2,3)
C
ITERATION = 0
AMXOFF = EPS
DMX = EPS
C=======================================================================
DO WHILE ( (AMXOFF/DMX.GT.EPS) .AND. (ITERATION.LT.MAXTRTN) )
DO I = 1 , NNODE
IF ( DABS(A(I,I)) .GT. DMX ) DMX = DABS(A(I,I))
END DO
C--------- FIND MAXIMUM IN OFF DIAGONAL ELEMENTS OF A(I,J) -------------
AMXOFF = 0.D0
DO I = 1 , NNODE-1
DO J = I+1 , NNODE
IF ( DABS(A(I,J)) .GT. AMXOFF ) THEN
AMXOFF = DABS(A(I,J))
K = I
L = J
END IF
END DO
END DO
IF ( AMXOFF .EQ. 0.D0 ) EXIT
C------------------ COMPUTATION OF ROTATIONAL ANGLE --------------------
C------------------ R IMPLIES RADIUS OF MOHR CIRCLE --------------------
C--------- CENTER IMPLIES CENTER COORDINATE OF MOHR CIRCLE -------------
ITERATION = ITERATION + 1
DIFFHALF = (A(K,K) - A(L,L))/2.D0
IF ( DIFFHALF .EQ. 0.D0 ) THEN
R = DABS(A(K,L))
COSINE = 1.D0/DSQRT(2.D0)
SINE = 1.D0/DSQRT(2.D0)
IF ( A(K,L) .LT. 0.D0 ) SINE = -SINE
IF ( R .EQ. 0.D0 ) THEN
COSINE = 1.D0
SINE = 0.D0
END IF
ELSE
R = DSQRT ( DIFFHALF**2 + A(K,L)**2 )
COSINE = DSQRT ( DABS(DIFFHALF)/(2.D0*R) + 0.5D0 )
SINE = A(K,L)/(2.D0*R*COSINE)
IF ( DIFFHALF .LT. 0.D0 ) COSINE = -COSINE
END IF
C------------------------- ORTHOGNALIZATION ----------------------------
IF ( R .GT. 0.D0 ) THEN
DO I = 1 , NNODE
AW = U(I,K)
AZ = U(I,L)
U(I,K) = AW*COSINE + AZ*SINE
U(I,L) = -AW*SINE + AZ*COSINE
IF ( ( I .NE. K ) .AND. ( I .NE. L ) ) THEN
AW = A(I,K)
AZ = A(I,L)
A(I,K) = AW*COSINE + AZ*SINE
A(I,L) = -AW*SINE + AZ*COSINE
A(K,I) = A(I,K)
A(L,I) = A(I,L)
END IF
END DO
A(K,L) = 0.D0
A(L,K) = 0.D0
CENTER = 0.5D0*(A(K,K)+A(L,L))
IF ( DIFFHALF .LT. 0.D0 ) R = -R
A(K,K) = CENTER + R
A(L,L) = CENTER - R
C----------------- LINES BELOW ALSO WORKS FINE.------------------------
C AKK = A(K,K)
C ALL = A(L,L)
C CROSS = 2.D0*A(K,L)*SINE*COSINE
C A(K,K) = AKK*COSINE**2 + ALL*SINE**2 + CROSS
C A(L,L) = AKK*SINE**2 + ALL*COSINE**2 - CROSS
C----------------------------------------------------------------------
END IF
END DO
C------------------------- EIGENVALUES = A(I,I)-------------------------
IC = ITERATION
TAU(1) = A(1,1)
TAU(2) = A(2,2)
TAU(3) = A(3,3)
C-------- MAKING TAU(1), TAU(2), AND TAU(3) DESCENDING ORDER -----------
IF ( TAU(2) .GT. TAU(1) ) THEN
TEMP = TAU(1)
TAU(1) = TAU(2)
TAU(2) = TEMP
DO J = 1 , 3
UTEMP(J) = U (J,1)
U (J,1) = U (J,2)
U (J,2) = UTEMP(J)
END DO
END IF
IF ( TAU(3) .GT. TAU(2) ) THEN
TEMP = TAU(2)
TAU(2) = TAU(3)
TAU(3) = TEMP
DO J = 1 , 3
UTEMP(J) = U (J,2)
U (J,2) = U (J,3)
U (J,3) = UTEMP(J)
END DO
END IF
IF ( TAU(2) .GT. TAU(1) ) THEN
TEMP = TAU(1)
TAU(1) = TAU(2)
TAU(2) = TEMP
DO J = 1 , 3
UTEMP(J) = U (J,1)
U (J,1) = U (J,2)
U (J,2) = UTEMP(J)
END DO
END IF
RETURN
END
C
C
SUBROUTINE NEWTON (EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, IC, TAU,U )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION TAU(3),DRV(3,3),DI(3,3),F(3),H(3,3),G(3,3),U(3,3)
DIMENSION TEMPU(3)
DATA MAXITA / 100 /
C-------------------------- STRESS INVARIANTS --------------------------
A = TXX+TYY+TZZ
B = TXX*TYY+TYY*TZZ+TZZ*TXX-TXY**2-TYZ**2-TXZ**2
C = TXX*TYY*TZZ+2.D0*TXY*TXZ*TYZ-
* TXY**2*TZZ-TYZ**2*TXX-TXZ**2*TYY
C---------------- ENERGY BY COMPUTAED GIVEN STRESSES -------------------
ENERGY0 = (TXX-TYY)**2+(TYY-TZZ)**2+(TZZ-TXX)**2+
* 6.D0*(TXY**2+TYZ**2+TXZ**2)
C------------------------ RESET EIGEN VECTORS --------------------------
DO I = 1 , 3
DO J = 1 , 3
U(I,J) = 0.D0
END DO
U(I,I) = 1.D0
END DO
C-------------------------- CASE OF A=B=C=0 ----------------------------
IF ( A.EQ.0.D0 .AND. B.EQ.0.D0 .AND. C.EQ.0.D0 ) THEN
IC = 0
TAU(1) = 0.D0
TAU(2) = 0.D0
TAU(3) = 0.D0
RETURN
END IF
C--------------------------- CASE OF B=C=0 -----------------------------
IC = -1
IF ( B .EQ. 0.D0 .AND. C .EQ. 0.D0 ) THEN
TAU(2) = 0.D0
IF ( A .GT. 0.D0 ) THEN
TAU(1) = A
TAU(3) = 0.D0
ELSE
TAU(1) = 0.D0
TAU(3) = A
END IF
RETURN
END IF
C---------------------------- CASE OF C=0 ------------------------------
IC = -2
IF ( C .EQ. 0.D0 ) THEN
TAU(1) = (A+DSQRT(A*A-4.D0*B))/2.D0
TAU(2) = 0.D0
TAU(3) = (A-DSQRT(A*A-4.D0*B))/2.D0
IF ( TAU(3) .GT. 0.D0 ) THEN
TAU(2) = TAU(3)
TAU(3) = 0.D0
END IF
END IF
C--------------- CASE OF P=Q=0 MEANING TRIPLE ROOTS---------------------
IC = -3
P = -(-A)**2/3.D0 + B
Q = 2.D0*(-A)**3/27.D0 -(-A)*B/3.D0 + (-C)
IF ( P.EQ.0.D0 .AND. Q.EQ.0.D0 ) THEN
TAU(1) = TXX
TAU(2) = TYY
TAU(3) = TZZ
RETURN
END IF
C=================== NEWTON-RAPHSON METHOD ITERATION ===================
C--------------------- IC = NUMBER OF ITERATIONS -----------------------
C---------------------- INITIAL GUESSES FOR TAU'S ----------------------
TAU(1) = DSQRT(2.D0) + A/3.D0
TAU(2) = DSQRT(3.D0) + B/TAU(1)
TAU(3) = DATAN(1.D0) + C/(TAU(1)*TAU(2))
C
DO I = 1 , MAXITA
IC = I
F(1) = TAU(1)+TAU(2)+TAU(3)-A
F(2) = TAU(1)*TAU(2)+TAU(2)*TAU(3)+TAU(3)*TAU(1)-B
F(3) = TAU(1)*TAU(2)*TAU(3)-C
C--------- DERIVATIVE MATRIX [DRV] ------
DRV(1,1) = 1.D0
DRV(1,2) = 1.D0
DRV(1,3) = 1.D0
DRV(2,1) = TAU(2)+TAU(3)
DRV(2,2) = TAU(1)+TAU(3)
DRV(2,3) = TAU(1)+TAU(2)
DRV(3,1) = TAU(2)*TAU(3)
DRV(3,2) = TAU(1)*TAU(3)
DRV(3,3) = TAU(1)*TAU(2)
C------- DETERMINANT OF [DRV]
DETJ = DRV(1,1) * ( DRV(2,2)*DRV(3,3)-DRV(2,3)*DRV(3,2) )
* - DRV(1,2) * ( DRV(2,1)*DRV(3,3)-DRV(2,3)*DRV(3,1) )
* + DRV(1,3) * ( DRV(2,1)*DRV(3,2)-DRV(2,2)*DRV(3,1) )
C------- [DI] = INVERSE OF THE MATRIX [DRV]
DI(1,1) = (DRV(2,2)*DRV(3,3) - DRV(2,3)*DRV(3,2))/DETJ
DI(2,1) = (DRV(2,3)*DRV(3,1) - DRV(2,1)*DRV(3,3))/DETJ
DI(3,1) = (DRV(2,1)*DRV(3,2) - DRV(2,2)*DRV(3,1))/DETJ
DI(1,2) = (DRV(1,3)*DRV(3,2) - DRV(1,2)*DRV(3,3))/DETJ
DI(2,2) = (DRV(1,1)*DRV(3,3) - DRV(1,3)*DRV(3,1))/DETJ
DI(3,2) = (DRV(1,2)*DRV(3,1) - DRV(1,1)*DRV(3,2))/DETJ
DI(1,3) = (DRV(1,2)*DRV(2,3) - DRV(1,3)*DRV(2,2))/DETJ
DI(2,3) = (DRV(1,3)*DRV(2,1) - DRV(1,1)*DRV(2,3))/DETJ
DI(3,3) = (DRV(1,1)*DRV(2,2) - DRV(1,2)*DRV(2,1))/DETJ
C-------- NEW VALUES FOR TAU'S
TAU(1) = TAU(1) - ( DI(1,1)*F(1)+DI(1,2)*F(2)+DI(1,3)*F(3) )
TAU(2) = TAU(2) - ( DI(2,1)*F(1)+DI(2,2)*F(2)+DI(2,3)*F(3) )
TAU(3) = TAU(3) - ( DI(3,1)*F(1)+DI(3,2)*F(2)+DI(3,3)*F(3) )
C--------- ENERGY BY COMPUTAED PRINCIPAL STRESSES ------
ENERGY1 = 2.D0*(TAU(1)+TAU(2)+TAU(3))**2-
* 6.D0*(TAU(1)*TAU(2)+TAU(2)*TAU(3)+TAU(3)*TAU(1))
C--------- DECISION --------
IF ( DABS(ENERGY1-ENERGY0)/ENERGY0 .LE. EPS ) EXIT
END DO
C--------------------------- NO CONVERGENCE ----------------------------
IF ( I .GE. MAXITA ) THEN
IC = -4
END IF
C===================== EIGEN VECTOR COMPUTATION ========================
CALL EVECTOR (EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, ICV, TAU, U )
IF ( IC .GT. 0 ) IC = IC + ICV
C============ MAKING TAU(i) AND U(i,j) DESCENDING ORDER ================
CALL ORDERING ( TAU, U )
RETURN
END
C
C
SUBROUTINE CARDANO (EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, IC, TAU, U )
C=========== SOLUTION OF CUBIC EQUATION BY CARDANO METHOD ==============
C-------------------- X**3 + Ax**2 + Bx +C = 0 -------------------------
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION TAU(3), U(3,3)
COMPLEX*16 T1, T2, ALPHA, X1, X2, X3, BETA
C------------------------ STRESS INVARIANTS ----------------------------
A = -(TXX+TYY+TZZ)
B = TXX*TYY+TYY*TZZ+TZZ*TXX-TXY**2-TYZ**2-TXZ**2
C = -(TXX*TYY*TZZ+2.D0*TXY*TXZ*TYZ-
* TXY**2*TZZ-TYZ**2*TXX-TXZ**2*TYY)
C------------------------- INITIAL SETTING -----------------------------
DO I = 1 , 3
DO J = 1 , 3
U(I,J) = 0.D0
END DO
U(I,I) = 1.D0
END DO
C------------------------- CASE OF A=B=C=0 -----------------------------
IF ( A.EQ.0.D0 .AND. B.EQ.0.D0 .AND. C.EQ.0.D0 ) THEN
IC = 0
TAU(1) = 0.D0
TAU(2) = 0.D0
TAU(3) = 0.D0
RETURN
END IF
C------------------------ VARIABLES P AND Q ----------------------------
P = -A**2/3.D0 + B
Q = 2.D0*A**3/27.D0 -A*B/3.D0 + C
C--------------- CASE OF P=Q=0 MEANING TRIPLE ROOTS---------------------
IF ( P.EQ.0.D0 .AND. Q.EQ.0.D0 ) THEN
IC = -3
TAU(1) = TXX
TAU(2) = TYY
TAU(3) = TZZ
RETURN
END IF
C------------------- COMPUTATION OF QUBIC ROOTS ------------------------
DISCR = Q**2 + (4.D0/27.D0)*P**3
BETA = (-DCMPLX(Q) - CDSQRT(DCMPLX(DISCR)))/2.0D0
ALPHA = BETA**(1.D0/3.D0)
T1 = DCMPLX(-0.5D0, +0.5D0*DSQRT(3.D0))
T2 = DCMPLX(-0.5D0, -0.5D0*DSQRT(3.D0))
X1 = ALPHA -DCMPLX(P)/(3.D0*ALPHA) -DCMPLX(A)/3.D0
X2 = ALPHA*T1-DCMPLX(P)/(3.D0*ALPHA*T1)-DCMPLX(A)/3.D0
X3 = ALPHA*T2-DCMPLX(P)/(3.D0*ALPHA*T2)-DCMPLX(A)/3.D0
C----------------- ALL PRINCIPAL STRESSES ARE REAL.---------------------
TAU(1) = DREAL(X1)
TAU(2) = DREAL(X2)
TAU(3) = DREAL(X3)
C====================== EIGEN VECTOR COMPUTATION =======================
CALL EVECTOR (EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, IC, TAU, U )
C========== MAKING TAU(i) ALONG WITH U(i,j) DESCENDING ORDER ===========
CALL ORDERING ( TAU, U )
RETURN
END
C
C
SUBROUTINE POWER (EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ,NITERA,TAU,U)
C---------- COMPUTATION OF PRINCIPAL STRESSES BY POWER METHOD ----------
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( MXITERA=99999 )
DIMENSION A(3,3), X(3),Y(3),U(3,3),TAU(3), UTEMP(3)
N = 3
C--------------- CASE OF P=Q=0 MEANING TRIPLE ROOTS---------------------
AA = -(TXX+TYY+TZZ)
B = TXX*TYY+TYY*TZZ+TZZ*TXX-TXY**2-TYZ**2-TXZ**2
C = -(TXX*TYY*TZZ+2.D0*TXY*TXZ*TYZ-
* TXY**2*TZZ-TYZ**2*TXX-TXZ**2*TYY)
P = -AA**2/3.D0 + B
Q = 2.D0*AA**3/27.D0 -AA*B/3.D0 + C
IF ( P.EQ.0.D0 .AND. Q.EQ.0.D0 ) THEN
TAU(1) = TXX
TAU(2) = TYY
TAU(3) = TZZ
DO I=1,3
DO J=1,3
U(I,J) = 0.0D0
END DO
U(I,I) = 1.0D0
END DO
NITERA = -3
RETURN
END IF
C----------------- SUBSTITUTION OF TAU-IJ INTO A(I,J) ------------------
A(1,1) = TXX
A(2,2) = TYY
A(3,3) = TZZ
A(1,2) = TXY
A(1,3) = TXZ
A(2,3) = TYZ
A(2,1) = A(1,2)
A(3,1) = A(1,3)
A(3,2) = A(2,3)
C----------------- CASE OF ALL STRESSES ARE ZERO. IC=0 ----------------
SUM = 0.D0
DO I = 1 , N
DO J = 1 , N
SUM = SUM + DABS(A(I,J))
END DO
END DO
IF ( SUM .EQ. 0.D0 ) THEN
TAU(1) = 0.D0
TAU(2) = 0.D0
TAU(3) = 0.D0
IC = 0
RETURN
END IF
C================ ADDING BIAS STRESS INTO NORMAL STRESSES ==============
TAUMIN = DABS(A(1,1))+DABS(A(1,2))+DABS(A(1,3))
TAUMIN = DMAX1 (TAUMIN,DABS(A(2,1))+DABS(A(2,2))+DABS(A(2,3)))
TAUMIN = DMAX1 (TAUMIN,DABS(A(3,1))+DABS(A(3,2))+DABS(A(3,3)))
A(1,1) = A(1,1) + 2.D0*TAUMIN
A(2,2) = A(2,2) + 2.D0*TAUMIN
A(3,3) = A(3,3) + 2.D0*TAUMIN
C============== PRINCIPAL STRESS COMPUTATION STARTS HERE ===============
NITERA = 0
DO MODE = 1 , 3
IF ( MODE .GT. 1 ) THEN
C-------------------- FOR SECOND AND THIRD MODES -----------------------
A(1,1) = A(1,1) - TAU(MODE-1)*U(1,MODE-1)*U(1,MODE-1)
A(1,2) = A(1,2) - TAU(MODE-1)*U(1,MODE-1)*U(2,MODE-1)
A(1,3) = A(1,3) - TAU(MODE-1)*U(1,MODE-1)*U(3,MODE-1)
A(2,1) = A(2,1) - TAU(MODE-1)*U(2,MODE-1)*U(1,MODE-1)
A(2,2) = A(2,2) - TAU(MODE-1)*U(2,MODE-1)*U(2,MODE-1)
A(2,3) = A(2,3) - TAU(MODE-1)*U(2,MODE-1)*U(3,MODE-1)
A(3,1) = A(3,1) - TAU(MODE-1)*U(3,MODE-1)*U(1,MODE-1)
A(3,2) = A(3,2) - TAU(MODE-1)*U(3,MODE-1)*U(2,MODE-1)
A(3,3) = A(3,3) - TAU(MODE-1)*U(3,MODE-1)*U(3,MODE-1)
END IF
C----------- INITIAL SETTING FOR ITERATION PROCEDURE -------------------
U(1,MODE) = 0.D0
U(2,MODE) = 0.D0
U(3,MODE) = 0.D0
U(MODE,MODE) = 1.D0
DIFF = 1.D0
C--------------- ITERATION FOR MODE=I STARTS HERE ----------------------
DO WHILE ( DIFF .GT. EPS .AND. NITERA .LE. MXITERA )
NITERA = NITERA + 1
C-------- COMPUTATION OF [A]{X}={Y}
Y(1) = A(1,1)*U(1,MODE) + A(1,2)*U(2,MODE) + A(1,3)*U(3,MODE)
Y(2) = A(2,1)*U(1,MODE) + A(2,2)*U(2,MODE) + A(2,3)*U(3,MODE)
Y(3) = A(3,1)*U(1,MODE) + A(3,2)*U(2,MODE) + A(3,3)*U(3,MODE)
C-------- COMPUTATION OF PRINCIPAL STRESSES
TAU(MODE) = DSQRT(Y(1)*Y(1)+Y(2)*Y(2)+Y(3)*Y(3))
C-------- COMPUTATION OF ERROR
IF ( TAU(MODE) .NE. 0.D0 ) THEN
X(1) = Y(1)/TAU(MODE)
X(2) = Y(2)/TAU(MODE)
X(3) = Y(3)/TAU(MODE)
DIFF=(DABS(U(1,MODE)-X(1))+DABS(U(2,MODE)-X(2))+
* DABS(U(3,MODE)-X(3)))/3
U(1,MODE) = X(1)
U(2,MODE) = X(2)
U(3,MODE) = X(3)
ELSE
DIFF = 0.D0
U(1,MODE) = 0.D0
U(2,MODE) = 0.D0
U(3,MODE) = 0.D0
U(MODE,MODE) = 1.D0
END IF
C----------------- END OF ITERATION FOR MODE=I -------------------------
END DO
C------------------------ END OF DO MODE -------------------------------
END DO
C============== SUBTRACT BIAS FROM NORMAL STRESSES =====================
TAU(1) = TAU(1) - 2.D0*TAUMIN
TAU(2) = TAU(2) - 2.D0*TAUMIN
TAU(3) = TAU(3) - 2.D0*TAUMIN
C======== MAKING TAU(1), TAU(2), AND TAU(3) DESCENDING ORDER ===========
CALL ORDERING ( TAU, U )
RETURN
END
C
C
SUBROUTINE BSECTION ( EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, IC,TAU,U )
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION P(3), TAU(3), U(3,3)
LOGICAL DCSN2, DCSN5
MAXITA = 1000
WRITE (2,*) '@BSECTION2'
WRITE (2,*) 'MODE BISECTION-ITERATION LAMBDA ERROR'
C--------
NEIGEN = 3
C-------- COMPUTATION OF LIMITS WHERE EIGENVALUES EXIST
TAUMAX = DABS(TXX)+DABS(TXY)+DABS(TXZ)
TAUMAX = DMAX1 (TAUMAX,DABS(TXY)+DABS(TYY)+DABS(TYZ))
TAUMAX = DMAX1 (TAUMAX,DABS(TXZ)+DABS(TYZ)+DABS(TZZ))
C-------- COMPUTATION OF IGENVALUES ----------
DO MODE = 1, NEIGEN
ITA = 0
BOTTOLMT =-TAUMAX*1.1D0
UPPERLMT = TAUMAX*1.2D0
FMAG = 1.D0
C------------ INITIAL SETTING FOR CONVERGENCE PARAMETERS ---------------
DCSN2 = .TRUE.
DCSN5 = .TRUE.
C*** BEGIN OF DO WHILE ***
DO WHILE ( DCSN2 .AND. DCSN5 )
DCSN2 = DABS(UPPERLMT-BOTTOLMT)/FMAG .GT. EPS
ITA = ITA + 1
DCSN5 = ITA .LT. MAXITA
FLMD = 0.5D0*(BOTTOLMT+UPPERLMT)
C------- FORMING STURM SEQUENCES
C------- P(0), P(1), P(2), P(3). NOTE: P(0)=1.
P(1) = TXX-FLMD
P(2) = (TXX-FLMD)*(TYY-FLMD)-TXY**2
POSI = (TXX-FLMD)*(TYY-FLMD)*(TZZ-FLMD)+2.D0*TXY*TYZ*TXZ
ANEG = (TYY-FLMD)*TXZ**2+(TXX-FLMD)*TYZ**2+(TZZ-FLMD)*TXY**2
P(3) = POSI - ANEG
C------- COUNT NUMBER OF SIGN-CHANGES IN P(I)
KOUNT = 0
PRODUCT = 1.D0
DO J = 1 , NEIGEN
IF ( PRODUCT*P(J) .LT. 0.D0 ) THEN
KOUNT = KOUNT + 1
PRODUCT = - PRODUCT
END IF
END DO
C------- NARROW DOWN LIMTS WHERE EIGENVALUE OF MODE EXISTS
IF ( KOUNT .LT. MODE ) THEN
BOTTOLMT = FLMD
ELSE
UPPERLMT = FLMD
END IF
FMAG = DMAX1 (DABS(UPPERLMT), DABS(BOTTOLMT))
IF ( FMAG .EQ. 0.D0 ) FMAG = 1.D0
END DO
C*** END OF DO WHILE ***
TAU(MODE) = FLMD
WRITE (2,*) MODE,ITA,FLMD,DABS(UPPERLMT-BOTTOLMT)
END DO
C====================== EIGEN VECTOR COMPUTATION =======================
CALL EVECTOR (EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, IC, TAU, U )
IC = IC + ITA
C========== MAKING TAU(i) ALONG WITH U(i,j) DESCENDING ORDER ===========
CALL ORDERING ( TAU, U )
RETURN
END
C
C
SUBROUTINE EVECTOR (EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, IC, TAU, U )
C========= SOLUTION OF EIGENVECTORS BY INVERSE POWER METHOD ============
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION TAU(3), F(3,3), G(3,3), U(3,3), TEMPU(3)
DATA MAXITA / 100 /
C
DO MODE = 1 , 3
F(1,1) = TXX-TAU(MODE)*(1.D0-EPS)
F(2,2) = TYY-TAU(MODE)*(1.D0-EPS)
F(3,3) = TZZ-TAU(MODE)*(1.D0-EPS)
F(1,2) = TXY
F(1,3) = TXZ
F(2,3) = TYZ
F(2,1) = F(1,2)
F(3,1) = F(1,3)
F(3,2) = F(2,3)
C----------- DETERMINANT OF [F]=[TAU(I,J)]-TAUI(1-EPS)[I] --------------
POSI = F(1,1)*F(2,2)*F(3,3)+F(2,1)*F(3,2)*F(1,3)+
* F(3,1)*F(1,2)*F(2,3)
ANEG = F(1,3)*F(2,2)*F(3,1)+F(1,1)*F(2,3)*F(3,2)+
* F(2,1)*F(1,2)*F(3,3)
DET = POSI - ANEG
C------- 1ST ROW OF INVERSE OF [F]
G(1,1) = (F(2,2)*F(3,3)-F(2,3)*F(3,2))/DET
G(1,2) = -(F(2,1)*F(3,3)-F(2,3)*F(3,1))/DET
G(1,3) = (F(2,1)*F(3,2)-F(2,2)*F(3,1))/DET
C------- 2ND ROW OF INVERSE OF [F]
G(2,1) = -(F(1,2)*F(3,3)-F(1,3)*F(3,2))/DET
G(2,2) = (F(1,1)*F(3,3)-F(1,3)*F(3,1))/DET
G(2,3) = -(F(1,1)*F(3,2)-F(1,2)*F(3,1))/DET
C------- 3RD ROW OF INVERSE OF [F]
G(3,1) = (F(1,2)*F(2,3)-F(1,3)*F(2,2))/DET
G(3,2) = -(F(1,1)*F(2,3)-F(1,3)*F(2,1))/DET
G(3,3) = (F(1,1)*F(2,2)-F(1,2)*F(2,1))/DET
C------- INVERSE OF [F] TIMES TAU*EPS
DO I = 1 , 3
DO J = 1 , 3
G(I,J) = TAU(MODE)*EPS*G(I,J)
END DO
END DO
C------------- INITIAL SETTING OF EIGEN VECTOR U(J,MODE) ---------------
U(1,MODE) = 0.D0
U(2,MODE) = 0.D0
U(3,MODE) = 0.D0
U(MODE,MODE) = 1.D0
C============================= ITERATION ===============================
VECL = 9999.D0
DO ITERATIN = 1 , MAXITA
TEMPU(1) = G(1,1)*U(1,MODE)+G(1,2)*U(2,MODE)+G(1,3)*U(3,MODE)
TEMPU(2) = G(2,1)*U(1,MODE)+G(2,2)*U(2,MODE)+G(2,3)*U(3,MODE)
TEMPU(3) = G(3,1)*U(1,MODE)+G(3,2)*U(2,MODE)+G(3,3)*U(3,MODE)
C--------------------------- VECTOR LENGTH -----------------------------
VECLENGH = DSQRT(TEMPU(1)**2+TEMPU(2)**2+TEMPU(3)**2)
C--------------------------- SUBSTITUTION ------------------------------
U(1,MODE) = TEMPU(1)/VECLENGH
U(2,MODE) = TEMPU(2)/VECLENGH
U(3,MODE) = TEMPU(3)/VECLENGH
C------------------------- CHECK CONVERGENCE ---------------------------
DELTA = DABS(VECL-VECLENGH)/DMAX1(VECL,VECLENGH)
IF ( DELTA .LE. EPS ) EXIT
VECL = VECLENGH
END DO
C------------------------ BELOW END OF DO MODE -------------------------
END DO
IC = ITERATIN
RETURN
END
C
C
SUBROUTINE ORDERING ( TAU, U )
C======== MAKING TAU(1), TAU(2), AND TAU(3) DESCENDING ORDER ===========
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION U(3,3),TAU(3), UTEMP(3)
IF ( TAU(2) .GT. TAU(1) ) THEN
TEMP = TAU(1)
TAU(1) = TAU(2)
TAU(2) = TEMP
DO J = 1 , 3
UTEMP(J) = U (J,1)
U (J,1) = U (J,2)
U (J,2) = UTEMP(J)
END DO
END IF
IF ( TAU(3) .GT. TAU(2) ) THEN
TEMP = TAU(2)
TAU(2) = TAU(3)
TAU(3) = TEMP
DO J = 1 , 3
UTEMP(J) = U (J,2)
U (J,2) = U (J,3)
U (J,3) = UTEMP(J)
END DO
END IF
IF ( TAU(2) .GT. TAU(1) ) THEN
TEMP = TAU(1)
TAU(1) = TAU(2)
TAU(2) = TEMP
DO J = 1 , 3
UTEMP(J) = U (J,1)
U (J,1) = U (J,2)
U (J,2) = UTEMP(J)
END DO
END IF
RETURN
END