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

平面四边形四节点等参单元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.55560

H

（2）=0.88890

H（3）=H

（1）

RSTG

（1）=-0.14830

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+DNIX*DNJY*DET*RH*SH

DYX=DYX+DNIY*DNJX*DET*RH*SH

DYY=DYY+DNIY*DNJY*DET*RH*SH

98CONTINUE

99CONTINUE

SKE（I1,J1）=DXX*D1+DYY*D3

SKE（I2,J2）=DYY*D1+DXX*D3

SKE（I1,J2）=DXY*D2+DYX*D3

SKE（I2,J1）=DYX*D2+DXY*D3

115CONTINUE

120CONTINUE

RETURN

END

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

SUBROUTINECONCR（NCP,R,JR）

DIMENSIONR（*）,JR（2,*）,XYZ

（2）

DO100I=1,NCP

READ（4,*）IP,PX,PY

XYZ

（1）=PX

XYZ

（2）=PY

DO95J=1,2

L=JR（J,IP）

IF（L.EQ.0）GOTO95

R（L）=R（L）+XYZ（J）

95CONTINUE

100CONTINUE

RETURN

END

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

SUBROUTINEBODYR（NBE,R,MEL,COOR,JR,AE）

DIMENSIONR（*）,MEL（5,*）,COOR（2,*）,JR（2,*）,

&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/CMN5/FUN（4）,PN（2,4）,XJAC（2,2）

COMMON/GAUSS/RSTG（3）,H（3）

H

（1）=1.0

H

（2）=1.0

RSTG

（1）=-0.96260

RSTG

（2）=-RSTG

（1）

DO10IE=1,NBE

DOI=1,8

RF（I）=0.00

ENDDO

cREAD（4,*）NEE

NEE=ie

NME=MEL（5,NEE）

GAMA=AE（3,NME）

DO75K=1,4

IEK=MEL（K,NEE）

iven（k）=iek

DO95M=1,2

JJ=2*（K-1）+M

NN（JJ）=JR（M,IEK）

95XYZ（M,K）=COOR（M,IEK）

75CONTINUE

DO99IS=1,2

S=RSTG（IS）

SH=H（IS）

DO98IR=1,2

RR=RSTG（IR）

RH=H（IR）

CALLFUN8（XYZ,RR,S,DET）

DO30I=1,4

J=2*I

RF（J）=RF（J）-FUN（I）*RH*SH*DET*GAMA

30CONTINUE

98CONTINUE

99CONTINUE

CALLASLOAD（R）

10CONTINUE

RETURN

END

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

SUBROUTINEFACER（iew,NSE,R,MEL,COOR,JR,WG）

DIMENSIONR（*）,MEL（5,*）,COOR（2,*）,JR（2,*）,wg（*）

*,XYZ（2,4）,iew（*）,PR

（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）

H

（1）=1.0

H

（2）=1.0

RSTG

（1）=-0.96260

RSTG

（2）=-RSTG

（1）

nwf=0

nnf=0

ir=wg

（1）+0.1

if（ir.eq.1）nwf=1

if（ir.eq.2）nnf=1

DO510IE=1,NSE

DOI=1,8

RF（I）=0.00

ENDDO

nee=iew（ie）

DO575K=1,4

IEK=MEL（K,NEE）

DO595M=1,2

JJ=2*（K-1）+M

NN（JJ）=JR（M,IEK）

595XYZ（M,K）=COOR（M,IEK）

575CONTINUE

IF（NWF.EQ.1）then

GAMA=WG

（2）

Z0=WG（3）

NSU=WG（4）+0.1

CALLSURLOD（NSU,XYZ,PR,Z0,GAMA,1）

endif

IF（NNF.EQ.1）then

q=WG

（2）

NSU=WG（4）+0.1

doj=1,2

PR（J）=q

enddo

CALLSURLOD（NSU,XYZ,PR,Z0,GAMA,2）

endif

CALLASLOAD（R）

510CONTINUE

RETURN

END

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

SUBROUTINESURLOD（NSU,XYZ,PR,Z0,GAMA,NSI）

DIMENSIONXYZ（2,*）,RST（3）,PR

（2）,KCRD（4）,KFACE（2,4）,

&FVAL（4）,NODES

（2）,FACT（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/CMN5/FUN（4）,PN（2,4）,XJAC（2,2）

