📄 spaceframebeam.cpp

📁 三维框架梁单元C++版本的源程序可以使结构力学的概念更上升一个层次
💻 CPP
📖 第 1 页 / 共 2 页
字号:
上一页 12
	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_CaElemTranFactorMatrixT[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_CaElemResultMatrix[loop](loop1,loop2);
			}	
			fout<<endl;
		}
		fout<<endl;	
		fout<<endl;	
	}
	fout.close();	
}

void CSpaceFrameBeam::FormElemMatrix()
{
	//调整好单元刚度矩阵
   m_CaElemOrigMatrix.SetSize(m_nElementNum);

	for(int loop=0;loop<m_nElementNum;loop++)
	{	
		//初始化,开设12x12的矩阵空间
		m_CaElemOrigMatrix[loop](12,12);

		//初始化,所有初值为0
		for(int loop1=0;loop1<12;loop1++)
		{
			for(int loop2=0;loop2<12;loop2++)
				m_CaElemOrigMatrix[loop](loop1,loop2) = 0;			
		}

		//再赋值，便于语句书写和校对！
		double E = m_aBeamElemProValue[loop].m_E;
		double G = m_aBeamElemProValue[loop].m_G;
		double A = m_aBeamElemProValue[loop].m_A;
		double L = m_aBeamElemProValue[loop].m_L;
		double IX = m_aBeamElemProValue[loop].m_Ix;
		double IY = m_aBeamElemProValue[loop].m_Iy;
		double IZ = m_aBeamElemProValue[loop].m_Iz;
		double FaiY = m_aBeamElemProValue[loop].m_Faiy;
		double FaiZ = m_aBeamElemProValue[loop].m_Faiz;
		//*
		m_CaElemOrigMatrix[loop](0,0) = E*A/L;
		m_CaElemOrigMatrix[loop](1,1) = 12*E*IZ/(L*L*L*(1+FaiY));
		m_CaElemOrigMatrix[loop](2,2) = 12*E*IY/(L*L*L*(1+FaiZ));
		m_CaElemOrigMatrix[loop](3,3) = G*IX/L;
		m_CaElemOrigMatrix[loop](4,4) = (4+FaiZ)*E*IY/((1+FaiZ)*L);
		//*	
		m_CaElemOrigMatrix[loop](4,2) = -6*E*IY/(L*L*(1+FaiZ));
		m_CaElemOrigMatrix[loop](2,4) = -6*E*IY/(L*L*(1+FaiZ));
		//*	
		m_CaElemOrigMatrix[loop](5,5) = (4+FaiY)*E*IZ/((1+FaiY)*L);
		//*	
		m_CaElemOrigMatrix[loop](5,1) = 6*E*IZ/(L*L*(1+FaiY));
		m_CaElemOrigMatrix[loop](1,5) = 6*E*IZ/(L*L*(1+FaiY));
		//*	
		m_CaElemOrigMatrix[loop](6,6) = E*A/L;
		//*	
		m_CaElemOrigMatrix[loop](6,0) = -E*A/L;
		m_CaElemOrigMatrix[loop](0,6) = -E*A/L;
		//*	
		m_CaElemOrigMatrix[loop](7,7) = 12*E*IZ/(L*L*L*(1+FaiY));
		//*	
		m_CaElemOrigMatrix[loop](7,1) = -12*E*IZ/(L*L*L*(1+FaiY));
		m_CaElemOrigMatrix[loop](1,7) = -12*E*IZ/(L*L*L*(1+FaiY));
		m_CaElemOrigMatrix[loop](7,5) = -6*E*IZ/(L*L*(1+FaiY));
		m_CaElemOrigMatrix[loop](5,7) = -6*E*IZ/(L*L*(1+FaiY));
		//*	
		m_CaElemOrigMatrix[loop](8,8) = 12*E*IY/(L*L*L*(1+FaiZ));
		//*	
		m_CaElemOrigMatrix[loop](8,2) = -12*E*IY/(L*L*L*(1+FaiZ));
		m_CaElemOrigMatrix[loop](2,8) = -12*E*IY/(L*L*L*(1+FaiZ));
		m_CaElemOrigMatrix[loop](8,4) = 6*E*IY/(L*L*(1+FaiZ));
		m_CaElemOrigMatrix[loop](4,8) = 6*E*IY/(L*L*(1+FaiZ));
		//*	
		m_CaElemOrigMatrix[loop](9,9) = G*IX/L;
		//*	
		m_CaElemOrigMatrix[loop](9,3) = -G*IX/L;
		m_CaElemOrigMatrix[loop](3,9) = -G*IX/L;
		//*	
		m_CaElemOrigMatrix[loop](10,10) = (4+FaiZ)*E*IY/((1+FaiZ)*L);
		//*	
		m_CaElemOrigMatrix[loop](10,2) = -6*E*IY/(L*L*(1+FaiZ));
		m_CaElemOrigMatrix[loop](2,10) = -6*E*IY/(L*L*(1+FaiZ));
		m_CaElemOrigMatrix[loop](10,4) = (2-FaiZ)*E*IY/(L*(1+FaiZ));
		m_CaElemOrigMatrix[loop](4,10) = (2-FaiZ)*E*IY/(L*(1+FaiZ));
		m_CaElemOrigMatrix[loop](10,8) = 6*E*IY/(L*L*(1+FaiZ));
		m_CaElemOrigMatrix[loop](8,10) = 6*E*IY/(L*L*(1+FaiZ));
		//*	
		m_CaElemOrigMatrix[loop](11,11) = (4+FaiY)*E*IZ/((1+FaiY)*L);
		//*	
		m_CaElemOrigMatrix[loop](11,1) = 6*E*IZ/(L*L*(1+FaiY));
		m_CaElemOrigMatrix[loop](1,11) = 6*E*IZ/(L*L*(1+FaiY));
		m_CaElemOrigMatrix[loop](11,5) = (2-FaiY)*E*IZ/((1+FaiY)*L);
		m_CaElemOrigMatrix[loop](5,11) = (2-FaiY)*E*IZ/((1+FaiY)*L);
		m_CaElemOrigMatrix[loop](11,7) = -6*E*IZ/(L*L*(1+FaiY));
		m_CaElemOrigMatrix[loop](7,11) = -6*E*IZ/(L*L*(1+FaiY));
	}

	//还不能考虑节点释放
	for(loop=0;loop<m_nElementNum;loop++)
	{}
}

