globalelement.cpp

三维框架梁单元C++版本的源程序 可以使结构力学的概念更上升一个层次

	}
	fout<<endl;
}

void CGlobalElement::CalcuInternalForce()
{
	int NODE1,NODE2;
	int m_nTemp;
	//节点位移向量
	CMatrix m_matNodeDisp(12,1,0);
	CMatrix m_matNodeDisp1(12,1,0);
	CMatrix m_matNodeDisp2(12,1,0);

	for(int loop1=0;loop1<m_apEle[0]->m_nElementNum;loop1++)
	{
		NODE1 = m_apEle[0]->m_aBeamElemProValue[loop1].m_nElemNode1;
		NODE2 = m_apEle[0]->m_aBeamElemProValue[loop1].m_nElemNode2;

		for(int loop2=0;loop2<(m_Node.m_nNodeNum - m_Node.m_nNode_K_Num);loop2++)
		{
			if(NODE1 == m_Node.m_aDOFIndex[loop2].m_aiNode)
			{
				m_nTemp = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_X;
				m_matNodeDisp(0,0) = m_adDisp2[m_nTemp];
				m_nTemp = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_Y;				
				m_matNodeDisp(1,0) = m_adDisp2[m_nTemp];
				m_nTemp = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_Z;
				m_matNodeDisp(2,0) = m_adDisp2[m_nTemp];
				m_nTemp = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_RX;
				m_matNodeDisp(3,0) = m_adDisp2[m_nTemp];
				m_nTemp = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_RY;
				m_matNodeDisp(4,0) = m_adDisp2[m_nTemp];
				m_nTemp = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_RZ;
				m_matNodeDisp(5,0) = m_adDisp2[m_nTemp];
			}
			if(NODE2 == m_Node.m_aDOFIndex[loop2].m_aiNode)
			{
				m_nTemp = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_X;
				m_matNodeDisp(6,0) = m_adDisp2[m_nTemp];
				m_nTemp = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_Y;
				m_matNodeDisp(7,0) = m_adDisp2[m_nTemp];
				m_nTemp = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_Z;
				m_matNodeDisp(8,0) = m_adDisp2[m_nTemp];
				m_nTemp = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_RX;
				m_matNodeDisp(9,0) = m_adDisp2[m_nTemp];
				m_nTemp = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_RY;
				m_matNodeDisp(10,0) = m_adDisp2[m_nTemp];
				m_nTemp = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_RZ;
				m_matNodeDisp(11,0) = m_adDisp2[m_nTemp];
			}
		}

		//求解单元内力
		CMatrix tempMK1,tempMK2;
		tempMK1 = m_apEle[0]->m_CaElemOrigMatrix[loop1];
		tempMK2 = m_apEle[0]->m_CaElemTranFactorMatrix[loop1];
		//
		m_matNodeDisp1 = tempMK2 * m_matNodeDisp;	
		m_matNodeDisp2 = tempMK1 * m_matNodeDisp1;  		

		//I节点*****
		//轴力
		m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemFX[0] += -m_matNodeDisp2(0,0);
		//剪力Y
		m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemShearY[0] += -m_matNodeDisp2(1,0);
		//剪力Z
		m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemShearZ[0] += -m_matNodeDisp2(2,0);		
		//扭距
		m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemT[0] += -m_matNodeDisp2(3,0);
		//弯距Y
		m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemMY[0] += -m_matNodeDisp2(4,0);
		//弯距Z		
		m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemMZ[0] += -m_matNodeDisp2(5,0);
		
		//J节点*****
		//轴力
		m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemFX[12] += m_matNodeDisp2(6,0);
		//剪力Y
		m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemShearY[12] += m_matNodeDisp2(7,0);
		//剪力Z	
		m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemShearZ[12] += m_matNodeDisp2(8,0);
		//扭距
		m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemT[12] += m_matNodeDisp2(9,0);
		//弯距Y
		m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemMY[12] += m_matNodeDisp2(10,0);
		//弯距Z
		m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemMZ[12] += m_matNodeDisp2(11,0);   		

		//节点力均分
		//轴力
		for(loop2=1;loop2<12;loop2++)
		m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemFX[loop2] = 
			m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemFX[12];
		//扭距
		for(loop2=1;loop2<12;loop2++)
		m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemT[loop2] = 
			m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemT[12];
		//剪力Y
		for(loop2=1;loop2<12;loop2++)
		m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemShearY[loop2] = 
			m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemShearY[0];
		//剪力Z
		for(loop2=1;loop2<12;loop2++)
		m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemShearZ[loop2] = 
			m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemShearZ[0];
		//弯距MY
		double MY1 = m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemMY[0];
		double MY2 = m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemMY[12];
		double MY3 = MY2 - MY1;
		for(loop2=1;loop2<12;loop2++)
			m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemMY[loop2] = MY1 + loop2*MY3/12;
		//弯距MZ
		double MZ1 = m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemMZ[0];
		double MZ2 = m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemMZ[12];
		double MZ3 = MZ2 - MZ1;
		for(loop2=1;loop2<12;loop2++)
			m_apEle[0]->m_aElemInForceDef_N[loop1].m_dElemMZ[loop2] = MZ1 + loop2*MZ3/12;			
	}


	//叠和	
	for(int loop=0;loop<m_apEle[0]->m_nElementNum;loop++)
	{
		for(int loop1=0;loop1<13;loop1++)
		{
		m_apEle[0]->m_aElemInForceDef[loop].m_dElemFX[loop1] += m_apEle[0]->m_aElemInForceDef_N[loop].m_dElemFX[0];
		m_apEle[0]->m_aElemInForceDef[loop].m_dElemShearY[loop1] += m_apEle[0]->m_aElemInForceDef_N[loop].m_dElemShearY[0];
		m_apEle[0]->m_aElemInForceDef[loop].m_dElemShearZ[loop1] += m_apEle[0]->m_aElemInForceDef_N[loop].m_dElemShearZ[0];
		m_apEle[0]->m_aElemInForceDef[loop].m_dElemT[loop1] += m_apEle[0]->m_aElemInForceDef_N[loop].m_dElemT[loop1];
		m_apEle[0]->m_aElemInForceDef[loop].m_dElemMY[loop1] += m_apEle[0]->m_aElemInForceDef_N[loop].m_dElemMY[loop1];
		m_apEle[0]->m_aElemInForceDef[loop].m_dElemMZ[loop1] += m_apEle[0]->m_aElemInForceDef_N[loop].m_dElemMZ[loop1];
		}
	}       

	//输出单元内力
	OutputInternalForce();
}

