spaceframebeam.cpp

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

#include "stdafx.h"
//#include "FORCE.h"
#include "SpaceFrameBeam.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CSpaceFrameBeam::CSpaceFrameBeam()
{

}

CSpaceFrameBeam::~CSpaceFrameBeam()
{

}

void CSpaceFrameBeam::ReadFile(CString m_cFilePath)
{
	ifstream fin;
	fin.open(m_cFilePath);

	// 临时参量	
	int m_nTempMat = 0; //临时材料号
	int m_nTempRC = 0;  //临时实常数号
	int m_N; 	
	char char1[100];
	bool m_bRUN = 0; // 控制运行的变量

	while(!fin.eof())
	{
		fin>>char1;	
    	if(!strcmp(char1,"ELEMENT_SPACEFRAME_START"))
			m_bRUN = 1;	

		while(m_bRUN)	
		{		
			fin>>char1;				
			if(!strcmp(char1,"ELEMENT_SPACEFRAME_END"))  
			{				
				m_bRUN = 0;
				break;			
			}
			//读取单元号
			if(!strcmp(char1,"EN"))
			{
				m_nElementNum++;
				m_aBeamElemProValue.SetSize(m_nElementNum);
				//单元序号
				fin>>char1;	
				m_N = atoi(char1);
				m_aBeamElemProValue[m_nElementNum-1].m_nElementNo = m_N;
				//单元节点1
				fin>>char1;	
				m_N = atoi(char1);
				m_aBeamElemProValue[m_nElementNum-1].m_nElemNode1 = m_N;
				//单元节点2
				fin>>char1;	
				m_N = atoi(char1);
				m_aBeamElemProValue[m_nElementNum-1].m_nElemNode2 = m_N;
				//单元节点3
				fin>>char1;	
				m_N = atoi(char1);
				m_aBeamElemProValue[m_nElementNum-1].m_nElemNode3 = m_N;			
				//材料号
				fin>>char1;	
				m_N = atoi(char1);
				m_aBeamElemProValue[m_nElementNum-1].m_nElementMat = m_N;
				//实常数号	
				fin>>char1;	
				m_N = atoi(char1);
				m_aBeamElemProValue[m_nElementNum-1].m_nElementRC = m_N;	
				//节点一释放
				fin>>char1;	
				m_N = atoi(char1);
				m_aBeamElemProValue[m_nElementNum-1].m_nEleRelease_I1 = m_N;
				fin>>char1;	
				m_N = atoi(char1);
				m_aBeamElemProValue[m_nElementNum-1].m_nEleRelease_I2 = m_N;
				fin>>char1;	
				m_N = atoi(char1);
				m_aBeamElemProValue[m_nElementNum-1].m_nEleRelease_I3 = m_N;
				fin>>char1;	
				m_N = atoi(char1);
				m_aBeamElemProValue[m_nElementNum-1].m_nEleRelease_I4 = m_N;
				fin>>char1;	
				m_N = atoi(char1);
				m_aBeamElemProValue[m_nElementNum-1].m_nEleRelease_I5 = m_N;
				fin>>char1;	
				m_N = atoi(char1);
				m_aBeamElemProValue[m_nElementNum-1].m_nEleRelease_I6 = m_N;
				//节点二释放
				fin>>char1;	
				m_N = atoi(char1);
				m_aBeamElemProValue[m_nElementNum-1].m_nEleRelease_J1 = m_N;
				fin>>char1;	
				m_N = atoi(char1);
				m_aBeamElemProValue[m_nElementNum-1].m_nEleRelease_J2 = m_N;
				fin>>char1;	
				m_N = atoi(char1);
				m_aBeamElemProValue[m_nElementNum-1].m_nEleRelease_J3 = m_N;
				fin>>char1;	
				m_N = atoi(char1);
				m_aBeamElemProValue[m_nElementNum-1].m_nEleRelease_J4 = m_N;
				fin>>char1;	
				m_N = atoi(char1);
				m_aBeamElemProValue[m_nElementNum-1].m_nEleRelease_J5 = m_N;
				fin>>char1;	
				m_N = atoi(char1);
				m_aBeamElemProValue[m_nElementNum-1].m_nEleRelease_J6 = m_N;
			}
		}
	}
	fin.close();	
}