COMMON/GAUSS/RSTG（3）,H（3）

DATAKCRD/1,1,2,2/

DATAKFACE/1,4,

*2,3,

*1,2,

*4,3/

DATAFVAL/-1.00,1.00,-1.00,1.00/

FACT

（1）=1.0

FACT

（2）=-1.0

FACT（3）=-1.0

FACT（4）=1.0

FACTNUS=FACT（NSU）

DOI=1,2

J=KFACE（I,NSU）

NODES（I）=J

ENDDO

IF（NSI.EQ.1）THEN

DOI=1,2

J=NODES（I）

Z=Z0-XYZ（2,J）

PR（I）=0.00

IF（Z.GT.0.00）PR（I）=Z*GAMA

ENDDO

ENDIF

ML=KCRD（NSU）

IF（ML.EQ.1）MM=2

IF（ML.EQ.2）MM=1

RST（ML）=FVAL（NSU）

DO70LX=1,2

RST（MM）=RSTG（LX）

CALLFUN8（XYZ,RST

（1）,RST

（2）,DET）

PXYZ=0.00

DO25I=1,2

J=NODES（I）

PXYZ=PXYZ+FUN（J）*PR（I）

25CONTINUE

A1=XJAC（MM,2）

A2=-XJAC（MM,1）

30DO60I=1,2

J=NODES（I）

K2=2*J

K1=K2-1

Q=PXYZ*FUN（J）*H（LX）*FACTNUS

RF（K1）=RF（K1）+Q*A1

RF（K2）=RF（K2）+Q*A2

60CONTINUE

70CONTINUE

RETURN

END

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

SUBROUTINEASLOAD（R）

DIMENSIONR（*）

COMMON/CMN1/NP,NE,NM,NR

COMMON/CMN3/RF（8）,SKE（8,8）,NN（8）

DO20I=1,8

L=NN（I）

IF（L.EQ.0）GOTO20

R（L）=R（L）+RF（I）

20CONTINUE

RETURN

END

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

SUBROUTINEDECOP（SK,MA）

DIMENSIONSK（*）,MA（*）

COMMON/CMN2/N,MX,NH

DO50I=2,N

L=I-MA（I）+MA（I-1）+1

K=I-1

L1=L+1

IF（L1.GT.K）GOTO30

DO20J=L1,K

IJ=MA（I）-I+J

M=J-MA（J）+MA（J-1）+1

IF（L.GT.M）M=L

MP=J-1

IF（M.GT.MP）GOTO20

DO10LP=M,MP

IP=MA（I）-I+LP

JP=MA（J）-J+LP

SK（IJ）=SK（IJ）-SK（IP）*SK（JP）

10CONTINUE

20CONTINUE

30IF（L.GT.K）GOTO50

DO40LP=L,K

IP=MA（I）-I+LP

LPP=MA（LP）

SK（IP）=SK（IP）/SK（LPP）

II=MA（I）

SK（II）=SK（II）-SK（IP）*SK（IP）*SK（LPP）

40CONTINUE

50CONTINUE

RETURN

END

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

SUBROUTINEFOBA（SK,MA,R）

DIMENSIONSK（*）,MA（*）,R（*）

COMMON/CMN2/N,MX,NH

DO10I=2,N

L=I-MA（I）+MA（I-1）+1

K=I-1

IF（L.GT.K）GOTO10

DO5LP=L,K

IP=MA（I）-I+LP

R（I）=R（I）-SK（IP）*R（LP）

5CONTINUE

10CONTINUE

DO20I=1,N

II=MA（I）

45R（I）=R（I）/SK（II）

20CONTINUE

DO30J1=2,N

I=2+N-J1

L=I-MA（I）+MA（I-1）+1

K=I-1

IF（L.GT.K）GOTO30

DO25J=L,K

IJ=MA（I）-I+J

55R（J）=R（J）-SK（IJ）*R（I）

25CONTINUE

30CONTINUE

RETURN

END

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

SUBROUTINESTRESS（COOR,MEL,JR,AE,R,STRE）

DIMENSIONXYZ（2,4）,DNX（2,4）,AE（4,*）,STRE（3,*）,

&COOR（2,*）,MEL（5,*）,JR（2,*）,RJAC（2,2）,SIG（