beam.cpp

一个计算悬臂梁的有限元vc源码

// Beam.cpp: implementation of the CBeam class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "afxtempl.h"
#include "fstream.h"
#include "math.h"
#include "Matrix.h"
#include "SparseMatrix.h"
#include "Node.h"
#include "BaseMaterial.h"
#include "Beam.h"

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

#define BEAM 1

CMatrix CBeam::m_matT=CMatrix(6,6);
CMatrix CBeam::m_matTT=CMatrix(6,6);
CMatrix CBeam::m_matKe=CMatrix(6,6);
CMatrix CBeam::m_matNodeDisp=CMatrix(6,1);
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CBeam::CBeam()
{

}

CBeam::~CBeam()
{

}

void CBeam::StiffAssemble(CSparseMatrix& smatGK)
{
	int loop,loop1,iBuf,aiDOFIndex[6];
	//int nFreeDOF;
	GetStiffness();
	for(loop=0;loop<2;loop++){
		aiDOFIndex[3*loop]=m_pNode->GetXDOFIndex(m_aiNode[loop]);
		aiDOFIndex[3*loop+1]=m_pNode->GetYDOFIndex(m_aiNode[loop]);
		aiDOFIndex[3*loop+2]=m_pNode->GetRDOFIndex(m_aiNode[loop]);
	}
	//nFreeDOF=m_pNode->GetFreeDOF();
	for(loop=0;loop<6;loop++){
		//if(aiDOFIndex[loop]<nFreeDOF){
			iBuf=aiDOFIndex[loop];
			smatGK(iBuf,iBuf)+=m_matKe(loop,loop);
			for(loop1=0;loop1<loop;loop1++){
				//if(aiDOFIndex[loop1]<nFreeDOF){
					smatGK(iBuf,aiDOFIndex[loop1])+=m_matKe(loop,loop1);
				//}
			}
		//}
	}
}

void CBeam::CalcuGKBandWidth(unsigned long* aiGKDiagAdd)
{
	int loop,iBuf,iMinDOFIndex,aiDOFIndex[6];

	for(loop=0;loop<2;loop++)
	{	aiDOFIndex[3*loop]=m_pNode->GetXDOFIndex(m_aiNode[loop]);
		aiDOFIndex[3*loop+1]=m_pNode->GetYDOFIndex(m_aiNode[loop]);
		aiDOFIndex[3*loop+2]=m_pNode->GetRDOFIndex(m_aiNode[loop]);
	}

	iMinDOFIndex=aiDOFIndex[0];
	for(loop=1;loop<6;loop++){
		if(iMinDOFIndex>aiDOFIndex[loop]) 
			iMinDOFIndex=aiDOFIndex[loop];
	}

	for(loop=0;loop<6;loop++){
		iBuf=aiDOFIndex[loop]-iMinDOFIndex+1;
		if(aiGKDiagAdd[aiDOFIndex[loop]]<iBuf)
			aiGKDiagAdd[aiDOFIndex[loop]]=iBuf;
	}
}
  
void CBeam::GetStiffness()
{
	double dBuf;
	double dE;
	CTypedPtrArray<CPtrArray,CBaseMaterial*>& apMaterial=*m_papMaterial;

	m_matT=0.0;
	m_matT(2,2)=m_matT(5,5)=1.0;	
	m_matT(0,0)=m_matT(1,1)=m_matT(3,3)=m_matT(4,4)=m_dCos;
	m_matT(0,1)=m_matT(3,4)=m_dSin;
	m_matT(1,0)=m_matT(4,3)=-m_dSin;
	m_matTT.Trans(m_matT);

	dE=apMaterial[m_iMaterialIndex]->GetE();

	m_matKe=0.0;
	dBuf=dE*m_dArea/m_dLength;
	m_matKe(0,0)=m_matKe(3,3)=dBuf;
	m_matKe(0,3)=m_matKe(3,0)=-dBuf;
	dBuf=m_dLength*m_dLength*m_dLength;
	dBuf=12.0*dE*m_dInertia/dBuf;
	m_matKe(1,1)=m_matKe(4,4)=dBuf;
	m_matKe(1,4)=m_matKe(4,1)=-dBuf;
	dBuf=m_dLength*m_dLength;
	dBuf=6.0*dE*m_dInertia/dBuf;
	m_matKe(1,2)=m_matKe(2,1)=m_matKe(1,5)
	=m_matKe(5,1)=dBuf;
	m_matKe(2,4)=m_matKe(4,2)=m_matKe(4,5)
	=m_matKe(5,4)=-dBuf;
	dBuf=4.0*dE*m_dInertia/m_dLength;
	m_matKe(2,2)=m_matKe(5,5)=dBuf;
	m_matKe(5,2)=m_matKe(2,5)=dBuf/2.0;

	m_matKe=m_matTT*m_matKe;
	m_matKe*=m_matT;
}

void CBeam::InternalForceInitial()
{
	for(int loop=0;loop<11;loop++)
		m_adM[loop]=m_adQ[loop]=m_adFy[loop]=0.0;
	m_dN=0.0;
}

