有限元编程的c++实现算例
有限元算例二维传热c++程序源代码4
//please turn to page 158 for more detailsprintf("计算基函数系数值...\n");for(id=0;id<nx*ny*2;id++){for(i=0;i<3;i++){if(i==0) j=1,k=2;else if(i==1) j=2,k=0;else if(i==2) j=0,k=1;pE[id].A=( (pE[id].nd[j].x-pE[id].nd[i].x)*(pE[id].nd[k].y-pE[id].nd[i].y)-(pE[id].nd[j].y-pE[id].nd[i].y)*(pE[id].nd[k].x-pE[id].nd[i].x) )/2.0;D=2.0*pE[id].A;pE[id].a[i]=( pE[id].nd[j].x*pE[id].nd[k].y- pE[id].nd[k].x*pE[id].nd[j].y )/D;pE[id].b[i]=( pE[id].nd[j].y-pE[id].nd[k].y )/D;pE[id].c[i]=( pE[id].nd[k].x-pE[id].nd[j].x )/D;}}printf("OK!\n");printf("计算单元有限元特征式系数矩阵...\n");int l,m;for(i=0;i<nx;i++) //计算单元有限元特征式系数矩阵for(j=0;j<ny;j++){for(l=0;l<3;l++) //for the first triangle in the rectanglefor(m=0;m<3;m++){pE[i*2+j*ny*2].Aij[l][m]=( pE[i*2+j*ny*2].b[l]*pE[i*2+j*ny*2].b[m] +pE[i*2+j*ny*2].c[l]*pE[i*2+j*ny*2].c[m] ) * pE[i*2+j*ny*2].A;}for(l=0;l<3;l++) //for the second triangle in the rectangle for(m=0;m<3;m++){pE[i*2+j*ny*2+1].Aij[l][m]=( pE[i*2+j*ny*2+1].b[l]*pE[i*2+j*ny*2+1].b[m] + pE[i*2+j*ny*2+1].c[l]*pE[i*2+j*ny*2+1].c[m] ) * pE[i*2+j*ny*2+1].A;}}printf("OK!\n");printf("计算积分值,填充到f函数向量数组...\n");static int js[2]={4,4}; //每一层积分区间均分为4个子区间int idx=0;for(i=0;i<nx;i++) //计算积分值,填充到f函数向量数组for(j=0;j<ny;j++){for(idx=0;idx<3;idx++) //for the first triangle in the rectangle。
有限元方法编程
有限元方法编程
【最新版】
目录
1.有限元方法概述
2.有限元方法的编程步骤
3.有限元方法的应用实例
4.总结
正文
一、有限元方法概述
有限元方法是一种数值分析方法，它通过将待求解的连续体划分为有限个小的、简单的子区域（单元），然后用这些单元的近似解描述整个连续体的行为。
这种方法主要用于求解偏微分方程，特别是在固体力学、流体力学、热传导等领域有着广泛的应用。
二、有限元方法的编程步骤
1.几何建模：首先需要对问题进行几何建模，即将实际问题转化为计算机可以处理的数学模型。
这包括对物体的边界、形状等进行描述。
2.网格划分：将整个模型划分为有限个小的单元，这些单元可以是四面体、六面体等，根据问题的实际情况和求解的需要来选择。
3.选择适当的有限元公式：根据问题的性质和求解的目标，选择合适的有限元公式来描述单元内的物理量，如应力、应变等。
4.组装方程：将所有单元的公式组合起来，得到整个模型的方程。
5.求解方程：通过数值方法（如迭代法）求解得到的方程组，得到模型的解。
6.后处理：对求解结果进行分析和处理，如绘制应力分布图、应变分
布图等。
三、有限元方法的应用实例
有限元方法在许多工程领域都有广泛的应用，如飞机设计、桥梁设计、汽车设计等。
例如，在飞机设计中，可以通过有限元方法求解机翼的应力分布，从而优化机翼的设计，提高飞行性能。
四、总结
有限元方法是一种强大的数值分析工具，它可以用于求解各种复杂的工程问题。
通过几何建模、网格划分、选择适当的有限元公式、组装方程、求解方程和后处理等步骤，可以得到问题的解。
三结点三角形有限单元程序设计C 语言
fclose(fp) ;
int NN2=NN*2; int NBW=NN2;
//半带宽 NBW
float GK[50][50];
float DU[50];
int *count=new int [NN];
float (*strain)[3]=new float [NN][3];
float (*stress)[3]=new float [NN][3];
int ngn[3]; for(i=0;i<3;i++)ngn[i]=NGN[NOE-1][i];
float cn[3][2]; for(i=0;i<3;i++) for(j=0;j<2;j++) cn[i][j]=CN[ngn[i]-1][j];
float a[3]={0,0,0}; float b[3]={0,0,0}; float c[3]={0,0,0};
ｖｏｉｄ ｅｌｅｓｔｒａｉｎ＿ｓｔｒｅｓｓ（
ｉｎｔＮＯＥ，ｉｎｔＮＧＮ［］［３］，ｆｌｏａｔＣＮ［］［２］，
ｆｌｏａｔＤＵ［５０］，ｆｌｏａｔＥ０，ｆｌｏａｔＵ０，ｆｌｏａｔｔ，
ｉｎｔＮＥ，ｉｎｔＮＮ，ｆｌｏａｔｅｌｅｓｔｒａｉｎ［３］，ｆｌｏａｔｅｌｅｓｔｒｅｓｓ［３］）；
＃ｅｎｄｉｆ
Ⅴ．节点应力－应变头文件 ＮｏｄｅＳｔｒａｉｎＳｔｒｅｓｓ．ｈ
fscanf(fp,"%s",&newlin); for(i=0;i<NN;i++)
for(j=0;j<2;j++) fscanf(fp, "%f", &CN[s",&newlin); int NRE=0; fscanf(fp, "%d", &NRE) ; int *restrict=new int [NRE]; for(i=0;i<NRE;i++)
Matlab 有限元法计算分析程序编写
MATLAB的使用方法
1) 最简单的计算器使用法 求[12＋2×(7－4)]÷32的算术运算结果 （1）用键盘在MATLAB指令窗中输入一下内容 (12+2*(7-4))/3^2 （2）在上述表达式输入完成后，按【Enter】键，该指令被执行 （3）在指令执行后，MATLAB指令窗中将显示一下内容 ans = 2 [说明] 加 + 减 乘 * 除 / 或 \ (这两个符号对于数组有不同的含义) 幂 ^ “ans”是answer的缩写，其含义是运算答案，它是MATLAB的一个默 认变量
material_number = fscanf( fid_in, '%d', 1 ) ; % read material number material = zeros( material_number, 3 ) ; for i=1:1:material_number nm = fscanf( fid_in, '%d', 1 ) ; material( i, : ) = fscanf( fid_in, '%f', [1,3] ) ; % read materials definition end bc_number = fscanf( fid_in, '%d', 1 ) ; % read boundary conditions number bc = zeros( bc_number, 3 ) ; for i=1:1:bc_number % read boundary condition definition bc( i, 1 ) = fscanf( fid_in, '%d', 1 ) ; bc( i, 2 ) = fscanf( fid_in, '%d', 1 ) ; bc( i, 3 ) = fscanf( fid_in, '%f', 1 ) ; end
有限元分析中《结构力学》矩阵位移法C语言程序(附例题)
程序：#include "stdafx.h"#include "stdio.h"#include "math.h"#include "stdlib.h"void main(){int loc[3][2]={0},ifix[6]={0};float area[3]={0.0},fint[3]={0.0},cx[4]={0.0},cy[4]={0.0},f[12]={0.0},fr[12]= {0.0},fe[3][6]={0.0};int nn,ne,nd,nfix;float ea;int i,j,k;FILE *shuru,*shuchu;shuru=fopen("shuru.dat","r");shuchu=fopen("shuchu.dat","w");fscanf(shuru,"%d%d%d%d%f",&nn,&ne,&nd,&nfix,&ea);fprintf(shuchu,"nn ne nd nfix e\n%d %d %d %d %f\n",nn,ne,nd,nfix,ea);i=0;while(i<=ne-1){fscanf(shuru,"%d%d%f%f",&loc[i][0],&loc[i][1],&area[i],&fint[i]);i++;}fprintf(shuchu,"element node1 node2 area fint\n");i=0;while(i<=ne-1){fprintf(shuchu,"%d %d %d %f %f\n",i+1,loc[i][0],loc[i][1],area[i],fint[i]);i++;}j=0;while(j<=nn-1){fscanf(shuru,"%f%f",&cx[j],&cy[j]);j++;}fprintf(shuchu,"node x-coord y-coord\n");j=0;while(j<=nn-1){fprintf(shuchu,"%d %f %f\n",j+1,cx[j],cy[j]);j++;}k=0;while(k<=nfix-1){fscanf(shuru,"%d",&ifix[k]);k++;}fprintf(shuchu,"ifix=");k=0;while(k<=nfix-1){fprintf(shuchu,"%d ",ifix[k]);k++;}fprintf(shuchu,"\n");void cst(int (*loc)[2],int *ifix,float *area,float *fint,float *cx,float *cy,float*f,float *fr,float (*fe)[6],FILE *shuru,FILE *shuchu,float ea);cst(loc,ifix,area,fint,cx,cy,f,fr,fe,shuru,shuchu,ea);fprintf(shuchu,"node x-disp y-disp thita\n");i=0;while(i<=3){fprintf(shuchu,"%d %f %f %f\n",i+1,f[3*i],f[3*i+1],f[3*i+2]);i++;}fprintf(shuchu,"reaction nodal forces from the equations\n");fprintf(shuchu,"node x-load y-load moment\n");i=0;while(i<=3){fprintf(shuchu,"%d %f %f %f\n",i+1,fr[3*i],fr[3*i+1],fr[3*i+2]);i++;}fprintf(shuchu,"element axi-f shear-q moment-m\n");i=0;while(i<=ne-1){fprintf(shuchu,"%d %f %f %f %f %f %f\n",i+1,fe[i][0],fe[i][1],fe[i][2],fe[i][3],fe [i][4],fe[i][5]);i++;}fclose(shuru);fclose(shuchu);}void cst(int (*loc)[2],int *ifix,float *area,float *fint,float *cx,float *cy,float*f,float *fr,float (*fe)[6],FILE *shuru,FILE *shuchu,float ea){int np,nvd;float p1[3][6]={0.0},p2[3][6]={0.0},gk[12][12]={0.0},gk1[12][12]={0.0},al[3]= {0.0},tt[3][6][6]={0.0},bkl[3][6][6]={0.0};float t[6][6]={0.0},css[3]={0.0},snn[3]={0.0},ek[6][6]={0.0},ekl[6][6]={0.0},ekk[3] [6][6]={0.0},xx[6]={0.0},ba[6][6]={0.0};int nn=4,ne=3,nd=12,nfix=6;int ii,jj,i,j,k,l,inode,nodei,idofn,nrows,nrowe,jnode,nodej,jdofn,ncols,ncole;int i1,i2,ie,ix;float x12,y12,q,eal,eil1,eil2,eil3;i=0;{for(;i<=ne-1;i++){i1=loc[i][0];i2=loc[i][1];x12=cx[i2-1]-cx[i1-1];y12=cy[i2-1]-cy[i1-1];al[i]=sqrt(pow(x12,2)+pow(y12,2));css[i]=x12/al[i];snn[i]=y12/al[i];}}fscanf(shuru,"%d%d",&np,&nvd);if(np!=0)i=0;for(;i<=np-1;i++){fscanf(shuru,"%d%f%f%f",&i,&f[3*i],&f[3*i+1],&f[3*i+2]); }if(nvd!=0)i=0;for(;i<=nvd-1;i++){fscanf(shuru,"%d%f",&ie,&q);i1=loc[ie-1][0];i2=loc[ie-1][1];p1[ie-1][1]=q*al[ie-1]/2;p1[ie-1][2]=q*al[ie-1]*al[ie-1]/12;p1[ie-1][4]=q*al[ie-1]/2;p1[ie-1][5]=-q*al[ie-1]*al[ie-1]/12;}i=0;for(;i<=nd-1;i++){j=0;for(;j<=nd-1;j++){gk[i][j]=0.0;}}for(i=0;i<=ne-1;i++){j=0;for(;j<=5;j++){k=0;for(;k<=5;k++){ekl[j][k]=0.0;ek[j][k]=0.0;t[j][k]=0.0;}}eal=ea*area[i]/al[i];eil1=ea*fint[i]/al[i];eil2=ea*fint[i]/(al[i]*al[i]);eil3=ea*fint[i]/(al[i]*al[i]*al[i]);ekl[0][0]=eal;ekl[1][1]=12.0*eil3;ekl[2][2]=4.0*eil1;ekl[3][3]=eal;ekl[4][4]=12.0*eil3;ekl[5][5]=4.0*eil1;ekl[2][1]=6.0*eil2;ekl[3][0]=-eal;ekl[4][1]=-12.0*eil3;ekl[4][2]=-6.0*eil2;ekl[5][1]=6.0*eil2;ekl[5][2]=2.0*eil1;ekl[5][4]=-6.0*eil2;ii=0;for(;ii<=4;ii++){jj=ii+1;for(;jj<=5;jj++){ekl[ii][jj]=ekl[jj][ii];}}k=0;for(;k<=5;k++){l=0;for(;l<=5;l++){{ekk[i][k][l]=ekl[k][l];fprintf(shuchu,"%d %d %d %f %f\n",i+1,k+1,l+1,ekl[k][l],ekk[i][k][l]);} }}t[0][0]=css[i];t[0][1]=-snn[i];t[1][0]=snn[i];t[1][1]=css[i];t[2][2]=1.0;t[3][3]=css[i];t[3][4]=-snn[i];t[4][3]=snn[i];t[4][4]=css[i];t[5][5]=1.0;j=0;for(;j<=5;j++){k=0;for(;k<=5;k++){tt[i][j][k]=t[j][k];p2[i][j]=p2[i][j]+t[j][k]*p1[i][k];}}ii=0;for(;ii<=5;ii++){j=0;for(;j<=5;j++){ba[ii][j]=0.0;k=0;for(;k<=5;k++){ba[ii][j]=ba[ii][j]+tt[i][ii][k]*ekl[k][j];}}}ek[ii][j]=0.0;ii=0;for(;ii<=5;ii++){j=0;for(;j<=5;j++){k=0;for(;k<=5;k++){ek[ii][j]=ek[ii][j]+ba[ii][k]*tt[i][j][k];}}}j=0;for(;j<=5;j++){ii=0;for(;ii<=5;ii++){fprintf(shuchu,"i,ii,j,ek,tt=%d %d %d %f %f\n",i+1,ii+1,j+1,ek[ii][j],tt[i][ii][j]); }}inode=0;while(inode<=1){nodei=loc[i][inode];idofn=0;while(idofn<=2){nrows=(nodei-1)*3+idofn;nrowe=inode*3+idofn;f[nrows]=f[nrows]+p2[i][nrowe];jnode=0;while(jnode<=1){nodej=loc[i][jnode];jdofn=0;while(jdofn<=2){ncols=(nodej-1)*3+jdofn;ncole=jnode*3+jdofn;gk[nrows][ncols]=gk[nrows][ncols]+ek[nrowe][ncole];jdofn++;}jnode++;}idofn++;}inode++;}}i=0;for(;i<=nd-1;i++){j=0;for(;j<=nd-1;j++){gk1[i][j]=gk[i][j];}}fprintf(shuchu,"nodal forces from applied loads\node x-load y-load moment\n"); i=0;for(;i<=nn-1;i++){fprintf(shuchu,"%d %f %f %f\n",i+1,f[3*i],f[3*i+1],f[3*i+2]);}i=0;for(;i<=nd-1;i++){j=0;for(;j<=nd-1;j++){fprintf(shuchu,"i,j,gk1 %d %d %f\n",i+1,j+1,gk[i][j]);}}i=0;for(;i<=nfix-1;i++){ix=ifix[i];gk[ix-1][ix-1]=gk[ix-1][ix-1]*1.0e20;}void gauss(float (*a)[12],float *b,int n);gauss(gk,f,nd);i=0;for(;i<=nd-1;i++){fr[i]=0.0;j=0;for(;j<=nd-1;j++){fr[i]=fr[i]+gk1[i][j]*f[j];fr[i]=fr[i]-f[i];}}for(i=0;i<=ne-1;i++){i1=loc[i][0];i2=loc[i][1];xx[0]=f[3*i1-3];xx[1]=f[3*i1-2];xx[2]=f[3*i1-1];xx[3]=f[3*i2-3];xx[4]=f[3*i2-2];xx[5]=f[3*i2-1];j=0;for(;j<=5;j++){k=0;for(;k<=5;k++){bkl[i][j][k]=0.0;l=0;for(;l<=5;l++){bkl[i][j][k]=bkl[i][j][k]+ekk[i][j][l]*tt[i][k][l];}}}j=0;for(;j<=5;j++){fe[i][j]=0.0;k=0;for(;k<=5;k++){fe[i][j]=fe[i][j]+bkl[i][j][k]*xx[k];}}j=0;for(;j<=5;j++){fe[i][j]=fe[i][j]-p1[i][j];}}}void gauss(float (*a)[12],float *b,int n){int i,i1,j,m;i=0;for(;i<=n-1;i++){i1=i+1;for(j=i1;j<=n-1;j++)a[i][j]=a[i][j]/a[i][i];b[i]=b[i]/a[i][i];a[i][i]=1.0;for(j=i1;j<=n-1;j++){for(m=i1;m<=n-1;m++)a[j][m]=a[j][m]-a[j][i]*a[i][m];b[j]=b[j]-a[j][i]*b[i];}}i=n-2;for(;i>=0;i--){j=i+1;for(;j<=n-1;j++){b[i]=b[i]-a[i][j]*b[j];}}}一、如图所示平面刚架的内力，各杆面积A=76.3cm2,惯性矩I=15760cm4，弹性模量E=2×105MPa程序：#include "stdafx.h"#include "stdio.h"#include "math.h"#include "stdlib.h"void main(){int loc[3][2]={0},ifix[6]={0};float area[3]={0.0},fint[3]={0.0},cx[4]={0.0},cy[4]={0.0},f[12]={0.0},fr[12]={0.0},fe[3][6]={0.0};int nn,ne,nd,nfix;float ea;int i,j,k;FILE *shuru,*shuchu;shuru=fopen("shuru.dat","r");shuchu=fopen("shuchu.dat","w");fscanf(shuru,"%d%d%d%d%f",&nn,&ne,&nd,&nfix,&ea);fprintf(shuchu,"nn ne nd nfix e\n%d %d %d %d %f\n",nn,ne,nd,nfix,ea);i=0;while(i<=ne-1){fscanf(shuru,"%d%d%f%f",&loc[i][0],&loc[i][1],&area[i],&fint[i]);fprintf(shuchu,"element node1 node2 area fint\n");i=0;while(i<=ne-1){fprintf(shuchu,"%d %d %d %f %f\n",i+1,loc[i][0],loc[i][1],area[i],fint[i]);i++;}j=0;while(j<=nn-1){fscanf(shuru,"%f%f",&cx[j],&cy[j]);j++;}fprintf(shuchu,"node x-coord y-coord\n");j=0;while(j<=nn-1){fprintf(shuchu,"%d %f %f\n",j+1,cx[j],cy[j]);j++;}k=0;while(k<=nfix-1){fscanf(shuru,"%d",&ifix[k]);k++;}fprintf(shuchu,"ifix=");k=0;while(k<=nfix-1){fprintf(shuchu,"%d ",ifix[k]);k++;}fprintf(shuchu,"\n");void cst(int (*loc)[2],int *ifix,float *area,float *fint,float *cx,float *cy,float*f,float *fr,float (*fe)[6],FILE *shuru,FILE *shuchu,float ea);cst(loc,ifix,area,fint,cx,cy,f,fr,fe,shuru,shuchu,ea);fprintf(shuchu,"node x-disp y-disp thita\n");i=0;while(i<=3){fprintf(shuchu,"%d %f %f %f\n",i+1,f[3*i],f[3*i+1],f[3*i+2]);i++;}fprintf(shuchu,"reaction nodal forces from the equations\n");fprintf(shuchu,"node x-load y-load moment\n");i=0;while(i<=3){fprintf(shuchu,"%d %f %f %f\n",i+1,fr[3*i],fr[3*i+1],fr[3*i+2]);i++;}fprintf(shuchu,"element axi-f shear-q moment-m\n");i=0;while(i<=ne-1){fprintf(shuchu,"%d %f %f %f %f %f %f\n",i+1,fe[i][0],fe[i][1],fe[i][2],fe[i][3],fe [i][4],fe[i][5]);fclose(shuru);fclose(shuchu);}void cst(int (*loc)[2],int *ifix,float *area,float *fint,float *cx,float *cy,float*f,float *fr,float (*fe)[6],FILE *shuru,FILE *shuchu,float ea){int np,nvd;float p1[3][6]={0.0},p2[3][6]={0.0},gk[12][12]={0.0},gk1[12][12]={0.0},al[3]= {0.0},tt[3][6][6]={0.0},bkl[3][6][6]={0.0};float t[6][6]={0.0},css[3]={0.0},snn[3]={0.0},ek[6][6]={0.0},ekl[6][6]={0.0},ekk[3] [6][6]={0.0},xx[6]={0.0},ba[6][6]={0.0};int nn=4,ne=3,nd=12,nfix=6;int ii,jj,i,j,k,l,inode,nodei,idofn,nrows,nrowe,jnode,nodej,jdofn,ncols,ncole;int i1,i2,ie,ix;float x12,y12,q,eal,eil1,eil2,eil3;i=0;{for(;i<=ne-1;i++){i1=loc[i][0];i2=loc[i][1];x12=cx[i2-1]-cx[i1-1];y12=cy[i2-1]-cy[i1-1];al[i]=sqrt(pow(x12,2)+pow(y12,2));css[i]=x12/al[i];snn[i]=y12/al[i];}}fscanf(shuru,"%d%d",&np,&nvd);if(np!=0)i=0;for(;i<=np-1;i++){fscanf(shuru,"%d%f%f%f",&i,&f[3*i],&f[3*i+1],&f[3*i+2]);}if(nvd!=0)i=0;for(;i<=nvd-1;i++){fscanf(shuru,"%d%f",&ie,&q);i1=loc[ie-1][0];i2=loc[ie-1][1];p1[ie-1][1]=q*al[ie-1]/2;p1[ie-1][2]=q*al[ie-1]*al[ie-1]/12;p1[ie-1][4]=q*al[ie-1]/2;p1[ie-1][5]=-q*al[ie-1]*al[ie-1]/12;}i=0;for(;i<=nd-1;i++){j=0;for(;j<=nd-1;j++){gk[i][j]=0.0;}}for(i=0;i<=ne-1;i++){j=0;for(;j<=5;j++){k=0;for(;k<=5;k++){ekl[j][k]=0.0;ek[j][k]=0.0;t[j][k]=0.0;}}eal=ea*area[i]/al[i];eil1=ea*fint[i]/al[i];eil2=ea*fint[i]/(al[i]*al[i]);eil3=ea*fint[i]/(al[i]*al[i]*al[i]); ekl[0][0]=eal;ekl[1][1]=12.0*eil3;ekl[2][2]=4.0*eil1;ekl[3][3]=eal;ekl[4][4]=12.0*eil3;ekl[5][5]=4.0*eil1;ekl[2][1]=6.0*eil2;ekl[3][0]=-eal;ekl[4][1]=-12.0*eil3;ekl[4][2]=-6.0*eil2;ekl[5][1]=6.0*eil2;ekl[5][2]=2.0*eil1;ekl[5][4]=-6.0*eil2;ii=0;for(;ii<=4;ii++){jj=ii+1;for(;jj<=5;jj++){ekl[ii][jj]=ekl[jj][ii];}}k=0;for(;k<=5;k++){l=0;for(;l<=5;l++){{ekk[i][k][l]=ekl[k][l];fprintf(shuchu,"%d %d %d %f %f\n",i+1,k+1,l+1,ekl[k][l],ekk[i][k][l]);} }}t[0][0]=css[i];t[0][1]=-snn[i];t[1][0]=snn[i];t[1][1]=css[i];t[2][2]=1.0;t[3][3]=css[i];t[3][4]=-snn[i];t[4][3]=snn[i];t[4][4]=css[i];t[5][5]=1.0;j=0;for(;j<=5;j++){k=0;for(;k<=5;k++){tt[i][j][k]=t[j][k];p2[i][j]=p2[i][j]+t[j][k]*p1[i][k];}}ii=0;for(;ii<=5;ii++){j=0;for(;j<=5;j++){ba[ii][j]=0.0;k=0;for(;k<=5;k++){ba[ii][j]=ba[ii][j]+tt[i][ii][k]*ekl[k][j];}}}ek[ii][j]=0.0;ii=0;for(;ii<=5;ii++){j=0;for(;j<=5;j++){k=0;for(;k<=5;k++){ek[ii][j]=ek[ii][j]+ba[ii][k]*tt[i][j][k];}}}j=0;for(;j<=5;j++){ii=0;for(;ii<=5;ii++){fprintf(shuchu,"i,ii,j,ek,tt=%d %d %d %f %f\n",i+1,ii+1,j+1,ek[ii][j],tt[i][ii][j]); }}inode=0;while(inode<=1){nodei=loc[i][inode];idofn=0;while(idofn<=2){nrows=(nodei-1)*3+idofn;nrowe=inode*3+idofn;f[nrows]=f[nrows]+p2[i][nrowe];jnode=0;while(jnode<=1){nodej=loc[i][jnode];jdofn=0;while(jdofn<=2){ncols=(nodej-1)*3+jdofn;ncole=jnode*3+jdofn;gk[nrows][ncols]=gk[nrows][ncols]+ek[nrowe][ncole];jdofn++;}jnode++;}idofn++;}inode++;}}i=0;for(;i<=nd-1;i++){j=0;for(;j<=nd-1;j++){gk1[i][j]=gk[i][j];}}fprintf(shuchu,"nodal forces from applied loads\node x-load y-load moment\n"); i=0;for(;i<=nn-1;i++){fprintf(shuchu,"%d %f %f %f\n",i+1,f[3*i],f[3*i+1],f[3*i+2]);}i=0;for(;i<=nd-1;i++){j=0;for(;j<=nd-1;j++){fprintf(shuchu,"i,j,gk1 %d %d %f\n",i+1,j+1,gk[i][j]); }}i=0;for(;i<=nfix-1;i++){ix=ifix[i];gk[ix-1][ix-1]=gk[ix-1][ix-1]*1.0e20;}void gauss(float (*a)[12],float *b,int n);gauss(gk,f,nd);i=0;for(;i<=nd-1;i++){fr[i]=0.0;j=0;for(;j<=nd-1;j++){fr[i]=fr[i]+gk1[i][j]*f[j];fr[i]=fr[i]-f[i];}}for(i=0;i<=ne-1;i++){i1=loc[i][0];i2=loc[i][1];xx[0]=f[3*i1-3];xx[1]=f[3*i1-2];xx[2]=f[3*i1-1];xx[3]=f[3*i2-3];xx[4]=f[3*i2-2];xx[5]=f[3*i2-1];j=0;for(;j<=5;j++){k=0;for(;k<=5;k++){bkl[i][j][k]=0.0;l=0;for(;l<=5;l++){bkl[i][j][k]=bkl[i][j][k]+ekk[i][j][l]*tt[i][k][l];}}}j=0;for(;j<=5;j++){fe[i][j]=0.0;k=0;for(;k<=5;k++){fe[i][j]=fe[i][j]+bkl[i][j][k]*xx[k]; }}j=0;for(;j<=5;j++){fe[i][j]=fe[i][j]-p1[i][j];}}}void gauss(float (*a)[12],float *b,int n) {int i,i1,j,m;i=0;for(;i<=n-1;i++){i1=i+1;for(j=i1;j<=n-1;j++)a[i][j]=a[i][j]/a[i][i];b[i]=b[i]/a[i][i];a[i][i]=1.0;for(j=i1;j<=n-1;j++){for(m=i1;m<=n-1;m++)a[j][m]=a[j][m]-a[j][i]*a[i][m];b[j]=b[j]-a[j][i]*b[i];}}i=n-2;for(;i>=0;i--){j=i+1;for(;j<=n-1;j++){b[i]=b[i]-a[i][j]*b[j];}}}程序输出：nn ne nd nfix e4 3 12 6 200000000.000000element node1 node2 area fint1 12 0.007630 0.0001582 3 1 0.007630 0.0001583 24 0.007630 0.000158node x-coord y-coord1 0.000000 5.0000002 6.400000 5.0000003 0.000000 0.0000004 9.600000 0.000000ifix=7 8 9 10 11 121 1 1 238437.500000 238437.5000001 12 0.000000 0.0000001 1 3 0.000000 0.0000001 1 4 -238437.500000 -238437.5000001 1 5 0.000000 0.0000001 1 6 0.000000 0.0000001 2 1 0.000000 0.0000001 2 2 1442.871094 1442.8710941 2 3 4617.187500 4617.1875001 2 4 0.000000 0.0000001 2 5 -1442.871094 -1442.8710941 2 6 4617.187500 4617.1875001 3 1 0.000000 0.0000001 32 4617.187500 4617.1875001 3 3 19700.000000 19700.0000001 3 4 0.000000 0.0000001 3 5 -4617.187500 -4617.1875001 3 6 9850.000000 9850.0000001 4 1 -238437.500000 -238437.5000001 42 0.000000 0.0000001 4 3 0.000000 0.0000001 4 4 238437.500000 238437.5000001 4 5 0.000000 0.0000001 4 6 0.000000 0.0000001 5 1 0.000000 0.0000001 52 -1442.871094 -1442.8710941 5 3 -4617.187500 -4617.1875001 5 4 0.000000 0.0000001 5 5 1442.871094 1442.8710941 5 6 -4617.187500 -4617.1875001 6 1 0.000000 0.0000001 62 4617.187500 4617.1875001 6 3 9850.000000 9850.0000001 6 4 0.000000 0.0000001 6 5 -4617.187500 -4617.1875001 6 6 19700.000000 19700.000000i,ii,j,ek,tt=1 1 1 238437.500000 1.000000 i,ii,j,ek,tt=1 2 1 0.000000 0.000000i,ii,j,ek,tt=1 3 1 0.000000 0.000000i,ii,j,ek,tt=1 4 1 -238437.500000 0.000000 i,ii,j,ek,tt=1 5 1 0.000000 0.000000i,ii,j,ek,tt=1 6 1 0.000000 0.000000i,ii,j,ek,tt=1 1 2 0.000000 0.000000i,ii,j,ek,tt=1 2 2 1442.871094 1.000000i,ii,j,ek,tt=1 3 2 4617.187500 0.000000 i,ii,j,ek,tt=1 4 2 0.000000 0.000000i,ii,j,ek,tt=1 5 2 -1442.871094 0.000000 i,ii,j,ek,tt=1 6 2 4617.187500 0.000000 i,ii,j,ek,tt=1 1 3 0.000000 0.000000i,ii,j,ek,tt=1 2 3 4617.187500 0.000000 i,ii,j,ek,tt=1 3 3 19700.000000 1.000000 i,ii,j,ek,tt=1 4 3 0.000000 0.000000i,ii,j,ek,tt=1 5 3 -4617.187500 0.000000 i,ii,j,ek,tt=1 6 3 9850.000000 0.000000 i,ii,j,ek,tt=1 1 4 -238437.500000 0.000000 i,ii,j,ek,tt=1 2 4 0.000000 0.000000i,ii,j,ek,tt=1 3 4 0.000000 0.000000i,ii,j,ek,tt=1 4 4 238437.500000 1.000000 i,ii,j,ek,tt=1 5 4 0.000000 0.000000i,ii,j,ek,tt=1 6 4 0.000000 0.000000i,ii,j,ek,tt=1 1 5 0.000000 0.000000i,ii,j,ek,tt=1 2 5 -1442.871094 0.000000 i,ii,j,ek,tt=1 3 5 -4617.187500 0.000000 i,ii,j,ek,tt=1 4 5 0.000000 0.000000i,ii,j,ek,tt=1 5 5 1442.871094 1.000000 i,ii,j,ek,tt=1 6 5 -4617.187500 0.000000 i,ii,j,ek,tt=1 1 6 0.000000 0.000000i,ii,j,ek,tt=1 2 6 4617.187500 0.000000 i,ii,j,ek,tt=1 3 6 9850.000000 0.000000 i,ii,j,ek,tt=1 4 6 0.000000 0.000000i,ii,j,ek,tt=1 5 6 -4617.187500 0.000000 i,ii,j,ek,tt=1 6 6 19700.000000 1.000000 2 1 1 305200.000000 305200.0000002 1 2 0.000000 0.0000002 13 0.000000 0.0000002 1 4 -305200.000000 -305200.0000002 1 5 0.000000 0.0000002 1 6 0.000000 0.0000002 2 1 0.000000 0.0000002 2 2 3025.919922 3025.9199222 23 7564.800293 7564.8002932 2 4 0.000000 0.0000002 2 5 -3025.919922 -3025.9199222 2 6 7564.800293 7564.8002932 3 1 0.000000 0.0000002 3 2 7564.800293 7564.8002932 3 3 25216.001953 25216.0019532 3 4 0.000000 0.0000002 3 5 -7564.800293 -7564.8002932 3 6 12608.000977 12608.0009772 4 1 -305200.000000 -305200.0000002 4 2 0.000000 0.0000002 43 0.000000 0.0000002 4 4 305200.000000 305200.0000002 4 5 0.000000 0.0000002 4 6 0.000000 0.0000002 5 1 0.000000 0.0000002 5 2 -3025.919922 -3025.9199222 53 -7564.800293 -7564.8002932 5 4 0.000000 0.0000002 5 5 3025.919922 3025.9199222 5 6 -7564.800293 -7564.8002932 6 1 0.000000 0.0000002 6 2 7564.800293 7564.8002932 63 12608.000977 12608.0009772 6 4 0.000000 0.0000002 6 5 -7564.800293 -7564.8002932 6 6 25216.001953 25216.001953i,ii,j,ek,tt=2 1 1 3025.919922 0.000000 i,ii,j,ek,tt=2 2 1 0.000000 1.000000i,ii,j,ek,tt=2 3 1 -7564.800293 0.000000 i,ii,j,ek,tt=2 4 1 -3025.919922 0.000000 i,ii,j,ek,tt=2 5 1 0.000000 0.000000i,ii,j,ek,tt=2 6 1 -7564.800293 0.000000 i,ii,j,ek,tt=2 1 2 0.000000 -1.000000i,ii,j,ek,tt=2 2 2 305200.000000 0.000000 i,ii,j,ek,tt=2 3 2 0.000000 0.000000i,ii,j,ek,tt=2 4 2 0.000000 0.000000i,ii,j,ek,tt=2 5 2 -305200.000000 0.000000 i,ii,j,ek,tt=2 6 2 0.000000 0.000000i,ii,j,ek,tt=2 1 3 -7564.800293 0.000000 i,ii,j,ek,tt=2 2 3 0.000000 0.000000i,ii,j,ek,tt=2 3 3 25216.001953 1.000000 i,ii,j,ek,tt=2 4 3 7564.800293 0.000000 i,ii,j,ek,tt=2 5 3 0.000000 0.000000i,ii,j,ek,tt=2 6 3 12608.000977 0.000000 i,ii,j,ek,tt=2 1 4 -3025.919922 0.000000 i,ii,j,ek,tt=2 2 4 0.000000 0.000000i,ii,j,ek,tt=2 3 4 7564.800293 0.000000 i,ii,j,ek,tt=2 4 4 3025.919922 0.000000 i,ii,j,ek,tt=2 5 4 0.000000 1.000000i,ii,j,ek,tt=2 6 4 7564.800293 0.000000i,ii,j,ek,tt=2 1 5 0.000000 0.000000i,ii,j,ek,tt=2 2 5 -305200.000000 0.000000 i,ii,j,ek,tt=2 3 5 0.000000 0.000000i,ii,j,ek,tt=2 4 5 0.000000 -1.000000i,ii,j,ek,tt=2 5 5 305200.000000 0.000000 i,ii,j,ek,tt=2 6 5 0.000000 0.000000i,ii,j,ek,tt=2 1 6 -7564.800293 0.000000 i,ii,j,ek,tt=2 2 6 0.000000 0.000000i,ii,j,ek,tt=2 3 6 12608.000977 0.000000 i,ii,j,ek,tt=2 4 6 7564.800293 0.000000 i,ii,j,ek,tt=2 5 6 0.000000 0.000000i,ii,j,ek,tt=2 6 6 25216.001953 1.000000 3 1 1 257061.218750 257061.2187503 1 2 0.000000 0.0000003 1 3 0.000000 0.0000003 14 -257061.218750 -257061.2187503 1 5 0.000000 0.0000003 1 6 0.000000 0.0000003 2 1 0.000000 0.0000003 2 2 1808.063232 1808.0632323 2 3 5366.628906 5366.6289063 24 0.000000 0.0000003 2 5 -1808.063232 -1808.0632323 2 6 5366.628906 5366.6289063 3 1 0.000000 0.0000003 3 2 5366.628906 5366.6289063 3 3 21238.716797 21238.7167973 34 0.000000 0.0000003 3 5 -5366.628906 -5366.6289063 3 6 10619.358398 10619.3583983 4 1 -257061.218750 -257061.2187503 4 2 0.000000 0.0000003 4 3 0.000000 0.0000003 4 4 257061.218750 257061.2187503 4 5 0.000000 0.0000003 4 6 0.000000 0.0000003 5 1 0.000000 0.0000003 5 2 -1808.063232 -1808.0632323 5 3 -5366.628906 -5366.6289063 54 0.000000 0.0000003 5 5 1808.063232 1808.0632323 5 6 -5366.628906 -5366.6289063 6 1 0.000000 0.0000003 6 2 5366.628906 5366.6289063 6 3 10619.358398 10619.3583983 64 0.000000 0.0000003 6 5 -5366.628906 -5366.6289063 6 6 21238.716797 21238.716797i,ii,j,ek,tt=3 1 1 75979.257813 0.539054 i,ii,j,ek,tt=3 2 1 -115892.476563 -0.842271 i,ii,j,ek,tt=3 3 1 4520.158203 0.000000i,ii,j,ek,tt=3 4 1 -75979.257813 0.000000 i,ii,j,ek,tt=3 5 1 115892.476563 0.000000 i,ii,j,ek,tt=3 6 1 4520.158203 0.000000i,ii,j,ek,tt=3 1 2 -115892.476563 0.842271 i,ii,j,ek,tt=3 2 2 182890.046875 0.539054 i,ii,j,ek,tt=3 3 2 2892.901367 0.000000i,ii,j,ek,tt=3 4 2 115892.476563 0.000000 i,ii,j,ek,tt=3 5 2 -182890.046875 0.000000 i,ii,j,ek,tt=3 6 2 2892.901367 0.000000i,ii,j,ek,tt=3 1 3 4520.158203 0.000000i,ii,j,ek,tt=3 2 3 2892.901367 0.000000i,ii,j,ek,tt=3 3 3 21238.716797 1.000000 i,ii,j,ek,tt=3 4 3 -4520.158203 0.000000 i,ii,j,ek,tt=3 5 3 -2892.901367 0.000000 i,ii,j,ek,tt=3 6 3 10619.358398 0.000000 i,ii,j,ek,tt=3 1 4 -75979.257813 0.000000 i,ii,j,ek,tt=3 2 4 115892.476563 0.000000 i,ii,j,ek,tt=3 3 4 -4520.158203 0.000000 i,ii,j,ek,tt=3 4 4 75979.257813 0.539054 i,ii,j,ek,tt=3 5 4 -115892.476563 -0.842271 i,ii,j,ek,tt=3 6 4 -4520.158203 0.000000 i,ii,j,ek,tt=3 1 5 115892.476563 0.000000 i,ii,j,ek,tt=3 2 5 -182890.046875 0.000000 i,ii,j,ek,tt=3 3 5 -2892.901367 0.000000 i,ii,j,ek,tt=3 4 5 -115892.476563 0.842271 i,ii,j,ek,tt=3 5 5 182890.046875 0.539054 i,ii,j,ek,tt=3 6 5 -2892.901367 0.000000 i,ii,j,ek,tt=3 1 6 4520.158203 0.000000i,ii,j,ek,tt=3 2 6 2892.901367 0.000000i,ii,j,ek,tt=3 3 6 10619.358398 0.000000 i,ii,j,ek,tt=3 4 6 -4520.158203 0.000000 i,ii,j,ek,tt=3 5 6 -2892.901367 0.000000 i,ii,j,ek,tt=3 6 6 21238.716797 1.000000 nodal forces from applied loadsode x-load y-load moment1 0.000000 -192.000000 -204.8000032 0.000000 -192.000000 204.8000033 0.000000 0.000000 0.0000004 0.000000 0.000000 0.000000 i,j,gk1 1 1 241463.421875i,j,gk1 1 2 0.000000i,j,gk1 1 3 7564.800293i,j,gk1 1 4 -238437.500000i,j,gk1 1 5 0.000000i,j,gk1 1 6 0.000000i,j,gk1 1 7 -3025.919922i,j,gk1 1 8 0.000000i,j,gk1 1 9 7564.800293i,j,gk1 1 10 0.000000i,j,gk1 1 11 0.000000i,j,gk1 1 12 0.000000i,j,gk1 2 1 0.000000i,j,gk1 2 2 306642.875000i,j,gk1 2 3 4617.187500i,j,gk1 2 4 0.000000i,j,gk1 2 5 -1442.871094i,j,gk1 2 6 4617.187500i,j,gk1 2 7 0.000000i,j,gk1 2 8 -305200.000000i,j,gk1 2 9 0.000000i,j,gk1 2 10 0.000000i,j,gk1 2 11 0.000000i,j,gk1 2 12 0.000000i,j,gk1 3 1 7564.800293i,j,gk1 3 2 4617.187500i,j,gk1 3 3 44916.000000i,j,gk1 3 4 0.000000i,j,gk1 3 5 -4617.187500i,j,gk1 3 6 9850.000000i,j,gk1 3 7 -7564.800293i,j,gk1 3 8 0.000000i,j,gk1 3 9 12608.000977i,j,gk1 3 10 0.000000i,j,gk1 3 11 0.000000i,j,gk1 3 12 0.000000i,j,gk1 4 1 -238437.500000i,j,gk1 4 2 0.000000i,j,gk1 4 3 0.000000i,j,gk1 4 4 314416.750000i,j,gk1 4 5 -115892.476563i,j,gk1 4 6 4520.158203i,j,gk1 4 8 0.000000i,j,gk1 4 9 0.000000i,j,gk1 4 10 -75979.257813 i,j,gk1 4 11 115892.476563 i,j,gk1 4 12 4520.158203 i,j,gk1 5 1 0.000000i,j,gk1 5 2 -1442.871094 i,j,gk1 5 3 -4617.187500 i,j,gk1 5 4 -115892.476563 i,j,gk1 5 5 184332.921875 i,j,gk1 5 6 -1724.286133 i,j,gk1 5 7 0.000000i,j,gk1 5 8 0.000000i,j,gk1 5 9 0.000000i,j,gk1 5 10 115892.476563 i,j,gk1 5 11 -182890.046875 i,j,gk1 5 12 2892.901367 i,j,gk1 6 1 0.000000i,j,gk1 6 2 4617.187500i,j,gk1 6 3 9850.000000i,j,gk1 6 4 4520.158203i,j,gk1 6 5 -1724.286133 i,j,gk1 6 6 40938.718750 i,j,gk1 6 7 0.000000i,j,gk1 6 8 0.000000i,j,gk1 6 9 0.000000i,j,gk1 6 10 -4520.158203 i,j,gk1 6 11 -2892.901367 i,j,gk1 6 12 10619.358398 i,j,gk1 7 1 -3025.919922 i,j,gk1 7 2 0.000000i,j,gk1 7 3 -7564.800293 i,j,gk1 7 4 0.000000i,j,gk1 7 5 0.000000i,j,gk1 7 6 0.000000i,j,gk1 7 7 3025.919922i,j,gk1 7 8 0.000000i,j,gk1 7 9 -7564.800293 i,j,gk1 7 10 0.000000i,j,gk1 7 11 0.000000i,j,gk1 7 12 0.000000i,j,gk1 8 1 0.000000i,j,gk1 8 2 -305200.000000i,j,gk1 8 4 0.000000i,j,gk1 8 5 0.000000i,j,gk1 8 6 0.000000i,j,gk1 8 7 0.000000i,j,gk1 8 8 305200.000000 i,j,gk1 8 9 0.000000i,j,gk1 8 10 0.000000i,j,gk1 8 11 0.000000i,j,gk1 8 12 0.000000i,j,gk1 9 1 7564.800293i,j,gk1 9 2 0.000000i,j,gk1 9 3 12608.000977i,j,gk1 9 4 0.000000i,j,gk1 9 5 0.000000i,j,gk1 9 6 0.000000i,j,gk1 9 7 -7564.800293i,j,gk1 9 8 0.000000i,j,gk1 9 9 25216.001953i,j,gk1 9 10 0.000000i,j,gk1 9 11 0.000000i,j,gk1 9 12 0.000000i,j,gk1 10 1 0.000000i,j,gk1 10 2 0.000000i,j,gk1 10 3 0.000000i,j,gk1 10 4 -75979.257813 i,j,gk1 10 5 115892.476563 i,j,gk1 10 6 -4520.158203 i,j,gk1 10 7 0.000000i,j,gk1 10 8 0.000000i,j,gk1 10 9 0.000000i,j,gk1 10 10 75979.257813 i,j,gk1 10 11 -115892.476563 i,j,gk1 10 12 -4520.158203 i,j,gk1 11 1 0.000000i,j,gk1 11 2 0.000000i,j,gk1 11 3 0.000000i,j,gk1 11 4 115892.476563 i,j,gk1 11 5 -182890.046875 i,j,gk1 11 6 -2892.901367 i,j,gk1 11 7 0.000000i,j,gk1 11 8 0.000000i,j,gk1 11 9 0.000000i,j,gk1 11 10 -115892.476563i,j,gk1 11 11 182890.046875i,j,gk1 11 12 -2892.901367i,j,gk1 12 1 0.000000i,j,gk1 12 2 0.000000i,j,gk1 12 3 0.000000i,j,gk1 12 4 4520.158203i,j,gk1 12 5 2892.901367i,j,gk1 12 6 10619.358398i,j,gk1 12 7 0.000000i,j,gk1 12 8 0.000000i,j,gk1 12 9 0.000000i,j,gk1 12 10 -4520.158203i,j,gk1 12 11 -2892.901367i,j,gk1 12 12 21238.716797node x-disp y-disp thita1 -0.020768 -0.000749 -0.0041792 -0.021164 -0.014385 0.0078233 -0.000000 -0.000000 0.0000004 0.000000 -0.000000 0.000000reaction nodal forces from the equationsnode x-load y-load moment1 0.249791 -191.991058 -204.7498172 0.253289 -191.827469 204.7061163 94.456879 228.500870 -209.7956244 -94.456657 155.499146 -54.196941element axi-f shear-q moment-m1 94.456924 228.500854 262.488831 -94.456924 155.499146 -28.8833622 228.500870 -94.456879 -209.795624 -228.500870 94.456879 -262.4888313 181.889847 -4.264181 28.883371 -181.889847 4.264181 -54.196941。
平面桁架有限元C#编程
1题目结构如图所示： 杆的弹性模量E 为200000Mpa ，横截面面积A 为3250mm 2。
图 1 桁架示意图2实验材料PC 机一台，Microsoft Visual Studio 软件，Ansys 软件。
3实验原理（1）桁架结构特点桁架结构中的桁架指的是桁架梁，是格一种梁式结构。
桁架结构常用于大跨度的厂房、展览馆、体育馆和桥梁等公共建筑中。
由于大多用于建筑的屋盖结构，桁架通常也被称作屋架。
结构上由光滑铰链连接，载荷只作用于节点处，只约束线位移，各杆只有轴向拉伸和压缩。
（2）平面桁架有限元分析1、单元分析局部坐标系中的干单元如图所示：图 2 局部坐标系中的杆单元以下公式描述了整体位移和局部位移之间的关系：U=Tu 其中U=[ U ix U iy U jx U jy ]，T=[cos θ−sin θ00sin θcos θ0000cos θ−sin θ00sin θcos θ],u=[u ix u iy u jx u jy ]U 和u 分别代表整体坐标系和局部坐标系XY 系和局部坐标系xy 下节点i 和节点j 的位移。
T 是变形从局部坐标转换到整体坐标系下的变换阵，类似的局部力和整体力也有以下关系：F=Tf其中F=[ F ixF iy F jx F jy ] ，是整体坐标系下施加在节点i 和j 上的力的分量而且f=[ f ix f iy f jx f jy ]，代表局部坐标系下施加在节点i和j上的分量。
在假设的二力杆条件下，杆只能沿着局部坐标系的x方向变形，内力也总是沿着局部坐标系x的方向，因此将y方向的位移设置为0，局部坐标系下内力和位移通过刚度矩阵有如下关系：[f ixf iyf jxf jy]=|k0−k00000−k0k00000|=[U ixU iyU jxU jy]这里k=k eq=AE/L，写成矩阵形式有：f=Ku将f和u替换成F和U有：T-1F=KT-1U将方程两边乘以T得到：F=TKT-1U其中T-1是变换矩阵T的逆矩阵，替换方程中的TKT-1和U矩阵的值，相乘后得到：[F ixF iy F jx F jy]= k[cos2θsinθcosθ−cos2θ−sinθcosθsinθcosθsin2θ−sinθcosθ−sin2θ−cos2θ−sinθcosθcos2θsinθcosθ−sinθcosθ−sin2θsinθcosθsin2θ][U ixU iyU jxU jy]上述方程代表了施加外力、单元刚度矩阵和任意单元节点的整体位移之间的关系。
有限单元法的编程实现，内附Fortran源代码
作业2：有限单元法的编程实现2.1有限元概述有限元分析的基本概念是将求解域离散为若干子域，并通过它们边界上的节点相互联结成为组合体，对每一单元假定一个合适的近似解，然后推导求解这个域总的满足条件，从而得到问题的解。
这个解只是近似解，因为实际问题被较简单的问题所代替。
由于大多数实际问题难以得到准确解，而有限元不仅计算精度高，而且能适应各种复杂形状，因而成为行之有效的工程分析手段。
对于不同物理性质和数学模型的问题，有限元求解法的基本步骤是相同的，只是具体公式推导和运算求解不同。
有限元求解问题的基本步骤通常为：2.1.1确定求解域及求解域的离散化，即子域的剖分根据实际问题近似确定求解域的物理性质和几何区域。
将求解域近似为具有不同有限大小和形状且彼此相连的有限个单元组成的离散域，习惯上称为有限元网络划分，求解域的离散化是有限元法的核心技术之一。
选择单元的形式，确定单元数，节点数及单元及节点的编号。
2.1.2单元分析一个具体的物理问题通常可以用一组包含问题状态变量边界条件的微分方程式表示，为适合有限元求解，通常将微分方程化为等价的泛函形式。
作为弹性力学微分方程的等效积分的形式，虚位移原理与虚应力原理分别是平衡方程与力的边界条件和几何方程与位移边界条件的等效积分形式。
在导出它们的的过程中都未涉及到物理方程所以它们适用于线弹性、非线性线弹性及弹塑性的问题。
现有的有限元计算多采用以位移为未知量的形式，可用虚位移原理来描述其平衡方程，其矩阵形式为：(u )u 0T T T V S f dV TdS σδεσδδ--=⎰⎰ (2.1) 我们需要对单元构造一个适合的近似解，即推导有限单元的格式，其中包括选择合理的单元坐标系，建立单元试函数，以某种方法给出单元各状态变量的离散关系，从而形成单元矩阵（结构力学中称刚度阵或柔度阵）。
这里以平面问题为例，单元位移与节点位移的关系表示为：[,,...]e i e e i i i j j a u u N a N N a Na ⎧⎫⎪⎪====⎨⎬⎪⎪⎩⎭∑ (2.2) 应变与节点位移的关系为：[][]e e e e i j i j Lu LNa L N N a B B a Ba ε===== (2.3)应力与节点位移的关系为：e e D DBa Sa σε=== (2.4)单元刚度矩阵的表达式为：T ee K B DBtdxdy Ω=⎰ (2.5) 单元等效节点载荷矩阵表示为T T e ee e ef s P P P N ftdxdy N TtdS ΩΩ=+=+⎰⎰ (2.6) 2.1.3 整体分析将单元总装形成离散域的总矩阵方程（联合方程组），反映对近似求解域的离散域的要求，即单元函数的连续性要满足一定的连续条件。
面向对象有限元程序框架设计及其在VC_中的实现
* 国家自 然 科 学 基 金 项 目 ( 90815010 ) ; 甘 肃 省 建 设 厅 科 技 项 目 ( JK 2009- 38) 。 第一作者: 罗加成, 男, 1982 年出生, 硕士研究生。 E - m ail: h uan vi 629@ 163. com 收稿日期: 2011- 04- 15
受复杂多样的单元类型、日益演化的本构模型、 繁琐的求解过程等诸多因素的限制, 有限元程序的 复杂性显而易见, 程序设计难以一蹴而就。因此, 构 建一个合理的有限元程序框架成为开发有限元软件 的重中之重, 合理的有限元程序框架体现了设计人 员对整个程序结构的宏观把握。然而, 有限元方法 应用领域及其理论的不断拓展, 对于程序的扩充性
面向对象技术的基本思想在于高度的数据抽象 和封装能力, 将其应用于有限元程序分析与设计, 改 变了以往面向过程有限元程序的结构, 对程序的扩 展性和代码重用性 产生了巨大的 推动作用。近年 来, 研究人员逐渐认识到面向对象技术应用于有限 元程序设计的优势, 使得面向对象有限元程序分析 与设计成为一种流行趋势[ 1-3] 。
2 面向对象的程序设计 2 1 分层框架设计
有限元分析中涉及到大量的数据, 以构建有限 元数据结构为中心, 可以方便地设计有限元程序框 架。有限元分析所需的主要数据包括: 1) 描述有限 元分析的整体数据, 如单元总数、节点总数、维数、材 料的种类数等; 2) 节点数据, 包括节点坐标、节点自 由度、节点力、节点位移等; 3) 单元数据, 包括单元类 型、单元属性数据、单元材料数据、高斯积分点数据 等; 4) 材料数据, 包括弹性模量、泊松比、容重等; 5) 荷载数据, 包括体力、集中力、面力、动力荷载; 6) 约 束数据, 包括约束节点、约束位移等。
有限元方法编程
有限元方法编程【原创版3篇】篇1 目录1.有限元方法概述2.有限元方法的编程步骤3.有限元方法的应用实例4.有限元方法的优缺点篇1正文一、有限元方法概述有限元方法是一种数值分析方法，它通过将待求解的连续体划分为有限个小的、简单的子区域（即有限元），从而将连续体问题转化为有限元上的离散问题。
这种方法可以大大简化问题的求解过程，并可以在计算机上进行高效的数值计算。
有限元方法被广泛应用于固体力学、流体力学、热传导、电磁场等领域。
二、有限元方法的编程步骤1.几何建模：首先需要对问题进行几何建模，即将问题的实际物理区域抽象为计算机可以处理的几何形状。
2.网格划分：将几何模型划分为有限个小的、简单的子区域，即有限元。
这一步需要考虑网格的密度和网格的类型，以保证求解的精度和效率。
3.选择合适的有限元公式：根据问题的性质和求解的目标，选择合适的有限元公式来描述问题的物理过程。
4.编写或选用求解器：根据所选公式，编写或选用相应的求解器，进行数值计算。
5.后处理：对计算结果进行处理，包括结果的可视化和结果的解析等。
三、有限元方法的应用实例有限元方法被广泛应用于各种工程问题中，例如飞机翼的强度分析、汽车底盘的振动分析、建筑物的抗震分析等。
四、有限元方法的优缺点优点：1.可以大大简化问题的求解过程，提高求解效率。
2.可以在计算机上进行高效的数值计算，便于进行结果的可视化和解析。
3.可以适用于各种复杂的几何形状和物理过程。
缺点：1.需要进行几何建模和网格划分，这需要耗费一定的时间和精力。
2.网格的选取对求解结果的精度和效率有重要影响，需要进行适当的选择。
篇2 目录1.有限元方法概述2.有限元方法编程的基本步骤3.有限元方法编程的实际应用4.有限元方法编程的挑战与未来发展篇2正文一、有限元方法概述有限元方法是一种数值分析方法，广泛应用于固体力学、流体力学、热传导等领域。
它的基本思想是将待求解的连续体划分为有限个小的、简单的子区域，即单元，然后用有限个简单的基本函数来近似描述每个单元的物理特性。
《有限元程序设计》课件
有限元程序设计的前景展望
广泛应用
随着计算机技术的不断发展，有 限元程序设计将在更多领域得到 广泛应用，为工程设计和科学研 究提供有力支持。
技术创新
未来有限元程序设计将不断涌现 出新的技术和方法，推动该领域 不断发展壮大。
国际化发展
随着国际化交流的加强，有限元 程序设计将实现国际化发展，推 动国际合作和共同进步。
求解
求解整体方程组得到近似解。
有限元方法的应用领域
01
02
03
04
结构力学
用于分析各种结构的力学行为 ，如桥梁、建筑、机械零件等
。
流体动力学
用于模拟流体在各种介质中的 流动行为，如流体动力学、渗
流等。
热传导
用于分析温度场在各种介质中 的分布和变化。
电磁场
用于分析电磁场在各种介质中 的分布和变化，如电磁场、电
磁波等。
02
有限元程序设计的关键技术
网格生成技术
网格生成技术是有限元分析中 的重要步骤，它涉及到将连续 的物理空间离散化为有限个小 的单元，以便进行数值计算。
网格的生成需要满足一定的规 则和条件，以保证计算的精度
和稳定性。
常见的网格生成方法包括结构 化网格、非结构化网格和自适 应网格等。
网格生成技术需要考虑的问题 包括网格大小、形状、方向和 连接方式等。
02
详细描述
弹性地基板的有限元分析是一 个二维问题，需要考虑复杂的 边界条件和非线性方程的求解 。通过将地基板划分为若干个 四边形单元，可以建立非线性 方程组进行求解。
03
计算过程
04
首先将地基板划分为若干个四边 形单元，然后根据每个单元的物 理性质和边界条件建立非线性方 程组。最后通过迭代方法求解非 线性方程组得到每个节点的位移 和应力。
Fortran语言编写及有限元结构程序
算例一计算简图及结果输出用平面刚架静力计算程序下图结构的内力。
各杆EA，EI相同。
已知：642EA=4.010KN,EI=1.610KN m⨯⨯•计算简图如下：（1）输入原始数据控制参数3,5,8,7,1,2（NE,NJ,N,NW,NPJ,NPF）结点坐标集结点未知量编号0.0,0.0,0,0,0 0.0,4.0,1,2,3 0.0,4.0,1,2,4 4.0,4.0,5,6,7 4.0,0.0,0,0,8单元杆端结点编号及单元EA、EI 1,2,4.0E+06,1.6E+04 3,4,4.0E+06,1.6E+04 5,4,4.0E+06,1.6E+04结点荷载7.0,-15.0非结点荷载1.0,2.0,2.0,-18.02.0,1.0,4.0,-25.0（2）输出结果NE= 3 NJ= 5 N= 8 NW= 7 NPJ= 1 NPF= 2 NODE X Y XX YY ZZ1 0.0000 0.0000 0 0 02 0.0000 4.0000 1 2 33 0.0000 4.0000 1 2 44 4.0000 4.000056 75 4.0000 0.0000 0 0 8ELEMENT NODE-I NODE-J EA EI1 12 0.400000E+07 0.160000E+052 3 4 0.400000E+07 0.160000E+053 54 0.400000E+07 0.160000E+05CODE PX-PY-PM7. -15.0000ELEMENT IND A Q1. 2. 2.0000 -18.00002. 1. 4.0000 -25.0000NODE U V CETA1 0.000000E+00 0.000000E+00 0.000000E+002 -0.221743E-02 -0.464619E-04 -0.139404E-023 -0.221743E-02 -0.464619E-04 0.357876E-024 -0.222472E-02 -0.535381E-04 -0.298554E-025 0.000000E+00 0.000000E+00 0.658499E-03ELEMENT N Q M1 N1= 46.4619 Q1= 10.7119 M1= -6.8477N2= -46.4619 Q2= 7.2881 M2= 0.00002 N1= 7.2881 Q1= 46.4619 M1= 0.0000N2= -7.2881 Q2= 53.5381 M2= 14.15233 N1= 53.5381 Q1= 7.2881 M1= 0.0000N2= -53.5381 Q2= -7.2881 M2= -29.1523算例二计算简图及结果输出用平面刚架静力计算程序下图结构的内力。
面向对象有限元分析程序设计及其VC++实现说明书
文章编号:1000_0887(2002)12_1283_06面向对象有限元分析程序设计及其V C++实现X马永其 冯 伟(上海大学,上海市应用数学和力学研究所,上海200072)(张禄坤推荐)摘要: 采用面向对象的方法确定了有限元分析过程的对象、对象间的关系及类库,并使用VC++语言,利用其MFC 类库实现了有限元分析类库及相应窗口图形化界面的程序体序# 程序系统具有良好的可靠性,可再用性和可扩充性# 为进一步开发大型、通用、功能性强的面向对象的有限元分析软件提供参考#关 键 词: 面向对象; 有限元; 程序; 设计; VC++中图分类号: TP311;O241 文献标识码: A引 言有限元分析程序涉及力学、应用数学和计算机科学3个不同学科的理论和方法,因而其编制工作十分复杂,且程序庞大易错# 因此,有限元软件的开发一直寻求从方法上得以改善# 面向对象程序设计方法的出现,给有限元软件的开发提供了新的途径# 国外学者于1990年开始应用面向对象方法进行有限元分析程序的设计[1,2]# 经过10余年的发展,目前已编制出能够进行非线性及流体问题有限元分析,并与CAD 或数据库等其它软件集成的有限元分析程序[3,4,5],但使用VC++语言的少见报道# 国内在此领域的研究目前正在起步,仅有个别文献涉足[6,7]# 本文将通过研究设计及编程实践,讨论应用面向对象的程序设计方法进行有限元程序设计的基本思想,及采用面向对象的VC++语言进行有限元分析程序编制的基本方法# 希望引起更多国内学者的兴趣和参与,进而开发出大型、通用、功能性强的面向对象的有限元工具软件#1 面向对象方法及VC++语言面向对象程序设计方法是计算机程序开发方法的一种变革# 其主要特征是:信息隐藏性;便于在概念上体现并行和分布式结构;便于软件的演化和扩充;增加程序设计的灵活性和可理解性[8]# 面向对象程序设计过程在形式上为认定类和组织类之间关系的过程,而实质为数据抽象和代码重用# 通常归结为3个步骤:识别对象;识别对象间的关系;识别类#1283应用数学和力学,第23卷第12期(2002年12月) Applied Mathematics and Mechanics应用数学和力学编委会编重庆出版社出版 X 收稿日期: 2001_10_12;修订日期: 2002_05_21作者简介: 马永其(1966)),男,宁夏人,博士(E _mail:mayongqi@)#1284马永其冯伟Microsoft Visual C++是建立在Windows95和WindowsNT32位应用程序上的强有力的复杂开发工具,是一个Windows集成开发环境(IDE:integrated development environment)[9]#其最大的优点是节约时间,不必编译、链接便可看到应用程序的式样,可用键盘或鼠标来完成可视化设计及应用程序的许多编程工作#因此它允许在较短的时间内开发很复杂的Windows程序#面向对象的编程基础是类,Visual C++中的MFC(Mic rosoft Foundation Class)类库提供了强大的面向对象编程接口,简化了可在屏幕上创建和操纵窗口、处理绘图、使程序连接到帮助文件、便利线程、管理内存等的复杂的API(Application Programming Interface)函数调用,使应用程序的编制可以用众多的可视化控件来实现#本文的程序是在Visual C++ 5.0环境中开发编制的#2程序的设计在有限元分析过程中,主要应用了结构、载荷、节点、单元、自由度、矩阵、材料、高斯积分点、边界条件、求解和辅助计算等物理概念#因此,根据面向对象程序设计方法可确定有限元分析过程的对象为:结构对象、载荷对象、节点对象、单元对象、自由度对象、矩阵对象、材料对象、高斯积分点对象、边界条件对象、求解对象和辅助计算对象等#单元是有限元分析过程中最为重要的物理概念,它是整个有限元分析过程的核心,贯穿于整个过程的始终#因此,在标识对象间关联的过程中,单元对象确定为最关键和重要的对象#其余对象,根据其与单元对象的关系可基本分为两类对象,一类是服务于单元对象的对象,另一类是被单元对象服务的对象#在两类对象的每一类中,根据所起作用的不同又可分为包容和利用两种#其中载荷对象、节点对象、材料对象、矩阵对象、高斯积分点对象是服务于单元对象的对象,而载荷对象、节点对象和材料对象是被单元对象包容的对象,矩阵对象、高斯积分点对象是被单元对象利用的对象#结构对象、求解对象和辅助计算对象是被单元对象服务的对象,结构对象包容单元对象,求解对象和辅助计算对象利用单元对象#边界条件对象、自由度对象间接作用于单元对象,边界条件对象被结构对象包容并被单元对象间接利用#自由度对象被节点对象包容并被单元对象间接包容#根据确定的有限元分析过程的对象和所标识的对象间的关联,便形成了一个由单元类、节点类、自由度类、载荷类、材料类、边界条件类、结构类、求解类以及矩阵类和高斯节点类等组成的有限元分析类库#具体设计以单元类为例,单元类要实现的任务是:1)单元材料系数矩阵的计算、位移_应变矩阵的计算、单元刚度矩阵的计算和单元载荷列阵的计算;2)单元内节点的管理#其成员数据主要有载荷、材料物性、高斯积分点、单元应变矩阵对象数组、单元刚度矩阵、单元载荷列阵和节点对象数组等#在对具体结构进行有限元分析时,将采用具体单元类型,如一维单元、二维单元、三维单元等#因此,视单元类为虚基类,采用protec ted关键字以便其子类能够存取有关数据,采用virtual关键字以实现多态性#这样就构筑了类之间的层次和体系结构,形成了继承关系#单元一旦确定后,利用单元的自然坐标可直接建立不同类型单元的插值函数,并可计算出相应的Jac obi矩阵与其行列式#因此,具体单元类除了继承单元类的所有数据成员和数据函数外,还要增加两个成员数据即插值函数矩阵和Jacobi矩阵与其行列式,成员函数也需增加操作这两个成员数据的相关函数,以便计算具体单元的刚度矩阵#设计程序的系统控制如图1所示#结构类接受用户界面输入的材料、载荷、边界条件和节点位置等信息后,发送消息并构造具体单元类、节点类、自由度类#然后由结构类、具体单元类、节点类和自由度类向求解类通讯并发送消息# 最后由结构类、具体单元类、节点类和求解类向辅助计算类通讯并发送消息# 从而完成一个结构的有限元分析过程#图1 类之间的关系及程序控制3 系统的实现根据以上有限元分析过程的面向对象程序设计,用VC++语言具体描述类库,以单元类为例,它是Visual C++的MFC 类库中的C Object 类的继承类,具体表示如下:class CElement:public CObject{protected:int noelement; //element numberint numnodes;//number of element nodes CMaterial *elmaterial;//pointer to element material CLoad *elload;//pointer to element load CGausspoint *gausspoint;//pointer to gausspoint class CMatrix *Dmat;//pointer to element constitutive matrix CObArray BmatArray;//a rray of e lement strain_displacement matrix CMatrix *kemat;//pointer to element stiffness matrix CMatrix *fevec;//pointer to element load ma trix int *lma rray;//pointer to equation numbe r of dof of element .s node CObArray nodearray;//a rray of e lement nodespublic:CElement();//construc tor CElement(int);//construc tor virtual ~CElement();//destructor 1285面向对象有限元分析程序设计及其VC++实现1286马永其冯伟virtual void setDma t()=0;//sets the element constitutive ma trixvirtual void setBmatArra y()=0;virtual void setkemat()=0;virtual void setfevec()=0;virtual void setnodearray()=0;virtual CMatrix*getDma t()=0;//returns the element constitutive matrixvirtual CObArray*getBmatArray()=0;virtual CMatrix*getkemat()=0;virtual CMatrix*getfe vec()=0;Virtual CObArray*getnodearray()=0;//,(other member data and functions)};利用Visual C++的App Wizard工具生成多文档界面(MDI)有限元分析系统的WINDWS应用程序框架,共生成6个类:About对话框的对话框类(C AboutDlg)、整个应用程序的CWinApp 类(C MyApp)、文档类(CMyDos)、视图类(CMyVie w)、框架类(CMainFra me)和子框架类(CChild-Frame)#加入有限元分析过程类库,并将有限元分析过程归纳、抽象、定义为一个问题类,由有限元分析的结构类、求解类和应力类组合而成,其成员数据为结构类对象、求解类对象和应力类对象#将问题类作为WI NDOWS应用程序框架中文档类的公共类型属性变量,加入其中,具体问题类和文档类属性变量描述如下:class CProblem:public CObject{private:CString name;//string of problem nameCStructure*structure;//pointer to structureCLinea rsolver*solver;//pointer to solverCStress*stress;//pointer to stresspublic:CProblem();//constructorCProblem(CString);//constructorvirtual~CProblem()://destructorvoid setT(CString);//sets the problem nameCString getT();//returns the problem namevoid setstructure(CStructure*);//sets the struc tureCStructure*getstruc ture();//returns the structurevoid setsolver(CLinearsolver*);//sets the solverCLinea rsolver*getsolver();//returns the solvervoid setstress(CStress*);//set the stressCStress*getstress();//returns the stress//,(other member func tions)};class CMyDoc:public CDocum ent{protected://c reate from serializa tion onlyCMyDoc();DECLARE -DYNCREATE(CMyDoc )//Attributespublic:CProblem m -problem;//attribute of CMyDoc problem //,(other member date and functions)};利用WI NDOWS 操作系统的消息循环机制,建立起相应的系统通讯体系,其作用方式通过图2所示#由MFC 类库中提供的丰富的for WI NDOW S 的窗口、菜单、对话框、属性表等类库,通过继承、修改和扩充编制所需的人机对话界面,如数据和分析等菜单项,结构对话框类(CStructure -Dialog)等窗口图形化界面#当用户要完成对一个结构进行有限元应力分析的计算过程时,只要将鼠标点到菜单(界面层)上,即发出消息将/菜单对话0激活,菜单对象产生反应# 通过完成对话框界面的人机交互后,发出消息至WI NDOW S 操作系统的消息队列中参加消息循环# /问题对象0(有限元分析层)获此消息后,产生相应的计算动作,完成有限元分析的计算过程#图2 系统消息流图4 小 结通过对有限元分析过程进行面向对象的程序设计,确立了相应的对象及对象间的关系,并建立了相应的类库,说明面向对象方法是一种易于用来进行科学和工程计算程序设计的强有力的程序设计方法# 采用VC++语言描述了有限元分析类库,通过继承、修改和扩充,利用MFC 类库方便地建立了面向对象的有限元分析及相应窗口图形化界面的程序体系,即保证了有限元分析类库的灵活使用,又保证了系统操作的清晰简便,从而提高了程序的可靠性,可再用性和可扩充性#本文的开发实践将为进一步开发大型、通用、功能性强的面向对象的有限元分析软件提供有价值参考#1287面向对象有限元分析程序设计及其VC++实现1288马永其冯伟[参考文献][1]Forde B W R,Foschi R O,Stiemer S F.Objec t_oriented finite elem ent analysis[J].Computer s andStr uctur es,1990,34(3):355)374.[2]Fenves G F.Object_oriented programming for engineering software development[J].Engin eer in gWith Com puter s,1990,6(1):1)15.[3]Archer G C,Fenves G,Thewalt C.A new objec t_oriented finite element analysis[J].Com puter s andStru ctur es,1999,70(1):63)75.[4]Robe rt Ian Mackie.An object_oriented approach to ca lculation control in finite element prog rams[J].Com puter s a nd Stru ctur es,2000,77(5):461)474.[5]YU Li_chao,Kuma r Ashok V.An object_oriented modular framework for implementing the finite ele-ment method[J].Computer s and Str uctur es,2001,79(9):919)928.[6]孔祥安.C++语言与面向对象有限元程序设计[M].成都:西南交通大学出版社,1995.[7]马永其,陈罕,李斯特.采用面向对象方法的有限元程序[J].北京化工大学学报,2000,27(3):51)55.[8]宛延凯.C++语言和面向对象程序设计[M].(第二版).北京:清华大学出版社,1998.[9]Kate Gregory.Visual C++5开发使用手册[M].康博创作室译,文珍审校.北京:机械工业出版社,1998.Object_Oriented Finite Eleme nt Analysisand Programming in VC++MA Yong_qi FENG W ei(Sha ngha i In stitut e of Applied Mathem atics and Mecha nics,Shan gha i Univer sity Shanghai200072,P R China)Abs tract:The design of finite element analysis program using objec t_oriented programm ing(OOP) techniques is presented.The objects,classes and the subclasses used in the programm ing are ex-plained.The system of classes library of finite element ana lysis program and Windows_type Graphica l User Interfaces by VC++and its MFC are developed.The reliability,reusability and extensibility of program are enhanced.It is a reference to develop the large_scale,versatile and powerful systems of object_oriented finite elem ent software.Key wo rds:object_oriented program ming;finite elem ent method;program;design;VC++。
有限元编程算例
有限元编程算例（Fortran）本程序通过Fortran 语言编写，程序在In tel Parallel Studio XE 2013 withVS2013中成功运行，程序为《计算力学》（龙述尧等编）一书中的源程序，仅作研究学习使用，省去了敲写的麻烦。
3.7.4算例例工夕设谦聲祇受沟布戴荷、如图3-3叙叮所示「假進£=」，泊松比世二①17•不计容重+厚度才=1 为平面应力问題，因对称取半边结构计真•結枸支承，跟兀划分，节点绸号如图3,16（b）JPfzii p试矗出及y = S m截啲的竖向位移图応=总mflg面的靳应力分布图，⑴（山图玉開燮均布栽荷的简支漂梁源程序为：!Page149C0MM0N/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,E0,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)0PEN(5,FILE='DATAIN')!OPEN(6,FILE='DATAOUT',STATUS='NEW')CALL DATAIF(IND.EQ.0)GOTO 10EO=EO/(1.0-UN*UN)UN=UN/(1.0-UN)10 CALL TOTSTICALL LOADCALL SUPPORCALL SOLVEQCALL STRESSPAUSE!STOPENDSUBROUTINE DATACOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)READ(5,*)NJ,NE,NZ,NDD,NPJ,INDNJ2=NJ*2NPJ1=NPJ+1READ(5,*)EO,UN,GAMA,TEREAD(5,*)((JM(I,J),J=1,3),I=1,NE)READ(5,*)((CJZ(I,J),J = 1,2),I=1,NJ)!Page150READ(5,*)(NZC(I),I=1,NZ)READ(5,*)((PJ(I,J),J=1,2),I=1,NPJ1)WRITE(6,10)(I,(CJZ(I,J),J=1,2),I=1,NJ)10 FORMA T(4X,2HNO,6X,1HX,6X,1HY/(I6,2X,F7.2,F7.2))RETURNENDSUBROUTINE ELEST(MEO,IASK)COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)IE=JM(MEO,1)JE=JM(MEO,2)ME=JM(MEO,3)CM=CJZ(JE,1)-CJZ(IE,1)BM=CJZ(IE,2)-CJZ(JE,2)CJ=CJZ(IE,1)-CJZ(ME,1)BJ=CJZ(ME,2)-CJZ(IE,2)AE=(BJ*CM-BM*CJ)/2.0IF(IASK.LE.1) GOTO 50DO 10 I=1,3DO 10 J=1,6B(I,J)=0.010 CONTINUEB(1,1)=-BJ-BMB(1,3)=BJB(1,5)=BMB(2,2)=-CJ-CMB(2,4)=CJB(2,6)=CMB(3,1)=B(2,2)B(3,2)=B(1,1)B(3,3)=B(2,4)B(3,4)=B(1,3)B(3,5)=B(2,6)!Page151B(3,6)=B(1,5)DO 20 I=1,3DO 20 J=1,6B(I,J)=B(I,J)/(2.0*AE)20 CONTINUED(1,1)=EO/(1.0-UN*UN)D(1,2)=EO*UN/(1.0-UN*UN)D(2,1)=D(1,2)D(2,2)=D(1,1)D(1,3)=0.0D(2,3)=0.0D(3,1)=0.0D(3,2)=0.0D(3,3)=EO/(2.0*(1.0+UN))DO 30 I=1,3DO 30 J=1,6S(I,J)=0.0DO 30 K=1,3S(I,J)=S(I,J)+D(I,K)*B(K,J) 30 CONTINUEIF(IASK.LE.2) GOTO 50DO 40 I=1,6DO 40 J=1,6EKE(I,J)=0.0DO 40 K=1,3!**********************************Exchange BAnd S***********************************************EKE(I,J)=EKE(I,J)+B(K,I)*S(K,J)*AE*TE40 CONTINUE50 CONTINUERETURNENDSUBROUTINE TOTSTICOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200) !Page152 DO 20 I=1,NJ2DO 20 J=1,NDDTKZ(I,J)=0.020 CONTINUE!*************Not Understanded*****************************DO 30 MEO=1,NECALL ELEST(MEO,3)DO 30 I=1,3DO 30 II=1,2LH=2*(I-1)+IILDH=2*(JM(MEO,I)-1)+IIDO 30 J=1,3DO 30 JJ=1,2L=2*(J-1)+JJLZ=2*(JM(MEO,J)-1)+JJLD=LZ-LDH+1IF(LD.LE.0) GOTO 30TKZ(LDH,LD)=TKZ(LDH,LD)+EKE(LH,L)30 CONTINUERETURNENDSUBROUTINE LOADCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)DO 10 I=1,NJ2P(I)=0.010 CONTINUEIF(NPJ.EQ.0) GOTO 30DO 20 I=1,NPJI1=I+1J=IFIX(PJ(I1,2))P(J)=PJ(I1,1)20 CONTINUE30 IF(GAMA.LE.0.0) GOTO 50!Page153DO 40 MEO=1,NECALL ELEST(MEO,1)PE=-GAMA*AE*TE/3.0IE=JM(MEO,1)JE=JM(MEO,2)ME=JM(MEO,3)P(2*IE)=P(2*IE)+PEP(2*JE)=P(2*JE)+PEP(2*ME)=P(2*ME)+PE40 CONTINUE 50 CONTINUERETURNENDSUBROUTINE SUPPORCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)DO 60 I=1,NZMZ=NZC(I)TKZ(MZ,1)=1.0DO 10 J=2,NDDTKZ(MZ,J)=0.010 CONTINUEIF(MZ-NDD)20,20,3020 JO=MZGOTO 4030 JO=NDD40 DO 50 J = 2,JOJ1=MZ-JTKZ(J1+1,J)=0.050 CONTINUEP(MZ)=0.060 CONTINUERETURNEND!Page154SUBROUTINE SOLVEQCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)NJ1=NJ2-1DO 50 K=1,NJ1IF(NJ2-K-NDD+1)10,10,2010 IM=NJ2GOTO 3020 IM=K+NDD-130 K1=K+1DO 50 I=K1,IML=I-K+1C=TKZ(K,L)/TKZ(K,1)LD1=NDD-L+1DO 40 J=1,LD1 M=J+I-K TKZ(I,J)=TKZ(I,J)-C*TKZ(K,M)40 CONTINUEP(I)=P(I)-C*P(K)50 CONTINUEP(NJ2)=P(NJ2)/TKZ(NJ2,1)DO 100 I1 = 1,NJ1I=NJ2-I1 !************************************************************************IF(NDD-NJ2+I-1)60,60,70 60 JO=NDD下面一行可能出错GOTO 8070 JO=NJ2-I+180 DO 90 J=2,JOLH=J+I-1P(I)=P(I)-TKZ(I,J)*P(LH)90 CONTINUEP(I)=P(I)/TKZ(I,1)100 CONTINUE!Page155WRITE(6,110)(I,P(2*I-1),P(2*I),I=1,NJ)!************************************************************************************110 FORMAT(2X,3HJD=,3X,2HU=,12X,2HV=/(I4,3X,F16.7,3X,F16.7))RETURNENDSUBROUTINE STRESSCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200) DIMENSION WY(6),YL(3)DO 60 MEO=1,NECALL ELEST(MEO,2)DO 10 I=1,3DO 10 J=1,2LH=2*(I-1)+JLDH=2*(JM(MEO,I)-1)+JWY(LH)=P(LDH)10 CONTINUEDO 20 I=1,3YL(I)=0.0DO 20 J=1,6YL(I)=YL(I)+S(I,J)*WY(J)20 CONTINUESIGX=YL(1)SIGY=YL(2)TOXY=YL(3)PYL=(SIGX+SIGY)/2.0SIG=(SIGX-SIGY)**2/4.0+TOXY*TOXYRYL=SQRT(SIG)SIG1=PYL+RYLSIG2=PYL-RYLIF(SIGY .EQ.SIG2) GOTO 30CETA1=TOXY/(SIGY-SIG2)CETA=90.0-57.29578*ATAN(CETA1)GOTO 40!Page15630 CETA=0.040 WRITE(6,50)MEO,SIGX,SIGY ,TOXY ,SIG1,SIG2,CETA50FORMA T(4X,2HE=,I3/2X,3HSX=,F11.3,3X,3HSY=,F11.3,3X,4HTAU=,F11.3/2X,3HS1=,F11.3,3X,3HS2=,F11. 3,3X,4HCET=,F11.3)!50FORMA T(4X,2HE=,I3/2X,3HSX=,Fll.3,3X,3HSY=,F11.3,3X,4HTAU=,F11.3/2X,3HSl=,Fll.3,3X,3HS2=,F11.3,3 X,4HCET=,F11.3)60 CONTINUERETURNEND输入文件为datain28,36,9,10,4,01,0.17,0,11,5,22,5,62,6,33,6,73,7,44,7,85,9,66,9,106,10,77,10,117,11,88,11,129,13,1010,13,1410,14,1111,14,1511,15,1212,15,1613,17,1414,17,1814,18,1515,18,1915,19,1616,19,2017,21,1818,21,2218,22,1919,22,2319,23,2020,23,2421,25,2222,25,2622,26,2323,26,2723,27,2424,27,280,61,62,63,60,51,52,53,50,41,42,43,40,31,32,33,30,21,22,23,20,11,12,13,10,01,02,03,07,15,23,31,39,47,49,50,55 0,0-5E4,2-10E4,4-10E4,6-5E4,8输出结果为：DATAOUTNO X Y1 0.00 6.002 1.00 6.003 2.00 6.004 3.00 6.005 0.00 5.006 1.00 5.007 2.00 5.008 3.00 5.009 0.00 4.0010 1.00 4.0011 2.00 4.0012 3.00 4.0013 0.00 3.0014 1.00 3.0015 2.00 3.0016 3.00 3.0017 0.00 2.0018 1.00 2.0019 2.00 2.0020 3.00 2.0021 0.00 1.0022 1.00 1.0023 2.00 1.0024 3.00 1.0025 0.00 0.0026 1.00 0.0027 2.00 0.0028 3.00 0.00JD= U= V=1 -29766.873 -1173917.7502 -14003.185 -1174018.8753 -3753.270 -1179518.1254 0.000 -1181719.7505 -26382.471 -1072681.5006 -10746.993 -1073615.0007 -2064.593 -1082360.7508 0.000 -1085873.2509 -13536.995 -964010.12510 3372.794 -970055.12511 7268.415 -989269.12512 0.000 -998401.81213 7816.581 -835383.43814 27176.234 -861713.93815 22063.230 -905726.12516 0.000 -927165.18817 29514.479 -665602.87518 53419.637 -747340.43819 34876.832 -839806.81220 0.000 -881219.12521 29580.273 -416288.71922 52944.918 -632601.12523 17504.195 -803765.68824 0.000 -859481.93825 0.000 0.00026 -120102.820 -583505.37527 -76202.375 -787347.18828 0.000 -829170.812E= 1SX= -1489.530 SY=-101489.383 TAU= -1489.531 S1= -1467.348 S2=-101511.562 CET= 179.147 E= 2SX= -1475.844 SY=-100654.875 TAU= -1790.500 S1= -1443.531 S2=-100687.188 CET= 178.966 E= 3SX= -7021.670 SY=-101597.672 TAU= -3741.688 S1= -6873.875 S2=-101745.469 CET= 177.738 E= 4SX= -8067.500 SY= -98528.750 TAU= -4459.156 S1= -7848.227 S2= -98748.023 CET= 177.185 E= 5SX= - 13143.328 SY= -99391.750 TAU= -1662.500S1= -13111.293 S2= -99423.781 CET= 178.896 E= 6SX= -14652.781 SY= -98337.500 TAU= -1501.062S1= -14625.867 S2= -98364.414 CET= 178.973E= 7SX= -2923.122 SY=-109168.297 TAU= -5888.469 S1= -2597.762 S2=-109493.656 CET= 176.837 E= 8SX= -716.078 SY=-103681.562 TAU= -8617.406 S1= 0.148 S2=-104397.789 CET= 175.249 E= 9SX= -9188.316 SY=-105121.867 TAU= -9771.594 S1= -8203.125 S2=-106107.062 CET= 174.243 E= 10SX= -12285.000 SY= -95180.250 TAU= -12199.594 S1= -10526.887 S2= -96938.359 CET= 171.799E= 11SX= -14170.516 SY= -95500.750 TAU= -5489.531 S1= -13801.664 S2= -95869.602 CET= 176.156E= 12SX= -22797.406 SY= -91347.000 TAU= -3902.844 S1= -22575.914 S2= -91568.492 CET= 176.752E= 13SX= -5104.269 SY=-129494.438 TAU= -11708.750 S1= -4011.727 S2=-130586.977 CET= 174.669 E= 14SX= 969.672 SY=-108176.375 TAU= -21424.750 S1= 5024.582 S2=-112231.281 CET= 169.283 E= 15SX= -14954.572 SY=-110883.469 TAU= - 18383.531 S1= -11552.273 S2=-114285.766 CET= 169.515E= 16SX= -19890.141 SY= -86924.312 TAU= -25131.188 S1= -11514.844 S2= -95299.609 CET= 161.569E= 17SX= -22109.688 SY= -87301.625 TAU= - 10225.406 S1= -20543.453 S2= -88867.859 CET= 171.292E= 18SX= -35190.453 SY= -77219.000 TAU= -9162.000 S1= -33280.023 S2= -79129.430 CET= 168.222E= 19SX= -9785.850 SY=-171444.172 TAU= -20524.969 S1= -7220.594 S2=-174009.422 CET= 172.876 E= 20SX= 4594.438 SY=-113592.375 TAU= -46145.688 S1= 20477.398 S2=-129475.336 CET= 161.007 E= 21SX= -25287.307 SY=-118672.312 TAU= -30023.750 S1= -16467.512 S2=-127492.109 CET= 163.629 E= 22SX= -30634.422 SY= -71127.188 TAU= -44991.469 S1= -1543.715 S2=-100217.891 CET= 147.114 E= 23SX= -34259.609 SY= -71743.438 TAU= -14637.906 S1= -29220.699 S2= -76782.344 CET= 161.005 E= 24SX= -43958.047 SY= -53418.938 TAU= - 17697.562 S1= -30369.627 S2= -67007.359 CET= 142.482 E= 25SX= -19028.160 SY=-252549.000 TAU= -34958.688 S1= -13907.055 S2=-257670.094 CET= 171.666 E= 26SX= 3973.812 SY=-114063.750 TAU= -92238.344 S1= 54459.047 S2=-164548.984 CET= 151.307 E= 27SX= -39180.809 SY=-121400.055 TAU= -39312.688 S1= -23409.074 S2=-137171.781 CET= 158.140 E= 28SX= -42804.766 SY= -43317.938 TAU= - 65723.062 S1= 22662.211 S2=-108784.914 CET= 135.112 E= 29SX= -42224.094 SY= -43219.188 TAU= - 10273.375 S1= -32436.225 S2= -53007.055 CET= 136.386 E= 30SX= -21830.422 SY= -25448.312 TAU= -23810.344 S1= 239.594 S2= -47518.328 CET= 137.172 E= 31SX= -48815.199 SY=-424587.344 TAU= -79800.078 S1= -32570.844 S2=-440831.688 CET= 168.494 E= 32SX=-132271.750 SY= -71582.000 TAU=-175409.250 S1= 76087.781 S2=-279941.531 CET= 130.093 E= 33SX= -45090.102 SY= -56761.105 TAU= 804.781 S1= -45034.867 S2= -56816.336 CET= 3.926 E= 34SX= 42332.711 SY= -9221.938 TAU= -47066.344 S1= 70218.328 S2= -37107.555 CET= 149.354 E= 35SX= -20899.344 SY= -19971.375 TAU= 16235.219 S1= -4193.512 S2= -36677.207 CET= 45.819E= 36SX= 73163.914 SY= -17873.250 TAU= -17873.344 S1= 76547.250 S2= -21256.586 CET= 169.281。
元编程的c实现算例修订稿
元编程的c实现算例公司标准化编码 [QQX96QT-XQQB89Q8-NQQJ6Q8-MQM9N]有限元编程的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 pi16.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];255. f0[6]=(g*c*c*c)*(4-3*c/l0)/(12*l0);256. }257. else258. {259. if(ind==2) {261. f0[1]=;262. f0[2]=(-(g*d*d)*(l0+2*c))/(l0*l0*l0);263. f0[3]=-(g*c*d*d)/(l0*l0);264. f0[4]=;265. f0[5]=(-g*c*c*(l0+2*d))/(l0*l0*l0);266. f0[6]=(g*c*c*d)/(l0*l0);267. }268. else269. {270. f0[1]=-(g*d/l0); f0[2]=;272. f0[3]=;273. f0[4]=-g*c/l0;274. f0[5]=;275. f0[6]=;276. }277. }278. }279.280. void zb(int e)284. {285. double ceta,co,si;286. int i,j;287. ceta=(gj[e]*pi)/180; co=cos(ceta); 289. si=sin(ceta);290. t[1][1]=co; t[1][2]=si;292. t[2][1]=-si;293. t[2][2]=co;294. t[3][3]=;295. for(i=1;i<=3;i++)296. {297. for(j=1;j<=3;j++) {299. t[i+3][j+3]=t[i][j];300. }301. }302. }303.304.305.306. void jdugd(int e)310. {311. double A0,l0,j0;312. int i;313. int j;314.315.316. A0=mj[e]; l0=gc[e]; j0=gx[e]; 320.321. for(i=0;i<=6;i++)322. for(j=0;j<=6;j++) kd[i][j]=;324.325. kd[1][1]=e0*A0/l0;326. kd[2][2]=12*e0*j0/pow(l0,3);327. kd[3][2]=6*e0*j0/pow(l0,2);328. kd[3][3]=4*e0*j0/l0;329. kd[4][1]=-kd[1][1];330. kd[4][4]=kd[1][1];331. kd[5][2]=-kd[2][2]; kd[5][3]=-kd[3][2];333. kd[5][5]=kd[2][2];334. kd[6][2]=kd[3][2];335. kd[6][3]=2*e0*j0/l0;336. kd[6][5]=-kd[3][2];337. kd[6][6]=kd[3][3];338.339. for(i=1;i<=6;i++)340. for(j=1;j<=i;j++) kd[j][i]=kd[i][j];342. }343.344.345. void dugd(int e)349. {350. int i,k,j,m;351. jdugd(e); zb(e); for(i=1;i<=6;i++)354. {355. for(j=1;j<=6;j++)356. {357. ke[i][j]=;358. for(k=1;k<=6;k++)359. {360. for(m=1;m<=6;m++)361. {362. ke[i][j]=ke[i][j]+t[k][i]*kd[k][m]*t[m][j]; } 364. }365. }366. }367. }368.369.370.372.373.374. </></>。
一个最基本的有限元计算程序
一个最基本的有限元计算程序有限元计算程序是一种数值计算方法，用于求解结构力学中的问题。
它将结构划分为有限个小单元，通过离散化方法，将结构的连续性问题转化为在节点上的离散化问题，然后利用数值方法求解，得到结构的应力、应变、位移等结果。
以下是一个最基本的有限元计算程序的设计和实现。
1.输入：用户需要输入结构的几何形状、材料属性、载荷和边界条件等信息。
-结构几何形状：可以通过输入结构的节点坐标和单元连接关系来描述结构的几何形状。
-材料属性：包括材料的弹性模量和泊松比等参数。
-载荷：可以输入结构上的节点力、边界面上的边界条件等，同时也可以输入分布载荷。
-边界条件：可以输入结构上的固约束条件，如支撑或固定。
2.网格划分：根据输入的节点坐标和单元连接关系，将结构划分为有限个小单元。
可以选择不同的划分方法，如三角形划分或四边形划分等。
3.单元刚度矩阵计算：对每个小单元，通过单元刚度矩阵的计算来建立整个结构的刚度矩阵。
单元刚度矩阵的计算需要根据材料属性和几何形状来求解。
4.结构总刚度矩阵组装：将每个小单元的刚度矩阵组装成整个结构的总刚度矩阵。
对于重叠节点，可以根据不同的组装方法来进行。
5.边界条件的处理：根据输入的边界条件，对总刚度矩阵进行边界条件的处理，将已知位移和力的约束转化为未知的位移和力的约束。
6.方程组的求解：利用数值方法（如高斯消元法、Cholesky分解法等）求解已经处理好的约束方程组，得到未知位移。
7.结果输出：输出计算结果，包括应力、应变、位移等。
以上是一个最基本的有限元计算程序的设计和实现过程。
在实际应用中，还可以对程序进行进一步的优化和改进，提高计算效率和准确性。
平板有限元C#编程
1题目一个钢板，其板厚为t (10mm)，中间开有一个小圆孔，在钢板两侧面受均布拉力的作用，p = 6MPa 。
材料的弹性模量E = 200GPa ，泊松比μ = 0.3。
钢板长度a （100mm ），宽度50mm （b ），圆孔直径10mm （d ），位置为钢板的几何中心。
试求小孔附近区域的应力结果。
p带孔薄板2分析 2.1问题简化该模型具有对称性，为减少总体刚度矩阵，于是只取整个模型的四分之一计算。
在简化后的模型的最下方的变上，x 方向的位移关于此边对称，y 方向上的位移关于此边反对称，而在最右方的边上，x 方向位移关于此边反对称，y 方向位移关于此边对称，并且左右两边载荷对称，因此，可在简化后的模型的下方的边上垂直放置滚动铰支座表示只有x 方向位移，在最右边与x 方向垂直放置滚动铰支座，表示只有y 方向位移。
简化后的模型2.2网格划分方法本题使用线性三角形单元分析，网格划分为均匀的三角形网格。
网格划分的方法为，将x1边m 等分，y1以及y2边n 分，再将圆弧m0等分，最后将x2边（m-m0）等分，并找到各个等分点的坐标。
再将x1边的等分点按照顺序分别与x2边和圆弧边的等分点进行连接，找出各条线段的n 等分点的坐标，然后依次连接各个等分点形成4边形网格，最后连接四边形网格的一条对角线形成三角形网格。
节点和单元编号的方法为从上到下从左至右，每一列编号完成后，立即跳至下一列的顶端，向下编号，直到编号完成。
2.3单元分析本问题为平面应力问题，在z 面上没有内力，其应力状态为：Txyxy=[]对于平面应力问题，胡克定律的通用形式为：x xx y yy 2xyxy10E 1011002或者写为：，其中E 是弹性模量，是泊松比。
对于线性三角形单元，用形函数和节点位移表示的位移变量为：i ix j jx k kx u S U S U S U i iyj jyk ky vS U S U S U将上式写为矩阵形式：ix iy i j k jxij k jykx kyU U S 0S 0S 0U u 0S 0S 0S U v U U然后将应变也用位移表示，并写成矩阵形式：ix iy x ijkjxy i j kjy xyiijjkkkx kyU U 0U 100U 2AU U 其中ij kY Y jk iY Y ki jY Y ik j X X ji kX X kj i X X ，进一步将上述方程写成：BU 。
有限元编程的c++实现算例
有限元编程的c++实现算例1. #include<>2. #include<>3.4.5. #define ne 3 #define nj 4 #define nz6 #define npj 0 #define npf1 #define nj3 12 #define dd6 #define e0 #definea0 #define i0 #define pi16.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()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. }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.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];255. f0[6]=(g*c*c*c)*(4-3*c/l0)/(12*l0);256. }257. else258. {259. if(ind==2) {261. f0[1]=;262. f0[2]=(-(g*d*d)*(l0+2*c))/(l0*l0*l0);263. f0[3]=-(g*c*d*d)/(l0*l0);264. f0[4]=;265. f0[5]=(-g*c*c*(l0+2*d))/(l0*l0*l0);266. f0[6]=(g*c*c*d)/(l0*l0);267. }268. else270. f0[1]=-(g*d/l0); f0[2]=;272. f0[3]=;273. f0[4]=-g*c/l0;274. f0[5]=;275. f0[6]=;276. }277. }278. }279.280. void zb(int e)284. {285. double ceta,co,si;286. int i,j;287. ceta=(gj[e]*pi)/180; co=cos(ceta); 289. si=sin(ceta);290. t[1][1]=co; t[1][2]=si;292. t[2][1]=-si;293. t[2][2]=co;294. t[3][3]=;295. for(i=1;i<=3;i++)296. {297. for(j=1;j<=3;j++) {299. t[i+3][j+3]=t[i][j];300. }301. }302. }303.304.305.306. void jdugd(int e)310. {311. double A0,l0,j0;312. int i;314.315.316. A0=mj[e]; l0=gc[e]; j0=gx[e];320.321. for(i=0;i<=6;i++)322. for(j=0;j<=6;j++) kd[i][j]=;324.325. kd[1][1]=e0*A0/l0;326. kd[2][2]=12*e0*j0/pow(l0,3);327. kd[3][2]=6*e0*j0/pow(l0,2);328. kd[3][3]=4*e0*j0/l0;329. kd[4][1]=-kd[1][1];330. kd[4][4]=kd[1][1];331. kd[5][2]=-kd[2][2]; kd[5][3]=-kd[3][2];333. kd[5][5]=kd[2][2];334. kd[6][2]=kd[3][2];335. kd[6][3]=2*e0*j0/l0;336. kd[6][5]=-kd[3][2];337. kd[6][6]=kd[3][3];338.339. for(i=1;i<=6;i++)340. for(j=1;j<=i;j++) kd[j][i]=kd[i][j];342. }343.344.345. void dugd(int e)349. {350. int i,k,j,m;351. jdugd(e); zb(e); for(i=1;i<=6;i++)354. {355. for(j=1;j<=6;j++)356. {357. ke[i][j]=;358. for(k=1;k<=6;k++)359. {360. for(m=1;m<=6;m++)361. {362. ke[i][j]=ke[i][j]+t[k][i]*kd[k][m]*t[m][j]; } 364. }365. }366. }367. }368.369.370.372.373.374. </></>。
有限元C++编程实践
基于有限元算法的编程实践学号：2011043010031姓名：廖校毅电子科技大学物理电子学院【摘要】本文就有限元算法在电磁场分析中的应用展开研究与实践，从电磁场的Maxwell方程出发,根据电磁场的边值问题及变分公式建立了有限元方程组。
具体在实践中，将这些知识运用到C++语言和Matlab中，并将这两种语言有机结合，编程并实现二维FEM。
程序最后通过计算含两种介质的电位槽电位分布来验证其正确性。
关键词: 有限元变分法C++ MatlabThe Programming Practice Based on The Finite Element Algorithm Student ID：2011043010031Name：Liao Xiaoyi University of Electronic Science and technology &School of PhysicalElectronicsAbstract In this paper, we take the application of finite element method in electromagnetic field analysis into research and practice. Starting from Maxwell equations of electromagnetic field, the electromagnetic field boundary value problems and variational formula established the finite element equations. In specific practice, this knowledge will be applied to the C++ language and Matlab, and the organic combination of two languages, programming and implementation of two-dimensional FEM. Finally, through the program to verify the validity of the calculation of potential distribution in channel potential containing two kinds of medium.Key words The finite element method The variational method C++ Matlab一、前言在数学中，有限元法（FEM，Finite Element Method）是一种为求得偏微分方程边值问题近似解的数值技术。
有限单元法程序代码
弹性力学平面问题有限元程序#define NODE_NUM 300#define ELE_NUM 200#define NODE2 NODE_NUM*2#include<stdlib.h>#include<stdio.h>#include<math.h>int nj,ne,nm,nz,npj,jm[ELE_NUM][4],nzc[60],nj2,bw,lxm;float ea,e,um,oul,te,xy[NODE_NUM][4],EA[5],p[NODE2],pj,F[ELE_NUM],etm[5][5]; float ek[7][7];double **K,*f;int BAND2(int n,int m, double **K,double *f);double **TwoArrayAlloc(int,int);void TwoArrayFree(double **);FILE *outf;/*-------------------------------------------*/double **TwoArrayAlloc(int r,int c){double *x,**y;int n;x=(double*)calloc(r*c,sizeof(double));y=(double**)calloc(r,sizeof(double*));for(n=0;n<=r-1;++n)y[n]=&x[c*n];}/*-------------------------------------------*/void TwoArrayFree(double **x){free(x[0]);free(x);}/*------------------------------------------*/int BAND2(int n,int m, double **K,double *f) {int i,j,t,ij,ji,it,jt,tm,m1;double s,w;m1=m+1;for(i=1;i<=n;i++){if(K[i-1][m1-1]<=0) return(0);w=0.0;if(i>m1) tm=i-m;else tm=1;for(j=tm;j<=i;j++){s=0.0;ij=j-i+m1;for(t=tm;t<=j-1;t++){jt=t-j+m1;s=s+K[i-1][it-1]*K[j-1][jt-1]/K[t-1][m1-1]; }K[i-1][ij-1]=K[i-1][ij-1]-s;if(j==i)f[i-1]=f[i-1]-w;elsew=w+K[i-1][ij-1]*f[j-1]/K[j-1][m1-1];}}for(i=n;i>=1;i--){s=0.0;if(i>n-m1) tm=n; else tm=i+m;for(j=i+1;j<=tm;j++){ji=i-j+m1;s=s+K[j-1][ji-1]*f[j-1];}f[i-1]=(f[i-1]-s)/K[i-1][m1-1];}return 1;}/*-----------------------------------------*/void input(){int jj,j,i,nh,nl;float dx,dy;FILE *infile;char name1[30],name2[30];R:printf("please enter data-filename\n");scanf("%s",name1);if((infile=fopen(name1,"r"))==NULL){ printf("the data-file not exit!");goto R;}printf("Please enter result-filename\n");scanf("%s",name2);outf=fopen(name2,"w");fscanf(infile,"%d%d%d%d%d%d5d",&nj,&ne,&nm,&nz,&npj,&lxm); for(i=1;i<=nj;i++){for(j=1;j<3;j++)fscanf(infile,"%f",&xy[i][j]);}for(i=1;i<=ne;i++){for(j=0;j<4;j++)fscanf(infile,"%d",&jm[i][j]);}for(i=1;i<=nm;i++){for(j=1;j<=4;j++)fscanf(infile,"%f",&etm[i][j]);}for(i=1;i<=nz;i++)fscanf(infile,"%d",&nzc[i]);nj2=nj*2;for(i=0;i<nj2;i++) p[i]=0.0;for(i=0;i<npj;i++){fscanf(infile,"%d%f",&jj,&pj);p[jj]=pj;}fprintf(outf,"The Num. Of Nodes: %3d\n",nj);fprintf(outf,"The Num. Of Elements.: %3d\n",ne);fprintf(outf,"The Num. Of Type Of Section Characteristic: %3d\n",nm); fprintf(outf,"The Num. Of Restriction: %3d\n",nz);fprintf(outf,"The Num. Of Nodal Loads: %3d\n",npj);fprintf(outf,"LXM= %d\n",lxm);fprintf(outf,"Coordinates x and y Of Nodes:\n");fprintf(outf," Node x y\n");for(i=1;i<=nj;i++){fprintf(outf,"%10d%10.2f%10.2f\n",i,xy[i][1],xy[i][2]);}fprintf(outf,"The Nodes Num. Of Mem:\n");fprintf(outf," Mem. Type Left Right\n");for(i=1;i<=ne;i++){fprintf(outf,"%10d%10d%10d%10d%10d\n",i,jm[i][0],jm[i][1],jm[i][2],jm[i][3]);}for(i=1;i<=nm;i++){fprintf(outf,"nm=%d E=%f,UM=%f,OUL=%f,T=%f\n",i,etm[i][1],etm[i][2],etm[i][3],etm[i][ 4]);}fprintf(outf,"The Position Of Retriction:");for(i=1;i<=nz;i++){fprintf(outf,"%d ",nzc[i]);}fprintf(outf,"\n===========loads of nodes==============\n");for(i=1;i<=nj2;i++){fprintf(outf,"%12d %12.2f\n",i,p[i]);}if(lxm>0){for(i=1;i<=nm;i++){etm[i][1]=etm[i][1]/(1-etm[i][2]*etm[i][2]); }}fclose(infile);}/*------------计算半带宽函数----------------*/ int bwidth(){int ie,i,j,min,max,jj,ib,lk,m,m1,m2,t,nn[8]; ib=0;for(ie=1;ie<=ne;ie++){m=abs(jm[ie][1]-jm[ie][2]);m1=abs(jm[ie][1]-jm[ie][3]);m2=abs(jm[ie][2]-jm[ie][3]);if(m>ib) ib=m;if(m1>ib) ib=m1;if(m2>ib)ib=m2;}bw=2*ib+2;fprintf(outf," bw=%d\n ",bw);return bw;}/*--------------形成单元刚度矩阵函数-----------*/void stiff(int ie){int i,j,k,m,n,i1,j1,i2,j2,ii,jj;float cc,bb,bc,cb,a1,u1,c[4],b[4];i=jm[ie][1];j=jm[ie][2];m=jm[ie][3]; /*取出单元节点号*/n=jm[ie][0];e=etm[n][1];um=etm[n][2];oul=etm[n][3];te=etm[n][4]; /*取出材料参数*/c[3]=xy[j][1]-xy[i][1];b[3]=xy[i][2]-xy[j][2];c[2]=xy[i][1]-xy[m][1];b[2]=xy[m][2]-xy[i][2];c[1]=-c[2]-c[3];b[1]=-b[2]-b[3]; /*计算[B]阵中的bi,ci,bj,cj,bm,cm*/ ea=(b[2]*c[3]-b[3]*c[2])/2.0;/*计算单元面积*/a1=e*te/4.0/(ea-um*um*ea);u1=(1.0-um)/2.0;for(ii=1;ii<=3;ii++) /*按2子块循环求出单刚计算单元面积*/ {for(jj=1;jj<=3;jj++){bb=b[ii]*b[jj];cc=c[ii]*c[jj];cb=c[ii]*b[jj];bc=b[ii]*c[jj];i2=2*ii;j2=2*jj;i1=i2-1;j1=j2-1;ek[i1][j1]=a1*(bb+u1*cc);ek[i1][j2]=a1*(um*bc+u1*cb);ek[i2][j1]=a1*(um*cb+u1*bc);ek[i2][j2]=a1*(cc+u1*bb);}}}/*-------集成总刚函数，存在刚下三角--------------------------*/int ekzk(int ie){int i1,j1,i,j,i2,j2,ii,jj,ji;for(i1=1;i1<=3;i1++){for(i2=1;i2<3;i2++){i=2*(i1-1)+i2;ii=2*(jm[ie][i1]-1)+i2;for(j1=1;j1<=3;j1++){for(j2=1;j2<3;j2++){j=2*(j1-1)+j2;jj=2*(jm[ie][j1]-1)+j2;ji=bw+jj-ii;if(ji<=bw)K[ii-1][ji-1]=K[ii-1][ji-1]+ek[i][j];}}}}}/*--------用0，1置换法进行约束处理-----------------------------*/ void restrict(){int i,j,i1,j1,ii;for(i=1;i<=nz;i++){j=nzc[i]-1;for(i1=0;i1<bw-1;i1++){K[j][i1]=00.00;ii=j+i1+1;if(ii<nj2)K[ii][bw-i1-2]=00.00;}K[j][bw-1]=1.00;f[j]=0;}}/*----------计算单元应力函数-----------------*/void force(){int i,j,ie,m,n,j1,j2,ii,jj;float aa,pyl,ryl,w[7],c[4],b[4],s1[4][7],yl[4],yl1[4];fprintf(outf,"==========stresses of elements ==========\n"); fprintf(outf,"NO SX SY ST S1 S2 CETA\n\n"); for(ie=1;ie<=ne;ie++){i=jm[ie][1];j=jm[ie][2];m=jm[ie][3];n=jm[ie][0];e=etm[n][1];um=etm[n][2];oul=etm[n][3];te=etm[n][4];c[3]=xy[j][1]-xy[i][1];b[3]=xy[i][2]-xy[j][2];c[2]=xy[i][1]-xy[m][1];b[2]=xy[m][2]-xy[i][2];c[1]=-c[2]-c[3];b[1]=-b[2]-b[3];ea=(b[2]*c[3]-b[3]*c[2])/2.0; aa=e/(2.0*ea*(1-um*um)); w[1]=f[2*i-2];w[2]=f[2*i-1];w[3]=f[2*j-2];w[4]=f[2*j-1];w[5]=f[2*m-2];w[6]=f[2*m-1];for(ii=1;ii<=3;ii++){j1=2*ii-1;j2=j1+1;s1[1][j1]=b[ii]*aa;s1[2][j1]=s1[1][j1]*um;s1[3][j2]=s1[1][j1]*(1-um)/2.0;s1[2][j2]=c[ii]*aa;s1[1][j2]=s1[2][j2]*um;s1[3][j1]=s1[2][j2]*(1-um)/2.0;}for(ii=1;ii<=3;ii++){yl[ii]=0.0;for(jj=1;jj<=6;jj++){yl[ii]=yl[ii]+s1[ii][jj]*w[jj];}}pyl=(yl[1]+yl[2])/2.0;ryl=sqrt(((yl[1]-yl[2])/2.0)*((yl[1]-yl[2])/2.0)+yl[3]*yl[3]); yl1[1]=pyl+ryl;yl1[2]=pyl-ryl;aa=abs(yl[2]-yl1[2]);if(abs(yl[2]-yl1[2])<6.1e-2){yl1[3]=0.0;}else yl1[3]=90.0-57.29578*atan(yl[3]/(yl[2]-yl1[2]));fprintf(outf,"%2d%10.2f%10.2f%10.2f%10.2f%10.2f%10.2f\n",ie,yl[1],yl[2],yl[3],yl1[1],yl1[2 ],yl1[3]);}}/*---------------------------------------------------------*/void main(){int ie,i,j,ii,jj,i1,j1;float p1,p2,p3;input();bw=bwidth();K=TwoArrayAlloc(nj2,bw);f=(double*)calloc(nj2,sizeof(double));if(f==NULL) exit(1);for(i=0;i<nj2;i++) f[i]=p[i+1];for(i=0;i<nj2;i++)for(j=0;j<bw;j++) K[i][j]=0;for(ie=1;ie<=ne;ie++){stiff(ie);ekzk(ie);}restrict();if(!BAND2(nj2,bw-1,K,f)) {printf("Failed !\n"); exit(0);}fprintf(outf,"=============== displacements of nodes ================\n");fprintf(outf,"No x y\n");for(ii=0;ii<nj;ii++){p1=f[2*ii];p2=f[2*ii+1];fprintf(outf,"%2d %14.5f %14.5f\n",ii+1,p1,p2);}force();}平面桁架静力分析有限元程序#define NODE_NUM 300#define ELE_NUM 200#include<stdlib.h>#include<stdio.h>#include<math.h>int nj,ne,nm,nz,npj,jm[ELE_NUM][3],nzc[60],nj2,bw;float xy[NODE_NUM][4],EA[5],p[NODE_NUM],pj,F[ELE_NUM];float ek[5][5];double **K,*f;int BAND2(int n,int m,double **K,double *f);double **TwoArrayAlloc(int,int);void TwoArrayFree(double**);FILE *outfile;/*-----------------------------*/double **TwoArrayAlloc(int r,int c){double *x,**y;int n;x=(double*)calloc(r*c,sizeof(double));y=(double**)calloc(r,sizeof(double*));for(n=0;n<=r-1;++n)y[n]=&x[c*n];return(y);}/*--------------------------------*/void TwoArrayFree(double **x){free(x[0]);free(x);}/*--------------------------------*/int BAND2(int n,int m,double **K,double *f) /*解方程函数*/ {int i,j,t,ij,ji,it,jt,tm,m1;double s,w;m1=m+1;for(i=1;i<=n;i++){if(K[i-1][m1-1]<=0) return(0);w=0.0;if(i>m1) tm=i-m;else tm=1;for(j=tm;j<=i;j++){s=0.0;ij=j-i+m1;for(t=tm;t<=j-1;t++){it=t-i+m1;jt=t-j+m1;s=s+K[i-1][it-1]*K[j-1][jt-1]/K[t-1][m1-1];}K[i-1][ij-1]=K[i-1][ij-1]-s;if(j==i)f[i-1]=f[i-1]-w;elsew=w+K[i-1][ij-1]*f[j-1]/K[j-1][m1-1];}}for(i=n;i>=1;i--){s=0.0;if(i>n-m1) tm=n; else tm=i+m;for(j=i+1;j<=tm;j++){ji=i-j+m1;s=s+K[j-1][ji-1]*f[j-1];}f[i-1]=(f[i-1]-s)/K[i-1][m1-1];}return 1;}/*--------------------------------------------*/void input()/*输入函数*/{int jj,j,i,nh,nl;float dx,dy;FILE *infile;infile=fopen("dat.txt","r");outfile=fopen("out.txt","w");fscanf(infile,"%d%d%d%d%d",&nj,&ne,&nm,&nz,&npj); /*读入控制数据*/fprintf(outfile,"The Num. Of Nodes: %3d\n",nj);/*输出节点数*/fprintf(outfile,"The Num. Of Mem.: %3d\n",ne);/*输出单元数*/fprintf(outfile,"The Num. Of Type Of Section Characteristic: %3d\n",nm);/*输出材料类型数*/fprintf(outfile,"The Num. Of Restriction: %3d\n",nz);/*输出荷载数*/fprintf(outfile,"The Num. Of Nodal Loads: %3d\n",npj);/*输出荷载数*/for(i=1;i<=nj;i++){for(j=1;j<3;j++)fscanf(infile,"%f",&xy[i][j]);/*读入节点坐标*/}fprintf(outfile,"Coordinates x and y Of Nodes:\n");fprintf(outfile," Node x y\n");for(i=1;i<=nj;i++){fprintf(outfile,"%10d%10.2f%10.2f\n",i,xy[i][1],xy[i][2]);/*输出节点坐标*/}for(i=1;i<=ne;i++){for(j=0;j<3;j++)fscanf(infile,"%d",&jm[i][j]);/*读入单元信息（材料类型、单元左右节点码）*/}fprintf(outfile,"The Nodes Num. Of Mem.:\n");fprintf(outfile," Mem. Type Left Right\n");for(i=1;i<=ne;i++){fprintf(outfile,"%10d%10d%10d%10d\n",i,jm[i][0],jm[i][1],jm[i][2]);/*输出单元信息*/}for(i=1;i<=nm;i++) fscanf(infile,"%f",&EA[i]);/*读入材料刚度*/for(i=1;i<=nm;i++){fprintf(outfile,"nm=%d EA=%f\n",i,EA[i]);/*输出材料刚度*/}for(i=1;i<=nz;i++) fscanf(infile,"%d",&nzc[i]);/*读入约束位置*/fprintf(outfile,"The Position Of Restriction: ");for(i=1;i<=nz;i++){fprintf(outfile,"%d",nzc[i]);/*输出约束位置*/}fprintf(outfile,"\n");nj2=nj*2;/*形成荷载*/for(i=0;i<nj2;i++) p[i]=0.0;for(i=0;i<npj;i++){fscanf(infile,"%d%f",&jj,&pj);/*读入荷载位置、荷载大小*/p[jj]=pj;}fclose(infile);}/*-------------------------------------------------------------------------*/int bwidth()/*求半带宽函数*/{int ie,i,j,min,max,jj,ib,lk,m,nn[8];ib=0;for(ie=1;ie<=ne;ie++){m=abs(jm[ie][1]-jm[ie][2]);if(m>ib) ib=m;/*找出单元左右节点码之差的最大值*/}return 2*ib+2;/*返回半带宽*/}/*-------------------------------------------------------------------------------*/ void stiff(int ie) /*形成单刚函数*/{int i,j,k,m,n;float dx,dy,dz,l,cx,cy,cz,ea;i=jm[ie][1];/*单元左节点号*/j=jm[ie][2];/*单元右节点号*/m=jm[ie][0];/*单元材料类型号*/dx=xy[j][1]-xy[i][1];/*单元左右节点横坐标之差*/dy=xy[j][2]-xy[i][2];/*单元左右节点纵坐标之差*/l=sqrt(dx*dx+dy*dy);/*单元长度*/cx=dx/l;/*求余弦*/cy=dy/l;/*求正弦*/ea=EA[m]/l;/*fprintf(outfile,"%d%10.2f\n",ie,l);*/ek[1][1]=ek[3][3]=cx*cx*ea;/*求单刚*/ek[2][2]=ek[4][4]=cy*cy*ea;ek[2][1]=ek[1][2]=ek[3][4]=ek[4][3]=cx*cy*ea;ek[4][1]=ek[1][4]=ek[3][2]=ek[2][3]=-cx*cy*ea;ek[1][3]=ek[3][1]=-cx*cx*ea;ek[2][4]=ek[4][2]=-cy*cy*ea;}/*-----------------------------------------------------------------------------*/ int ekzk(int ie)/*集成总刚度*/{int i1,j1,i,j,i2,j2,ii,jj,ji;for(i1=1;i1<=2;i1++){for(i2=1;i2<=2;i2++){i=2*(i1-1)+i2;ii=2*(jm[ie][i1]-1)+i2;for(j1=1;j1<=2;j1++){for(j2=1;j2<=2;j2++){j=2*(j1-1)+j2;jj=2*(jm[ie][j1]-1)+j2;ji=bw+jj-ii;if(ji<=bw) K[ii-1][ji-1]=K[ii-1][ji-1]+ek[i][j];}}}}}/*---------------------------------------------------------------------------*/ void force()/*求单元轴力函数*/{int i,j,ie,m;float dx,dy,dz,l,cx,cy,cz,ea,w[7];fprintf(outfile,"============dyzl==================\n"); for(ie=1;ie<=ne;ie++){i=jm[ie][1];j=jm[ie][2];m=jm[ie][0];w[1]=f[2*i-2];w[2]=f[2*i-1];w[3]=f[2*j-2];w[4]=f[2*j-1];dx=xy[j][1]-xy[i][1];dy=xy[j][2]-xy[i][2];l=sqrt(dx*dx+dy*dy);cx=dx/l;cy=dy/l;ea=EA[m]/l;dx=w[3]-w[1];/*左右节点水平位移之差*/dy=w[4]-w[2];/*左右节点竖直位移之差*/l=ea*(cx*dx+cy*dy);/*求单元轴力*/F[ie]=1;fprintf(outfile,"%d%10.4f\n",ie,l);/*输出单元号、单元轴力*/}}/*--------------------------------------------------------------------------------------*/void main(){int ie,i,j,ii,jj,i1,j1;float p1,p2,p3;input();bw=bwidth();K=TwoArrayAlloc(nj2,bw);f=(double*)calloc(nj2,sizeof(double));if(f==NULL) exit(1);for(i=0;i<nj2;i++) f[i]=p[i+1];for(ie=1;ie<=ne;ie++){stiff(ie);ekzk(ie);}for(i=1;i<=nz;i++)/*约束处理（后处理的0、1置换法）*/{j=nzc[i]-1;for(i1=0;i1<bw-1;i1++){K[j][i1]=00.00;/*将有约束处的行置0*/ii=j+i1+1;if(ii<nj2) K[ii][bw-i1-2]=00.00;/*将有约束处的列置0*/}K[j][bw-1]=1.00;/*将对角线元素置1*/f[j]=0;}/*fprintf(outfile,"========================zk======================\n");for(i=0;i<nj2;i++){for(j=0;j<bw;j++){fprintf(outfile,"%6.2f",K[i][j]);}fprintf(outfile,"\n");}*/fprintf(outfile,"============jdhz=================\n");for(i=0;i<nj2;i++){fprintf(outfile,"%d %6.2f\n",i+1,f[i]);/*输出节点力*/}if(!BAND2(nj2,bw-1,K,f)){printf("Failed!\n"); exit(0);}fprintf(outfile,"=========jdwy==========\n");fprintf(outfile,"xy\n");for(ii=0;ii<nj;ii++){p1=f[2*ii];p2=f[2*ii+1];fprintf(outfile,"%2d %14.5f %14.5f\n",ii+1,p1,p2);/*输出节点位移*/ }force();fclose(outfile);}平面刚架静力分析有限元程序#include<io.h>#include"string.h"#include"stdlib.h"#include"stdio.h"#include"math.h"void hbw();void sncs(int nel);void fix(int np);void trmat();void fis(int nel);void fpj();void force();void stiff();void addsm();void restr();void matmul();void soleq();void outdis();float sm[6][6],tsm[300][30],res[60][2],pj[300],t[6][6],d[6][6],r[300], fo[6],foj[6],pf[200][4],x[100],y[100],ae[10][3],sl,sn,cs,eal,eil;int nj,ne,nt,nr,npj,npf,nn,mbw,jel[100][2],nae[100],is[6];FILE *infile,*outfile;/************主函数**************/void main(){char name1[30],name2[30];int i,j,nel,np;printf("please enter data-filename\n");scanf("%s",name1);printf("please enter result-filename\n");scanf("%s",name2);if((infile=fopen(name1,"r"))!=NULL){fscanf(infile,"%d%d%d%d%d%d",&nj,&ne,&nt,&nr,&npj,&npf);for(i=0;i<ne;i++) fscanf(infile,"%d%d%d",&jel[i][0],&jel[i][1],&nae[i]); for(i=0;i<nj;i++) fscanf(infile,"%f%f",&x[i],&y[i]);for(i=0;i<nt;i++) fscanf(infile,"%f%f%f",&ae[i][0],&ae[i][1],&ae[i][2]); for(i=0;i<nr;i++) fscanf(infile,"%f%f",&res[i][0],&res[i][1]);}else{printf("the data-file not exit!");exit(1);}nn=3*nj;outfile=fopen(name2,"w");if(outfile==NULL){printf("the result-file not exist!");exit(1);}fprintf(outfile,"statis analysis of plane frame\n");fprintf(outfile,"input data\n");fprintf(outfile,"************\n");fprintf(outfile,"control data\n");fprintf(outfile,"the num.of node:%3d\n",nj);fprintf(outfile,"the num.of mem:%3d\n",ne);fprintf(outfile,"the num.of type of section characteristic:%3d\n",nt); fprintf(outfile,"the num.of restricted degrees of freedom:%3d\n",nr); fprintf(outfile,"the num.of nodal loads:%3d\n",npj);fprintf(outfile,"the num.of non-nodal loads:%3d\n",npf);fprintf(outfile,"the num.of nodal degrees of freedom:%3d\n",nn); fprintf(outfile,"information of mem.\n");fprintf(outfile,"mem. start node end node type\n");for(i=0;i<ne;i++)fprintf(outfile,"%5d%10d%10d%10d\n",i+1,jel[i][0],jel[i][1],nae[i]); fprintf(outfile,"coordinates x and y of nodes\n");fprintf(outfile,"node x y\n");for(i=0;i<nj;i++)/*第二次输入*/fprintf(outfile,"%10d%10.2f%10.2f\n",i+1,x[i],y[i]);fprintf(outfile,"information of cross section each mem.\n");fprintf(outfile,"typeaeramoment of interia elastic modulus\n");for(i=0;i<nt;i++)fprintf(outfile,"%8d%15.5f%15.5f%20.2f\n",i+1,ae[i][0],ae[i][1],ae[i][2]); fprintf(outfile,"information restriction\n");fprintf(outfile,"rester.-norestr.-disp.-norestr.-disp.\n");for(i=0;i<nr;i++)fprintf(outfile,"%10d%19.3f%19.3f\n",i+1,res[i][0],res[i][1]);hbw();for(i=0;i<nn;i++) pj[i]=0.0;if(npj==0) goto aa;fprintf(outfile,"nodal loads\n");fprintf(outfile,"nodepxpyzm\n");for(i=0;i<npj;i++){ int no;float px,py,zm;fscanf(infile,"%d%f%f%f",&no,&px,&py,&zm);fprintf(outfile,"%10d%10.2f%10.2f%10.2f\n",no,px,py,zm);pj[3*no-3]=px; pj[3*no-2]=py; pj[3*no-1]=zm;}aa:for(i=0;i<nn;i++) r[i]=0.0;for(i=0;i<nr;i++){int ni;ni=floor(res[i][0]+0.1)-1;if(pj[ni]!=0.0)r[ni]=r[ni]-pj[ni];}if(npf==0) goto bb;fprintf(outfile,"non-nodal loads\n");fprintf(outfile,"loads leng.supp.to load mem type\n");for(i=0;i<npf;i++){fscanf(infile,"%f%f%f%f%f",&pf[i][0],&pf[i][1],&pf[i][2],&pf[i][3]);fprintf(outfile,"%15.2f%15.2f%15.2f%15.2f\n",pf[i][0],pf[i][1],pf[i][2],pf[i][3]); }for(np=0;np<npf;np++){nel=floor(pf[np][2]+0.1);sncs(nel-1);fix(np);trmat();fis(nel-1);fpj();}bb:for(i=0;i<nn;i++)for(j=0;j<mbw;j++) tsm[i][j]=0.0; for(nel=0;nel<ne;nel++){sncs(nel);trmat();stiff();matmul();fis(nel);addsm();}restr();soleq();outdis();force();}/*===================求最大半带宽========================*/ void hbw(){int nel;mbw=0;for(nel=0;nel<ne;nel++){int ma=abs(jel[nel][0]-jel[nel][1]);if(mbw<ma) mbw=ma;}mbw=3*(mbw+1);fprintf(outfile,"Half_Bandwidth Mbw:%5d\n",mbw);}/*矩阵相乘*/void matmul(){int i,j,k;for(i=0;i<6;i++)for(j=0;j<6;j++) d[i][j]=0.0;for(i=0;i<6;i++)for(j=0;j<6;j++)for(k=0;k<6;k++) d[i][j]=d[i][j]+t[k][i]*sm[k][j];for(i=0;i<6;i++)for(j=0;j<6;j++) sm[i][j]=0.0;for(i=0;i<6;i++)for(j=0;j<6;j++)for(k=0;k<6;k++)sm[i][j]=sm[i][j]+d[i][k]*t[k][j];}/*求单元常数*/void sncs(int nel){int ii,jj,k;float xi,yi,xj,yj;ii=jel[nel][0];jj=jel[nel][1];xi=x[ii-1];xj=x[jj-1];yi=y[ii-1];yj=y[jj-1];sl=sqrt((xj-xi)*(xj-xi)+(yj-yi)*(yj-yi)); sn=(yj-yi)/sl;cs=(xj-xi)/sl;k=nae[nel];eal=ae[k-1][0]*ae[k-1][2]/sl;eil=ae[k-1][1]*ae[k-1][2]/sl;}/*求单元坐标转换矩阵*/void trmat(){int i,j;for(i=0;i<6;i++)for(j=0;j<6;j++)t[i][j]=0.0;t[0][0]=cs;t[0][1]=sn;t[1][0]=-sn;t[1][1]=cs;t[2][2]=1.0;for(i=0;i<3;i++)for(j=0;j<3;j++) t[i+3][j+3]=t[i][j];}/*求单元固端力*/void fix(int np){float w,c,c1,c2,cc,cc1,cc2;int i,im;w=pf[np][0];c=pf[np][1];c1=c/sl; c2=c1*c1; cc=sl-c;cc1=cc/sl; cc2=cc1*cc1;for(i=0;i<6;i++) fo[i]=0.0;im=floor(pf[np][3]+0.1);if(im==1){fo[1]=-w*c*(1.0-c2+c2*c1/2.0);fo[2]=-w*c*c*(6.0-8.0*c1+3*c2)/12.0;fo[4]=-w*c-fo[1];fo[5]=w*c1*c*c*(4.0-3.0*c1)/12.0;}else if(im==2){fo[1]=-w*(1+2*c1)*cc2;fo[2]=-w*c*c2;fo[4]=-w*(1+2*cc1)*c2;fo[5]=w*cc*c2;}else if(im==3){fo[1]=6*w*c1*cc1/sl;fo[2]=w*cc1*(2-3*cc1);fo[4]=-fo[2];fo[5]=w*c1*(2-3*c1);}else if(im==4){fo[1]=0.25*w*c*(2.-3.*c2+1.6*c2*c1); fo[2]=-w*c*c*(2.-3.*c1+1.2*c2)/6.0; fo[4]=-w*c/2.0-fo[2];fo[5]=0.25*w*c*c*c1*(1.-0.8*c1);}else if(im==5){fo[0]=-w*cc1;fo[3]=-w*c1;}else if(im==6){fo[0]=-w*c*(1.-0.5*c1);fo[3]=-0.5*w*c*c1;}}/*形成总结点荷载向量*/void fpj(){ int i,j,k,ii,jj;for(k=0;k<6;k++)foj[k]=0.0;for(i=0;i<6;i++)for(j=0;j<6;j++) foj[i]=foj[i]+t[j][i]*fo[j]; for(ii=0;ii<6;ii++) pj[is[ii]]=pj[is[ii]]-foj[ii]; }/*形成单元刚度矩阵*/void stiff(){ int i,j;float s1,s2;for(i=0;i<6;i++)for(j=0;j<6;j++)sm[i][j]=0.0; sm[0][0]=eal;sm[3][0]=-eal; sm[3][3]=eal;s1=12.*eil/(sl*sl);s2=6.*eil/sl;sm[1][1]=s1; sm[4][1]=-s1;sm[4][4]=s1; sm[2][1]=s2;sm[5][1]=s2; sm[4][2]=-s2;sm[5][4]=-s2; sm[2][2]=4*eil;sm[5][5]=4*eil; sm[5][2]=2*eil;for(i=0;i<6;i++)for(j=0;j<i;j++)sm[j][i]=sm[i][j];}/*由单元位移分量码L形成总刚位移分量码IS（L）*/ void fis(int nel){ int i,j;for(i=0;i<2;i++)for(j=0;j<3;j++) is[3*i+j]=3*(jel[nel][i]-1)+j; }/*形成结构原始总刚度矩阵*/void addsm(){ int i,j;for(i=0;i<6;i++)for(j=0;j<6;j++){ int kc;kc=is[j]-is[i];if(kc>=0) tsm[is[i]][kc]=tsm[is[i]][kc]+sm[i][j];}}/*约束处理*/void restr(){ int i;for(i=0;i<nr;i++){ int ni;ni=floor(res[i][0]+0.1); tsm[ni-1][0]=1.0e25;pj[ni-1]=tsm[ni-1][0]*res[i][1]; }}/*解线性方程组*/void soleq(){ float c1;int k,ni,im,i,m,j,nm,jm;for(k=1;k<nn;k++){ if(fabs(tsm[k-1][0])<=0.000001){ fprintf(outfile,"*****singularity in row %4d*****\n",k+1); goto fin;}ni=k+mbw-1; im=ni;if(ni>nn) im=nn;for(i=k+1;i<im+1;i++){ m=i-k+1;c1=tsm[k-1][m-1]/tsm[k-1][0];for(j=1;j<mbw-m+2;j++)tsm[i-1][j-1]=tsm[i-1][j-1]-c1*tsm[k-1][j+i-k-1];pj[i-1]=pj[i-1]-c1*pj[k-1]; }}if(fabs(tsm[nn-1][0])<=0.000001){ fprintf(outfile,"******singularity in row %4d******\n",nn); goto fin; }pj[nn-1]=pj[nn-1]/tsm[nn-1][0];for(i=nn-1;i>0;i--){ nm=nn-i+1;jm=nm;if(nm>mbw)jm=mbw;for(j=2;j<jm+1;j++)pj[i-1]=pj[i-1]-tsm[i-1][j-1]*pj[j+i-2];pj[i-1]=pj[i-1]/tsm[i-1][0]; }fin:return;}/*输出位移向量*/void outdis(){ int i;fprintf(outfile,"output results\n");fprintf(outfile,"******\n");fprintf(outfile,"nodal displacements\n");fprintf(outfile," node u v ceta\n");for(i=0;i<nj;i++) fprintf(outfile,"%10d%15.6f%15.6f%15.6f\n",i+1,pj[3*i],pj[3*i+1],pj[3*i+2]); }/*求单元杆端力、支座反力或结合点力*/void force (){ float dj[6],f[6],fj[6],dd[6];int nel,np,i,j,ip;fprintf(outfile,"Mem.N1 Q1 M1 N2 Q2 M2 \n");for(nel=0;nel<ne;nel++){sncs(nel);trmat();stiff();fis(nel);for(i=0;i<6;i++) dj[i]=pj[is[i]];for(i=0;i<6;i++) {dd[i]=0.0;f[i]=0.0;}for(i=0;i<6;i++)for(j=0;j<6;j++) dd[i]=dd[i]+t[i][j]*dj[j];for(i=0;i<6;i++)for(j=0;j<6;j++) f[i]=f[i]+sm[i][j]*dd[j];if(npf==0) goto a;for(np=0;np<npf;np++){ip=floor(pf[np][2]+0.1)-1;if(ip==nel){ fix(np);for(j=0;j<6;j++) f[j]=f[j]+fo[j];}}a:fprintf(outfile,"%5d%11.2f%11.2f%11.2f%11.2f%11.2f%11.2f\n",nel+1,f[0],f[1],f[2],f[3],f[ 4],f[5]);for(i=0;i<6;i++)fj[i]=0.0;for(i=0;i<6;i++)fj[i]=0.0;for(i=0;i<6;i++)for(j=0;j<6;j++)fj[i]=fj[i]+t[j][i]*f[j];for(i=0;i<6;i++) r[is[i]]=r[is[i]]+fj[i]; }fprintf(outfile,"Nodal Reaction\n");fprintf(outfile,"Node RX RY RM\n");for(i=0;i<nj;i++)fprintf(outfile,"%10d%15.4f%15.4f%15.4f\n",i+1,r[3*i],r[3*i+1],r[3*i+2]);}。
在VC＋＋编程平台上实现面向对象有限元程序
在VC＋＋编程平台上实现面向对象有限元程序
曹骥; 阮红河; 等
【期刊名称】《《工程设计CAD与智能建筑》》
【年(卷),期】2002(000)001
【摘要】如何在面向对象有限元理论基础上实现一个实用的有限元软件系统是本文的研究重点。
本文首先在对编程语言、类库以及编程平台的选择作了一番讨论后,确定用VC++编程平台实现面向对象有限元程序。
然后探讨了有限元程序的面向对象编程、MFC类库的使用以及VC++编辑器的使用等具体问题。
文末给出了一个面向对象有限元程序实例(OOFEA)。
【总页数】4页(P14-17)
【作者】曹骥; 阮红河; 等
【作者单位】同济大学建筑工程系
【正文语种】中文
【中图分类】TP311.1
【相关文献】
1.VC＋＋面向对象编程初探 [J], 康芊
2.面向对象通用优化程序在VC++平台上的实现 [J], 鲍弋海;曹骥;袁勇
3.地下工程有限元程序面向对象的设计与实现 [J], 李晓军;朱合华
4.面向对象有限元程序设计及其VC++与Matlab混合编程实现 [J], 史贵才;葛修润
5.利用VC＋＋实现有限元程序设计 [J], 王居林;周彦;高立名
因版权原因，仅展示原文概要，查看原文内容请购买。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
有限元编程的c++实现算例1.#include<stdio.h>2.#include<math.h>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.#definepi16.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. { //把节点载荷送入P58. 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. { //求固端反力F066. hz=i;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++)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; //第一列置1121. for(j=2;j<=dd;j++)122. {123. kz[z][j]=0.0; //行置0124. }125. if((z!=1))126. {127. if(z>dd)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.138.139.140. for(k=1;k<=nj3-1;k++)141. {142. if(nj3>k+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(dd>nj3-i+1)168. j0=nj3-i+1;169. elsej0=dd; //求最大列码j0170. for(j=2;j<=j0;j++)171. {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("NJ U V CETA \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("E N Q M \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++) //求由节点位移引起的单元节点力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++) //按非节点载荷数循环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.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]=(g*c*c*d)/(l0*l0);267. }268. else269. {270. f0[1]=-(g*d/l0); //纵向集中力的固端反力271. f0[2]=0.0;272. f0[3]=0.0;273. f0[4]=-g*c/l0;274. f0[5]=0.0;275. f0[6]=0.0;276. }278.}279.280.//************************************************************ 281.//<功能：构成坐标变换矩阵>282.//************************************************************ 283.voidzb(inte)284.{285. doubleceta,co,si;286. inti,j;287. ceta=(gj[e]*pi)/180; //角度变弧度288. co=cos(ceta);289. si=sin(ceta);290. t[1][1]=co; //计算T右上角元素291. t[1][2]=si;292. t[2][1]=-si;293. t[2][2]=co;294. t[3][3]=1.0;295. for(i=1;i<=3;i++)296. {297. for(j=1;j<=3;j++) //计算T的左下角元素298. {299. t[i+3][j+3]=t[i][j];300. }301. }302.}303.304.305.306.//*****************************************************307.//<计算局部坐标下单元刚度矩阵kd[][]>308.//*****************************************************309.voidjdugd(inte)310.{311. doubleA0,l0,j0;313. intj;314.315.316. A0=mj[e]; //面积317. l0=gc[e]; //杆长318. j0=gx[e]; //惯性钜319.320.321. for(i=0;i<=6;i++)322. for(j=0;j<=6;j++) //kd清0323. kd[i][j]=0.0;324.325. kd[1][1]=e0*A0/l0;326. kd[2][2]=12*e0*j0/pow(l0,3);327. kd[3][2]=6*e0*j0/pow(l0,2);328. kd[3][3]=4*e0*j0/l0;329. kd[4][1]=-kd[1][1];330. kd[4][4]=kd[1][1];331. kd[5][2]=-kd[2][2]; //计算kd左下角各元素332. kd[5][3]=-kd[3][2];333. kd[5][5]=kd[2][2];334. kd[6][2]=kd[3][2];335. kd[6][3]=2*e0*j0/l0;336. kd[6][5]=-kd[3][2];337. kd[6][6]=kd[3][3];338.339. for(i=1;i<=6;i++)340. for(j=1;j<=i;j++) //将kd左下角元素按对称原则送到右下角 341. kd[j][i]=kd[i][j];342.}343.344.345.//********************************************************** 346.//<计算整体坐标下单元刚度矩阵ke[][]>348.voiddugd(inte)349.{350. inti,k,j,m;351. jdugd(e); //计算局部单元刚度矩阵kd352. zb(e); //计算坐标变换矩阵T353. for(i=1;i<=6;i++)354. {355. for(j=1;j<=6;j++)356. {357. ke[i][j]=0.0;358. for(k=1;k<=6;k++)359. {360. for(m=1;m<=6;m++)361. {362. ke[i][j]=ke[i][j]+t[k][i]*kd[k][m]*t[m][j]; //计算刚度坐标内单元刚度矩阵ke 363. }364. }365. }366. }367.}368.369.370.//**程序结束**371.372.373.374.</math.h></stdio.h>。