element.cpp

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

// Element.cpp: implementation of the CElement class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "FiniteElement.h"
#include "Element.h"
#include "math.h"

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

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

CElement::CElement(int id)
{
	m_element_id=id;
	Initial();
}

CElement::CElement()
{
	m_element_id=0;
	Initial();
}

CElement::~CElement()
{

}

double CElement::GetArea()
{
	double xij,yij,xik,yik;
	xij=m_pNode[1]->m_nodexy[0]-m_pNode[0]->m_nodexy[0];
	yij=m_pNode[1]->m_nodexy[1]-m_pNode[0]->m_nodexy[1];
	xik=m_pNode[2]->m_nodexy[0]-m_pNode[0]->m_nodexy[0];
	yik=m_pNode[2]->m_nodexy[1]-m_pNode[0]->m_nodexy[1];
	m_area=(xij*yik-xik*yij)*0.5;
	return m_area;
}

void CElement::SetNode(CNode* pNode1, CNode* pNode2, CNode* pNode3)
{
	m_pNode[0]=pNode1;
	m_pNode[1]=pNode2;
	m_pNode[2]=pNode3;
	GetArea();
}

void CElement::SetMaterialInfo(double young,double possion,double thickness,double weight)
{
	m_young=young;
	m_possion=possion;
	m_thickness=thickness;
	m_weight=weight;
}

