main.cpp

结构矩阵分析

  
  

 //跨内集中荷载的等效节点力的计算
struct equivalent_jiedian_force_focus *focus_equal=new struct equivalent_jiedian_force_focus[number_of_focus_load];
double a,b;
 for(i=0;i<number_of_focus_load;i++)
 {
	 a=(*(changdu+((focus_load+i)->serial_of_element-1)))*((focus_load+i)->local_of_force)*(-1);
	 b=(*(changdu+((focus_load+i)->serial_of_element-1)))-a;
(focus_equal+i)->Nx1_x=(-1)*0.5*((focus_load+i)->fx)*(*(cos+((focus_load+i)->serial_of_element-1)))*(*(cos+((focus_load+i)->serial_of_element-1)));
(focus_equal+i)->Nx1_y=(-1)*0.5*((focus_load+i)->fx)*(*(cos+((focus_load+i)->serial_of_element-1)))*(*(sin+((focus_load+i)->serial_of_element-1)));
(focus_equal+i)->Nx2_x=0.5*((focus_load+i)->fx)*(*(cos+((focus_load+i)->serial_of_element-1)))*(*(cos+((focus_load+i)->serial_of_element-1)));
(focus_equal+i)->Nx2_y=0.5*((focus_load+i)->fx)*(*(cos+((focus_load+i)->serial_of_element-1)))*(*(sin+((focus_load+i)->serial_of_element-1)));
(focus_equal+i)->Qx1_x=(-1)*((focus_load+i)->fx)*(*(sin+((focus_load+i)->serial_of_element-1)))*pow(b,2)*((*(changdu+((focus_load+i)->serial_of_element-1)))+a*2)/pow(*(changdu+((focus_load+i)->serial_of_element-1)),3)*(*(sin+((focus_load+i)->serial_of_element-1)));
(focus_equal+i)->Qx1_y=(-1)*((focus_load+i)->fx)*(*(sin+((focus_load+i)->serial_of_element-1)))*pow(b,2)*((*(changdu+((focus_load+i)->serial_of_element-1)))+a*2)/pow(*(changdu+((focus_load+i)->serial_of_element-1)),3)*(*(cos+((focus_load+i)->serial_of_element-1)));
(focus_equal+i)->Qx2_x=((focus_load+i)->fx)*(*(sin+((focus_load+i)->serial_of_element-1)))*pow(a,2)*(*(changdu+((focus_load+i)->serial_of_element-1))+b*2)/pow(*(changdu+((focus_load+i)->serial_of_element-1)),3)*(*(sin+((focus_load+i)->serial_of_element-1)));
(focus_equal+i)->Qx2_y=((focus_load+i)->fx)*(*(sin+((focus_load+i)->serial_of_element-1)))*pow(a,2)*(*(changdu+((focus_load+i)->serial_of_element-1))+b*2)/pow(*(changdu+((focus_load+i)->serial_of_element-1)),3)*(*(cos+((focus_load+i)->serial_of_element-1)));
(focus_equal+i)->Mx1=((focus_load+i)->fx)*(*(sin+((focus_load+i)->serial_of_element-1)))*a*pow(b,2)/pow(*(changdu+((focus_load+i)->serial_of_element-1)),2);
(focus_equal+i)->Mx2=(-1)*((focus_load+i)->fx)*(*(sin+((focus_load+i)->serial_of_element-1)))*b*pow(a,2)/pow(*(changdu+((focus_load+i)->serial_of_element-1)),2);

(focus_equal+i)->Ny1_x=0.5*(-1)*((focus_load+i)->fy)*(*(sin+((focus_load+i)->serial_of_element-1)))*(*(cos+((focus_load+i)->serial_of_element-1)));
(focus_equal+i)->Ny1_y=0.5*(-1)*((focus_load+i)->fy)*(*(sin+((focus_load+i)->serial_of_element-1)))*(*(sin+((focus_load+i)->serial_of_element-1)));
(focus_equal+i)->Ny2_x=0.5*(-1)*((focus_load+i)->fy)*(*(sin+((focus_load+i)->serial_of_element-1)))*(*(cos+((focus_load+i)->serial_of_element-1)));
(focus_equal+i)->Ny2_y=0.5*(-1)*((focus_load+i)->fy)*(*(sin+((focus_load+i)->serial_of_element-1)))*(*(sin+((focus_load+i)->serial_of_element-1)));
(focus_equal+i)->Qy1_x=(-1)*((focus_load+i)->fy)*(*(cos+((focus_load+i)->serial_of_element-1)))*pow(b,2)*((*(changdu+((focus_load+i)->serial_of_element-1)))+a*2)/pow(*(changdu+((focus_load+i)->serial_of_element-1)),3)*(*(sin+((focus_load+i)->serial_of_element-1)));
(focus_equal+i)->Qy1_y=((focus_load+i)->fy)*(*(cos+((focus_load+i)->serial_of_element-1)))*pow(b,2)*((*(changdu+((focus_load+i)->serial_of_element-1)))+a*2)/pow(*(changdu+((focus_load+i)->serial_of_element-1)),3)*(*(cos+((focus_load+i)->serial_of_element-1)));
(focus_equal+i)->Qy2_x=(-1)*((focus_load+i)->fy)*(*(cos+((focus_load+i)->serial_of_element-1)))*pow(a,2)*((*(changdu+((focus_load+i)->serial_of_element-1)))+b*2)/pow(*(changdu+((focus_load+i)->serial_of_element-1)),3)*(*(sin+((focus_load+i)->serial_of_element-1)));
(focus_equal+i)->Qy2_y=((focus_load+i)->fy)*(*(cos+((focus_load+i)->serial_of_element-1)))*pow(a,2)*((*(changdu+((focus_load+i)->serial_of_element-1)))+b*2)/pow(*(changdu+((focus_load+i)->serial_of_element-1)),3)*(*(cos+((focus_load+i)->serial_of_element-1)));
(focus_equal+i)->My1=((focus_load+i)->fy)*(*(cos+((focus_load+i)->serial_of_element-1)))*a*pow(b,2)/pow(*(changdu+((focus_load+i)->serial_of_element-1)),2);
(focus_equal+i)->My2=(-1)*((focus_load+i)->fy)*(*(cos+((focus_load+i)->serial_of_element-1)))*b*pow(a,2)/pow(*(changdu+((focus_load+i)->serial_of_element-1)),2);

//若考虑空间结构，还要在以下的程序段中加入杆和Z方向的夹角的正余弦值
/*(focus_equal+i)->Nz1=(-1)*((focus_load+i)->fz);
(focus_equal+i)->Nz2=(focus_load+i)->fz;
(focus_equal+i)->Qz1=(-1)*((focus_load+i)->fz)*pow(b,2)*((*(changdu+((focus_load+i)->serial_of_element-1)))+a*2)/pow(*(changdu+((focus_load+i)->serial_of_element-1)),3);
(focus_equal+i)->Qz2=((focus_load+i)->fz)*pow(a,2)*(*(changdu+((focus_load+i)->serial_of_element-1))+b*2)/pow(*(changdu+((focus_load+i)->serial_of_element-1)),3);
(focus_equal+i)->Mz1=((focus_load+i)->fz)*a*pow(b,2)/pow(*(changdu+((focus_load+i)->serial_of_element-1)),2);
(focus_equal+i)->Mz2=(-1)*((focus_load+i)->fz)*b*pow(a,2)/pow(*(changdu+((focus_load+i)->serial_of_element-1)),2);*/
 }

 
 //以下输出例子的等效节点力