void CGlobalElement::OutputInternalForce()
{
		//输出单元内力
		ofstream fout;
		fout.open("InternalForce.twm");
		
		fout.width(16);
		fout<<"Element No";
		fout.width(16);
		fout<<"Section No";
		fout.width(16);
		fout<<"FX";
		fout.width(16);
		fout<<"ShearY";
		fout.width(16);
		fout<<"ShearZ";
		fout.width(16);
		fout<<"MX";
		fout.width(16);
		fout<<"MY";
		fout.width(16);
		fout<<"MZ"<<endl;

		for(int loop=0;loop<m_apEle[0]->m_nElementNum;loop++)
		{			
			for(int loop1=0;loop1<13;loop1++)
			{
				fout.width(16);
				fout<<m_apEle[0]->m_aBeamElemProValue[loop].m_nElemNode1;
				fout.width(16);
				fout<<loop1+1;
				fout.width(16);
				fout<<m_apEle[0]->m_aElemInForceDef[loop].m_dElemFX[loop1];
				fout.width(16);
				fout<<m_apEle[0]->m_aElemInForceDef[loop].m_dElemShearY[loop1];
				fout.width(16);
				fout<<m_apEle[0]->m_aElemInForceDef[loop].m_dElemShearZ[loop1];
				fout.width(16);
				fout<<m_apEle[0]->m_aElemInForceDef[loop].m_dElemT[loop1];
				fout.width(16);
				fout<<m_apEle[0]->m_aElemInForceDef[loop].m_dElemMY[loop1];
				fout.width(16);
				fout<<m_apEle[0]->m_aElemInForceDef[loop].m_dElemMZ[loop1]<<endl;	
			}
			fout<<endl;
		}

		fout.close();
}

