计算结构力学程序设计源代码（指针流）

#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>
#include <windows.h>

/* 全局变量声明 */
int TNN;	//节点总数		
int NFIN;	//固定节点数	
int NFRN;	//可动节点数
int NOR;	//杆件数
int NOL;	//载荷总数
int NOS;	//截面数

float* XCN;		//节点x方向的坐标
float* YCN;		//节点y方向的坐标
int* BNR;		//杆件始端节点号
int* ENR;		//杆件末端节点号
float* EEE;		//弹性模量E
float* AAA;		//横截面积A
float* JJJ;		//主惯性矩J

int* NRL;		//加载杆件号
int* KOL;		//载荷类型
float* VOL;		//载荷大小
float* DLB;		//载荷距始端的距离

int* NRS;		//藏面所在的杆件编号
float* DSB;		//截面距始端的距离

float** LCS;	//杆件的长度、倾角余弦和倾角正弦
float* DON;		//存放节点位移
float* IFS;		//存放截面内力、

/*函数声明*/
void PGRead();//从文本文件中读取数据
float** PGBuildTotalStif();//组集总刚度阵
void PGLCosSin(int k);//求杆件的长度，倾角余弦和正弦
int* PGI0J0(int rod);//节点对号函数
void PGBuildUnitStif(int k, int flag, float us[3][3]);//求单根杆件的单元刚度阵
void PGBuildRD(int k, int flag, float rd[3][3]);//求单根杆件局部系下刚度阵
float* PGBuildLoadVector();//组集载荷向量
void PGReactionForce(int i, float rfb[3], float rde[3]);//计算支座反力
void PGCholesky(float** A, float* b, float* x, int n);//改进的平方根法求解方程组
void PGInternalForce(int ii, int k, float xq);//计算截面内力
void PGCtlInternalForce(int jm, int i, float rfe[3], float tinf[3]);//计算简支梁截面处内力
void PGDispRodEndForce(int k, float tref[3]);//计算由节点位移产生的截面内力
void PGaaa();//格式输出函数
void PGPrint();//结果输出函数


int main()
{
	SetConsoleOutputCP(65001);
	char value = 0;	//存放用户输入的字符
	int i;			//循环控制变量
	float** kk, * pp;//kk为指向总刚度阵的指针，pp为指向载荷向量的指针
	printf("欢迎使用平面刚架求解器 \n");
	printf("\n请在根目录中的参数输入.txt中输入结构数据，保存后返回该界面，按任意键继续操作。\n");
	printf("需注意:\n1.节点编号时先编固定节点再编可动节点，杆件的小号端为始端，大号端为末端。\n");
	printf("2.弹性模量单位为千牛/平方米， 荷载单位为千牛。\n");
	printf("3.荷载类型7的“荷载距始端”距离对应“线膨胀系数”，类型8对应“线膨胀系数/杆件厚度”。\n");
	printf("\n是否开始进行结构节点位移和指定截面内力进行计算?(Y/N): ");
	scanf_s("%c", &value);//用户根据提示信息选择输入字符
	if (value == 'Y' || value == 'y')//用户选择'y'或'Y'则进行以下计算
	{
		PGRead();//从数据文件读入数据
		kk = PGBuildTotalStif();//组集总刚度阵并将其指针赋给kk
		printf("f",kk[0][0]);
		pp = PGBuildLoadVector();//组集载荷向量并将其指针赋给pp
		PGCholesky(kk, pp, DON, 3 * NFRN);//改进的平方根法求解，结果存放在DON中
		for (i = 0; i < NOS; i++)//对每一截面逐个求内力
			PGInternalForce(3 * i, NRS[i], DSB[i]);//求局部系截面内力，结果存放在IFS中

		printf("\n");
		PGaaa();//输出间隔符号
		printf("\n\n");
		PGPrint();//结构参数以及计算结果的输出
	}
	printf("\n 感谢您的使用，再见! \n\n\n");
	getchar();
	system("pause");
	return 0;
}

