youixanyuan.cpp

有限元算法。平面绗架结构的有限元分析

 
    //用高斯主列消去法解方程
    //计算增广矩阵的变形矩阵
	for(Ti=0;Ti<NJ*3-1;Ti++)
	{
		max=fabs(gd_kz[Ti][Ti]);
		i_max=Ti;
		for(Ti2=Ti;Ti2<NJ*3;Ti2++)
		{
			if(fabs(gd_kz[Ti2][Ti])>max)
			{
				i_max=Ti2-1;
				max=gd_kz[Ti2][Ti];				
			}
		}
       
		if(i_max!=Ti)
		{
			for(i=Ti;i<NJ*3+1;i++)
			{
				zhongjian=gd_kz[Ti][i];
				gd_kz[Ti][i]=gd_kz[i_max][i];
				gd_kz[i_max][i]=zhongjian;
			}
		}
		
		for(j=Ti+1;j<NJ*3;j++)
		{	
			double shang=gd_kz[j][Ti]/gd_kz[Ti][Ti];
			for(k=0;k<NJ*3+1;k++)
				gd_kz[j][k]=-shang*gd_kz[Ti][k]+gd_kz[j][k];
		}	
	}    

	//反代求值
	double chengji=0;
	weiyi[NJ*3-1]=gd_kz[NJ*3-1][NJ*3]/gd_kz[NJ*3-1][NJ*3-1];
	for(i=NJ*3-2;i>=0;i--)
	{
		for(j=i+1;j<NJ*3;j++)
			chengji=chengji+gd_kz[i][j]*weiyi[j];
		weiyi[i]=(gd_kz[i][NJ*3]-chengji)/gd_kz[i][i];
		chengji=0;
	}
}











//********************
//引入支承条件的函数
//********************
void zhichengyinru()
{
	int i,Ti;
	int weizhi;
	fangchengshu=NJ*3-NZ;
	for(i=0;i<NZ;i++)
	{
		weizhi=zhichengweizhi[i]-1;
		for(Ti=0;Ti<NJ*3;Ti++)
			zhengtigangdu[weizhi][Ti]=0;
        for(Ti=0;Ti<NJ*3;Ti++)
			zhengtigangdu[Ti][weizhi]=0;
		zhengtigangdu[weizhi][weizhi]=1;
        F[weizhi]=0;
	}
}




