有限元编程算例.docx
《有限元编程算例.docx》由会员分享，可在线阅读，更多相关《有限元编程算例.docx（18页珍藏版）》请在冰豆网上搜索。
有限元编程算例
有限元编程算例（Fortran）
本程序通过Fortran语言编写，程序在IntelParallelStudioXE2013withVS2013中成功运行，程序为《计算力学》（龙述尧等编）一书中的源程序，仅作研究学习使用，省去了敲写的麻烦。
源程序为：
!
Page149
COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE
COMMON/X2/JM（100,3）,NZC（50）,CJZ（100,2）,PJ（100,2）,B（3,6）,D（3,3）,S（3,6）,TKZ（200,20）,EKE（6,6）,P（200）
OPEN（5,FILE='DATAIN'）
!
OPEN（6,FILE='DATAOUT',STATUS='NEW'）
CALLDATA
IF（IND.EQ.0）GOTO10
EO=EO/（1.0-UN*UN）
UN=UN/（1.0-UN）
10CALLTOTSTI
CALLLOAD
CALLSUPPOR
CALLSOLVEQ
CALLSTRESS
PAUSE
!
STOP
END
SUBROUTINEDATA
COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE
COMMON/X2/JM（100,3）,NZC（50）,CJZ（100,2）,PJ（100,2）,B（3,6）,D（3,3）,S（3,6）,TKZ（200,20）,EKE（6,6）,P（200）
READ（5,*）NJ,NE,NZ,NDD,NPJ,IND
NJ2=NJ*2
NPJ1=NPJ+1
READ（5,*）EO,UN,GAMA,TE
READ（5,*）（（JM（I,J）,J=1,3）,I=1,NE）
READ（5,*）（（CJZ（I,J）,J=1,2）,I=1,NJ）
!
Page150
READ（5,*）（NZC（I）,I=1,NZ）
READ（5,*）（（PJ（I,J）,J=1,2）,I=1,NPJ1）
WRITE（6,10）（I,（CJZ（I,J）,J=1,2）,I=1,NJ）
10FORMAT（4X,2HNO,6X,1HX,6X,1HY/（I6,2X,F7.2,F7.2））
RETURN
END
SUBROUTINEELEST（MEO,IASK）
COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE
COMMON/X2/JM（100,3）,NZC（50）,CJZ（100,2）,PJ（100,2）,B（3,6）,D（3,3）,S（3,6）,TKZ（200,20）,EKE（6,6）,P（200）
IE=JM（MEO,1）
JE=JM（MEO,2）
ME=JM（MEO,3）
CM=CJZ（JE,1）-CJZ（IE,1）
BM=CJZ（IE,2）-CJZ（JE,2）
CJ=CJZ（IE,1）-CJZ（ME,1）
BJ=CJZ（ME,2）-CJZ（IE,2）
AE=（BJ*CM-BM*CJ）/2.0
IF（IASK.LE.1）GOTO50
DO10I=1,3
DO10J=1,6
B（I,J）=0.0
10CONTINUE
B（1,1）=-BJ-BM
B（1,3）=BJ
B（1,5）=BM
B（2,2）=-CJ-CM
B（2,4）=CJ
B（2,6）=CM
B（3,1）=B（2,2）
B（3,2）=B（1,1）
B（3,3）=B（2,4）
B（3,4）=B（1,3）
B（3,5）=B（2,6）
!
Page151
B（3,6）=B（1,5）
DO20I=1,3
DO20J=1,6
B（I,J）=B（I,J）/（2.0*AE）
20CONTINUE
D（1,1）=EO/（1.0-UN*UN）
D（1,2）=EO*UN/（1.0-UN*UN）
D（2,1）=D（1,2）
D（2,2）=D（1,1）
D（1,3）=0.0
D（2,3）=0.0
D（3,1）=0.0
D（3,2）=0.0
D（3,3）=EO/（2.0*（1.0+UN））
DO30I=1,3
DO30J=1,6
S（I,J）=0.0
DO30K=1,3
S（I,J）=S（I,J）+D（I,K）*B（K,J）
30CONTINUE
IF（IASK.LE.2）GOTO50
DO40I=1,6
DO40J=1,6
EKE（I,J）=0.0
DO40K=1,3
!
**********************************ExchangeBAndS***********************************************
EKE（I,J）=EKE（I,J）+B（K,I）*S（K,J）*AE*TE
40CONTINUE
50CONTINUE
RETURN
END
SUBROUTINETOTSTI
COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE
COMMON/X2/JM（100,3）,NZC（50）,CJZ（100,2）,PJ（100,2）,B（3,6）,D（3,3）,S（3,6）,TKZ（200,20）,EKE（6,6）,P（200）
!
Page152
DO20I=1,NJ2
DO20J=1,NDD
TKZ（I,J）=0.0
20CONTINUE
!
*************NotUnderstanded*****************************
DO30MEO=1,NE
CALLELEST（MEO,3）
DO30I=1,3
DO30II=1,2
LH=2*（I-1）+II
LDH=2*（JM（MEO,I）-1）+II
DO30J=1,3
DO30JJ=1,2
L=2*（J-1）+JJ
LZ=2*（JM（MEO,J）-1）+JJ
LD=LZ-LDH+1
IF（LD.LE.0）GOTO30
TKZ（LDH,LD）=TKZ（LDH,LD）+EKE（LH,L）
30CONTINUE
RETURN
END
SUBROUTINELOAD
COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE
COMMON/X2/JM（100,3）,NZC（50）,CJZ（100,2）,PJ（100,2）,B（3,6）,D（3,3）,S（3,6）,TKZ（200,20）,EKE（6,6）,P（200）
DO10I=1,NJ2
P（I）=0.0
10CONTINUE
IF（NPJ.EQ.0）GOTO30
DO20I=1,NPJ
I1=I+1
J=IFIX（PJ（I1,2））
P（J）=PJ（I1,1）
20CONTINUE
30IF（GAMA.LE.0.0）GOTO50
!
Page153
DO40MEO=1,NE
CALLELEST（MEO,1）
PE=-GAMA*AE*TE/3.0
IE=JM（MEO,1）
JE=JM（MEO,2）
ME=JM（MEO,3）
P（2*IE）=P（2*IE）+PE
P（2*JE）=P（2*JE）+PE
P（2*ME）=P（2*ME）+PE
40CONTINUE
50CONTINUE
RETURN
END
SUBROUTINESUPPOR
COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE
COMMON/X2/JM（100,3）,NZC（50）,CJZ（100,2）,PJ（100,2）,B（3,6）,D（3,3）,S（3,6）,TKZ（200,20）,EKE（6,6）,P（200）
DO60I=1,NZ
MZ=NZC（I）
TKZ（MZ,1）=1.0
DO10J=2,NDD
TKZ（MZ,J）=0.0
10CONTINUE
IF（MZ-NDD）20,20,30
20JO=MZ
GOTO40
30JO=NDD
40DO50J=2,JO
J1=MZ-J
TKZ（J1+1,J）=0.0
50CONTINUE
P（MZ）=0.0
60CONTINUE
RETURN
END
!
Page154
SUBROUTINESOLVEQ
COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE
COMMON/X2/JM（100,3）,NZC（50）,CJZ（100,2）,PJ（100,2）,B（3,6）,D（3,3）,S（3,6）,TKZ（200,20）,EKE（6,6）,P（200）
NJ1=NJ2-1
DO50K=1,NJ1
IF（NJ2-K-NDD+1）10,10,20
10IM=NJ2
GOTO30
20IM=K+NDD-1
30K1=K+1
DO50I=K1,IM
L=I-K+1
C=TKZ（K,L）/TKZ（K,1）
LD1=NDD-L+1
DO40J=1,LD1
M=J+I-K
TKZ（I,J）=TKZ（I,J）-C*TKZ（K,M）
40CONTINUE
P（I）=P（I）-C*P（K）
50CONTINUE
P（NJ2）=P（NJ2）/TKZ（NJ2,1）
DO100I1=1,NJ1
I=NJ2-I1
!
************************************************************************下面一行可能出错
IF（NDD-NJ2+I-1）60,60,70
60JO=NDD
GOTO80
70JO=NJ2-I+1
80DO90J=2,JO
LH=J+I-1
P（I）=P（I）-TKZ（I,J）*P（LH）
90CONTINUE
P（I）=P（I）/TKZ（I,1）
100CONTINUE
!
Page155
WRITE（6,110）（I,P（2*I-1）,P（2*I）,I=1,NJ）
!
************************************************************************************
110FORMAT（2X,3HJD=,3X,2HU=,12X,2HV=/（I4,3X,F16.7,3X,F16.7））
RETURN
END
SUBROUTINESTRESS
COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE
COMMON/X2/JM（100,3）,NZC（50）,CJZ（100,2）,PJ（100,2）,B（3,6）,D（3,3）,S（3,6）,TKZ（200,20）,EKE（6,6）,P（200）
DIMENSIONWY（6）,YL（3）
DO60MEO=1,NE
CALLELEST（MEO,2）
DO10I=1,3
DO10J=1,2
LH=2*（I-1）+J
LDH=2*（JM（MEO,I）-1）+J
WY（LH）=P（LDH）
10
升级会员 