void PGRead()
{
	FILE* fp;//定义文件指针
	char c;//存放临时的字符型数据
	int i, j;//循环控制变量
	fopen_s(&fp,"F:\\FastWorkLand\\C\\structural mechanics\\ddd.txt", "r");//为读取数据打开文本文件
	fseek(fp, 15L, 0);//将fp所指位置从初始位置向后移动15个字节
	fscanf_s(fp, "%d", &TNN,1);//读取fp指向的整形数据，存放在TNN中
	fseek(fp, 18L, 1);//将fp所指位置从当前位置向后移动17个字节
	fscanf_s(fp, "%d", &NFIN,1);//读取fp指向的整形数据，存放在NFIN中
	fseek(fp, 15L, 1);//将fp所指位置从当前位置向后移动17个字节
	fscanf_s(fp, "%d", &NOR,1);//读取fp指向的整形数据，存放在NOR中
	fseek(fp, 15L, 1);//将fp所指位置从当前位置向后移动17个字节
	fscanf_s(fp, "%d", &NOL,1);//读取fp指向的整形数据，存放在NOL中
	fseek(fp, 15L, 1);//将fp所指位置从当前位置向后移动17个字节
	fscanf_s(fp, "%d", &NOS,1);//读取fp指向的整形数据，存放在NOS中
	fseek(fp, 2L, 1);//将fp所指位置从当前位置向后移动2个字节
	NFRN = TNN - NFIN;//计算可动节点数
	XCN = (float*)calloc(TNN, sizeof(float));//为XCN分配TNN个长度等于float变量的内存空间，下同
	YCN = (float*)calloc(TNN, sizeof(float));
	BNR = (int*)calloc(NOR, sizeof(int));
	ENR = (int*)calloc(NOR, sizeof(int));
	EEE = (float*)calloc(NOR, sizeof(float));
	AAA = (float*)calloc(NOR, sizeof(float));
	JJJ = (float*)calloc(NOR, sizeof(float));
	NRL = (int*)calloc(NOL, sizeof(int));
	KOL = (int*)calloc(NOL, sizeof(int));
	VOL = (float*)calloc(NOL, sizeof(float));
	DLB = (float*)calloc(NOL, sizeof(float));
	NRS = (int*)calloc(NOS, sizeof(int));
	DSB = (float*)calloc(NOS, sizeof(float));
	IFS = (float*)calloc(3 * NOS, sizeof(float));
	DON = (float*)calloc(3 * NFRN, sizeof(float));
	for (i = 0; i < 13; i++)//分别读取13组数据存放在13个数组变量中
	{
		fseek(fp, 24L, 1);//将fp所指位置从当前位置向后移动15个字节
		j = 0;//数组指标，从每组数据的第一个数据开始
		do
		{
			switch (i)//用switch语句控制对各组数据的读取
			{
			case 0:fscanf_s(fp, "%f", &XCN[j],1); break;//i=0时读取fp指向的浮点型数据存放在指标为j的XCN数组中，跳出switch语句，下同
			case 1:fscanf_s(fp, "%f", &YCN[j], 1); break;
			case 2:fscanf_s(fp, "%d", &BNR[j], 1); break;
			case 3:fscanf_s(fp, "%d", &ENR[j], 1); break;
			case 4:fscanf_s(fp, "%f", &EEE[j], 1); break;
			case 5:fscanf_s(fp, "%f", &AAA[j], 1); break;
			case 6:fscanf_s(fp, "%f", &JJJ[j], 1); break;
			case 7:fscanf_s(fp, "%d", &NRL[j], 1); break;
			case 8:fscanf_s(fp, "%d", &KOL[j], 1); break;
			case 9:fscanf_s(fp, "%f", &VOL[j], 1); break;
			case 10:fscanf_s(fp, "%f", &DLB[j], 1); break;
			case 11:fscanf_s(fp, "%d", &NRS[j], 1); break;
			case 12:fscanf_s(fp, "%f", &DSB[j], 1); break;
			}
			fscanf_s(fp, "%c", &c, 1);//读取每个数据后的逗号或换行符
			j++;//数组指标自加
		} while (c != '\n');//若读取的数据后面不是换行符则继续读取
	}
}