//*************************
//计算等效总节点载荷的函数
//*************************
void zongzaihe()	
{

	int i,Ti;
	double U1=0,V1=0,M1=0,U2=0,V2=0,M2=0;
	double c,G,l,d;
	double Ff[6];
	double T[6][6],T1[6][6];          //T为坐标转换矩阵,T1为T的转置矩阵
	for(i=0;i<NJ*3;i++)
		F[i]=0;
	for(i=0;i<6;i++)
		Ff[i]=0;
    for( i=0;i<NE;i++)
	{
		for(int jl=0;jl<6;jl++)
			dy_weiyi[i][jl]=0;
	}
	for(i=0;i<NPF;i++)
	{
		int dym=fjzhzuoyongdanyuan[i];
		U1=V1=M1=U2=V2=M2=0;
		for(int t=0;t<6;t++)
			Ff[t]=0;
		c=fjzhzuoyongdian[i];
		l=changdu[(fjzhzuoyongdanyuan[i]-1)];
		d=l-c;
		G=fjzhzhi[i];
		
		if(fjzhleixing[i]==1)
		{
			V1=-0.5*G*c*(2-2*c*c/(l*l)+c*c*c/(l*l*l));
			V2=-G*c-V1;
			M1=G*c*c*(6-8*c/l+3*c*c/(l*l))/12;
			M2=-G*c*c*c*(4-3*c/l)/(12*l);
			U1=U2=0;		
		}		
		else if(fjzhleixing[i]==2)
		{
			V1=-G*(l+2*c)*d*d/(l*l*l);
			V2=-G*(l+2*d)*c*c/(l*l*l);
			M1=G*c*d*d/(l*l);
			M2=-G*c*c*d/(l*l);
			U1=U2=0;

		}
		else if(fjzhleixing[i]==3)
		{
			U1=-G*d/l;
			U2=-G*c/l;
			M1=M2=V1=V2=0;
		}
		else if(fjzhleixing[i]==4)
		{
			U1=U2=0;
			V1=-6*G*c*d/(l*l*l);
			V2=-V1;
			M1=G*d*(2-3*d/l)/l;
			M2=G*c*(2-3*c/l)/l;
		}
		else if(fjzhleixing[i]==5)
		{
			V1=-G*sin(fjzhjiaodu[i])*(l+2*c)*d*d/(l*l*l);
			V2=-G*sin(fjzhjiaodu[i])*(l+2*d)*c*c/(l*l*l);
			M1=G*sin(fjzhjiaodu[i])*c*d*d/(l*l);
			M2=-G*sin(fjzhjiaodu[i])*c*c*d/(l*l);
			U1=U2=0;
			U1=-G*cos(fjzhjiaodu[i])*d/l+U1;
			U2=-G*cos(fjzhjiaodu[i])*c/l+U2;	        
		}
	    

        /*
		cout<<"操你妈的;操你妈的;操你妈的;操你妈的;操你妈的:"<<endl;
		cout<<U1<<endl;
		cout<<V1<<endl;
		cout<<M1<<endl;
		cout<<U2<<endl;
		cout<<V2<<endl;
		cout<<M2<<endl;
		cout<<endl;
        */

		//记录竿的固端内力
		dy_neili[dym-1][0]=U1;
		dy_neili[dym-1][1]=V1;
		dy_neili[dym-1][2]=M1;
		dy_neili[dym-1][3]=U2;
		dy_neili[dym-1][4]=V2;
		dy_neili[dym-1][5]=M2;
		
		//cout<<"你输入的载荷类型不正确,非节点载荷类型:1-均布,2-垂直集中,3-平行集中,4-力偶,5-角度集中.输入其他值都是不正确的"<<endl;
	
        //计算坐标转换矩阵T及其逆阵T1
		js_T_T1( fjzhzuoyongdanyuan[i]-1, T,T1);
		//cout<<"单元序号:"<<fjzhzuoyongdanyuan[i]-1<<endl;
        //cout<<"角度:"<<jiaodu[fjzhzuoyongdanyuan[i]-1]<<endl;
		//cout<<"操："<<endl;
		for(Ti=0;Ti<6;Ti++)
			Ff[Ti]=-(T1[Ti][0]*U1+T1[Ti][1]*V1+T1[Ti][2]*M1+T1[Ti][3]*U2+T1[Ti][4]*V2+T1[Ti][5]*M2);		
		
		//for(int h=0;h<6;h++)
		//	cout<<Ff[h]<<endl;
		//cout<<endl;
		
		    //************************
		    //将载荷转换到整体坐标系中
		    //************************
		    
			F[(dym_jdm[dym-1][0]-1)*3]=F[(dym_jdm[dym-1][0]-1)*3]+Ff[0];
			F[(dym_jdm[dym-1][0]-1)*3+1]=F[(dym_jdm[dym-1][0]-1)*3+1]+Ff[1];
			F[(dym_jdm[dym-1][0]-1)*3+2]=F[(dym_jdm[dym-1][0]-1)*3+2]+Ff[2];
			F[(dym_jdm[dym-1][1]-1)*3]=F[(dym_jdm[dym-1][1]-1)*3]+Ff[3];
			F[(dym_jdm[dym-1][1]-1)*3+1]=F[(dym_jdm[dym-1][1]-1)*3+1]+Ff[4];
			F[(dym_jdm[dym-1][1]-1)*3+2]=F[(dym_jdm[dym-1][1]-1)*3+2]+Ff[5];
			/*
			F[(dym_jdm[fjzhzuoyongdanyuan[i]-1][0]-1)*3]=F[(dym_jdm[fjzhzuoyongdanyuan[i]-1][0]-1)*3]+Ff[0];
			F[(dym_jdm[fjzhzuoyongdanyuan[i]-1][0]-1)*3+1]=F[(dym_jdm[fjzhzuoyongdanyuan[i]-1][0]-1)*3+1]+Ff[1];
			F[(dym_jdm[fjzhzuoyongdanyuan[i]-1][0]-1)*3+2]=F[(dym_jdm[fjzhzuoyongdanyuan[i]-1][0]-1)*3+2]+Ff[2];
			F[(dym_jdm[fjzhzuoyongdanyuan[i]-1][1]-1)*3]=F[(dym_jdm[fjzhzuoyongdanyuan[i]-1][1]-1)*3]+Ff[3];
			F[(dym_jdm[fjzhzuoyongdanyuan[i]-1][1]-1)*3+1]=F[(dym_jdm[fjzhzuoyongdanyuan[i]-1][1]-1)*3+1]+Ff[4];
			F[(dym_jdm[fjzhzuoyongdanyuan[i]-1][1]-1)*3+2]=F[(dym_jdm[fjzhzuoyongdanyuan[i]-1][1]-1)*3+2]+Ff[5];
			*/
	}
    

	for(i=0;i<NPJ;i++)
	{
		F[(jdzhzuoyongdian[i]-1)*3]=F[(jdzhzuoyongdian[i]-1)*3]+jiedianzaihe[i][0];
		F[(jdzhzuoyongdian[i]-1)*3+1]=F[(jdzhzuoyongdian[i]-1)*3+1]+jiedianzaihe[i][1];
        F[(jdzhzuoyongdian[i]-1)*3+2]=F[(jdzhzuoyongdian[i]-1)*3+2]+jiedianzaihe[i][2];
	}

}
                     