/*for(i=0;i<number_of_spread_load;i++)//输出着整体坐标系x方向均布力的等效节点力
 {
	 cout<<"局部坐标系下的x方向均布力在单元起始端的等效轴力在整体坐标系x轴方向的分力是"
		 <<(spread_equal+i)->Nx1_x<<endl
		 <<"局部坐标系下的x方向均布力在单元起始端的等效轴力在整体坐标系x轴方向的分力是"
		 <<(spread_equal+i)->Nx1_y<<endl
		 <<"局部坐标系下的x方向均布力在单元终末端的等效轴力在整体坐标系x轴方向的分力是"
		 <<(spread_equal+i)->Nx2_x<<endl
		 <<"局部坐标系下的x方向均布力在单元终末端的等效轴力在整体坐标系y轴方向的分力是"
		 <<(spread_equal+i)->Nx2_y<<endl
		 <<"局部坐标系下的x方向均布力的单元起始端的等效切向力在整体坐标系x轴方向的分力是"
		 <<(spread_equal+i)->Qx1_x<<endl
		 <<"局部坐标系下的x方向均布力的单元起始端的等效切向力在整体坐标系y轴方向的分力是"
		 <<(spread_equal+i)->Qx1_y<<endl
		 <<"局部坐标系下的x方向均布力的单元终末端的等效切向力在整体坐标系x轴方向的分力是"
		 <<(spread_equal+i)->Qx2_x<<endl
         <<"局部坐标系下的x方向均布力的单元终末端的等效切向力在整体坐标系y轴方向的分力是"
		 <<(spread_equal+i)->Qx2_y<<endl
         <<"局部坐标系下的x方向均布力的单元起始端的等效弯矩是"
		 <<(spread_equal+i)->Mx1<<endl
		 <<"局部坐标系下的x方向均布力的单元终末端的等效弯矩是"
		 <<(spread_equal+i)->Mx2<<endl;
 }
for(i=0;i<number_of_spread_load;i++)//输出整体坐标系y方向均布力的等效节点力
 {
	 cout<<"局部坐标系下的y方向均布力的单元起始端的等效轴力在整体坐标系x轴方向的分力是"
		 <<(spread_equal+i)->Ny1_x<<endl
         <<"局部坐标系下的y方向均布力的单元起始端的等效轴力在整体坐标系y轴方向的分力是"
		 <<(spread_equal+i)->Ny1_y<<endl
		 <<"局部坐标系下的y方向均布力的单元终末端的等效轴力在整体坐标系x轴方向的分力是"
		 <<(spread_equal+i)->Ny2_x<<endl
		 <<"局部坐标系下的y方向均布力的单元终末端的等效轴力在整体坐标系y轴方向的分力是"
		 <<(spread_equal+i)->Ny2_y<<endl
		 <<"局部坐标系下的y方向均布力的单元起始端的等效切向力在整体坐标系x轴方向的分力是"
		 <<(spread_equal+i)->Qy1_x<<endl
		 <<"局部坐标系下的y方向均布力的单元起始端的等效切向力在整体坐标系y轴方向的分力是"
		 <<(spread_equal+i)->Qy1_y<<endl
		 <<"局部坐标系下的y方向均布力的单元终末端的等效切向力在整体坐标系x轴方向的分力是"
		 <<(spread_equal+i)->Qy2_x<<endl
		 <<"局部坐标系下的y方向均布力的单元终末端的等效切向力在整体坐标系y轴方向的分力是"
		 <<(spread_equal+i)->Qy2_y<<endl
		 <<"局部坐标系下的y方向均布力的单元起始端的等效弯矩是"
		 <<(spread_equal+i)->My1<<endl
		 <<"局部坐标系下的y方向均布力的单元终末端的等效弯矩是"
		 <<(spread_equal+i)->My2<<endl;
 }
for(i=0;i<number_of_focus_load;i++)//输出着整体坐标系x方向集中力的等效节点力
 {
	 cout<<"局部坐标系下的x方向集中力的单元起始端的等效轴力在整体坐标系x轴方向的分力是"
		 <<(focus_equal+i)->Nx1_x<<endl
		 <<"局部坐标系下的x方向集中力的单元起始端的等效轴力在整体坐标系y轴方向的分力是"
		 <<(focus_equal+i)->Nx1_y<<endl
		 <<"局部坐标系下的x方向集中力的单元终末端的等效轴力在整体坐标系x轴方向的分力是"
		 <<(focus_equal+i)->Nx2_x<<endl
		 <<"局部坐标系下的x方向集中力的单元终末端的等效轴力在整体坐标系y轴方向的分力是"
		 <<(focus_equal+i)->Nx2_y<<endl
		 <<"局部坐标系下的x方向集中力的单元起始端的等效切向力在整体坐标系x轴方向的分力是"
		 <<(focus_equal+i)->Qx1_x<<endl
		 <<"局部坐标系下的x方向集中力的单元起始端的等效切向力在整体坐标系y轴方向的分力是"
		 <<(focus_equal+i)->Qx1_y<<endl
         <<"局部坐标系下的x方向集中力的单元终末端的等效切向力在整体坐标系x轴方向的分力是"
		 <<(focus_equal+i)->Qx2_x<<endl
		 <<"局部坐标系下的x方向集中力的单元终末端的等效切向力在整体坐标系y轴方向的分力是"
		 <<(focus_equal+i)->Qx2_y<<endl
		 <<"局部坐标系下的x方向集中力的单元起始端的等效弯矩是"
		 <<(focus_equal+i)->Mx1<<endl
         <<"局部坐标系下的x方向集中力的单元终末端的等效弯矩是"
		 <<(focus_equal+i)->Mx2<<endl;
 }
for(i=0;i<number_of_focus_load;i++)//输出着整体坐标系y方向集中力的等效节点力
 {
	 cout<<"局部坐标系下的y方向集中力的单元起始端的等效轴力在整体坐标系x轴方向的分力是"
		 <<(focus_equal+i)->Ny1_x<<endl
		 <<"局部坐标系下的y方向集中力的单元起始端的等效轴力在整体坐标系y轴方向的分力是"
		 <<(focus_equal+i)->Ny1_y<<endl
		 <<"局部坐标系下的y方向集中力的单元终末端的等效轴力在整体坐标系x轴方向的分力是"
		 <<(focus_equal+i)->Ny2_x<<endl
		 <<"局部坐标系下的y方向集中力的单元终末端的等效轴力在整体坐标系y轴方向的分力是"
		 <<(focus_equal+i)->Ny2_y<<endl
		 <<"局部坐标系下的y方向集中力的单元起始端的等效切向力在整体坐标系x轴方向的分力是"
		 <<(focus_equal+i)->Qy1_x<<endl
		 <<"局部坐标系下的y方向集中力的单元起始端的等效切向力在整体坐标系y轴方向的分力是"
		 <<(focus_equal+i)->Qy1_y<<endl
         <<"局部坐标系下的y方向集中力的单元终末端的等效切向力在整体坐标系x轴方向的分力是"
		 <<(focus_equal+i)->Qy2_x<<endl
		 <<"局部坐标系下的y方向集中力的单元终末端的等效切向力在整体坐标系y轴方向的分力是"
		 <<(focus_equal+i)->Qy2_y<<endl
		 <<"局部坐标系下的y方向集中力的单元起始端的等效弯矩是"
		 <<(focus_equal+i)->My1<<endl
         <<"局部坐标系下的y方向集中力的单元终末端的等效弯矩是"
		 <<(focus_equal+i)->My2<<endl;
 }*/




