有限元编程的c++实现算例

1、有限元编程的c+实现算例/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 /弹性模量

2、 13. #define a0 0.008 /截面积 14. #define i0 1.22E-4 /单元惯性距 15. #define pi 3.141592654 16. 17. 18. int jmne+13=0,0,0,0,1,2,0,2,3,0,4,3; /*gghjghg*/ 19. double gcne+1=0.0,1.0,2.0,1.0; 20. double gjne+1=0.0,90.0,0.0,90.0; 21. double mjne+1=0.0,a0,a0,a0; 22. double gxne+1=0.0,i0,i0,i0; 23. int zcnz+1=0,1,

3、2,3,10,11,12; 24. double pjnpj+13=0.0,0.0,0.0; 25. double pfnpf+15=0,0,0,0,0,0,-20,1.0,2.0,2.0; 26. double kznj3+1dd+1,pnj3+1; 27. double pe7,f7,f07,t77; 28. double ke77,kd77; 29. 30. 31. /*kz整体刚度矩阵 32. /*ke整体坐标下的单元刚度矩阵 33. /*kd局部坐标下的单位刚度矩阵 34. /*t坐标变换矩阵 35. 36. /*这是函数声明 37. void jdugd(int); 38. voi

4、d zb(int); 39. void gdnl(int); 40. void dugd(int); 41. 42. 43. /*主程序开始 44. void main() 45. 46. int i,j,k,e,dh,h,ii,jj,hz,al,bl,m,l,dl,zl,z,j0; 47. double cl,wy7; 48. int im,in,jn; 49. 50. /* 51. / 52. /* 53. 54. if(npj0) 55. 56. for(i=1;i0) 63. 64. for(i=1;i=npf;i+) 65. /求固端反力F0 66. hz=i; 67. gdnl(h

5、z); 68. e=(int)pfhz3; 69. zb(e); /求单元号码 70. for(j=1;j=6;j+) /求坐标变换矩阵T 71. 72. pej=0.0; 73. for(k=1;k=6;k+) /求等效节点载荷 74. 75. pej=pej-tkj*f0k; 76. 77. 78. al=jme1; 79. bl=jme2; 80. p3*al-2=p3*al-2+pe1; /将等效节点载荷送到P中 81. p3*al-1=p3*al-1+pe2; 82. p3*al=p3*al+pe3; 83. p3*bl-2=p3*bl-2+pe4; 84. p3*bl-1=p3*b

6、l-1+pe5; 85. p3*bl=p3*bl+pe6; 86. 87. 88. 89. 90. /* 91. / 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*(jmei-1)+ii; /该元素在KZ中的行码 101. for(j=1;j=2;j+) 102. 103. for(jj=1;jj0) 109. kzdhdl

7、=kzdhdl+kehl; /刚度集成 110. 111. 112. 113. 114. 115. 116. /*引入边界条件* 117. for(i=1;i=nz;i+) /按支撑循环 118. 119. z=zci; /支撑对应的位移数 120. kzzl=1.0; /第一列置1 121. for(j=2;jdd) 128. j0=dd; 129. else if(z=dd) 130. j0=z; /列（45度斜线）置0 131. for(j=2;j=j0;j+) 132. kzz-j+1j=0.0; 133. 134. pz=0.0; /P置0 135. 136. 137. 138. 1

8、39. 140. for(k=1;kk+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=kzkl/kzk1; /修改KZ 151. jn=dd-l+1; 152. for(j=1;j=1;i-) 166. 167. if(ddnj3-i+1) 168. j0=nj3-i+1; 169. else j0=dd; /求最大列码j0 170. for(j=2;j=j0;j+) 171. 17

9、2. h=j+i-1; 173. pi=pi-kzij*ph; 174. 175. pi=pi/kzi1; /求其他位移分量 176. 177. printf(n); 178. printf(_n); 179. printf(NJ U V CETA n); /输出位移 180. for(i=1;i=nj;i+) 181. 182. printf( %-5d %14.11f %14.11f %14.11fn,i,p3*i-2,p3*i-1,p3*i); 183. 184. printf(_n); 185. /*根据E的值输出相应E单元的N,Q,M(A,B)的结果* 186. printf(E N

10、 Q M n); 187. /*计算轴力N,剪力Q,弯矩M* 188. for(e=1;e=ne;e+) /按单元循环 189. 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*(jmei-1)+ii; /给出整体坐标下单元节点位移 198. wyh=pdh; 199. 200. 201. for(i=1;i=6;i+) 202. 203. fi=0.0; 204. for(j

11、=1;j=6;j+) 205. 206. for(k=1;k0) 213. 214. for(i=1;i=npf;i+) /按非节点载荷数循环 215. if(pfi3=e) /找到荷载所在的单元 216. 217. hz=i; 218. gdnl(hz); /求固端反力 219. for(j=1;j=6;j+) /将固端反力累加 220. 221. fj=fj+f0j; 222. 223. 224. 225. printf(%-3d(A) %9.5f %9.5f %9.5fn,e,f1,f2,f3); /输出单元A（i）端内力 226. printf( (B) %9.5f %9.5f %9.

12、5fn,f4,f5,f6); /输出单元B（i）端内力 227. 228. return; 229. 230. /*主程序结束* 231. 232. /* 233. / 234. /* 235. 236. void gdnl(int hz) 237. 238. int ind,e; 239. double g,c,l0,d; 240. 241. 242. g=pfhz1; /载荷值 243. c=pfhz2; /载荷位置 244. e=(int)pfhz3; /作用单元 245. ind=(int)pfhz4; /载荷类型 246. l0=gce; /杆长 247. d=l0-c; 248.

13、if(ind=1) 249. 250. f01=0.0; 251. f02=-(g*c*(2-2*c*c/(l0*l0)+(c*c*c)/(l0*l0*l0)/2; /均布载荷的固端反力 252. f03=-(g*c*c)*(6-8*c/l0+3*c*c/(l0*l0)/12; 253. f04=0.0; 254. f05=-g*c-f02; 255. f06=(g*c*c*c)*(4-3*c/l0)/(12*l0); 256. 257. else 258. 259. if(ind=2) /横向集中力的固端反力 260. 261. f01=0.0; 262. f02=(-(g*d*d)*(l0

14、+2*c)/(l0*l0*l0); 263. f03=-(g*c*d*d)/(l0*l0); 264. f04=0.0; 265. f05=(-g*c*c*(l0+2*d)/(l0*l0*l0); 266. f06=(g*c*c*d)/(l0*l0); 267. 268. else 269. 270. f01=-(g*d/l0); /纵向集中力的固端反力 271. f02=0.0; 272. f03=0.0; 273. f04=-g*c/l0; 274. f05=0.0; 275. f06=0.0; 276. 277. 278. 279. 280. /* 281. / 282. /* 283.

15、 void zb(int e) 284. 285. double ceta,co,si; 286. int i,j; 287. ceta=(gje*pi)/180; /角度变弧度 288. co=cos(ceta); 289. si=sin(ceta); 290. t11=co; /计算T右上角元素 291. t12=si; 292. t21=-si; 293. t22=co; 294. t33=1.0; 295. for(i=1;i=3;i+) 296. 297. for(j=1;j=3;j+) /计算T的左下角元素 298. 299. ti+3j+3=tij; 300. 301. 302.

16、 303. 304. 305. 306. /* 307. / 308. /* 309. void jdugd(int e) 310. 311. double A0,l0,j0; 312. int i; 313. int j; 314. 315. 316. A0=mje; /面积 317. l0=gce; /杆长 318. j0=gxe; /惯性钜 319. 320. 321. for(i=0;i=6;i+) 322. for(j=0;j=6;j+) /kd清0 323. kdij=0.0; 324. 325. kd11=e0*A0/l0; 326. kd22=12*e0*j0/pow(l0,3

17、); 327. kd32=6*e0*j0/pow(l0,2); 328. kd33=4*e0*j0/l0; 329. kd41=-kd11; 330. kd44=kd11; 331. kd52=-kd22; /计算kd左下角各元素 332. kd53=-kd32; 333. kd55=kd22; 334. kd62=kd32; 335. kd63=2*e0*j0/l0; 336. kd65=-kd32; 337. kd66=kd33; 338. 339. for(i=1;i=6;i+) 340. for(j=1;j=i;j+) /将kd左下角元素按对称原则送到右下角 341. kdji=kdij; 342. 343. 344. 345. /* 346. / 347. /* 348. void dugd(int e) 349. 350. int i,k,j,m; 351. jdugd(e