truss.cpp

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

// Truss.cpp: implementation of the CTruss 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 "Truss.h"

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

#define TRUSS 0

CMatrix CTruss::m_matT=CMatrix(4,4);
CMatrix CTruss::m_matTT=CMatrix(4,4);
CMatrix CTruss::m_matKe=CMatrix(4,4);
CMatrix CTruss::m_matNodeDisp=CMatrix(4,1);

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

CTruss::CTruss()
{
}

CTruss::~CTruss()
{

}

void CTruss::GetStiffness()
{
	double dBuf,dBuf1;
	double dE;
	CTypedPtrArray<CPtrArray,CBaseMaterial*>& apMaterial=*m_papMaterial;
	
	dE=apMaterial[m_iMaterialIndex]->GetE();
	dBuf=dE*m_dArea/m_dLength;
	dBuf1=m_dCos*m_dCos*dBuf;
	m_matKe(0,0)=m_matKe(2,2)=dBuf1;
	m_matKe(0,2)=m_matKe(2,0)=-dBuf1;
	dBuf1=m_dSin*m_dSin*dBuf;
	m_matKe(1,1)=m_matKe(3,3)=dBuf1;
	m_matKe(1,3)=m_matKe(3,1)=-dBuf1;
	dBuf1=m_dCos*m_dSin*dBuf;
	m_matKe(0,1)=m_matKe(1,0)=m_matKe(2,3)=m_matKe(3,2)=dBuf1;
	m_matKe(0,3)=m_matKe(3,0)=m_matKe(1,2)=m_matKe(2,1)=-dBuf1;

}

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

	for(loop=0;loop<2;loop++)
	{	aiDOFIndex[2*loop]=m_pNode->GetXDOFIndex(m_aiNode[loop]);
		aiDOFIndex[2*loop+1]=m_pNode->GetYDOFIndex(m_aiNode[loop]);
	}
	
	iMinDOFIndex=aiDOFIndex[0];
	for(loop=1;loop<4;loop++){
		if(iMinDOFIndex>aiDOFIndex[loop]) 
			iMinDOFIndex=aiDOFIndex[loop];
	}
	for(loop=0;loop<4;loop++){
		iBuf=aiDOFIndex[loop]-iMinDOFIndex+1;
		if(aiGKDiagAdd[aiDOFIndex[loop]]<iBuf)
			aiGKDiagAdd[aiDOFIndex[loop]]=iBuf;
	}
}

void CTruss::StiffAssemble(CSparseMatrix& smatGK)
{
	int loop,loop1,iBuf,aiDOFIndex[4];
	//int nFreeDOF;
	GetStiffness();
	for(loop=0;loop<2;loop++){
		aiDOFIndex[2*loop]=m_pNode->GetXDOFIndex(m_aiNode[loop]);
		aiDOFIndex[2*loop+1]=m_pNode->GetYDOFIndex(m_aiNode[loop]);
	}
	//nFreeDOF=m_pNode->GetFreeDOF();
	for(loop=0;loop<4;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  CTruss::CalcuInternalForce(const double* adNodeDisp)
{
	int loop,iBuf;

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

	GetTransferMatrixGCToLC(m_matT);
	m_matNodeDisp=m_matT*m_matNodeDisp;
	m_dN=m_matNodeDisp(2,0);	
}


void CTruss::InternalForceInitial()
{
	m_dN=0;
}

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

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

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

void CTruss::ReadNode(int& iCurCharPos, CString &sData)
{
	double dX,dY,dX1,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 CTruss::ReadParameter(int &iCurCharPos, CString &sData)
{
	m_dArea=ReadDouble(iCurCharPos,sData);
	m_iMaterialIndex=ReadInt(iCurCharPos,sData);
}

int CTruss::GetElementType()
{
	return TRUSS;
}

void CTruss::OutputInternalForce(ofstream &fout)
{
	fout<<"N = "<<m_dN<<endl;
}

void CTruss::OutputParameter(ofstream &fout)
{
	fout<<"Truss  面积:"<<m_dArea;
	fout<<"  材料类型编号:"<<m_iMaterialIndex;
	fout<<"  结点号:"<<m_aiNode[0]<<"  "<<m_aiNode[1]<<endl;
}

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

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

	dDensity=apMaterial[m_iMaterialIndex]->GetDensity();
	m_matKe=0.0;
	dBuf=dDensity*m_dArea*m_dLength;
	m_matKe(0,0)=m_matKe(1,1)=m_matKe(2,2)
				=m_matKe(3,3)=dBuf/3.0;
	m_matKe(0,2)=m_matKe(2,0)=m_matKe(1,3)
				=m_matKe(3,1)=dBuf/6.0;

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

void CTruss::MassAssemble(CSparseMatrix& smatMass)
{
	int loop,loop1,iBuf,aiDOFIndex[4];
	//int nFreeDOF;
	GetMassMatrix();
	for(loop=0;loop<2;loop++){
		aiDOFIndex[2*loop]=m_pNode->GetXDOFIndex(m_aiNode[loop]);
		aiDOFIndex[2*loop+1]=m_pNode->GetYDOFIndex(m_aiNode[loop]);
	}
	//nFreeDOF=m_pNode->GetFreeDOF();
	for(loop=0;loop<4;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);
				//}
			}
		//}
	}
}

void CTruss::GetAxialForce(double &dN)
{
	dN=m_dN;
}