void CElement::Initial()
{
	m_young=0.0;
	m_possion=0.0;
	m_thickness=0.0;
	m_weight=0.0;
	//-----------------------------
	 m_sigmaX=0.0;			//X应力分量
	 m_sigmaY=0.0;			//Y应力分量
	 m_taoXY=0.0;			//剪切应力分量
	 m_averageStress=0.0;	//平均应力
	 m_radiusStress=0.0;	//应力圆半径
	 m_maxMainStress=0.0;	//最大主应力
	 m_minMainStress=0.0;	//最小主应力
	 m_mainPlaneAngle=0.0;	//主平面角	

}
//===========================================================
void CElement::GetMatrixD(double** matrixD)
{
	int i,j;
	for(i=0;i<3;i++)
	{
		for(j=0;j<3;j++)
		{
			matrixD[i][j]=0.0;
		}
	}
	matrixD[0][0]=m_young/(1.0-m_possion*m_possion);
	matrixD[0][1]=m_possion*matrixD[0][0];
	matrixD[1][0]=matrixD[0][1];
	matrixD[1][1]=matrixD[0][0];
	matrixD[2][2]=0.5*m_young/(1.0+m_possion);
}
//============================================================
void CElement::GetMatrixB(double** matrixB)
{
	int i,j;
	for(i=0;i<3;i++)
	{
		for(j=0;j<6;j++)
		{
			matrixB[i][j]=0.0;
		}
	}
	matrixB[0][0]=m_pNode[1]->m_nodexy[1]-m_pNode[2]->m_nodexy[1];
	matrixB[0][2]=m_pNode[2]->m_nodexy[1]-m_pNode[0]->m_nodexy[1];
	matrixB[0][4]=m_pNode[0]->m_nodexy[1]-m_pNode[1]->m_nodexy[1];
	matrixB[1][1]=m_pNode[2]->m_nodexy[0]-m_pNode[1]->m_nodexy[0];
	matrixB[1][3]=m_pNode[0]->m_nodexy[0]-m_pNode[2]->m_nodexy[0];
	matrixB[1][5]=m_pNode[1]->m_nodexy[0]-m_pNode[0]->m_nodexy[0];
	matrixB[2][0]=matrixB[1][1];matrixB[2][1]=matrixB[0][0];
	matrixB[2][2]=matrixB[1][3];matrixB[2][3]=matrixB[0][2];
	matrixB[2][4]=matrixB[1][5];matrixB[2][5]=matrixB[0][4];
	for(i=0;i<3;i++)
	{
		for(j=0;j<6;j++)
		{
			matrixB[i][j]=0.5/m_area*matrixB[i][j];
		}
	}
}
//===================================================================
void CElement::GetMatrixS(double **matrixS)
{
	int i,j,k;
	double** matrixD;
	double** matrixB;
	//-----------------------------
	matrixD=new double*[3];//开辟内存区
	for(i=0;i<3;i++)
	{
		matrixD[i]=new double[3];
	}
	matrixB=new double*[3];//
	for(i=0;i<3;i++)
	{
		matrixB[i]=new double[6];
	}
	//-----------------------------
	GetMatrixD((double**)matrixD);
	GetMatrixB((double**)matrixB);
	for(i=0;i<3;i++)
	{
		for(j=0;j<6;j++)
		{
			matrixS[i][j]=0.0;
			for(k=0;k<3;k++)
			{
				matrixS[i][j]=matrixS[i][j]+matrixD[i][k]*matrixB[k][j];
			}
		}
	}

	//----------delete new data---------------
	for(i=0;i<3;i++)
	{
		delete matrixD[i];
	}
	delete matrixD;
	for(i=0;i<3;i++)
	{
		delete matrixB[i];
	}
	delete matrixB;
}
//======================================================================
void CElement::GetElementStiffMatrix(double **elementStiffMatrix)
{
	int i,j,k;
	double** matrixS;
	double** matrixB;
	//-----------------------------
	matrixS=new double*[3];//开辟内存区
	for(i=0;i<3;i++)
	{
		matrixS[i]=new double[6];
	}
	matrixB=new double*[3];//
	for(i=0;i<3;i++)
	{
		matrixB[i]=new double[6];
	}
	//----------------------------
	//----------------------------
	GetMatrixS(matrixS);
	GetMatrixB(matrixB);
	for(i=0;i<6;i++)
	{
		for(j=0;j<6;j++)
		{
			elementStiffMatrix[i][j]=0.0;
			for(k=0;k<3;k++)
			{
				elementStiffMatrix[i][j]=elementStiffMatrix[i][j]+
					matrixS[k][i]*matrixB[k][j]*m_area*m_thickness;
			}
		}
	}
	//----------delete new data---------------
	for(i=0;i<3;i++)
	{
		delete matrixS[i];
	}
	delete matrixS;
	for(i=0;i<3;i++)
	{
		delete matrixB[i];
	}
	delete matrixB;

}
//================================================================= 
void CElement::GetStress()
{
	int i,j;
	double** matrixS;
	double temp[3]={0.0,0.0,0.0};
	double disp[6];
	//----------------------
	matrixS=new double*[3];//开辟内存区
	for(i=0;i<3;i++)
	{
		matrixS[i]=new double[6];
	}
	//----------------------
	for(i=0;i<3;i++)
	{
		disp[i*2]=m_pNode[i]->m_displace[0];
		disp[i*2+1]=m_pNode[i]->m_displace[1];;
	}
	GetMatrixS(matrixS);
	for(i=0;i<3;i++)
	{
		for(j=0;j<6;j++)
		{
			temp[i]=temp[i]+matrixS[i][j]*disp[j];
		}
	}
	//----------------------------------
	m_sigmaX=temp[0];
	m_sigmaY=temp[1];
	m_taoXY=temp[2];
	m_averageStress=(m_sigmaX+m_sigmaY)/2;
	double dbuf;
	dbuf=(m_sigmaX-m_sigmaY)*(m_sigmaX-m_sigmaY)/4+m_taoXY*m_taoXY;
	m_radiusStress=sqrt(dbuf);
	m_maxMainStress=m_averageStress+m_radiusStress;
	m_minMainStress=m_averageStress-m_radiusStress;
	if(m_sigmaY==m_minMainStress)
	{
		m_mainPlaneAngle=0.0;
	}
	else 
	{
		m_mainPlaneAngle=90-57.29578*atan(m_taoXY/(m_sigmaY-m_minMainStress));
	}
	//----------------------
	for(i=0;i<3;i++)		//回收内存区
	{
		delete matrixS[i];
	}
	delete matrixS;
}