//***********************************************
//计算总刚度的函数,存放于zhengtigangdu[][]数组中
//***********************************************
void zonggang()
{
	double gd[6][6],gd2[6][6];        // gd1为中间变量.
	double T[6][6],T1[6][6];          //T为坐标转换矩阵,T1为T的转置矩阵
	double danyuan_gangdu[6][6];
    int Ti,Ti2,Ti3;
	double chengji=0;                 //中间变量

    for(Ti=0;Ti<NJ*3;Ti++)            //先将整体刚度付0
	{
		for(Ti2=0;Ti2<NJ*3;Ti2++)
			zhengtigangdu[Ti][Ti2]=0;
	}


	for(int i=0;i<NE;i++)
	{
		//***********************
    	//局部坐标系中的单元刚度
    	//***********************
		dy_gangdu(i,danyuan_gangdu);

		//********************
		//计算单元的转换矩阵T
		//********************
		
		js_T_T1( i, T,T1);
        
        //****************************************
		//计算单元在整体坐标系中的单刚
		//单元在整体坐标系中的刚度存在gd[]数组中.
        //****************************************
				
        for(Ti=0;Ti<6;Ti++)
		{
			for(Ti2=0;Ti2<6;Ti2++)
			{
				gd[Ti][Ti2]=0;
				gd2[Ti][Ti2]=0;
			}

		}
        chengji=0;

		for(Ti=0;Ti<6;Ti++)
		{
			for(Ti2=0;Ti2<6;Ti2++)
			{
				for(Ti3=0;Ti3<6;Ti3++)
					chengji=chengji+T1[Ti][Ti3]*danyuan_gangdu[Ti3][Ti2];    //改过
				gd2[Ti][Ti2]=chengji;
				chengji=0;
			}
			
		}
        chengji=0;

		for(Ti=0;Ti<6;Ti++)
		{
			for(Ti2=0;Ti2<6;Ti2++)
			{
				for(Ti3=0;Ti3<6;Ti3++)
					chengji=chengji+gd2[Ti][Ti3]*T[Ti3][Ti2];       //改过
                gd[Ti][Ti2]=chengji;
				chengji=0;
			}			
		}
		chengji=0;

        //*******************************************
		//将每个单元刚度集成到整体刚度中计算整体刚度
		//*******************************************

        
		int qianjiedian,houjiedian;
		qianjiedian=dym_jdm[i][0];
		houjiedian=dym_jdm[i][1];
		for(Ti=0;Ti<3;Ti++)
		{
			for(Ti2=0;Ti2<3;Ti2++)
				zhengtigangdu[(qianjiedian-1)*3+Ti][(qianjiedian-1)*3+Ti2]=zhengtigangdu[(qianjiedian-1)*3+Ti][(qianjiedian-1)*3+Ti2]+gd[Ti][Ti2];
		}
		for(Ti=0;Ti<3;Ti++)
		{
			for(Ti2=0;Ti2<3;Ti2++)
				zhengtigangdu[(houjiedian-1)*3+Ti][(houjiedian-1)*3+Ti2]=zhengtigangdu[(houjiedian-1)*3+Ti][(houjiedian-1)*3+Ti2]+gd[Ti+3][Ti2+3];
		}
		for(Ti=0;Ti<3;Ti++)
		{
			for(Ti2=0;Ti2<3;Ti2++)
				zhengtigangdu[(qianjiedian-1)*3+Ti][(houjiedian-1)*3+Ti2]=zhengtigangdu[(qianjiedian-1)*3+Ti][(houjiedian-1)*3+Ti2]+gd[Ti][Ti2+3];
		}
		for(Ti=0;Ti<3;Ti++)
		{
			for(Ti2=0;Ti2<3;Ti2++)
				zhengtigangdu[(houjiedian-1)*3+Ti][(qianjiedian-1)*3+Ti2]=zhengtigangdu[(houjiedian-1)*3+Ti][(qianjiedian-1)*3+Ti2]+gd[Ti+3][Ti2];
		}
     }
}

//**************************
//输入数据的函数
//**************************
void input()
{
	char queding;
	char qdxuhao;
	int i_qdxuhao;
	//************
	//基本参数输入
	//************
	cout<<"请输入基本参数:单元数,节点数,支持数,节点载荷数,非节点载荷数;"<<endl<<endl;
	cout<<"单元数:";
	cin>>NE;
	cout<<"节点数:";
	cin>>NJ;
	cout<<"支持数:";
	cin>>NZ;
	cout<<"有节点载荷作用的节点数:";
	cin>>NPJ;
	cout<<"非节点载荷数:";
	cin>>NPF;
	cout<<endl<<endl;
	cout<<"单元数:"<<NE<<" 节点数:"<<NJ<<" 支持数:"<<NZ<<" 有节点载荷作用的节点数:"<<NPJ<<" 非节点载荷数:"<<NPF<<endl;
    cout<<"以上数据是否有错?确定有错,请输入Y,按其他任意键被视为无错."<<endl;
	cin>>queding;
	cout<<endl;
	while(queding=='y')
	{
		cout<<"你要修改第几个数据?请输入号数:";
		cin>>qdxuhao;
		switch(qdxuhao)
		{
		case '1':
			cout<<"重新输入单元数:";
	        cin>>NE;
			break;
		case '2':
			cout<<"重新输入节点数:";
	        cin>>NJ;
			break;
		case '3':
			cout<<"重新输入支持数:";
	        cin>>NZ;
			break;
		case '4':
			cout<<"重新输入有节点载荷作用的节点数:";
	        cin>>NPJ;
			break;