void CGlobalElement::CalcuSupportReaction()
{
	//AfxMessageBox(_T("CalcuSupportReaction()"),MB_OK,1);
	//临时变量		
	int nTotalDOF = m_Node.m_nTotalDOF;
	int nFreeDOF = m_Node.m_nFreeDOF;
	int Num = nTotalDOF - nFreeDOF;
	double dBuf;	

	//定义并且初始化	
	double *SupportReaction = new double [nTotalDOF];
	for(int loop=0;loop<nTotalDOF;loop++)
	{
		SupportReaction[loop] = 0;	
	}

	//输出支座反力过程数据和结果数据
	ofstream fout;
	fout.open("SupportReaction.twm");

	/*
	fout<<"输出整体坐标下节点位移向量"<<endl;
	for(loop=0;loop<nTotalDOF;loop++)
	{
		fout<<m_adDisp2[loop]<<endl;
	}
	fout<<endl;
	fout<<endl;
	fout<<endl;

	fout<<"输出整体坐标下外荷载向量"<<endl;
	for(loop=0;loop<nTotalDOF;loop++)
	{
		fout<<LoadVector[loop]<<endl;
	}
	fout<<endl;
	fout<<endl;
	fout<<endl;	   */


	//计算支座反力
	for(loop=0;loop<nTotalDOF;loop++)
	{
		for(int loop1=0;loop1<nTotalDOF;loop1++)
		{			
			dBuf = m_apEle[0]->m_smatGK(loop,loop1);	
			SupportReaction[loop] += dBuf*m_adDisp2[loop1]; 			
		}		
	}

	//支座反力 ＝ 总体力的向量 － 外荷载向量
	for(loop=0;loop<nTotalDOF;loop++)
	{
		SupportReaction[loop] = SupportReaction[loop] - LoadVector[loop];
		//规范数据，没有无穷小数据
		if( (SupportReaction[loop]<0.000000001) && (SupportReaction[loop]>-0.000000001) )
			SupportReaction[loop] = 0;
	}
	
	/*
	//fout<<"总体座标下支座反力的向量"<<endl;
	for(loop=0;loop<nTotalDOF;loop++)
	{
		fout<<SupportReaction[loop]<<endl;
	}
	fout<<endl;
	fout<<endl;
	fout<<endl;	
	*/


	//生成支座反力（过程变量）

	//初始化
	double **Reaction = new double *[m_Node.m_nNodeNum - m_Node.m_nNode_K_Num];
	for(loop=0;loop<(m_Node.m_nNodeNum - m_Node.m_nNode_K_Num);loop++)
		Reaction[loop]  = new double [7];
	for(loop=0;loop<(m_Node.m_nNodeNum - m_Node.m_nNode_K_Num);loop++)
	{
		for(int loop1=0;loop1<7;loop1++)
			Reaction[loop][loop1] = 0;
	}

	//Reaction具体赋值
	for(loop=0;loop<(m_Node.m_nNodeNum - m_Node.m_nNode_K_Num);loop++)
	{
		Reaction[loop][0] = int(  m_Node.m_aDOFIndex[loop].m_aiNode  );	
		Reaction[loop][1] = SupportReaction[m_Node.m_aDOFIndex[loop].m_aiDOFIndex_X];
		Reaction[loop][2] = SupportReaction[m_Node.m_aDOFIndex[loop].m_aiDOFIndex_Y];
		Reaction[loop][3] = SupportReaction[m_Node.m_aDOFIndex[loop].m_aiDOFIndex_Z];
		Reaction[loop][4] = SupportReaction[m_Node.m_aDOFIndex[loop].m_aiDOFIndex_RX];
		Reaction[loop][5] = SupportReaction[m_Node.m_aDOFIndex[loop].m_aiDOFIndex_RY];
		Reaction[loop][6] = SupportReaction[m_Node.m_aDOFIndex[loop].m_aiDOFIndex_RZ];
	}

	//由Reaction向m_apEle[0]->m_nSupportReaction转化
	int num = 0;
	double m_temp;
	for(loop=0;loop<(m_Node.m_nNodeNum - m_Node.m_nNode_K_Num);loop++)
	{
		m_temp = 0;
		m_temp += (Reaction[loop][1])*(Reaction[loop][1]);
		m_temp += (Reaction[loop][2])*(Reaction[loop][2]);
		m_temp += (Reaction[loop][3])*(Reaction[loop][3]);
		m_temp += (Reaction[loop][4])*(Reaction[loop][4]);
		m_temp += (Reaction[loop][5])*(Reaction[loop][5]);
		m_temp += (Reaction[loop][6])*(Reaction[loop][6]);

		if(m_temp>0.000000001)
		{
			m_apEle[0]->m_nSupportReaction[num].m_nNodeNo = Reaction[loop][0];
			m_apEle[0]->m_nSupportReaction[num].m_dSR_X =   Reaction[loop][1];
			m_apEle[0]->m_nSupportReaction[num].m_dSR_Y =   Reaction[loop][2];
			m_apEle[0]->m_nSupportReaction[num].m_dSR_Z =   Reaction[loop][3];
			m_apEle[0]->m_nSupportReaction[num].m_dSR_RX =  Reaction[loop][4];
			m_apEle[0]->m_nSupportReaction[num].m_dSR_RY =  Reaction[loop][5];
			m_apEle[0]->m_nSupportReaction[num].m_dSR_RZ =  Reaction[loop][6];
			num++;
		}			
	}	

	fout<<"支座反力："<<endl;
	fout.width(15);
	fout<<"节点号";
	fout.width(15);
	fout<<"FX";
	fout.width(15);
	fout<<"FY";
	fout.width(15);
	fout<<"FZ";
	fout.width(15);
	fout<<"MX";
	fout.width(15);
	fout<<"MY";
	fout.width(15);
	fout<<"MZ"<<endl;
	for(loop=0;loop<m_Node.m_nNodeFixNum;loop++)
	//for(loop=0;loop<(m_Node.m_nNodeNum - m_Node.m_nNode_K_Num);loop++)
	{
		fout.width(15);
		//fout<<num;
		fout<<m_apEle[0]->m_nSupportReaction[loop].m_nNodeNo;
		//fout<<Reaction[loop][0];
		fout.width(15);
		fout<<m_apEle[0]->m_nSupportReaction[loop].m_dSR_X;
		//fout<<Reaction[loop][1];
		fout.width(15);
		fout<<m_apEle[0]->m_nSupportReaction[loop].m_dSR_Y;
		//fout<<Reaction[loop][2];
		fout.width(15);
		fout<<m_apEle[0]->m_nSupportReaction[loop].m_dSR_Z;
		//fout<<Reaction[loop][3];
		fout.width(15);
		fout<<m_apEle[0]->m_nSupportReaction[loop].m_dSR_RX;
		//fout<<Reaction[loop][4];
		fout.width(15);
		fout<<m_apEle[0]->m_nSupportReaction[loop].m_dSR_RY;
		//fout<<Reaction[loop][5];
		fout.width(15);
		fout<<m_apEle[0]->m_nSupportReaction[loop].m_dSR_RZ<<endl;
		//fout<<Reaction[loop][6]<<endl;
	}
	fout.close();
	
	//删除过程变量
	delete Reaction;
	delete SupportReaction;
}