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

有限元编程的c++实现算例单元数6.#definenj4//节点数7.#definenz6//支撑数8.#definenpj0//节点载荷数9.#definenpf1//非节点载荷数10.#definenj312//节点位移总数11.#definedd6//半带宽12.#definee02.1E8//弹性模量13.#definea00.008//截面积14.#definei01.22E-4//单元惯性距15.#definepi3.14159265416.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.//功能：形成矩阵P52.//***********************************************53.54.if(npj0)55.{56.for(i=1;i=npj;i++)57.{//把节点载荷送入P58.j=pj[i][2];59.p[j]=pj[i][1];60.}61.}62.if(npf0)63.{64.for(i=1;i=npf;i++)65.{//求固端反力F066.hz=i;67.gdnl(hz);68.e=(int)pf[hz][3];69.zb(e);//求单元号码70.for(j=1;j=6;j++)//求坐标变换矩阵T71.{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);//求整体坐标系中的单元刚度矩阵ke95.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(dl0)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;//第一列置1121.for(j=2;j=dd;j++)122.{123.kz[z][j]=0.0;//行置0124.}125.if((z!=1))126.{127.if(zdd)128.j0=dd;129.elseif(z=dd)130.j0=z;//列（45度斜线）置0131.for(j=2;j=j0;j++)132.kz[z-j+1][j]=0.0;133.}134.p[z]=0.0;//P置0135.}136.137.138.139.140.for(k=1;k=nj3-1;k++)141.{142.if(nj3k+dd-1)//求最大行码143.im=k+dd-1;144.elseif(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];//修改KZ151.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];//修改P158.}159.}160.161.162.163.164.p[nj3]=p[nj3]/kz[nj3][1];//求最后一个位移分量165.for(i=nj3-1;i=1;i--)166.{167.if(ddnj3-i+1)168.j0=nj3-i+1;169.elsej0=dd;//求最大列码j0170.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];//求其他位移分量176.}177.printf(\n);178.printf(_____________________________________________________________\n);179.printf(NJUVCETA\n);//输出位移180.for(i=1;i=nj;i++)181.{182.printf(%-5d%14.11f%14.11f%14.11f\n,i,p[3*i-2],p[3*i-1],p[3*i]);183.}184.printf(_____________________________________________________________\n);185.//*根据E的值输出相应E单元的N,Q,M(A,B)的结果**186.printf(ENQM\n);187.//*计算轴力N,剪力Q,弯矩M*188.for(e=1;e=ne;e++)//按单元循环189.{190.jdugd(e);//求局部单元刚度矩阵kd191.zb(e);//求坐标变换矩阵T192.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(npf0)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];222.}223.}224.}225.printf(%-3d(A)%9.5f%9.5f%9.5f\n,e,f[1],f[2],f[3]);//输出单元A（i）端内力226.printf((B)%9.5f%9.5f%9.5f\n,f[4],f[5],f[6]);//输出单元B（i）端内力227.}228.return;229.}230.//**主程序结束**231.232.//************************************************************233.//功能：将非节点载荷下的杆端力计算出来存入f0[]234.//************************************************************235.236.voidgdnl(inthz)237.{238.intind,e;239.doubleg,c,l0,d;240.241.242.g=pf[hz][1];//载荷值243.c=pf[hz][2];//载荷位置244.e=(int)pf[hz][3];//作用单元245.ind=(int)pf[hz][4];//载荷类型246.l0=gc[e];//杆长247.d=l0-c;248.if(ind==1)249.{250.f0[1]=0.0;251.f0[2]=-(g*c*(2-2*c*c/(l0*l0)+(c*c*c)/(l0*l0*l0)))/2;//均布载荷的固端反力252.f0[3]=-(g*c*c)*(6-8*c/l0+3*c*c/(l0*l0))/12;253.f0[4]=0.0;254.f0[5]=-g*c-f0[2];255.f0[6]=(g*c*c*c)*(4-3*c/l0)/(12*l0);256.}257.else258.{259.if(ind==2)//横向集中力的固端反力260.{261.f0[1]=0.0;262.f0[2]=(-(g*d*d)*(l0+2*c))/(l0*l0*l0);263.f0[3]=-(g*c*d*d)/(l0*l0);264.f0[4]=0.0;265.f0[5]=(-g*c*c*(l0+2*d))/(l0*l0*l0);266.f0[6]