平面四边形四节点等参单元Fortran源程序

第第页平面四边形四节点等参单元Fortran源程序

C************************************************

C*FINITEELEMENTPROGRAM*C*FORTwoDIMENSIONALELASticityPROBLEM*C*WITH4NODE*C************************************************PROGRAMELASTICITYcharacter*32dat,cch

DIMENSIONSK(80000),COOR(2,300),AE(4,11),MEL(5,200),&WG(4),JR(2,300),MA(600),R(600),iew(30),STRE(3,200)COMMON/CMN1/NP,NE,NM,NRCOMMON/CMN2/N,MX,NH

COMMON/CMN3/RF(8),SKE(8,8),NN(8)

WRITE(*,*)'PLEASEENTERINPUTFILENAME'READ(*,'(A)')DAT

OPEN(4,FILE=dat,STATUS='OLD')

OPEN(7,FILE='OUT',STATUS='UNKNOWN')READ(4,*)NP,NE,NM,NR

WRITE(7,'(A,I6)')'NUMBEROFNODENP=',npWRITE(7,'(A,I6)')'NUMBEROFELEMENTNE=',neWRITE(7,'(A,I6)')'NUMBEROFMATERIALNM=',nmWRITE(7,'(A,I6)')'NUMBEROFsurportingNC=',NrCALLINPUT(JR,COOR,AE,MEL)CALLCBAND(MA,JR,MEL)DOI=1,NHSK(I)=0.0enddo

CALLSK0(SK,MEL,COOR,JR,MA,AE)doI=1,NR(I)=0.0enddopause'aaa'stop

READ(4,*)NCP,NBE,iz

WRITE(*,'(5i8)')NCP,NBE,izWRITE(7,'(5i8)')NCP,NBE,iz

IF(NCP.GT.0)CALLCONCR(NCP,R,JR)

IF(NBE.GT.0)CALLBODYR(NBE,R,MEL,COOR,JR,AE)IF(iz.GT.0)thendojj=1,iz

READ(4,*)Js,nse,(WG(I),I=1,4)read(4,*)(iew(m),m=1,nse)

CALLFACER(iew,NSE,R,MEL,COOR,JR,WG)enddoendif

CALLDECOP(SK,MA)CALLFOBA(SK,MA,R)CALLOUTDISP(NP,R,JR)

CALLSTRESS(COOR,MEL,JR,AE,R,STRE)

WRITE(7,'(A)')'PROGRAMSAFFHASBEENENDED'WRITE(*,'(A)')'PROGRAMSAFFHASBEENENDED'STOPcRETURNEND

C*********************************************SUBROUTINEINPUT(JR,COOR,AE,MEL)

DIMENSIONJR(2,*),COOR(2,*),AE(4,*),MEL(5,*)COMMON/CMN1/NP,NE,NM,NRCOMMON/CMN2/N,MX,NHDO70I=1,NPREAD(4,*)IP,X,YCOOR(1,IP)=XCOOR(2,IP)=Y70CONTINUEDO11J=1,NE

READ(4,*)NEE,NME,(MEL(I,NEE),I=1,4)MEL(5,NEE)=NME11CONTINUEDO10I=1,NPDO10J=1,210JR(J,I)=1

DO20I=1,NR

READ(4,*)IP,IX,IYJR(1,IP)=IXJR(2,IP)=IY20CONTINUEN=0

DO30I=1,NPDO30J=1,2

IF(JR(J,I))30,30,2525N=N+1JR(J,I)=N30CONTINUEDO55J=1,NM

READ(4,*)JJ,(AE(I,JJ),I=1,4)WRITE(*,910)JJ,(AE(I,JJ),I=1,4)55CONTINUE

910FORMAT(/20X,'MATERIALPROPERTIES'/(3X,I5,4(1x,E8.3)))RETURN

END

C**********************************************SUBROUTINECBAND(MA,JR,MEL)

DIMENSIONMA(*),JR(2,*),MEL(5,*),NN(8)COMMON/CMN1/NP,NE,NM,NRCOMMON/CMN2/N,MX,NHDO65I=1,N65MA(I)=0

DO90IE=1,NEDO75K=1,4IEK=MEL(K,IE)DO95M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)95CONTINUE75CONTINUEL=N

DO80I=1,2*4NNI=NN(I)

IF(NNI.EQ.0)GOTO80IF(NNI.LT.L)L=NNI80CONTINUEDO85M=1,2*4JP=NN(M)

IF(JP.EQ.0)GOTO85JPL=JP-L+1

IF(JPL.GT.MA(JP))MA(JP)=JPL85CONTINUE90CONTINUEMX=0MA(1)=1DO10I=2,N

IF(MA(I).GT.MX)MX=MA(I)MA(I)=MA(I)+MA(I-1)10CONTINUENH=MA(N)

WRITE(7,'(A,I8)')'TOTALDEGREESOFFREEDOMN=',NWRITE(7,'(A,I8)')'MAX-SEMI-BANDWIDTHMX=',MXWRITE(7,'(A,I8)')'TOTAL-STORAGENH=',NH500FORMAT(/5X,'FREEDOMN='

*,I5,3X,'SEMI-BANDWI.MX=',I5,3X,*'STORAGENH=',I7)RETURNEND

C**********************************************SUBROUTINESK0(SK,MEL,COOR,JR,MA,AE)

DIMENSIONSK(*),MEL(5,*),COOR(2,*),JR(2,*),MA(*),*AE(4,*),XYZ(2,4),iven(4)COMMON/CMN1/NP,NE,NM,NRCOMMON/CMN2/N,MX,NH

COMMON/CMN3/RF(8),SKE(8,8),NN(8)COMMON/CMN4/NEE,NMECOMMON/GAUSS/RSTG(3),H(3)H(1)=0.5555555555555560H(2)=0.8888888888888890H(3)=H(1)

RSTG(1)=-0.7745966692414830RSTG(2)=0.00

RSTG(3)=-RSTG(1)DO10IE=1,NENEE=IE

NME=MEL(5,IE)DO75K=1,4IEK=MEL(K,IE)iven(k)=IEKDO95M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)

95XYZ(M,K)=COOR(M,IEK)75CONTINUE

CALLSTIF(XYZ,AE,iven)DO60I=1,8DO60J=1,8II=NN(I)JJ=NN(J)

IF((JJ.EQ.0).OR.(II.LT.JJ))GOTO60JN=MA(II)-(II-JJ)

SK(JN)=SK(JN)+SKE(I,J)60CONTINUE70CONTINUE

write(7,1111)((ske(i,j),j=1,8),i=1,8)1111format(2x,8f12.2)10CONTINUERETURNEND

C*********************************************SUBROUTINESTIF(XYZ,AE,iven)

DIMENSIONAE(4,*),DNX(2,4),XYZ(2,*),iven(*),

*RJAC(2,2)

COMMON/CMN1/NP,NE,NM,NRCOMMON/CMN2/N,MX,NH

COMMON/CMN3/RF(8),SKE(8,8),NN(8)COMMON/CMN4/NEE,NMECOMMON/GAUSS/RSTG(3),H(3)DO40I=1,8RF(I)=0.00DO30J=1,8SKE(I,J)=0.0030CONTINUE40CONTINUEE=AE(1,NME)U=AE(2,NME)

GAMA=AE(3,NME)

D1=E*(1.00-U)/((1.00+U)*(1.00-2.00*U))D