有限元分析中《结构力学》矩阵位移法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];