CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C USER-ELEMENT for HYDROGEL with coupled deformation and mass diffusion C Jiaping Zhang and Hanqing Jiang C Department of Mechanical and Aerospace Engineering C Arizona State University C Email: hanqing.jiang@asu.edu C Tel: 480-965-1483 C September, 2008 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS, 1 PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME, 2 KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF, 3 LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,PERIOD) C C This is a 8-node brick element. At each note, there are four degrees of freedom: C Ux, Uy, Uz, and C, where C is the chemical potential C Finite deformation is considered. C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C DIMENSION RHS(MLVARX,*),AMATRX(NDOFEL,NDOFEL),PROPS(*), 1 SVARS(NSVARS),ENERGY(8),COORDS(MCRD,NNODE),U(NDOFEL), 2 DU(MLVARX,*),V(NDOFEL),A(NDOFEL),TIME(2),PARAMS(*), 3 JDLTYP(MDLOAD,*),ADLMAG(MDLOAD,*),DDLMAG(MDLOAD,*), 4 PREDEF(2,NPREDF,NNODE),LFLAGS(*),JPROPS(*) C----- The above is the standard input in the USER-ELEMENT in ABAQUS. C NNODE---Number of nodes on the element. This value is defined by the NODES parameter C on the *USER ELEMENT option. C MCRD is the dimension of the *notes and is defined by the COORDINATES parameter on the C *USER ELEMENT option. C NDOFEL---Number of degrees of freedom in the element. PARAMETER (PI=3.14159265358979323846) PARAMETER (IPNML=8) C----- IPNML is the number of Gaussian integration points DIMENSION DIS_T(3,NNODE),DIS_C(3,NNODE),DIS_DELTA(3*NNODE) DIMENSION CON_INC_T(NNODE),CON_INC_C(NNODE),CON_INC_DELTA(NNODE) DIMENSION XYZ(3,NNODE),CON_0(NNODE) DIMENSION POINTS(IPNML,4),XI(3),DSTRAIN(6),DGIT_RATE_T(3,3) DIMENSION STRESS_T(6),STRAIN_T(6),DG_T(3,3),FLUX_T(3) DIMENSION SF(NNODE),DSFDX(NNODE,3),DGIT_C(3,3) DIMENSION BL(6,3*NNODE),BNL(9,3*NNODE),B2(3,NNODE),DG_C(3,3) DIMENSION DGI_T(3,3),DGIT_T(3,3),DGIT_RATE_C(3,3) DIMENSION TN1(3,3*NNODE),TN2(1,NNODE),SMATRIX(9,9) DIMENSION STRAIN_C(6),CON_GRAD_C(3),DDSDDE(6,6),DJDM(3,3) DIMENSION FLUX_C(3),STRESS_C(6) DIMENSION ARRAY1(3*NNODE,6),ARRAY2(3*NNODE,9),ARRAY3(NNODE,3) DIMENSION STLT(3*NNODE,3*NNODE),STL(3*NNODE,3*NNODE) DIMENSION STNLT(3*NNODE,3*NNODE),STNL(3*NNODE,3*NNODE) DIMENSION ZMATRIXT(NNODE,3*NNODE),ZMATRIX(NNODE,3*NNODE) DIMENSION H1MATRIXT(NNODE,NNODE),H1MATRIX(NNODE,NNODE) DIMENSION H3MATRIXT(NNODE,3*NNODE),H3MATRIX(NNODE,3*NNODE) DIMENSION H4MATRIXT(NNODE,3*NNODE),H4MATRIX(NNODE,3*NNODE) DIMENSION P1T(3*NNODE,1),P1(3*NNODE,1) DIMENSION P2T(3*NNODE,1),P2(3*NNODE,1) DIMENSION P3T(3*NNODE,1),P3(3*NNODE,1) DIMENSION P4T(NNODE,1),P4(NNODE,1) DIMENSION P5T(NNODE,1),P5(NNODE,1) DIMENSION P6T(NNODE,1),P6(NNODE,1) DIMENSION P7T(NNODE,1),P7(NNODE,1) DIMENSION P8T(NNODE,1),P8(NNODE,1) DIMENSION P9T(NNODE,1),P9(NNODE,1) DIMENSION DAMPING(4*NNODE,4*NNODE),STIFF(4*NNODE,4*NNODE) DIMENSION BYF(3,1),TRACTION(6,3,1),SINJECTION(6,1,1) DIMENSION EIGEN1(4*NNODE),EIGEN2(4*NNODE) DIMENSION AMATRX_COPY(NDOFEL,NDOFEL) DIMENSION TSTRESSF(9,1),DDSDDEF(9,6) DIMENSION DDSDDMU(6,1),H5MATRIX(3*NNODE,NNODE) DIMENSION H5MATRIXT(3*NNODE,NNODE) C XYZ for initial coordinates at t=0 C CON_0 for initial chemical potential at t=0 C Depending on specific problems, these two matrixes should be adjusted correspondingly DO J=1,NNODE DO I=1,3 c XYZ(I,J)=COORDS(I,J) XYZ(I,J)=COORDS(I,J)/3.2150215080650055063 ENDDO c CON_0(J)=PREDEF(1,1,J) CON_0(J)=0. ENDDO ICOUNT=0 DO INODE=1,NNODE DO IMCRD=1,3 ICOUNT=ICOUNT+1 IPOSN=IMCRD+4*(INODE-1) DISPT=U(IPOSN)-DU(IPOSN,1) DIS_T(IMCRD,INODE)=DISPT+coords(imcrd,inode) & -xyz(imcrd,inode) DIS_C(IMCRD,INODE)=U(IPOSN)+coords(imcrd,inode) & -xyz(imcrd,inode) DIS_DELTA(ICOUNT)=DU(IPOSN,1) ENDDO CON_INC_C(INODE)=U(4*INODE)-CON_0(INODE) CON_INC_T(INODE)=U(4*INODE)-DU(4*INODE,1)-PREDEF(1,1,INODE) CON_INC_DELTA(INODE)=DU(4*INODE,1)+PREDEF(1,1,INODE) & -CON_0(INODE) ENDDO C Read body and surface traction. BYF=0.D0 TRACTION=0.d0 BINJECTION=0.D0 SINJECTION=0.D0 DO I=1,MDLOAD INDEX=JDLTYP(I,1) IF(INDEX.LT.20.AND.INDEX.GT.10) THEN I1=INDEX-10 BYF(I1,1)=ADLMAG(I,1) ENDIF IF(INDEX.GT.200.AND.INDEX.LT.300) THEN I1=(INDEX-200)/10 I2=MOD((INDEX-200),10) TRACTION(I1,I2,1)=ADLMAG(I,1) ENDIF IF(INDEX.EQ.3) BINJECTION=ADLMAG(I,1) IF(INDEX.GT.40.AND.INDEX.LT.50) THEN I1=INDEX-40 SINJECTION(I1,1,1)=ADLMAG(I,1) ENDIF IF(INDEX.LT.-10.AND.INDEX.GT.-20) THEN C Non-uniform distributed body force, using keywords UnNU in *dload, and index=-n C I1=ABS(INDEX)-10 C BYF(I1,1)=ADLMAG(I,1) ENDIF IF(INDEX.GT.-300.AND.INDEX.LT.-200) THEN C Non-uniform distributed surface traction, using keywords UnNU in *dload, and index=-n C I1=(ABS(INDEX)-200)/10 C I2=MOD((ABS(INDEX)-200),10) C TRACTION(I1,I2,1)=ADLMAG(I,1) ENDIF C IF(INDEX.EQ.-3) BINJECTION=ADLMAG(I,1) IF(INDEX.GT.-50.AND.INDEX.LT.-40) THEN C I1=ABS(INDEX)-40 C SINJECTION(I1,1,1)=ADLMAG(I,1) ENDIF ENDDO C C----- Set initial values (zero!) for the stiffness matrix and the right-hand-side term and other matrix DO K1=1,NDOFEL DO KRHS=1,NRHS RHS(K1,KRHS) = 0.0D0 END DO DO K2=1,NDOFEL AMATRX(K2,K1)=0.0D0 END DO END DO STL=0.D0 STNL=0.D0 ZMATRIX=0.D0 H1MATRIX=0.D0 H3MATRIX=0.D0 H4MATRIX=0.D0 H5MATRIX=0.D0 P1=0.D0 P2=0.D0 P3=0.D0 P4=0.D0 P5=0.D0 P6=0.D0 P7=0.D0 P8=0.D0 P9=0.D0 C C----- Gaussian integration points CALL GET_POINTS(IPNML,POINTS) DO IPOINT=1,IPNML C read state variables that were stored at time increment time=t DO I=1,6 STRESS_T(I)=SVARS((IPOINT-1)*32+I) STRAIN_T(I)=SVARS((IPOINT-1)*32+6+I) ENDDO ICOUNT=0 DO I=1,3 DO J=1,3 ICOUNT=ICOUNT+1 DG_T(I,J)=SVARS(ICOUNT+(IPOINT-1)*32+12) ENDDO ENDDO ICOUNT=0 DO I=1,3 DO J=1,3 ICOUNT=ICOUNT+1 DGIT_RATE_T(I,J)=SVARS(ICOUNT+(IPOINT-1)*32+21) ENDDO ENDDO DETDG_T=SVARS((IPOINT-1)*32+31) DETDG_RATE_T=SVARS((IPOINT-1)*32+32) DGI_T(1,1)= (DG_T(2,2)*DG_T(3,3)-DG_T(2,3)*DG_T(3,2))/DETDG_T DGI_T(1,2)=-(DG_T(1,2)*DG_T(3,3)-DG_T(1,3)*DG_T(3,2))/DETDG_T DGI_T(1,3)=-(DG_T(1,3)*DG_T(2,2)-DG_T(1,2)*DG_T(2,3))/DETDG_T DGI_T(2,1)=-(DG_T(2,1)*DG_T(3,3)-DG_T(2,3)*DG_T(3,1))/DETDG_T DGI_T(2,2)= (DG_T(1,1)*DG_T(3,3)-DG_T(1,3)*DG_T(3,1))/DETDG_T DGI_T(2,3)=-(DG_T(1,1)*DG_T(2,3)-DG_T(1,3)*DG_T(2,1))/DETDG_T DGI_T(3,1)=-(DG_T(2,2)*DG_T(3,1)-DG_T(2,1)*DG_T(3,2))/DETDG_T DGI_T(3,2)=-(DG_T(1,1)*DG_T(3,2)-DG_T(1,2)*DG_T(3,1))/DETDG_T DGI_T(3,3)= (DG_T(1,1)*DG_T(2,2)-DG_T(1,2)*DG_T(2,1))/DETDG_T C---- DGIT_T is the transpose of DGI_T DO I=1,3 DO J=1,3 DGIT_T(I,J)=DGI_T(J,I) ENDDO ENDDO XI=POINTS(IPOINT,1:3) WEIGHT=POINTS(IPOINT,4) CALL GRADIENT(NNODE,XYZ,DIS_T,DIS_C,XI,DETJAC,SF,DSFDX, & BL,BNL,B2,DG_C,DGIT_C,DETDG_C,DETDG_T,DGIT_T) DETDG_RATE_C=(DETDG_C-DETDG_T)/DTIME IF(DTIME.LE.1E-10) DETDG_RATE_C=0.D0 DGIT_RATE_C(:,:)=(DGIT_C(:,:)-DGIT_T(:,:))/DTIME IF(DTIME.LE.1E-10) DGIT_RATE_C=0.D0 CALL SHAPEFUNCTION(NNODE,SF,TN1,TN2) TEMPY=0.D0 DO KK=1,NNODE TEMPY=TEMPY+SF(KK)*XYZ(2,KK) ENDDO CON_GAUSS_C=0.D0 DO I=1,NNODE CON_GAUSS_C=CON_GAUSS_C+TN2(1,I)*(CON_INC_C(I)+CON_0(I)) ENDDO CON_INC_GAUSS_C=0.D0 DO I=1,NNODE CON_INC_GAUSS_C=CON_INC_GAUSS_C+TN2(1,I)*CON_INC_DELTA(I) ENDDO CON_RATE_C=CON_INC_GAUSS_C/DTIME IF(DABS(CON_INC_GAUSS_C).LT.1E-10) CON_RATE_C=0.D0 IF(DTIME.LE.1E-10) CON_RATE_C=0.D0 CALL FORM_SMATRIX(SMATRIX,STRESS_T) DSTRAIN=0.D0 DO I=1,6 DO J=1,3*NNODE DSTRAIN(I)=DSTRAIN(I)+BL(I,J)*DIS_DELTA(J) ENDDO ENDDO STRAIN_C(:)=STRAIN_T(:)+DSTRAIN(:) CON_GRAD_C=0.D0 DO I=1,3 DO J=1,NNODE CON_GRAD_C(I)=CON_GRAD_C(I)+B2(I,J) & *(CON_INC_C(J)+CON_0(J)) ENDDO ENDDO CALL DMATRIX(NPROPS,DETDG_C,CON_GRAD_C,CON_GAUSS_C,PROPS,DG_C, & DDSDDE,DJDM,STRESS_C,FLUX_C,CHEMPOTEN, & DDSDDEF,TSTRESSF,DDSDDMU) CALL TTMULTIPLY(BL,DDSDDE,ARRAY1,6,3*NNODE,6) CALL TMULTIPLY(ARRAY1,BL,STLT,3*NNODE,6,3*NNODE) STLT(:,:)=STLT(:,:)*WEIGHT*DETJAC STL=STL+STLT CALL TTMULTIPLY(B2,DJDM,ARRAY3,3,NNODE,3) CALL TMULTIPLY(ARRAY3,B2,H1MATRIXT,NNODE,3,NNODE) H1MATRIXT(:,:)=H1MATRIXT(:,:)*WEIGHT*DETJAC H1MATRIX=H1MATRIX-H1MATRIXT if(lflags(1).eq.71.OR.LFLAGS(1).EQ.1.OR.LFLAGS(1).EQ.2) then ZMATRIX=0.D0 H3MATRIX=0.D0 H4MATRIX=0.D0 P4=0.D0 endif if(lflags(1).eq.72.or.lflags(1).eq.73) then CALL FORM_ZMATRIX(NNODE,DSFDX,TN2,DGIT_T,DETDG_T,ZMATRIXT) ZMATRIXT(:,:)=ZMATRIXT(:,:)*WEIGHT*DETJAC ZMATRIX=ZMATRIX+ZMATRIXT H3MATRIX(:,:)=ZMATRIX(:,:)/DETDG_T*DETDG_RATE_T CALL FORM_H4MATRIX(NNODE,DSFDX,TN2,DETDG_T,DGIT_RATE_T & ,H4MATRIXT) H4MATRIXT(:,:)=H4MATRIXT(:,:)*WEIGHT*DETJAC H4MATRIX=H4MATRIX+H4MATRIXT CALL FORM_H5MATRIX(NNODE,BL,DDSDDMU,TN2,H5MATRIXT) H5MATRIXT(:,:)=H5MATRIXT(:,:)*WEIGHT*DETJAC H5MATRIX=H5MATRIX+H5MATRIXT P4T(:,1)=TN2(1,:)*DETDG_RATE_C*WEIGHT*DETJAC P4=P4+P4T ENDIF CALL TTMULTIPLY(BNL,TSTRESSF,P1T,9,3*NNODE,1) P1T(:,:)=P1T(:,:)*WEIGHT*DETJAC P1=P1+P1T CALL TTMULTIPLY(TN1,BYF,P2T,3,3*NNODE,1) P2T(:,:)=P2T(:,:)*WEIGHT*DETJAC P2=P2+P2T CALL SURFACEINT(IPOINT,NNODE,XI,XYZ,TRACTION,P3T,3,WEIGHT) P3=P3+P3T CALL TTMULTIPLY(B2,FLUX_C,P5T,3,NNODE,1) P5T(:,:)=P5T(:,:)*WEIGHT*DETJAC P5=P5+P5T P6T(:,1)=TN2(1,:)*BINJECTION*WEIGHT*DETJAC P6=P6+P6T P7T(:,1)=TN2(1,:)*(1.D0+CON_GAUSS_C)*WEIGHT*DETJAC P7=P7+P7T P8T(:,1)=TN2(1,:)*DETDG_C*WEIGHT*DETJAC P8=P8+P8T CALL SURFACEINT(IPOINT,NNODE,XI,XYZ,SINJECTION,P9T,1,WEIGHT) P9=P9+P9T DO I=1,6 SVARS((IPOINT-1)*32+I)=STRESS_C(I) SVARS((IPOINT-1)*32+I+6)=STRAIN_C(I) ENDDO ICOUNT=0 DO I=1,3 DO J=1,3 ICOUNT=ICOUNT+1 SVARS((IPOINT-1)*32+ICOUNT+12)=DG_C(I,J) ENDDO ENDDO ICOUNT=0 DO I=1,3 DO J=1,3 ICOUNT=ICOUNT+1 SVARS((IPOINT-1)*32+ICOUNT+21)=DGIT_RATE_C(I,J) ENDDO ENDDO SVARS((IPOINT-1)*32+31)=DETDG_C SVARS((IPOINT-1)*32+32)=DETDG_RATE_C ENDDO DO I=1,NNODE RHS(4*I-3,1)= -P1(3*I-2,1)+P2(3*I-2,1)+P3(3*I-2,1) RHS(4*I-2,1)= -P1(3*I-1,1)+P2(3*I-1,1)+P3(3*I-1,1) RHS(4*I-1,1)= -P1(3*I,1)+P2(3*I,1)+P3(3*I,1) RHS(4*I,1) = -P4(I,1)+P5(I,1)+P6(I,1)-P9(I,1) ENDDO DAMPING=0.D0 DO I=1,4*NNODE DO J=1,4*NNODE IF(MOD(I,4).EQ.0.AND.MOD(J,4).NE.0)THEN K=J/4*3+MOD(J,4) DAMPING(I,J)=0.5*ZMATRIX(I/4,K) ENDIF IF(MOD(I,4).NE.0.AND.MOD(J,4).EQ.0)THEN K=I/4*3+MOD(I,4) DAMPING(I,J)=0.5*ZMATRIX(J/4,K) ENDIF ENDDO ENDDO STIFF=0.D0 DO I=1,4*NNODE DO J=1,4*NNODE IF(MOD(I,4).EQ.0.AND.MOD(J,4).EQ.0) THEN STIFF(I,J) = -H1MATRIX(I/4,J/4) ENDIF IF(MOD(I,4).EQ.0.AND.MOD(J,4).NE.0)THEN K=J/4*3+MOD(J,4) STIFF(I,J)=-H5MATRIX(K,I/4) ENDIF IF(MOD(I,4).NE.0.AND.MOD(J,4).EQ.0)THEN K=I/4*3+MOD(I,4) STIFF(I,J)=H5MATRIX(K,J/4) ENDIF IF(MOD(I,4).NE.0.AND.MOD(J,4).NE.0)THEN K=I/4*3+MOD(I,4) L=J/4*3+MOD(J,4) STIFF(I,J)=STL(K,L) ENDIF ENDDO ENDDO if(lflags(1).eq.71.OR.LFLAGS(1).EQ.1.OR.LFLAGS(1).EQ.2) then AMATRX(:,:)=STIFF(:,:) endif if(lflags(1).eq.72.or.lflags(1).eq.73) then IF(DTIME.LE.1E-10) DTIME=1E-10 AMATRX(:,:)=STIFF(:,:)+DAMPING(:,:)/dtime endif C Modification of the stiffness matrix EPS=1E-6 888 AMATRX_COPY(:,:)=AMATRX(:,:) CALL SYMQR(AMATRX_COPY,EIGEN1,EIGEN2,0.D0,4*NNODE,4*NNODE,EPS, & .FALSE.,.FALSE.,.FALSE.,II) SMALLEST=1000000 DO I=1,4*NNODE IF(EIGEN1(I).LT.SMALLEST) SMALLEST=EIGEN1(I) ENDDO IF(SMALLEST.LT.0.D0) THEN SHIFT=-1.D0*SMALLEST+0.0001 ELSE GOTO 999 ENDIF do i=1,4*nnode amatrx(i,i)=amatrx(i,i)+SHIFT enddo GOTO 888 999 CONTINUE RETURN END C C---- This subroutine calculates DDSDDE, DJDM and stress and flux at time t+delta SUBROUTINE DMATRIX(NPROPS,DETDG,C_GRAD,CONCEN,PROPS,DG, & DDSDDE,DJDM,TSTRESS,FLUX,CHEMPOTEN,DDSDDEF,TSTRESSF & ,DDSDDMU) IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION C_GRAD(3),PROPS(NPROPS),DG(3,3),DDSDDE(6,6),DJDM(3,3) DIMENSION TSTRESS(6),FLUX(3),CC(3,3),CCI(3,3),DDSDDEF(9,6) DIMENSION DGI(3,3),DGIT(3,3),TSTRESSF(9),DDSDDMU(6,1) REAL*8 NV,CHI,DD,INV1,INV3 NV =PROPS(1) CHI=PROPS(2) DD =PROPS(3) C CC is the right Cauchy-Born tensor CC=0.D0 DO I=1,3 DO J=1,3 DO K=1,3 CC(I,J)=CC(I,J)+DG(K,I)*DG(K,J) ENDDO ENDDO ENDDO DETCC=CC(1,1)*CC(2,2)*CC(3,3)-CC(1,1)*CC(2,3)*CC(3,2) & -CC(2,1)*CC(1,2)*CC(3,3)+CC(2,1)*CC(1,3)*CC(3,2) & +CC(3,1)*CC(1,2)*CC(2,3)-CC(3,1)*CC(1,3)*CC(2,2) INV3=DETCC INV1=CC(1,1)+CC(2,2)+CC(3,3) IF(INV3-1.D0.LE.1E-10) INV3=1.D0+1E-10 C---- CCI is the inverse of CC CCI(1,1)= (CC(2,2)*CC(3,3)-CC(2,3)*CC(3,2))/DETCC CCI(1,2)=-(CC(1,2)*CC(3,3)-CC(1,3)*CC(3,2))/DETCC CCI(1,3)=-(CC(1,3)*CC(2,2)-CC(1,2)*CC(2,3))/DETCC CCI(2,1)=-(CC(2,1)*CC(3,3)-CC(2,3)*CC(3,1))/DETCC CCI(2,2)= (CC(1,1)*CC(3,3)-CC(1,3)*CC(3,1))/DETCC CCI(2,3)=-(CC(1,1)*CC(2,3)-CC(1,3)*CC(2,1))/DETCC CCI(3,1)=-(CC(2,2)*CC(3,1)-CC(2,1)*CC(3,2))/DETCC CCI(3,2)=-(CC(1,1)*CC(3,2)-CC(1,2)*CC(3,1))/DETCC CCI(3,3)= (CC(1,1)*CC(2,2)-CC(1,2)*CC(2,1))/DETCC detdg=dsqrt(inv3) DDSDDE=0.D0 if(detdg.lt.1.d0) then ddsdde=0.d0 tstress=0.d0 goto 1000 endif PAR1=-NV-DETDG*DLOG(DETDG/(DETDG-1.D0))+1.D0 & +CHI/DETDG-CONCEN*DETDG PAR2=DETDG*DLOG(DETDG/(DETDG-1.D0))-DETDG/(DETDG-1.D0) & +CHI/DETDG+CONCEN*DETDG DDSDDE(1,1)=PAR1*2*CCI(1,1)*CCI(1,1)-PAR2*CCI(1,1)*CCI(1,1) DDSDDE(1,2)=PAR1*2*CCI(1,2)*CCI(1,2)-PAR2*CCI(1,1)*CCI(2,2) DDSDDE(1,3)=PAR1*2*CCI(1,3)*CCI(1,3)-PAR2*CCI(1,1)*CCI(3,3) DDSDDE(1,4)=PAR1*2*CCI(1,1)*CCI(1,2)-PAR2*CCI(1,1)*CCI(1,2) DDSDDE(1,5)=PAR1*2*CCI(1,1)*CCI(1,3)-PAR2*CCI(1,1)*CCI(1,3) DDSDDE(1,6)=PAR1*2*CCI(1,2)*CCI(1,3)-PAR2*CCI(1,1)*CCI(2,3) DDSDDE(2,1)=DDSDDE(1,2) DDSDDE(2,2)=PAR1*2*CCI(2,2)*CCI(2,2)-PAR2*CCI(2,2)*CCI(2,2) DDSDDE(2,3)=PAR1*2*CCI(2,3)*CCI(2,3)-PAR2*CCI(2,2)*CCI(3,3) DDSDDE(2,4)=PAR1*2*CCI(2,1)*CCI(2,2)-PAR2*CCI(2,2)*CCI(1,2) DDSDDE(2,5)=PAR1*2*CCI(2,1)*CCI(2,3)-PAR2*CCI(2,2)*CCI(1,3) DDSDDE(2,6)=PAR1*2*CCI(2,2)*CCI(2,3)-PAR2*CCI(2,2)*CCI(2,3) DDSDDE(3,1)=DDSDDE(1,3) DDSDDE(3,2)=DDSDDE(2,3) DDSDDE(3,3)=PAR1*2*CCI(3,3)*CCI(3,3)-PAR2*CCI(3,3)*CCI(3,3) DDSDDE(3,4)=PAR1*2*CCI(3,1)*CCI(3,2)-PAR2*CCI(3,3)*CCI(1,2) DDSDDE(3,5)=PAR1*2*CCI(3,1)*CCI(3,3)-PAR2*CCI(3,3)*CCI(1,3) DDSDDE(3,6)=PAR1*2*CCI(3,2)*CCI(3,3)-PAR2*CCI(3,3)*CCI(2,3) DDSDDE(4,1)=DDSDDE(1,4) DDSDDE(4,2)=DDSDDE(2,4) DDSDDE(4,3)=DDSDDE(3,4) DDSDDE(4,4)=PAR1*(CCI(1,1)*CCI(2,2)+CCI(1,2)*CCI(1,2)) & -PAR2*CCI(1,2)*CCI(1,2) DDSDDE(4,5)=PAR1*(CCI(1,1)*CCI(2,3)+CCI(1,3)*CCI(2,1)) & -PAR2*CCI(1,2)*CCI(1,3) DDSDDE(4,6)=PAR1*(CCI(1,2)*CCI(2,3)+CCI(1,3)*CCI(2,2)) & -PAR2*CCI(1,2)*CCI(2,3) DDSDDE(5,1)=DDSDDE(1,5) DDSDDE(5,2)=DDSDDE(2,5) DDSDDE(5,3)=DDSDDE(3,5) DDSDDE(5,4)=DDSDDE(4,5) DDSDDE(5,5)=PAR1*(CCI(1,1)*CCI(3,3)+CCI(1,3)*CCI(1,3)) & -PAR2*CCI(1,3)*CCI(1,3) DDSDDE(5,6)=PAR1*(CCI(1,2)*CCI(3,3)+CCI(1,3)*CCI(3,2)) & -PAR2*CCI(1,3)*CCI(2,3) DDSDDE(6,1)=DDSDDE(1,6) DDSDDE(6,2)=DDSDDE(2,6) DDSDDE(6,3)=DDSDDE(3,6) DDSDDE(6,4)=DDSDDE(4,6) DDSDDE(6,5)=DDSDDE(5,6) DDSDDE(6,6)=PAR1*(CCI(2,2)*CCI(3,3)+CCI(2,3)*CCI(2,3)) & -PAR2*CCI(2,3)*CCI(2,3) DDSDDMU(1,1)=-DETDG*CCI(1,1) DDSDDMU(2,1)=-DETDG*CCI(2,2) DDSDDMU(3,1)=-DETDG*CCI(3,3) DDSDDMU(4,1)=-DETDG*CCI(1,2) DDSDDMU(5,1)=-DETDG*CCI(1,3) DDSDDMU(6,1)=-DETDG*CCI(2,3) DDSDDEF=0.D0 DO J=1,6 DDSDDEF(1,J)=DG(1,1)*DDSDDE(1,J)+DG(1,2)*DDSDDE(4,J) & +DG(1,3)*DDSDDE(5,J) DDSDDEF(2,J)=DG(1,1)*DDSDDE(4,J)+DG(1,2)*DDSDDE(2,J) & +DG(1,3)*DDSDDE(6,J) DDSDDEF(3,J)=DG(1,1)*DDSDDE(5,J)+DG(1,2)*DDSDDE(6,J) & +DG(1,3)*DDSDDE(3,J) DDSDDEF(4,J)=DG(2,1)*DDSDDE(1,J)+DG(2,2)*DDSDDE(4,J) & +DG(2,3)*DDSDDE(5,J) DDSDDEF(5,J)=DG(2,1)*DDSDDE(4,J)+DG(2,2)*DDSDDE(2,J) & +DG(2,3)*DDSDDE(6,J) DDSDDEF(6,J)=DG(2,1)*DDSDDE(5,J)+DG(2,2)*DDSDDE(6,J) & +DG(2,3)*DDSDDE(3,J) DDSDDEF(7,J)=DG(3,1)*DDSDDE(1,J)+DG(3,2)*DDSDDE(4,J) & +DG(3,3)*DDSDDE(5,J) DDSDDEF(8,J)=DG(3,1)*DDSDDE(4,J)+DG(3,2)*DDSDDE(2,J) & +DG(3,3)*DDSDDE(6,J) DDSDDEF(9,J)=DG(3,1)*DDSDDE(5,J)+DG(3,2)*DDSDDE(6,J) & +DG(3,3)*DDSDDE(3,J) ENDDO C---- TSTRESS is the stress at t+delta PAR3=DETDG*DLOG(DETDG/(DETDG-1.D0))-1.D0 & -CHI/DETDG+CONCEN*DETDG TSTRESS(1)= NV*(1.D0-CCI(1,1))-PAR3*CCI(1,1) TSTRESS(2)= NV*(1.D0-CCI(2,2))-PAR3*CCI(2,2) TSTRESS(3)= NV*(1.D0-CCI(3,3))-PAR3*CCI(3,3) TSTRESS(4)=-NV*CCI(1,2)-PAR3*CCI(1,2) TSTRESS(5)=-NV*CCI(1,3)-PAR3*CCI(1,3) TSTRESS(6)=-NV*CCI(2,3)-PAR3*CCI(2,3) C TSTRESSF(1)=DG(1,1)*TSTRESS(1)+DG(1,2)*TSTRESS(4) & +DG(1,3)*TSTRESS(5) TSTRESSF(2)=DG(1,1)*TSTRESS(4)+DG(1,2)*TSTRESS(2) & +DG(1,3)*TSTRESS(6) TSTRESSF(3)=DG(1,1)*TSTRESS(5)+DG(1,2)*TSTRESS(6) & +DG(1,3)*TSTRESS(3) TSTRESSF(4)=DG(2,1)*TSTRESS(1)+DG(2,2)*TSTRESS(4) & +DG(2,3)*TSTRESS(5) TSTRESSF(5)=DG(2,1)*TSTRESS(4)+DG(2,2)*TSTRESS(2) & +DG(2,3)*TSTRESS(6) TSTRESSF(6)=DG(2,1)*TSTRESS(5)+DG(2,2)*TSTRESS(6) & +DG(2,3)*TSTRESS(3) TSTRESSF(7)=DG(3,1)*TSTRESS(1)+DG(3,2)*TSTRESS(4) & +DG(3,3)*TSTRESS(5) TSTRESSF(8)=DG(3,1)*TSTRESS(4)+DG(3,2)*TSTRESS(2) & +DG(3,3)*TSTRESS(6) TSTRESSF(9)=DG(3,1)*TSTRESS(5)+DG(3,2)*TSTRESS(6) & +DG(3,3)*TSTRESS(3) C---- DGI is the inverse of DG 1000 DGI(1,1)= (DG(2,2)*DG(3,3)-DG(2,3)*DG(3,2))/DETDG DGI(1,2)=-(DG(1,2)*DG(3,3)-DG(1,3)*DG(3,2))/DETDG DGI(1,3)=-(DG(1,3)*DG(2,2)-DG(1,2)*DG(2,3))/DETDG DGI(2,1)=-(DG(2,1)*DG(3,3)-DG(2,3)*DG(3,1))/DETDG DGI(2,2)= (DG(1,1)*DG(3,3)-DG(1,3)*DG(3,1))/DETDG DGI(2,3)=-(DG(1,1)*DG(2,3)-DG(1,3)*DG(2,1))/DETDG DGI(3,1)=-(DG(2,2)*DG(3,1)-DG(2,1)*DG(3,2))/DETDG DGI(3,2)=-(DG(1,1)*DG(3,2)-DG(1,2)*DG(3,1))/DETDG DGI(3,3)= (DG(1,1)*DG(2,2)-DG(1,2)*DG(2,1))/DETDG C---- DGIT is the transpose of DGI DO I=1,3 DO J=1,3 DGIT(I,J)=DGI(J,I) ENDDO ENDDO C---- DJDM is the mobility tensor M DJDM=0.D0 DO I=1,3 DO J=1,3 DO K=1,3 DJDM(I,J)=DJDM(I,J)+DD*(DETDG-1.D0) & *DGIT(K,I)*DGIT(K,J) ENDDO ENDDO ENDDO FLUX=0.D0 DO I=1,3 DO J=1,3 FLUX(I)=FLUX(I)-DJDM(I,J)*C_GRAD(J) ENDDO ENDDO RETURN END C C C---- This subroutine gets the Gaussian integration points SUBROUTINE GET_POINTS(IP,POINTS) IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION POINTS(IP,4) POINTS(1,1)= 1.D0/DSQRT(3.D0) POINTS(1,2)=-1.D0/DSQRT(3.D0) POINTS(1,3)=-1.D0/DSQRT(3.D0) POINTS(1,4)= 1.D0 POINTS(2,1)= 1.D0/DSQRT(3.D0) POINTS(2,2)= 1.D0/DSQRT(3.D0) POINTS(2,3)=-1.D0/DSQRT(3.D0) POINTS(2,4)= 1.D0 POINTS(3,1)=-1.D0/DSQRT(3.D0) POINTS(3,2)= 1.D0/DSQRT(3.D0) POINTS(3,3)=-1.D0/DSQRT(3.D0) POINTS(3,4)= 1.D0 POINTS(4,1)=-1.D0/DSQRT(3.D0) POINTS(4,2)=-1.D0/DSQRT(3.D0) POINTS(4,3)=-1.D0/DSQRT(3.D0) POINTS(4,4)= 1.D0 POINTS(5,1)= 1.D0/DSQRT(3.D0) POINTS(5,2)=-1.D0/DSQRT(3.D0) POINTS(5,3)= 1.D0/DSQRT(3.D0) POINTS(5,4)= 1.D0 POINTS(6,1)= 1.D0/DSQRT(3.D0) POINTS(6,2)= 1.D0/DSQRT(3.D0) POINTS(6,3)= 1.D0/DSQRT(3.D0) POINTS(6,4)= 1.D0 POINTS(7,1)=-1.D0/DSQRT(3.D0) POINTS(7,2)= 1.D0/DSQRT(3.D0) POINTS(7,3)= 1.D0/DSQRT(3.D0) POINTS(7,4)= 1.D0 POINTS(8,1)=-1.D0/DSQRT(3.D0) POINTS(8,2)=-1.D0/DSQRT(3.D0) POINTS(8,3)= 1.D0/DSQRT(3.D0) POINTS(8,4)= 1.D0 RETURN END C C C C This subroutine calculates all gradient operations, such as Jacobian, deformation C gradient and B matrix SUBROUTINE GRADIENT(NNODE,XYZ,DIS_T,DIS_C,XI,DETJAC,SF,DSFDX, & BL,BNL,B2,DG,DGIT,DETDG,DETDG_T,DGIT_T) IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION XYZ(3,NNODE),DIS_T(3,NNODE),DIS_C(3,NNODE),XI(3) DIMENSION SF(NNODE),DSFDX(NNODE,3),BL(6,3*NNODE),BNL(9,3*NNODE) DIMENSION B2(3,NNODE),DG(3,3),DGI(3,3),DGIT(3,3) DIMENSION DSFDXI(NNODE,3),FJAC(3,3),FJACI(3,3) DIMENSION BL0(6,3*NNODE),BL1(6,3*NNODE),SL(3,3),CC(3,3) DIMENSION DG_T(3,3),DGIT_T(3,3),DGI_T(3,3) C---- XYZ is the INITIAL positions of nodes, i.e., COORS = X C---- SF is the shape function SF(1)=(1.D0+XI(1))*(1.D0-XI(2))*(1.D0-XI(3))*0.125 SF(2)=(1.D0+XI(1))*(1.D0+XI(2))*(1.D0-XI(3))*0.125 SF(3)=(1.D0-XI(1))*(1.D0+XI(2))*(1.D0-XI(3))*0.125 SF(4)=(1.D0-XI(1))*(1.D0-XI(2))*(1.D0-XI(3))*0.125 SF(5)=(1.D0+XI(1))*(1.D0-XI(2))*(1.D0+XI(3))*0.125 SF(6)=(1.D0+XI(1))*(1.D0+XI(2))*(1.D0+XI(3))*0.125 SF(7)=(1.D0-XI(1))*(1.D0+XI(2))*(1.D0+XI(3))*0.125 SF(8)=(1.D0-XI(1))*(1.D0-XI(2))*(1.D0+XI(3))*0.125 C---- DSFDXI = Partial SF_I/Partial XI_J; SF=1..NNODE, XI=1..3 DSFDXI(1,1)= (1.D0-XI(2))*(1.D0-XI(3))*0.125 DSFDXI(1,2)=-(1.D0+XI(1))*(1.D0-XI(3))*0.125 DSFDXI(1,3)=-(1.D0+XI(1))*(1.D0-XI(2))*0.125 DSFDXI(2,1)= (1.D0+XI(2))*(1.D0-XI(3))*0.125 DSFDXI(2,2)= (1.D0+XI(1))*(1.D0-XI(3))*0.125 DSFDXI(2,3)=-(1.D0+XI(1))*(1.D0+XI(2))*0.125 DSFDXI(3,1)=-(1.D0+XI(2))*(1.D0-XI(3))*0.125 DSFDXI(3,2)= (1.D0-XI(1))*(1.D0-XI(3))*0.125 DSFDXI(3,3)=-(1.D0-XI(1))*(1.D0+XI(2))*0.125 DSFDXI(4,1)=-(1.D0-XI(2))*(1.D0-XI(3))*0.125 DSFDXI(4,2)=-(1.D0-XI(1))*(1.D0-XI(3))*0.125 DSFDXI(4,3)=-(1.D0-XI(1))*(1.D0-XI(2))*0.125 DSFDXI(5,1)= (1.D0-XI(2))*(1.D0+XI(3))*0.125 DSFDXI(5,2)=-(1.D0+XI(1))*(1.D0+XI(3))*0.125 DSFDXI(5,3)= (1.D0+XI(1))*(1.D0-XI(2))*0.125 DSFDXI(6,1)= (1.D0+XI(2))*(1.D0+XI(3))*0.125 DSFDXI(6,2)= (1.D0+XI(1))*(1.D0+XI(3))*0.125 DSFDXI(6,3)= (1.D0+XI(1))*(1.D0+XI(2))*0.125 DSFDXI(7,1)=-(1.D0+XI(2))*(1.D0+XI(3))*0.125 DSFDXI(7,2)= (1.D0-XI(1))*(1.D0+XI(3))*0.125 DSFDXI(7,3)= (1.D0-XI(1))*(1.D0+XI(2))*0.125 DSFDXI(8,1)=-(1.D0-XI(2))*(1.D0+XI(3))*0.125 DSFDXI(8,2)=-(1.D0-XI(1))*(1.D0+XI(3))*0.125 DSFDXI(8,3)= (1.D0-XI(1))*(1.D0-XI(2))*0.125 C DEFINE JACOBIAN MATRIX FJAC=0.0D0 DO I=1,3 DO J=1,3 DO K=1,NNODE FJAC(I,J)=FJAC(I,J)+DSFDXI(K,I)*XYZ(J,K) ENDDO ENDDO ENDDO C DEFINE DETERMINANT OF JACOB DETJAC=FJAC(1,1)*FJAC(2,2)*FJAC(3,3)-FJAC(1,1)*FJAC(2,3)*FJAC(3,2) & -FJAC(2,1)*FJAC(1,2)*FJAC(3,3)+FJAC(2,1)*FJAC(1,3)*FJAC(3,2) & +FJAC(3,1)*FJAC(1,2)*FJAC(2,3)-FJAC(3,1)*FJAC(1,3)*FJAC(2,2) C---- FJACI is the inverse of FJAC FJACI(1,1)= (FJAC(2,2)*FJAC(3,3)-FJAC(2,3)*FJAC(3,2))/DETJAC FJACI(1,2)=-(FJAC(1,2)*FJAC(3,3)-FJAC(1,3)*FJAC(3,2))/DETJAC FJACI(1,3)=-(FJAC(1,3)*FJAC(2,2)-FJAC(1,2)*FJAC(2,3))/DETJAC FJACI(2,1)=-(FJAC(2,1)*FJAC(3,3)-FJAC(2,3)*FJAC(3,1))/DETJAC FJACI(2,2)= (FJAC(1,1)*FJAC(3,3)-FJAC(1,3)*FJAC(3,1))/DETJAC FJACI(2,3)=-(FJAC(1,1)*FJAC(2,3)-FJAC(1,3)*FJAC(2,1))/DETJAC FJACI(3,1)=-(FJAC(2,2)*FJAC(3,1)-FJAC(2,1)*FJAC(3,2))/DETJAC FJACI(3,2)=-(FJAC(1,1)*FJAC(3,2)-FJAC(1,2)*FJAC(3,1))/DETJAC FJACI(3,3)= (FJAC(1,1)*FJAC(2,2)-FJAC(1,2)*FJAC(2,1))/DETJAC C GET DSFDX:PATRIAL SF_I/PARTIAL X_J DSFDX=0.D0 DO I=1,NNODE DO J=1,3 DO K=1,3 DSFDX(I,J)=DSFDX(I,J)+FJACI(J,K)*DSFDXI(I,K) ENDDO ENDDO ENDDO C---- BL0 is the linear geometrical matrix. C---- The corresponding strains are e11, e22, e33, 2e12, 2e13, 2e23. BL0=0.D0 DO I=1,NNODE BL0(1,3*I-2)=DSFDX(I,1) BL0(2,3*I-1)=DSFDX(I,2) BL0(3,3*I) =DSFDX(I,3) BL0(4,3*I-2)=DSFDX(I,2) BL0(4,3*I-1)=DSFDX(I,1) BL0(5,3*I-2)=DSFDX(I,3) BL0(5,3*I) =DSFDX(I,1) BL0(6,3*I-1)=DSFDX(I,3) BL0(6,3*I) =DSFDX(I,2) ENDDO C DG --- deformation gradient at time t+delta DG=0.D0 DO I=1,NNODE DG(1,1)=DG(1,1)+BL0(1,3*I-2)*DIS_C(1,I) DG(1,2)=DG(1,2)+BL0(4,3*I-2)*DIS_C(1,I) DG(1,3)=DG(1,3)+BL0(5,3*I-2)*DIS_C(1,I) DG(2,1)=DG(2,1)+BL0(4,3*I-1)*DIS_C(2,I) DG(2,2)=DG(2,2)+BL0(2,3*I-1)*DIS_C(2,I) DG(2,3)=DG(2,3)+BL0(6,3*I-1)*DIS_C(2,I) DG(3,1)=DG(3,1)+BL0(5,3*I) *DIS_C(3,I) DG(3,2)=DG(3,2)+BL0(6,3*I) *DIS_C(3,I) DG(3,3)=DG(3,3)+BL0(3,3*I) *DIS_C(3,I) ENDDO DG(1,1)=1.D0+DG(1,1) DG(2,2)=1.D0+DG(2,2) DG(3,3)=1.D0+DG(3,3) CC=0.D0 DO I=1,3 DO J=1,3 DO K=1,3 CC(I,J)=CC(I,J)+DG(K,I)*DG(K,J) ENDDO ENDDO ENDDO DETCC=CC(1,1)*CC(2,2)*CC(3,3)-CC(1,1)*CC(2,3)*CC(3,2) & -CC(2,1)*CC(1,2)*CC(3,3)+CC(2,1)*CC(1,3)*CC(3,2) & +CC(3,1)*CC(1,2)*CC(2,3)-CC(3,1)*CC(1,3)*CC(2,2) DETDG=DSQRT(DETCC) C---- DGI is the inverse of DG DGI(1,1)= (DG(2,2)*DG(3,3)-DG(2,3)*DG(3,2))/DETDG DGI(1,2)=-(DG(1,2)*DG(3,3)-DG(1,3)*DG(3,2))/DETDG DGI(1,3)=-(DG(1,3)*DG(2,2)-DG(1,2)*DG(2,3))/DETDG DGI(2,1)=-(DG(2,1)*DG(3,3)-DG(2,3)*DG(3,1))/DETDG DGI(2,2)= (DG(1,1)*DG(3,3)-DG(1,3)*DG(3,1))/DETDG DGI(2,3)=-(DG(1,1)*DG(2,3)-DG(1,3)*DG(2,1))/DETDG DGI(3,1)=-(DG(2,2)*DG(3,1)-DG(2,1)*DG(3,2))/DETDG DGI(3,2)=-(DG(1,1)*DG(3,2)-DG(1,2)*DG(3,1))/DETDG DGI(3,3)= (DG(1,1)*DG(2,2)-DG(1,2)*DG(2,1))/DETDG C---- DGIT is the transpose of DGI DO I=1,3 DO J=1,3 DGIT(I,J)=DGI(J,I) ENDDO ENDDO DG_T=0.D0 DO I=1,NNODE DG_T(1,1)=DG_T(1,1)+BL0(1,3*I-2)*DIS_T(1,I) DG_T(1,2)=DG_T(1,2)+BL0(4,3*I-2)*DIS_T(1,I) DG_T(1,3)=DG_T(1,3)+BL0(5,3*I-2)*DIS_T(1,I) DG_T(2,1)=DG_T(2,1)+BL0(4,3*I-1)*DIS_T(2,I) DG_T(2,2)=DG_T(2,2)+BL0(2,3*I-1)*DIS_T(2,I) DG_T(2,3)=DG_T(2,3)+BL0(6,3*I-1)*DIS_T(2,I) DG_T(3,1)=DG_T(3,1)+BL0(5,3*I) *DIS_T(3,I) DG_T(3,2)=DG_T(3,2)+BL0(6,3*I) *DIS_T(3,I) DG_T(3,3)=DG_T(3,3)+BL0(3,3*I) *DIS_T(3,I) ENDDO DG_T(1,1)=1.D0+DG_T(1,1) DG_T(2,2)=1.D0+DG_T(2,2) DG_T(3,3)=1.D0+DG_T(3,3) DETDG_T=DG_T(1,1)*DG_T(2,2)*DG_T(3,3)-DG_T(1,1)*DG_T(2,3) & *DG_T(3,2)-DG_T(2,1)*DG_T(1,2)*DG_T(3,3)+DG_T(2,1)*DG_T(1,3) & *DG_T(3,2)+DG_T(3,1)*DG_T(1,2)*DG_T(2,3)-DG_T(3,1)*DG_T(1,3) & *DG_T(2,2) C---- DGI is the inverse of DG DGI_T(1,1)= (DG_T(2,2)*DG_T(3,3)-DG_T(2,3)*DG_T(3,2))/DETDG_T DGI_T(1,2)=-(DG_T(1,2)*DG_T(3,3)-DG_T(1,3)*DG_T(3,2))/DETDG_T DGI_T(1,3)=-(DG_T(1,3)*DG_T(2,2)-DG_T(1,2)*DG_T(2,3))/DETDG_T DGI_T(2,1)=-(DG_T(2,1)*DG_T(3,3)-DG_T(2,3)*DG_T(3,1))/DETDG_T DGI_T(2,2)= (DG_T(1,1)*DG_T(3,3)-DG_T(1,3)*DG_T(3,1))/DETDG_T DGI_T(2,3)=-(DG_T(1,1)*DG_T(2,3)-DG_T(1,3)*DG_T(2,1))/DETDG_T DGI_T(3,1)=-(DG_T(2,2)*DG_T(3,1)-DG_T(2,1)*DG_T(3,2))/DETDG_T DGI_T(3,2)=-(DG_T(1,1)*DG_T(3,2)-DG_T(1,2)*DG_T(3,1))/DETDG_T DGI_T(3,3)= (DG_T(1,1)*DG_T(2,2)-DG_T(1,2)*DG_T(2,1))/DETDG_T C---- DGIT is the transpose of DGI DO I=1,3 DO J=1,3 DGIT_T(I,J)=DGI_T(J,I) ENDDO ENDDO C SL is the partial (0_t U) / partial (X) SL=0.D0 DO I=1,3 DO J=1,3 DO INODE=1,NNODE SL(I,J)=SL(I,J)+DIS_T(I,INODE)*DSFDX(INODE,J) ENDDO ENDDO ENDDO C---- BL1 is the linear geometrical matrix. C---- The corresponding strains are e11, e22, e33, 2e12, 2e13, 2e23. BL1=0.D0 DO I=1,NNODE BL1(1,3*I-2)=SL(1,1)*DSFDX(I,1) BL1(1,3*I-1)=SL(2,1)*DSFDX(I,1) BL1(1,3*I) =SL(3,1)*DSFDX(I,1) BL1(2,3*I-2)=SL(1,2)*DSFDX(I,2) BL1(2,3*I-1)=SL(2,2)*DSFDX(I,2) BL1(2,3*I) =SL(3,2)*DSFDX(I,2) BL1(3,3*I-2)=SL(1,3)*DSFDX(I,3) BL1(3,3*I-1)=SL(2,3)*DSFDX(I,3) BL1(3,3*I) =SL(3,3)*DSFDX(I,3) BL1(4,3*I-2)=SL(1,1)*DSFDX(I,2)+SL(1,2)*DSFDX(I,1) BL1(4,3*I-1)=SL(2,1)*DSFDX(I,2)+SL(2,2)*DSFDX(I,1) BL1(4,3*I) =SL(3,1)*DSFDX(I,2)+SL(3,2)*DSFDX(I,1) BL1(5,3*I-2)=SL(1,1)*DSFDX(I,3)+SL(1,3)*DSFDX(I,1) BL1(5,3*I-1)=SL(2,1)*DSFDX(I,3)+SL(2,3)*DSFDX(I,1) BL1(5,3*I) =SL(3,1)*DSFDX(I,3)+SL(3,3)*DSFDX(I,1) BL1(6,3*I-2)=SL(1,2)*DSFDX(I,3)+SL(1,3)*DSFDX(I,2) BL1(6,3*I-1)=SL(2,2)*DSFDX(I,3)+SL(2,3)*DSFDX(I,2) BL1(6,3*I) =SL(3,2)*DSFDX(I,3)+SL(3,3)*DSFDX(I,2) ENDDO BL=BL0+BL1 C BNL is the non-linear geometrical matrix BNL=0.D0 DO I=1,NNODE BNL(1,3*I-2) = DSFDX(I,1) BNL(2,3*I-2) = DSFDX(I,2) BNL(3,3*I-2) = DSFDX(I,3) BNL(4,3*I-1) = DSFDX(I,1) BNL(5,3*I-1) = DSFDX(I,2) BNL(6,3*I-1) = DSFDX(I,3) BNL(7,3*I) = DSFDX(I,1) BNL(8,3*I) = DSFDX(I,2) BNL(9,3*I) = DSFDX(I,3) ENDDO C C B2 is the geometrical function for nable C DO I=1,NNODE DO J=1,3 B2(J,I)=DSFDX(I,J) ENDDO ENDDO RETURN END C C C C---- This subroutine calculates the shape function SUBROUTINE SHAPEFUNCTION(NNODE,SF,TN1,TN2) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SF(NNODE),TN1(3,3*NNODE),TN2(1,NNODE) TN1=0.D0 TN2=0.D0 DO I=1,NNODE TN1(1,3*I-2)=SF(I) TN1(2,3*I-1)=SF(I) TN1(3,3*I) =SF(I) TN2(1,I) =SF(I) ENDDO RETURN END C C C C---- This subrontine calculates the stress matrix for finite deformation SUBROUTINE FORM_SMATRIX(SMATRIX,STRESS_T) IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION SMATRIX(9,9),STRESS_T(6) C SMATRIX is the second PK stress matrix associated with non-linear geometrical matrix BNL C STRESS_T is the stress vector at time t. The components are T11, T22, T33, T12, T13, T23 SMATRIX = 0.D0 SMATRIX(1,1) = STRESS_T(1) SMATRIX(1,2) = STRESS_T(4) SMATRIX(1,3) = STRESS_T(5) SMATRIX(2,1) = STRESS_T(4) SMATRIX(2,2) = STRESS_T(2) SMATRIX(2,3) = STRESS_T(6) SMATRIX(3,1) = STRESS_T(5) SMATRIX(3,2) = STRESS_T(6) SMATRIX(3,3) = STRESS_T(3) DO I=1,3 DO J=1,3 SMATRIX(I+3,J+3)=SMATRIX(I,J) SMATRIX(I+6,J+6)=SMATRIX(I,J) ENDDO ENDDO RETURN END C C C C SUBROUTINE TMULTIPLY(A,B,C,L,M,N) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(L,M),B(M,N),C(L,N) C C=A*B C=0.D0 DO I=1,L DO J=1,N DO K=1,M C(I,J)=C(I,J)+A(I,K)*B(K,J) ENDDO ENDDO ENDDO RETURN END SUBROUTINE TTMULTIPLY(A,B,C,L,M,N) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(L,M),B(L,N),C(M,N) C C=Transpose(A)*B C=0.D0 DO I=1,M DO J=1,N DO K=1,L C(I,J)=C(I,J)+A(K,I)*B(K,J) ENDDO ENDDO ENDDO RETURN END SUBROUTINE FORM_H4MATRIX(NNODE,DSFDX,TN2,DETDG_T,DGIT_RATE_T & ,H4MATRIX) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION DSFDX(NNODE,3),TN2(1,NNODE),A(1,3*NNODE) DIMENSION H4MATRIX(NNODE,3*NNODE),DGIT_RATE_T(3,3) DO I=1,NNODE A(1,3*I-2) = DGIT_RATE_T(1,1)*DSFDX(I,1)+DGIT_RATE_T(1,2) & *DSFDX(I,2) +DGIT_RATE_T(1,3)*DSFDX(I,3) A(1,3*I-1) = DGIT_RATE_T(2,1)*DSFDX(I,1)+DGIT_RATE_T(2,2) & *DSFDX(I,2) +DGIT_RATE_T(2,3)*DSFDX(I,3) A(1,3*I) = DGIT_RATE_T(3,1)*DSFDX(I,1)+DGIT_RATE_T(3,2) & *DSFDX(I,2) +DGIT_RATE_T(3,3)*DSFDX(I,3) ENDDO CALL TTMULTIPLY(TN2,A,H4MATRIX,1,NNODE,3*NNODE) H4MATRIX(:,:)=H4MATRIX(:,:)*DETDG_T RETURN END SUBROUTINE FORM_ZMATRIX(NNODE,DSFDX,TN2,DGIT_T,DETDG_T,ZMATRIX) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION DSFDX(NNODE,3),TN2(1,NNODE),A(1,3*NNODE) DIMENSION ZMATRIX(NNODE,3*NNODE),DGIT_T(3,3) DO I=1,NNODE A(1,3*I-2) = DGIT_T(1,1)*DSFDX(I,1)+DGIT_T(1,2)*DSFDX(I,2) & +DGIT_T(1,3)*DSFDX(I,3) A(1,3*I-1) = DGIT_T(2,1)*DSFDX(I,1)+DGIT_T(2,2)*DSFDX(I,2) & +DGIT_T(2,3)*DSFDX(I,3) A(1,3*I) = DGIT_T(3,1)*DSFDX(I,1)+DGIT_T(3,2)*DSFDX(I,2) & +DGIT_T(3,3)*DSFDX(I,3) ENDDO CALL TTMULTIPLY(TN2,A,ZMATRIX,1,NNODE,3*NNODE) ZMATRIX(:,:)=ZMATRIX(:,:)*DETDG_T RETURN END SUBROUTINE SURFACEINT(IPOINT,NNODE,XI,COORDS,VECTOR, & FORCE,NDIM,WEIGHT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION XI(3),COORDS(3,NNODE),VECTOR(6,NDIM,1),AA(3) DIMENSION FORCE(NDIM*NNODE,1),TN1S(NDIM,NDIM*NNODE) IF(IPOINT.EQ.1) THEN NSURFACE1=1 NSURFACE2=3 NSURFACE3=6 ENDIF IF(IPOINT.EQ.2) THEN NSURFACE1=1 NSURFACE2=3 NSURFACE3=4 ENDIF IF(IPOINT.EQ.3) THEN NSURFACE1=1 NSURFACE2=4 NSURFACE3=5 ENDIF IF(IPOINT.EQ.4) THEN NSURFACE1=1 NSURFACE2=5 NSURFACE3=6 ENDIF IF(IPOINT.EQ.5) THEN NSURFACE1=2 NSURFACE2=3 NSURFACE3=6 ENDIF IF(IPOINT.EQ.6) THEN NSURFACE1=2 NSURFACE2=3 NSURFACE3=4 ENDIF IF(IPOINT.EQ.7) THEN NSURFACE1=2 NSURFACE2=4 NSURFACE3=5 ENDIF IF(IPOINT.EQ.8) THEN NSURFACE1=2 NSURFACE2=5 NSURFACE3=6 ENDIF FORCE=0.D0 CALL SURFACESHAPEFUNCTION(NNODE,NDIM,XI,NSURFACE1,TN1S,AA,COORDS) DO I=1,3*NNODE DO K=1,NDIM FORCE(I,1)=FORCE(I,1)+VECTOR(NSURFACE1,K,1)*TN1S(K,I) & *WEIGHT*AA(1) ENDDO ENDDO CALL SURFACESHAPEFUNCTION(NNODE,NDIM,XI,NSURFACE2,TN1S,AA,COORDS) DO I=1,3*NNODE DO K=1,NDIM FORCE(I,1)=FORCE(I,1)+VECTOR(NSURFACE2,K,1)*TN1S(K,I) & *WEIGHT*AA(2) ENDDO ENDDO CALL SURFACESHAPEFUNCTION(NNODE,NDIM,XI,NSURFACE3,TN1S,AA,COORDS) DO I=1,3*NNODE DO K=1,NDIM FORCE(I,1)=FORCE(I,1)+VECTOR(NSURFACE3,K,1)*TN1S(K,I) & *WEIGHT*AA(3) ENDDO ENDDO RETURN END SUBROUTINE SURFACESHAPEFUNCTION(NNODE,NDIM,XI,NSURFACE,TN1S,AA, & COORDS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SF(NNODE),TN1S(NDIM,NDIM*NNODE),XI(3),XYZ(3) DIMENSION FJAC(3,3),AA(3),COORDS(3,NNODE),DSFDXI(NNODE,3) XYZ=XI IF(NSURFACE.EQ.1) XYZ(3)=-1.D0 IF(NSURFACE.EQ.2) XYZ(3)= 1.D0 IF(NSURFACE.EQ.3) XYZ(1)= 1.D0 IF(NSURFACE.EQ.4) XYZ(2)= 1.D0 IF(NSURFACE.EQ.5) XYZ(1)=-1.D0 IF(NSURFACE.EQ.6) XYZ(2)=-1.D0 SF(1)=(1.D0+XYZ(1))*(1.D0-XYZ(2))*(1.D0-XYZ(3))*0.125 SF(2)=(1.D0+XYZ(1))*(1.D0+XYZ(2))*(1.D0-XYZ(3))*0.125 SF(3)=(1.D0-XYZ(1))*(1.D0+XYZ(2))*(1.D0-XYZ(3))*0.125 SF(4)=(1.D0-XYZ(1))*(1.D0-XYZ(2))*(1.D0-XYZ(3))*0.125 SF(5)=(1.D0+XYZ(1))*(1.D0-XYZ(2))*(1.D0+XYZ(3))*0.125 SF(6)=(1.D0+XYZ(1))*(1.D0+XYZ(2))*(1.D0+XYZ(3))*0.125 SF(7)=(1.D0-XYZ(1))*(1.D0+XYZ(2))*(1.D0+XYZ(3))*0.125 SF(8)=(1.D0-XYZ(1))*(1.D0-XYZ(2))*(1.D0+XYZ(3))*0.125 TN1S=0.D0 IF(NDIM.EQ.3) THEN DO I=1,NNODE TN1S(1,3*I-2)=SF(I) TN1S(2,3*I-1)=SF(I) TN1S(3,3*I) =SF(I) ENDDO ELSE DO I=1,NNODE TN1S(1,I)=SF(I) ENDDO ENDIF C DEFINE PARTIAL SHAPEFUNCTION TO XI AND ELTA; DSFDXI(1,1)= (1.D0-XYZ(2))*(1.D0-XYZ(3))*0.125 DSFDXI(1,2)=-(1.D0+XYZ(1))*(1.D0-XYZ(3))*0.125 DSFDXI(1,3)=-(1.D0+XYZ(1))*(1.D0-XYZ(2))*0.125 DSFDXI(2,1)= (1.D0+XYZ(2))*(1.D0-XYZ(3))*0.125 DSFDXI(2,2)= (1.D0+XYZ(1))*(1.D0-XYZ(3))*0.125 DSFDXI(2,3)=-(1.D0+XYZ(1))*(1.D0+XYZ(2))*0.125 DSFDXI(3,1)=-(1.D0+XYZ(2))*(1.D0-XYZ(3))*0.125 DSFDXI(3,2)= (1.D0-XYZ(1))*(1.D0-XYZ(3))*0.125 DSFDXI(3,3)=-(1.D0-XYZ(1))*(1.D0+XYZ(2))*0.125 DSFDXI(4,1)=-(1.D0-XYZ(2))*(1.D0-XYZ(3))*0.125 DSFDXI(4,2)=-(1.D0-XYZ(1))*(1.D0-XYZ(3))*0.125 DSFDXI(4,3)=-(1.D0-XYZ(1))*(1.D0-XYZ(2))*0.125 DSFDXI(5,1)= (1.D0-XYZ(2))*(1.D0+XYZ(3))*0.125 DSFDXI(5,2)=-(1.D0+XYZ(1))*(1.D0+XYZ(3))*0.125 DSFDXI(5,3)= (1.D0+XYZ(1))*(1.D0-XYZ(2))*0.125 DSFDXI(6,1)= (1.D0+XYZ(2))*(1.D0+XYZ(3))*0.125 DSFDXI(6,2)= (1.D0+XYZ(1))*(1.D0+XYZ(3))*0.125 DSFDXI(6,3)= (1.D0+XYZ(1))*(1.D0+XYZ(2))*0.125 DSFDXI(7,1)=-(1.D0+XYZ(2))*(1.D0+XYZ(3))*0.125 DSFDXI(7,2)= (1.D0-XYZ(1))*(1.D0+XYZ(3))*0.125 DSFDXI(7,3)= (1.D0-XYZ(1))*(1.D0+XYZ(2))*0.125 DSFDXI(8,1)=-(1.D0-XYZ(2))*(1.D0+XYZ(3))*0.125 DSFDXI(8,2)=-(1.D0-XYZ(1))*(1.D0+XYZ(3))*0.125 DSFDXI(8,3)= (1.D0-XYZ(1))*(1.D0-XYZ(2))*0.125 C DEFINE JACOBIAN MATRIX FJAC=0.0D0 DO I=1,3 DO J=1,3 DO K=1,NNODE FJAC(I,J)=FJAC(I,J)+DSFDXI(K,I)*COORDS(J,K) ENDDO ENDDO ENDDO AA(1)=((FJAC(1,2)*FJAC(2,3)-FJAC(2,2)*FJAC(1,3))**2.D0 & +(FJAC(1,1)*FJAC(2,3)-FJAC(2,1)*FJAC(1,3))**2.D0 & +(FJAC(1,1)*FJAC(2,2)-FJAC(2,1)*FJAC(1,2))**2.D0)**0.5 AA(2)=((FJAC(2,2)*FJAC(3,3)-FJAC(3,2)*FJAC(2,3))**2.D0 & +(FJAC(2,1)*FJAC(3,3)-FJAC(3,1)*FJAC(2,3))**2.D0 & +(FJAC(2,1)*FJAC(3,2)-FJAC(3,1)*FJAC(2,2))**2.D0)**0.5 AA(3)=((FJAC(1,2)*FJAC(3,3)-FJAC(3,2)*FJAC(1,3))**2.D0 & +(FJAC(1,1)*FJAC(3,3)-FJAC(3,1)*FJAC(1,3))**2.D0 & +(FJAC(1,1)*FJAC(3,2)-FJAC(3,1)*FJAC(1,2))**2.D0)**0.5 RETURN END SUBROUTINE SYMQR (A,D,E,K0,N,NA,EPS,ABSCNV,VEC,TRD,FAIL) c********************************************************************* C C EXPLANATION OF THE PARAMETERS IN THE CALLING SEQUENCE: C C A A DOUBLE DIMENSIONED ARRAY. IF THE MATRIX IS NOT C INITIALLY TRIDIAGONAL, IT IS CONTAINED IN THE LOWER C TRIANGLE OF A. IF EIGENVECTORS ARE NOT REQUESTED, C THE LOWER TRIANGLE OF A IS DESTROYED WHILE THE C ELEMENTS ABOVE THE DIAGONAL ARE LEFT UNDISTURBED. C IF EIGENVECTORS ARE REQUESTED, THEY ARE RETURNED IN THE C COLUMNS OF A. C C D A SINGLY SUBSCRIPTED ARRAY. IF THE MATRIX IS C INITIALLY TRIDIAGONAL, D CONTAINS ITS DIAGONAL C ELEMENTS. ON RETURN, D CONTAINS THE EIGENVALUES OF C THE MATRIX. C C E A SINGLY SUBSCRIPTED ARRAY. IF THE MATRIX IS C INITIALLY TRIDIAGONAL, E CONTAINS ITS OFF-DIAGONAL C ELEMENTS. UPON RETURN, E(I) CONTAINS THE NUMBER OF C ITERATIONS REQUIRED TO COMPUTE THE APPROXIMATE C EIGENVALUE D(I). C C K0 A REAL*8 VARIABLE CONTAINING AN INITIAL ORIGIN SHIFT TO C BE USED UNTIL THE COMPUTED SHIFTS SETTLE DOWN. C C N AN INTEGER VARIABLE CONTAINING THE FIRST DIMENSION C OF THE ARRAY A. C C EPS A REAL*8 VARIABLE CONTAINING A CONVERGENCE TOLERANCE. C C ABSCNV A LOGICAL VARIABLE CONTAINING THE VALUE .TRUE. IF C THE ABSOLUTE CONVERGENCE CRITERION IS TO BE USED C OR THE VALUE .FALSE. IF THE RELATIVE CRITERION C IS TO BE USED. C C VEC A LOGICAL VARIABLE CONTAINING THE VALUE .TRUE. IF C EIGENVECTORS ARE TO BE COMPUTED AND RETURNED IN C THE ARRAY A AND OTHERWISE CONTAINING THE VALUE C .FALSE.. C C TRD A LOGICAL VARIABLE CONTAINING THE VALUE .TRUE. C IF THE MATRIX IS TRIDIAGONAL AND LOCATED IN THE ARRAYS C D AND E AND OTHERWISE CONTAINING THE VALUE .FALSE.. C C FAIL AN INTEGER VARIABLE CONTAINING AN ERROR SIGNAL. C ON RETURN, THE EIGENVALUES IN D(FAIL+1), ..., D(N) C AND THEIR CORRESPONDING EIGENVECTORS MAY BE PRESUMED C ACCURATE. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER N INTEGER NA DIMENSION A(NA,N),D(N),E(N) LOGICAL ABSCNV REAL*8 C REAL*8 CB REAL*8 CC REAL*8 CD REAL*8 CON REAL*8 EPS INTEGER FAIL INTEGER I INTEGER I1 INTEGER J REAL*8 K REAL*8 K0 REAL*8 K1 REAL*8 K2 INTEGER L INTEGER L1 INTEGER LL INTEGER LL1 REAL*8 MAX REAL*8 NINF INTEGER NL INTEGER NM1 INTEGER NM2 INTEGER NNL REAL*8 NORM INTEGER NU INTEGER NUM1 REAL*8 P REAL*8 PP REAL*8 Q REAL*8 QQ REAL*8 R INTEGER RETURN REAL*8 S REAL*8 S2 LOGICAL SHFT REAL*8 SUM REAL*8 SUM1 REAL*8 TEMP REAL*8 TEST REAL*8 TITTER LOGICAL TRD LOGICAL VEC TITTER = 50.0E+00 NM1 = N - 1 NM2 = N - 2 NINF = 0.0E+00 FAIL = 0 D=0.D0 C C SIGNAL ERROR IF N IS NOT POSITIVE. C IF ( N .GT. 0 ) GO TO 1 FAIL = -1 RETURN C C SPECIAL TREATMENT FOR A MATRIX OF ORDER ONE. C 1 IF ( N .GT. 1 ) GO TO 5 IF ( .NOT. TRD ) D(1) = A(1,1) IF ( VEC ) A(1,1) = 1.0E+00 FAIL = 0 RETURN C C IF THE MATRIX IS TRIDIAGONAL, SKIP THE REDUCTION.C C 5 IF ( TRD ) GO TO 100 IF ( N .EQ. 2 ) GO TO 80 C C REDUCE THE MATRIX TO TRIDIAGONAL FORM BY HOUSEHOLDER'S METHOD. C DO 70 L = 1, NM2 L1 = L + 1 D(L) = A(L,L) MAX = 0.0E+00 DO 10 I = L1, N 10 MAX = DMAX1 ( MAX, DABS ( A(I,L) ) ) IF ( MAX .NE. 0.0 ) GO TO 13 E(L) = 0.0E+00 A(L,L) = 1.0E+00 GO TO 70 13 SUM = 0.0E+00 DO 17 I = L1, N A(I,L) = A(I,L) / MAX 17 SUM = SUM + A(I,L)**2 S2 = SUM S2 = DSQRT ( S2 ) IF ( A(L1,L) .LT. 0.0 ) S2 = -S2 E(L) = -S2 * MAX A(L1,L) = A(L1,L) + S2 A(L,L) = S2 * A(L1,L) SUM1 = 0.0E+00 DO 50 I = L1, N SUM = 0.0E+00 DO 20 J = L1, I 20 SUM = SUM + A(I,J) * A(J,L) IF ( I .EQ. N ) GO TO 40 I1 = I + 1 DO 30 J = I1, N 30 SUM = SUM + A(J,L) * A(J,I) 40 E(I) = SUM / A(L,L) 50 SUM1 = SUM1 + A(I,L) * E(I) CON = 0.5E+00 * SUM1 / A(L,L) DO 60 I = L1, N E(I) = E(I) - CON * A(I,L) DO 60 J = L1, I 60 A(I,J) = A(I,J) - A(I,L) * E(J) - A(J,L) * E(I) 70 CONTINUE 80 D(NM1) = A(NM1,NM1) D(N) = A(N,N) E(NM1) = A(N,NM1) C C IF EIGENVECTORS ARE REQUIRED, INITIALIZE A. C 100 IF ( .NOT. VEC ) GO TO 180 C C IF THE MATRIX WAS TRIDIAGONAL, SET A EQUAL TO THE IDENTITY MATRIX. C IF ( .NOT. TRD .AND. N .NE. 2 ) GO TO 130 DO 120 I = 1, N DO 110 J = 1, N 110 A(I,J) = 0.0E+00 120 A(I,I) = 1.0E+00 GO TO 180 C C IF THE MATRIX WAS NOT TRIDIAGONAL, MULTIPLY OUT THE C TRANSFORMATIONS OBTAINED IN THE HOUSEHOLDER REDUCTION.C C 130 A(N,N) = 1.0E+00 A(NM1,NM1) = 1.0E+00 A(NM1,N) = 0.0E+00 DO 170 L = 1, NM2 LL = NM2 - L + 1 LL1 = LL + 1 DO 140 I = LL1, N SUM = 0.0E+00 DO 135 J = LL1, N 135 SUM = SUM + A(J,LL) * A(J,I) 140 A(LL,I) = SUM / A(LL,LL) DO 150 I = LL1, N DO 150 J = LL1, N 150 A(I,J) = A(I,J) - A(I,LL) * A(LL,J) DO 160 I = LL1, N A(I,LL) = 0.0E+00 160 A(LL,I) = 0.0E+00 170 A(LL,LL) = 1.0E+00 C C IF AN ABSOLUTE CONVERGENCE CRITERION IS REQUESTED, C ( ABSCNV = .TRUE. ), COMPUTE THE INFINITY NORM OF THE MATRIX. C 180 IF ( .NOT. ABSCNV ) GO TO 200 NINF = DMAX1 ( & NINF, DABS ( D(1) ) + DABS ( E(1) ) + DABS ( E(I-1) ) ) IF ( N .EQ. 2 ) GO TO 200 DO 190 I = 2, NM1 190 NINF = DMAX1 ( NINF, & DABS ( D(I) ) + DABS ( E(I) ) + DABS ( E(I-1) ) ) C C START THE QR ITERATION. C 200 NU = N NUM1 = N - 1 SHFT = .FALSE. K1 = K0 TEST = NINF * EPS E(N) = 0.0E+00 C C CHECK FOR CONVERGENCE AND LOCATE THE SUBMATRIX IN WHICH THE C QR STEP IS TO BE PERFORMED. C 210 DO 220 NNL = 1, NUM1 NL = NUM1 - NNL + 1 IF ( .NOT. ABSCNV ) & TEST = EPS * DMIN1 ( DABS ( D(NL) ), DABS ( D(NL+1) ) ) IF ( DABS ( E(NL) ) .LE. TEST ) GO TO 230 220 CONTINUE GO TO 240 230 E(NL) = 0.0E+00 NL = NL + 1 IF ( NL .NE. NU ) GO TO 240 IF ( NUM1 .EQ. 1 ) RETURN NU = NUM1 NUM1 = NU - 1 GO TO 210 240 E(NU) = E(NU) + 1.0E+00 IF ( E(NU) .LE. TITTER ) GO TO 250 FAIL = NU RETURN C C CALCULATE THE SHIFT. C 250 CB = ( D(NUM1) - D(NU) ) / 2.0E+00 MAX = DMAX1 ( DABS ( CB ), DABS ( E(NUM1) ) ) CB = CB / MAX CC = ( E(NUM1) / MAX )**2 CD = DSQRT ( CB**2 + CC ) IF ( CB .NE. 0.0E+00 ) CD = SIGN ( CD, CB ) K2 = D(NU) - MAX * CC / ( CB + CD ) IF ( SHFT ) GO TO 270 IF ( DABS ( K2 - K1 ) .LT. 0.5E+00 * DABS ( K2 ) ) GO TO 260 K1 = K2 K = K0 GO TO 300 260 SHFT = .TRUE. 270 K = K2 C C PERFORM ONE QR STEP WITH SHIFT K ON ROWS AND COLUMNS C NL THROUGH NU. C 300 P = D(NL) - K Q = E(NL) CALL SINCOS ( P, Q, C, S, NORM ) 310 DO 380 I = NL, NUM1 C C IF REQUIRED, ROTATE THE EIGENVECTORS. C IF ( .NOT. VEC ) GO TO 330 DO 320 J = 1, N TEMP = C * A(J,I) + S * A(J,I+1) A(J,I+1) = -S * A(J,I) + C*A(J,I+1) 320 A(J,I) = TEMP C C PERFORM THE SIMILARITY TRANSFORMATION AND CALCULATE THE NEXT C ROTATION. C 330 D(I) = C * D(I) + S * E(I) TEMP = C * E(I) + S * D(I+1) D(I+1) = -S * E(I) + C * D(I+1) E(I) = -S * K D(I) = C * D(I) + S * TEMP IF ( I .EQ. NUM1 ) GO TO 380 IF ( ABS ( S ) .GT. ABS ( C ) ) GO TO 350 R = S / C D(I+1) = -S * E(I) + C * D(I+1) P = D(I+1) - K Q = C * E(I+1) CALL SINCOS ( P, Q, C, S, NORM ) 340 E(I) = R * NORM E(I+1) = Q GO TO 380 350 P = C * E(I) + S * D(I+1) Q = S * E(I+1) D(I+1) = C * P / S + K E(I+1) = C * E(I+1) CALL SINCOS ( P, Q, C, S, NORM ) 360 E(I) = NORM 380 CONTINUE TEMP = C * E(NUM1) + S * D(NU) D(NU) = -S * E(NUM1) + C * D(NU) E(NUM1) = TEMP GO TO 210 END SUBROUTINE SINCOS ( P, Q, C, S, NORM ) c*********************************************************************72 C C SINCOS CALCULATES THE ROTATION CORRESPONDING TO THE VECTOR (P,Q). C IMPLICIT NONE REAL*8 C REAL*8 NORM REAL*8 P REAL*8 PP REAL*8 Q REAL*8 QQ REAL*8 S 500 PP = ABS ( P ) QQ = ABS ( Q ) IF ( QQ .GT. PP ) GO TO 510 NORM = PP * DSQRT ( 1.0E+00 + ( QQ / PP )**2 ) GO TO 520 510 IF ( QQ .EQ. 0.0E+00 ) GO TO 530 NORM = QQ * DSQRT ( 1.0E+00 + ( PP / QQ )**2 ) 520 C = P / NORM S = Q / NORM RETURN 530 C = 1.0E+00 S = 0.0E+00 NORM = 0.0E+00 RETURN END SUBROUTINE FORM_H5MATRIX(NNODE,BL,DDSDDMU,TN2,H5MATRIXT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION BL(6,3*NNODE),DDSDDMU(6,1),TN2(1,NNODE) DIMENSION H5MATRIXT(3*NNODE,NNODE),TEMP(3*NNODE,1) CALL TTMULTIPLY(BL,DDSDDMU,TEMP,6,3*NNODE,1) CALL TMULTIPLY(TEMP,TN2,H5MATRIXT,3*NNODE,1,NNODE) RETURN END