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

有限元编程的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.

30.

31. //**kz[][]—整体刚度矩阵

32. //**ke[][]—整体坐标下的单元刚度矩阵

33. //**kd[][]—局部坐标下的单位刚度矩阵

34. //**t[][]—坐标变换矩阵

35.

36. //**这是函数声明

37. void jdugd(int);

38. void 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,wy[7];

48.int im,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;

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];//将等效节点载荷送到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++)

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))