void CBeam::CalcuInternalForce(const double* adNodeDisp)
{
	int loop,iBuf;
	double dx,dDx,dMi,dMj,dBuf,dBuf1;
	double dKxi;

	for(loop=0;loop<2;loop++){
		iBuf=m_pNode->GetXDOFIndex(m_aiNode[loop]);
		m_matNodeDisp(3*loop,0)=adNodeDisp[iBuf];
		iBuf=m_pNode->GetYDOFIndex(m_aiNode[loop]);
		m_matNodeDisp(3*loop+1,0)=adNodeDisp[iBuf];
		iBuf=m_pNode->GetRDOFIndex(m_aiNode[loop]);
		m_matNodeDisp(3*loop+2,0)=adNodeDisp[iBuf];
	}
	
	GetStiffness();
	m_matNodeDisp=m_matKe*m_matNodeDisp;
	
	GetTransferMatrixGCToLC(m_matT);
	m_matNodeDisp=m_matT*m_matNodeDisp;
	
	m_dN=-m_matNodeDisp(0,0);

	dMi=-m_matNodeDisp(2,0);
	dMj=m_matNodeDisp(5,0);
	m_adM[0]+=dMi;
	m_adM[10]+=dMj;
	m_adQ[0]-=m_matNodeDisp(1,0);
	m_adQ[10]-=m_matNodeDisp(1,0);

	dBuf=(dMi-dMj)/10.0;
	dx=dDx=m_dLength/10.0;
	for(loop=1;loop<10;loop++){
		m_adM[loop]+=dMi-dBuf*(double)loop;
		m_adQ[loop]-=m_matNodeDisp(1,0);
	}

	for(loop=0;loop<2;loop++){
		iBuf=m_pNode->GetXDOFIndex(m_aiNode[loop]);
		m_matNodeDisp(3*loop,0)=adNodeDisp[iBuf];
		iBuf=m_pNode->GetYDOFIndex(m_aiNode[loop]);
		m_matNodeDisp(3*loop+1,0)=adNodeDisp[iBuf];
		iBuf=m_pNode->GetRDOFIndex(m_aiNode[loop]);
		m_matNodeDisp(3*loop+2,0)=adNodeDisp[iBuf];
	}
	GetTransferMatrixGCToLC(m_matT);
	m_matNodeDisp=m_matT*m_matNodeDisp;

	m_adFy[0]-=m_matNodeDisp(1,0);
	m_adFy[10]-=m_matNodeDisp(4,0);
	for(loop=1;loop<10;loop++){
		dKxi=loop/10.0;
		dBuf=dKxi*dKxi;
		dBuf1=dKxi*dKxi*dKxi;
		m_adFy[loop]-=(1.0-3.0*dBuf+2.0*dBuf1)*m_matNodeDisp(1,0)
					  +(dKxi-2.0*dBuf+dBuf1)*m_dLength*m_matNodeDisp(2,0)
					  +(3.0*dBuf-2.0*dBuf1)*m_matNodeDisp(4,0)
					  +(dBuf1-dBuf)*m_dLength*m_matNodeDisp(5,0);
	}
	m_dU=m_matNodeDisp(0,0);
	m_dU1=m_matNodeDisp(3,0);
}

void CBeam::GetNode(int *aiNode, int& nEleNode)
{
	aiNode[0]=m_aiNode[0];
	aiNode[1]=m_aiNode[1];
	nEleNode=2;
}	

void CBeam::GetTransferMatrixGCToLC(CMatrix &matT)
{
	matT=0.0;
	matT(2,2)=matT(5,5)=1.0;	
	matT(0,0)=matT(1,1)=matT(3,3)=matT(4,4)=m_dCos;
	matT(0,1)=matT(3,4)=m_dSin;
	matT(1,0)=matT(4,3)=-m_dSin;
}

double CBeam::GetLength()
{
	return m_dLength;
}

void CBeam::SetInitialSectionForce(double *adQInit, double *adMInit)
{
	int loop;
	for(loop=0;loop<11;loop++){
		m_adQ[loop]+=adQInit[loop];
		m_adM[loop]+=adMInit[loop];
	}
}

void CBeam::SetInitialSectionDef(double *adFyInit)
{
	for(int loop=0;loop<11;loop++){
		m_adFy[loop]+=adFyInit[loop];
	}
}

void CBeam::ReadNode(int& iCurCharPos, CString &sData)
{
	double dX,dX1,dY,dY1;

	m_aiNode[0]=ReadInt(iCurCharPos,sData);
	m_aiNode[1]=ReadInt(iCurCharPos,sData);

	dX=m_pNode->GetX(m_aiNode[0]);	
	dY=m_pNode->GetY(m_aiNode[0]);	
	dX1=m_pNode->GetX(m_aiNode[1]);	
	dY1=m_pNode->GetY(m_aiNode[1]);
	m_dLength=sqrt((dY1-dY)*(dY1-dY)+(dX1-dX)*(dX1-dX));
	m_dCos=(dX1-dX)/m_dLength;
	m_dSin=(dY1-dY)/m_dLength;
}

