globalelement.cpp

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

				for(int loop3=0;loop3<6;loop3++)
				{
					x = loop3*L/12;
					m_apEle[0]->m_aElemDeflection[loop1].m_dDelf_Y[loop3] += 
						m_dDY*L*L*x*x*(x/L)*(x/L)/(24*E*IZ);
					m_apEle[0]->m_aElemDeflection[loop1].m_dDelf_Z[loop3] += 
						m_dDZ*L*L*x*x*(x/L)*(x/L)/(24*E*IY);
				}
				for(loop3=6;loop3<13;loop3++)
				{
					x = loop3*L/12;
					m_apEle[0]->m_aElemDeflection[loop1].m_dDelf_Y[loop3] += 
						m_dDY*L*L*x*x*(1-(x/L))*(1-(x/L))/(24*E*IZ);
					m_apEle[0]->m_aElemDeflection[loop1].m_dDelf_Z[loop3] += 
						m_dDZ*L*L*x*x*(1-(x/L))*(1-(x/L))/(24*E*IY);
				}
			}			
		}		
	}	


	//CASE2 初始化
	m_adCase2_ElemDis_LoadVector.SetSize(m_Node.m_nTotalDOF);
	for(loop=0;loop<m_Node.m_nTotalDOF;loop++)
		m_adCase2_ElemDis_LoadVector[loop] = 0;
	//临时荷载向量（局部座标->整体座标）
	for(loop=0;loop<12;loop++)
	{
		TempCM_DisLoadVector(loop,0) = 0;	
		TempCM_DisLoadVectorG(loop,0) = 0;
	}

	//具体赋值
	for(loop=0;loop<m_Load.m_nElemLoadNum_Case2_Dis;loop++)
	{
		//受荷单元号
		m_nTemp = m_Load.m_adElemLoadValue_Case2_Dis[loop].m_nElementNO;
		//单元对应节点
		for(int loop1=0;loop1<m_apEle[0]->m_nElementNum;loop1++)
		{
			if(m_nTemp == m_apEle[0]->m_aBeamElemProValue[loop1].m_nElementNo)
			{
				m_nNode1 = m_apEle[0]->m_aBeamElemProValue[loop1].m_nElemNode1;
				m_nNode2 = m_apEle[0]->m_aBeamElemProValue[loop1].m_nElemNode2;
				L = m_apEle[0]->m_aBeamElemProValue[loop1].m_L;
				//转化矩阵
				CMatrix cm_TT = m_apEle[0]->m_CaElemTranFactorMatrixT[loop1];				
				//生成荷载向量
				m_dDX = m_Load.m_adElemLoadValue_Case2_Dis[loop].m_dElemForceX * L/2;
				m_dDY = m_Load.m_adElemLoadValue_Case2_Dis[loop].m_dElemForceY * L/2;
				m_dDZ = m_Load.m_adElemLoadValue_Case2_Dis[loop].m_dElemForceZ * L/2;
				m_dDMX = 0;
				m_dDMY = m_Load.m_adElemLoadValue_Case2_Dis[loop].m_dElemForceZ *L*L/12;
				m_dDMZ = m_Load.m_adElemLoadValue_Case2_Dis[loop].m_dElemForceY *L*L/12;
				//局部荷载向量
				TempCM_DisLoadVector(0,0) = m_dDX;
				TempCM_DisLoadVector(1,0) = m_dDY;
				TempCM_DisLoadVector(2,0) = m_dDZ;
				TempCM_DisLoadVector(3,0) = m_dDMX;
				TempCM_DisLoadVector(4,0) = -m_dDMY;
				TempCM_DisLoadVector(5,0) = m_dDMZ;
				//
				TempCM_DisLoadVector(6,0) = m_dDX;
				TempCM_DisLoadVector(7,0) = m_dDY;
				TempCM_DisLoadVector(8,0) = m_dDZ;
				TempCM_DisLoadVector(9,0) = -m_dDMX;
				TempCM_DisLoadVector(10,0) = m_dDMY;
				TempCM_DisLoadVector(11,0) = -m_dDMZ;
				//整体荷载向量
				TempCM_DisLoadVectorG = cm_TT * TempCM_DisLoadVector;

				for(int loop2=0;loop2<(m_Node.m_nNodeNum - m_Node.m_nNode_K_Num);loop2++)
				{
					if(m_nNode1 == m_Node.m_aDOFIndex[loop2].m_aiNode)
					{
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_X;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase2_ElemDis_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(0,0);
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_Y;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase2_ElemDis_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(1,0);
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_Z;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase2_ElemDis_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(2,0);	
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_RX;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase2_ElemDis_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(3,0);
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_RY;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase2_ElemDis_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(4,0);
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_RZ;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase2_ElemDis_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(5,0);	
					}
					if(m_nNode2 == m_Node.m_aDOFIndex[loop2].m_aiNode)
					{
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_X;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase2_ElemDis_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(6,0);
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_Y;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase2_ElemDis_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(7,0);
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_Z;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase2_ElemDis_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(8,0);	
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_RX;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase2_ElemDis_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(9,0);
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_RY;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase2_ElemDis_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(10,0);
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_RZ;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase2_ElemDis_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(11,0);	
					}
				}

				//调整轴力内力
				for(int loop3=0;loop3<13;loop3++)
				{
					double x = L*loop3/12;
					Q = m_Load.m_adElemLoadValue_Case2_Dis[loop].m_dElemForceX;
					double RA = -Q*L/2;
					//
					m_apEle[0]->m_aElemInForceDef[loop1].m_dElemFX[loop3] += (RA + Q*(L-x));
				}				
				//调整剪力Y内力
				for(loop3=0;loop3<13;loop3++)
				{					
					double x = L*loop3/12;
					Q = m_Load.m_adElemLoadValue_Case2_Dis[loop].m_dElemForceY;
					double RA = -Q*L/2;
					//					
                    m_apEle[0]->m_aElemInForceDef[loop1].m_dElemShearY[loop3] += (RA + Q*(L-x));						
				}
				
				//调整剪力Z内力 
				for(loop3=0;loop3<13;loop3++)
				{					
					double x = L*loop3/12;
					Q = m_Load.m_adElemLoadValue_Case2_Dis[loop].m_dElemForceZ;
					double RA = -Q*L/2;
					//					
                    m_apEle[0]->m_aElemInForceDef[loop1].m_dElemShearZ[loop3] += (RA + Q*(L-x));						
				}

				//调整MX内力（不要做！）

				//调整MY内力
				for(loop3=0;loop3<13;loop3++)
				{	
					double x = L*loop3/12;
					Q = m_Load.m_adElemLoadValue_Case2_Dis[loop].m_dElemForceZ;
					MA = -Q * L*L/12;
					RA = Q * L/2;
					//
					m_apEle[0]->m_aElemInForceDef[loop1].m_dElemMY[loop3] += (MA + RA*x -Q*x*x/2);	
				}

				//调整MZ内力
				for(loop3=0;loop3<13;loop3++)
				{	
					double x = L*loop3/12;
					Q = m_Load.m_adElemLoadValue_Case2_Dis[loop].m_dElemForceY;
					MA = -Q * L*L/12;
					RA = Q * L/2;
					//
					m_apEle[0]->m_aElemInForceDef[loop1].m_dElemMZ[loop3] += -(MA + RA*x -Q*x*x/2);	
				}
			}			
		}		
	}
}

