structure.cpp

自编的有限元方法解决弹性力学平面问题的软件。有图形显示

// Structure.cpp: implementation of the CStructure class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "FiniteElement.h"
#include "Structure.h"
#include "Element.h"
#include "Node.h"

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

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

CStructure::CStructure()
{
	Initial();
}

CStructure::~CStructure()
{
	Initial();
}
void CStructure::Initial()
{

	m_title=_T("弹性力学平面问题有限元分析");
	m_element_num=0;		//结构单元总数
	m_node_num=0;			//结构节点总数
	m_supportnode_num=0;	//支承节点数
	m_loadnode_num=0;		//有荷载节点数
	m_question_id=0;		//问题类型码
	//材料性质参数
	m_young=0.0;			//杨氏弹性模量	
	m_possion=0.0;		//泊松比
	m_thickness=0.0;		//厚度
	m_weight=0.0;		//容重
	//
	m_halfband=0;
	//元素对象
	for(int i=0;i<m_nodeArray.GetSize();i++)
	{
		delete (CNode*)(m_nodeArray.GetAt(i));
		m_nodeArray.RemoveAll();
	}
	for(int j=0;j<m_elementArray.GetSize();j++)
	{
		delete (CElement*)(m_elementArray.GetAt(j));
		m_elementArray.RemoveAll();
	}

}
/////////////////////////////////////////////////////////////////////

