平面四边形四节点等参单元Fortran源程序.docx
《平面四边形四节点等参单元Fortran源程序.docx》由会员分享，可在线阅读，更多相关《平面四边形四节点等参单元Fortran源程序.docx（21页珍藏版）》请在冰豆网上搜索。
平面四边形四节点等参单元Fortran源程序
C************************************************
C*FINITEELEMENTPROGRAM*
C*FORTwoDIMENSIONALELASticityPROBLEM*
C*WITH4NODE*
C************************************************
PROGRAMELASTICITY
character*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,NR
COMMON/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）'）'NUMBEROFNODE---------------------NP=',np
WRITE（7,'（A,I6）'）'NUMBEROFELEMENT------------------NE=',ne
WRITE（7,'（A,I6）'）'NUMBEROFMATERIAL-----------------NM=',nm
WRITE（7,'（A,I6）'）'NUMBEROFsurporting---------------NC=',Nr
CALLINPUT（JR,COOR,AE,MEL）
CALLCBAND（MA,JR,MEL）
DOI=1,NH
SK（I）=0.0
enddo
CALLSK0（SK,MEL,COOR,JR,MA,AE）
doI=1,N
R（I）=0.0
enddo
pause'aaa'
stop
READ（4,*）NCP,NBE,iz
WRITE（*,'（5i8）'）NCP,NBE,iz
WRITE（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）then
dojj=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）
enddo
endif
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'
STOP
cRETURN
END
C*********************************************
SUBROUTINEINPUT（JR,COOR,AE,MEL）
DIMENSIONJR（2,*）,COOR（2,*）,AE（4,*）,MEL（5,*）
COMMON/CMN1/NP,NE,NM,NR
COMMON/CMN2/N,MX,NH
DO70I=1,NP
READ（4,*）IP,X,Y
COOR（1,IP）=X
COOR（2,IP）=Y
70CONTINUE
DO11J=1,NE
READ（4,*）NEE,NME,（MEL（I,NEE）,I=1,4）
MEL（5,NEE）=NME
11CONTINUE
DO10I=1,NP
DO10J=1,2
10JR（J,I）=1
DO20I=1,NR
READ（4,*）IP,IX,IY
JR（1,IP）=IX
JR（2,IP）=IY
20CONTINUE
N=0
DO30I=1,NP
DO30J=1,2
IF（JR（J,I））30,30,25
25N=N+1
JR（J,I）=N
30CONTINUE
DO55J=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,NR
COMMON/CMN2/N,MX,NH
DO65I=1,N
65MA（I）=0
DO90IE=1,NE
DO75K=1,4
IEK=MEL（K,IE）
DO95M=1,2
JJ=2*（K-1）+M
NN（JJ）=JR（M,IEK）
95CONTINUE
75CONTINUE
L=N
DO80I=1,2*4
NNI=NN（I）
IF（NNI.EQ.0）GOTO80
IF（NNI.LT.L）L=NNI
80CONTINUE
DO85M=1,2*4
JP=NN（M）
IF（JP.EQ.0）GOTO85
JPL=JP-L+1
IF（JPL.GT.MA（JP））MA（JP）=JPL
85CONTINUE
90CONTINUE
MX=0
MA
（1）=1
DO10I=2,N
IF（MA（I）.GT.MX）MX=MA（I）
MA（I）=MA（I）+MA（I-1）
10CONTINUE
NH=MA（N）
WRITE（7,'（A,I8）'）'TOTALDEGREESOFFREEDOM-----------N=',N
WRITE（7,'（A,I8）'）'MAX-SEMI-BANDWIDTH-----------------MX=',MX
WRITE（7,'（A,I8）'）'TOTAL-STORAGE----------------------NH=',NH
500FORMAT（/5X,'FREEDOMN='
*,I5,3X,'SEMI-BANDWI.MX=',I5,3X,
*'STORAGENH=',I7）
RETURN
END
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,NR
COMMON/CMN2/N,MX,NH
COMMON/CMN3/RF（8）,SKE（8,8）,NN（8）
COMMON/CMN4/NEE,NME
COMMON/GAUSS/RSTG（3）,H（3）
H
（1）=0.5555555555555560
H
（2）=0.8888888888888890
H（3）=H
（1）
RSTG
（1）=-0.7745966692414830
RSTG
（2）=0.00
RSTG（3）=-RSTG
（1）
DO10IE=1,NE
NEE=IE
NME=MEL（5,IE）
DO75K=1,4
IEK=MEL（K,IE）
iven（k）=IEK
DO95M=1,2
JJ=2*（K-1）+M
NN（JJ）=JR（M,IEK）
95XYZ（M,K）=COOR（M,IEK）
75CONTINUE
CALLSTIF（XYZ,AE,iven）
DO60I=1,8
DO60J=1,8
II=NN（I）
JJ=NN（J）
IF（（JJ.EQ.0）.OR.（II.LT.JJ））GOTO60
JN=MA（II）-（II-JJ）
SK（JN）=SK（JN）+SKE（I,J）
60CONTINUE
70CONTINUE
write（7,1111）（（ske（i,j）,j=1,8）,i=1,8）
1111format（2x,8f12.2）
10CONTINUE
RETURN
END
C*********************************************
SUBROUTINESTIF（XYZ,AE,iven）
DIMENSIONAE（4,*）,DNX（2,4）,XYZ（2,*）,iven（*）,
*RJAC（2,2）
COMMON/CMN1/NP,NE,NM,NR
COMMON/CMN2/N,MX,NH
COMMON/CMN3/RF（8）,SKE（8,8）,NN（8）
COMMON/CMN4/NEE,NME
COMMON/GAUSS/RSTG（3）,H（3）
DO40I=1,8
RF（I）=0.00
DO30J=1,8
SKE（I,J）=0.00
30CONTINUE
40CONTINUE
E=AE（1,NME）
U=AE（2,NME）
GAMA=AE（3,NME）
D1=E*（1.00-U）/（（1.00+U）*（1.00-2.00*U））
D2=E*U/（（1.00+U）*（1.00-2.00*U））
D3=E*0.50/（1.00+U）
DO120I=1,4
II=2*（I-1）
I1=II+1
I2=II+2
DO115J=1,4
JJ=2*（J-1）
J1=JJ+1
J2=JJ+2
DXX=0
DXY=0
DYX=0
DYY=0
DO99IS=1,3
S=RSTG（IS）
SH=H（IS）
DO98IR=1,3
R=RSTG（IR）
RH=H（IR）
CALLFDNX（XYZ,DNX,DET,R,S,RJAC,iven,NEE）
DNIX=DNX（1,I）
DNIY=DNX（2,I）
DNJX=DNX（1,J）
DNJY=DNX（2,J）
DXX=DXX+DNIX*DNJX*DET*RH*SH
DXY=DXY
升级会员 