平面刚架有限元程序C语言完整版

//计算平面钢架的程序
//程序开始

#include
#include


#define NE 2
#define NJ 3                   //定义并输入基本参数
#define NZ 6
#define NPJ 3
#define NPF 2
#define NJ3 9
#define DD 9
#define E0 3.0000E8            // 定义并输入常数
#define A0 0.5
#define I0 4.1667E-2
#define PI 3.141592654


// 这是输入参数的初始化和定义全局变量

int jm[NE+1][3]={{000}{012}{031}};
double gc[NE+1]={0.05.05.0};
double gj[NE+1]={0.00.090.0};
double mj[NE+1]={0.0A0A0};
double gx[NE+1]={0.0I0I0};
int zc[NZ+1]={0456789};
double pj[NPJ+1][3]={{0.00.00.0}{0.06.01.0}{0.02.02.0}{0.0-5.03.0}};
double pf[NPF+1][5]={{00000}{0.0-4.85.01.01.0}{0.0-8.02.52.02.0}};
double kz[NJ3+1][NJ3+1]p[NJ3+1];
double pe[7]f[7]f0[7]t[7][7];
double ke[7][7]kd[7][7];
//  kz[][]---整体刚度矩阵
//  ke[][]---整体坐标系下的单元刚度矩阵
//  kd[][]---局部坐标系下的单元刚度矩阵
//  t[][] ---坐标板换矩阵


//*****函数声明*****
void jdugd（int）;
void zb（int）;
void gdnl（int）;
void dugd（int）;


//***主函数开始*****
void main（void）
{
     FILE *fp;
	int ijkedhhiijjhza1b1mldlzlzj0;
	double clwy[7];
	int IMINjn;


	//<功能：形成矩阵p>
	if（NPJ>0）
	{
		for（i=1;i<=NPJ;i++）
		{
			j=（int）pj[i][2];
			p[j]=pj[i][1];
		}
	}
	if（NPF>0）
	{
		for（i=1;i<=NPF;i++）
		{
			hz=i;  //可否省略此行代码，下一行直接写gdnl（i），这样也是将实参i传递给形参hz
			gdnl（hz）;
			e=（int）pf[hz][3];
			zb（e）;
			for（j=1;j<=6;j++）
			{
				pe[j]=0.0;
				for（k=1;k<=6;k++）
				{
					pe[j]=pe[j]-t[k][j]*f0[k];    //用的是转秩矩阵     zb（）函数出来的是直接矩阵   在这先让列不动行动  相当于乘的转秩
				}
			}
			a1=jm[e][1];
			b1=jm[e][2];
			p[3*a1-2]=p[3*a1-2]+pe[1];
			p[3*a1-1]=p[3*a1-1]+pe[2];
			p[3*a1]=p[3*a1]+pe[3];
			p[3*b1-2]=p[3*b1-2]+pe[4];
			p[3*b1-1]=p[3*b1-1]+pe[5];
			p[3*b1]=p[3*b1]+pe[6];

		}
	}






//*********************************************
	//<功能：生成整体刚度矩阵ke[][]>
	for（e=1;e<=NE;e++）                            //这里的总刚KE实际是已经处理过的带状矩阵KZ*了!
	{

		dugd（e）;
		for（i=1;i<=2;i++）
		{
			for（ii=1;ii<=3;ii++）
			{
				h=3*（i-1）+ii;
				dh=3*（jm[e][i]-1）+ii;
				
			
		
				for（j=1;j<=2;j++）
				{
					for（jj=1;jj<=3;jj++）
					{
						l=3*（j-1）+jj;
						zl=3*（jm[e][j]-1）+jj;
						dl=zl-dh+1;    
						if（dl>0）
							kz[dh][dl]=kz[dh][dl]+ke[h][l];
					}
				}
			}
		}
	}

	//****引入边界条件*******
	for（i=1;i<=NZ;i++）    //按支撑循环，即约束总数（假设共一个固定端支座和一个固定铰支座，则总数仍为6?）
	{
		z=zc[i];          //取出支撑对应的位移数（如固定铰支座对应节点码为3，则z=7 8 9）放入变量z中
		kz[z][1]=1.0;     //Z行第一个元素置为1.0，注意这里的KZ已是处理后的带状矩阵了，即对应于原总刚中的Z行中对角线处的元素置一了!
		for（j=2;j<=DD;j++）  //这个循环是将带状矩阵中Z行的其他元素全部置为零，由于KZ与KZ*中行元素是不变的，即KZ中的r行元素在KZ*中
			kz[z][j]=0.0;   //仍在r行。即这个循环的本质是将KZ总刚中Z行除对角线位置元素外的所有元素置为零!
		if（（z!=1））
		{
			if（z>DD）j0=DD;    //先通过一个判断语句，得出正确的列数J0，要用这个列数来实现Z行斜向上45°所有元素置零，这个操作等同于
			else if（z<=DD）j0=z;  //将总刚KZ中的Z列除对角线元素外的所有元素置零!
			for（j=2;j<=j0;j++）    //j从2到j0代表了斜向上45°各元素的列变化规律
				kz[z-j+1][j]=0.0;  //代表了斜向上45°各元素的行变化规律
		}
		p[z]=0.0;   //将位移约束代码所对应的荷载列阵处的元素置为零                          
	}



	/*高斯