allprogram.cpp

有限元程序。用c语言开发的三角形有限元程序。可以计算偎依

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
void main()
{
	int i,j,m,p,k;
	int n,s;
	int judge;
    struct em
	{
		int nd1;
		int nd2;
		int nd3;
	}//定义单元
	*em1;
	struct node
	{
		float x;
		float y;
	};//定义节点
	typedef struct node nd;
	nd *node1;
	float E,v,b[3],c[3],B[3][6],D[3][3],area,t;
	double **K,*F,*d;
	FILE *ifp,*ofp;
		/**********************************************************************/
	if((ifp=fopen("data.txt","r"))==NULL)
	{printf("File open fail\n");
	exit(0);
	}
	fscanf(ifp,"%f%f%f%d%d",&E,&v,&t,&n,&s);
	fscanf(ifp,"%d",&judge);
    em1=(struct em *)malloc(sizeof(struct em)*n);
    node1=(nd *)malloc(sizeof(nd)*s);
	K=(double **)malloc(sizeof(double*)*2*s);
	for(i=0;i<2*s;i++)
		K[i]=(double *)malloc(sizeof(double)*2*s);
	F=(double *)malloc(sizeof(double)*2*s);
	d=(double *)malloc(sizeof(double)*2*s);
	for(i=0;i<n;i++)
	fscanf(ifp,"%d%d%d",&em1[i].nd1,&em1[i].nd2,&em1[i].nd3);
	for(i=0;i<s;i++)
	fscanf(ifp,"%f%f",&node1[i].x,&node1[i].y);
	for(i=0;i<2*s;i++)

	fscanf(ifp,"%lf",&F[i]);
	for(i=0;i<2*s;i++)
	fscanf(ifp,"%lf",&d[i]);
	/*------------------------------------------------*/
		ofp=fopen("output.txt","w");	
		fprintf(ofp,"弹性模量E为:%f\t泊松比v为:%f\t厚度为:%f\n",E,v,t);
			fprintf(ofp,"单元数n为：%d\t节点数s为：%d\n",n,s);
			fprintf(ofp,"单元及其节点号\n");
			for(i=0;i<n;i++)
				{fprintf(ofp,"%d,%d,%d",em1[i].nd1,em1[i].nd2,em1[i].nd3);
					fprintf(ofp,"\n");
				}
				fprintf(ofp,"节点及其坐标\n");
				for(i=0;i<s;i++)
					{fprintf(ofp,"%12f,%12f",node1[i].x,node1[i].y);
						fprintf(ofp,"\n");
					}//以上是原始数据的读入和输出
	if(judge<0)//区分问题性质
	{
		E=E/(1-v*v);
		v=v*v/(1-v);
	}//下面给弹性矩阵赋值
		D[0][0]=D[1][1]=E/(1-v*v),D[0][1]=D[1][0]=E*v/(1-v*v),D[0][2]=D[1][2]=D[2][0]=D[2][1]=0,D[2][2]=E/(2*(1+v));
  
	for(i=0;i<3;i++)
      for(j=0;j<6;j++)
       B[i][j]=0;
	  for(i=0;i<2*s;i++)
		  for(j=0;j<2*s;j++)
			  K[i][j]=0;         //总刚赋零             

	for(p=0;p<n;p++)
	{	
		i=em1[p].nd1-1;
        j=em1[p].nd2-1;
        m=em1[p].nd3-1;

        b[0]=node1[j].y-node1[m].y;b[1]=node1[m].y-node1[i].y;b[2]=node1[i].y-node1[j].y;
	    c[0]=-node1[j].x+node1[m].x;c[1]=-node1[m].x+node1[i].x;c[2]=-node1[i].x+node1[j].x;
		area=(b[1]*c[2]-c[1]*b[2]);
	    B[0][0]=b[0]/area;B[0][2]=b[1]/area;B[0][4]=b[2]/area;
        B[1][1]=c[0]/area;B[1][3]=c[1]/area;B[1][5]=c[2]/area;
        B[2][0]=c[0]/area;B[2][1]=b[0]/area;B[2][2]=c[1]/area;
        B[2][3]=b[1]/area;B[2][4]=c[2]/area;B[2][5]=b[2]/area;
    float B1[6][3],C[6][3],ek[6][6];
	int r,q;
			for(r=0;r<3;r++)
				for(q=0;q<6;q++)
					B1[q][r]=B[r][q];
	//矩阵相乘求单刚
for(r=0;r<6;r++)
 for(q=0;q<3;q++)
    {C[r][q]=0;
     for(k=0;k<3;k++)
     C[r][q]+=B1[r][k]*D[k][q];
 }
	 for(r=0;r<6;r++)
		 for(q=0;q<6;q++)
		 {ek[r][q]=0;
		 for(k=0;k<3;k++)
			 ek[r][q]+=(C[r][k]*B[k][q]);
		 }
		 for(r=0;r<6;r++)
			 for(q=0;q<6;q++)
			 {
				 ek[r][q]=ek[r][q]*area*t/2;
			 }
             for(r=-1;r<1;r++)
				 for(q=-1;q<1;q++)
				 {//组合总刚
					 K[2*i+r+1][2*i+q+1]+=ek[1+r][1+q];K[2*i+r+1][2*j+q+1]+=ek[1+r][3+q];K[2*i+r+1][2*m+q+1]+=ek[1+r][5+q];
			         K[2*j+r+1][2*i+q+1]+=ek[3+r][1+q];K[2*j+r+1][2*j+q+1]+=ek[3+r][3+q];K[2*j+r+1][2*m+q+1]+=ek[3+r][5+q];
			         K[2*m+r+1][2*i+q+1]+=ek[5+r][1+q];K[2*m+r+1][2*j+q+1]+=ek[5+r][3+q];K[2*m+r+1][2*m+q+1]+=ek[5+r][5+q];
				 }
	}
	/*****************************************************************************/

	for(i=0;i<2*s;i++)
		for(j=0;j<2*s;j++)
			if(i==j&&d[i]==0)
				K[i][j]=9.e20;//对角元素大数法
  /*---------------------------------------------------------------*/
			//下面是高斯消元法求线性方程组的解
	int h;
	double *y,*temp,te,sum;
  y=(double*)malloc(sizeof(double)*2*s);
  for(i=0;i<2*s-1;i++)
   {
    k=i;
    for(j=i;j<2*s;j++)
     {
      if(fabs(K[k][i])<fabs(K[j][i]))
       k=j;
     }
    temp=K[i];K[i]=K[k];K[k]=temp;
    te=F[i];F[i]=F[k];F[k]=te;
    for(j=i+1;j<2*s;j++)
{ 
F[j]=F[j]-F[i]*(K[j][i]/K[i][i]);
       for(h=2*s-1;h>=i;h--)
        K[j][h]-=K[i][h]*K[j][i]/K[i][i];}
     }
  for(i=2*s-1;i>=0;i--)
   {
    sum=F[i];
    for(j=i+1;j<2*s;j++)
     sum-=y[j]*(K[i][j]);
    y[i]=sum/K[i][i];
   }
  fprintf(ofp,"\n\n");
  for(i=0;i<2*s;i++)
  {
	  fprintf(ofp,"y[%d]=%g\t",i+1,y[i]);
  }
free(K);free(y);free(F);
fclose(ifp);
fclose(ofp);
}