void CSpaceFrameBeam::TranElemMatrix()
{
	//******************************定义转换矩阵******************************
	//设定大小
	m_CaElemTranFactorMatrix.SetSize(m_nElementNum);
	for(int loop=0;loop<m_nElementNum;loop++)
	{
		//初始化,开设12x12的矩阵空间
		m_CaElemTranFactorMatrix[loop](12,12);

		//初始化,所有初值为0
		for(int loop1=0;loop1<12;loop1++)
		{
			for(int loop2=0;loop2<12;loop2++)
				m_CaElemTranFactorMatrix[loop](loop1,loop2) = 0;			
		}

		//**********************************************************************		
		double a1 = m_aBeamElemProValue[loop].m_dX2 - m_aBeamElemProValue[loop].m_dX1;
		double a2 = m_aBeamElemProValue[loop].m_dY2 - m_aBeamElemProValue[loop].m_dY1;
		double a3 = m_aBeamElemProValue[loop].m_dZ2 - m_aBeamElemProValue[loop].m_dZ1;
		double t11 = a1/m_aBeamElemProValue[loop].m_L;
		double t12 = a2/m_aBeamElemProValue[loop].m_L;
		double t13 = a3/m_aBeamElemProValue[loop].m_L;

		double b1 = m_aBeamElemProValue[loop].m_dX3 - m_aBeamElemProValue[loop].m_dX1;
		double b2 = m_aBeamElemProValue[loop].m_dY3 - m_aBeamElemProValue[loop].m_dY1;
		double b3 = m_aBeamElemProValue[loop].m_dZ3 - m_aBeamElemProValue[loop].m_dZ1; 
		double s1 = a2*b3 - a3*b2;
		double s2 = a3*b1 - a1*b3;
		double s3 = a1*b2 - a2*b1;
		double u1 = s2*a3 - s3*a2;
		double u2 = s3*a1 - s1*a3;
		double u3 = s1*a2 - s2*a1;
		double u = sqrt(pow(u1,2) + pow(u2,2) + pow(u3,2));		

		double t21 = u1/u;
		double t22 = u2/u;
		double t23 = u3/u;

		double t31 = t12*t23 - t13*t22;
		double t32 = t13*t21 - t11*t23;
		double t33 = t11*t22 - t12*t21;
		//**********************************************************************

		//*****
		m_CaElemTranFactorMatrix[loop](0,0) = t11;
		m_CaElemTranFactorMatrix[loop](0,1) = t12;
		m_CaElemTranFactorMatrix[loop](0,2) = t13;
		//*
		m_CaElemTranFactorMatrix[loop](1,0) = t21;
		m_CaElemTranFactorMatrix[loop](1,1) = t22;
		m_CaElemTranFactorMatrix[loop](1,2) = t23;
		//*
		m_CaElemTranFactorMatrix[loop](2,0) = t31;
		m_CaElemTranFactorMatrix[loop](2,1) = t32;
		m_CaElemTranFactorMatrix[loop](2,2) = t33;

		//*****
		m_CaElemTranFactorMatrix[loop](3,3) = t11;
		m_CaElemTranFactorMatrix[loop](3,4) = t12;
		m_CaElemTranFactorMatrix[loop](3,5) = t13;
		//*
		m_CaElemTranFactorMatrix[loop](4,3) = t21;
		m_CaElemTranFactorMatrix[loop](4,4) = t22;
		m_CaElemTranFactorMatrix[loop](4,5) = t23;
		//*
		m_CaElemTranFactorMatrix[loop](5,3) = t31;
		m_CaElemTranFactorMatrix[loop](5,4) = t32;
		m_CaElemTranFactorMatrix[loop](5,5) = t33;

		//*****
		m_CaElemTranFactorMatrix[loop](6,6) = t11;
		m_CaElemTranFactorMatrix[loop](6,7) = t12;
		m_CaElemTranFactorMatrix[loop](6,8) = t13;
		//*
		m_CaElemTranFactorMatrix[loop](7,6) = t21;
		m_CaElemTranFactorMatrix[loop](7,7) = t22;
		m_CaElemTranFactorMatrix[loop](7,8) = t23;
		//*
		m_CaElemTranFactorMatrix[loop](8,6) = t31;
		m_CaElemTranFactorMatrix[loop](8,7) = t32;
		m_CaElemTranFactorMatrix[loop](8,8) = t33;

		//*****
		m_CaElemTranFactorMatrix[loop](9,9) = t11;
		m_CaElemTranFactorMatrix[loop](9,10) = t12;
		m_CaElemTranFactorMatrix[loop](9,11) = t13;
		//*
		m_CaElemTranFactorMatrix[loop](10,9) = t21;
		m_CaElemTranFactorMatrix[loop](10,10) = t22;
		m_CaElemTranFactorMatrix[loop](10,11) = t23;
		//*
		m_CaElemTranFactorMatrix[loop](11,9) = t31;
		m_CaElemTranFactorMatrix[loop](11,10) = t32;
		m_CaElemTranFactorMatrix[loop](11,11) = t33;		
	}

	//******************************求转换矩阵的转置******************************
	m_CaElemTranFactorMatrixT.SetSize(m_nElementNum);	
	for(loop=0;loop<m_nElementNum;loop++)
	{
		m_CaElemTranFactorMatrixT[loop].Trans(m_CaElemTranFactorMatrix[loop]);
	}

	
	//******************************求整体坐标系下的单元刚度矩阵******************************
	m_CaElemResultMatrix.SetSize(m_nElementNum);	
	
	for(loop=0;loop<m_nElementNum;loop++)
	{
		//临时矩阵
		CMatrix TempMat(12,12,0);
		//矩阵类似乎只有初始化以后才能试用
		m_CaElemResultMatrix[loop] = TempMat;

		for(int loop1=0;loop1<12;loop1++)
		{
			for(int loop2=0;loop2<12;loop2++)
			{
				for(int loop3=0;loop3<12;loop3++)
				{
					TempMat(loop1,loop2) +=
						m_CaElemTranFactorMatrixT[loop](loop1,loop3) *
						m_CaElemOrigMatrix[loop](loop3,loop2);
				}
			}
		}
		for(loop1=0;loop1<12;loop1++)
		{
			for(int loop2=0;loop2<12;loop2++)
			{
				for(int loop3=0;loop3<12;loop3++)
				{
					m_CaElemResultMatrix[loop](loop1,loop2) +=
						TempMat(loop1,loop3) * m_CaElemTranFactorMatrix[loop](loop3,loop2);
				}
			}	
		}    
	}

}

