有限元编程的C++实现算例
read.pudn./downloads76/doc/fileformat/290377/ganjian.cpp__.htm
1. #include
2. #include
3.
4.
5. #define ne 3 //单元数
6. #define nj 4 //节点数
7. #define nz 6 //支撑数
8. #define npj 0 //节点载荷数
9. #define npf 1 //非节点载荷数
10. #define nj3 12 //节点位移总数
11. #define dd 6 //半带宽
12. #define e0 2.1E8 //弹性模量
13. #define a0 0.008 //截面积
14. #define i0 1.22E-4 //单元惯性距
15. #define pi 3.141592654
16.
17.
18. int jm[ne+1][3]={{0,0,0},{0,1,2},{0,2,3},{0,4,3}}; /*gghjghg*/
19. double gc[ne+1]={0.0,1.0,2.0,1.0};
20. double gj[ne+1]={0.0,90.0,0.0,90.0};
21. double mj[ne+1]={0.0,a0,a0,a0};
22. double gx[ne+1]={0.0,i0,i0,i0};
23. int zc[nz+1]={0,1,2,3,10,11,12};
24. double pj[npj+1][3]={{0.0,0.0,0.0}};
25. double pf[npf+1][5]={{0,0,0,0,0},{0,-20,1.0,2.0,2.0}};
26. double kz[nj3+1][dd+1],p[nj3+1];
27. double pe[7],f[7],f0[7],t[7][7];
28. double ke[7][7],kd[7][7];
29.
61. }
30.
31. //**kz[][]一整体刚度矩阵
32. //**ke[][]一整体坐标下的单元刚度矩阵 33. //**kd[][]一局部坐标下的单位刚度矩阵 34. //**皿]一坐标变换矩阵 35.
35. //**这是函数声明 36. void jdugd(int); 37. void zb(int); 38. void gdnl(int); 39. void dugd(int); 41. 42.
43. //**主程序开始 44. void main() 45. {
int i,j,k,e,dh,h,ii,jj,hz,al,bl,m,l,dl,zl,z,j0;
49.
50. //*********************************************** 51. //<功能：形成矩阵P>
52. //*********************************************** 53.
53. if(npj>0) 54. { 55. for(i=1;i<=npj;i++) 56. {
57. j=pj[i][2]; 58. p[j]=pj[i][1];
59.
}
47. double cl,wy[7]; 48. int im,in,jn; 46. 〃把节点载荷送入P
62. if(npf>0)
63.
64. for(i=1;i<=npf;i++)
65. //求固端反力F0
66. hz=i;
67. gdnl(hz);
68. e=(int)pf[hz][3];
69. zb(e); //求单元
70. for(j=1;j<=6;j++) //求坐标变换矩阵T
71.
72. pe[j]=0.0;
73. for(k=1;k<=6;k++) //求等效节点载荷
74.
75. pe[j]=pe[j]-t[k][j]*f0[k];
76.
77.
78. al=jm[e][1];
79. bl=jm[e][2];
80. p[3*al-2]=p[3*al-2]+pe[1]; //将等效节点载荷送到
81. p[3*al-1]=p[3*al-1]+pe[2];
82. p[3*al]=p[3*al]+pe[3];
83. p[3*bl-2]=p[3*bl-2]+pe[4];
84. p[3*bl-1]=p[3*bl-1]+pe[5];
85. p[3*bl]=p[3*bl]+pe[6];
86.
87.
88.
89.
90. //*************************************************
91. //<功能：生成整体刚度矩阵kz[][]>
92. for(e=1;e<=ne;e++) //按单元循环
93. (
94. dugd(e); //求整体坐标系中的单元刚度矩阵ke
95. for(i=1;i<=2;i++) //对行码循环
96. (
97. for(ii=1;ii<=3;ii++)
98. (
99. h=3*(i-1)+ii; 〃元素在ke中的行码
100. dh=3*(jm[e][i]-1)+ii; 〃该元素在KZ中的行码
101. for(j=1;j<=2;j++)
102. (
103. for(jj=1;jj<=3;jj++) //对列码循环
104. (
105. l=3*(j-1)+jj; 〃元素在ke中的列码
106. zl=3*(jm[e][j]-1)+jj; 〃该元素在KZ中的行码
107. dl=zl-dh+1; 〃该元素在KZ*中的行码
108. if(dl>0)
109. kz[dh][dl]=kz[dh][dl]+ke[h][l]; //刚度集成
110. }
111. }
112. }
113. }
114. }
115.
116. //**引入边界条件**
117. for(i=1;i<=nz;i++) //按支撑循环
118. (
119. z=zc[i]; //支撑对应的位移数
120. kz[z][l]=1.0; 〃第一列置1
121. for(j=2;j<=dd;j++)
122. (
123. kz[z][j]=0.0; //行置0
124. }
125. if((z!=1))
126. (
127. if(z>dd)
128. j0=dd;
129. else if(z<=dd)
130. j0=z; 〃列（45度斜线）置0 131. for(j=2;j<=j0;j++)
132. kz[z-j+1][j]=0.0;
133. }
134. p[z]=0.0; //P 置0
135. }
136.
137.
138.
139.
140. for(k=1;k<=nj3-1;k++)
141. (
142. if(nj3>k+dd-1) //求最大行码
143. im=k+dd-1;
144. else if(nj3<=k+dd-1)
145. im=nj3;
146. in=k+1;
147. for(i=in;i<=im;i++)
148. (
149. l=i-k+1;
150. cl=kz[k][l]/kz[k][1]; 〃修改KZ
151. jn=dd-l+1;
152. for(j=1;j<=jn;j++)
153. (
154. m=j+i-k;
155. kz[i][j]=kz[i][j]-cl*kz[k][m];
156. }
157. P[i]=P[i]-cl*p[k]; // 修改P
158. 159. 160. 161. 162. 163.
164. 165. 166. 167. 168. 169. 170. 171. 172. 173. 174. 175. 176. 177. 178. 179. 180. 181. 182. 183. 184. 185. 186. 187. 188. 189. p[nj3]=p[nj3]/kz[nj3][1];
for(i=nj3-1;i>=1;i--)
(
if(dd>nj3-i+1)
j0=nj3-i+1;
else j0=dd;
for(j=2;j<=j0;j++)
(
h=j+i-1;
p[i]=p[i]-kz[i][j]*p[h];
)
p[i]=p[i]/kz[i][1];
)
printf("\n");
printf(" ___________
//求最后一个位移分量
//求最大列码j0
//求其他位移分量
,\n")
;
printf("NJ U V CETA \n"); // 输出位移for(i=1;i<=nj;i++)
(
printf(" %-5d %14.11f %14.11f %14.11f\n",i,p[3*i-2],p[3*i-1],p[3*i]);
)
printf("\n");
//*根据E的值输出相应E单元的N,Q,M(A,B)的结果**
printf("E N Q M \n");
//*计算轴力N,剪力Q,弯矩M*
for(e=1;e<=ne;e++) // 按单元循环
190. jdugd(e); 〃求局部单元刚度矩阵kd
191. zb(e); 〃求坐标变换矩阵T
192. for(i=1;i<=2;i++)
193. {
194. for(ii=1;ii<=3;ii++)
195. {
196. h=3*(i-1)+ii;
197. dh=3*(jm[e][i]-1)+ii; //给出整体坐标下单元节点位移198. wy[h]=p[dh];
199. }
200. }
201. for(i=1;i<=6;i++)
202. {
203. f[i]=0.0;
204. for(j=1;j<=6;j++)
205. {
206. for(k=1;k<=6;k++) 〃求由节点位移引起的单元节点力207. {
208. f[i]=f[i]+kd[i][j]*t[j][k]*wy[k];
209. }
210. }
211. }
212. if(npf>0)
213. {
214. for(i=1;i<=npf;i++) //按非节点载荷数循环
215. if(pf[i][3]==e) //找到荷载所在的单元
216. {
217. hz=i;
218. gdnl(hz); //求固端反力
219. for(j=1;j<=6;j++) //将固端反力累加
220. {
221. f[j]=f[j]+f0[j];