int CStructure::DataInput(CString &sData)
{
	int i,j;//循环变量
	int id;	//单元结点标志
	int tag1,tag2,tag3;//结点码
	double dbuf1=0.0,dbuf2=0.0;//双精度缓冲区
	int iCurCharPos;		//当前文件光标位置
	CString sBuf;			//字符串缓冲区
	Initial();	
	//-----------------------------------------------------
	sBuf="结点总数\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return 0;
	m_node_num=ReadInt(iCurCharPos,sData);
	sBuf="单元总数\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return 0;
	m_element_num=ReadInt(iCurCharPos,sData);
	sBuf="支承结点数\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return 0;
	m_supportnode_num=ReadInt(iCurCharPos,sData);
	sBuf="结点荷载数\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return 0;
	m_loadnode_num=ReadInt(iCurCharPos,sData);
	sBuf="问题类型码\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return 0;
	m_question_id=ReadInt(iCurCharPos,sData);
	//------------------------------------------------------
	sBuf="弹性模量\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return 0;
	m_young=ReadDouble(iCurCharPos,sData);
	sBuf="泊松比\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return 0;
	m_possion=ReadDouble(iCurCharPos,sData);
	sBuf="单元厚度\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return 0;
	m_thickness=ReadDouble(iCurCharPos,sData);
	sBuf="容重\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return 0;
	m_weight=ReadDouble(iCurCharPos,sData);
	//----------------------给结点和单元分配内存空间--------
	m_nodeArray.SetSize(m_node_num);
	m_elementArray.SetSize(m_element_num);
	for(i=0;i<m_node_num;i++)
	{
		m_nodeArray.SetAt(i,new CNode(i+1));
	}
	for(j=0;j<m_element_num;j++)
	{
		m_elementArray.SetAt(j,new CElement(j+1));
		((CElement*)(m_elementArray.GetAt(j)))->SetMaterialInfo(
			m_young,m_possion,m_thickness,m_weight);
	}
	//---------读入结点坐标数组-----------------------------
	sBuf="结点坐标数组\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return 0;
	for(j=0;j<m_node_num;j++){
		id=ReadInt(iCurCharPos,sData);
		dbuf1=ReadDouble(iCurCharPos,sData);
		dbuf2=ReadDouble(iCurCharPos,sData);
		((CNode*)(m_nodeArray.GetAt(id-1)))->SetNodeXY(dbuf1,dbuf2);
	}
	//---------读入单元结点数组-----------------------------
	sBuf="单元结点码数据\0";//在这里浪费了N个小时。要仔细啊！
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return 0;
	for(j=0;j<m_element_num;j++){
		id=ReadInt(iCurCharPos,sData);
		tag1=ReadInt(iCurCharPos,sData);	
		tag2=ReadInt(iCurCharPos,sData);
		tag3=ReadInt(iCurCharPos,sData);
		((CElement*)(m_elementArray.GetAt(id-1)))->SetNode(
			(CNode*)(m_nodeArray.GetAt(tag1-1)),
			(CNode*)(m_nodeArray.GetAt(tag2-1)),
			(CNode*)(m_nodeArray.GetAt(tag3-1)));
	}
	//----------读取支承结点信息----------------------------
	sBuf="支承结点数组\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return 0;
	for(i=0;i<m_supportnode_num;i++){
		id=ReadInt(iCurCharPos,sData);
		tag1=ReadInt(iCurCharPos,sData);
		tag2=ReadInt(iCurCharPos,sData);
		((CNode*)(m_nodeArray.GetAt(id-1)))->SetSupportXY(tag1,tag2);
	}
	//----------读取结点荷载信息----------------------------
	sBuf="结点荷载数组\0";
	if((iCurCharPos=FindKeyWordsPos(sBuf,sData))<0) return 0;
	for(i=0;i<m_loadnode_num;i++){
		id=ReadInt(iCurCharPos,sData);
		dbuf1=ReadDouble(iCurCharPos,sData);
		dbuf2=ReadDouble(iCurCharPos,sData);
		((CNode*)(m_nodeArray.GetAt(id-1)))->SetLoadXY(dbuf1,dbuf2);
	}
	//------------------------------------------------------
	GetHalfBand();//要利用单元和节点中的信息，必须在它们生成之后调用
	m_aStiffMatrix=new double*[m_node_num*2];//刚度矩阵
	for(i=0;i<m_node_num*2;i++)
	{
		m_aStiffMatrix[i]=new double[m_halfband];
	}
	m_aLoadVector=new double[m_node_num*2];//荷载向量
	//读取数据成功，返回1
	return 1;
}
//======================================================================
int CStructure::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 CStructure::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 );
}
//================================================================
int CStructure::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;
}
//=================================================
int CStructure::GetHalfBand()
{
	int i,j,tempmax;
	CNode* pNode[3];
	CElement* pElement;
	m_halfband=0;
	tempmax=0;
	for(i=0;i<m_element_num;i++)
	{
		pElement=(CElement*)(m_elementArray.GetAt(i));
		for(j=0;j<3;j++)
		{
			pNode[j]=pElement->m_pNode[j];
		}
		tempmax=max(abs(pNode[0]->m_node_id-pNode[1]->m_node_id),
			max(abs(pNode[1]->m_node_id-pNode[2]->m_node_id),
			abs(pNode[2]->m_node_id-pNode[0]->m_node_id)));
		if(tempmax>m_halfband)m_halfband=tempmax;
	}
	m_halfband=2*(m_halfband+1);
	return m_halfband;
}
//==================================================
//注意：只有用new 生成的二维数组才能用**方式传递
void CStructure::GetStiffMatrix(double **stiffmatrix)
{
	int ie,i,j,ii,jj;
	int ih,il,idh,idl,l;
	double** elestiffmatrix;
	//-----------------------------------------
	elestiffmatrix=new double*[6];//刚度矩阵
	for(i=0;i<6;i++)
	{
		elestiffmatrix[i]=new double[6];
	}
	//-------------------------------------------
	CElement* tempElement;
	for(i=0;i<m_node_num*2;i++)
	{
		for(j=0;j<m_halfband;j++)
		{
			stiffmatrix[i][j]=0.0;
		}
	}
	//------------------------------
	for(ie=0;ie<m_element_num;ie++)
	{	
		tempElement=(CElement*)(m_elementArray.GetAt(ie));
		tempElement->GetElementStiffMatrix(elestiffmatrix);
		for(i=0;i<3;i++)
		{
			for(ii=0;ii<2;ii++)
			{
				ih=2*i+ii;
				idh=2*(tempElement->m_pNode[i]->m_node_id-1)+ii;
				for(j=0;j<3;j++)
				{
					for(jj=0;jj<2;jj++)
					{
						l=2*j+jj;
						il=2*(tempElement->m_pNode[j]->m_node_id-1)+jj;
						idl=il-idh;
						if(idl>=0)
						{
							stiffmatrix[idh][idl]=stiffmatrix[idh][idl]+elestiffmatrix[ih][l];
						}
					}
				}
			}
		}
	}
	//-----------删除生成的临时数组elementMatrix[][]---------
	for(i=0;i<6;i++)
	{
		delete elestiffmatrix[i];
	}
	delete elestiffmatrix;
}
//================================================================
void CStructure::GetLoadVector(double* loadvector)
{
	int i,ie;
	int id;
	CNode* tempNode;
	CElement* tempElement;
	double weight;
	for(i=0;i<m_node_num;i++)
	{
		loadvector[i]=0.0;
	}
	if(m_loadnode_num!=0)
	{
		for(i=0;i<m_node_num;i++)
		{
			tempNode=(CNode*)(m_nodeArray.GetAt(i));
			loadvector[2*i]=tempNode->m_loadinfo[0];
			loadvector[2*i+1]=tempNode->m_loadinfo[1];
		}
	}
	if(m_weight!=0.0)
	{
		for(ie=0;ie<m_element_num;ie++)
		{
			tempElement=(CElement*)(m_elementArray.GetAt(ie));
			weight=tempElement->m_area*(-1)*tempElement->m_weight*tempElement->m_thickness/3;
			for(i=0;i<3;i++)
			{
				id=tempElement->m_pNode[i]->m_node_id;
				loadvector[2*id-1]=loadvector[2*id-1]+weight;
			}
		}
	}
}
//==========================================================
void CStructure::AnalysisStructure()
{
	GetStiffMatrix(m_aStiffMatrix);	//生成刚度矩阵
	SetSupport(m_aStiffMatrix);		//引入支承条件
	GetLoadVector(m_aLoadVector);	//生成荷载向量
	SolveEquationGroup(m_aStiffMatrix,m_aLoadVector);//解方程组
	GetStress();//求得所要求的数据，并存到结点和单元中去，在视图中可以调出显示。
}
//==========================================================
void CStructure::SolveEquationGroup(double **matrixKE, double *vectorP)
{
	int i,j,k;
	int m,im,iw;
	int nj1,nj2;
	int k1,l,i1,jo;
	double coef;
	nj2=m_node_num*2;
	nj1=nj2-1;
	for(k=1;k<=nj1;k++)
	{
		if(nj2>k+m_halfband-1)
			im=k+m_halfband-1;
		else 
			im=nj2;
		k1=k+1;
		for(i=k1;i<=im;i++)
		{
			l=i-k+1;
			coef=matrixKE[k-1][l-1]/matrixKE[k-1][1-1];
			iw=m_halfband-l+1;
			for(j=1;j<=iw;j++)
			{
				m=j+i-k;
				matrixKE[i-1][j-1]=matrixKE[i-1][j-1]-coef*matrixKE[k-1][m-1];
				vectorP[i-1]=vectorP[i-1]-matrixKE[i-1][j-1]*vectorP[k-1];
			}
		}
	}
	vectorP[nj2-1]=vectorP[nj2-1]/matrixKE[nj2-1][1-1];
	for(i1=1;i1<=nj1;i1++)
	{
		i=nj2-i1;
		if(m_halfband>nj2-i+1)jo=nj2-i+1;
		else 
			jo=m_halfband;
		for(j=2;j<=jo;j++)
		{
			k=j+i-1;
			vectorP[i-1]=vectorP[i-1]-matrixKE[i-1][j-1]*vectorP[k-1];
		}
		vectorP[i-1]=vectorP[i-1]/matrixKE[i-1][1-1];
	}
}
//==========================================================
void CStructure::GetStress()
{
	int i;
	CElement* pEle;
	CNode* pNode;
	//---------------将位移存到节点中去--------------
	for(i=0;i<m_node_num;i++)
	{
		pNode=(CNode*)(m_nodeArray.GetAt(i));
		pNode->m_displace[0]=m_aLoadVector[2*i];
		pNode->m_displace[1]=m_aLoadVector[2*i+1];
	}
	//---------------分析单元的应力情况--------------
	for(i=0;i<m_element_num;i++)
	{
		pEle=(CElement*)(m_elementArray.GetAt(i));		
		pEle->GetStress();
	}
}

void CStructure::SetSupport(double **stiffmatrix)
{
	int i,j,k;
	int ii,jo;
	CNode* curNode;
	for(i=0;i<m_node_num;i++)
	{
		curNode=(CNode*)(m_nodeArray.GetAt(i));
		if(curNode->m_supportinfo[0]==1||curNode->m_supportinfo[1]==1)
		{
			for(j=0;j<2;j++)
			{
				if(curNode->m_supportinfo[j]==1)
				{
					ii=2*(curNode->m_node_id-1)+j;//确定整体行号
					stiffmatrix[ii][0]=1.0f;
					for(k=1;k<m_halfband;k++)
					{
						stiffmatrix[ii][k]=0.0f;
					}
					//-------------修改斜线上的元素-----------
					if(ii>m_halfband)jo=m_halfband;
					else jo=ii;
					for(k=1;k<jo;k++)
					{
						stiffmatrix[ii-k][k]=0.0f;
					}	
					m_aLoadVector[ii]=0.0f;//修改右端向量
				}
			}
		}
	}
}