void CSpaceFrameBeam::GetStiffness()
{	
	FillElementData();	
	FormElemMatrix();
	TranElemMatrix();	
	OutputParameter();
	InternalForceInitial();  //单元内力位移初始化
}

void CSpaceFrameBeam::StiffAssemble()
{
	CalcuGKBandWidth();
	CalcuGKDiagAddress();
	SparseMatrixStiffAssemble();
}

void CSpaceFrameBeam::FillElementData()
{
	//确定材料取值
	//临时变量
	int m_N1,m_N2;
	
	//材料类型的数量
	//实在费劲！
	CTypedPtrArray <CPtrArray,CBaseMaterial*> &apMaterial = *m_papMaterial;
	m_N1 = apMaterial[0]->m_nMatNum;
	for(int loop=0;loop<m_nElementNum;loop++)
	{		
		m_N2 = m_aBeamElemProValue[loop].m_nElementMat;		
		for(int loop1=0;loop1<m_N1;loop1++)
		{
			if(m_N2 == apMaterial[0]->m_aElasticMatPro[loop1].m_nNo)
			{
				//确定弹性模量
				m_aBeamElemProValue[loop].m_E = apMaterial[0]->m_aElasticMatPro[loop1].m_dEO[0];
				//确定剪切模量
				m_aBeamElemProValue[loop].m_G = apMaterial[0]->m_aElasticMatPro[loop1].m_dEO[4];
				//确定泊松比
				m_aBeamElemProValue[loop].m_Miu = apMaterial[0]->m_aElasticMatPro[loop1].m_dMiu;
			}
		} 
	}

	//确定实常数
	CTypedPtrArray <CPtrArray,CBaseRealConstant*> &apBaseRealConstant = *m_papBaseRealConstant;
	m_N1 = apBaseRealConstant[0]->m_nRealConstantNum;
	for(loop=0;loop<m_nElementNum;loop++)
	{	
		m_N2 =  m_aBeamElemProValue[loop].m_nElementRC;		
		for(int loop1=0;loop1<m_N1;loop1++)
		{
			if(m_N2 == apBaseRealConstant[0]->SpaceFrameRCValue[loop1].m_nRealConsNo)
			{
				//确定截面积
				m_aBeamElemProValue[loop].m_A = apBaseRealConstant[0]->SpaceFrameRCValue[loop1].m_dArea;	
				//确定Ix
				m_aBeamElemProValue[loop].m_Ix = apBaseRealConstant[0]->SpaceFrameRCValue[loop1].m_dIx;	
				//确定Iy
				m_aBeamElemProValue[loop].m_Iy = apBaseRealConstant[0]->SpaceFrameRCValue[loop1].m_dIy;
				//确定Iz
				m_aBeamElemProValue[loop].m_Iz = apBaseRealConstant[0]->SpaceFrameRCValue[loop1].m_dIz;
				//确定m_Faiy 的面积
				m_aBeamElemProValue[loop].m_Faiy = apBaseRealConstant[0]->SpaceFrameRCValue[loop1].m_dFaiy;
				//确定m_Faiz 的面积
				m_aBeamElemProValue[loop].m_Faiz = apBaseRealConstant[0]->SpaceFrameRCValue[loop1].m_dFaiz;
			}
		}
	}

	//确定单元长度
	for(loop=0;loop<m_nElementNum;loop++)
	{
		m_N1 = m_pNode->m_nNodeNum;
		for(int loop1=0;loop1<m_N1;loop1++)
		{
			if(m_aBeamElemProValue[loop].m_nElemNode1 == m_pNode->m_CaNodeCoor[loop1].m_aiNode)
			{
				m_aBeamElemProValue[loop].m_dX1 = m_pNode->m_CaNodeCoor[loop1].m_adX;
				m_aBeamElemProValue[loop].m_dY1 = m_pNode->m_CaNodeCoor[loop1].m_adY;
				m_aBeamElemProValue[loop].m_dZ1 = m_pNode->m_CaNodeCoor[loop1].m_adZ;
			}
			if(m_aBeamElemProValue[loop].m_nElemNode2 == m_pNode->m_CaNodeCoor[loop1].m_aiNode)
			{
				m_aBeamElemProValue[loop].m_dX2 = m_pNode->m_CaNodeCoor[loop1].m_adX;
				m_aBeamElemProValue[loop].m_dY2 = m_pNode->m_CaNodeCoor[loop1].m_adY;
				m_aBeamElemProValue[loop].m_dZ2 = m_pNode->m_CaNodeCoor[loop1].m_adZ;
			}	
			if(m_aBeamElemProValue[loop].m_nElemNode3 == m_pNode->m_CaNodeCoor[loop1].m_aiNode)
			{
				m_aBeamElemProValue[loop].m_dX3 = m_pNode->m_CaNodeCoor[loop1].m_adX;
				m_aBeamElemProValue[loop].m_dY3 = m_pNode->m_CaNodeCoor[loop1].m_adY;
				m_aBeamElemProValue[loop].m_dZ3 = m_pNode->m_CaNodeCoor[loop1].m_adZ;
			}	
		}
		m_aBeamElemProValue[loop].m_L = sqrt( pow((m_aBeamElemProValue[loop].m_dX1 - m_aBeamElemProValue[loop].m_dX2),2) 
			                             + pow((m_aBeamElemProValue[loop].m_dY1 - m_aBeamElemProValue[loop].m_dY2),2) 
			                             + pow((m_aBeamElemProValue[loop].m_dZ1 - m_aBeamElemProValue[loop].m_dZ2),2) );
	}


	//确定剪切影响系数
	for(loop=0;loop<m_nElementNum;loop++)
	{	
		double EE = m_aBeamElemProValue[loop].m_E;
		double GG = m_aBeamElemProValue[loop].m_G;
		double IZZ = m_aBeamElemProValue[loop].m_Iz;
		double IYY = m_aBeamElemProValue[loop].m_Iy;		
		double LL = m_aBeamElemProValue[loop].m_L;
		double Ay = m_aBeamElemProValue[loop].m_Faiy;
		double Az = m_aBeamElemProValue[loop].m_Faiz;
		double AA = m_aBeamElemProValue[loop].m_A;		

		if(Ay>0)
			m_aBeamElemProValue[loop].m_Faiy = 12*EE*IZZ/(GG*Ay*LL*LL);
		else
			m_aBeamElemProValue[loop].m_Faiy = 0;

		if(Az>0)
			m_aBeamElemProValue[loop].m_Faiz = 12*EE*IYY/(GG*Az*LL*LL);	
		else
			m_aBeamElemProValue[loop].m_Faiz = 0;
	}
	
}


