有限元编程的c++实现算例
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
有限元编程的c++实现算例
1.#include
2.#include
3.
4.
5.#definene3 //单元数
6.#definenj4 //节点数
7.#definenz6 //支撑数
8.#definenpj0 //节点载荷数
9.#definenpf1 //非节点载荷数
10.#definenj312 //节点位移总数
11.#definedd6 //半带宽
12.#definee02.1E8 //弹性模量
13.#definea00.008 //截面积
14.#definei01.22E-4 //单元惯性距
15.#definepi
16.
17.
18.intjm[ne+1][3]={{0,0,0},{0,1,2},{0,2,3},{0,4,3}}; /*gghjghg*/
19.doublegc[ne+1]={0.0,1.0,2.0,1.0};
20.doublegj[ne+1]={0.0,90.0,0.0,90.0};
21.doublemj[ne+1]={0.0,a0,a0,a0};
22.doublegx[ne+1]={0.0,i0,i0,i0};
23.intzc[nz+1]={0,1,2,3,10,11,12};
24.doublepj[npj+1][3]={{0.0,0.0,0.0}};
25.doublepf[npf+1][5]={{0,0,0,0,0},{0,-20,1.0,2.0,2.0}};
26.doublekz[nj3+1][dd+1],p[nj3+1];
27.doublepe[7],f[7],f0[7],t[7][7];
28.doubleke[7][7],kd[7][7];
29.
30.
31.//**kz[][]—整体刚度矩阵
32.//**ke[][]—整体坐标下的单元刚度矩阵
33.//**kd[][]—局部坐标下的单位刚度矩阵
34.//**t[][]—坐标变换矩阵
35.
36.//**这是函数声明
37.voidjdugd(int);
38.voidzb(int);
39.voidgdnl(int);
40.voiddugd(int);
41.
42.
43.//**主程序开始
44.voidmain()
45.{
46. inti,j,k,e,dh,h,ii,jj,hz,al,bl,m,l,dl,zl,z,j0;
47. doublecl,wy[7];
48. intim,in,jn;
49.
50.//***********************************************
51.//<功能：形成矩阵P>
52.//***********************************************
53.
54. if(npj>0)
55. {
56. for(i=1;i<=npj;i++)
57. { //把节点载荷送入P
58. j=pj[i][2];
59. p[j]=pj[i][1];
60. }
61. }
62. if(npf>0)
63. {
64. for(i=1;i<=npf;i++)
65. { //求固端反力F0
66. hz=i;
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]; //将等效节点载荷送到P中
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++)