void CBeam::ReadParameter(int &iCurCharPos, CString &sData)
{
	m_dArea=ReadDouble(iCurCharPos,sData);
	m_dInertia=ReadDouble(iCurCharPos,sData);
	m_iMaterialIndex=ReadInt(iCurCharPos,sData);
}

double CBeam::GetInertia()
{
	return m_dInertia;
}

int CBeam::GetElementType()
{
	return BEAM;
}

void CBeam::OutputInternalForce(ofstream &fout)
{
	int loop;

	fout<<"N = "<<m_dN<<endl;
	fout<<"M:  ";
	for(loop=0;loop<11;loop++){
		fout<<m_adM[loop]<<"  ";
	}
	fout<<endl;

	fout<<"Q:  ";	
	for(loop=0;loop<11;loop++){
		fout<<m_adQ[loop]<<"  ";
	}
	fout<<endl;

	fout<<"Fy:  ";	
	for(loop=0;loop<11;loop++){
		fout<<m_adFy[loop]<<"  ";
	}
	fout<<endl;
}	

void CBeam::OutputParameter(ofstream &fout)
{
	fout<<"Beam  面积:"<<m_dArea<<"  惯性矩:"<<m_dInertia;
	fout<<"  材料类型编号:"<<m_iMaterialIndex;
	fout<<"  结点号:"<<m_aiNode[0]<<"  "<<m_aiNode[1]<<endl;
}

void CBeam::GetSectionForce(double& dN,double *adQ, double *adM)
{
	dN=m_dN;
	for(int loop=0;loop<11;loop++){
		adM[loop]=m_adM[loop];
		adQ[loop]=m_adQ[loop];
	}
}


void CBeam::GetSectionDef(double *adFy,double&dU,double&dU1)
{
	for(int loop=0;loop<11;loop++)
		adFy[loop]=m_adFy[loop];

	dU=m_dU;
	dU1=m_dU1;
}

void CBeam::GetMassMatrix()
{
	double dBuf;
	double dDensity;
	CTypedPtrArray<CPtrArray,CBaseMaterial*>& apMaterial=*m_papMaterial;

	m_matT=0.0;
	m_matT(2,2)=m_matT(5,5)=1.0;	
	m_matT(0,0)=m_matT(1,1)=m_matT(3,3)=m_matT(4,4)=m_dCos;
	m_matT(0,1)=m_matT(3,4)=m_dSin;
	m_matT(1,0)=m_matT(4,3)=-m_dSin;
	m_matTT.Trans(m_matT);

	dDensity=apMaterial[m_iMaterialIndex]->GetDensity();

	m_matKe=0.0;
	dBuf=dDensity*m_dArea*m_dLength/420.0;
	m_matKe(0,0)=m_matKe(3,3)=140.0*dBuf;
	m_matKe(0,3)=m_matKe(3,0)=70.0*dBuf;
	m_matKe(1,1)=m_matKe(4,4)=156.0*dBuf;
	m_matKe(1,4)=m_matKe(4,1)=54.0*dBuf;
	m_matKe(1,2)=m_matKe(2,1)=22.0*dBuf*m_dLength;
	m_matKe(1,5)=m_matKe(5,1)=-13.0*dBuf*m_dLength;
	m_matKe(2,2)=m_matKe(5,5)=4.0*dBuf*m_dLength*m_dLength;
	m_matKe(2,4)=m_matKe(4,2)=13.0*dBuf*m_dLength;
	m_matKe(2,5)=m_matKe(5,2)=-3.0*dBuf*m_dLength*m_dLength;
	m_matKe(4,5)=m_matKe(5,4)=-22.0*dBuf*m_dLength;

	m_matKe=m_matTT*m_matKe;
	m_matKe*=m_matT;
}

void CBeam::MassAssemble(CSparseMatrix& smatMass)
{
	int loop,loop1,iBuf,aiDOFIndex[6];
	//int nFreeDOF;
	GetMassMatrix();
	for(loop=0;loop<2;loop++){
		aiDOFIndex[3*loop]=m_pNode->GetXDOFIndex(m_aiNode[loop]);
		aiDOFIndex[3*loop+1]=m_pNode->GetYDOFIndex(m_aiNode[loop]);
		aiDOFIndex[3*loop+2]=m_pNode->GetRDOFIndex(m_aiNode[loop]);
	}
	//nFreeDOF=m_pNode->GetFreeDOF();
	for(loop=0;loop<6;loop++){
		//if(aiDOFIndex[loop]<nFreeDOF){
			iBuf=aiDOFIndex[loop];
			smatMass(iBuf,iBuf)+=m_matKe(loop,loop);
			for(loop1=0;loop1<loop;loop1++){
				//if(aiDOFIndex[loop1]<nFreeDOF){
					smatMass(iBuf,aiDOFIndex[loop1])+=m_matKe(loop,loop1);
				//}
			}
		//}
	}
}