PEATP

CCC PROGRAM FEATPC====================================================================C THIS PROGRAM CAN SOLVE THE ELASTIC PROBLEM OF PLANE STRESS,C PLANE STRAIN,AXISYMMETRIC SOLID AND MINDLIN PLATEC====================================================================C 7 SUBROUTINES(ALLOCAT,INPUT,ASSEM,STATIC_SOLVE,STRESS,DYNAM,C EIGEN)ARE CALLED BY THE MAIN PROGRAM.C BESIDES,ANOTHER 18 SUBROUTINES ARE CALLED BY ABOVE 7 SUBROUTINES.C======================================================================= IMPLICIT REAL*8(A-H,O-Z) CHARACTER*5 FILENAMEDIMENSION IZ(300000),AR(15000000)C !IZ-MAXIMUN SPACE FOR DYNAMIC INTEGE ARRAYC !AR-MAXIMUM SPACE FOR DYNAMIC REAL ARRAY COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI COMMON/COMN/NFIX,NPC,GRAV COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA COMMON/COM3/MND2,NUMPT2 COMMON/ELEM/NODE,INTX,INTY COMMON/DYN/HUV,OMEGA,CC1,CC2,TT,DT,ALPHA,DELTA print*,'===> Please Input the Title of the Problem to be Solved!' READ(*,'(a5)')FILENAMEOPEN(5,FILE=FILENAME//'.inp',STATUS='OLD') OPEN(6,FILE=FILENAME//'.dat',STATUS='UNKNOWN') OPEN(10,FILE=FILENAME//'.mkp',STATUS='UNKNOWN') OPEN(14,FILE=FILENAME//'.dis',STATUS='UNKNOWN') OPEN(15,FILE=FILENAME//'.str',STATUS='UNKNOWN')C*********************************************************************************************************************************C ALLOCATE STORAGE SPACE FOR THE ARRAYS OF FE MODELC***************************************************************************************************** CALL ALLOCAT(M1,M2,M3,M4,M5,M8,M11,M12,N1,N2,N3,N4,N5,N6,N7, $ N8,N9,N10,N11,N12,N13,N14,N15,N16,N17,N18,N19,N20,N21,N22, $ N23,N24,N25,N26,N27,N31,N32,N33,N34,N35,N36,N37,N38)C***************************************************************************************************C INPUT ALL THE DATA OF FE MODEL AT APPOINTED ADDRESSC**************************************************************************************************** CALL INPUT(IZ(M2),IZ(M3),IZ(M4),AR(N1),AR(N2),AR(N3),AR(N4))C****************************************************************************************************C ASSEMBLE THE ELEMENT MATRIXES TO FORM GLOBAL MATRIXES:[M],[K],{P}C******************************************************************************************************* CALL ASSEM(AR(N1),IZ(M2),AR(N2),IZ(M3),AR(N3),IZ(M4), $ AR(N4),AR(N11),AR(N12),AR(N13),AR(N14), $ IZ(M11),AR(N15),AR(N16),AR(N17))C********************************************************************************************************C GET NODAL DISPLACEMENTS BY SOLVING EQUATION[K]{U}={P}C FOR STATIC PROBLEM(MSOLVE=1)C******************************************************************************************************** IF(MSOLV.EQ.1)THEN CALL STATIC_SOLVE(AR(N12),AR(N13),AR(N14))C*********************************************************************************************************C COMPUTE STRESSES AT THE INTEGRATION POINTS AND NODES C*********************************************************************************************************** CALL STRESS(AR(N1),IZ(M2),AR(N2),AR(N14),AR(N15),IZ(M11), $ AR(N18),AR(N9),AR(N10),IZ(M8),AR(N8)) ENDIFC****************************************************************************************************C SOLVE DYNAMIC RESPONSE PROBLEM BY CENTRAL-DIFFERENCE METHOD C (MSOLVE=2) AND BY NEWMARK METHOD(MSOLV=3)C******************************************************************************************************* IF(MSOLV.EQ.2.OR.MSOLV.EQ.3) $CALL DYNAM(AR(N11),AR(N12),AR(N13),AR(N21),AR(N22), $ AR(N23),AR(N24),AR(N25),AR(N26))C************************************************************************************************************C SOLVE DYNAMIC CHARACTER PROBLEM BY THE INVERSE ITERATION METHOD C (MSOLV=4) AND BY THE SUBSPACE ITERATION METHOD(MSOLV=5)C************************************************************************************************************ IF(MSOLV.EQ.4.OR.MSOLV.EQ.5) $CALL EIGEN(AR(N11),AR(N12),AR(N31),AR(N32),AR(N33), $ AR(N34),AR(N35),AR(N36),AR(N37)) STOP ENDC=========================================================================C=======================SUB:1============================================= SUBROUTINE ALLOCAT(M1,M2,M3,M4,M5,M8,M11,M12,N1,N2,N3,N4,N5, $ N6,N7,N8,N9,N10,N11,N12,N13,N14,N15,N16,N17,N18,N19,N20, $ N21,N22,N23,N24,N25,N26,N27,N31,N32,N33,N34,N35,N36,N37,N38) C***********************************************************************************************************C INPUT BASIC PAPRMETERS FROM FILE'IN_DAT'C ALLOCAT DYANMICAL STORAGE SPACEC*********************************************************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA COMMON/COMN/NFIX,NPC,GRAV COMMON/COM3/MND2,NUMPT2C !MND-MAXIMAL NODE NUMBER IN ALL ELEMENTSC !NUMEL-NEMBER OF GLOBAL ELEMENTS C !NUMPT-NUMBER OF GLOBAL NODESC !MBAND-HALF BANDWIDTH(INCLUDING DIAGONAL ELEMENT) READ(5,*)MND,NUMEL,NUMPT,MBAND READ(5,*)NFIX,NPC,MPROB,MSOLVC !NFIX-NUMBER OF NODES SUBJECTER TO CONSTRIANTC !NPC-NUBMER OF NODES SUBJECTED TO EQUIVALENT LOADC==========================================================C MPROB=1-PLANE STRESS PROBLEM, MPROB=2-PLANE STRAIN PROBLEMC MPROB=3-AXISYMMETRIC PROBLEM, MPROB=4-MINDLIN PLATE PROBLEMC------------------------------------------------------------------------------------------------------------------------------C MSOLV=1-STATIC ANALYSISC MSOLV=2-DYNAMIC RESPONSE ANALYSIS BY CENTRAL DIFFERENCE METHODC MSOLV=3-DYNAMIC RESPONSE ANALYSIS BY NEWMARK METHODC MSOLV=4-DYNAMIC CHARACTER ANALYSIS BY INVERSE ITERATION METHODC MSOLV=5-DYNAMIC CHARACTER ANALYSIS BY SUBSPACE ITERATION METHODC================================================================== IF(MSOLV.NE.4.OR.MSOLV.NE.5)READ(5,*)NMATI,GRAV,MTYPE IF(MSOLV.EQ.4.OR.MSOLV.EQ.5)READ(5,*)NMATI,GRAV,MTYPE,NVAC !NMATI-KIND OF MATERIALSC !GRAV-GRAVITY ACCELERATIONC !NVA-NUMBER OF EIGENVALUESC----------------------------------------------------------------------------------------------------------------------------------------C MTYPE-CONTROL KEY FOR OUTPUT RESULTS C MTYPE=0-OUTPUT RESULTS INGLUDE GLOBAL MATRIXES AND STRESS AT C GAUSS POINTSC MTYPE=1-OUTPUT RESULTS INCLUDE GLOBAL MATRIXES C MTYPE=2-OUTPUT RESULTS INCLUDE STRESSES AT GAUSS POINTSC------------------------------------------------------------------------------------------------------------------------------------------- IF(MPROB.EQ.1.OR.MPROB.EQ.2) THEN NF=2 !NF-NUMBER OF NODAL FREEDOMS(DISPLACEMENT COMPONENTS) NFSTR=3 !NFSTR-NUMBER OF STRESS COMPONENTS ENDIF IF(MPROB.EQ.3)THEN NF=2 NFSTR=4 ENDIF IF(MPROB.EQ.4)THEN NF=3 NFSTR=5 ENDIFC*************************************************************************************************************************C ALLOCATE STORAGE SPACE FOR INPUT DATAC********************************************************************************************************************** M1=1 M2=M1 M3=M2+NUMEL*14 !IELEM(NUMEL,14)-CODE OF MATERIAL AND NODES M4=M3+NFIX*4 !IFIXD(NFIX,4)-CODE OF NODES HAMING CONSTRIANT M5=M4+NPC*4 !ILOAD(NPC,4)-CODE OF NODES HAVING LOAD N1=1 N2=N1+NMATI*4 !VAMATI(NMAT1,4)-PARAMETERS OF MATERIALS N3=N2+NUMPT*2 !VCOOD(NUMPT,2)-GLOBAL NODE COORDINATES N4=N3+NFIX*3 !VFIXD(NFIX,3)-CONSTRIANED DISPLACEMENTS N5=N4+NPC*3 !VLOAD(NPC,3)-VALUES OF EQUIVALENT LOADC****************************************************************************************************************************C ALLOCATE STORAGE SPACE FOR ELEMENTAL MATRIXES AND GLOBAL MATRIXESC****************************************************************************************************************************** MND2=MND*NF !MND2-NUMBER OF FREEDOMS IN A ELEMENT NUMPT2=NUMPT*NF !NUMPT2-NUMBER OF GLOBAL FREEDOMS M8=M5 M11=M8+NUMPT !IADD(NUMPT)-USED FOR NODAL STRESS M12=M11+MND !IEL(MND)-NODE CODE IN A ELEMENT N6=N5 N7=N6+NFSTR*MND2 !VSG(NFSTR,MND2)-ELEMENT STRESS MATRIX AT !INTEGRATION POINTS N8=N7+NFSTR*MND2 !VSN(NFSTR,MND2)-ELEMENT STRESS AT NODES N9=N8+NUMPT*NFSTR !SSN(NUMPT,NFSTR)-STRESSES AT NODES N10=N9+9*NFSTR !SSS(9,NFSTR)-STRESSES AT INTERGARION POINTSN11=N10+4*NFSTR !VSS(4,NFSTR)-STRESSES AT DESIGNATED !INTEGRATION POINTS` N12=N11+NUMPT2*MBAND !GMM(NUMPT2,MBAND)-GLOBAL MASS MATRIX N13=N12+NUMPT2*MBAND !GKK(NUMPT2,MBAND)-GLOBAL STIFENESS MATRIX N14=N13+NUMPT2 !GP(NUMPT2)-GLOBAL LOAD VECTOR N15=N14+NUMPT2 !GU(NUMPT2)-GLOBAL DISPLACEMENT VECTOR N16=N15+MND*2 !VXY(MND,2)-NODE COORDINATE IN ELEMENT N17=N16+MND2*MND2 !VMM(MND2,MND2)-ELEMENT MASS MATRIX N18=N17+MND2*MND2 !VKK(MND2,MND2)-ELEMENT STIFFNESS MATRIX N19=N18+MND2 !VU(MND2)-NODE DISPLACEMENTS IN A ELEMENT N20=N19+4 !VMAE(4)-MATERIAL PARAMETERS IN A ELEMENTC******************************************************************************************************************************C ALLOCAT STORAGE SPACE FOR DYNAMIC RESPONSE ANALYSISC********************************************************************************************************************************* IF(MSOLV.EQ.2.OR.MSOLV.EQ.3)THEN N21=N20 N22=N21+NUMPT2*MBAND !DAMP(NUMPT2,MBAND)-DAMPLING MATRIX N23=N22+NUMPT2 !U0(NUMPT2)-INITIAL DISPLACEMENT VECTOR N24=N23+NUMPT2 !V0(NUMPT2)-INITIAL VELOCITY VECTOR N25=N24+NUMPT2 !A(NUMPT2)-INITIAL ACCELERATION VECTOR N26=N25+NUMPT2*MBAND !AW(NUMPT2,MBAND)-WORKING ARRAY N27=N26+NUMPT2 !B(NUMPT2)-WORKING ARRAY ENDIFC*************************************************************************************************************************************C ALLOCAT STORAGE SPACE FOR DYNAMIC CHARACTER ANALYSISC*************************************************************************************************************************************** IF(MSOLV.EQ.4.OR.MSOLV.EQ.5)THEN IF(NVA.EQ.0)THEN WRITE(*,11) 11 FORMAT('PLEASE READ THE NUMBERS OF EIGENVALUE-NVA=') READ(*,*)NVA ENDIF N31=N20 N32=N31+NUMPT2*NVA !AA(NUMPT2,NVA)-INITIAL ITERATION VECTOR N33=N32+NUMPT2*NVA !BB(NUMPT2,NVA)-WORKING ARRAY N34=N33+NVA*NVA !GM(NVA,NVA)-MASS MATRIX IN SUBSPACE N35=N34+NVA*NVA !GK(NVA,NVA)-STIFFNESS MATRIX IN SUBSPACE N36=N35+NVA*NVA !V(NVA,NVA)-EIGENVECTORS IN SUBSPACE N37=N36+NVA !W1(NVA)-WORKING ARRAY N38=N37+NVA !W2(NVA)-EIGENVALUES IN SUBSPACE ENDIF RETURN ENDC=======================SUB:2========================================================= SUBROUTINE INPUT(IELEM,IFIXD,ILOAD,VMATI,VCOOD,VFIXD,VLOAD)C*****************************************************************************************************************************C INPUT ALL THE INFORMATION OF FE MODELC******************************************************************************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI COMMON/COMN/NFIX,NPC,GRAV COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA COMMON/DYN/HUV,OMEGA,CC1,CC2,TT,DT,ALPHA,DELTA DIMENSION IELEM(NUMEL,14),IFIXD(NFIX,4),ILOAD(NPC,4), $ VCOOD(NUMPT,2),VFIXD(NFIX,3),VLOAD(NPC,3), $ VMATI(NMATI,5) WRITE(*,101) 101 FORMAT(/6X,'## IMPUT DATA FROM FILETO MEMORY #')C======================================================================================== !INPUT NODAL COORDINATES READ(5,*)(II,(VCOOD(I,J),J=1,2),I=1,NUMPT) !INPUT ELEMENT INFORMATION READ(5,*)(II,(IELEM(I,J),J=1,4+MND),I=1,NUMEL) !IMPUT CONSTRAIN INFORMATION READ(5,*)(II,(IFIXD(I,J),J=1,NF+1),(VFIXD(I,J),J=1,NF),I=1,NFIX) IF(NPC.GT.0)THEN !INPUT EQUIVALENT LOAD AT NODES READ(5,*)(II,(ILOAD(I,J),J=1,NF+1),(VLOAD(I,J),J=1,NF),I=1,NPC) ENDIF !INPUT MATERIAL PARAMETERS READ(5,*)(II,(VMATI(I,J),J=1,3),I=1,NMATI) IF(MSOLV.EQ.2.OR.MSOLV.EQ.3) THEN READ(5,*) READ(5,*) READ(5,*) HUV,FREQ,CC1,CC2,TT,DT,ALPHA,DELTAC !HUV-INPUT CONTROL OF INITIAL DISPLACEMENT AND VELOCITYC (HUV=1:INPUT DISPLACEMENT;HUV=2:INPUT VELOCITY;C HUV=3:INPUT BOTH OF THEM)ADC !OMEGA-FREQUENCE OF LO C !CC1,CC2-CONSTANTS TO COMPUTE DAMPLING MATRIXSC !TT-TOTAL TIMEC !DT-TIME STEP LENGTHC !ALPHA,DELTA-CONSTANTS USED IN THE NEWMARK METHOD ENDIFC***************************************************************************************C OUTPUT ABOVE INPUT DATAC**************************************************************************** WRITE(*,102)102 FORMAT(/5X,'%% OUTPUT INPUT-DATA TO %%') WRITE(6,10)10 FORMAT(/8X,'MAXIMAL-ELEM-NODES,ELEMENTS,NODES,BANDWIDTH') WRITE(6,11) MND,NUMEL,NUMPT,MBAND11 FORMAT(10X,I10,I11,I9,I9) WRITE(6,12)12 FORMAT(/8X,'FIXED-NODES,EQUIVALENT-LOADS,MATERIAL-KINDS,', $ 2X,'GRAVITY') WRITE(6,13) NFIX,NPC,NMATI,GRAV13 FORMAT(11X,I9,2(3X,I9),2X,F16.5) WRITE(6,14)14 FORMAT(/8X,'PROBLEM-KIND,SOLVING-KIND,OUTPUT-KEY') WRITE(6,15) MPROB,MSOLV,MTYPE15 FORMAT(11X,I9,2(3X,I9)) WRITE(6,16)16 FORMAT(/8X,'NODAL-FREEDOMS,STRESS-COMPONENTS') WRITE(6,17) NF,NFSTR17 FORMAT(11X,I9,3X,I10) WRITE(6,18) ! OUTPUT NODAL COORDIATES18 FORMAT(/8X,'NODAL COORDINATES'/9X'NO.',15X,'X-',11X,'Y-') DO 19 I=1,NUMPT19 WRITE(6,20)I,(VCOOD(I,J),J=1,2)20 FORMAT(5X,I5,5X,3F15.6) WRITE(6,21)(I,I=1,9) ! OUTPUT ELEMENT INFORMATION 21 FORMAT(/8X,'ELEMENT INFORMATION'/3X,'NP.',1X,'NODES',1X, $ 'MATERIAL',1X,'INTX',1X,'INTY',1Z,9(2X,2HN-,I1)) DO 22 I=1,NUMEL22 WRITE(6,23) I,(IELEM(I,J),J=1,MND+4)23 FORMAT(1X,I5,2I5,3X,2I5,3X,9(1X,I4)) WRITE(6,24) ! OUTPUT XONSTRAIN INFORMATION24 FORMAT(/8X,'CONSRAINT INFORMATION ON NODES'/ $ 11X,'NO.',2X,'NODE NO.',2X,'X-',1X,'Y-',1X,'Z-', $ 5X,'X-VALUE',2X,'Y-VALUE',2X,'Z-VALUE') DO 25 I=1,NFIX25 WRITE(6,26) I,(IFIXD(I,J),J=1,4),(VFIXD(I,J),J=1,3)26 FORMAT(10X,I4,4X,I3,3X,3I3,3X,3E10.3) IF(NPC.NE.0) THEN WRITE(6,27) ! OUTPUT EQUIVALENT LOAD AT NODES27 FORMAT(/8X,'EQUIVALENT LOAD ON NODES'/12X,'NO.',1X, $ 'NODE NO.',2X,'X-',1X,'Y-',1X,'Z-',5X,'X-VALUE',2X, $ 'Y-VALUE',2X,'Z-VALUE') DO 28 I=1,NPC28 WRITE(6,29) I,(ILOAD(I,J),J=1,4),(VLOAD(I,J),J=1,3)29 FORMAT(10X,I4,4X,I3,4X,3I3,3X,3E10.3) ENDIF WRITE(6,33)33 FORMAT(/8X,'MATERIAL PARAMETERS'/5X,'NO.'3X,'E(MODULUS)',2X, $ 'V(POISSON RATIO)',2X,'DENS(DENSITY)',2X,'TH(THICKNESS)') DO 34 I=1,NMATI34 WRITE(6,35) I,(VMATI(I,J),J=1,4)35 FORMAT(5X,I2,4E14.3)C*******C********OUTPUT THE INPUT PARAMETERS FOR DYNAMIC REPONSE COMPUTATIONC********** IF(MSOLV.EQ.2.OR.MSOLV.EQ.3) $WRITE(6,40) HUV,OMEGA,CC1,CC2,TT,DT,ALPHA,DELTA40 FORMAT(/8X,'PARAMETERS FOR DYNAMIC RESPONSE COMPUTATION:'/ $ 5X,'INPUT CONTROL FOR INITIAL DISPLACE AND VELOCITY-HUV=',I3/ $ 5X,'FREQUENCE OF LOAD-OMEGA=',E10.3/ $ 5X,'DAMPING COEFFICIENTS-CC1=',E10.3,8X,'CC2=',E10.3/ $ 5X,'TOTAL TIME-TT=',E12.6,10X,'STEP LENGTH-DT=',E12.6/ $ 5X,'PARAMETERS OF NEWMARK-ALPHA=',E10.3,8X,'DELTA=',E10.3)C**********C*******OUTPUT THE INPUT PARAMETER FOR CHARACTER VALUE COMPUTATIONC********* IF(MSOLV.EQ.4.OR.MSOLV.EQ.5) WRITE(6,45) NVA45 FORMAT(/8X,'PARAMERER FOR DYNAMIC CHARACTER COMPUTATION:'/ $ 10X,'NUMBERS OF EIGENVALUES-NVA=',I3) WRITE(6,301)301 FORMAT(/'C************************************************', $ '******************************',/) RETURN ENDC=========================SUB:3===================================== SUBROUTINE ASSEM(VMATI,IELEM,VCOOD,IFIXD,VFIXD, $ ILOAD,VLOAD,GMM,GKK,GP,GU,IEL,VXY,VMM,VKK)C*******************************************************************C ASSEMBLE THE GLOBAL MATRIXES:[M],[K],AND{P}C CALL ELEMENT-MATRIXC***************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI COMMON/COMN/NFIX,NPC,GRAV COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA COMMON/COM3/MND2,NUMPT2 COMMON/ELEM/NODE,INTX,INTY DIMENSION IELEM(NUMEL,MND+4),VCOOD(NUMPT,2),IFIXD(NFIX,4), $ VFIXD(NFIX,3),ILOAD(NPC,4),VLOAD(NPC,3), $ GMM(NUMPT2,MBAND),GKK(NUMPT2,MBAND),GP(NUMPT2), $ GU(NUMPT2),VXY(MND,2),IEL(MND),VMATI(NMATI,4), $ VMAE(4),VMM(MND2,MND2),VKK(MND2,MND2) WRITE(*,101)101 FORMAT(/5X,'## ASSEMBLE GLOBAL MATRIX[GKK],[GMM]##')C*********C**********GLEAR THE MEMORY FOR GLOBAL MATRIX:[M]AND[K]C******** DO 10 I=1,NUMPT2 DO 10 J=1,MBAND GKK(I,J)=0.0 GMM(I,J)=0.010 CONTINUEC*******C*********LOOP OVER EACH ELEMENTC********* DO 320 IE=1,NUMEL !NUMEL-NUMBER OF ELEMENTSC********C*********FORM INFORMATION OF EACH ELEMENT FROM THE INPUT DATAC******* DO 11 I=1,MND !MND-MAXIMUM NODE NUMBER IN ALL ELEMENTS IEL(I)=IELEM(IE,I+4) DO 11 J=1,2 VXY(I,J)=0.0 IF(IEL(I).GT.0) VXY(I,J)=VCOOD(IEL(I),J)11 CONTINUE NODE=IELEM(IE,1) INTX=IELEM(IE,3) INTY=IELEM(IE,4) IMATI=IELEM(IE,2) DO 13 J=1,4 VMAE(J)=VMATI(IMATI,J)13 CONTINUEC****************************************************************************C COMPUTE ELEMENTAL MASS MATRIX[VMM] AND STIFFNESS MATRIX[VKK]C FOR THE ELEMENTS WITH 3-6,4-8,9 NODESC********************************************************************************C CALL ELEMENT_MATRIX(IE,VXY,IEL,VMM,VKK,VMAE)CC*************************************************************************C ASSEMBLE ELEMENT MATRIXES TO FORM THE GLOBAL MASS MATRIX[GMM]C AND GLOBAL STIFFNESS MATRIX[GKK]C***************************************************************************** DO 20 I=1,MND DO 20 J=1,MND DO 25 II=1,NF DO 25 JJ=1,NF IH=NF*(I-1)+II IV=NF*(J-1)+JJ IHH=NF*(IEL(I)-1)+II IVV=NF*(IEL(J)-1)+JJ IVV=IVV-IHH+1 IF(IHH.GT.0.AND.IVV.GT.0) $ GKK(IHH,IVV)=GKK(IHH,IVV)+VKK(IH,IV) IF(IHH.GT.0.AND.IVV.GT.0) $ GMM(IHH,IVV)=GMM(IHH,IVV)+VMM(IH,IV)25 CONTINUE20 CONTINUE320 CONTINUE WRITE(*,102)102 FORMAT(/6X,'## DRAW BOUNDARY COMNDITIONS TO FORM [GP]##')C******************************************************************************C FORM THE GRAVITY LOADING [P]C******************************************************************************* DO 30 I=1,NUMPT DO 30 J=1,NF GP(NF*(I-1)+J)=0.0 GU(NF*(I-1)+J)=0.0 IF(J.EQ.NF) GU(NF*(I-1)+J)=GRAV !GRAV-GRAVITY ACCELERATION30 CONTINUE DO 31 I=1,NUMPT2 GP(I)=GMM(I,1)*GU(I) DO 31 K=I+1,I+MBAND-1 IF(K.LE.NUMPT2) GP(I)=GP(I)+GMM(I,K-I+1)*GU(K) IF(2*I-K.GE.1) GP(I)=GP(I)+GMM(2*I-K,K-I+1)*GU(2*I-K)31 CONTINUE IF(GRAV.NE.0.0) THEN WRITE(6,*)'LOAD UNDER GRAVITY:' DO 32 I=1,NUMPT32 WRITE(6,74) I,(GP(NF*(I-1)+J),J=1,NF) WRITE(6,301) ENDIFC*********************************************************************************C ADD THE NODAL LOAD VECTOR[P]C******************************************************************************** DO 40 I=1,NPC DO 40 J=1,NF IF(ILOAD(I,J+1).NE.0) THEN II=NF*(ILOAD(I,1)-1)+J GP(II)=GP(II)+VLOAD(I,J) ENDIF40 CONTINUEC************************************************************************************C DRAW THE NODAL CONSTRAINTC*********************************************************************************** IF(MSOLV.EQ.1) THEN DO 42 I=1,NFIX DO 42 J=1,NF IF(IFIXD(I,J+1).NE.0) THEN II=NF*(IFIXD(I,1)-1)+J GU(II)=VFIXD(I,J)*1E30 GKK(II,1)=GKK(II,1)*1E30 IF(GKK(II,1).LE.1E-10) GKK(II,1)=0. ENDIF42 CONTINUE ENDIF IF(MSOLV.GT.1) THEN DO 50 I=1,NFIX DO 50 J=1,NF IF(IFIXD(I,J+1).NE.0) THEN II=NF*(IFIXD(I,1)-1)+J GU(II)=0.0 GKK(II,1)=1.0 GMM(II,1)=1.0 IF(MSOLV.EQ.4.OR.MSOLV.EQ.5) GMM(II,1)=0.0 DO 60 K=2,MBAND GKK(II,K)=0.0 GMM(II,K)=0.0 IF(II.LT.K) GOTO 60 GKK(II-K+1,K)=0.0 GMM(II-K+1,K)=0.060 CONTINUE ENDIF50 CONTINUE ENDIFC****************************************************************************************C OUTPUT THE MATRIXES:[M],[K] AND{P} TO FILE 'OUT_MKP'C (WHEN MTYPE=0 OR MTYPE=1)C*************************************************************************************** CLOSE(10,STATUS='DELETE') IF(MTYPE.EQ.0.OR.MTYPE.EQ.1) THEN WRITE(*,105)105 FORMAT(/6X,'%% OUTPUT GLOBAL MATRIX INTO%%')C*******C****************OUTPUT THE GLOBAL MASS MATRIX:[M]C******** OPEN(10,FILE='OUT_MKP',STATUS='UNKNOWN') WRITE(10,301) WRITE(10,*)'GLOBAL MASS MATRIX(GMM):' DO 70 I=1,NUMPT DO 70 J=1,NF II=NF*(I-1)+J WRITE(10,71) I,II,(GMM(II,K),K=1,MBAND) 70 CONTINUE71 FORMAT(2I5,20(3X,2E15.6))C***********C******************OUTPUT THE GLOBAL STIFFNESS MATRIX[K]C********* WRITE(10,301) WRITE(10,*) 'GLOBAL STIFFNESS MATRIX (GKK):' DO 72 I=1,NUMPT DO 72 J=1,NF II=NF*(I-1)+J WRITE(10,71) I,II,(GKK(II,K),K=1,MBAND)72 CONTINUEC*********C******************OUTPUT THE GLOBAL LOADING{P}C******** WRITE(10,301) WRITE(10,*)'GLOBAL LOADING(GP):' DO 73 I=1,NUMPT WRITE(10,74) I,(GP(NF*(I-1)+J),J=1,NF)73 CONTINUE74 FORMAT(5X,I6,4X,3E15.6) WRITE(10,301)301 FORMAT(/3X,'********************', $ '***************************',/) CLOSE(10) ENDIF RETURN ENDC=========================SUB:3-1========================================================== SUBROUTINE ELEMENT_MATRIX(IE,VXY,IEL,VMM,VKK,VMAE)C***************************************************************************************C FORM THE ELEMENT MATRIXES:[M],[K] AND[S]C CALL SUBROUTINES OF ELEMENT_VD,GAUSS_INTEGRAATION ORC HAMMER_INTEGATION,SHAPE_QUADRANGLE_8 OR SHAPE_QUADRANGLE_9C OR SHAPE_TRIANGLE,ELEMENT_VB,ELEMENT_JOCABIC******************************************************************************************* IMPLICIT REAL*8(A-H,O-Z) COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA COMMON/COM3/MND2,NUMPT2 COMMON/ELEM/NODE,INTX,INTY DIMENSION VXY(MND,2),IEL(MND),VMAE(4),VMM(MND2,MND2), $ VKK(MND2,MND2),VSG(NFSTR,MND2) DIMENSION VN(9),VDN(3,9),VD0(3,9),VD(5,5),VB(5,27)C*******C**********FORM THE[D] MATRIX ACCORDING TO TYPE OF THE PROBLEMC******** CALL ELEMENT_VD(MPROB,VD,VMAE) !VD-ELASTIC MATRIXC********C**********CLEAR THE MEMORY FOR ELEMENT MATRIXES[M]AND[K]C******** DO 10 I=1,MND2 DO 10 J=1,MND2 VMM(I,J)=0.0 VKK(I,J)=0.010 CONTINUEC**********************************************************************************************C COMPUTE SHAPE FUNCTION VN AND ITS LOCAL DERIVATIVE VDN AT EACHC INTEGRATOPM POINTS IN ELEMENTSC****************************************************************************************** IF(NODE.EQ.3.OR.NODE.EQ.6) INTY=1 DO 302 I=1,INTX !INTX-INTEGRATION POINTS IN X DIRECTION DO 302 J=1,INTY !INTY-INTEGRATION POINTS IN Y DIRECTION IF(NODE.EQ.4.OR.NODE.EQ.8) THEN CALL GAUSS_INTEGRATION(INTX,INTY,I,J,X,Y,WXY) CALL SHAPE_QUADRANGLE_8(NODE,X,Y,IEL,VN,VDN)C !VN-SHAPE FUNCTIONC !VDN-LOCAL DERIVATE OF VN ENDIF IF(NODE.EQ.9) THEN CALL GAUSS_INTEGRATION(INTX,INTY,I,J,X,Y,WXY) CALL SHAPE_QUADRANGLE_9(NODE,X,Y,IEL,VN,VDN) ENDIF IF(NODE.EQ.3.OR.NODE.EQ.6) THEN CALL HAMMER_INTEGRATION(INTX,I,X,Y,Z,WXY) CALL SHAPE_TRIANGLE(NODE,X,Y,Z,IEL,VN,VDN) ENDIFC******C*********************FORM THE JACOBI MAXTRIX[J] AT INTEGRATION POINTSC****** CALL ELEMENT_JACOBI(MND,VXY,VDN,SJ,VD0)C !SJ=|J|C !VD0-GLOBAL DERIVATIVE OF VNC !VDN-LOCAL DERIVATE OF VN IF(SJ.LE.0.0) THEN WRITE(*,99)IE,I,J,SJ99 FORMAT(/3X,'***SJ.LE.0.0 IN ELEMENT=',I4,3X,'INTX=',I2, $ 3X,'INTY=',I2,3X,'SJ=',E10.4) ENDIFC***C**********FORM THE[B]MATRIX AT INTEGRATION POINTSC*** CALL ELEMENT_VB(MPROB,MND,VXY,VN,VD0,VB,SR) !VB-ELEMENT STRAIN MATRIX !SR-PARAMETER OF INTEGRATIONC***C*****************FORM THE ELEMENT STRESS MATRIX[S]=[D]*[B]C*** DO 20 II=1,NFSTR DO 20 JJ=1,MND2 VSG(II,JJ)=0.0C !VSG-ELEMENT STRESS MATRIX AT INTEGRATION POINT DO 20 KK=1,NFSTR VSG(II,JJ)=VSG(II,JJ)+VD(II,KK)*VB(KK,JJ)20 CONTINUEC***C*****************FORM ELEMENT STIFFNESS MATRIX:[K]=[B]*[S]C*** DO 22 II=1,MND2 DO 22 JJ=1,MND2 DO 22 KK=1,NFSTR VKK(II,JJ)=VKK(II,JJ)+VB(KK,II)*VSG(KK,JJ)*WXY*SJ*SRC !VKK-ELEMENT STIFFNESS MATRIX22 CONTINUE C***C******************FORM THE ELEMENT MASS MATRIX[M]C*** DENS=VMAE(3) DO 25 II=1,MND DO 25 JJ=1,MND DO 25 KK=1,NF II1=NF*(II-1)+KK JJ1=NF*(JJ-1)+KK VMM(II1,JJ1)=VMM(II1,JJ1)+DENS*VN(II)*VN(JJ)*WXY*SJ*SRC !VMM-ELEMENT MASS MATRIX25 CONTINUE302 CONTINUE RETURN ENDC================================SUB3-1-1================================== SUBROUTINE GAUSS_INTEGRATION(INTX,INTY,I,J,X,Y,WXY)C*******************************************************************C GET THE INFORMATION OF GAUSS INTEGRATION POINTC**************************************************************8 IMPLICIT REAL*8(A-H,O-Z) DIMENSION GXY(3,3),WG(3,3)C***C********************GAUSS INTEGRATION CONSTANTS FOR 1,2 AND 3POINTSC*** GXY(1,1)=0.0 WG(1,1)=2.0 GXY(1,2)=-0.577350269189626 GXY(2,2)=0.577350269189626 WG(1,2)=1.0 WG(2,2)=1.0 GXY(1,3)=-0.774596669241483 GXY(2,3)=0.0 GXY(3,3)=0.774596669241483 WG(1,3)=0.555555555555556 WG(2,3)=0.888888888888889 WG(3,3)=0.555555555555556C***C*******************GET PARAMETERS OF INTEGRATION POINTC*** X=GXY(I,INTX) Y=GXY(J,INTY) WXY=WG(I,INTX)*WG(J,INTY) RETURN ENDC=======================SUB:3-1-2====================================================== SUBROUTINE SHAPE_QUADRANGLE_8(NODE,X,Y,IEL,VN,VDN)C***********************************************************************C COMPUTE SHAPE FUNCTION OF QUADRANGLE ELEMENT AT INTEGRATION POINTC******************************************************************************* IMPLICIT REAL*8(A-H,O-Z) DIMENSION IEL(NODE),VN(9),VDN(3,9)C***C******************CLEAR THE MEMORY FOR SHAPE FUNCTION AND ITS DERIVATIVESC*** DO 10 I=1,9 VN(I)=0.0 !VN-SHAPE FUNCTION DO 10 J=1,3 VDN(J,I)=0.0 !VDN-LOCAL DERIVATIVE OF VN10 CONTINUEC***C******************SET FUNCTION VALUE FOR QUADRANDGLE ELEMENT OF 4 NODESC*** VN(1)=(1-X)*(1-Y)/4 VN(2)=(1+X)*(1-Y)/4 VN(3)=(1+X)*(1+Y)/4 VN(4)=(1-X)*(1+Y)/4 VDN(1,1)=-(1-Y)/4 VDN(1,2)=(1-Y)/4 VDN(1,3)=(1+Y)/4 VDN(1,4)=-(1+Y)/4 VDN(2,1)=-(1-X)/4 VDN(2,2)=-(1+X)/4 VDN(2,3)=(1+X)/4 VDN(2,4)=(1-X)/4C***C*****************SET FUNCTION VALUE FOR QUADRANGLE ELEMENT OF THE 5-8 NODESC*** IF(NODE.EQ.8) THEN VN(5)=(1-X*X)*(1-Y)/2 VN(6)=(1-Y*Y )*(1+X)/2 VN(7)=(1-X*X)*(1+Y)/2 VN(8)=(1-Y*Y)*(1-X)/2 VDN(1,5)=(-2*X)*(1-Y)/2 VDN(1,6)=(1-Y*Y)*(+1)/2 VDN(1,7)=(-2*X)*(1+Y)/2 VDN(1,8)=(1-Y*Y)*(-1)/2 VDN(2,5)=(1-X*X)*(-1)/2 VDN(2,6)=(-2*Y)*(1+X)/2 VDN(2,7)=(1-X*X)*(+1)/2 VDN(2,8)=(-2*Y)*(1-X)/2 DO 30 I=1,4 IF(IEL(4+I).EQ.0) VN(4+I)=0.0 !IEL-NODE CODE IN A ELEMENT IF(IEL(4+I).EQ.0) VDN(1,4+I)=0.0 IF(IEL(4+I).EQ.0) VDN(2,4+I)=0.030 CONTINUE VN(1)=VN(1)-(VN(5)+VN(8))/2 VN(2)=VN(2)-(VN(5)+VN(6))/2 VN(3)=VN(3)-(VN(6)+VN(7))/2 VN(4)=VN(4)-(VN(7)+VN(8))/2 DO 40 I=1,2 VDN(I,1)=VDN(I,1)-(VDN(I,5)+VDN(I,8))/2 VDN(I,2)=VDN(I,2)-(VDN(I,5)+VDN(I,6))/2 VDN(I,3)=VDN(I,3)-(VDN(I,6)+VDN(I,7))/2 VDN(I,4)=VDN(I,4)-(VDN(I,7)+VDN(I,8))/240 CONTINUE ENDIF RETURN ENDC=========================SUB:3-1-3================================================================= SUBROUTINE SHAPE_QUADRANGLE_9(NODE,X,Y,IEL,VN,VDN)C***********************************************************************************C COMPUTE SHAPE FUNCTION OF QUADRNGLE ELEMENT ON INTEGRATION POINTC************************************************************************************ IMPLICIT REAL*8(A-H,O-Z) DIMENSION IEL(NODE),VN(9),VDN(3,9) DIMENSION IX(9),IY(9),VL0X(3),VL1X(3),VL0Y(3),VL1Y(3) DATA IX/1,-1,-1,1,0,-1,0,1,0/ DATA IY/1,1,-1,-1,1,0,-1,0,0/C***C*****************CLEAR THE MEMORY FOR SHAPE FUNCTION AND ITS DERIVATIVESC*** DO 10 I=1,9 VN(I)=0.0 !VN-SHAPE FUNCTION DO 10 J=1,3 VDN(J,I)=0.0 !VDN-LOCAL DERIVATIVE OF VN10 CONTINUEC***C******************SET THE VALUE OF LAGRANGEFUNCTION AND ITS DERIVATIVESC*** VL0X(1)=0.5*X*(X-1) VL0X(2)=1-X*X VL0X(3)=0.5*X*(X+1) VL1X(1)=X-0.5 VL1X(2)=-2*X VL1X(3)=X+0.5 VL0Y(1)=0.5*Y*(Y-1) VL0Y(2)=1-Y*Y VL0Y(3)=0.5*Y*(Y+1) VL1Y(1)=Y-0.5 VL1Y(2)=-2*Y VL1Y(3)=Y+0.5C***C*******************SET SHAPE FUNCTION FOR QUADRANGLE ELEMENT OF 9 NODESC*** DO 20 I=1,9 IIX=IX(I)+2 IIY=IY(I)+2 VN(I)=VL0X(IIX)*VL0Y(IIY) VDN(1,I)=VL1X(IIX)*VL0Y(IIY) VDN(2,I)=VL0X(IIX)*VL1Y(IIY)20 CONTINUE IF(IEL(1).EQ.0) RETURN ENDC============================SUB:3-1-4================================================ SUBROUTINE HAMMER_INTEGRATION(INTX,I,X,Y,Z,WX)C***************************************************************************C GET THE INFORMATION OF HAMMER INTEGRATION POINTC*************************************************************************** IMPLICIT REAL*8(A-H,O-Z) DIMENSION HXY(3,5),WH(5),INDEX1(7)C !HXY-COORDINATE OF HAMMER INTEGRATION POINTC !WH-WEIGHT OF HAMMER INTEGRATION POINT INDEX1(1)=1 INDEX1(3)=2 INDEX1(4)=3 INDEX1(7)=4C***C***********INTEGRATION CONSTANTS OF ONE POINTC*** HXY(1,1)=1.0/3.0 HXY(2,1)=1.0/3.0 HXY(3,1)=1.0/3.0 WH(1)=1.0C***C***********INTEGRATION CONSTANTS OF 3 POINTSC*** HXY(1,2)=2.0/3.0 HXY(2,2)=1.0/6.0 HXY(3,2)=1.0/6.0 WH(2)=1.0/3.0C***C**********INTEGRATION CONSTANTS OF 4 POINTSC*** HXY(1,3)=0.6 HXY(2,3)=0.2 HXY(3,3)=0.2 WH(3)=25.0/48.0C***C**********INTEGRATION CONSTANTS OF 7 POINTS C*** A1=0.0597158717 B1=0.4701420641 A2=0.7974269853 B2=0.1012865073 HXY(1,4)=A1 HXY(2,4)=B1 HXY(3,4)=B1 WH(4)=0.1323941527 HXY(1,5)=A2 HXY(2,5)=B2 HXY(3,5)=B2 WH(5)=0.1259391805C***C**********GET PARAMETERS OF INTEGRATION POINTSC*** X=HXY(I+0-(I-1)/3*3,INDEX1(INTX)) Y=HXY(I+2-(I+1)/3*3,INDEX1(INTX)) Z=HXY(I+1-(I+0)/3*3,INDEX1(INTX)) WX=WH(INDEX1(INTX))/2.0CCCCC IF(INTX.EQ.7.AND.I.GE.4) THEN J=I-3 X=HXY(J+0-(J-1)/3*3,5) Y=HXY(J+2-(J+1)/3*3,5) Z=HXY(J+1-(J+0)/3*3,5) WX=WH(5)/2.0 ENDIFCCCCC IF(INTX.EQ.4.AND.I.EQ.4) THEN X=HXY(I+0-(I-1)/3*3,1) Y=HXY(I+1-(I+0)/3*3,1) Z=HXY(I+2-(I+1)/3*3,1) WX=-27.0/48.0/2.0 ENDIFCCCCC IF(INTX.EQ.7.AND.I.EQ.7) THEN X=HXY(1,1) Y=HXY(2,1) Z=HXY(3,1) WX=0.9/8.0 ENDIF RETURN ENDC==========================SUB:3-1-5========================================== SUBROUTINE SHAPE_TRIANGLE(NODE,X,Y,Z,IEL,VN,VDN)C*********************************************************************C COMPUTE SHAPE FUNCTION OF TRIANGLE ELEMENT AT INTEGRATION POINTC********************************************************************************* IMPLICIT REAL*8(A-H,O-Z) DIMENSION IEL(NODE),VN(9),VDN(3,9)C***C******************CLEAR THE MEMORY FOR SHAPE FUNCTION AND ITS DERIVATIVESC*** DO 10 I=1,9 VN(I)=0.0 !VN-SHAPE FUNCTION DO 10 J=1,3 VDN(J,I)=0.0 !VDN-LOCAL DERIVATIVE OF VN10 CONTINUEC***C******************SET THE FUNCTION VALUE FOR 3-NODE TRIANGLE ELEMENTC*** VN(1)=X VN(2)=Y VN(3)=Z VDN(1,1)=1.0 VDN(1,2)=0.0 VDN(1,3)=0.0 VDN(2,1)=0.0 VDN(2,2)=1.0 VDN(2,3)=0.0 VDN(3,1)=0.0 VDN(3,2)=0.0 VDN(3,3)=1.0C***C******************SET THE FUNCTION VALUE FOR TRIANGLE ELEMENT OF 4-6NODESC*** IF(NODE.EQ.6) THEN VN(4)=4*X*Y VN(5)=4*Y*Z VN(6)=4*Z*X VDN(1,4)=4*Y VDN(1,5)=0.0 VDN(1,6)=4*Z VDN(2,4)=4*X VDN(2,5)=4*Z VDN(2,6)=0.0 VDN(3,4)=0.0 VDN(3,5)=4*Y VDN(3,6)=4*X DO 20 I=1,3 IF(IEL(3+I).EQ.0) VN(3+I)=0.0 IF(IEL(3+I).EQ.0) VDN(1,3+I)=0.0 IF(IEL(3+I).EQ.0) VDN(2,3+I)=0.0 IF(IEL(3+I).EQ.0) VDN(3,3+I)=0.020 CONTINUE VN(1)=VN(1)-(VN(4)+VN(6))/2 VN(2)=VN(2)-(VN(4)+VN(5))/2 VN(3)=VN(3)-(VN(5)+VN(6))/2 DO 30 I=1,3 VDN(I,1)=VDN(I,1)-(VDN(I,4)+VDN(I,6))/2 VDN(I,2)=VDN(I,2)-(VDN(I,4)+VDN(I,5))/2 VDN(I,3)=VDN(I,3)-(VDN(I,5)+VDN(I,6))/230 CONTINUE ENDIF DO 40 I=1,NODE VDN(1,I)=VDN(1,I)-VDN(3,I) VDN(2,I)=VDN(2,I)-VDN(3,I)40 CONTINUE RETURN ENDC==============================SUB:3-1-6================================================= SUBROUTINE ELEMENT_VD(MPROB,VD,VMAE)C************************************************************************************C**************FORM ELEMENT ELASTIC MATRIX[D] ACCORDING TO TYPE OF PROBLEMC*********************************************************************************** IMPLICIT REAL*8(A-H,O-Z) DIMENSION VD(5,5),VMAE(4) E=VMAE(1) !E-ELASTIC MODULUS V=VMAE(2) !V-POSSION'S RATIOC***C**************CLEAR MEMORY FOR THE MATRIX:[D]C*** DO 30 I=1,5 DO 30 J=1,5 VD(I,J)=0.0 !VD-ELASTICMATRIX30 CONTINUEC***C**************COMPUTE[D] MATRIX FOR PLANE STRESS OR PLANE STRAIN PROBLEMC*** IF(MPROB.EQ.1.OR.MPROB.EQ.2) THEN IF(MPROB.EQ.2) E=E/(1-V*V) IF(MPROB.EQ.2) V=V/(1-V) D0=E/(1-V*V) VD(1,1)=D0 !MPROB=1--PLANE STRESS VD(2,2)=D0 VD(3,3)=D0*(1-V)/2 !MPROB=2--PLANE STRAIN VD(1,2)=D0*V VD(2,1)=D0*V ENDIFC***C***************COMPUTE [D]MATRIX FOR AXISYMMETRIC PROBLEMC*** IF(MPROB.EQ.3) THEN D0=E*(1-V)/(1+V)/(1-2*V) VD(1,1)=D0 VD(2,2)=D0 VD(3,3)=D0*(1-2*V)/2/(1-V) VD(4,4)=D0 VD(2,1)=D0*V/(1-V) !MPROB=3-AXISYMMETRIC VD(1,2)=D0*V/(1-V) VD(4,1)=D0*V/(1-V) VD(1,4)=D0*V/(1-V) VD(4,2)=D0*V/(1-V) VD(2,4)=D0*V/(1-V) ENDIFC***C****************COMPUTE[D]MATRIX FOR MINDLIN PLATE PROBLEMC*** IF(MPROB.EQ.4) THEN TH=VMAE(4) D0=E*TH*TH*TH/12/(1-V*V) VD(1,1)=D0 VD(2,2)=D0 !MPROB=4-MINDLIN PLATE VD(3,3)=D0*(1-V)/2 VD(1,2)=D0*V VD(2,1)=D0*V VD(4,4)=E/2/(1+V)*TH/(6.0/5.0) VD(5,5)=E/2/(1+V)*TH/(6.0/5.0) ENDIF RETURN ENDC============= ==================SUB:3-1-7=============================================== SUBROUTINE ELEMENT_JACOBI(MND,VXY,VDN,SJ,VD0)C*************************************************************************C GET THE DETERMINANT OF JACOBI MATRIX AND CARTESIAN DERIVATIVESC******************************************************************************** IMPLICIT REAL*8(A-H,O-Z) DIMENSION VXY(MND,2),VDN(3,9),VD0(3,9) DIMENSION VJJ(2,2),VJ1(2,2)C******C*******************FORM JACOBI MATRIXC****** DO 11 II=1,2 DO 11 JJ=1,2 VJJ(II,JJ)=0.0 !VJJ-JACOBI MATRIX[J] DO 11 KK=1,MND VJJ(II,JJ)=VJJ(II,JJ)+VDN(II,KK)*VXY(KK,JJ)11 CONTINUEC***C********************FORM THE DETERMINANT OF JACOBI MATRIXC*** SJ=VJJ(1,1)*VJJ(2,2)-VJJ(1,2)*VJJ(2,1) !SJ-|J|C***C*********************COMPUTE THE INVERSE OF JACOBI MATRIXC*** VJ1(1,1)=+VJJ(2,2)/SJ VJ1(1,2)=-VJJ(1,2)/SJ !VJ1-INVERSE OF [J] VJ1(2,1)=-VJJ(2,1)/SJ VJ1(2,2)=+VJJ(1,1)/SJC***C*********************COMPUTE THE CARTESIAN DERIVATIVES OF SHAPE FUNCTIONC*** DO 51 II=1,2 DO 51 JJ=1,9 VD0(II,JJ)=0.0 !VD0-GLOBAL DERIVATIVE OF VN DO 51 KK=1,2 VD0(II,JJ)=VD0(II,JJ)+VJ1(II,KK)*VDN(KK,JJ)51 CONTINUE RETURN ENDC===========================SUB:3-1-8======================================= SUBROUTINE ELEMENT_VB(MPROB,MND,VXY,VN,VD0,VB,SR)C**********************************************************************C FORM ELEMENT STRAIN MATRIX[B]ACCODING TO TYPE OF PROBLEMC*********************************************************************** IMPLICIT REAL*8(A-H,O-Z) DIMENSION VXY(MND,2),VN(9),VD0(3,9),VB(5,27) !VXY-单元节点坐标 !VN-形函数 !VD0-整体形函数导数 !VB-应变矩阵C***C**********************CLEAR MEMORY FOR THE MATRIX[B]C*** DO 30 JJ=1,5 DO 30 II=1,27 VB(JJ,II)=0.030 CONTINUEC***C***********************COMPUTE[B] MATRIX FOR PLANE STRESS OR PLANE STRAIN PROBLEMC*** IF(MPROB.EQ.1.OR.MPROB.EQ.2) THEN SR=1.0 DO 40 II=1,9 VB(1,(II-1)*2+1)=VD0(1,II) !MPROB=1-PLANE STRESS VB(2,(II-1)*2+2)=VD0(2,II) VB(3,(II-1)*2+1)=VD0(2,II) !MPROB=2-PLANE STRAIN VB(3,(II-1)*2+2)=VD0(1,II) 40 CONTINUE ENDIFC***C************************COMPUTE[B]MATRIX FOR AXISYMMETRIX PROBLEMC*** IF(MPROB.EQ.3) THENSR=0.0DO 45 II=1,MNDSR=SR+VXY(II,1)*VN(II)45 CONTINUE !MPROB=3-AXISYMMETRICDO 50 II=1,9 VB(1,(II-1)*2+1)=VD0(1,II) VB(2,(II-1)*2+2)=VD0(2,II) VB(3,(II-1)*2+1)=VD0(2,II) VB(3,(II-1)*2+2)=VD0(1,II) IF(ABS(SR).LT.1E-20) VB(4,(II-1)*2+1)=0.0IF(ABS(SR).GE.1E-20) VB(4,(II-1)*2+1)=VN(II)/SR50 CONTINUESR=SR*2*3.14159265ENDIFC***C***********************COMPUTE[B]MATRIC FOR THE MINDLIN PLATE PROBLEMC*** IF(MPROB.EQ.4) THEN SR=1.0 !MPROB=4-MINDLIN PLATE DO 70 II=1,9 VB(1,(II-1)*3+1)=-VD0(1,II) VB(2,(II-1)*3+2)=-VD0(2,II) VB(3,(II-1)*3+1)=-VD0(2,II) VB(3,(II-1)*3+2)=-VD0(1,II) VB(4,(II-1)*3+1)=-VN(II) VB(4,(II-1)*3+3)=-VD0(1,II) VB(5,(II-1)*3+2)=-VN(II) VB(5,(II-1)*3+3)=-VD0(2,II)70 CONTINUE ENDIF RETURN ENDC=============================SUB:4============================================================C********************************************************************************************C SOLVE STATIC PROBLEM TO GET THE NODE DISPLACEMENTSC CALL SUBROUTINES :DECOMPOS AND BACKSUBSC*************************************************************************************************** SUBROUTINE STATIC_SOLVE(GKK,GP,GU) IMPLICIT REAL*8(A-H,O-Z) COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA COMMON/COM3/MND2,NUMPT2 COMMON/ELEM/NODE,INTX,INTY DIMENSION GKK(NUMPT2,MBAND),GP(NUMPT2),GU(NUMPT2) WRITE(*,101)101 FORMAT(/5X,'## SOLVE EQUATION [GKK]{GU}={GP}##')C***C*******************TRIANGLE DECOMPOSITION OF THE MATRIX:[GKK]C*** C***C********************SUBSTITUTE BACK TO GET THE VECTOR {GU}C*** CALL GXIAO(NUMPT2,MBAND,GKK,GP) DO 30 I=1,NUMPT2 GU(I)=GP(I)30 CONTINUEC***************************************************************C OUTPUT THE NODAL DISPLACEMENTC**************************************************************** WRITE(*,103)103 FORMAT(/6X,'#OUTPUT NODAL DISPLACEMENT TO FILE#') IF(MPROB.EQ.4) THEN WRITE(14,40) ELSE WRITE(14,39) ENDIF39 FORMAT(10X,'NODAL DISPLACEMENTS'/2X,'NO.OF NODES',10X, $ 'X-',15X,'Y-')40 FORMAT(10X,'NODAL DISPLANCEMENTS'/2X,'NO.OF NODES',5X, $ 'THETA-X',7X,'THETA-Y',10X,'W-Z') DO 41 I=1,NUMPT IF(NF.EQ.2) WRITE(14,42) I,(GU(NF*(I-1)+J),J=1,NF) IF(NF.EQ.3) WRITE(14,43) I,(GU(NF*(I-1)+J),J=1,NF) 41 CONTINUE42 FORMAT(2X,I5,4X,2E18.8)43 FORMAT(2X,I5,4X,3E16.8) CLOSE(14) RETURN ENDC=============================SUB:4-1============================== SUBROUTINE GXIAO(NUMPT2,MBAND,ARRAY1,ARRAY2)IMPLICIT REAL*8(A-H,O-Z)DIMENSION ARRAY1(NUMPT2,MBAND), ARRAY2(NUMPT2)DO 20 K=1,NUMPT2-1IF(NUMPT2.GT.K+MBAND-1)THENIM=K+MBAND-1ELSEIM=NUMPT2END IFDO 20 I=K+1,IMIL=I-K+1C=ARRAY1(K,IL)/ARRAY1(K,1)DO 10 J=1,MBAND-IL+1M=J+I-K10 ARRAY1(I,J)=ARRAY1(I,J)-C*ARRAY1(K,M)ARRAY2(I)=ARRAY2(I)-C*ARRAY2(K)20 CONTINUE ARRAY2(NUMPT2)=ARRAY2(NUMPT2)/ARRAY1(NUMPT2,1)DO 40 I=NUMPT2-1,1,-1IF(MBAND.GT.NUMPT2-I+1)THENJO=NUMPT2-I+1ELSEJO=MBANDEND IFDO 30 J=2,JOIH=J+I-130 ARRAY2(I)=ARRAY2(I)-ARRAY1(I,J)*ARRAY2(IH)40 ARRAY2(I)=ARRAY2(I)/ARRAY1(I,1)50 CONTINUERETURNENDC=============================SUB:4-2============================== SUBROUTINE DECOMPOS(NUMPT2,MBAND,ARRAY)C******************************************************************************************C TRIANGLE DECOMPOSITION OF MATRIXC********************************************************************************************* IMPLICIT REAL*8(A-H,O-Z)DIMENSION ARRAY(NUMPT2,MBAND)DO 14 M=1,NUMPT2-1MM=M+MBAND-1IF(MM.GT.NUMPT2)MM=NUMPT2DO 13 I=M+1,MMMI=I+MBAND-1IF(MI.GT.NUMPT2)MI=NUMPT2DO 10 J=I,MIARRAY(I,J-I+1)=ARRAY(I,J-I+1)-ARRAY(M,I-M+1)*ARRAY(M,J-M+1) $ /ARRAY(M,1)10 CONTINUE13 CONTINUE14 CONTINUE RETURNENDC=============================SUB:4-3=========================================================== SUBROUTINE BACKSUBS(NUMPT2,MBAND,ARRAY1,ARRAY2)C***********************************************************************************************C BACKSUBSTITUTION OF TRIANGLE DECOMPOSITIONC********************************************************************************************** IMPLICIT REAL*8(A-H,O-Z)DIMENSION ARRAY1(NUMPT2,MBAND),ARRAY2(NUMPT2)DO 11 M=1,NUMPT2-1MM=M+MBAND-1IF(MM.GT.NUMPT2)MM=NUMPT2DO 11 I=M+1,MMMI=I+MBAND-1IF(MI.GT.NUMPT2)MI=NUMPT2DO 11 J=1,MIARRAY2(I)=ARRAY2(I)-ARRAY1(M,I-M+1)*ARRAY2(M)/ARRAY1(M,1)11 CONTINUEDO 20 I=NUMPT2,1,-1MI=I+MBAND-1IF(MI.GT.NUMPT2)MI=NUMPT2DO 21 J=I+1,MIARRAY2(I)=(ARRAY2(I)-ARRAY1(I,J-I+1)*ARRAY2(J))/ARRAY1(I,1)21 CONTINUEARRAY2(I)=ARRAY2(I)/ARRAY1(I,1)20 CONTINUERETURNENDC===========================SUB:5============================================================== SUBROUTINE STRESS(VMATI,IELEM,VCOOD,GU,VXY,IEL,VU,SSS, $ VSS,IADD,SSN)C*******************************************************************************C COMPUTE STRESSES AT THE INTEGRATION POINTS AND NODES C CALL SUBROUTINE STRESS_MATRIXC*********************************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA COMMON/COM3/MND2,NUMPT2 COMMON/ELEM/NODE,INTX,INTY DIMENSION IELEM(NUMEL,MND+4),VCOOD(NUMPT,2),GU(NUMPT2), $ VXY(MND,2),IEL(MND),VU(MND2),VSS(4,NFSTR), $ SSS(9,NFSTR),IADD(NUMPT),VMATI(NMATI,4), $ SSN(NUMPT,NFSTR),VSNN(9,NFSTR),VMAE(4)C DO 20 I=1,NUMPT IADD(I)=0 DO 20 J=1,NFSTR SSN(I,J)=0.020 CONTINUE WRITE(*,101)101 FORMAT(/6X,'# OUTPUT ELEMENT STRESS TO FILE#') IF(MTYPE.EQ.0.OR.MTYPE.EQ.2) THEN IF(MPROB.EQ.1.OR.MPROB.EQ.2) WRITE(15,21) IF(MPROB.EQ.3) WRITE(15,22) IF(MPROB.EQ.4) WRITE(15,23) ENDIF21 FORMAT(10X,'STRESS ON EACH INTEGRATION POINT'/'ELEM-NO.',1X, $ 'INTEG-NO.',5X,'SINGMA-X',7X,'SIGMA-Y',9X,'TAO-XY')22 FORMAT(10X,'STRESS ON EACH INTEGRATION POINT'/'ELEM-NO.',1X, $ 'INTEG-NO.',3X,'SIGMA-X',8X,'SIGMA-Y',8X,'TAO-XY',10X, $ 'ROTATION')23 FORMAT(10X,'STRESS ON EACH INTEGRATION POINT'/'ELEM-NO.',1X, $ 'INTEG.',5X,'M-X',10X,'M-Y',9X,'M-XY',7X,'TAO-XZ', $ 7X,'TAO-YZ')C***C*********FORM ALL THE ELEMENT INFORMATION FROM THE INPUT DATA DO 320 IE=1,NUMELDO 30 I=I,MND IEL(I)=IELEM(IE,I+4) DO 30 J=1,2 VXY(I,J)=0.0 IF(IEL(I).GT.0) VXY(I,J)=VCOOD(IEL(I),J)30 CONTINUE NODE=IELEM(IE,1) !NODE-NUMBER OF NODEES IN THE ELEMENT INTX=IELEM(IE,3) !INTX-INTEGERATION POINTS IN X-DIRECEION INTY=IELEM(IE,4) !INTY-INTEGERATION POINTS IN Y-DIRECEION IF(NODE.EQ.3.OR.NODE.EQ.6) INTY=1 INTXY=INTX*INTY IMATI=IELEM(IE,2) DO 31 J=1,431 VMAE(J)=VMATI(IMATI,J) DO 40 I=1,MND DO 40 J=1,NF IF (IEL(I).GT.0) VU(NF*(I-1)+J)=GU(NF*(IELEM(IE,I+4)-1)+J)40 CONTINUEC***C******************COMPUTE STRESS AT THE INTEGRATION POINTSC*** CALL STRESS_MATRIX(IE,VXY,IEL,VMAE,VU,SSS,VSS)C***C******************OUTPUT STRESS AT THE INTEGRATION POINTC*** IF(MTYPE.EQ.0.OR.MTYPE.EQ.2)THEN DO 70 I=1,INTXY IF(MPROB.EQ.1.OR.MPROB.EQ.2) $ WRITE(15,71) IE,I,(SSS(I,J),J=1,NFSTR) IF(MPROB.EQ.3) WRITE(15,72) IE,I,(SSS(I,J),J=1,NFSTR) IF(MPROB.EQ.4) WRITE(15,73) IE,I,(SSS(I,J),J=1,NFSTR)70 CONTINUE ENDIF71 FORMAT(1X,I5,1X,I5,1X,3(2X,E15.6))72 FORMAT(1X,I5,1X,I5,1X,4(1X,E15.6))73 FORMAT(1X,I5,1X,I4,1X,5E13.6)C***C*****************COMPUTE NODAL STRESS FROM STRESS MATRIXBY ELEMENT SMOOTHINGC*** DO 83 J=1,NFSTR IF(NODE.EQ.3.OR.NODE.EQ.6) THEN A1=5.0/3.0 B1=1.0/3.0 VSNN(1,J)=VSS(1,J)*A1-VSS(2,J)*B1-VSS(3,J)*B1 VSNN(2,J)=-VSS(1,J)*B1+VSS(2,J)*A1-VSS(3,J)*B1 VSNN(3,J)=-VSS(1,J)*B1-VSS(2,J)*B1+VSS(3,J)*A1 VSNN(4,J)=(VSNN(1,J)+VSNN(2,J))/2.0 VSNN(5,J)=(VSNN(2,J)+VSNN(3,J))/2.0 VSNN(6,J)=(VSNN(3,J)+VSNN(1,J))/2.0 ELSE A=1.8660254 B=-0.5 C=0.1339746 VSNN(1,J)=VSS(1,J)*A+VSS(2,J)*B+VSS(3,J)*C+VSS(4,J)*B VSNN(2,J)=VSS(1,J)*B+VSS(2,J)*A+VSS(3,J)*B+VSS(4,J)*C VSNN(3,J)=VSS(1,J)*C+VSS(2,J)*B+VSS(3,J)*A+VSS(4,J)*B VSNN(4,J)=VSS(1,J)*B+VSS(2,J)*C+VSS(3,J)*B+VSS(4,J)*A VSNN(5,J)=(VSNN(1,J)+VSNN(2,J))/2 VSNN(6,J)=(VSNN(2,J)+VSNN(3,J))/2 VSNN(7,J)=(VSNN(3,J)+VSNN(4,J))/2 VSNN(8,J)=(VSNN(4,J)+VSNN(1,J))/2 VSNN(9,J)=(VSNN(5,J)+VSNN(7,J))/2 ENDIF83 CONTINUEC***C*****************COMPUTE THE AVERAGE VALUE OF NODAL STRESSESC***DO 89 I=1,MNDIEL(I)=IELEM(IE,I+4)IH=IEL(I) IF(IH.GT.0) THEN IADD(IH)=IADD(IH)+1 DO 90 J=1,NFSTR SSN(IH,J)=(SSN(IH,J)*(IADD(IH)-1)+VSNN(I,J))/IADD(IH) 90 CONTINUE ENDIF89 CONTINUE320 CONTINUEC***C******************OUTPUT STRESS AT THE NODESC*** IF(MPROB.EQ.1.OR.MPROB.EQ.2) WRITE(15,91) IF(MPROB.EQ.3) WRITE(15,92) IF(MPROB.EQ.4) WRITE(15.93)91 FORMAT(/10X,'STRESSES AT EACH NODE'/2X,'NODE-NO.',8X, $ 'SIGMA-X',10X,'SIGMA-Y',10X,'TAO-XY',10X,'ROTATION')92 FORMAT(/10X,'STRESSES AT EACH NODE'/2X,'NODE-NO.',8X, $ 'SIGMA-X',10X,'SIGMA-Y',10X,'TAO-XY',10X,'ROTATION')93 FORMAT(/10X,'STRESSES AT EACH NODES'/2X,'NODE-NO.',5X,'M-X', $ 10X,'M-Y',10X,'M-XY',8X,'TAO-XZ',9X,'TAO-YZ') DO 95 I=1,NUMPT IF(MPROB.EQ.1.OR.MPROB.EQ.2) $ WRITE(15,96) I,(SSN(I,J),J=1,NFSTR) IF(MPROB.EQ.3) WRITE(15,97) I,(SSN(I,J),J=1,NFSTR) IF(MPROB.EQ.4) WRITE(15,98) I,(SSN(I,J),J=1,NFSTR)95 CONTINUE96 FORMAT(2X,I5,2X,3(2X,E15.6))97 FORMAT(2X,I5,2X,4(2X,E15.6))98 FORMAT(1X,I5,1X,5E14.6) CLOSE(15) RETURN ENDC==============================SUB:5-1============================================= SUBROUTINE STRESS_MATRIX(IE,VXY,IEL,VMAE,VU,SSS,VSS,GU)C***********************************************************************************C FORM THE ELEMENT STRESS MATRIXES:[VSG]AND[VSN]C******************************************************************************* IMPLICIT REAL*8(A-H,O-Z) COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA COMMON/COM3/MND2,NUMPT2 COMMON/ELEM/NODE,INTX,INTY DIMENSION VXY(MND,2),IEL(MND),VMAE(4),VU(MND2),VSS(4,NFSTR), $ SSS(9,NFSTR),VSG(NFSTR,MND2),VSN(NFSTR,MND2),GU(NUMPT2) DIMENSION VN(9),VDN(3,9),VD0(3,9),VD(5,5),VB(5,27)C***C*******************FORM THE[D]MATRIX ACCORDING TO TYPE OF THE PROBLEMC*** CALL ELEMENT_VD(MPROB,VD,VMAE) !VD-ELASTIC MATRIXC************************************************************************************C LOOP OVER THE NUMERICAL INTEGRATION POINTS(INTX,INTY)C************************************************************************************* IF(NODE.EQ.3.OR.NODE.EQ.6) INTY=1 DO 32 I=1,INTX DO 32 J=1,INTY IF(NODE.EQ.4.OR.NODE.EQ.8) THEN CALL GAUSS_INTEGRATION(INTX,INTY,I,J,X,Y,WXY) CALL SHAPE_QUADRANGLE_8(NODE,X,Y,IEL,VN,VDN) ENDIF !VN-SHAPE FUNCTION IF(NODE.EQ.9) THEN !VU-单元节点位移 CALL GAUSS_INTEGRATION (INTX,INTY,I,J,X,Y,WXY) CALL SHAPE_QUADRANGLE_9(NODE,X,Y,IEL,VN,VDN) ENDIF !VDN-LOCAL DERIVATIVE OF[VN] IF(NODE.EQ.3.OR.NODE.EQ.6)THENCALL HAMMER_INTEGRATION(INTX,I,X,Y,Z,WXY)CALL SHAPE_TRIANGLE(NODE,X,Y,Z,IEL,VN,VDN)ENDIFCALL ELEMENT_JACOBI(MND,VXY,VDN,SJ,VD0) !SJ-|J| !VD0-GLOBAL DERIVATIVE OF[VN] CALL ELEMENT_VB(MPROB,MND,VXY,VN,VD0,VB,SR) !VB-ELEMENT STRAIN MATRIX !SR-PARAMETER OF INTEGRATIONC***C********************FORM ELEMENT STRESS MATRIX[VSG]AT INTEGRATION POINTSC*** ([VSG]=[D]*[B]) DO 20 II=1,NFSTR DO 20 JJ=1,MND2 VSG(II,JJ)=0.0 DO 22 KK=1,NFSTR VSG(II,JJ)=VSG(II,JJ)+VD(II,KK)*VB(KK,JJ)22 CONTINUE20 CONTINUE k=(I-1)*INTY+J DO 30 II=1,NFSTRZ=0.0 IF(INTX.EQ.3.AND.INTY.EQ.1) VSS(K,II)=0.0 IF(INTX.EQ.2.AND.INTY.EQ.2) VSS(K,II)=0.0DO 30 JJ=1,MND2 Z=Z+VSG(II,JJ)*VU(JJ)SSS(K,II)=Z IF(INTX.EQ.3.AND.INTY.EQ.1) VSS(K,II)=SSS(K,II) IF(INTX.EQ.2.AND.INTY.EQ.2) VSS(K,II)=SSS(K,II)30 CONTINUE32 CONTINUEC*********************************************************************************************C FORM ELEMENT STESS MATRIX[VSN] AT OPTIMAL STRESS POINTSC ( [VSN]=[D]*[B])C************************************************************************************************** IF(NODE.EQ.3.OR.NODE.EQ.6) THEN INTR=3 INTS=1 IF(INTX.EQ.3.AND.INTY.EQ.1) GOTO 46 ELSE INTR=2 INTS=2 IF(INTX.EQ.2.AND.INTY.EQ.2) GOTO 46 ENDIF DO 45 I=1,INTR DO 45 J=1,INTS IF(NODE.EQ.4.OR.NODE.EQ.8) THEN CALL GAUSS_INTEGRATION(INTR,INTS,I,J,X,Y,WXY) CALL SHAPE_QUADRANGLE_8(NODE,X,Y,IEL,VN,VDN) ENDIF IF(NODE.EQ.9) THEN CALL GAUSS_INTEGRATION(INTR,INTS,I,J,X,Y,WXY) CALL SHAPE_QUADRANGLE_9(NODE,X,Y,IEL,VN,VDN) ENDIF IF(NODE.EQ.3.OR.NODE.EQ.6) THEN CALL HAMMER_INTEGRATION(INTR,I,X,Y,Z,WXY) CALL SHAPE_TRIANGLE(NODE,X,Y,Z,IEL,VN,VDN) ENDIF CALL ELEMENT_JACOBI(MND,VXY,VDN,SJ,VD0) CALL ELEMENT_VB(MPROB,MND,VXY,VN,VD0,VB,SR) DO 42 II=1,NFSTR DO 42 JJ=1,MND2 VSN(II,JJ)=0.0 DO 42 KK=1,NFSTR VSN(II,JJ)=VSN(II,JJ)+VD(II,KK)*VB(KK,JJ)42 CONTINUE K=(I-1)*INTS+J DO 43 II=1,NFSTR VSS(K,II)=0.0 DO 43 JJ=1,MND2 VSS(K,II)=VSS(K,II)+VSN(II,JJ)*VU(JJ) 43 CONTINUE45 CONTINUE46 CONTINUE RETURN ENDC==============================SUB:6============================================ SUBROUTINE DYNAM(GMM,GKK,GP,DAMP,U0,V0,A,AW,B)C*******************************************************************************C SOLVE THE DYNAMIC PROBLEM BY CENTRAL DIFFERENCE METHODC OR NEWMARK METHODC CALL SUBROUTINES:DECOMPOS,BACKSUBS,CENTER OR NEWMARKC******************************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA COMMON/COM3/MND2,NUMPT2 COMMON/DYN/HUV,OMEGA,CC1,CC2,TT,DT,ALPHA,DELTA DIMENSION GMM(NUMPT2,MBAND),GKK(NUMPT2,MBAND),GP(NUMPT2), $ A(NUMPT2),U0(NUMPT2),V0(NUMPT2),AW(NUMPT2,MBAND), $ B(NUMPT2),DAMP(NUMPT2,MBAND)C***C*********************COMPUTE DAMPLING MATRIX(DAMP)C*** DO 50 I=1,NUMPT2 DO 50 J=1,MBAND DAMP(I,J)=CC1*GMM(I,J)+CC2*GKK(I,J)50 CONTINUEC***C**********************FORM INITIAL DISPLACEMENT(U0) AND VELOCITY(V0)C*** DO 60 I=1,NUMPT2 U0(I)=0.0 V0(I)=0.060 CONTINUE IF(HUV.EQ.1.OR.HUV.EQ.3) THEN READ(5,*) READ(5,*)(U0(I),I=1,NUMPT2) ENDIF IF(HUV.EQ.2.OR.HUV.EQ.3) THEN READ(5,*) READ(5,*)(V0(I),I=1,NUMPT2) ENDIF WRITE(6,*)'INITIAL DISPLACEMENTS-U0:' WRITE(6,62)(U0(I),I=1,NUMPT2) WRITE(6,*)'INITIAL VELOCITIES-V0:' WRITE(6,62)(V0(I),I=1,NUMPT2)62 FORMAT(2X,4E15.6) C**************************************************************************************C COMPUTE INITIAL ACCELERATION-A(NUMPT2)C*********************************************************************************** DO 71 I=1,NUMPT271 A(I)=GP(I) DO 70 I=1,NUMPT2 IK=I-MBAND+1 IF(IK.LT.1) IK=1 DO 70 J=IK,I A(I)=A(I)-GKK(J,I-J+1)*U0(J)-DAMP(I,J-I+1)*V0(J)70 CONTINUE DO 72 I=1,NUMPT2IK=I+MBAND-1IF(IK.GT.NUMPT2) IK=NUMPT2DO 72 J=I+1,IKA(I)=A(I)-GKK(I,J-I+1)*U0(J)-DAMP(I,J-1)*V0(J)72 CONTINUE DO 73 I=1,NUMPT2 DO 73 J=1,MBAND AW(I,J)=GMM(I,J)73 CONTINUE CALL DECOMPOS(NUMPT2,MBAND,AW) CALL BACKSUBS(NUMPT2,MBAND,AW,A) !^^^^^^^^^C***C*******************BY USING GENTRAL DIFFERENCE METHOD WHEN MSOLV=2C*** IF(MSOLV.EQ.2) CALL CENTER(GMM,GKK,DAMP,GP,U0,V0,A,AW,B) RETURN ENDC===========================SUB:6-1===================================================== SUBROUTINE CENTER(GMM,GKK,DAMP,GP,U0,V0,A,AW,B)C********************************************************************************C SOLVE DYNAMIC RESPONSE BY CENTRAL DIFFERENCE METHODC CALL SUBROUTINES:DECOMPOS,BACKSUBSC********************************************************************************* IMPLICIT REAL*8(A-H,O-Z) COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA COMMON/COM3/MND2,NUMPT2 COMMON/DYN/HUV,OMEGA,CC1,CC2,TT,DT,ALPHA,DELTA DIMENSION GMM(NUMPT2,MBAND),GKK(NUMPT2,MBAND), $ DAMP(NUMPT2,MBAND),GP(NUMPT2),U0(NUMPT2), $ V0(NUMPT2),AW(NUMPT2,MBAND),A(NUMPT2),B(NUMPT2) CLOSE(30,STATUS='DELETE') OPEN(30,FILE='OUT_CEN',STATUS='UNKNOWN') WRITE(30,*)'DYNAMIC RESPONSE RESULT BY CENTRAL DIFFERENCE' WRITE(30,*)'TOTAL TIME=',TT WRITE(*,101)101 FORMAT(/6X,'#SOLVE DYNAMIC RESPONSE BY CENTRAL DIFFERENCE') WRITE(*,102)102 FORMAT(/6X,'%OUTPUT DYNAMIC DISPLACEMENT IN FILE')C***C***********************INITIAL COMPUTATIONSC*** C0=1./(DT*DT) C1=0.5/DT C2=2.*C0 C3=1./C2 DO 30 I=1,NUMPT2 A(I)=U0(I)-DT*V0(I)+C3*A(I)30 CONTINUE DO 40 I=1,NUMPT2 B(I)=GP(I) DO 40 J=1,MBAND AW(I,J)=GMM(I,J)40 GMM(I,J)=C0*GMM(I,J)+C1*DAMP(I,J)C***C************************COMPUTATIONS FOR EACH TIME STEPC*** CALL DECOMPOS(NUMPT2,MBAND,GMM) DO 300 Y=0,TT,DT IF(OMEGA.GT.0.0) AP=SIN(OMEGA*Y) IF(OMEGA.LT.0.0) AP=COS(OMEGA*Y) DO 41 I=1,NUMPT2 GP(I)=B(I) IF(OMEGA.NE.0.0) GP(I)=GP(I)*AP41 CONTINUE DO 50 I=1,NUMPT2 IK=I-MBAND+1 IF(IK.LT.1) IK=1 DO 50 J=IK,I GP(I)=GP(I)-(GKK(J,I-J+1)-C2*AW(J,I-J+1))*U0(J)- $ (C0*AW(J,I-J+1)-C1*DAMP(J,I-J+1))*A(J)50 CONTINUE DO 60 I=1,NUMPT2 IK=I+MBAND-1 IF(IK.GT.NUMPT2) IK=NUMPT2 DO 60 J=I+1,IK GP(I)=GP(I)-(GKK(I,J-I+1)-C2*AW(I,J-I+1))*U0(J)- $ (C0*AW(I,J-I+1)-C1*DAMP(I,J-I+1))*A(J)60 CONTINUE CALL BACKSUBS(NUMPT2,MBAND,GMM,GP) NSTEP=Y/DT+1 WRITE(30,61) NSTEP,DT,Y+DT IF(MPROB.EQ.4) THEN WRITE(30,73) WRITE(30,74)(GP(I),I=1,NUMPT2) ELSE WRITE(30,71) WRITE(30,72)(GP(I),I=1,NUMPT2) ENDIF DO 70 I=1,NUMPT2 A(I)=U0(I) U0(I)=GP(I)70 CONTINUE300 CONTINUE61 FORMAT(2X,'NO.OF STEP=',I5,3X,'STEP LENGTH=',E15.6,3X, $ 'AT TIME=',E15.6)71 FORMAT(2X,'DISPLACEMENT:'/2(13X,'X-',13X,'Y-'))72 FORMAT(4X,4E16.8)73 FORMAT(2X,'DISPLACEMENT:'/11X,'THETA-X',11X,'THETA-Y', $ 12X,'W-Z')74 FORMAT(6X,3E16.8) CLOSE(30) RETURN ENDC=============================SUB:6-2=============================================== SUBROUTINE NEWMARK(GMM,GKK,DAMP,GP,U0,V0,A,AW,B)C***************************************************************************C SOLVE DYNAMIC RESPONSE BY NEWMARK METHODC CALL SUBROUTINES:DECOMPOS,BACKSUBSC*************************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA COMMON/COM3/MND2,NUMPT2 COMMON/DYN/HUV,OMEGA,CC1,CC2,TT,DT,ALPHA,DELTA DIMENSION GMM(NUMPT2,MBAND),GKK(NUMPT2,MBAND), $ AW(NUMPT2,MBAND),DAMP(NUMPT2,MBAND),GP(NUMPT2), $ U0(NUMPT2),V0(NUMPT2),A(NUMPT2),B(NUMPT2) CLOSE(31,STATUS='DELETE') OPEN(31,FILE='OUT_NMK',STATUS='UNKNOWN') WRITE(31,*)'DYNAMIC RESPONSE RESULT BY NEWMARK METHOD' WRITE(31,*)'TOTAL TIME=',TT WRITE(*,101)101 FORMAT(/6X,'# SOLVE DYNAMIC RESPONSE BY NEWMARK METHOD##') WRITE(*,102)102 FORMAT(/6X,'% OUTPUT DYNAMIC DISPLACEMENT IN FILE')C***C*****************************INITIAL COMPUTATIONSC*** C0=1.0/(ALPHA*DT*DT) C1=DELTA/(ALPHA*DT) C2=1.0/(ALPHA*DT) C3=0.5/ALPHA-1.0 C4=DELTA/ALPHA-1.0 C5=DT/2.0*(DELTA/ALPHA-2.0) C6=DT*(1.0-DELTA) C7=DELTA*DTC************************************************************************************C COMPUTE K'=K+C0*M+C1*CC************************************************************************************ DO 40 I=1,NUMPT2 B(I)=GP(I) DO 40 J=1,MBAND AW(I,J)=GKK(I,J)40 GKK(I,J)=GKK(I,J)+C0*GMM(I,J)+C1*DAMP(I,J)C***********************************************************************************C TRIANGLE DECOMPOSITION OF THE MATRIX:[GKK]C*********************************************************************************** CALL DECOMPOS(NUMPT2,MBAND,GKK)C***C*************************COMPUTATIONS FOR EACH TIME STEPC*** DO 300 Y=0,TT,DT IF(OMEGA.GT.0.0) AP=SIN(OMEGA*Y) IF(OMEGA.LT.0.0) AP=COS(OMEGA*Y) DO 41 I=1,NUMPT2 GP(I)=B(I) IF(OMEGA.NE.0.0) GP(I)=GP(I)*AP41 CONTINUE DO 50 I=1,NUMPT2 IK=I-MBAND+1 IF(IK.LT.1) IK=1 DO 50 J=IK,I GP(I)=GP(I)+GMM(J,I-J+1)*(C0*U0(J)+C2*V0(J)+C3*A(J))+ $ DAMP(J,I-J+1)*(C1*U0(J)+C4*V0(J)+C5*A(J))50 CONTINUE DO 60 I=1,NUMPT2 IK=I+MBAND-1 IF(IK.GT.NUMPT2) IK=NUMPT2 DO 60 J=I+1,IK GP(I)=GP(I)+GMM(J,I-J+1)*(C0*U0(J)+C2*V0(J)+C3*A(J))+ $ DAMP(J,I-J+1)*(C1*U0(J)+C4*V0(J)+C5*A(J))60 CONTINUE CALL BACKSUBS(NUMPT2,MBAND,GKK,GP) NSTEP=Y/DT+1 WRITE(31,61) NSTEP,DT,Y+DT IF(MPROB.EQ.4) THEN WRITE(31,73) WRITE(31,74) (GP(I),I=1,NUMPT2) ELSE WRITE(31,71) WRITE(31,72) (GP(I),I=1,NUMPT2) ENDIF DO 70 I=1,NUMPT2 A(I)=C0*(GP(I)-U0(I))-C2*V0(I)-C3*A(I) V0(I)=V0(I)+C6*A(I)+C7*A(I) U0(I)=GP(I)70 CONTINUE300 CONTINUE61 FORMAT(2X,'NO.OF STEP=',I5,3X,'STEP LENGTH=',E15.6,3X, $ 'AT TIME=',E15.6)71 FORMAT(2X,'DIAPLACEMENT:'/2(13X,'X-',13X,'Y-'))72 FORMAT(4X,4E16.8)73 FORMAT(2X,'DISPLACEMENT:'/11X'THETA-X',11X,'THETA-Y', $ 12X,'W-Z')74 FORMAT(6X,3E16.8) CLOSE(31) RETURN ENDC===================================SUB:7============================================== SUBROUTINE EIGEN(GMM,GKK,AA,BB,GM,GK,V,W1,W2)C*************************************************************************************C SOLVE EIGENVALUE PROBLEM BY INVERSE OR SUBSPACE METHODSC CALL SUBROUTINE INVERSE OR SUBROUTINE SUBSPACEC*************************************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA COMMON/COM3/MND2,NUMPT2 DIMENSION GMM(NUMPT2,MBAND),GKK(NUMPT2,MBAND),AA(NUMPT2,NVA), $ BB(NUMPT2,NVA),GM(NVA,NVA),GK(NVA,NVA),V(NVA,NVA), $ W1(NVA),W2(NVA)C***C***********************BY USING INVERSE METHOD WHEN MSOLV=4C*** IF(MSOLV.EQ.4) CALL INVERSE(GMM,GKK,AA,BB)C***C***********************BY USING SUBSPACE METHOD WHEN MSOLV=5C*** IF(MSOLV.EQ.5) CALL SUBSPACE(GMM,GKK,AA,BB,GM,GK,V,W1,W2) RETURN ENDC===============================SUB:7-1================================================== SUBROUTINE INVERSE(GMM,GKK,AA,BB)C******************************************************************************************C SOLVE EIGENVALUE BY INVERSE ITERATION METHODC CALL SUBROUTINES:DECOMPOS,BACKSUBSC**************************************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA COMMON/COM3/MND2,NUMPT2 DIMENSION GKK(NUMPT2,MBAND),GMM(NUMPT2,MBAND),AA(NUMPT2,NVA), $ BB(NUMPT2) CLOSE(32,STATUS='DELETE') OPEN(32,FILE='OUT_VERS',STATUS='UNKNOWN') WRITE(32,*)'DYNAMIC CHARACTER RESUITS BY INVERSE METHOD' WRITE(*,101)101 FORMAT(/6X,'# SOLVE EIGENVALUE BY INVERSE METHOD#') WRITE(*,102)102 FORMAT(/6X,'% OUTPUT EIGENVALUE AND MODE IN FILE')C***C**************************FORM INITIAL ITERATION VECTORS A(I,J)C*** DO 5 I=1,NUMPT2 DO 5 J=1,NVA K=I+JAA(I,J)=RAN(K)5 CONTINUEC***C**************************TRIANGLE DECOMPOSITION OF STIFFNESS MATRIXC*** CALL DECOMPOS(NUMPT2,MBAND,GKK) DO 100 II=1,NVA W1=-1.0 GK=0.0 GM=0.0C***C*************************COMPUTE Y=MXC***333 CONTINUE DO 8 I=1,NUMPT28 BB(I)=0.0 DO 10 I=1,NUMPT2 IK=I-MBAND+1 IF(IK.LT.1) IK=1 DO 10 J=IK,I BB(I)=BB(I)+GMM(J,I-J+1)*AA(J,II)10 CONTINUE DO 12 I=1,NUMPT2 IK=I+MBAND-1 IF(IK.GT.NUMPT2) IK=NUMPT2 DO 12 J=I+1,IK BB(I)=BB(I)+GMM(I,J-I+1)*AA(J,II)12 CONTINUEC***C***********************SOLVE KX=YC***111 CONTINUE NSTEP=NSTEP+1 IF(NSTEP.GT.500) THEN WRITE(*,*)'NO.=',II,'STEP=',NSTEP WRITE(32,*)'NO.=',II,'STEP=',NSTEP RETURN ENDIF DO 20 I=1,NUMPT220 AA(I,II)=BB(I) CALL BACKSUBS(NUMPT2,MBAND,GKK,AA(I,II))C***C**********************COMPUTE W**2=K'/M'C*** DO 24 I=1,NUMPT2 GK=GK+AA(I,II)*BB(I) BB(I)=0.024 CONTINUE DO 25 I=1,NUMPT2 IK=I-MBAND+1 IF(IK.LT.1) IK=1 DO 25 J=IK,I BB(I)=BB(I)+GMM(I,J-I+1)*AA(J,II)25 CONTINUE DO 26 I=NUMPT2,1,-1 IK=I+MBAND-1 IF(IK.GT.NUMPT2) IK=NUMPT2 DO 26 J=I+1,IK BB(I)=BB(I)+GMM(I,J-I+1)*AA(J,II)26 CONTINUE DO 27 I=1,NUMPT227 GM=GM+AA(I,II)*BB(I) W2=GK/GM W1=ABS((W1-W2)/W2) IF(W1.GT.1.0E-6) THEN W1=W2C***C*********************GRAM-SCHMIDT ORTHOGONIGATIONC*** IF(II.EQ.1) GOTO 222 IF(II.NE.1) THEN DO 30 JJ=1,II-1 S=0.0 DO 40 I=1,NUMPT2 IK=I-MBAND+1 IF(IK.LT.1) IK=1 DO 40 J=IK,I S=S+AA(I,JJ)*GMM(J,I-J+1)*AA(J,II)40 CONTINUE DO 50 I=1,NUMPT2 IK=I+MBAND-1 IF(IK.GT.NUMPT2) IK=NUMPT2 DO 50 J=I+1,IK S=S+AA(I,JJ)*GMM(I,J-I+1)*AA(J,II)50 CONTINUE DO 51 K=1,NUMPT2 AA(K,II)=AA(K,II)-AA(K,JJ)*S51 CONTINUE30 CONTINUE S=0.0 DO 52 I=1,NUMPT2 IK=I-MBAND+1 IF(IK.LT.1) IK=1 DO 52 J=IK,I S=S+AA(I,II)*GMM(J,I-J+1)*AA(J,II)52 CONTINUE DO 54 I=1,NUMPT2 IK=I+MBAND-1 IF(IK.GT.NUMPT2) IK=NUMPT2 DO 54 J=I+1,IK S=S+AA(I,II)*GMM(I,J-I+1)*AA(J,II)54 CONTINUE DO 55 K=1,NUMPT2 AA(K,II)=AA(K,II)/S55 CONTINUE GK=0.0 GM=0.0 GOTO 333 ENDIF222 DO 80 I=1,NUMPT2 W2=SQRT(GM) BB(I)=BB(I)/W280 CONTINUE GK=0.0 GM=0.0 GOTO 111 ELSEC***C***************************COMPUTE FREQUENCY,PERIOD,EIGENVECTORC*** WW=SQRT(W2) PD=2*3.1415926/WW WRITE(32,93) II,NSTEP,WW,PD DO 90 I=1,NUMPT2 W2=SQRT(GM) AA(I,II)=AA(I,II)/W290 CONTINUE IF(MPROB.EQ.4) THEN WRITE(32,92) WRITE(32,95)(AA(I,II),I=1,NUMPT2) ELSE WRITE(32,91) WRITE(32,94)(AA(I,II),I=1,NUMPT2) ENDIF ENDIF NSTEP=0100 CONTINUE91 FORMAT(2X,'VIBRATION MODE:'/2(13X,'X-',13X,'Y-'))92 FORMAT(2X,'VIBRATION MODE:'/11X,'THETA-X',11X,'THETA-Y', $ 12X,'W-Z')93 FORMAT('***',2X,'NO.OF EIGENVALUE=',I5,4X,'ITERATION TIMES=', $ I5/5X,'FREQUENCE=',F16.4,4X,'PERION='E16.6)94 FORMAT(2X,4E16.8)95 FORMAT(4X,3E16.8) CLOSE(32) RETURN ENDC=======================================SUB:7-2===================================== SUBROUTINE SUBSPACE(GMM,GKK,AA,BB,GM,GK,V,W1,W2)C********************************************************************************C SOLVE EIGENVALUE BY SUBSPACE ITERATION METHODC CALL SUBROUTINES:DECOMPOS,BACKSUBS,JOCOBI,ARRANGEC************************************************************************************** IMPLICIT REAL*8(A-H,O-Z) COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA COMMON/COM3/MND2,NUMPT2 DIMENSION GMM(NUMPT2,MBAND),GKK(NUMPT2,MBAND),AA(NUMPT2,NVA), $ BB(NUMPT2,NVA),GM(NVA,NVA),GK(NVA,NVA),V(NVA,NVA), $ W2(NVA),GMK(NUMPT2),W1(NVA) CLOSE(33,STATUS='DELETE') OPEN(33,FILE='OUT_SUBS',STATUS='UNKNOWN') WRITE(*,101) WRITE(33,*)'RESULTS OF EIGENVALUES BY SUBSPACE METHOD:'101 FORMAT(/6X,'#SOLVE EIGENVALUE BY SUBSPACE METHOD#') WRITE(*,102)102 FORMAT(/6X,'% OUTPUT EIGENVALUE AND MODE INTO FILE')C***C***************************FORM INITIAL ITERATION VECTORSC*** DO 5 I=1,NUMPT2 GMK(I)=GMM(I,1)/GKK(I,1) AA(I,1)=1.05 CONTINUE DO 7 I=1,NUMPT2 DO 7 J=2,NVA7 AA(I,J)=0.0 DO 8 J=2,NVA QQ=0.0 DO 9 I=1,NUMPT2 IF(GMK(I).GT.QQ) THEN QQ=GMK(I) K=I ENDIF9 CONTINUE GMK(K)=0.0 AA(K,J)=1.08 CONTINUEC***C************************TRIANGLE DECOMPOSTITION OF STIFFNESS MATRIXC*** CALL DECOMPOS(NUMPT2,MBAND,GKK)C***C************************COMPUTE Y=MXC*** DO 10 K=1,NVA DO 11 I=1,NUMPT2 IK=I-MBAND+1 IF(IK.LT.1) IK=1 DO 11 J=IK,I BB(I,K)=BB(I,K)+GMM(J,I-J+1)*AA(J,K)11 CONTINUE DO 12 I=1,NUMPT2 IK=I+MBAND-1 IF(IK.GT.NUMPT2) IK=NUMPT2 DO 12 J=I+1,IK BB(I,K)=BB(I,K)+GMM(I,J-I+1)*AA(J,K)12 CONTINUEC***C**********************SOLVING EQUATION KX=YC*** DO 20 I=1,NUMPT220 AA(I,K)=BB(I,K)10 CONTINUE NSTEP=0111 CONTINUE NSTEP=NSTEP+1 DO 30 J=1,NUMPT2 DO 30 K=1,NVA30 BB(J,K)=AA(J,K) DO 32 K=1,NVA CALL BACKSUBS(NUMPT2,MBAND,GKK,AA(1,K))32 CONTINUEC***C***********************COMPUTE K1=XT*YC*** DO 40 I=1,NVA DO 40 J=1,NVA GK(I,J)=0.0 GK(I,J)=GK(I,J)+AA(K,I)*BB(K,J)40 CONTINUEC***C***********************COMPUTE Y1=M*X1C*** DO 45 K=1,NVA DO 50 I=1,NUMPT2 IK=I-MBAND+1 IF(IK.LT.1) IK=1 DO 50 J=IK,I BB(I,K)=BB(I,K)+GMM(J,I-J+1)*AA(J,K)50 CONTINUE DO 60 I=1,NUMPT2 IK=I+MBAND-1 IF(IK.GT.NUMPT2) IK=NUMPT2 DO 60 J=I+1,IK BB(I,K)=BB(I,K)+GMM(J,I-J+1)*AA(J,K)60 CONTINUE45 CONTINUEC***C*******************COMPUTE M1=XT*Y1C*** DO 70 I=1,NVA DO 70 J=1,NVA GM(I,J)=0.0 DO 70 K=1,NUMPT2 GM(I,J)=GM(I,J)+AA(K,I)*BB(K,J)70 CONTINUEC***********************************************************************************C SOLVE GENERAL EIGENVALUE PROBLEMC*********************************************************************************** CALL JACOBI(GK,GM,V,W2,NVA) CALL ARRANGE(W2,V,NVA)C***C**********************CHECK IF THE CONVERGENT CONDITION IS SATISFIEDC*** IF(NVA.EQ.NUMPT2) GOTO 222 DO 80 J=1,NVA IF(ABS((W2(J)-W1(J))/W2(J)).GT.1.E-4) THEN DO 82 I=1,NUMPT2 DO 82 K=1,NVA AA(I,K)=0.0 DO 82 M=1,NVA AA(I,K)=AA(I,K)+BB(I,M)*V(M,K)82 CONTINUE DO 85 K=1,NVA85 W1(K)=W2(K) GOTO 111 ENDIF80 CONTINUE222 CONTINUE DO 88 I=1,NUMPT2 DO 88 J=1,NVA BB(I,J)=0.0 DO 88 K=1,NVA BB(I,K)=BB(I,K)+AA(I,K)*V(K,J)88 CONTINUE DO 90 J=1,NVA WW=SQRT(W2(J)) PD=2*3.1415926/WW WRITE(33,91) J,NSTEP,WW,PD IF(MPROB.EQ.4) THEN WRITE(33,93) WRITE(33,95)(BB(I,J),I=1,NUMPT2) ELSE WRITE(33,92) WRITE(33,94)(BB(I,J),I=1,NUMPT2) ENDIF90 CONTINUE91 FORMAT('***',2X,'NO.OF EIGENVALUE=',I5,4X,'ITERATIN TIMES=', $ I5/5X,'FREQUENCY=',F16.4,4X,'PERIOD=',E16.6)92 FORMAT(2X,'VIBRATION MODE:'/2(13X,'X-',13X,'Y-'))93 FORMAT(2X,'VIBRATION MODE:'/11X,'THETA-X',11X,'THETA-Y', $ 12X,'W-Z')94 FORMAT(2X,4E16.8)95 FORMAT(5X,3E16.8) CLOSE(33) RETURN ENDC===================================SUB:7-2-1======================================== SUBROUTINE JACOBI(GK,GM,V,EIGV,N)C*******************************************************************************C SOLVE EIGENVALUE BY JACOBI METHODC********************************************************************************* IMPLICIT REAL*8(A-H,O-Z) DIMENSION GK(N,N),GM(N,N),V(N,N),EIGV(N),D(N) RTOL=1.E-12 NSMAX=15 DO 10 I=1,N D(I)=0.0 D(I)=GK(I,I)/GM(I,I)10 EIGV(I)=D(I) DO 30 I=1,N DO 20 J=1,N20 V(I,J)=0.030 V(I,I)=1.0 NSWEEP=0NN=N-140 NSWEEP=NSWEEP+1 EPS=(0.01**NSWEEP)**2 DO 110 J=1,NN JJ=J+1 DO 110 K=JJ,N EPTOLA=(GK(J,K)*GK(J,K))/(GK(J,J)*GK(K,K)) EPTOLB=(GM(J,K)*GM(J,K))/(GM(J,J)*GM(K,K)) IF((EPTOLA.LT.EPS).AND.(EPTOLB.LT.EPS)) GOTO 110 AKK=GK(K,K)*GM(J,K)-GM(K,K)*GK(J,K) AJJ=GK(J,J)*GM(J,K)-GM(J,J)*GK(J,K) AB=GK(J,J)*GM(K,K)-GK(K,K)*GM(J,J) CHECK=(AB*AB+4.0*AKK*AJJ)/4.0 IF(CHECK) 50,51,5150 STOP 22251 SQCH=SQRT(CHECK) D1=AB/2.+SQCH D2=AB/2.-SQCH DEN=D1 IF(ABS(D2).GT.ABS(D1)) DEN=D2 IF(DEN) 55,57,5557 CA=0.0 CG=-GK(J,K)/GK(K,K)GOTO 6055 CA=AKK/DENCG=-AJJ/DEN60 IF(N-2) 61,90,6161 JP1=J+1 JM1=J-1 KP1=K+1 KM1=K-1 IF(JM1-1) 63,62,6262 DO 68 I=1,JM1 AJ=GK(I,J) BJ=GM(I,J) AK=GK(I,K) BK=GM(I,K) GK(I,J)=AJ+CG*AK GM(I,J)=BJ+CG*BK GK(I,K)=AK+CA*AJ68 GM(I,K)=BK+CA*BJ63 IF(KP1-N) 64,64,6664 DO 65 I=KP1,N AJ=GK(J,I) BJ=GM(J,I) AK=GK(K,I) BK=GM(K,I) GK(J,I)=AJ+CG*AK GM(J,I)=BJ+CG*BK GK(K,I)=AK+CA*AJ65 GM(K,I)=BK+CA*BJ66 IF(JP1-KM1) 70,70,9070 DO 80 I=JP1,KM1 AJ=GK(J,I) BJ=GM(J,I) AK=GK(I,K) BK=GM(I,K) GK(J,I)=AJ+CG*AK GM(J,I)=BJ+CG*BK GK(I,K)=AK+CA*AJ80 GM(I,K)=BK+CA*BJ90 AK=GK(K,K) BK=GM(K,K) GK(K,K)=AK+2.*CA*GK(J,K)+CA*CA*GK(J,J) GM(K,K)=BK+2.*CA*GM(J,K)+CA*CA*GM(J,J) GK(J,J)=GK(J,J)+2.*CG*GK(J,K)+CG*CG*AK GM(J,J)=GM(J,J)+2.*CG*GM(J,K)+CG*CG*BK GK(J,K)=0.0 GM(J,K)=0.0 DO 91 I=1,N XJ=V(I,J) XK=V(I,K) V(I,J)=XJ+CG*XK91 V(I,K)=XK+CA*XJ110 CONTINUE DO 92 I=1,N IF(GK(I,I).GT.0.0.AND.GM(I,I).GT.0.0) GOTO 92 STOP 33392 EIGV(I)=GK(I,I)/GM(I,I) DO 93 I=1,N TOL=RTOL*D(I) DIF=ABS(EIGV(I)-D(I)) IF(DIF.GT.TOL) GOTO 9793 CONTINUE EPS=TOL**2 DO 94 J=1,NNJJ=J+1DO 94 K=JJ,N EPSA=(GK(J,K)*GK(J,K))/(GK(J,J)*GK(K,K)) EPSB=(GM(J,K)*GK(J,K))/(GM(J,J)*GM(K,K)) IF((EPSA.LT.EPS).AND.(EPSB.LT.EPS)) GOTO 94 GOTO 9794 CONTINUE DO 95 I=1,N DO 95 J=1,N GK(J,I)=GK(I,J)95 GM(J,I)=GM(I,J) DO 96 J=1,N BB=SQRT(GM(J,J)) DO 96 K=1,N96 V(K,J)=V(K,J)/BB RETURN97 DO 98 I=1,N98 D(I)=EIGV(I) IF(NSWEEP.LT.NSMAX) GOTO 40 GOTO 94RETURN ENDC===========================SUB:7-2-2==================================================== SUBROUTINE ARRANGE(W,V,N)C*************************************************************************************C ARRANGE EIGENVALUE INTO ORDERC************************************************************************************ IMPLICIT REAL*8(A-H,O-Z) DIMENSION W(N),V(N,N) DO 13 I=1,N-1 K=I P=W(I) DO 11 J=I+1,N IF(W(J).LT.P) THEN K=J P=W(J)ENDIF11 CONTINUE IF(K.NE.I) THEN W(K)=W(I) W(I)=P DO 12 J=1,N P=V(J,I) V(J,I)=V(J,K) V(J,K)=P12 CONTINUE ENDIF13 CONTINUE RETURN END