float** PGBuildTotalStif()
{
	float** kk, us[3][3];//kk为总刚度阵，us为单元刚度阵分块]
	int i, j, k, m = 0, n = 0, t[2], * p;//i,j,m,n为循环控制变量，k为杆件程序编号，t为普通整形数组，p为存放杆端节点对号位置的数组指针
	kk = (float**)calloc(3 * NFRN, sizeof(float*));//以下三行语句为kk申请二维内存空间
	for (i = 0; i < 3 * NFRN; i++)
		*(kk + i) = (float*)calloc(3 * NFRN, sizeof(float));
	
	for (i = 0; i < 3 * NFRN; i++)					 //对kk总刚清零
		for (j = 0; j < 3 * NFRN; j++)
			kk[i][j] = 0;
	LCS = (float**)calloc(3, sizeof(float*));		//对LCS申请二维内存空间
	for (i = 0; i < 3; i++)
		*(LCS + i) = (float*)calloc(NOR, sizeof(float));
			
	for (k = 0; k < NOR; k++)					    //k为杆件号，对干循环，组装总刚度阵
	{
		PGLCosSin(k + 1);	//计算该杆几何参数
		
		p = PGI0J0(k + 1);	//计算k杆件端点对号位置，存入p数组中
		for (i = 0; i < 2; i++)
			t[i] = p[i];	//防止指针p丢失，赋值给t
		for (i = 0; i < 2; i++)	//叠加主对角分块位置
		{
			if (t[i] >= 0)	//可动节点叠加，否则不叠加
			{
				PGBuildUnitStif(k, 1 + i, us);	
				//计算程序杆件编号为k的杆件的单刚分块，存入us
				for (m = 0; m < 3; m++)
					for (n = 0; n < 3; n++)
						kk[t[i] + m][t[i] + n] += us[m][n];	//对us中的9个元素按相应位置进行叠加
			}
		}
		if (t[0] >= 0 && t[1] >= 0)//符合条件说明两端点均为可动节点并进行叠加，否则不叠加
		{
			for (i = 0; i < 2; i++)//对杆件两端点对应的非主对角分块位置进行叠加
			{
				PGBuildUnitStif(k, 3 + i, us);	//计算程序杆件编号为k的杆件的单元刚度阵分块，存入us
				for (m = 0; m < 3; m++)
					for (n = 0; n < 3; n++)
						kk[t[i] + m][t[1 - i] + n] += us[m][n];	//对us中的9个元素按相应位置进行叠加
			}
						
		}
	}
	return kk;	//返回总刚的数组指针
}
void PGLCosSin(int k)	//k为杆号
{
	int i, j;
	k--;		//kk自减，为程序杆号
	i = BNR[k] - 1;	//i存放杆始节点对应程序中的数组指标
	j = ENR[k] - 1;	//j存放杆终节点对应程序中的数组指标
	LCS[1][k] = XCN[j] - XCN[i];	//杆件始末端节点横坐标差
	LCS[2][k] = YCN[j] - YCN[i];	//杆件始末端节点纵坐标差
	LCS[0][k] = sqrt(LCS[1][k] * LCS[1][k] + LCS[2][k] * LCS[2][k]);  //杆件长度
	LCS[1][k] = LCS[1][k] / LCS[0][k];    //杆件余弦值
	LCS[2][k] = LCS[2][k] / LCS[0][k];	  //杆件正弦值
	printf("hello world!");
}
int* PGI0J0(int rod)	//rod为杆件编号
{
	int bl, br, ij[2];
	bl = BNR[rod - 1];	//bl存放杆始节点号
	br = ENR[rod - 1];	//br存放杆末节点号
	ij[0] = 3 * (bl - NFIN - 1);	//将始端节点在总刚的对应位置存入ij中
	ij[1] = 3 * (br - NFIN - 1);	//将末端节点在总刚的对应位置存入ij中
	printf("hello world!4");
	return ij;	//返回ij指针
}
void PGBuildUnitStif(int k, int flag, float us[3][3])	//k为杆件号，flag为分块号，us为单刚分块
{
	int i, j, m = 0;	//循环量
	float rd[3][3] = { 0 }, t[3][3] = { 0 }, tt[3][3] = { 0 }, c[3][3] = { 0 };	//rd为局部系刚度阵，t为坐标转换阵，tt为其转制
	PGBuildRD(k, flag, rd);	//根据flag计算杆件（k+1）得rd阵相应分块阵
	
	t[0][0] = tt[0][0] = t[1][1] = tt[1][1] = LCS[1][k];	//计算坐标转换阵
	t[1][0] = tt[0][1] = LCS[2][k];			//	cos a	-sin a   0	
	t[0][1] = tt[1][0] = -1 * LCS[2][k];	//	sin a	cos a    0
	t[2][2] = tt[2][2] = 1;					//	  0	      0      1	
	for (i = 0; i < 3; i++)		//对us清零
		for (j = 0; j < 3; j++)
			us[i][j] = 0;
	for (i = 0; i < 3; i++)    //计算t*rd，存入c
		for (j = 0; j < 3; j++)
			for (m = 0; m < 3; m++)
				c[i][j] += t[i][m] * rd[m][j];
	for (i = 0; i < 3; i++)		//计算rd*tt，存入us
		for (j = 0; j < 3; j++)
			for (m = 0; m < 3; m++)
				us[i][j] += c[i][m] * tt[m][j];
			printf("hello world!5");
}
void PGBuildRD(int k, int flag, float rd[3][3])//k为杆件程序编号，flag为矩阵分块标号
{
	float a, b, c, d, e;
	a = EEE[k] * AAA[k] / LCS[0][k];	//EA/L
	d = 4 * EEE[k] * JJJ[k] / LCS[0][k];	//4EJ/L
	c = d / 2 * 3 / LCS[0][k];			//6EJ/(L*L)
	b = c * 2 / LCS[0][k];				//12EJ/(L*L*L)
	e = d / 2;								//2EJ/L
	switch (flag)//根据flag计算相应分块阵
	{
	case 1:rd[0][0] = a; rd[1][1] = b; rd[1][2] = rd[2][1] = -c; rd[2][2] = d; break;  //k11
	case 2:rd[0][0] = a; rd[1][1] = b; rd[1][2] = rd[2][1] = c; rd[2][2] = d; break;  //k22
	case 3:rd[0][0] = -a; rd[1][1] = -b; rd[1][2] = -c; rd[2][1] = c; rd[2][2] = e; break;  //k12
	case 4:rd[0][0] = -a; rd[1][1] = -b; rd[1][2] = c; rd[2][1] = -c; rd[2][2] = e; break;  //k21
	printf("hello world!6");
	}
}
float* PGBuildLoadVector()
{
	int i, j, u, k, rod, * p;		//ij为循环量，rod为杆号
	float rf[2][3], * pp;	//rf存杆件始末端支反力分量，pp为载荷向量
	pp = (float*)calloc(3 * NFRN, sizeof(float));		//为pp分配3*NFRN个长度等于float变量的内存空间
	for (i = 0; i < 3 * NFRN; i++)	//载荷向量清零
		pp[i] = 0;
	for (u = 0; u < NOR; u++)
	{
		KOL[NOL + u] = 1;
		VOL[NOL + u] = 0;
		NRL[NOL + u] = u + 1;
		DLB[NOL + u] = LCS[0][u] / 2;

	}
	for (i = 0; i < NOL + NOR; i++)		//对每一载荷循环
	{
		rod = NRL[i] - 1;	//将NRL中存放的杆件实际编号-1即为程序编号
		for (j = 0; j < 3; j++)		//对rf清零
			rf[0][j] = rf[1][j] = 0;
		PGReactionForce(i, rf[0], rf[1]);    //计算固支梁第i个载荷始末端支反力，存入rf[0]与rf[1]中
		p = PGI0J0(NRL[i]);					//对号位置存入p
		for (j = 0; j < 2; j++)					//对始末节点的相反数进行叠加
			if (p[j] >= 0)					//符合条件为自由节点，叠加，反正固定不叠加
			{
				pp[p[j]] -= rf[j][0] * LCS[1][rod] - rf[j][1] * LCS[2][rod];
				pp[p[j] + 1] -= rf[j][0] * LCS[2][rod] + rf[j][1] * LCS[1][rod];
				pp[p[j] + 2] -= rf[j][2];

			}
			printf("hello world!7");
	}


	/*	for (j = 0; j < NOR;j++)//考虑杆自重计算,对每一根杆循环
		{

			for (j = 0; j < 3; j++)		//对rf清零
				rf[0][j] = rf[1][j] = 0;
			PGReactionForce(j, rf[0], rf[1]);//计算固支梁每根杆自重产生始末端点支反力，分别存放在rf[0],rf[1]中
			p = PGI0J0(j);
				for (k = 0; k < 2; k++)					//对始末节点的相反数进行叠加
					if (p[k] >= 0)					//符合条件为自由节点，叠加，反正固定不叠加
					{
						pp[p[k]] -= rf[k][0] * LCS[1][j] - rf[k][1] * LCS[2][j];
						pp[p[k] + 1] -= rf[k][0] * LCS[2][j] + rf[k][1] * LCS[1][j];
						pp[p[k] + 2] -= rf[k][2];

					}


		}*/
	return pp;
}
void PGReactionForce(int i, float rfb[3], float rfe[3])//i为载荷编号，rfb、rfe分别为始末端支反力分量
{
	float ra, rb, t;
	int rod = NRL[i] - 1;//rod为杆件程序编号
	ra = DLB[i] / LCS[0][rod];//ra=x(q)/L
	rb = 1 - ra;//rb=1-x(q)/L
	switch (KOL[i])//根据KOL[i]计算相应载荷类型的支反力
	{
	case 1:t = rb * rb; rfb[1] = -VOL[i] * (1 + 2 * ra) * t; rfe[1] = -VOL[i] - rfb[1]; rfb[2] = VOL[i] * DLB[i] * t; rfe[2] = -VOL[i] * DLB[i] * rb * ra; break;//竖直集中力
	case 2:t = VOL[i] * DLB[i]; rfb[1] = -t * (1 - ra * ra + 0.5 * ra * ra * ra); rfe[1] = -t - rfb[1]; rfb[2] = t * DLB[i] * (6 - 8 * ra + 3 * ra * ra) / 12; rfe[2] = -t * ra * DLB[i] * (4 - 3 * ra) / 12; break;//竖直均布载荷
	case 3:rfb[0] = -VOL[i] * rb; rfe[0] = -VOL[i] * ra; break;//轴向集中力
	case 4:t = VOL[i] * DLB[i]; rfe[0] = -t * ra / 2; rfb[0] = -t - rfe[0]; break;//轴向均布载荷
	case 5:t = VOL[i] * DLB[i]; rfb[1] = -0.25 * t * (2 - 3 * ra * ra + 1.6 * ra * ra * ra); rfe[1] = -t / 2 - rfb[1]; rfb[2] = t * DLB[i] * (2 - 3 * ra + 1.2 * ra * ra) / 6; rfe[2] = -t * ra * DLB[i] * (1 - 0.8 * ra) / 4; break;//竖直三角形分布载荷
	case 6:rfb[1] = 6 * VOL[i] * rb * ra / LCS[0][rod]; rfe[1] = -rfb[1]; rfb[2] = -VOL[i] * rb * (2 - 3 * rb); rfe[2] = -VOL[i] * ra * (2 - 3 * ra); break;//集中弯矩
	case 7:rfb[0] = VOL[i] * DLB[i] * EEE[rod] * AAA[rod]; rfe[0] = -rfb[0]; break;//均匀温升
	case 8:rfb[2] = VOL[i] * DLB[i] * EEE[rod] * JJJ[rod] * 2; rfe[2] = -rfb[2]; break;//上升温下降温
	}
	printf("hello world!8");
}
void PGCholesky(float** A, float* b, float* x, int n)//A为对称系数阵，b为常数向量，x为未知数向量，n为维数
{
	
	int i, j, k;//循环控制变量
	float s, ** L, * D, * y;//s为中间变量，L、D为分解矩阵，y为中间向量
	L = (float**)calloc(n, sizeof(float*));//以下三行语句为kk申请二维内存空间
	for (i = 0; i < n; i++)
		*(L + i) = (float*)calloc(n, sizeof(float));
	D = (float*)calloc(n, sizeof(float));//为D申请n个长度等于float变量的内存空间
	y = (float*)calloc(n, sizeof(float));//为y申请n个长度等于float变量的内存空间
	for (i = 0; i < n; i++)
		L[i][i] = 1;//L初始化
	/*将A分解为LDL(t)*/
	printf("%f",A[0][0]);
	D[0] = A[0][0];
	for (i = 1; i < n; i++)
	{
		for (j = 0; j < i; j++)
		{
			s = 0;
			for (k = 0; k < j; k++)
			{
				s = s + L[i][k] * D[k] * L[j][k];
			}
			L[i][j] = (A[i][j] - s) / D[j];

		}
		s = 0;
		for (k = 0; k < i; k++)
			s = s + L[i][k] * L[i][k] * D[k];
		D[i] = A[i][i] - s;
	}
	/*由Ly=b求解y*/
	y[0] = b[0];
	for (i = 1; i < n; i++)
	{
		s = 0;
		for (k = 0; k < i; k++)
			s = s + L[i][k] * y[k];
		y[i] = b[i] - s;
	}
	/*由DL(t)x=y求解x*/
	x[n - 1] = y[n - 1] / D[n - 1];
	for (i = n - 2; i >= 0; i--)
	{
		s = 0;
		for (k = i + 1; k < n; k++)
			s = s + L[k][i] * x[k];
		x[i] = y[i] / D[i] - s;
	}
}
void PGInternalForce(int ii, int k, float xp)//ii为截面对号位置，k为杆件实际编号，xp为截面与始端距离
{
	int i, j;//i,j为循环控制变量
	float rfb[3], rfe[3], tf[3];//rfb，rfe分别为始末端支反力分量
	for (i = 0; i < NOL; i++)//对每一个载荷进行循环
	{
		if (NRL[i] == k)//若该载荷施加于杆件k
		{
			for (j = 0; j < 3; j++)//对rfb，rfe，tinf清零
				rfb[j] = rfe[j] = tf[j] = 0;
			PGReactionForce(i, rfb, rfe);//计算固支梁第i个载荷始末端点支反力，分别存放在rfb、rfe中
			IFS[ii] += rfe[0];//以下三行语句根据末端支反力求出截面内力并叠加到对应位置
			IFS[ii + 1] -= rfe[1];
			IFS[ii + 2] = IFS[ii + 2] + rfe[1] * (LCS[0][k - 1] - xp) - rfe[2];//
			PGCtlInternalForce(ii / 3, i, rfe, tf);//对应悬臂梁在载荷下的截面内力计算
			for (j = 0; j < 3; j++)//将悬臂梁的截面内力分量按相应位置进行叠加
				IFS[ii + j] += tf[j];
		}
	}
	PGDispRodEndForce(k, tf);//计算由节点位移产生的杆端力，并存放在数组tf中
	IFS[ii + 0] -= tf[0];//以下三行语句根据公式计算杆端力引起的截面内力并按相应位置进行叠加
	IFS[ii + 1] += tf[1];
	IFS[ii + 2] += tf[2] + tf[1] * xp;
	printf("hello world!10");
}
void PGCtlInternalForce(int jm, int i, float rfe[3], float tinf[3])//jm为截面程序的编号，i为载荷程序的编号，rfe为末端支反力，tinf为截面内力分量
{
	float t = DLB[i] - DSB[jm], r = DSB[jm] / DLB[i];//t=x(q)-x(p),r=x(p)/x(q)
	switch (KOL[i])//根据KOL[i]计算相应载荷类型对应悬臂梁的截面内力
	{
	case 1:if (DSB[jm] < DLB[i]) { tinf[1] -= VOL[i]; tinf[2] += VOL[i] * t; }break;//若符合条件则按照公式进行计算，否则分量为0跳过计算过程，下同
	case 2:if (DSB[jm] < DLB[i]) { tinf[1] -= VOL[i] * t; tinf[2] += 0.5 * VOL[i] * t * t; }break;
	case 3:if (DSB[jm] < DLB[i])  tinf[0] += VOL[i]; break;
	case 4:if (DSB[jm] < DLB[i])  tinf[0] += VOL[i] * t; break;
	case 5:if (DSB[jm] < DLB[i]) { tinf[1] -= VOL[i] * (1 + r) * t / 2; tinf[2] += VOL[i] * t * t * (2 + r) / 6; }break;
	case 6:if (DSB[jm] < DLB[i])  tinf[2] += VOL[i]; break;
	case 7:break;//温度载荷不再悬臂梁上产生内力，下同
	case 8:break;


	}
	printf("hello world!11");
}
void PGDispRodEndForce(int k, float tref[3])//k为杆件实际编号tref存放杆端力分量
{
	int* p, i, j, m, n, t[2];//i,j,m,n为循环控制变量，t为普通整型数组
	float rd[3][3] = { 0 }, tt[3][3] = { 0 }, rdb[3][3], temp;//tt为坐标转换矩阵的转置阵，rdb，temp为临时储存变量
	for (i = 0; i < 3; i++)//对tinf清零
		tref[i] = 0;
	tt[0][0] = tt[1][1] = LCS[1][k - 1];//以下四行语句构造坐标转换矩阵的转置阵
	tt[0][1] = LCS[2][k - 1];
	tt[1][0] = -1 * LCS[2][k - 1];
	tt[2][2] = 1;
	p = PGI0J0(k);//计算实际编号为k的杆件端点对号位置，并存放在p指向的数组中
	for (i = 0; i < 2; i++)
		t[i] = p[i];//将p指向的数值赋给t，避免指针p丢失
	for (i = 0; i < 2; i++)//分别对始末端进行计算叠加
	{
		if (t[i] >= 0)//符合条件为自由节点，执行计算，否则为固定节点
		{
			for (j = 0; j < 3; j++)//对以下三行语句rdb清零
				for (m = 0; m < 3; m++)
					rdb[j][m] = 0;
			PGBuildRD(k - 1, 2 * i + 1, rd);//根据标号计算杆件k的rd阵的相应分块阵
			for (j = 0; j < 3; j++)//以下四行语句计算矩阵rd与矩阵tt相乘，结果存放在矩阵rdb中
				for (m = 0; m < 3; m++)
					for (n = 0; n < 3; n++)
						rdb[j][m] += rd[j][n] * tt[n][m];
			for (j = 0; j < 3; j++)//对三个分量循环
			{
				temp = 0;//临时变量归0
				for (m = 0; m < 3; m++)//对以下两行语句计算矩阵rdb与向量DON相乘，结果存放在矩阵temp中
					temp += rdb[j][m] * DON[t[i] + m];
				tref[j] += temp;//将计算出的杆端力分量叠加到相应位置
			}
		}
		else //固定节点的位移为0；矩阵相乘后仍为0
			for (j = 0; j < 3; j++)
				tref[j] += 0;//将0叠加到tref中
	}
	printf("hello world!12");
}
void PGaaa()
{
	printf("******************************************************************");
}
void PGPrint()
{
	int i, j; //循环控制变量
	printf(" 平面刚架结构计算\n");
	PGaaa();
	printf("\n 节点总数=%d\t\t固定节点数=%d\t\t可动节点数=%d\n 杆件数=%d\t\t荷载总数=%d\t\t截面总数=%d\n", TNN, NFIN, NFRN, NOR, NOL, NOS);
	PGaaa();
	printf("\n 节点号     X坐标       Y坐标\n");
	for (i = 1; i <= TNN; i++)
		printf("   %d\t   %5.7f\t%5.7f\n", i, XCN[i - 1], YCN[i - 1]);
	PGaaa();
	printf("\n 杆件号    左节点      右节点    弹性模量E    截面积A    主惯性距J\n");
	for (i = 1; i <= NOR; i++)
		printf("   %d\t     %d\t\t %d\t %5.0f     %5.3f    %5.8f\n", i, BNR[i - 1], ENR[i - 1], EEE[i - 1], AAA[i - 1], JJJ[i - 1]);
	PGaaa();
	printf("\n 荷载号   所在杆件号   荷载类型   荷载大小   距左端距离\n");
	for (i = 0; i < NOL; i++)
		printf("   %d\t     %d\t         %d\t  %5.2f      %5.5f\n", i + 1, NRL[i], KOL[i], VOL[i], DLB[i]);
	PGaaa();
	printf("\n 截面号   所在杆件号   距左端距离\n");
	for (i = 1; i <= NOS; i++)
		printf("   %d\t      %d\t       %5.7f\n", i, NRS[i - 1], DSB[i - 1]);
	PGaaa();
	printf("\n 节点号   位移X分量    位移Y分量     转角\n");
	for (i = NFIN + 1, j = 0; i <= TNN; i++, j++)
		printf("   %d\t  %5.7f   %5.7f   %5.7f\n", i, DON[3 * j], DON[3 * j + 1], DON[3 * j + 2]);
	PGaaa();
	printf("\n 截面号   内力X分量    内力Y分量     弯矩\n");
	for (i = 1; i <= NOS; i++)
		printf("   %d\t  %5.7f    %5.7f    %5.7f\n", i, IFS[3 * i - 3], IFS[3 * i - 2], IFS[3 * i - 1]);
	PGaaa();
}