子程序subroutine vumat
subroutine vumat(
C Read only -
1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
2 stepTime, totalTime, dt, cmname, coordMp, charLength,
3 props, density, strainInc, relSpinInc,
4 tempOld, stretchOld, defgradOld, fieldOld,
3 stressOld, stateOld, enerInternOld, enerInelasOld,
6 tempNew, stretchNew, defgradNew, fieldNew,
C Write only -
5 stressNew, stateNew, enerInternNew, enerInelasNew )
include 'vaba_param.inc'
C
C MODEL FOR THE SHEAR DEFORMATION AN
D CRACK OF POL YMERS
C SHEAR MODEL: SPRING+BURGERS ELEMENT+PLASTIC ELEMENT(D-P & EYRING)
C CRACK MODEL: S11>S11_CR WHEN SKK>0
C
C WRITTEN BY ZHANG CHUNYU(CHUNYU@https://www.sodocs.net/doc/972674948.html,.SG)
C
C All arrays dimensioned by (*) are not used in this algorithm
dimension props(nprops), density(nblock),
1 coordMp(nblock,*),
2 charLength(*), strainInc(nblock,ndir+nshr),
3 relSpinInc(*), tempOld(*),
4 stretchOld(*), defgradOld(nblock,ndir+nshr),
5 fieldOld(*), stressOld(nblock,ndir+nshr),
6 stateOld(nblock,nstatev), enerInternOld(nblock),
7 enerInelasOld(nblock), tempNew(*),
8 stretchNew(*), defgradNew(nblock,ndir+nshr), fieldNew(*),
9 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
1 enerInternNew(nblock), enerInelasNew(nblock)
C strain components stored as state variables
dimension eelas(ndir+nshr),eplas(ndir+nshr),deplas(ndir+nshr)
dimension evisco(ndir+nshr),devisco(ndir+nshr),veint(ndir+nshr)
dimension effstrn(ndir+nshr),ps(ndir),an(ndir,ndir),s(ndir+nshr)
dimension dr(3,3),f_old(6),f_new(6)
data newton,toler/10,1.0E-6/
C
character*80 cmname
C
parameter( zero = 0., one = 1., two = 2., three = 3.,
1 third = one/three, half = .5, twoThirds = two/three,
2 threeHalfs = 1.5 )
C -----------------------------------------------------------
C PROPS(1) - E
C PROPS(2) - NU
C PROPS(3) - E1
C PROPS(4) - ETA1
C PROPS(5) - ETA0
C PROPS(6) - PHAI
C PROPS(7) - PHAI2 FOR NON-ASSOCIATIVE FLOW
C PROPS(8) - B FOR EYRING MODEL
C PROPS(9) - C1
C PROPS(10)- C2
C PROPS(11)- M
C PROPS(12)- SCRAZE
C PROPS(13)- ZETA0
C PROPS(14)- CRAZE STRAIN LIMIT
C PROPS(15)- SHEAR STRAIN LIMIT
C PROPS(16)- Y0
C PROPS(17)... PLASTIC CURVE
C CALLS AHAR
D FOR CURV
E O
F SYIELD VS. PEEQ
C -----------------------------------------------------------
C
if (ndir .NE. 3) then
write(6,1)
1 format(//,30X,'***ERROR - THIS VUMA T MAY ONL Y BE USED FOR ', 1 'ELEMENTS WITH THREE DIRECT STRESS COMPONENTS')
endif
C
C ELASTIC PROPERTIES
ntens=ndir+nshr
bniu=props(2)
if(bniu.GT.0.4999.AND.bniu.LT.0.5001) bniu=0.499
C PARAMETERS FOR THE KELVIN UNIT
phai=props(6)
phai2=props(7)
B=props(8)
c1=props(9)
c2=props(10)
c3=three*bniu/(one+bniu)
bm=props(11)
scraze=props(12)
zeta0=props(13)
craze_cr=props(14)
shear_cr=props(15)
c loop the material block
do 350 i = 1,nblock
c_factor=stateOld(i,4*ntens+4)
if(c_factor .EQ. 0) c_factor=1.0
q_factor=c_factor
e0=props(1)
e1=props(3)
eta1=props(4)
eta0=props(5)
e0=e0*q_factor
bk=e0/(one-two*bniu)/three
g=e0/(one+bniu)/two
e1=e1*q_factor
g1=e1/(2.0*(1.0+bniu))
eta1=eta1*q_factor
eta0=eta0*q_factor
coef1=1.0+eta1/eta0+1.0/4.0*g1/eta0*dt
coef2=g1+2.0*eta1/dt
coef3=1.0/4.0*g1/eta0*dt
coef=1.0+g*coef1/coef2
C RECOVER AN
D ROTAT
E DEVIATORIC SHEAR STRAIN & TOTAL STRESS c CALL ROTSIG(STATEV(1),DROT,EELAS,2,NDI,NSHR)
c CALL ROTSIG(STA TEV(NTENS+1),DROT,EPLAS,2,NDI,NSHR)
c CALL ROTSIG(STATEV(2*NTENS+3),DROT,EVISCO,2,NDI,NSHR)
do 10 j=1,ntens
eelas(j)=stateOld(i,j)
eplas(j)=stateOld(i,j+ntens)
evisco(j)=stateOld(i,j+3+2*ntens)
veint(j)=stateOld(i,j+3+3*ntens)
deplas(j)=0.0
f_old(j)=defgradOld(i,j)
f_new(j)=defgradNew(i,j)
10 continue
c get rotation tensor an
d rotat
e old strain tensor
call getdr(f_old,f_new,dr)
call ROTTENSOR(eelas,dr,3,3)
call ROTTENSOR(eplas,dr,3,3)
call ROTTENSOR(evisco,dr,3,3)
eqplas_shr=stateOld(i,1+2*ntens)
eqplas_crz=stateOld(i,2+2*ntens)
vplas=stateOld(i,3+2*ntens)
C CALCULATE DEVIATORIC STRESS
oldstresskk=0.0
dstrankk=0.0
strankk=0.0
do 20 j=1,ndir
oldstresskk=oldstresskk+stressOld(i,j)
dstrankk=dstrankk+strainInc(i,j)
strankk=strankk+eelas(j)+eplas(j)+evisco(j)+strainInc(i,j)
20 continue
C CALCULATE EFFECTIVE STRAIN
do 30 j=1,ndir
effstrn(j)=eelas(j)+(strainInc(i,j)-dstrankk/3.0)-
1 (0.5*(coef1+2.0*coef3)*(stressOld(i,j)-oldstresskk/3.0)+
1 g1/eta0*veint(j)-2.0*g1*evisco(j))/coef2
30 continue
do 31 j=ndir+1,ntens
effstrn(j)=eelas(j)+strainInc(i,j)/2.0-
1 (0.5*(coef1+2.0*coef3)*stressOld(i,j)+
1 g1/eta0*veint(j)-2.0*g1*evisco(j))/coef2
31 continue
C UPDATE
D DEVIATORIC STRESS
do 40 j=1,ntens
stressNew(i,j)=2.0*g/coef*effstrn(j)
40 continue
C CALCULATE BULK STRESS & PRESSURE
stresskk=oldstresskk+3.0*bk*dstrankk
p=-stresskk/3.0
C UPDATE DIRECT STRESS
do 45 j=1,ndir
stressNew(i,j)=stressNew(i,j)-p
45 continue
C DETERMINE SHEAR YILE
D OR CRAZ
E YIELD
do 50 j=1,ntens
s(j)=stressNew(i,j)
50 continue
C PRINCIPAL STRESS AND
call PRINCIPAL(s,ps,an,ndir,nshr,1.0E-6)
ps1=ps(1)
nmax=1
if(ps(2) .GT. ps1) then
ps1=ps(2)
nmax=2
endif
if(ps(3) .GT. ps1) then
ps1=ps(3)
nmax=3
endif
bms=-p
cs=c1+c2/bms+c3*bms
if (bms .GT. 0 .AND. ps1 .GT. 0 .AND. ps1 .GT. cs) then
flow=one
else
flow=zero
endif
C CRAZE FLOW,calculate deplas
if (flow .EQ. one .AND. c_factor .GT. 1.0E-4) then
do 60 j=1,ndir
deplas(j)=zeta0*(ps1/scraze/q_factor)**(1.0/bm)*dt*
1 an(j,nmax)**2
60 continue
do 61 j=ndir+1,ntens
n1=j-ndir
n2=j-ndir+1
if(n2 .GT. ndir) then
n2=n2-ndir
endif
deplas(j)=zeta0*(ps1/scraze/q_factor)**(1.0/bm)*dt*
1 an(n1,nmax)*an(n2,nmax)
61 continue
eqplas_crz=eqplas_crz+zeta0*(ps1/scraze/q_factor)**(1.0/bm)*dt if(eqplas_crz .GT. craze_cr .AND. c_factor .GT. zero) then
c_factor=c_factor*0.1*craze_cr/eqplas_crz
endif
c update stress
dvplaskk=deplas(1)+deplas(2)+deplas(3)
vplas=vplas+dvplaskk
stresskk=oldstresskk+3.0*bk/q_factor*c_factor*
1 (dstrankk-dvplaskk)
do 70 j=1,ndir
deplas(j)=deplas(j)-dvplaskk/3.0
stressNew(i,j)=2.0*g/q_factor*c_factor*
1 (effstrn(j)-deplas(j))/coef
stressNew(i,j)=stressNew(i,j)+stresskk/3.0
70 continue
do 71 j=ndir+1,ntens
stressNew(i,j)=2.0*g/q_factor*c_factor*
1 (effstrn(j)-deplas(j))/coef
71 continue
endif
C SHEAR YIEL
D FLOW,CALCULAT
E DEPLAS
c test crack-->shear yielding
c if crack an
d shear yielding can not coexist,delet
e the following the command c flow=zero
if(nprops.GT.15 .AND. props(16).GT. 0.0 .AND. flow .EQ. zero) then
C
C MISES STRESS
C
smises=(s(1)-s(2))*(s(1)-s(2)) +
1 (s(2)-s(3))*(s(2)-s(3)) +
1 (s(3)-s(1))*(s(3)-s(1))
do 90 j=ndir+1,ntens
smises=smises+6.0*s(j)*s(j)
90 CONTINUE
smises=sqrt(smises/2.0)
C SHEAR STRAIN RATE
effrate=0.0
do 95 j=1,ndir
effrate=effrate+(strainInc(i,j)-dstrankk/3.0)*
1 (strainInc(i,j)-dstrankk/3.0)
95 continue
do 100 j=ndir+1,ntens
effrate=effrate+2.0*strainInc(i,j)*strainInc(i,j)
100 CONTINUE
effrate=sqrt(effrate/2.0)/dt
if (effrate .LT. 1.0E-6) then
effrate=1.0E-6
endif
C HARDENING CURVE, GET YIEL
D STRESS
C CALCULATE HOW MANY POINTS FOR THE PLASTIC CURVE
C COEF5 DEPQ=COEF5*LAMTA
coef5=sqrt(1.0+phai2*phai2/2.0)
NV ALUE=nprops/2-7
CALL AHARD(syiel0,HARD,eqplas_shr*coef5,props(16),NVALUE)
C DETERMINE IF ACTIVEL Y YIELDING
C
effyield=syiel0+B*log10(effrate)+phai*p
if (smises .GT. (1+toler)*effyield ) then
C CALCULATE EQUIV ALENT STRAIN
equie=0.0
do 110 j=1,ndir
equie=equie+effstrn(j)*effstrn(j)
110 continue
do 111 j=ndir+1,ntens
equie=equie+2*effstrn(j)*effstrn(j)
111 continue
equie=sqrt(2.0/3.0*equie)
C
C SOLVE FOR EQUIV STRESS, NEWTON ITERATION
C Maybe unnecessary for explicit solver
syield=syiel0
deqpl=0.0
do 130 j=1,newton
rhs=3.0*g*(equie-deqpl)-coef*(syield+B*log10(effrate)
1 -phai*bk*(strankk-vplas-phai2*deqpl))
deqpl=deqpl+rhs/(3.0*g+coef*(HARD*coef5+bk*phai*phai2)) call AHARD(syield,HARD,(eqplas_shr+deqpl)*coef5,props(16),NV ALUE)
if(abs(rhs).LT.toler*effyield) goto 140
130 continue
write(6,2) newton
2 format(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT ', 1 'CONVERGE AFTER ',I3,' ITERATIONS')
140 continue
stresskk=3.0*bk*(strankk-vplas-phai2*deqpl)
p=-stresskk/3.0
effyield=syield+phai*p+B*log10(effrate)
coef4=3.0*g*deqpl/effyield+coef
do 150 j=1,ndir
stressNew(i,j)=2.0*g*effstrn(j)/coef4
deplas(j)=3.0*stressNew(i,j)*deqpl/(2.0*effyield)
stressNew(i,j)=stressNew(i,j)+stresskk/3.0
150 continue
do 160 j=ndir+1,ntens
stressNew(i,j)=2.0*g*effstrn(j)/coef4
deplas(j)=3.0*stressNew(i,j)*deqpl/(2.0*effyield)
160 continue
eqplas_shr=eqplas_shr+deqpl
vplas=vplas+phai2*deqpl
if(eqplas_shr .GT. shear_cr .AND. c_factor .GT. zero) then c_factor=c_factor*0.1*shear_cr/eqplas_shr
endif
endif
endif
C UPDATE STATE V ARIABLES
C CALCULATE DEVISCO
do 260 j=1,ndir
devisco(j)=(0.5*coef1*(stressOld(i,j)-oldstresskk/3.0+
1 stressNew(i,j)-stresskk/3.0)+coef3*(stressOld(i,j)-
1 oldstresskk/3.0)+g1/eta0*veint(j)-2*g1*evisco(j))/coef2
260 continue
do 265 j=ndir+1,ntens
devisco(j)=(0.5*coef1*(stressOld(i,j)+stressNew(i,j))+
1 coef3*stressOld(i,j)+g1/eta0*veint(j)-
1 2*g1*evisco(j))/coef2
265 continue
C STORE STATE V ARIABLE ARRAY
do 300 j=1,ndir
stateNew(i,j)=eelas(j)+strainInc(i,j)-dstrankk/3.0-devisco(j)-
1 deplas(j)
stateNew(i,j+ntens)=eplas(j)+deplas(j)
stateNew(i,j+3+2*ntens)=evisco(j)+devisco(j)
stateNew(i,j+3+3*ntens)=veint(j)+0.5*dt*(stressOld(i,j)-
1 oldstresskk/3.0+stressNew(i,j)-stresskk/3.0)
300 continue
do 310 j=ndir+1,ntens
stateNew(i,j)=eelas(j)+strainInc(i,j)/2.0-devisco(j)-deplas(j) stateNew(i,j+ntens)=eplas(j)+deplas(j)
stateNew(i,j+3+2*ntens)=evisco(j)+devisco(j)
stateNew(i,j+3+3*ntens)=veint(j)+0.5*dt*(stressOld(i,j)+
1 +stressNew(i,j))
310 continue
stateNew(i,1+2*ntens)=eqplas_shr
stateNew(i,2+2*ntens)=eqplas_crz
stateNew(i,3+2*ntens)=vplas
stateNew(i,4*ntens+4)=c_factor
if(c_factor .LT. 1.0E-4) then
stateNew(i,4*ntens+5)=zero
else
stateNew(i,4*ntens+5)=one
endif
350 continue
C
return
end
C
c CALCULATE EQUIV ALENT STRESS
SUBROUTINE AHARD(SYIELD,HARD,EQPLAS,TABLE,NV ALUE)
C
INCLUDE 'ABA_PARAM.INC'
DIMENSION TABLE(2,NV ALUE)
C
C SET YIEL
D STRESS TO LAST V ALU
E O
F TABLE, HARDENIN
G TO ZERO
SYIELD=TABLE(1,NV ALUE)
HARD=0.0
C
C IF MORE THAN ONE ENTRY, SEARCH TABLE
C
IF(NV ALUE.GT.1) THEN
DO 10 K1=1,NV ALUE-1
EQPL1=TABLE(2,K1+1)
IF(EQPLAS.LT.EQPL1) THEN
EQPL0=TABLE(2,K1)
IF(EQPL1.LE.EQPL0) THEN
WRITE(6,1)
1 FORMAT(//,30X,'***ERROR - PLASTIC STRAIN MUST BE ',
1 'ENTERED IN ASCENDING ORDER')
C CALL XIT
ENDIF
C
C CURRENT YIEL
D STRESS AND HARDENING
C
DEQPL=EQPL1-EQPL0
SYIEL0=TABLE(1,K1)
SYIEL1=TABLE(1,K1+1)
DSYIEL=SYIEL1-SYIEL0
HARD=DSYIEL/DEQPL
SYIELD=SYIEL0+(EQPLAS-EQPL0)*HARD
GOTO 20
ENDIF
10 CONTINUE
20 CONTINUE
ENDIF
RETURN
END
C ROUNTINE TO CALCULATE EIGENV ALUES AN
D EIGENVECTORS
SUBROUTINE PRINCIPAL(S,PS,V,NDIR,NSHR,TL)
C IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION A(NDIR,NDIR),V(NDIR,NDIR),S(NDIR+NSHR),PS(NDIR)
NEQ=NDIR
C EIGENV ALUE SOLUTION BY JACOBI METHO
D -
C A - MATRIX (ANY RANK) TO BE SOLVE
D ---
C EIGENV ALUES ON DIAGONAL
C V - MA TRIX OF EIGENVECTORS PRODUCED
C TL- NUMBER OF SIGNIFICANT FIGURES
C---- INITIALIZATION -----------------------
ZERO = 0.0D0
SUM = ZERO
TOL = ABS(TL)
A(1,1)=S(1)
A(2,2)=S(2)
A(3,3)=S(3)
A(1,2)=S(NDIR+1)
A(2,1)=S(NDIR+1)
IF(NSHR .GT. 1) THEN
A(1,3)=S(NDIR+3)
A(3,1)=S(NDIR+3)
A(2,3)=S(NDIR+2)
A(3,2)=S(NDIR+2)
ENDIF
C---- SET INITIAL EIGENVECTORS -------------
DO 200 I=1,NEQ
DO 190 J=1,NEQ
IF (TL.GT.ZERO) V(I,J) = ZERO
190 SUM = SUM + ABS(A(I,J))
IF (TL.GT.ZERO) V(I,I) = 1.0
200 CONTINUE
C---- CHECK FOR TRIVIAL PROBLEM -----------
IF (NEQ.EQ.1) RETURN
IF (SUM.LE.ZERO) RETURN
SUM = SUM/FLOAT(NEQ*NEQ)
C-------------------------------------------
C---- REDUCE MA TRIX TO DIAGONAL ------------
C-------------------------------------------
400 SSUM = ZERO
AMAX = ZERO
DO 700 J=2,NEQ
IH = J-1
DO 700 I=1,IH
C---- CHECK IF A(I,J) IS TO BE REDUCED -----
AA = ABS(A(I,J))
IF (AA.GT.AMAX) AMAX = AA
SSUM = SSUM + AA
IF (AA.LT.0.1*AMAX) GO TO 700
C---- CALCULATE ROTA TION ANGLE ----------
AA=ATAN2(2.0*A(I,J),A(I,I)-A(J,J))/2.0
SI = SIN(AA)
CO = COS(AA)
C---- MODIFY "I" AND "J" COLUMNS --------
DO 500 K=1,NEQ
TT = A(K,I)
A(K,I) = CO*TT + SI*A(K,J)
A(K,J) = -SI*TT + CO*A(K,J)
TT = V(K,I)
V(K,I) = CO*TT + SI*V(K,J)
500 V(K,J) = -SI*TT + CO*V(K,J)
C---- MODIFY DIAGONAL TERMS -------------
A(I,I) = CO*A(I,I) + SI*A(J,I)
A(J,J) =-SI*A(I,J) + CO*A(J,J)
A(I,J) = ZERO
C---- MAKE "A" MATRIX SYMMETRICAL -------
DO 600 K=1,NEQ
A(I,K) = A(K,I)
A(J,K) = A(K,J)
600 CONTINUE
C---- A(I,J) MADE ZERO BY ROTATION ------
700 CONTINUE
C---- CHECK FOR CONVERGENCE -------------
IF(ABS(SSUM)/SUM .GT.TOL)GO TO 400
DO 900 I=1,NDIR
PS(I)=A(I,I)
900 CONTINUE
RETURN
END
C**************************************************************************** SUBROUTINE ROTTENSOR(TENSOR,DR,NDIR,NSHR)
C ROTATE A TENSOR: TENSOR=DR.TENSOR.trans(DR)
C TENSOR IS STORE
D IN A VETOR(NDIR+NSHR)
C****************************************************************************
c IMPLICIT REAL*8(A-H,O-Z)
DIMENSION TENSOR(NDIR+NSHR),DR(NDIR,NDIR)
DIMENSION STENSOR(3,3),TEMP(3,3),TDR(3,3)
CALL TRANS(DR,TDR)
CALL GETU(TENSOR,STENSOR,NDIR,NSHR)
CALL MUL(DR,STENSOR,TEMP,NDIR)
CALL MUL(TEMP,TDR,STENSOR,NDIR)
C REARRANGE AN
D STOR
E IN A VECTOR
TENSOR(1)=STENSOR(1,1)
TENSOR(2)=STENSOR(2,2)
TENSOR(3)=STENSOR(3,3)
TENSOR(4)=STENSOR(1,2)
TENSOR(5)=STENSOR(2,3)
TENSOR(6)=STENSOR(1,3)
RETURN
END
C**************************************************************************** SUBROUTINE GETDR(FO,FN,DR)
C CALCULATE DELTA_R DF=DR.DU
C DU=TRANSPOSE(DU) DR.TRANSPOSE(DR)=I
C**************************************************************************** c IMPLICIT REAL*8(A-H,O-Z)
DIMENSION FOLD(6),FNEW(6),DR(3,3)
DIMENSION FO(3,3),FN(3,3),FOINV(3,3),DF(3,3),DC(3,3),TRANSDF(3,3) DIMENSION V(3,3),DU(3,3),DUINV(3,3)
CALL GETF(FOLD,FO,3,3)
CALL GETF(FNEW,FN,3,3)
C FOLD^-1
CALL KMINV(FO,FOINV,DET)
C DF=FNEW.FOLD^-1
CALL MUL(FN,FOINV,DF,3)
C TRANSPOSE(DF)
CALL TRANS(DF,TRANSDF)
C CALCULATE TRANSPOSE(DF).DF
CALL MUL(TRANSDF,DF,DC,3)
C CALCULATE THE EIGENV ALUES AN
D EIGENVECTORS OF DC
CALL EIGEN(DC,V,3,1.0e-6)
C CALCULATE DU
DO 20 I=1,3
DO 10 J=1,3
DU(I,J)=SQRT(DC(1,1))*V(I,1)*V(J,1)+
1 SQRT(DC(2,2))*V(I,2)*V(J,2)+
1 SQRT(DC(3,3))*V(I,3)*V(J,3)
10 CONTINUE
20 CONTINUE
C CLACULATE DU^-1
CALL KMINV(DU,DUINV,DET)
C DR=DF.DU^-1
CALL MUL(DF,DUINV,DR,3)
RETURN
END
C**************************************************************************** SUBROUTINE EIGEN(A,V,NEQ,TL)
c IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION A(NEQ,NEQ),V(NEQ,NEQ)
C EIGENV ALUE SOLUTION BY JACOBI METHO
D -
C A --MATRIX (ANY RANK) TO BE SOLVE
D ---
C EIGENV ALUES ON DIAGONAL
C V - MA TRIX OF EIGENVECTORS PRODUCED
C TL- NUMBER OF SIGNIFICANT FIGURES
C---- INITIALIZATION -----------------------
ZERO = 0.0D0
SUM = ZERO
TOL = ABS(TL)
C---- SET INITIAL EIGENVECTORS -------------
DO 200 I=1,NEQ
DO 190 J=1,NEQ
IF (TL.GT.ZERO) V(I,J) = ZERO
190 SUM = SUM + ABS(A(I,J))
IF (TL.GT.ZERO) V(I,I) = 1.0
200 CONTINUE
C---- CHECK FOR TRIVIAL PROBLEM -----------
IF (NEQ.EQ.1) RETURN
IF (SUM.LE.ZERO) RETURN
SUM = SUM/FLOAT(NEQ*NEQ)
C-------------------------------------------
C---- REDUCE MA TRIX TO DIAGONAL ------------
C-------------------------------------------
400 SSUM = ZERO
AMAX = ZERO
DO 700 J=2,NEQ
IH = J -1
DO 700 I=1,IH
C---- CHECK IF A(I,J) IS TO BE REDUCED -----
AA = ABS(A(I,J))
IF (AA.GT.AMAX) AMAX = AA
SSUM = SSUM + AA
IF (AA.LT.0.1*AMAX) GO TO 700
C---- CALCULATE ROTA TION ANGLE ----------
AA=ATAN2(2.0*A(I,J),A(I,I)-A(J,J))/2.0
SI = SIN(AA)
CO = COS(AA)
C---- MODIFY "I" AND "J" COLUMNS --------
DO 500 K=1,NEQ
TT = A(K,I)
A(K,I) = CO*TT + SI*A(K,J)
A(K,J) = -SI*TT + CO*A(K,J)
TT = V(K,I)
V(K,I) = CO*TT + SI*V(K,J)
500 V(K,J) = -SI*TT + CO*V(K,J)
C---- MODIFY DIAGONAL TERMS -------------
A(I,I) = CO*A(I,I) + SI*A(J,I)
A(J,J) =-SI*A(I,J) + CO*A(J,J)
A(I,J) = ZERO
C---- MAKE "A" MATRIX SYMMETRICAL -------
DO 600 K=1,NEQ
A(I,K) = A(K,I)
A(J,K) = A(K,J)
600 CONTINUE
C---- A(I,J) MADE ZERO BY ROTATION ------
700 CONTINUE
C---- CHECK FOR CONVERGENCE -------------
IF(ABS(SSUM)/SUM .GT.TOL)GO TO 400
RETURN
END
C****************************************************************************
C**************************************************************************** SUBROUTINE KMINV(A,AINV,DET_AINV)
C This subroutine calculates the inverse of a {3 x 3} matrix and the
C determinant of the inverse
C**************************************************************************** c IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(3,3), AINV(3,3)
PARAMETER(ZERO=0.D0, ONE=1.D0)
DET_A = A(1,1)*(A(2,2)*A(3,3) - A(3,2)*A(2,3)) -
+ A(2,1)*(A(1,2)*A(3,3) - A(3,2)*A(1,3)) +
+ A(3,1)*(A(1,2)*A(2,3) - A(2,2)*A(1,3))
IF (DET_A .LE. ZERO) THEN
WRITE(80,*) 'WARNING: DET OF MAT IS ZERO/NEGATIVE !!'
ENDIF
DET_AINV = ONE/DET_A
AINV(1,1) = DET_AINV*(A(2,2)*A(3,3) - A(3,2)*A(2,3))
AINV(1,2) = DET_AINV*(A(3,2)*A(1,3) - A(1,2)*A(3,3))
AINV(1,3) = DET_AINV*(A(1,2)*A(2,3) - A(2,2)*A(1,3))
AINV(2,1) = DET_AINV*(A(3,1)*A(2,3) - A(2,1)*A(3,3))
AINV(2,2) = DET_AINV*(A(1,1)*A(3,3) - A(3,1)*A(1,3))
AINV(2,3) = DET_AINV*(A(2,1)*A(1,3) - A(1,1)*A(2,3))
AINV(3,1) = DET_AINV*(A(2,1)*A(3,2) - A(3,1)*A(2,2))
AINV(3,2) = DET_AINV*(A(3,1)*A(1,2) - A(1,1)*A(3,2))
AINV(3,3) = DET_AINV*(A(1,1)*A(2,2) - A(2,1)*A(1,2))
RETURN
END
C******************************************************************* SUBROUTINE MUL(A,B,C,N)
C MULTIPLICATION OF SQUARE MA TRIX A,B. THE RSULTE IS STORE
D IN C C N IS TH
E DIMENSION O
F THE MATRIX
C*******************************************************************
c IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(N,N),B(N,N),C(N,N)
DO 30 I=1,N
DO 20 J=1,N
C(I,J)=0
DO 10 K=1,N
C(I,J)=C(I,J)+A(I,K)*B(K,J)
10 CONTINUE
20 CONTINUE
30 CONTINUE
RETURN
END
C**********************************************************************
SUBROUTINE TRANS(A,A TRANS)
C THIS SUBROUTINE CALCULATES THE TRANSPOSE OF AN 3 BY 3
C MATRIX [A], AN
D PLACES TH
E RESULT IN A TRANS.
C**********************************************************************
c REAL*8 A(3,3),ATRANS(3,3)
dimension A(3,3),A TRANS(3,3)
DO 1 I=1,3
DO 2 J=1,3
ATRANS(J,I) = A(I,J)
2 CONTINUE
1 CONTINUE
RETURN
END
C**********************************************************************
C**************************************************************** SUBROUTINE GETF(F,SF,NDIR,NSHR)
C REARRANGE DEFORMA TION GRADIENT STORE
D IN A ONE-DIMENSION VETOR
C TO A SQUARE MA TRIX
C F(NDIR+2*NSHR) SF(NIDR,NDIR)
C****************************************************************
c IMPLICIT REAL*8(A-H,O-Z)
DIMENSION F(NDIR+2*NSHR),SF(NDIR,NDIR)
IF (NSHR .LT. 3) THEN
WRITE(80,*) 'WARNING: CAN ONL Y BE USED IN 3D MODEL '
ENDIF
SF(1,1)=F(1)
SF(1,2)=F(4)
SF(1,3)=F(9)
SF(2,1)=F(7)
SF(2,2)=F(2)
SF(2,3)=F(5)
SF(3,1)=F(6)
SF(3,2)=F(8)
SF(3,3)=F(3)
RETURN
END
C**************************************************************** SUBROUTINE GETU(U,SU,NDIR,NSHR)
C REARRANGE STRETCH TENSOR STORE
D IN A ONE-DIMENSION VETOR C TO A SQUAR
E MA TRIX
C U(NDIR+NSHR) SU(NIDR,NDIR)
C****************************************************************
c IMPLICIT REAL*8(A-H,O-Z)
DIMENSION U(NDIR+NSHR),SU(NDIR,NDIR)
IF (NSHR .LT. 3) THEN
WRITE(80,*) 'WARNING: CAN ONL Y BE USED IN 3D MODEL '
ENDIF
SU(1,1)=U(1)
SU(1,2)=U(4)
SU(1,3)=U(6)
SU(2,1)=U(4)
SU(2,2)=U(2)
SU(2,3)=U(5)
SU(3,1)=U(6)
SU(3,2)=U(5)
SU(3,3)=U(3)
RETURN
END