globalelement.cpp

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

// GlobalElement.cpp: implementation of the CGlobalElement class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "math.h"
#include "GlobalElement.h"

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

#define TRUSS 0
#define BEAM 1
#define BEAM_HINGE 2
#define PLANE_8NODE_ELE 3

#define MATERIAL_A 0
#define MATERIAL_B 1

#define TWO_DOF		0
#define THREE_DOF	1

#define NODE_LOAD 0
#define BEAM_LOAD_1 1
#define BEAM_LOAD_2 2
#define BEAM_LOAD_3 3
#define BEAM_LOAD_4 4
#define BEAM_LOAD_5 5

#define DISPLAY_ELEMENT 1
#define DISPLAY_MOMENT 2
#define DISPLAY_SHEAR 3
#define DISPLAY_DEFORMATION 4
#define DISPLAY_TSTRESS 5
#define DISPLAY_CSTRESS 6
#define DISPLAY_AXIAL_FORCE 10

#define CONSISTENT_MASS 0
#define LUMP_MASS 1

#define DYNAMIC_INSTANT_LOAD 0
#define DYNAMIC_SEISMIC_LOAD 1

#define PI 3.1415926

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

CGlobalElement::CGlobalElement()
{
	m_nMode=3;
	m_iCurLoadGroup=0;
	m_nKindMaterial=2;
	m_iMassMatrixType=CONSISTENT_MASS;
	m_iDynamicLoadType=DYNAMIC_INSTANT_LOAD;
	m_adDampRatio.SetSize(m_nMode);
	m_adDampRatio[0]=m_adDampRatio[1]=0.05;
	m_dResponseDuration=10.0;
	m_dLoadDuration=0.3;
	m_dTimeStep=0.01;
	m_dPeakAcc=1.0;
	m_sGroundAccFile="C:\\";
	m_dNewmarkGama=0.5;
	m_dNewmarkBeta=0.25;
	m_dWilsonXita=1.4;
}

CGlobalElement::~CGlobalElement()
{
	Initial();
}

void CGlobalElement::Initial()
{
	int loop,loop1,iBuf,iBuf1;
	m_bFlagForceAnalyzed=false;
	m_bFlagModeAnalyzed=false;
	m_bFlagDynamicAnalyzed=false;
	iBuf=m_apEle.GetSize();
	for(loop=0;loop<iBuf;loop++)
		delete m_apEle[loop];
	m_apEle.RemoveAll();

	if(iBuf!=0){
		delete m_adDisp;
		delete m_adLoadVector;
	}

	iBuf=m_apMaterial.GetSize();
	for(loop=0;loop<iBuf;loop++)
		delete m_apMaterial[loop];
	m_apMaterial.RemoveAll();

	iBuf=m_apLoadGroup.GetSize();
	for(loop=0;loop<iBuf;loop++){
		iBuf1=m_apLoadGroup[loop]->m_apLoad.GetSize();
		for(loop1=0;loop1<iBuf1;loop1++){
			delete m_apLoadGroup[loop]->m_apLoad[loop1];
		}
		m_apLoadGroup[loop]->m_apLoad.RemoveAll();
		delete m_apLoadGroup[loop];
	}
	m_apLoadGroup.RemoveAll();

	iBuf=m_LumpMass.m_apLoad.GetSize();
	for(loop1=0;loop1<iBuf;loop1++){
		delete m_LumpMass.m_apLoad[loop1];
	}
	m_LumpMass.m_apLoad.RemoveAll();
	m_Node.Initial();
	m_iCurLoadGroup=0;
}

