有限元编程的c++实现算例 1. #include<> 2. #include<> 3. 4. 5. #define ne 3 #define nj 4 #define nz 6 #define npj 0 #define npf 1 #define nj3 12 #define dd 6 #define e0 #define a0 #define i0 #define pi 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]={,,,}; 20. double gj[ne+1]={,,,}; 21. double mj[ne+1]={,a0,a0,a0}; 22. double gx[ne+1]={,i0,i0,i0}; 23. int zc[nz+1]={0,1,2,3,10,11,12}; 24. double pj[npj+1][3]={{,,}}; 25. double pf[npf+1][5]={{0,0,0,0,0},{0,-20,,,}}; 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. 36. void jdugd(int); 38. void zb(int); 39. void gdnl(int); 40. void dugd(int); 41. 42. 43. 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. 54. if(npj>0) 55. { 56. for(i=1;i<=npj;i++) 57. { 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. { hz=i; 67. gdnl(hz); 68. e=(int)pf[hz][3]; 69. zb(e); for(j=1;j<=6;j++) { 72. pe[j]=; 73. for(k=1;k<=6;k++) { 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[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. for(e=1;e<=ne;e++) { 94. dugd(e); for(i=1;i<=2;i++) { 97. for(ii=1;ii<=3;ii++) 98. { 99. h=3*(i-1)+ii; dh=3*(jm[e][i]-1)+ii; for(j=1;j<=2;j++) 102. { 103. for(jj=1;jj<=3;jj++) { 105. l=3*(j-1)+jj; zl=3*(jm[e][j]-1)+jj; dl=zl-dh+1; if(dl>0) 109. kz[dh][dl]=kz[dh][dl]+ke[h][l]; } 111. } 112. } 113. } 114. } 115. 116. for(i=1;i<=nz;i++) { 119. z=zc[i]; kz[z][l]=; for(j=2;j<=dd;j++) 122. { 123. kz[z][j]=; } 125. if((z!=1)) 126. { 127. if(z>dd) 128. j0=dd; 129. else if(z<=dd) 130. j0=z; for(j=2;j<=j0;j++) 132. kz[z-j+1][j]=; 133. } 134. p[z]=; } 136. 137. 138. 139. 140. for(k=1;k<=nj3-1;k++) 141. { 142. if(nj3>k+dd-1) 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]; 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]; } 159. } 160. 161. 162. 163. 164. p[nj3]=p[nj3]/kz[nj3][1]; for(i=nj3-1;i>=1;i--) 166. { 167. if(dd>nj3-i+1) 168. j0=nj3-i+1; 169. else j0=dd; for(j=2;j<=j0;j++) 171. { 172. h=j+i-1; 173. p[i]=p[i]-kz[i][j]*p[h]; 174. } 175. p[i]=p[i]/kz[i][1]; } 177. printf("\n"); 178. printf("_____________________________________________________________\n"); 179. printf("NJ U V CETA \n"); for(i=1;i<=nj;i++) 181. { 182. printf(" %-5d % % %\n",i,p[3*i-2],p[3*i-1],p[3*i]); 183. } 184. printf("_____________________________________________________________\n"); 185. printf("E N Q M \n"); 187. for(e=1;e<=ne;e++) { 190. jdugd(e); zb(e); 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; wy[h]=p[dh]; 199. } 200. } 201. for(i=1;i<=6;i++) 202. { 203. f[i]=; 204. for(j=1;j<=6;j++) 205. { 206. for(k=1;k<=6;k++) { 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++) if(pf[i][3]==e) { 217. hz=i; 218. gdnl(hz); for(j=1;j<=6;j++) { 221. f[j]=f[j]+f0[j]; 222. } 223. } 224. } 225. printf("%-3d(A) % % %\n",e,f[1],f[2],f[3]); printf(" (B) % % %\n",f[4],f[5],f[6]); } 228. return; 229. } 230. 232. 236. void gdnl(int hz) 237. { 238. int ind,e; 239. double g,c,l0,d; 240. 241. 242. g=pf[hz][1]; c=pf[hz][2]; e=(int)pf[hz][3]; ind=(int)pf[hz][4]; l0=gc[e]; d=l0-c; 248. if(ind==1) 249. { 250. f0[1]=; 251. f0[2]=-(g*c*(2-2*c*c/(l0*l0)+(c*c*c)/(l0*l0*l0)))/2; f0[3]=-(g*c*c)*(6-8*c/l0+3*c*c/(l0*l0))/12; 253. f0[4]=; 254. f0[5]=-g*c-f0[2];