void CSpaceFrameBeam::OutputParameter()
{	
	ofstream fout;
	//********************************************	
	fout.open("CSpaceFrameBeam Class Test.tst",ios::out);
	fout<<"单元的数量   "<<m_nElementNum<<endl;

	fout.width(12);
	fout<<"单元号";
	fout.width(12);
	fout<<"节点一";
	fout.width(12);
	fout<<"节点二";
	fout.width(12);
	fout<<"节点三";
	fout.width(12);
	fout<<"材料号";
	fout.width(12);
	fout<<"实常数号"<<endl;
	for(int loop=0;loop<m_nElementNum;loop++)
	{
		fout.width(12);
		fout<<m_aBeamElemProValue[loop].m_nElementNo;
		fout.width(12);
		fout<<m_aBeamElemProValue[loop].m_nElemNode1;
		fout.width(12);
		fout<<m_aBeamElemProValue[loop].m_nElemNode2;
		fout.width(12);
		fout<<m_aBeamElemProValue[loop].m_nElemNode3;

		fout.width(12);
		fout<<m_aBeamElemProValue[loop].m_nElementMat;
		fout.width(12);
		fout<<m_aBeamElemProValue[loop].m_nElementRC<<endl;
	}
	fout<<endl;
	fout<<endl;

	fout.width(14);
	fout<<"更详细的数据："<<endl;
	fout.width(8);
	fout<<"单元号";
	fout.width(14);
	fout<<"弹性模量";
	fout.width(14);
	fout<<"剪切模量";
	fout.width(14);
	fout<<"长度";
	fout.width(14);
	fout<<"截面积";
	fout.width(14);
	fout<<"Ix";
	fout.width(14);
	fout<<"Iy";
	fout.width(14);
	fout<<"Iz";
	fout.width(14);
	fout<<"FaiY";
	fout.width(14);
	fout<<"FaiZ"<<endl;
	for(loop=0;loop<m_nElementNum;loop++)
	{
		fout.width(8);
		fout<<m_aBeamElemProValue[loop].m_nElementNo;
		fout.width(14);
		fout<<m_aBeamElemProValue[loop].m_E;
		fout.width(14);
		fout<<m_aBeamElemProValue[loop].m_G;
		fout.width(14);
		fout<<m_aBeamElemProValue[loop].m_L;
		fout.width(14);
		fout<<m_aBeamElemProValue[loop].m_A;
		fout.width(14);
		fout<<m_aBeamElemProValue[loop].m_Ix;
		fout.width(14);
		fout<<m_aBeamElemProValue[loop].m_Iy;
		fout.width(14);
		fout<<m_aBeamElemProValue[loop].m_Iz;
		fout.width(14);
		fout<<m_aBeamElemProValue[loop].m_Faiy;
		fout.width(14);
		fout<<m_aBeamElemProValue[loop].m_Faiz<<endl;
	}
	fout<<endl;

	fout.width(15);
	fout<<"单元号";
	fout.width(15);
	fout<<"节点1";
	fout.width(15);
	fout<<"X座标";
	fout.width(15);
	fout<<"Y座标";
	fout.width(15);
	fout<<"Z座标"<<endl;
	for(loop=0;loop<m_nElementNum;loop++)
	{
		fout.width(15);
		fout<<m_aBeamElemProValue[loop].m_nElementNo;
		fout.width(15);
		fout<<m_aBeamElemProValue[loop].m_nElemNode1;
		fout.width(15);
		fout<<m_aBeamElemProValue[loop].m_dX1;
		fout.width(15);
		fout<<m_aBeamElemProValue[loop].m_dY1;
		fout.width(15);
		fout<<m_aBeamElemProValue[loop].m_dZ1<<endl;
	}	

	fout.width(15);
	fout<<"单元号";
	fout.width(15);
	fout<<"节点2";
	fout.width(15);
	fout<<"X座标";
	fout.width(15);
	fout<<"Y座标";
	fout.width(15);
	fout<<"Z座标"<<endl;
	for(loop=0;loop<m_nElementNum;loop++)
	{
		fout.width(15);
		fout<<m_aBeamElemProValue[loop].m_nElementNo;
		fout.width(15);
		fout<<m_aBeamElemProValue[loop].m_nElemNode2;
		fout.width(15);
		fout<<m_aBeamElemProValue[loop].m_dX2;
		fout.width(15);
		fout<<m_aBeamElemProValue[loop].m_dY2;
		fout.width(15);
		fout<<m_aBeamElemProValue[loop].m_dZ2<<endl;
	}

	fout.width(15);
	fout<<"单元号";
	fout.width(15);
	fout<<"节点3";
	fout.width(15);
	fout<<"X座标";
	fout.width(15);
	fout<<"Y座标";
	fout.width(15);
	fout<<"Z座标"<<endl;
	for(loop=0;loop<m_nElementNum;loop++)
	{
		fout.width(15);
		fout<<m_aBeamElemProValue[loop].m_nElementNo;
		fout.width(15);
		fout<<m_aBeamElemProValue[loop].m_nElemNode3;
		fout.width(15);
		fout<<m_aBeamElemProValue[loop].m_dX3;
		fout.width(15);
		fout<<m_aBeamElemProValue[loop].m_dY3;
		fout.width(15);
		fout<<m_aBeamElemProValue[loop].m_dZ3<<endl;
	}	
	fout<<endl;
	
	fout<<"单元节点释放情况："<<endl;
	fout.width(10);
	fout<<"单元号：";
	fout.width(10);
	fout<<"I_UX";
	fout.width(10);
	fout<<"I_UY";
	fout.width(10);
	fout<<"I_UZ";
	fout.width(10);
	fout<<"I_ROTX";
	fout.width(10);
	fout<<"I_ROTY";
	fout.width(10);
	fout<<"I_ROTZ";
	fout.width(10);
	fout<<"J_UX";
	fout.width(10);
	fout<<"J_UY";
	fout.width(10);
	fout<<"J_UZ";
	fout.width(10);
	fout<<"J_ROTX";
	fout.width(10);
	fout<<"J_ROTY";
	fout.width(10);
	fout<<"J_ROTZ"<<endl;
	for(loop=0;loop<m_nElementNum;loop++)
	{
		fout.width(10);		
		fout<<m_aBeamElemProValue[loop].m_nElementNo;
		fout.width(10);
		fout<<m_aBeamElemProValue[loop].m_nEleRelease_I1;
		fout.width(10);
		fout<<m_aBeamElemProValue[loop].m_nEleRelease_I2;
		fout.width(10);
		fout<<m_aBeamElemProValue[loop].m_nEleRelease_I3;
		fout.width(10);
		fout<<m_aBeamElemProValue[loop].m_nEleRelease_I4;
		fout.width(10);
		fout<<m_aBeamElemProValue[loop].m_nEleRelease_I5;
		fout.width(10);
		fout<<m_aBeamElemProValue[loop].m_nEleRelease_I6;
		fout.width(10);
		fout<<m_aBeamElemProValue[loop].m_nEleRelease_J1;
		fout.width(10);
		fout<<m_aBeamElemProValue[loop].m_nEleRelease_J2;
		fout.width(10);
		fout<<m_aBeamElemProValue[loop].m_nEleRelease_J3;
		fout.width(10);
		fout<<m_aBeamElemProValue[loop].m_nEleRelease_J4;
		fout.width(10);
		fout<<m_aBeamElemProValue[loop].m_nEleRelease_J5;
		fout.width(10);
		fout<<m_aBeamElemProValue[loop].m_nEleRelease_J6<<endl;
	}
	fout<<endl;
	fout<<endl;

	fout<<"矩阵计算数据  "<<endl;
	for(loop=0;loop<m_nElementNum;loop++)
	{
		for(int loop1=0;loop1<12;loop1++)
		{
			for(int loop2=0;loop2<12;loop2++)
			{
				fout.width(10);
				fout<<m_CaElemOrigMatrix[loop](loop1,loop2);
			}	
			fout<<endl;
		}		
		fout<<endl;	
		fout<<endl;	
	}

	fout<<endl;
	fout<<"单元座标系－>整体座标系 转换矩阵  "<<endl;
	for(loop=0;loop<m_nElementNum;loop++)
	{
        for(int loop1=0;loop1<12;loop1++)
		{
			for(int loop2=0;loop2<12;loop2++)
			{
				fout.width(10);
				fout<<m_CaElemTranFactorMatrix[loop](loop1,loop2);
			}	
			fout<<endl;
		}
		fout<<endl;	
		fout<<endl;	
	}

	fout<<endl;
	fout<<"单元座标系－>整体座标系 转换矩阵 转置  "<<endl;