structure.cpp

用VC写得计算电磁学的有限元分析程序

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

#include "stdafx.h"
#include "FiniteElec.h"
#include "Structure.h"
#include "SolveEquation.h"

#include "Node.h"
#include "Element.h"

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

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

CStructure::CStructure()
{
	Init();
}

CStructure::~CStructure()
{
	Init();
}

void CStructure::Init()
{
	m_strExamName=_T("无数据");
	m_strStructureData=_T("");
	m_nodecount=0;
	m_elementcount=0;
	m_halfband=0;
	//----------------
	m_done=0;
	m_solved=0;	//未求解
	m_aLoadVector=NULL;
	m_aStiffMatrix=NULL;
	//----------------
	CNode* pNode;
	CElement* pElement;
	m_model.Init();
//	if(m_model.m_nodecount==0)return;//已删除
	//删除结点对象
	for(int i=0;i<m_nodeArray.GetSize();i++)
	{
		pNode=(CNode*)(m_nodeArray.GetAt(i));
		delete pNode;
		m_nodeArray.RemoveAll();
	}
	//删除单元对象
	for(i=0;i<m_elementArray.GetSize();i++)
	{
		pElement=(CElement*)(m_elementArray.GetAt(i));
		delete pElement;
		m_elementArray.RemoveAll();
	}
	//初始化方程组
	InitEquation(0);//nodecount=0
}
///////////////////////////////////////
//算半带宽
int CStructure::GetHalfBand()
{
	if(m_halfband!=0)return m_halfband;//已求

	int i,j,tempmax;
	CNode* pNode[3];
	CElement* pElement;
	m_halfband=0;
	tempmax=0;
	for(i=0;i<m_elementcount;i++)
	{
		pElement=(CElement*)(m_elementArray.GetAt(i));
		for(j=0;j<3;j++)
		{
			if((pNode[j]=pElement->m_pNodes[j])==NULL)return 0;
		}
		tempmax=max(abs(pNode[0]->m_nodeid-pNode[1]->m_nodeid),
			max(abs(pNode[1]->m_nodeid-pNode[2]->m_nodeid),
			abs(pNode[2]->m_nodeid-pNode[0]->m_nodeid)));
		if(tempmax>m_halfband)m_halfband=tempmax;
	}
	m_halfband=(m_halfband+1);
	return m_halfband;
}
////////////////////////////////////////
//求总刚度矩阵
void CStructure::GetStiffMatrix(double **stiffmatrix)
{
	int i,j,k;//循环变量
	int iz,jz;//单刚的元素在总刚中编号
	int ib,jb;//单刚的元素在半带中的编号
	double** elestiffmatrix;
	CElement* pElement;
	//-----------------------------------------
	elestiffmatrix=new double*[3];//临时单元刚度矩阵
	for(i=0;i<3;i++)
	{
		elestiffmatrix[i]=new double[3];
	}
	//==================================
	//==对各个单元扫描一遍-------
	for(i=0;i<m_elementcount;i++)
	{
		pElement=(CElement*)(m_elementArray.GetAt(i));
		pElement->GetElementStiffMatrix(elestiffmatrix);
		//对单刚进行扫描
		for(j=0;j<3;j++)
		{
			for(k=0;k<3;k++)
			{
				//单元编号第一位是1,而数组是0
				iz=pElement->m_pNodes[j]->m_nodeid-1;
				jz=pElement->m_pNodes[k]->m_nodeid-1;
				m_aStiffMatrixAll[iz][jz]=m_aStiffMatrixAll[iz][jz]+elestiffmatrix[j][k];
				if(iz>jz)continue;//下半三角不用存储
				ib=iz;//半带宽中行号不变
				jb=jz-ib;//半带宽中列号
				stiffmatrix[ib][jb]=stiffmatrix[ib][jb]+elestiffmatrix[j][k];
			}
		}	
	}
	//==================================
	//-----------删除生成的临时数组elementMatrix[][]
	for(i=0;i<3;i++)
	{
		delete elestiffmatrix[i];
	}
	delete elestiffmatrix;
	///////////////////////
/*	for(i=0;i<m_nodecount;i++)
	{
		for(j=0;j<m_halfband;j++)
		{
			m_tempArray[i][j]=m_aStiffMatrix[i][j];
		}
	}*/

}
//////////////////////
//生成荷载向量
void CStructure::GetLoadVector(double *loadvector)
{
	int i;
	CNode* pNode;
	//结各个结点扫描一遍
	for(i=0;i<m_nodecount;i++)
	{
		pNode=(CNode*)(m_nodeArray.GetAt(i));
		loadvector[i]=pNode->m_nodeValue;
	}
}
//////////////////////
//引入边界条件，相当于固定支撑
void CStructure::SetSupport(double **stiffmatrix)
{
	//用于处理边界电位为0的情况
	int i;
	int j;
	CNode* pNode;
	//结各个结点扫描一遍
	for(i=0;i<m_nodecount;i++)
	{
		pNode=(CNode*)(m_nodeArray.GetAt(i));
		//如果为2，是零电位点，对方程进行处理。
		if(pNode->m_nodetype==0)
		{
			//无载荷电压
			
		}
		else if(pNode->m_nodetype==1)
		{
			//有载荷电压
			//---------全刚度矩阵---------
			for(j=0;j<m_nodecount;j++)
			{
				if(i!=j)
				{
					m_aStiffMatrixAll[i][j]=0.0f;
					//m_aStiffMatrixAll[j][i]=0.0f;
				}
				else
				{//对角线元素化为1
					m_aStiffMatrixAll[i][j]=1.0f;
				}
			}
			//---------半带宽--------------
			stiffmatrix[i][0]=1.0;
			for(j=1;j<m_halfband;j++)
			{
				stiffmatrix[i][j]=0.0f;
			}
			//斜线上元素化零
		/*	for(j=0;j<((i<m_halfband)?i:m_halfband);j++)
			{
				if(i==0)break;
				stiffmatrix[i-j-1][j+1]=0.0f;
			}
		*/	//-----------------------------
		}
		else if(pNode->m_nodetype>100)
		{//重复结点处理
			m_aStiffMatrixAll[i][i]=1.0f;
			stiffmatrix[i][0]=1.0;
			m_aLoadVector[i]=0.0f;
		}
	}
}
////////////////////////////////////////
//求解半带宽方程组
int 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_nodecount;
	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];
	}
	return 1;
}