void CGlobalElement::ElemCen_LoadVectorAssemble()
{
	//临时变量
	int m_nTemp = 0;
	int m_nTemp1 = 0;
	int m_nNode1 = 0;
	int m_nNode2 = 0;
	double L = 0;
	double X = 0;
	double RA = 0;
	double RB = 0;
	double MA = 0;
	double MB = 0;

	double P = 0;
	//
	//临时荷载向量（局部座标->整体座标）
	CMatrix mat_TT;
	CMatrix TempCM_DisLoadVector(12,1,0);	
	CMatrix TempCM_DisLoadVectorG(12,1,0);	

	//初始化  CASE 1
	m_adCase1_ElemCen_LoadVector.SetSize(m_Node.m_nTotalDOF);
	for(int loop=0;loop<m_Node.m_nTotalDOF;loop++)
		m_adCase1_ElemCen_LoadVector[loop] = 0;

	//具体赋值
	for(loop=0;loop<m_Load.m_nElemLoadNum_Case1_Cen;loop++)
	{
		//受荷单元号
		m_nTemp = m_Load.m_adElemLoadValue_Case1_Cen[loop].m_nElementNO;
		//作用点偏移
		X = m_Load.m_adElemLoadValue_Case1_Cen[loop].m_dDistance;
		if(X <= 0.0)
		{
			AfxMessageBox("单元荷载输入有误",MB_OK,1);
			break;
		}
		
		//单元对应节点
		for(int loop1=0;loop1<m_apEle[0]->m_nElementNum;loop1++)
		{
			if(m_nTemp == m_apEle[0]->m_aBeamElemProValue[loop1].m_nElementNo)
			{
				m_nNode1 = m_apEle[0]->m_aBeamElemProValue[loop1].m_nElemNode1;
				m_nNode2 = m_apEle[0]->m_aBeamElemProValue[loop1].m_nElemNode2;
				L = m_apEle[0]->m_aBeamElemProValue[loop1].m_L;
				if( X > L)
				{
					AfxMessageBox("单元荷载输入有误",MB_OK,1);
					break;
				}
				//转换矩阵
				mat_TT = m_apEle[0]->m_CaElemTranFactorMatrixT[loop1];
				//生成荷载向量					
				P = m_Load.m_adElemLoadValue_Case1_Cen[loop].m_dElemForceX;
				TempCM_DisLoadVector(0,0) = P/2;
				TempCM_DisLoadVector(6,0) = P/2;
				P = m_Load.m_adElemLoadValue_Case1_Cen[loop].m_dElemForceY;
				double P1 = P;
				TempCM_DisLoadVector(1,0) = P*(L-X)*(L-X)*(L+2*X)/(L*L*L);
				TempCM_DisLoadVector(7,0) = P*X*X*(L+2*(L-X))/(L*L*L);
				P = m_Load.m_adElemLoadValue_Case1_Cen[loop].m_dElemForceZ;
				double P2 = P;
				TempCM_DisLoadVector(2,0) = P*(L-X)*(L-X)*(L+2*X)/(L*L*L);
				TempCM_DisLoadVector(8,0) = P*X*X*(L+2*(L-X))/(L*L*L);
				//
				P = m_Load.m_adElemLoadValue_Case1_Cen[loop].m_dElemMX;
				TempCM_DisLoadVector(3,0) = P*(L-X)/L;
				TempCM_DisLoadVector(9,0) = P*X/L;
				//
				P = m_Load.m_adElemLoadValue_Case1_Cen[loop].m_dElemMY;
				TempCM_DisLoadVector(4,0) = -6*P*X*(L-X)/(L*L*L) - P2*X*(L-X)*(L-X)/(L*L);
				TempCM_DisLoadVector(10,0) = 6*P*X*(L-X)/(L*L*L) + P2*X*X*(L-X)/(L*L);
				P = m_Load.m_adElemLoadValue_Case1_Cen[loop].m_dElemMZ;
				TempCM_DisLoadVector(5,0) = 6*P*X*(L-X)/(L*L*L) + P1*X*(L-X)*(L-X)/(L*L);
				TempCM_DisLoadVector(11,0) = -6*P*X*(L-X)/(L*L*L) - P1*X*X*(L-X)/(L*L);			
				
				//整体荷载向量
				TempCM_DisLoadVectorG = mat_TT * TempCM_DisLoadVector;				

				for(int loop2=0;loop2<(m_Node.m_nNodeNum - m_Node.m_nNode_K_Num);loop2++)
				{
					if(m_nNode1 == m_Node.m_aDOFIndex[loop2].m_aiNode)
					{
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_X;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase1_ElemCen_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(0,0);
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_Y;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase1_ElemCen_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(1,0);
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_Z;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase1_ElemCen_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(2,0);
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_RX;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase1_ElemCen_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(3,0);
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_RY;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase1_ElemCen_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(4,0);
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_RZ;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase1_ElemCen_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(5,0);						
					}
					if(m_nNode2 == m_Node.m_aDOFIndex[loop2].m_aiNode)
					{
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_X;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase1_ElemCen_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(6,0);
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_Y;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase1_ElemCen_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(7,0);
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_Z;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase1_ElemCen_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(8,0);
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_RX;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase1_ElemCen_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(9,0);
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_RY;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase1_ElemCen_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(10,0);
						m_nTemp1 = m_Node.m_aDOFIndex[loop2].m_aiDOFIndex_RZ;
						//if(m_nTemp1<m_Node.m_nFreeDOF)	
							m_adCase1_ElemCen_LoadVector[m_nTemp1] += TempCM_DisLoadVectorG(11,0);							
					}					
				}

				//修改单元的挠度
				L = m_apEle[0]->m_aBeamElemProValue[loop1].m_L;
				X = m_Load.m_adElemLoadValue_Case1_Cen[loop].m_dDistance;
				double PY = m_Load.m_adElemLoadValue_Case1_Cen[loop].m_dElemForceY;
				double PZ = m_Load.m_adElemLoadValue_Case1_Cen[loop].m_dElemForceZ;
				double E = m_apEle[0]->m_aBeamElemProValue[loop1].m_E;
				double IY = m_apEle[0]->m_aBeamElemProValue[loop1].m_Iy;
				double IZ = m_apEle[0]->m_aBeamElemProValue[loop1].m_Iz;
				for(int loop3=0;loop3<13;loop3++)
				{
					double x = loop3*L/12;
					if(x <= X)
					{
						m_apEle[0]->m_aElemDeflection[loop1].m_dDelf_Y[loop3] += 
							PY*(L-X)*(L-X)*x*x*( 3*X - (1+2*X/L)*x) /(6*E*IZ*L*L);
						m_apEle[0]->m_aElemDeflection[loop1].m_dDelf_Z[loop3] += 
							PZ*(L-X)*(L-X)*x*x*( 3*X - (1+2*X/L)*x) /(6*E*IY*L*L);
					}
					else
					{
						m_apEle[0]->m_aElemDeflection[loop1].m_dDelf_Y[loop3] += 
							-(PY*X*X*(L-x)*(L-x)*(X - (1+2*((L-X)/L))*x))/(6*E*IZ*L*L);
						m_apEle[0]->m_aElemDeflection[loop1].m_dDelf_Z[loop3] += 
							-(PZ*X*X*(L-x)*(L-x)*(X - (1+2*((L-X)/L))*x))/(6*E*IY*L*L);
					}
				}			
			}
		}
	}

	//初始化 CASE 2
	m_adCase2_ElemCen_LoadVector.SetSize(m_Node.m_nTotalDOF);
	for(loop=0;loop<m_Node.m_nTotalDOF;loop++)
		m_adCase2_ElemCen_LoadVector[loop] = 0;

	//具体赋值
	for(loop=0;loop<m_Load.m_nElemLoadNum_Case2_Cen;loop++)
	{
		//受荷单元号
		m_nTemp = m_Load.m_adElemLoadValue_Case2_Cen[loop].m_nElementNO;
		//作用点偏移
		X = m_Load.m_adElemLoadValue_Case2_Cen[loop].m_dDistance;
		if(X <= 0.0)
		{
			AfxMessageBox("单元荷载输入有误",MB_OK,1);
			break;
		}
		
		//单元对应节点
		for(int loop1=0;loop1<m_apEle[0]->m_nElementNum;loop1++)
		{
			if(m_nTemp == m_apEle[0]->m_aBeamElemProValue[loop1].m_nElementNo)
			{
				m_nNode1 = m_apEle[0]->m_aBeamElemProValue[loop1].m_nElemNode1;
				m_nNode2 = m_apEle[0]->m_aBeamElemProValue[loop1].m_nElemNode2;
				L = m_apEle[0]->m_aBeamElemProValue[loop1].m_L;
				if( X > L)
				{
					AfxMessageBox("单元荷载输入有误",MB_OK,1);