CCC
PROGRAM FEATP
C====================================================================
C THIS PROGRAM CAN SOLVE THE ELASTIC PROBLEM OF PLANE STRESS,
C PLANE STRAIN,AXISYMMETRIC SOLID AND MINDLIN PLATE
C====================================================================
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 FILENAME
DIMENSION IZ(300000),AR(15000000)
C !IZ-MAXIMUN SPACE FOR DYNAMIC INTEGE ARRAY
C !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)')FILENAME
OPEN(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 MODEL
C*****************************************************************************************************
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 ADDRESS
C****************************************************************************************************
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))
ENDIF
C****************************************************************************************************
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
END
C=========================================================================
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 SPACE
C***********************************************************************************************************
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,NUMPT2
C !MND-MAXIMAL NODE NUMBER IN ALL ELEMENTS
C !NUMEL-NEMBER OF GLOBAL ELEMENTS
C !NUMPT-NUMBER OF GLOBAL NODES
C !MBAND-HALF BANDWIDTH(INCLUDING DIAGONAL ELEMENT)
READ(5,*)MND,NUMEL,NUMPT,MBAND
READ(5,*)NFIX,NPC,MPROB,MSOLV
C !NFIX-NUMBER OF NODES SUBJECTER TO CONSTRIANT
C !NPC-NUBMER OF NODES SUBJECTED TO EQUIVALENT LOAD
C==========================================================
C MPROB=1-PLANE STRESS PROBLEM, MPROB=2-PLANE STRAIN PROBLEM
C MPROB=3-AXISYMMETRIC PROBLEM, MPROB=4-MINDLIN PLATE PROBLEM
C------------------------------------------------------------------------------------------------------------------------------
C MSOLV=1-STATIC ANALYSIS
C MSOLV=2-DYNAMIC RESPONSE ANALYSIS BY CENTRAL DIFFERENCE METHOD
C MSOLV=3-DYNAMIC RESPONSE ANALYSIS BY NEWMARK METHOD
C MSOLV=4-DYNAMIC CHARACTER ANALYSIS BY INVERSE ITERATION METHOD
C MSOLV=5-DYNAMIC CHARACTER ANALYSIS BY SUBSPACE ITERATION METHOD
C==================================================================
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,NVA
C !NMATI-KIND OF MATERIALS
C !GRAV-GRAVITY ACCELERATION
C !NVA-NUMBER OF EIGENVALUES
C----------------------------------------------------------------------------------------------------------------------------------------
C MTYPE-CONTROL KEY FOR OUTPUT RESULTS
C MTYPE=0-OUTPUT RESULTS INGLUDE GLOBAL MATRIXES AND STRESS AT
C GAUSS POINTS
C MTYPE=1-OUTPUT RESULTS INCLUDE GLOBAL MATRIXES
C MTYPE=2-OUTPUT RESULTS INCLUDE STRESSES AT GAUSS POINTS
C-------------------------------------------------------------------------------------------------------------------------------------------
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
ENDIF
C*************************************************************************************************************************
C ALLOCATE STORAGE SPACE FOR INPUT DATA
C**********************************************************************************************************************
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 LOAD
C****************************************************************************************************************************
C ALLOCATE STORAGE SPACE FOR ELEMENTAL MATRIXES AND GLOBAL MATRIXES
C******************************************************************************************************************************
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 POINTS
N11=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 ELEMENT
C******************************************************************************************************************************
C ALLOCAT STORAGE SPACE FOR DYNAMIC RESPONSE ANALYSIS
C*********************************************************************************************************************************
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
ENDIF
C*************************************************************************************************************************************
C ALLOCAT STORAGE SPACE FOR DYNAMIC CHARACTER ANALYSIS
C***************************************************************************************************************************************
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
END
C=======================SUB:2=========================================================
SUBROUTINE INPUT(IELEM,IFIXD,ILOAD,VMATI,VCOOD,VFIXD,VLOAD)
C*****************************************************************************************************************************
C INPUT ALL THE INFORMATION OF FE MODEL
C********************************************************************************************************************************
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 FILE<IN_DAT>TO 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,DELTA
C !HUV-INPUT CONTROL OF INITIAL DISPLACEMENT AND VELOCITY
C (HUV=1:INPUT DISPLACEMENT;HUV=2:INPUT VELOCITY;
C HUV=3:INPUT BOTH OF THEM)AD
C !OMEGA-FREQUENCE OF LO
C !CC1,CC2-CONSTANTS TO COMPUTE DAMPLING MATRIXS
C !TT-TOTAL TIME
C !DT-TIME STEP LENGTH
C !ALPHA,DELTA-CONSTANTS USED IN THE NEWMARK METHOD
ENDIF
C***************************************************************************************
C OUTPUT ABOVE INPUT DATA
C****************************************************************************
WRITE(*,102)
102 FORMAT(/5X,'%% OUTPUT INPUT-DATA TO<OUT_DAT> %%')
WRITE(6,10)
10 FORMAT(/8X,'MAXIMAL-ELEM-NODES,ELEMENTS,NODES,BANDWIDTH')
WRITE(6,11) MND,NUMEL,NUMPT,MBAND
11 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,GRAV
13 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,MTYPE
15 FORMAT(11X,I9,2(3X,I9))
WRITE(6,16)
16 FORMAT(/8X,'NODAL-FREEDOMS,STRESS-COMPONENTS')
WRITE(6,17) NF,NFSTR
17 FORMAT(11X,I9,3X,I10)
WRITE(6,18) ! OUTPUT NODAL COORDIATES
18 FORMAT(/8X,'NODAL COORDINATES'/9X'NO.',15X,'X-',11X,'Y-')
DO 19 I=1,NUMPT
19 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,NUMEL
22 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 INFORMATION
24 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,NFIX
25 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 NODES
27 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,NPC
28 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,NMATI
34 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 COMPUTATION
C**********
IF(MSOLV.EQ.2.OR.MSOLV.EQ.3)
$WRITE(6,40) HUV,OMEGA,CC1,CC2,TT,DT,ALPHA,DELTA
40 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 COMPUTATION
C*********
IF(MSOLV.EQ.4.OR.MSOLV.EQ.5) WRITE(6,45) NVA
45 FORMAT(/8X,'PARAMERER FOR DYNAMIC CHARACTER COMPUTATION:'/
$ 10X,'NUMBERS OF EIGENVALUES-NVA=',I3)
WRITE(6,301)
301 FORMAT(/'C************************************************',
$ '******************************',/)
RETURN
END
C=========================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-MATRIX
C*****************************************************************
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.0
10 CONTINUE
C*******
C*********LOOP OVER EACH ELEMENT
C*********
DO 320 IE=1,NUMEL !NUMEL-NUMBER OF ELEMENTS
C********
C*********FORM INFORMATION OF EACH ELEMENT FROM THE INPUT DATA
C*******
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 CONTINUE
C****************************************************************************
C COMPUTE ELEMENTAL MASS MATRIX[VMM] AND STIFFNESS MATRIX[VKK]
C FOR THE ELEMENTS WITH 3-6,4-8,9 NODES
C********************************************************************************
C
CALL ELEMENT_MATRIX(IE,VXY,IEL,VMM,VKK,VMAE)
C
C*************************************************************************
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 CONTINUE
20 CONTINUE
320 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 ACCELERATION
30 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,NUMPT
32 WRITE(6,74) I,(GP(NF*(I-1)+J),J=1,NF)
WRITE(6,301)
ENDIF
C*********************************************************************************
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)
ENDIF
40 CONTINUE
C************************************************************************************
C DRAW THE NODAL CONSTRAINT
C***********************************************************************************
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.
ENDIF
42 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.0
60 CONTINUE
ENDIF
50 CONTINUE
ENDIF
C****************************************************************************************
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<OUT_MKP>%%')
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 CONTINUE
71 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 CONTINUE
C*********
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 CONTINUE
74 FORMAT(5X,I6,4X,3E15.6)
WRITE(10,301)
301 FORMAT(/3X,'********************',
$ '***************************',/)
CLOSE(10)
ENDIF
RETURN
END
C=========================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 OR
C HAMMER_INTEGATION,SHAPE_QUADRANGLE_8 OR SHAPE_QUADRANGLE_9
C OR SHAPE_TRIANGLE,ELEMENT_VB,ELEMENT_JOCABI
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),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 PROBLEM
C********
CALL ELEMENT_VD(MPROB,VD,VMAE) !VD-ELASTIC MATRIX
C********
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.0
10 CONTINUE
C**********************************************************************************************
C COMPUTE SHAPE FUNCTION VN AND ITS LOCAL DERIVATIVE VDN AT EACH
C INTEGRATOPM POINTS IN ELEMENTS
C******************************************************************************************
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 FUNCTION
C !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)
ENDIF
C******
C*********************FORM THE JACOBI MAXTRIX[J] AT INTEGRATION POINTS
C******
CALL ELEMENT_JACOBI(MND,VXY,VDN,SJ,VD0)
C !SJ=|J|
C !VD0-GLOBAL DERIVATIVE OF VN
C !VDN-LOCAL DERIVATE OF VN
IF(SJ.LE.0.0) THEN
WRITE(*,99)IE,I,J,SJ
99 FORMAT(/3X,'***SJ.LE.0.0 IN ELEMENT=',I4,3X,'INTX=',I2,
$ 3X,'INTY=',I2,3X,'SJ=',E10.4)
ENDIF
C***
C**********FORM THE[B]MATRIX AT INTEGRATION POINTS
C***
CALL ELEMENT_VB(MPROB,MND,VXY,VN,VD0,VB,SR)
!VB-ELEMENT STRAIN MATRIX
!SR-PARAMETER OF INTEGRATION
C***
C*****************FORM THE ELEMENT STRESS MATRIX[S]=[D]*[B]
C***
DO 20 II=1,NFSTR
DO 20 JJ=1,MND2
VSG(II,JJ)=0.0
C !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 CONTINUE
C***
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*SR
C !VKK-ELEMENT STIFFNESS MATRIX
22 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*SR
C !VMM-ELEMENT MASS MATRIX
25 CONTINUE
302 CONTINUE
RETURN
END
C================================SUB3-1-1==================================
SUBROUTINE GAUSS_INTEGRATION(INTX,INTY,I,J,X,Y,WXY)
C*******************************************************************
C GET THE INFORMATION OF GAUSS INTEGRATION POINT
C**************************************************************8
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION GXY(3,3),WG(3,3)
C***
C********************GAUSS INTEGRATION CONSTANTS FOR 1,2 AND 3POINTS
C***
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.555555555555556
C***
C*******************GET PARAMETERS OF INTEGRATION POINT
C***
X=GXY(I,INTX)
Y=GXY(J,INTY)
WXY=WG(I,INTX)*WG(J,INTY)
RETURN
END
C=======================SUB:3-1-2======================================================
SUBROUTINE SHAPE_QUADRANGLE_8(NODE,X,Y,IEL,VN,VDN)
C***********************************************************************
C COMPUTE SHAPE FUNCTION OF QUADRANGLE ELEMENT AT INTEGRATION POINT
C*******************************************************************************
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 DERIVATIVES
C***
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 VN
10 CONTINUE
C***
C******************SET FUNCTION VALUE FOR QUADRANDGLE ELEMENT OF 4 NODES
C***
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)/4
C***
C*****************SET FUNCTION VALUE FOR QUADRANGLE ELEMENT OF THE 5-8 NODES
C***
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.0
30 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))/2
40 CONTINUE
ENDIF
RETURN
END
C=========================SUB:3-1-3=================================================================
SUBROUTINE SHAPE_QUADRANGLE_9(NODE,X,Y,IEL,VN,VDN)
C***********************************************************************************
C COMPUTE SHAPE FUNCTION OF QUADRNGLE ELEMENT ON INTEGRATION POINT
C************************************************************************************
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 DERIVATIVES
C***
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 VN
10 CONTINUE
C***
C******************SET THE VALUE OF LAGRANGEFUNCTION AND ITS DERIVATIVES
C***
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.5
C***
C*******************SET SHAPE FUNCTION FOR QUADRANGLE ELEMENT OF 9 NODES
C***
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
END
C============================SUB:3-1-4================================================
SUBROUTINE HAMMER_INTEGRATION(INTX,I,X,Y,Z,WX)
C***************************************************************************
C GET THE INFORMATION OF HAMMER INTEGRATION POINT
C***************************************************************************
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION HXY(3,5),WH(5),INDEX1(7)
C !HXY-COORDINATE OF HAMMER INTEGRATION POINT
C !WH-WEIGHT OF HAMMER INTEGRATION POINT
INDEX1(1)=1
INDEX1(3)=2
INDEX1(4)=3
INDEX1(7)=4
C***
C***********INTEGRATION CONSTANTS OF ONE POINT
C***
HXY(1,1)=1.0/3.0
HXY(2,1)=1.0/3.0
HXY(3,1)=1.0/3.0
WH(1)=1.0
C***
C***********INTEGRATION CONSTANTS OF 3 POINTS
C***
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.0
C***
C**********INTEGRATION CONSTANTS OF 4 POINTS
C***
HXY(1,3)=0.6
HXY(2,3)=0.2
HXY(3,3)=0.2
WH(3)=25.0/48.0
C***
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.1259391805
C***
C**********GET PARAMETERS OF INTEGRATION POINTS
C***
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.0
CCCCC
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
ENDIF
CCCCC
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
ENDIF
CCCCC
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
END
C==========================SUB:3-1-5==========================================
SUBROUTINE SHAPE_TRIANGLE(NODE,X,Y,Z,IEL,VN,VDN)
C*********************************************************************
C COMPUTE SHAPE FUNCTION OF TRIANGLE ELEMENT AT INTEGRATION POINT
C*********************************************************************************
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 DERIVATIVES
C***
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 VN
10 CONTINUE
C***
C******************SET THE FUNCTION VALUE FOR 3-NODE TRIANGLE ELEMENT
C***
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.0
C***
C******************SET THE FUNCTION VALUE FOR TRIANGLE ELEMENT OF 4-6NODES
C***
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.0
20 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))/2
30 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
END
C==============================SUB:3-1-6=================================================
SUBROUTINE ELEMENT_VD(MPROB,VD,VMAE)
C************************************************************************************
C**************FORM ELEMENT ELASTIC MATRIX[D] ACCORDING TO TYPE OF PROBLEM
C***********************************************************************************
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 RATIO
C***
C**************CLEAR MEMORY FOR THE MATRIX:[D]
C***
DO 30 I=1,5
DO 30 J=1,5
VD(I,J)=0.0 !VD-ELASTICMATRIX
30 CONTINUE
C***
C**************COMPUTE[D] MATRIX FOR PLANE STRESS OR PLANE STRAIN PROBLEM
C***
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
ENDIF
C***
C***************COMPUTE [D]MATRIX FOR AXISYMMETRIC PROBLEM
C***
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)
ENDIF
C***
C****************COMPUTE[D]MATRIX FOR MINDLIN PLATE PROBLEM
C***
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
END
C============= ==================SUB:3-1-7===============================================
SUBROUTINE ELEMENT_JACOBI(MND,VXY,VDN,SJ,VD0)
C*************************************************************************
C GET THE DETERMINANT OF JACOBI MATRIX AND CARTESIAN DERIVATIVES
C********************************************************************************
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 MATRIX
C******
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 CONTINUE
C***
C********************FORM THE DETERMINANT OF JACOBI MATRIX
C***
SJ=VJJ(1,1)*VJJ(2,2)-VJJ(1,2)*VJJ(2,1) !SJ-|J|
C***
C*********************COMPUTE THE INVERSE OF JACOBI MATRIX
C***
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)/SJ
C***
C*********************COMPUTE THE CARTESIAN DERIVATIVES OF SHAPE FUNCTION
C***
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
END
C===========================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 PROBLEM
C***********************************************************************
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.0
30 CONTINUE
C***
C***********************COMPUTE[B] MATRIX FOR PLANE STRESS OR PLANE STRAIN PROBLEM
C***
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
ENDIF
C***
C************************COMPUTE[B]MATRIX FOR AXISYMMETRIX PROBLEM
C***
IF(MPROB.EQ.3) THEN
SR=0.0
DO 45 II=1,MND
SR=SR+VXY(II,1)*VN(II)
45 CONTINUE !MPROB=3-AXISYMMETRIC
DO 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.0
IF(ABS(SR).GE.1E-20) VB(4,(II-1)*2+1)=VN(II)/SR
50 CONTINUE
SR=SR*2*3.14159265
ENDIF
C***
C***********************COMPUTE[B]MATRIC FOR THE MINDLIN PLATE PROBLEM
C***
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
END
C=============================SUB:4============================================================
C********************************************************************************************
C SOLVE STATIC PROBLEM TO GET THE NODE DISPLACEMENTS
C CALL SUBROUTINES :DECOMPOS AND BACKSUBS
C***************************************************************************************************
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 CONTINUE
C***************************************************************
C OUTPUT THE NODAL DISPLACEMENT
C****************************************************************
WRITE(*,103)
103 FORMAT(/6X,'#OUTPUT NODAL DISPLACEMENT TO FILE<OUT_DIS>#')
IF(MPROB.EQ.4) THEN
WRITE(14,40)
ELSE
WRITE(14,39)
ENDIF
39 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 CONTINUE
42 FORMAT(2X,I5,4X,2E18.8)
43 FORMAT(2X,I5,4X,3E16.8)
CLOSE(14)
RETURN
END
C=============================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-1
IF(NUMPT2.GT.K+MBAND-1)THEN
IM=K+MBAND-1
ELSE
IM=NUMPT2
END IF
DO 20 I=K+1,IM
IL=I-K+1
C=ARRAY1(K,IL)/ARRAY1(K,1)
DO 10 J=1,MBAND-IL+1
M=J+I-K
10 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,-1
IF(MBAND.GT.NUMPT2-I+1)THEN
JO=NUMPT2-I+1
ELSE
JO=MBAND
END IF
DO 30 J=2,JO
IH=J+I-1
30 ARRAY2(I)=ARRAY2(I)-ARRAY1(I,J)*ARRAY2(IH)
40 ARRAY2(I)=ARRAY2(I)/ARRAY1(I,1)
50 CONTINUE
RETURN
END
C=============================SUB:4-2==============================
SUBROUTINE DECOMPOS(NUMPT2,MBAND,ARRAY)
C******************************************************************************************
C TRIANGLE DECOMPOSITION OF MATRIX
C*********************************************************************************************
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION ARRAY(NUMPT2,MBAND)
DO 14 M=1,NUMPT2-1
MM=M+MBAND-1
IF(MM.GT.NUMPT2)MM=NUMPT2
DO 13 I=M+1,MM
MI=I+MBAND-1
IF(MI.GT.NUMPT2)MI=NUMPT2
DO 10 J=I,MI
ARRAY(I,J-I+1)=ARRAY(I,J-I+1)-ARRAY(M,I-M+1)*ARRAY(M,J-M+1)
$ /ARRAY(M,1)
10 CONTINUE
13 CONTINUE
14 CONTINUE
RETURN
END
C=============================SUB:4-3===========================================================
SUBROUTINE BACKSUBS(NUMPT2,MBAND,ARRAY1,ARRAY2)
C***********************************************************************************************
C BACKSUBSTITUTION OF TRIANGLE DECOMPOSITION
C**********************************************************************************************
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION ARRAY1(NUMPT2,MBAND),ARRAY2(NUMPT2)
DO 11 M=1,NUMPT2-1
MM=M+MBAND-1
IF(MM.GT.NUMPT2)MM=NUMPT2
DO 11 I=M+1,MM
MI=I+MBAND-1
IF(MI.GT.NUMPT2)MI=NUMPT2
DO 11 J=1,MI
ARRAY2(I)=ARRAY2(I)-ARRAY1(M,I-M+1)*ARRAY2(M)/ARRAY1(M,1)
11 CONTINUE
DO 20 I=NUMPT2,1,-1
MI=I+MBAND-1
IF(MI.GT.NUMPT2)MI=NUMPT2
DO 21 J=I+1,MI
ARRAY2(I)=(ARRAY2(I)-ARRAY1(I,J-I+1)*ARRAY2(J))/ARRAY1(I,1)
21 CONTINUE
ARRAY2(I)=ARRAY2(I)/ARRAY1(I,1)
20 CONTINUE
RETURN
END
C===========================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_MATRIX
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 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.0
20 CONTINUE
WRITE(*,101)
101 FORMAT(/6X,'# OUTPUT ELEMENT STRESS TO FILE<OUT_STR>#')
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)
ENDIF
21 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,NUMEL
DO 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,4
31 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 CONTINUE
C***
C******************COMPUTE STRESS AT THE INTEGRATION POINTS
C***
CALL STRESS_MATRIX(IE,VXY,IEL,VMAE,VU,SSS,VSS)
C***
C******************OUTPUT STRESS AT THE INTEGRATION POINT
C***
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
ENDIF
71 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 MATRIX<VSN>BY ELEMENT SMOOTHING
C***
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
ENDIF
83 CONTINUE
C***
C*****************COMPUTE THE AVERAGE VALUE OF NODAL STRESSES
C***
DO 89 I=1,MND
IEL(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
ENDIF
89 CONTINUE
320 CONTINUE
C***
C******************OUTPUT STRESS AT THE NODES
C***
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 CONTINUE
96 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
END
C==============================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 PROBLEM
C***
CALL ELEMENT_VD(MPROB,VD,VMAE) !VD-ELASTIC MATRIX
C************************************************************************************
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)THEN
CALL HAMMER_INTEGRATION(INTX,I,X,Y,Z,WXY)
CALL SHAPE_TRIANGLE(NODE,X,Y,Z,IEL,VN,VDN)
ENDIF
CALL 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 INTEGRATION
C***
C********************FORM ELEMENT STRESS MATRIX[VSG]AT INTEGRATION POINTS
C*** ([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 CONTINUE
20 CONTINUE
k=(I-1)*INTY+J
DO 30 II=1,NFSTR
Z=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.0
DO 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 CONTINUE
32 CONTINUE
C*********************************************************************************************
C FORM ELEMENT STESS MATRIX[VSN] AT OPTIMAL STRESS POINTS
C ( [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 CONTINUE
45 CONTINUE
46 CONTINUE
RETURN
END
C==============================SUB:6============================================
SUBROUTINE DYNAM(GMM,GKK,GP,DAMP,U0,V0,A,AW,B)
C*******************************************************************************
C SOLVE THE DYNAMIC PROBLEM BY CENTRAL DIFFERENCE METHOD
C OR NEWMARK METHOD
C CALL SUBROUTINES:DECOMPOS,BACKSUBS,CENTER OR NEWMARK
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/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 CONTINUE
C***
C**********************FORM INITIAL DISPLACEMENT(U0) AND VELOCITY(V0)
C***
DO 60 I=1,NUMPT2
U0(I)=0.0
V0(I)=0.0
60 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,NUMPT2
71 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,NUMPT2
IK=I+MBAND-1
IF(IK.GT.NUMPT2) IK=NUMPT2
DO 72 J=I+1,IK
A(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=2
C***
IF(MSOLV.EQ.2) CALL CENTER(GMM,GKK,DAMP,GP,U0,V0,A,AW,B)
RETURN
END
C===========================SUB:6-1=====================================================
SUBROUTINE CENTER(GMM,GKK,DAMP,GP,U0,V0,A,AW,B)
C********************************************************************************
C SOLVE DYNAMIC RESPONSE BY CENTRAL DIFFERENCE METHOD
C CALL SUBROUTINES:DECOMPOS,BACKSUBS
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/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<OUT_PUT>')
C***
C***********************INITIAL COMPUTATIONS
C***
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 STEP
C***
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)*AP
41 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 CONTINUE
300 CONTINUE
61 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
END
C=============================SUB:6-2===============================================
SUBROUTINE NEWMARK(GMM,GKK,DAMP,GP,U0,V0,A,AW,B)
C***************************************************************************
C SOLVE DYNAMIC RESPONSE BY NEWMARK METHOD
C CALL SUBROUTINES:DECOMPOS,BACKSUBS
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/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<OUT_NMK>')
C***
C*****************************INITIAL COMPUTATIONS
C***
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*DT
C************************************************************************************
C COMPUTE K'=K+C0*M+C1*C
C************************************************************************************
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 STEP
C***
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)*AP
41 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 CONTINUE
300 CONTINUE
61 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
END
C===================================SUB:7==============================================
SUBROUTINE EIGEN(GMM,GKK,AA,BB,GM,GK,V,W1,W2)
C*************************************************************************************
C SOLVE EIGENVALUE PROBLEM BY INVERSE OR SUBSPACE METHODS
C CALL SUBROUTINE INVERSE OR SUBROUTINE SUBSPACE
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
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=4
C***
IF(MSOLV.EQ.4) CALL INVERSE(GMM,GKK,AA,BB)
C***
C***********************BY USING SUBSPACE METHOD WHEN MSOLV=5
C***
IF(MSOLV.EQ.5) CALL SUBSPACE(GMM,GKK,AA,BB,GM,GK,V,W1,W2)
RETURN
END
C===============================SUB:7-1==================================================
SUBROUTINE INVERSE(GMM,GKK,AA,BB)
C******************************************************************************************
C SOLVE EIGENVALUE BY INVERSE ITERATION METHOD
C CALL SUBROUTINES:DECOMPOS,BACKSUBS
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
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<OUT_VERS>')
C***
C**************************FORM INITIAL ITERATION VECTORS A(I,J)
C***
DO 5 I=1,NUMPT2
DO 5 J=1,NVA
K=I+J
AA(I,J)=RAN(K)
5 CONTINUE
C***
C**************************TRIANGLE DECOMPOSITION OF STIFFNESS MATRIX
C***
CALL DECOMPOS(NUMPT2,MBAND,GKK)
DO 100 II=1,NVA
W1=-1.0
GK=0.0
GM=0.0
C***
C*************************COMPUTE Y=MX
C***
333 CONTINUE
DO 8 I=1,NUMPT2
8 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 CONTINUE
C***
C***********************SOLVE KX=Y
C***
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,NUMPT2
20 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.0
24 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,NUMPT2
27 GM=GM+AA(I,II)*BB(I)
W2=GK/GM
W1=ABS((W1-W2)/W2)
IF(W1.GT.1.0E-6) THEN
W1=W2
C***
C*********************GRAM-SCHMIDT ORTHOGONIGATION
C***
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)*S
51 CONTINUE
30 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)/S
55 CONTINUE
GK=0.0
GM=0.0
GOTO 333
ENDIF
222 DO 80 I=1,NUMPT2
W2=SQRT(GM)
BB(I)=BB(I)/W2
80 CONTINUE
GK=0.0
GM=0.0
GOTO 111
ELSE
C***
C***************************COMPUTE FREQUENCY,PERIOD,EIGENVECTOR
C***
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)/W2
90 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=0
100 CONTINUE
91 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
END
C=======================================SUB:7-2=====================================
SUBROUTINE SUBSPACE(GMM,GKK,AA,BB,GM,GK,V,W1,W2)
C********************************************************************************
C SOLVE EIGENVALUE BY SUBSPACE ITERATION METHOD
C CALL SUBROUTINES:DECOMPOS,BACKSUBS,JOCOBI,ARRANGE
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
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<OUT_SUBS>')
C***
C***************************FORM INITIAL ITERATION VECTORS
C***
DO 5 I=1,NUMPT2
GMK(I)=GMM(I,1)/GKK(I,1)
AA(I,1)=1.0
5 CONTINUE
DO 7 I=1,NUMPT2
DO 7 J=2,NVA
7 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
ENDIF
9 CONTINUE
GMK(K)=0.0
AA(K,J)=1.0
8 CONTINUE
C***
C************************TRIANGLE DECOMPOSTITION OF STIFFNESS MATRIX
C***
CALL DECOMPOS(NUMPT2,MBAND,GKK)
C***
C************************COMPUTE Y=MX
C***
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 CONTINUE
C***
C**********************SOLVING EQUATION KX=Y
C***
DO 20 I=1,NUMPT2
20 AA(I,K)=BB(I,K)
10 CONTINUE
NSTEP=0
111 CONTINUE
NSTEP=NSTEP+1
DO 30 J=1,NUMPT2
DO 30 K=1,NVA
30 BB(J,K)=AA(J,K)
DO 32 K=1,NVA
CALL BACKSUBS(NUMPT2,MBAND,GKK,AA(1,K))
32 CONTINUE
C***
C***********************COMPUTE K1=XT*Y
C***
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 CONTINUE
C***
C***********************COMPUTE Y1=M*X1
C***
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 CONTINUE
45 CONTINUE
C***
C*******************COMPUTE M1=XT*Y1
C***
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 CONTINUE
C***********************************************************************************
C SOLVE GENERAL EIGENVALUE PROBLEM
C***********************************************************************************
CALL JACOBI(GK,GM,V,W2,NVA)
CALL ARRANGE(W2,V,NVA)
C***
C**********************CHECK IF THE CONVERGENT CONDITION IS SATISFIED
C***
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,NVA
85 W1(K)=W2(K)
GOTO 111
ENDIF
80 CONTINUE
222 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)
ENDIF
90 CONTINUE
91 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
END
C===================================SUB:7-2-1========================================
SUBROUTINE JACOBI(GK,GM,V,EIGV,N)
C*******************************************************************************
C SOLVE EIGENVALUE BY JACOBI METHOD
C*********************************************************************************
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,N
20 V(I,J)=0.0
30 V(I,I)=1.0
NSWEEP=0
NN=N-1
40 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,51
50 STOP 222
51 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,55
57 CA=0.0
CG=-GK(J,K)/GK(K,K)
GOTO 60
55 CA=AKK/DEN
CG=-AJJ/DEN
60 IF(N-2) 61,90,61
61 JP1=J+1
JM1=J-1
KP1=K+1
KM1=K-1
IF(JM1-1) 63,62,62
62 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*AJ
68 GM(I,K)=BK+CA*BJ
63 IF(KP1-N) 64,64,66
64 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*AJ
65 GM(K,I)=BK+CA*BJ
66 IF(JP1-KM1) 70,70,90
70 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*AJ
80 GM(I,K)=BK+CA*BJ
90 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*XK
91 V(I,K)=XK+CA*XJ
110 CONTINUE
DO 92 I=1,N
IF(GK(I,I).GT.0.0.AND.GM(I,I).GT.0.0) GOTO 92
STOP 333
92 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 97
93 CONTINUE
EPS=TOL**2
DO 94 J=1,NN
JJ=J+1
DO 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 97
94 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,N
96 V(K,J)=V(K,J)/BB
RETURN
97 DO 98 I=1,N
98 D(I)=EIGV(I)
IF(NSWEEP.LT.NSMAX) GOTO 40
GOTO 94
RETURN
END
C===========================SUB:7-2-2====================================================
SUBROUTINE ARRANGE(W,V,N)
C*************************************************************************************
C ARRANGE EIGENVALUE INTO ORDER
C************************************************************************************
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)
ENDIF
11 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)=P
12 CONTINUE
ENDIF
13 CONTINUE
RETURN
END
文章标题
最新推荐文章于 2022-05-12 11:32:07 发布
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