void CGlobalElement::DataInput(CString &sData)
{
	int loop,loop1,iBuf;
	int iCurCharPos;
	int nLoadGroup,nLoad;
	double dBuf;
	int iEle,iEleType,iNode,iNodeType;
	int iMaterialType;
	int nConstrainedNode,iIsFixedX,iIsFixedY,iIsFixedR;
	int iMaterialIndex;
	CString sBuf;
	
	Initial();

	CBaseElement::m_pNode=&m_Node;
	CLoad::m_papEle=&m_apEle;
	CLoad::m_pNode=&m_Node;
	CLoad::m_papMaterial=&m_apMaterial;
	CBaseElement::m_papMaterial=&m_apMaterial;

	sBuf="num. of node\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return;
	m_nNode=ReadInt(iCurCharPos,sData);
	m_Node.SetNum(m_nNode);

	sBuf="node type\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return;
	for(loop=0;loop<m_nNode;loop++){
		iNode=ReadInt(iCurCharPos,sData);
		iNodeType=ReadInt(iCurCharPos,sData);
		m_Node.SetType(iNode,iNodeType);
	}

	sBuf="coor\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return;
	for(loop=0;loop<m_nNode;loop++){
		iNode=ReadInt(iCurCharPos,sData);
		dBuf=ReadDouble(iCurCharPos,sData);
		m_Node.SetNodeX(iNode,dBuf);
		dBuf=ReadDouble(iCurCharPos,sData);
		m_Node.SetNodeY(iNode,dBuf);
	}
	
	sBuf="num. constrained node\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return;
	nConstrainedNode=ReadInt(iCurCharPos,sData);
	for(loop=0;loop<nConstrainedNode;loop++){
		iNode=ReadInt(iCurCharPos,sData);
		iIsFixedX=ReadInt(iCurCharPos,sData);
		iIsFixedY=ReadInt(iCurCharPos,sData);
		iIsFixedR=ReadInt(iCurCharPos,sData);
		m_Node.SetConstrainedNodeDescription(iNode,iIsFixedX,iIsFixedY,iIsFixedR);
	}

	sBuf="num. of element\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return;
	m_nEle=ReadInt(iCurCharPos,sData);
	m_apEle.SetSize(m_nEle);
	
	sBuf="element type\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return;
	for(loop=0;loop<m_nEle;loop++){
		iEle=ReadInt(iCurCharPos,sData);
		iEleType=ReadInt(iCurCharPos,sData);
		switch(iEleType){
		case TRUSS:
			m_apEle[iEle]=new CTruss();
			break;
		case BEAM:
			m_apEle[iEle]=new CBeam();
			break;
		case BEAM_HINGE:
			m_apEle[iEle]=new CBeamHinge();
			break;
		case PLANE_8NODE_ELE:
			m_apEle[iEle]=new CPlane8NodeEle();
			break;
		}
	}

	sBuf="element node\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return;

	for(loop=0;loop<m_nEle;loop++){
		iEle=ReadInt(iCurCharPos,sData);
		m_apEle[iEle]->ReadNode(iCurCharPos,sData);
	}

	sBuf="element parameter\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return;
	for(loop=0;loop<m_nEle;loop++){
		iEle=ReadInt(iCurCharPos,sData);
		m_apEle[iEle]->ReadParameter(iCurCharPos,sData);
	}

	sBuf="num. of material\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return;
	iBuf=ReadInt(iCurCharPos,sData);
	m_apMaterial.SetSize(iBuf);
	sBuf="material parameter\0";
	for(loop=0;loop<iBuf;loop++){
		iMaterialIndex=ReadInt(iCurCharPos,sData);
		iMaterialType=ReadInt(iCurCharPos,sData);
		if(iMaterialIndex>=iBuf||iMaterialType>=m_nKindMaterial){
			sBuf="Error in material index!";
			AfxMessageBox( sBuf, MB_OK, 0 );
			return;
		}
		switch(iMaterialType){
		case MATERIAL_A:
			m_apMaterial[iMaterialIndex]=new CMaterialA();
			break;
		case MATERIAL_B:
			m_apMaterial[iMaterialIndex]=new CMaterialB();
			break;
		}
		m_apMaterial[loop]->ReadParameter(iCurCharPos,sData);
	}
	
	sBuf="num. of load group\0";
	if((iCurCharPos=sData.Find(sBuf))>=0){
		iCurCharPos+=sBuf.GetLength();
		nLoadGroup=ReadInt(iCurCharPos,sData);		
		for(loop=0;loop<nLoadGroup;loop++){
			m_apLoadGroup.Add(new CLoadGroup());
			sBuf.Format("load group %d",loop);
			if((iCurCharPos=sData.Find(sBuf))>=0){
				iCurCharPos+=sBuf.GetLength();
				nLoad=ReadInt(iCurCharPos,sData);
				for(loop1=0;loop1<nLoad;loop1++){
					m_apLoadGroup[loop]->m_apLoad.Add(new CLoad());
					m_apLoadGroup[loop]->m_apLoad[loop1]->ReadParameter(iCurCharPos,sData);
				}
			}
		}
	}

	int nLumpMass;
	sBuf="lump mass\0";
	if((iCurCharPos=sData.Find(sBuf))>=0){
		nLumpMass=ReadInt(iCurCharPos,sData);
		for(loop=0;loop<nLumpMass;loop++){
			m_LumpMass.m_apLoad.Add(new CLoad());
			m_LumpMass.m_apLoad[loop]->ReadMassParameter(iCurCharPos,sData);
		}
	}

	m_Node.CalcuDOFIndex();
	iBuf=m_Node.GetTotalDOF();
	m_adDisp=new double[iBuf];
	m_adLoadVector=new double[iBuf];
	CalcuGKDiagAddress();	
}

int CGlobalElement::ReadInt(int &iCurCharPos, const CString &sData)
{
	int nCount,nLength;
	char cBuf,sBuf[50];

	nLength=sData.GetLength();
	cBuf=sData[iCurCharPos];
	while(cBuf<(char)43||cBuf>(char)57||cBuf==(char)44||cBuf==(char)47)
	{
		iCurCharPos++;
		if(iCurCharPos==nLength) return -1;
		cBuf=sData[iCurCharPos];
	}
	
	nCount=0;
	while((cBuf>=(char)48&&cBuf<=(char)57)||cBuf==(char)43||cBuf==(char)45)
	{
		sBuf[nCount]=cBuf;
		iCurCharPos++;
		cBuf=sData[iCurCharPos];
		nCount++;
	}
	sBuf[nCount]='\0';
	return atoi(sBuf);
}

