Python和科学计算认证群组 - 讨论区
标题：Fortran的一个例题
2011年04月12日 星期二 12:38
以前没用过Fortran，今天忽然要写个东西。
这是《有限元法基础与程序设计》书中P26-P30页一个例题的程序，平时也没真正接触过Fortran，手头找的的书竟然是1990年谭浩强出版的fortran的书..，昨天刚说，今晚就用
这段代码基本没看懂，不知道什么意思，还要根据这个来编写P31页，例2-2的程序，希望有接触过的一定帮帮忙....
DIMENSION LOC(4,3),CX(6),CY(6),IFIX(6),F(12),GK(12,12),STRES(4,3),BAK(4,3,6)
COMMON NN,NE,ND,NFIX,E,ANU,T,GM,NTYPE
OPEN(5,FILE='E1.DAT',STATUS='OLD')
OPEN(6,FILE='OUT1.')
READ(5,*) NN,NE,ND,NFIX,E,ANU,T,GM,NTYPE
WRITE(6,105) NN,NE,ND,NFIX,E,ANU,T,GM,NTYPE
105 FORMAT(2X,'NN NE ND NFIX E ANU T GM NTYPE'1/4I4,E10.4,2F7.3,E10.4,I3)
READ(5,*)(LOC(I,1),LOC(I,2),LOC(I,3),I=1,NE)
READ(5,*)(CX(J),CY(J),J=1,NN)
READ(5,*)(IFIX(K),K=1,NFIX)
DO 10 I=1,ND
10 F(I)=0.0
F(2)=-1.0
CALL CST(LOC,CX,CY,IFIX,F,GK,STRES,BAK)
WRITE(6,110)
110 FORMAT(/4X,'NODE',5X,'X-DISP',8X,'Y-DISP')
WRITE(6,115)(I,F(2*I-1),F(2*I),I=1,NN)
115 FORMAT(2X,I5,2E15.6)
WRITE(6,120)
120 FORMAT(/4X,'ELEMENT',4X,'X-STR',8X,'Y-STR',8X,'XY-STR')
WRITE(6,125)(I,(STRES(I,J),J=1,3),I=1,NE)
125 FORMAT(2X,I4,3E15.6)
STOP
END
SUBROUTINE CST(LOC,CX,CY,IFIX,F,GK,STRES,BAK)
DIMENSION LOC(NE,3),CX(NN),CY(NN),IFIX(NFIX),F(ND),GK(ND,ND),STRES(NE,3),D(3,3),BB(3,6),EK(6,6),XX(6),BE(3),CE(3),BA(3,6),BAK(NE,3,6)
COMMON NN,NE,ND,NFIX,E,ANU,T,GM,NTYPE
DO 10 I=1,ND
DO 10 J=1,ND
10 GK(I,J)=0.0
DO 20 II=1,3
DO 20 JJ=1,3
20 D(II,JJ)=0.0
IF(NTYPE.EQ.1)GO TO 30
E=E/(1.0-ANU**2)
ANU=ANU/(1.0-ANU)
30 S=E/(1.0-ANU**2)
D(1,1)=S
D(1,2)=S*ANU
D(2,2)=S
D(2,1)=D(1,2)
D(3,3)=0.5*S*(1.0-ANU)
DO 100 I=1,NE
DO 40 II=1,3
DO 40 JJ=1,6
40 BB(II,JJ)=0.
I1=LOC(I,1)
I2=LOC(I,2)
I3=LOC(I,3)
BE(1)=CY(I2)-CY(I3)
BE(2)=CY(I3)-CY(I1)
BE(3)=CY(I1)-CY(I2)
CE(1)=CX(I3)-CX(I2)
CE(2)=CX(I1)-CX(I3)
CE(3)=CX(I2)-CX(I1)
S2=CX(I1)*BE(1)+CX(I2)*BE(2)+CX(I3)*BE(3)
DO 50 II=1,3
L=2*II
MM=L-1
BB(1,MM)=BE(II)/S2
BB(2,L)=CE(II)/S2
BB(3,MM)=BB(2,L)
50 BB(3,L)=BB(1,MM)
DO 60 K=1,3
DO 60 L=1,6
BA(K,L)=0.0
DO 60 MM=1,3
BA(K,L)=BA(K,L)+D(K,MM)*BB(MM,L)
60 BAK(I,K,L)=BA(K,L)
IF(GM.EQ.0.0) GO TO 65
DO 70 INODE=1,3
NODEI=LOC(I,INODE)
J2=NODEI*2
70 F(J2)=F(J2)-T*GM*(0.5*S2)/3.
65 CONTINUE
DO 75 K=1,6
DO 75 L=1,6
B1=0.0
DO 80 MM=1,3
80 B1=B1+BB(MM,K)*BA(MM,L)
EK(K,L)=0.5*S2*B1*T
WRITE(6,140) I,K,L,EK(K,L)
140 FORMAT(1X,'I K L,EK',3I4,E12.5)
75 CONTINUE
DO 85 INODE=1,3
NODEI=LOC(I,INODE)
DO 85 IDOFN=1,2
NROWS=(NODEI-1)*2+IDOFN
NROWE=(INODE-1)*2+IDOFN
DO 85 JNODE=1,3
NODEJ=LOC(I,JNODE)
DO 85 JDOFN=1,2
NCOLS=(NODEJ-1)*2+JDOFN
NCOLE=(JNODE-1)*2+JDOFN
85 GK(NROWS,NCOLS)=GK(NROWS,NCOLS)+EK(NROWE,NCOLE)
100 CONTINUE
WRITE(6,160)
160 FORMAT(/4X,'NODE',5X,'X-LOAD',8X,'Y-LOAD')
WRITE(6,165)(I,F(2*I-1),F(2*I),I=1,NN)
165 FORMAT(2X,I5,2E15.6)
WRITE(6,170) ((I,J,GK(I,J),J=1,ND),I=1,ND)
170 FORMAT(1X,'IJ,GK',2I4,E12.5,2X,2I4,E12.5)
DO 90 I=1,NFIX
IX=IFIX(I)
90 GK(IX,IX)=GK(IX,IX)*1.0E15
CALL GAUSS(GK,F,ND)
DO 95 I=1,NE
DO 95 J=1,3
XX(1)=F(2*LOC(I,1)-1)
XX(2)=F(2*LOC(I,1))
XX(3)=F(2*LOC(I,2)-1)
XX(4)=F(2*LOC(I,2))
XX(5)=F(2*LOC(I,3)-1)
XX(6)=F(2*LOC(I,3))
DO 95 K=1,6
95 STRES(I,J)=STRES(I,J)+BAK(I,J,K)*XX(K)
RETURN
END
SUBROUTINE GAUSS(A,B,N)
DIMENSION A(N,N),B(N)
DO 1 I=1,N
I1=I+1
DO 10 J=I1,N
10 A(I,J)=A(I,J)/A(I,I)
B(I)=B(I)/A(I,I)
A(I,I)=1.0
DO 20 J=I1,N
DO 30 M=I1,N
30 A(J,M)=A(J,M)-A(J,I)*A(I,M)
20 B(J)=B(J)-A(J,I)*B(I)
1 CONTINUE
DO 40 I=N-1,1,-1
DO 50 J=I+1,N
50 B(I)=B(I)-A(I,J)*B(J)
40 CONTINUE
RETURN
END
2011年04月12日 星期二 14:07
刚才有几处错误，编译过去了。
要自己写一个E1.DAT,运行后会输出OUT，这段程序相当于一个公式。
不同的题除了要改E1.DAT,还得该这个公式，代码太不友好了，
不知道这个例题怎么改
不过在试把这段代码，用F2py转一下
2011年04月12日 星期二 15:58
1 -0.235102E+03 0.196394E+01 0.179446E+02
2 -0.157505E+03 -0.457531E+00 -0.277912E+02
3 -0.157708E+03 -0.167055E+01 0.223379E+02
4 -0.791844E+02 0.230318E+01 -0.250900E+02
5 -0.798636E+02 -0.176540E+01 0.249486E+02
6 -0.132146E+01 0.505946E+01 -0.229986E+02
7 -0.204386E+01 0.733318E+00 0.267462E+02
8 0.769114E+02 0.761440E+01 -0.214570E+02
9 0.765514E+02 0.545907E+01 0.285467E+02
10 0.157258E+03 0.933106E+01 -0.201619E+02
11 0.157851E+03 0.128792E+02 0.312907E+02
12 0.244167E+03 0.868761E+01 -0.193180E+02
13 -0.236666E+03 -0.136773E+02 0.267469E+02
14 -0.154222E+03 -0.306593E+01 -0.241972E+02
15 -0.155452E+03 -0.104341E+02 0.277387E+02
16 -0.763528E+02 -0.114424E+01 -0.209798E+02
17 -0.774398E+02 -0.765246E+01 0.291528E+02
18 0.699715E+00 0.135491E+01 -0.189443E+02
19 -0.411388E+00 -0.529869E+01 0.300478E+02
20 0.781434E+02 0.375863E+01 -0.183050E+02
21 0.768310E+02 -0.410058E+01 0.297594E+02
22 0.156609E+03 0.552584E+01 -0.194234E+02
23 0.154641E+03 -0.626357E+01 0.263445E+02
24 0.233628E+03 0.666478E+01 -0.229558E+02
25 -0.230218E+03 -0.157564E+02 0.284959E+02
26 -0.149664E+03 -0.413972E+01 -0.215026E+02
27 -0.151253E+03 -0.136530E+02 0.311182E+02
28 -0.735740E+02 -0.258980E+01 -0.170711E+02
29 -0.750273E+02 -0.112916E+02 0.325446E+02
30 0.136307E+01 -0.286377E+00 -0.148682E+02
31 -0.144916E+00 -0.931559E+01 0.325238E+02
32 0.761457E+02 0.223235E+01 -0.149260E+02
33 0.744032E+02 -0.820170E+01 0.306906E+02
34 0.151287E+03 0.455364E+01 -0.173275E+02
35 0.149287E+03 -0.742116E+01 0.271067E+02
36 0.227401E+03 0.635939E+01 -0.218262E+02
37 -0.219642E+03 -0.158262E+02 0.285531E+02
38 -0.143023E+03 -0.488657E+01 -0.189837E+02
39 -0.144621E+03 -0.144541E+02 0.329691E+02
40 -0.703245E+02 -0.337098E+01 -0.132118E+02
41 -0.718512E+02 -0.125129E+02 0.350499E+02
42 0.121490E+01 -0.102348E+01 -0.104710E+02
43 -0.375519E+00 -0.105469E+02 0.348550E+02
44 0.725290E+02 0.165516E+01 -0.106971E+02
45 0.707674E+02 -0.889362E+01 0.324217E+02
46 0.144533E+03 0.420781E+01 -0.138029E+02
47 0.142525E+03 -0.781444E+01 0.279097E+02
48 0.218276E+03 0.617208E+01 -0.196537E+02
49 -0.205910E+03 -0.153885E+02 0.281363E+02
50 -0.134250E+03 -0.564081E+01 -0.162380E+02
51 -0.135694E+03 -0.142880E+02 0.341638E+02
52 -0.661368E+02 -0.399002E+01 -0.905710E+01
53 -0.675595E+02 -0.125088E+02 0.370773E+02
54 0.931111E+00 -0.148364E+01 -0.563873E+01
55 -0.583456E+00 -0.105548E+02 0.369772E+02
56 0.679125E+02 0.133486E+01 -0.589620E+01
57 0.662080E+02 -0.887227E+01 0.339661E+02
58 0.135772E+03 0.395353E+01 -0.970998E+01
59 0.133786E+03 -0.793730E+01 0.281517E+02
60 0.205529E+03 0.592302E+01 -0.170046E+02
61 -0.189029E+03 -0.147892E+02 0.274934E+02
62 -0.123305E+03 -0.648990E+01 -0.130908E+02
63 -0.124515E+03 -0.137390E+02 0.350624E+02
64 -0.608717E+02 -0.459789E+01 -0.443387E+01
65 -0.621098E+02 -0.120120E+02 0.387876E+02
66 0.641255E+00 -0.183853E+01 -0.362068E+00
67 -0.743642E+00 -0.101310E+02 0.387458E+02
68 0.622501E+02 0.114850E+01 -0.745073E+00
69 0.606225E+02 -0.859628E+01 0.350604E+02
70 0.124860E+03 0.377640E+01 -0.536556E+01
71 0.122914E+03 -0.787483E+01 0.278320E+02
72 0.189290E+03 0.563025E+01 -0.140676E+02
73 -0.168677E+03 -0.140670E+02 0.266583E+02
74 -0.110099E+03 -0.748490E+01 -0.936801E+01
75 -0.111015E+03 -0.129719E+02 0.358043E+02
76 -0.545907E+02 -0.519900E+01 0.847212E+00
77 -0.555976E+02 -0.112292E+02 0.402604E+02
78 0.177752E+00 -0.201940E+01 0.543351E+01
79 -0.105870E+01 -0.942377E+01 0.401243E+02
80 0.553727E+02 0.123115E+01 0.464918E+01
81 0.538165E+02 -0.808689E+01 0.356223E+02
82 0.111826E+03 0.385724E+01 -0.102116E+01
83 0.109908E+03 -0.762658E+01 0.269669E+02
84 0.169941E+03 0.543703E+01 -0.110806E+02
85 -0.144267E+03 -0.131910E+02 0.256406E+02
86 -0.944957E+02 -0.872626E+01 -0.476782E+01
87 -0.950347E+02 -0.119529E+02 0.364830E+02
88 -0.474860E+02 -0.586206E+01 0.710692E+01
89 -0.481935E+02 -0.100975E+02 0.415677E+02
90 -0.852685E+00 -0.209350E+01 0.119740E+02
91 -0.189477E+01 -0.833232E+01 0.411070E+02
92 0.469251E+02 0.161314E+01 0.102872E+02
93 0.454487E+02 -0.722821E+01 0.355690E+02
94 0.967355E+02 0.437142E+01 0.298641E+01
95 0.948184E+02 -0.710838E+01 0.255107E+02
96 0.148296E+03 0.557954E+01 -0.856926E+01
97 -0.114873E+03 -0.121291E+02 0.243627E+02
98 -0.763154E+02 -0.104014E+02 0.111397E+01
99 -0.763625E+02 -0.106838E+02 0.370552E+02
100 -0.398448E+02 -0.714367E+01 0.146968E+02
101 -0.400958E+02 -0.864672E+01 0.426810E+02
102 -0.307591E+01 -0.306545E+01 0.196600E+02
103 -0.370669E+01 -0.684243E+01 0.416778E+02
104 0.361817E+02 0.122671E+01 0.165065E+02
105 0.349848E+02 -0.594059E+01 0.348419E+02
106 0.794619E+02 0.477417E+01 0.646860E+01
107 0.776257E+02 -0.622222E+01 0.234206E+02
108 0.126019E+03 0.627213E+01 -0.759995E+01
109 -0.797547E+02 -0.109761E+02 0.222985E+02
110 -0.555194E+02 -0.127711E+02 0.819923E+01
111 -0.550022E+02 -0.967492E+01 0.368195E+02
112 -0.318093E+02 -0.109996E+02 0.232564E+02
113 -0.312707E+02 -0.777395E+01 0.430957E+02
114 -0.698591E+01 -0.870597E+01 0.288131E+02
115 -0.651803E+01 -0.590415E+01 0.417690E+02
116 0.220097E+02 -0.474860E+01 0.245896E+02
117 0.219974E+02 -0.482300E+01 0.337007E+02
118 0.589974E+02 0.796125E+00 0.106848E+02
119 0.580153E+02 -0.508432E+01 0.210524E+02
120 0.105840E+03 0.602229E+01 -0.938215E+01
121 -0.410355E+02 -0.103524E+02 0.176122E+02
122 -0.327636E+02 -0.154937E+02 0.136073E+02
123 -0.320189E+02 -0.110343E+02 0.330478E+02
124 -0.226070E+02 -0.212029E+02 0.292610E+02
125 -0.209083E+02 -0.110309E+02 0.405333E+02
126 -0.114652E+02 -0.275266E+02 0.378563E+02
127 -0.851256E+01 -0.984602E+01 0.404412E+02
128 0.407786E+01 -0.311056E+02 0.367410E+02
129 0.798273E+01 -0.772316E+01 0.330494E+02
130 0.323347E+02 -0.273295E+02 0.220567E+02
131 0.359554E+02 -0.564836E+01 0.206355E+02
132 0.889595E+02 -0.110236E+02 -0.994203E+01
133 -0.105990E+02 -0.117920E+02 0.668375E+01
134 -0.104442E+02 -0.144557E+02 0.706642E+01
135 -0.112543E+02 -0.193069E+02 0.197522E+02
136 -0.103426E+02 -0.323346E+02 0.215326E+02
137 -0.948441E+01 -0.271958E+02 0.269498E+02
138 -0.113690E+02 -0.646560E+02 0.347504E+02
139 -0.604854E+01 -0.327966E+02 0.297715E+02
140 -0.114945E+02 -0.107671E+03 0.463622E+02
141 0.104993E+01 -0.325541E+02 0.294520E+02
142 -0.141647E+01 -0.158214E+03 0.533252E+02
143 0.212767E+02 -0.223266E+02 0.291644E+02
144 0.601269E+02 -0.213207E+03 0.400848E+02
总算得出了结果..
这是144个单元在三个不同方向上的单元应力。
大概是就是平面的一个矩形，用一行48个，一共3行的三角形划分了144个单元，（可能是这样..）得出
了上边的结果我怎么画图，比较直观的表达出这144个单元的位置，以及应力大小呢
2011年04月12日 星期二 17:26
画出矩形，三个方向的应力nomalize以后分别代表rgb～用颜色就可以表达了。。
2011年04月12日 星期二 18:02
恩，这个知道,画过了。
主要是fortran皮毛不懂，不知道这144个怎么分布的，144/2 ，144/3，144/4，144/6，144/8，试了下
画出来的都不符合物理现象..
直接用仿真软件画了个图，应付了
2011年04月12日 星期二 20:01
看这个数组的定义GK(12,12)，应该是12x12的图像。
Zeuux © 2024
京ICP备05028076号