///////////////////////////////////////////////////////////////////////////////////////////
//																						///
//		组合总刚和总节点力向量--team 4 徐文权,胡徐强、陈维良、张林，王神莲				///
//																			6-2-2004	///
///////////////////////////////////////////////////////////////////////////////////////////
//考虑空间和时间效率
//一个单元的所占的空间很少,但总刚的所要的空间很大,考虑总刚的优化.
//考虑刚度矩阵的对称性,当整体自由度neq已知,那么只存储K上半个矩阵.
//所要的空间为(neq+1)*neq/2,K的类型定义为double*(一维数组)
	int e;
	double *K,*F;//为了总体刚度K，F的规模不受限制，定义为double*
	int numK=(neq+1)*neq/2;
	K=new double[numK];
	F=new double[neq];
	Element<6> *Elem;
	Elem=new Element<6>[nElem];
//初始化Ke,F
	for(i=0;i<neq;++i){
		F[i]=0;
	}
	for(i=0;i<numK;++i){
		K[i]=0;
	}
	for(e=0;e<nElem;e++){
		Elem[e].LM[0]=(e_df+e)->x1;
		Elem[e].LM[1]=(e_df+e)->y1;
		Elem[e].LM[2]=(e_df+e)->m1;
		Elem[e].LM[3]=(e_df+e)->x2;
		Elem[e].LM[4]=(e_df+e)->y2;
		Elem[e].LM[5]=(e_df+e)->m2;
	}
	for(e=0;e<nElem;e++){
		if((unit+e)->C=='T'){
			for(i=0;i<6;i++){
				Elem[e].Fe[i]=0.0;
			}
		}
		else{
			Elem[e].Fe[0]=(focus_equal+e)->Nx1_x+(focus_equal+e)->Qx1_x+(spread_equal+e)->Nx1_x+(spread_equal+e)->Qx1_x;
			Elem[e].Fe[1]=(focus_equal+e)->Nx1_y+(focus_equal+e)->Qx1_y+(spread_equal+e)->Nx1_y+(spread_equal+e)->Qx1_y;
			Elem[e].Fe[2]=(focus_equal+e)->Mx1+(focus_equal+e)->Mx1+(spread_equal+e)->Mx1+(spread_equal+e)->Mx1;
			Elem[e].Fe[3]=(focus_equal+e)->Nx2_x+(focus_equal+e)->Qx2_x+(spread_equal+e)->Nx2_x+(spread_equal+e)->Qx2_x;
			Elem[e].Fe[4]=(focus_equal+e)->Nx2_y+(focus_equal+e)->Qx2_y+(spread_equal+e)->Nx2_y+(spread_equal+e)->Qx2_y;
			Elem[e].Fe[5]=(focus_equal+e)->Mx2+(focus_equal+e)->Mx2+(spread_equal+e)->Mx2+(spread_equal+e)->Mx2;
		}
	}
	int nd=6,ki,kj,kk,j;
	for(e=0;e<nElem;++e){//for all element
		for(i=0;i<nd;i++){
			if(Elem[e].LM[i]>0){
				ki=Elem[e].LM[i];
				F[ki-1]+=Elem[e].Fe[i];
				for(j=i;j<nd;j++){
					if(Elem[e].LM[j]>0){
						kj=Elem[e].LM[j];
						if(ki<=kj)//当ki<kj是为下三角
						{	
							if(ki==0)kk=kj;
							else kk=(2*neq-ki+2)*(ki-1)/2+kj-ki;//原K[ki][kj]在现在K中的位置为kk
							if((unit+e)->C=='T')
								K[kk]+=(unit+e)->gan_matrix[i][j];
							else
								K[kk]+=(unit+e)->gangjia_matrix[i][j];
						}
					}
				}
			}
		}
	}
	for(i=0;i<neq;i++){
		j=(load+i)->serial_of_node;
		node_df tmp=*(df+j);
		F[tmp.x]+=(load+i)->fx;
		F[tmp.y]+=(load+i)->fy;
		F[tmp.m]+=(load+i)->mz;
	}
//输出刚度矩阵
	cout<<"总刚:"<<endl;
	for(j=0;j<neq;j++){
		cout.width(2);cout<<K[j]<<" ";
	}
	cout<<endl;
	for(i=1;i<neq;++i){
		cout.width(2);cout<<K[i]<<" ";
		for(j=1;j<i;++j){
			kk=(2*neq-j+1)*j/2+i-j;
			cout.width(2);cout<<K[kk]<<" ";
		}
		for(j=i;j<neq;++j){
			kk=(2*neq-i+1)*i/2+j-i;
			cout.width(2);cout<<K[kk]<<" ";
		}
		cout<<endl;
	}
	cout<<"总组合力:"<<endl;
	for(i=0;i<neq;i++){
		cout<<F[i]<<endl;
	}
	int l=1;//工况
	LDLTsolve(K,neq,F,l);
////////////////////////////////////////////////////////////////////////////////////////////////////
	//最后释放内存
	delete[] K,F,Elem,load,df,focus_equal,spread_equal;
}


//////////////////////////////////////////////////////////////////////////////////////

void main()
{
	work1_work2_work3();
	cout<<endl;
	
}