void CSpaceFrameBeam::CalcuGKBandWidth()
{
	//带宽数组初始化
	m_aGKBandWith.SetSize(m_pNode->m_nTotalDOF);	
	for(int ii=0;ii<m_pNode->m_nTotalDOF;ii++)	
		m_aGKBandWith[ii] = 1; //带宽至少为1

	//临时存储空间，存储一个单元、两个节点的12个自由度的序列编号
	CUIntArray aiDOFIndex;
	aiDOFIndex.SetSize(12);	

	//按单元循环
	for(int loop=0;loop<m_nElementNum;loop++)
	{
		//结点数量
		int NodeNum = m_pNode->m_nNodeNum - m_pNode->m_nNode_K_Num;
		//节点1、2
		int Node1 = m_aBeamElemProValue[loop].m_nElemNode1;
		int Node2 = m_aBeamElemProValue[loop].m_nElemNode2;
		//匹配结点数据		
		for(int loop1=0;loop1<NodeNum;loop1++)
		{			
			if(Node1 == m_pNode->m_aDOFIndex[loop1].m_aiNode)
			{
				aiDOFIndex[0] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_X;
				aiDOFIndex[1] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_Y;
				aiDOFIndex[2] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_Z;
				aiDOFIndex[3] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_RX;
				aiDOFIndex[4] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_RY;
				aiDOFIndex[5] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_RZ;
			}	
			if(Node2 == m_pNode->m_aDOFIndex[loop1].m_aiNode)
			{
				aiDOFIndex[6] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_X;
				aiDOFIndex[7] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_Y;
				aiDOFIndex[8] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_Z;
				aiDOFIndex[9] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_RX;
				aiDOFIndex[10] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_RY;
				aiDOFIndex[11] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_RZ;
			}	
		}		

		//求最小自由度编号
		int iMinDOFIndex = m_pNode->m_nTotalDOF;
		for(int loop2=0;loop2<12;loop2++)
		{
			if(iMinDOFIndex > aiDOFIndex[loop2])
				iMinDOFIndex = aiDOFIndex[loop2];			
		}		

		//求矩阵带宽
		for(int loop3=0;loop3<12;loop3++)
		{
			int iBuf = aiDOFIndex[loop3] - iMinDOFIndex + 1;			
			if(m_aGKBandWith[aiDOFIndex[loop3]] < iBuf)
				m_aGKBandWith[aiDOFIndex[loop3]] = iBuf;			
		}		
	}

	//********************************************
	ofstream fout;
	fout.open("测试带宽数组.tst");
	for(int loop5=0;loop5<m_pNode->m_nTotalDOF;loop5++)
		fout<<m_aGKBandWith[loop5]<<endl;
	fout.close();   	
}


