计算结构力学，平行桁架杆件轴力计算源程序

此题算例来源于科学出版社《计算结构力学》，阎军，杨春秋编著中第三章，课后习题第3-2题，仅供参考！！！

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <windows.h>
/*========================
**************************
-----全局变量定义区域-----
**************************
========================*/
int n,m,nn,nc;
int ihl[100],ihr[100];
float clxy[100],t[100],h[100],ea[100];
int reg;
int i0,rd,j01;

float x[100],y[100];
float r[100][100],c[100][100];
float v[100],dp[100];
void i0j0(int k)	//计算对号指示数
{
	int i,j;
	i=ihl[k-1];
	j=ihr[k-1];
	i0=2*(i-nc-1);
	j01=2*(j-nc-1);
}

void ch(int k)	//计算杆长方向余弦
{
	int i,j;
	i=ihl[k-1];
	j=ihr[k-1];
	clxy[1]=x[j-1]-x[i-1];
	clxy[2]=y[j-1]-y[i-1];
	
	clxy[0]=sqrt(clxy[1]*clxy[1]+clxy[2]*clxy[2]);
	clxy[1]=clxy[1]/clxy[0];
	clxy[2]=clxy[2]/clxy[0];
}

void stif(int k)	//计算轴力杆总体系中的单刚
{
	int i,j;
	ch(k);
	t[0]=clxy[1];
	t[1]=clxy[2];
	rd=ea[k-1]/clxy[0];
	for(i=0;i<2;i++)
		for(j=0;j<2;j++)
			c[i][j]=t[i]*t[j]*rd;
}

void dor()	//将平面桁架总体刚度方阵储存于r[nn][nn]中
{	
	int i,j,k;
	for(i=0;i<nn;i++)
		for(j=0;j<nn;j++)
			r[i][j]=0.0;
			for(k=1;k<=m;k++)
			{
				i0j0(k);
				stif(k);
				if(i0>=0)
					for(i=i0+1;i<=i0+2;i++)
					{
						for(j=i0+1;j<=i0+2;j++)
							r[i-1][j-1]+=c[i-i0-1][j-i0-1];
						for(j=j01+1;j<=j01+2;j++)
						{
							r[i-1][j-1]-=c[i-i0-1][j-j01-1];
							r[j-1][i-1]=r[i-1][j-1];
						}
					}
					for(i=j01+1;i<=j01+2;i++)
						
						{
							for(j=j01+1;j<=j01+2;j++)
								r[i-1][j-1]+=c[i-j01-1][j-j01-1];
						}
			}
}

void doh()	//计算杆件轴力
{
int i,k;
	for(k=1;k<=m;k++)
	{
		i0j0(k);
			for(i=0;i<2;i++)
				{
					if(i0<0)
						v[i]=dp[j01+i];
					
					else
						v[i]=dp[j01+i]-dp[i0+i];
					
				}
		stif(k);
		h[k-1]=0.0;
		for(i=1;i<=2;i++)
			h[k-1]=h[k-1]+t[i-1]*v[i-1]*rd;
	}
}

void choldlt()	//改进的平方根法
{
	int i,j,k;
	double s;
	for(i=1;i<nn;i++)
		{
	for(j=0;j<i;j++)
	{
		s=0.0;
		for(k=0;k<j;k++)
		{
			s+=r[i][k]*r[j][k]*r[k][k];
		}
		r[i][j]=(r[i][j]-s)/r[j][j];
	}
	s=0;
	for(k=0;k<i;k++)
	{
		s+=r[i][k]*r[i][k]*r[k][k];
	}
	r[i][i]-=s;
		}
}

void trildlt()	//计算节点位移
{
	int i,j;
	double s;
	for(i=1;i<nn;i++){
		s=0;
		for(j=0;j<i;j++)
		{
			s+=r[i][j]*dp[j];
		}
		dp[i]=dp[i]-s;
	}
	dp[nn-1]/=r[nn-1][nn-1];
	for(i=nn-2;i>=0;i--)
	{
		s=0;
		for(j=i+1;j<nn;j++)
		{
			s+=r[j][i]*dp[j];
		}
		dp[i]=dp[i]/r[i][i]-s;
		
	}
}

int main()	//主函数
{ 
	SetConsoleOutputCP(65001);
	printf("请输入总节点数，杆件数，固定节点数\n");

	scanf("%d%d%d",&n,&m,&nc);
	nn=2*(n-nc);
	printf("请输入节点坐标\n");
	for(reg=0;reg<n;reg++)
	{
		scanf("%f%f",&x[reg],&y[reg]);
	}

	printf("请输入左右节点情况\n");
	for(reg=0;reg<m;reg++)
	{
		scanf("%d%d",&ihl[reg],&ihr[reg]);
	}
	printf("请输入抗拉刚度\n");
	for(reg=0;reg<m;reg++)
	{
		scanf("%f",&ea[reg]);
	}
	printf("正在计算...");
	dor();
	choldlt();
	printf("请输入载荷x,y分量\n");
		for(reg=0;reg<nn;reg++)
	{
		scanf("%f",&dp[reg]);
	}
	trildlt();
	doh();
	for(reg=0;reg<m;reg++)
	{
		printf("杆件%d的内力为：%f\n",reg+1,h[reg]);
	}
printf("计算成功！！！Hello World\n");
	getchar();
	system("pause");
	return 0;
}