double CGlobalElement::ReadDouble(int &iCurCharPos, const CString &sData)
{
	char cBuf,sBuf[50],*sStopBuf;
	int nCount,nLength;
	
	nLength=sData.GetLength();
	cBuf=sData[iCurCharPos];
	while(cBuf<(char)43||cBuf>(char)57||cBuf==(char)44||cBuf==(char)47)
	{
		iCurCharPos++;
		if(iCurCharPos==nLength) return -1;
		cBuf=sData[iCurCharPos];
	}
	
	nCount=0;
	while((cBuf>=(char)48&&cBuf<=(char)57)||
		cBuf==(char)43||cBuf==(char)45||cBuf==(char)46||
		cBuf==(char)69||cBuf==(char)101)
	{
		sBuf[nCount]=cBuf;
		iCurCharPos++;
		if(iCurCharPos==nLength) break;
		cBuf=sData[iCurCharPos];
		nCount++;
	}
	sBuf[nCount]='\0';
	return strtod( sBuf, &sStopBuf );
}

void CGlobalElement::Loading()
{
	int loop,nLoadGroup;
	int nTotalDOF,nFreeDOF;
	CString sBuf;
	
	m_bFlagForceAnalyzed=true;
	nLoadGroup=m_apLoadGroup.GetSize();
	if(m_iCurLoadGroup==nLoadGroup) return;

	nTotalDOF=m_Node.GetTotalDOF();
	nFreeDOF=m_Node.GetFreeDOF();

	for(loop=0;loop<m_nEle;loop++)
		m_apEle[loop]->InternalForceInitial();

	m_smatGK.ElementBufRealloc();
	m_smatGK=0.0;
	for(loop=0;loop<m_nEle;loop++){
		m_apEle[loop]->StiffAssemble(m_smatGK);
	}

	GetLoadVector();
	for(loop=0;loop<nFreeDOF;loop++)
		m_adDisp[loop]=m_adLoadVector[loop];

	if(!m_smatGK.LdltSolve(nFreeDOF,m_adDisp)) return;
	
	m_Node.SetCurCoordinate(m_adDisp);
	for(loop=0;loop<m_nEle;loop++){
		m_apEle[loop]->CalcuInternalForce(m_adDisp);
	}
	CalcuSupportReaction();
	m_smatGK.ElementBufEmpty();
	m_iCurLoadGroup++;
}

void CGlobalElement::GetLoadVector()
{
	int loop,loop1,nLoad;
	int nTotalDOF,nFreeDOF;
	double dBuf;
	
	nTotalDOF=m_Node.GetTotalDOF();
	nFreeDOF=m_Node.GetFreeDOF();

	for(loop=0;loop<nTotalDOF;loop++){
		m_adLoadVector[loop]=0.0;
		m_adDisp[loop]=0.0;
	}
	
	nLoad=m_apLoadGroup[m_iCurLoadGroup]->m_apLoad.GetSize();
	for(loop=0;loop<nLoad;loop++){
		m_apLoadGroup[m_iCurLoadGroup]->m_apLoad[loop]
			->LoadVectorAssemble(m_adLoadVector,m_adDisp);
	}

	for(loop=nFreeDOF;loop<nTotalDOF;loop++){
		if(m_adDisp[loop]!=0.0){
			for(loop1=0;loop1<nFreeDOF;loop1++){
				dBuf=m_smatGK(loop,loop1);
				if(dBuf!=0.0){
					m_adLoadVector[loop1]-=dBuf*m_adDisp[loop];
				}
			}
		}
	}
}

int CGlobalElement::FindKeyWordsPos(CString &sKeyWords, CString &sData)
{
	int	iCharPos;
	CString sBuf;
	if((iCharPos=sData.Find(sKeyWords))<0){
		sBuf.Format("Can not find the data: %s",sKeyWords);
		AfxMessageBox( sBuf, MB_OK, 0 );
	}
	else{
		iCharPos+=sKeyWords.GetLength();
	}
	return iCharPos;
}

void CGlobalElement::CalcuGKDiagAddress()
{
	int loop;
	int nTotalDOF;
	unsigned long *aiGKDiagAdd;
		
	nTotalDOF=m_Node.GetTotalDOF();
	aiGKDiagAdd=new unsigned long [nTotalDOF];
	
	for(loop=1;loop<nTotalDOF;loop++){
		aiGKDiagAdd[loop]=1;
	}
	for(loop=0;loop<m_nEle;loop++){
		m_apEle[loop]->CalcuGKBandWidth(aiGKDiagAdd);
	}
	aiGKDiagAdd[0]=0;	
	for(loop=1;loop<nTotalDOF;loop++)
		aiGKDiagAdd[loop]+=aiGKDiagAdd[loop-1];
	m_smatGK.Realloc(nTotalDOF,aiGKDiagAdd);
	delete aiGKDiagAdd;
}

void CGlobalElement::OutputData(ofstream& fout)
{