void CSpaceFrameBeam::CalcuGKDiagAddress()
{
	//对角地址实例化
	aiGKDiagAdd = new unsigned long [m_pNode->m_nTotalDOF]; 
	//初始化	
	for(int loop=0;loop<m_pNode->m_nTotalDOF;loop++)
		aiGKDiagAdd[loop] = 0;
	//具体赋值
	for(loop=1;loop<m_pNode->m_nTotalDOF;loop++)
	{
		int iBuf = 0;		
		for(int loop1=0;loop1<=loop;loop1++)
		{			
			iBuf += m_aGKBandWith[loop1];
		}
		aiGKDiagAdd[loop] = iBuf - 1; 
	} 

	 //********************************************	
	ofstream fout;
	fout.open("测试对角地址数组.tst");
	for(loop=0;loop<m_pNode->m_nTotalDOF;loop++)
		fout<<aiGKDiagAdd[loop]<<endl;
	fout.close();
}

void CSpaceFrameBeam::SparseMatrixStiffAssemble()
{
	//初始化
	m_smatGK.Realloc(m_pNode->m_nTotalDOF,aiGKDiagAdd);		
	delete aiGKDiagAdd;
	//具体赋值
	m_smatGK.ElementBufRealloc();
	m_smatGK=0.0;

	//储存自由度编号的临时数组
	CUIntArray DOFIndex;  
	DOFIndex.SetSize(12); 	

    for(int loop=0;loop<m_nElementNum;loop++)
	{
		//节点1、2
		int Node1 = m_aBeamElemProValue[loop].m_nElemNode1;
		int Node2 = m_aBeamElemProValue[loop].m_nElemNode2;

		//临时数组初始化
		for(int ii=0;ii<12;ii++)
			DOFIndex[ii] = 0;	

		//临时数组赋值
		for(int loop1=0;loop1<(m_pNode->m_nNodeNum - m_pNode->m_nNode_K_Num);loop1++)
		{
			if(Node1 == m_pNode->m_aDOFIndex[loop1].m_aiNode)
			{
				DOFIndex[0] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_X;
				DOFIndex[1] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_Y;
				DOFIndex[2] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_Z;	
				DOFIndex[3] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_RX;	
				DOFIndex[4] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_RY;	
				DOFIndex[5] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_RZ;				
			}
			if(Node2 == m_pNode->m_aDOFIndex[loop1].m_aiNode)
			{	
				DOFIndex[6] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_X;	
				DOFIndex[7] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_Y;
				DOFIndex[8] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_Z;
				DOFIndex[9] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_RX;	
				DOFIndex[10] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_RY;	
				DOFIndex[11] = m_pNode->m_aDOFIndex[loop1].m_aiDOFIndex_RZ;				
			}
		}		

		for(int loop2=0;loop2<12;loop2++)
		{
			int iBuf = DOFIndex[loop2];
			//对角元素加入总刚矩阵
			m_smatGK(iBuf,iBuf) += m_CaElemResultMatrix[loop](loop2,loop2);
				
			//下三角非对称元素加入总刚矩阵
			for(int loop3=0;loop3<loop2;loop3++)			
				m_smatGK(iBuf,DOFIndex[loop3]) += m_CaElemResultMatrix[loop](loop2,loop3);						
			
		}
	}

	//********************************************	
	ofstream fout;
	fout.open("Sparse Matrix.tst");
	for(int ii=0;ii<m_pNode->m_nTotalDOF;ii++)
	{
		for(int jj=0;jj<m_pNode->m_nTotalDOF;jj++)
		{
			fout.width(10);
			fout<<m_smatGK(ii,jj)<<"  ";
		}
		fout<<endl;
	}
	fout.close();	
}