//////////////////////
//===================================
//==进行有限元分析-主运行模块-主线
int CStructure::AnalysisStructure()
{
	if(m_solved)return 1;	//已求解
	if(m_done!=1)return 0;//数据体未准备好
	if(GetHalfBand()==0)return 0;		//获取半带宽
	InitEquation(1);		//初始化方程组
	GetStiffMatrix(m_aStiffMatrix);	//生成刚度矩阵
	GetLoadVector(m_aLoadVector);	//生成荷载向量
	SetSupport(m_aStiffMatrix);		//引入支承条件
	//SolveEquationGroup(m_aStiffMatrix,m_aLoadVector);//半带宽消元法
	CSolveEquation::Gauss(m_aStiffMatrixAll,m_aLoadVector,m_nodecount);
	//保存结果电压值到结点中去。
	int i;
	CNode* pNode;
	for(i=0;i<m_nodecount;i++)
	{
		pNode=(CNode*)(m_nodeArray.GetAt(i));
		pNode->m_nodeValue=m_aLoadVector[i];
	}
	m_solved=1;	//成功解决问题
	return 1;
}
//////////////////////////////////////
///
void CStructure::Serialize(CArchive &ar)
{
	int loop;
	int i,j;
	long nChar;
	char cBuf;
	CNode* pNode;
	CString temstr;
	if(ar.IsStoring())
	{
		if(!m_solved)//没有结算，返回
		{
	//		AfxMessageBox("还没有进行求解，不能输出结果。请先求解!　");
			return;
		}
		//保存数据
		///////////////////////////////////////////
		ar.WriteString("+++++++++++++++++++++++\r\n");