结构分析程序设计.docx
《结构分析程序设计.docx》由会员分享，可在线阅读，更多相关《结构分析程序设计.docx（18页珍藏版）》请在冰豆网上搜索。
结构分析程序设计
#include"stdio.h"
#include"math.h"
main（）
{
FILE*fp;
intp,q,a[10],b[10],i,j,k,m,g,n,z,zd[10],u[10],v[10];
floats,t,cx,cy,sum,fx[10],fy[10],d[10],d1[100][100]={0},x[10],y[10],E[10],A[10],K[100][100]={0};
floatk1[10][6][6]={0},T[10][6][6]={0},T1[10][6][6]={0},T2[10][6][6]={0},k2[10][6][6]={0};
floatl[10],fxi[10],fxj[10],fyi[10],fyj[10];
fp=fopen（"in1.txt","r"）;
fscanf（fp,"%d%d%d",&n,&g,&z）;//输入结点数、杆件数及支承数（支承数是指有几个结点有支撑）
for（i=0;i{
fscanf（fp,"%d%d%d",&zd[i],&u[i],&v[i]）;
}//输入支承信息（依次输入水平方向、竖直方向的支承信息，有为‘1’，无为‘0’）
for（i=0;i{
fscanf（fp,"%d%d",&a[i],&b[i]）;
}//输入杆件的两端点
for（i=0;i{
fscanf（fp,"%f%f",&x[i],&y[i]）;
}//输入结点坐标
for（i=0;i{
fscanf（fp,"%f%f",&E[i],&A[i]）;
}//输入杆件信息
for（i=0;i{
fscanf（fp,"%f%f",&fx[i],&fy[i]）;
}//输入各个节点外力情况
fclose（fp）;
for（i=0;i{
j=a[i]-1;
k=b[i]-1;
t=（（x[j]-x[k]）*（x[j]-x[k]））+（（y[j]-y[k]）*（y[j]-y[k]））;
l[i]=sqrt（t）;
}
for（i=0;i{
k1[i][0][0]=E[i]*A[i]/l[i];
k1[i][0][2]=（-1）*E[i]*A[i]/l[i];
k1[i][2][0]=（-1）*E[i]*A[i]/l[i];
k1[i][2][2]=E[i]*A[i]/l[i];
}
/*形成局部坐标系下的单元刚度矩阵*/
/*fopen（"out3.txt","w"）;
fprintf（fp,"局部坐标系下的单元刚度矩阵\n"）;
for（i=0;i{
for（j=0;j<4;j++）
{
for（k=0;k<4;k++）
{
fprintf（fp,"%.2f",k1[i][j][k]）;
}
fprintf（fp,"\n"）;
}
fprintf（fp,"\n\n"）;
}
fclose（fp）;
*/
for（i=0;i{
p=a[i]-1;
q=b[i]-1;
cx=（x[q]-x[p]）/l[i];
cy=（y[q]-y[p]）/l[i];
for（j=0;j<4;j++）
{
T[i][j][j]=cx;
}
T[i][0][1]=cy;
T[i][1][0]=（-1）*cy;
T[i][2][3]=cy;
T[i][3][2]=（-1）*cy;
}
/*形成坐标转换矩阵*/
/*fopen（"out3.txt","w"）;
fprintf（fp,"坐标转换矩阵\n"）;
for（i=0;i{
for（j=0;j<4;j++）
{
for（k=0;k<4;k++）
{
fprintf（fp,"%.2f",T[i][j][k]）;
}
fprintf（fp,"\n"）;
}
fprintf（fp,"\n\n"）;
}
fclose（fp）;
*/
for（i=0;i{
p=a[i]-1;
q=b[i]-1;
cx=（x[q]-x[p]）/l[i];
cy=（y[q]-y[p]）/l[i];
for（j=0;j<4;j++）
{
T1[i][j][j]=cx;
}
T1[i][0][1]=（-1）*cy;
T1[i][1][0]=cy;
T1[i][2][3]=（-1）*cy;
T1[i][3][2]=cy;
}
/*形成坐标转换矩阵的逆矩阵*/
/*fopen（"out3.txt","w"）;
fprintf（fp,"坐标转换矩阵的逆矩阵\n"）;
for（i=0;i{
for（j=0;j<4;j++）
{
for（k=0;k<4;k++）
{
fprintf（fp,"%.2f",T1[i][j][k]）;
}
fprintf（fp,"\n"）;
}
fprintf（fp,"\n\n"）;
}
fclose（fp）;
*/
for（i=0;i{
for（j=0;j<4;j++）
{
for（k=0;k<4;k++）
{
sum=0;
for（m=0;m<4;m++）
{
sum=sum+T1[i][j][m]*k1[i][m][k];
}
T2[i][j][k]=sum;
}
}
}
for（i=0;i{
for（j=0;j<4;j++）
{
for（k=0;k<4;k++）
{
sum=0;
for（m=0;m<4;m++）
{
sum=sum+T2[i][j][m]*T[i][m][k];
}
k2[i][j][k]=sum;
}
}
}
/*形成整体坐标系下的单元刚度矩阵*/
/*fopen（"out3.txt","w"）;
fprintf（fp,"整体坐标系下的单元刚度矩阵\n"）;
for（i=0;i{
for（j=0;j<4;j++）
{
for（k=0;k<4;k++）
{
fprintf（fp,"%.2f",k2[i][j][k]）;
}
fprintf（fp,"\n"）;
}
fprintf（fp,"\n\n"）;
}
fclose（fp）;
*/
for（i=0;i{
p=2*（a[i]-1）;
q=2*（b[i]-1）;
K[p][p]=K[p][p]+k2[i][0][0];
K[p][p+1]=K[p][p+1]+k2[i][0][1];
K[p+1][p]=K[p+1][p]+k2[i][1][0];
K[p+1][p+1]=K[p+1][p+1]+k2[i][1][1];
K[p][q]=K[p][q]+k2[i][0][2];
K[p][q+1]=K[p][q+1]+k2[i][0][3];
K[p+1][q]=K[p+1][q]+k2[i][1][2];
K[p+1][q+1]=K[p+1][q+1]+k2[i][1][3];
K[q][p]=K[q][p]+k2[i][2][0];
K[q][p+1]=K[q][p+1]+k2[i][2][1];
K[q+1][p]=K[q+1][p]+k2[i][3][0];
K[q+1][p+1]=K[q+1][p+1]+k2[i][3][1];
K[q][q]=K[q][q]+k2[i][2][2];
K[q][q+1]=K[q][q+1]+k2[i][2][3];
K[q+1][q]=K[q+1][q]+k2[i][3][2];
K[q+1][q+1]=K[q+1][q+1]+k2[i][3][3];
}
/*形成总刚矩阵*/
for（i=0;i{
K[2*i][2*n]=fx[i];
K[2*i+1][2*n]=fy[i];
}
/*形成增广矩阵*/
/*fopen（"out3.txt","w"）;
fprintf（fp,"形成增广矩阵\n"）;
for（i=0;i<2*n;i++）
{
for（j=0;j<2*n+1;j++）
{
fprintf（fp,"%.3f",K[i][j]）;
}
fprintf（fp,"\n\n"）;
}
fclose（fp）;
*/
for（i=0;i{
if（u[i]==1）
{
for（j=0;j<=2*n;j++）
{
p=2*（zd[i]-1）;
K[p][j]=0;
}
for（j=0;j<2*n;j++）
{
q=2*（zd[i]-1）;
K[j][q]=0;
}
p=2*（zd[i]-1）;
q=2*（zd[i]-1）;
K[p][q]=1;
}
if（v[i]==1）
{
for（j=0;j<=2*n;j++）
{
p=2*zd[i]-1;
K[p][j]=0;
}
for（j=0;j<2*n;j++）
{
q=2*zd[i]-1;
K[j][q]=0;
}
p=2*zd[i]-1;
q=2*zd[i]-1;
K[p][q]=1;
}
}
/*引入支承条件*/
/*fopen（"out3.txt","w"）;
fprintf（fp,"引入支承条件\n"）;
for（i=0;i<2*n;i++）
{
for（j=0;j<2*n+1;j++）
{
fprintf（fp,"%.3f",K[i][j]）;
}
fprintf（fp,"\n\n"）;
}
fclose（fp）;
*/
for（k=0;k<2*n;k++）
{
for（;K[k][k]==0;）
{
for（j=k;j<=2*n;j++）
{
t=K[k][j];
K[k][j]=K[k+1][j];
K[k+1][j]=t;
}
}
/*判断第一个系数是否为0*/
s=K[k][k];
for（i=k,j=k;j<=2*n;j++）
{
K[i][j]=K[i][j]/s;
}
/*将对角线上的书化为1*/
for（i=k+1;i<2*n;i++）
{
s=K[i][k];
for（j=k;j<=2*n;j++）
{
K[i][j]=K[i][j]-s*K[k][j];
}
}
/*将第一行下面的的每一行的第一个数化为0*/
}
/*形成上三角矩阵*/
for（i=2*n-1;i>=0;i--）
{
sum=0;
for（m=2*n-1;m>i;m--）
{
sum=sum+K[i][m]*d[m];
}
/*每个系数与对应的X相乘，得到除X[l]项以外的其他项和*/
d[i]=K[i][2*n]-sum;
/*求d[l]*/
}
/*用高斯消元法求d[i]*/
/*fopen（"out3.txt","w"）;
fprintf（fp,"求解之后的增广矩阵\n"）;
for（i=0;i<2*n;i++）
{
for（j=0;j<2*n+1;j++）
{
fprintf（fp,"%.2f",K[i][j]）;
}
fprintf（fp,"\n\n"）;
}
fclose（fp）;
*/
/*fopen（"out2.txt","w"）;
fprintf（fp,"整体坐标系下的d[i]\n"）;
for（i=0;i<2*n;i++）
{
fprintf（fp,"d[%d]=%.4f\n",i+1,d[i]）;
}
fclose（fp）;
*/
for（i=0;i{
for（j=0;j<4;j++）
{
k=0;
sum=0;
p=2*（a[i]-1）;
q=2*（b[i]-1）;
sum=sum+T[i][j][k]*d[p];
sum=sum+T[i][j][k+1]*d[p+1];
sum=sum+T[i][j][k+2]*d[q];
sum=sum+T[i][j][k+3]*d[q+1];
d1[i][j]=sum;
}
}
/*求局部坐标系下的d[i]*/
/*fopen（"out2.txt","w"）;
fprintf（fp,"局部坐标系下的d[i]\n"）;
for（i=0;i{
for（j=0;j<4;j++）
{
fprintf（fp,"d[%d][%d]=%.4f",i+1,j+1,d1[i][j]）;
}
fprintf（fp,"\n"）;
}
fclose（fp）;
*/
for（i=0;i{
sum=0;
for（k=0;k<4;k++）
{
sum=sum+k1[i][0][k]*d1[i][k];
}
fxi[i]=sum;
sum=0;
for（k=0;k<4;k++）
{
sum=sum+k1[i][1][k]*d1[i][k];
}
fxj[i]=sum;
sum=0;
for（k=0;k<4;k++）
{
sum=sum+k1[i][2][k]*d1[i][k];
}
fyi[i]=sum;
sum=0;
for（k=0;k<4;k++）
{
sum=sum+k1[i][3][k]*d1[i][k];
}
fyj[i]=sum;
}
/*求各个杆件的轴力*/
fopen（"out1.txt","w"）;
fprintf（fp,"各个杆件的轴力\n"）;
for（i=0;i{
fprintf（fp,"fxi[%d]=%.7f\n",i+1,fxi[i]）;
fprintf（fp,"fxj[%d]=%.7f\n",i+1,fxj[i]）;
fprintf（fp,"fyi[%d]=%.7f\n",i+1,fyi[i]）;
fprintf（fp,"fyj[%d]=%.7f\n",i+1,fyj[i]）;
}
fclose（fp）;
/*输出结果*/
}
//实例
/*输入：
442
111201
13232434
00400444
112112121112
00005-1000
结果：
坐标转换矩阵
0.001.000.000.00
-1.000.000.000.00
0.000.000.001.00
0.000.00-1.000.00
-0.710.710.000.00
-0.71-0.710.000.00
0.000.00-0.710.71
0.000.00-0.71-0.71
0.001.000.000.00
-1.000.000.000.00
0.000.000.001.00
0.000.00-1.000.00
1.000.000.000.00
0.001.000.000.00
0.000.001.000.00
0.000.000.001.00
整体坐标系下的单元刚度矩阵
0.000.000.000.00
0.003.000.00-3.00
0.000.000.000.00
0.00-3.000.003.00
1.06-1.06-1.061.06
-1.061.061.06-1.06
-1.061.061.06-1.06
1.06-1.06-1.061.06
0.000.000.000.00
0.003.000.00-3.00
0.000.000.000.00
0.00-3.000.003.00
3.000.00-3.000.00
0.000.000.000.00
-3.000.003.000.00
0.000.000.000.00
将对应位置的/**/去掉之后即可得到相应结果
将内容复制到vc++6.0下之后，在程序的同一个文件夹下建立in1.txt输入要求输入内容，再建对应空白的out1.txt文件，运行之后可在out1.txt看见结果*/
以上内容是自己所编，有不足之处，还请各位多多包涵！
升级会员 