void CSpaceFrameBeam::InternalForceInitial()
{
	//初始化单元内力
	m_aElemInForceDef.SetSize(m_nElementNum);
	for(int loop=0;loop<m_nElementNum;loop++)
	{
		m_aElemInForceDef[loop].m_nElemNo = m_aBeamElemProValue[loop].m_nElementNo;	
		for(int loop1=0;loop1<13;loop1++)
		{
			m_aElemInForceDef[loop].m_dElemFX[loop1] = 0;
			m_aElemInForceDef[loop].m_dElemShearY[loop1] = 0;
			m_aElemInForceDef[loop].m_dElemShearZ[loop1] = 0;
			m_aElemInForceDef[loop].m_dElemT[loop1] = 0;
			m_aElemInForceDef[loop].m_dElemMY[loop1] = 0;
			m_aElemInForceDef[loop].m_dElemMZ[loop1] = 0;			
		}
	}

	m_aElemInForceDef_N.SetSize(m_nElementNum);
	for(loop=0;loop<m_nElementNum;loop++)
	{
		m_aElemInForceDef_N[loop].m_nElemNo = m_aBeamElemProValue[loop].m_nElementNo;	
		for(int loop1=0;loop1<13;loop1++)
		{
			m_aElemInForceDef_N[loop].m_dElemFX[loop1] = 0;
			m_aElemInForceDef_N[loop].m_dElemShearY[loop1] = 0;
			m_aElemInForceDef_N[loop].m_dElemShearZ[loop1] = 0;
			m_aElemInForceDef_N[loop].m_dElemT[loop1] = 0;
			m_aElemInForceDef_N[loop].m_dElemMY[loop1] = 0;
			m_aElemInForceDef_N[loop].m_dElemMZ[loop1] = 0;	
		}
	}

	//初始化单元挠度
	m_aElemDeflection.SetSize(m_nElementNum);
	m_aElemDeflection_N.SetSize(m_nElementNum);
	for(loop=0;loop<m_nElementNum;loop++)
	{
		m_aElemDeflection[loop].m_nElemNo = m_aBeamElemProValue[loop].m_nElementNo;	
		m_aElemDeflection_N[loop].m_nElemNo = m_aBeamElemProValue[loop].m_nElementNo;	
		for(int loop1=0;loop1<13;loop1++)
		{
			m_aElemDeflection[loop].m_dDelf_X[loop1] = 0;
			m_aElemDeflection[loop].m_dDelf_Y[loop1] = 0;
			m_aElemDeflection[loop].m_dDelf_Z[loop1] = 0;
			//
			m_aElemDeflection_N[loop].m_dDelf_X[loop1] = 0;
			m_aElemDeflection_N[loop].m_dDelf_Y[loop1] = 0;
			m_aElemDeflection_N[loop].m_dDelf_Z[loop1] = 0;
		}
	}

	//初始化支座反力
	m_nSupportReaction.SetSize(m_pNode->m_nNodeFixNum);
	for(loop=0;loop<m_pNode->m_nNodeFixNum;loop++)
	{
		m_nSupportReaction[loop].m_nNodeNo = m_pNode->m_CaNodeFix[loop].m_aiConstrainedNode;
		m_nSupportReaction[loop].m_dSR_X = 0;
		m_nSupportReaction[loop].m_dSR_Y = 0;
		m_nSupportReaction[loop].m_dSR_Z = 0;	
		m_nSupportReaction[loop].m_dSR_RX = 0;
		m_nSupportReaction[loop].m_dSR_RY = 0;
		m_nSupportReaction[loop].m_dSR_RZ = 0;	
	}
}
上一页 12
⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -