femunit.cpp

用vc写的有限元计算代码

//
// 有限元单元类
//
// 黄玉洲 2008.08.01

#include "StdAfx.h"
#include "FemUnit.h"
#include "FEM.h"
#include <math.h>

CFemUnit::CFemUnit(void)
{
	m_nUnitNo = 0;
	m_nDielectric = 1;
}

CFemUnit::CFemUnit(int nUnitNo, CFemNode *pFemNode1, CFemNode *pFemNode2, CFemNode *pFemNode3)
{
	this->m_pNode[0] = pFemNode1;
	this->m_pNode[1] = pFemNode2;
	this->m_pNode[2] = pFemNode3;
	m_nUnitNo = nUnitNo;
	m_nDielectric = 1;
	MakeCoefficientMatrix();
}

CFemUnit::CFemUnit(const CFemUnit & femUnit)
{
	*this = femUnit;
}

CFemUnit::~CFemUnit(void)
{
	//delete[] m_pNode;
}

CFemUnit & CFemUnit:: operator =(const CFemUnit &femUnit)
{
	/*memcpy( m_pNode, femUnit.m_pNode, sizeof(m_pNode) );
	memcpy( m_CoefMatrix, femUnit.m_CoefMatrix, sizeof(m_CoefMatrix) );
	m_nUnitNo = femUnit.m_nUnitNo;*/

	// 在没有智能指针的情况下可以做以下操作
	memcpy(this, &femUnit, sizeof(CFemUnit));

	return *this;
}


// 生成系数矩阵
void CFemUnit::MakeCoefficientMatrix(void)
{
	double x1 = m_pNode[0]->GetX();
	double x2 = m_pNode[1]->GetX();
	double x3 = m_pNode[2]->GetX();
	double y1 = m_pNode[0]->GetY();
	double y2 = m_pNode[1]->GetY();
	double y3 = m_pNode[2]->GetY();

	double A = ( (x2-x1)*(y3-y1) - (x3-x1)*(y2-y1) ) / 2.0;
	double nCoef = 2.2 / (m_nDielectric * ( 4*A ));

	m_CoefMatrix[0][0] = ( (y2-y3)*(y2-y3) + (x2-x3)*(x2-x3) ) * nCoef;
	m_CoefMatrix[0][1] = ( (y2-y3)*(y3-y1) + (x2-x3)*(x3-x1) ) * nCoef;
	m_CoefMatrix[0][2] = ( (y2-y3)*(y1-y2) + (x2-x3)*(x1-x2) ) * nCoef;
	m_CoefMatrix[1][0] = m_CoefMatrix[0][1];
	m_CoefMatrix[1][1] = ( (y3-y1)*(y3-y1) + (x3-x1)*(x3-x1) ) * nCoef;
	m_CoefMatrix[1][2] = ( (y3-y1)*(y1-y2) + (x3-x1)*(x1-x2) ) * nCoef;
	m_CoefMatrix[2][0] = m_CoefMatrix[0][2];
	m_CoefMatrix[2][1] = m_CoefMatrix[1][2];
	m_CoefMatrix[2][2] = ( (y1-y2)*(y1-y2) + (x1-x2)*(x1-x2) ) * nCoef;
}

// 取得节点的系数
double CFemUnit::GetNodeCoefficient(int nNodeIndex1, int nNodeIndex2)
{
	int nRealNode1 = -1;
	int nRealNode2 = -1;
	double nNodeCoef = 0.0;

	// 根据节点编号判断节点的索引值
	for(int i=0; i<3; i++)
	{
		if( m_pNode[i]->GetNodeNo() == nNodeIndex1 )
			nRealNode1 = i;
		if( m_pNode[i]->GetNodeNo() == nNodeIndex2 )
			nRealNode2 = i;
	}

	// 如果节点和本单元相关，则返回节点的矩阵系数
	if( nRealNode1 > -1 && nRealNode2 > -1)
	{
		nNodeCoef = m_CoefMatrix[nRealNode1][nRealNode2];
	}

	return nNodeCoef;
}

// 显示
void CFemUnit::Show(CDC* pDC, double nZoom)
{
	double nMaxX, nMinX;
	double nMaxY, nMinY;

	double nMaxVal = 100;
	double nMinVal = 0;

	nMinX = nMaxX = m_pNode[0]->GetX();
	nMinY = nMaxY = m_pNode[0]->GetY();
	for(int i=1; i<3; i++)
	{
		if( m_pNode[i]->GetX() > nMaxX) 
			nMaxX = m_pNode[i]->GetX();
		else if( m_pNode[i]->GetX() < nMinX) 
			nMinX = m_pNode[i]->GetX();

		if( m_pNode[i]->GetY() > nMaxY) 
			nMaxY = m_pNode[i]->GetY();
		else if( m_pNode[i]->GetY() < nMinY) 
			nMinY = m_pNode[i]->GetY();
	}

	nMinX *= nZoom;
	nMaxX *= nZoom;
	nMinY *= nZoom;
	nMaxY *= nZoom;

	for(double x=nMinX; x <= nMaxX; x++)
		for(double y=nMinY; y<= nMaxY; y++)
		{
			if( IsInUnit(x/nZoom, y/nZoom) )
			{
				double nVal = GetPosValue(x/nZoom, y/nZoom);

				double nClrG = 255.0 * nVal / (nMaxVal - nMinVal);
				double nClrB = 255.0 - 255.0 * nVal / (nMaxVal - nMinVal);

				pDC->SetPixel(CPoint((int)x, (int)y), RGB( 0, nClrG,nClrB));
			}
		}

	
}
// 显示电压值
void CFemUnit::ShowNodeValue(CDC* pDC, double nZoom)
{
	// 显示电压值
	CString str;
	str.Format(_T("%0.3f"), m_pNode[0]->GetNodeValue() );
	pDC->TextOut( m_pNode[0]->GetX()*nZoom, m_pNode[0]->GetY()*nZoom, str);

	str.Format(_T("%0.3f"), m_pNode[1]->GetNodeValue() );
	pDC->TextOut( m_pNode[1]->GetX()*nZoom, m_pNode[1]->GetY()*nZoom, str);

	str.Format(_T("%0.3f"), m_pNode[2]->GetNodeValue() );
	pDC->TextOut( m_pNode[2]->GetX()*nZoom, m_pNode[2]->GetY()*nZoom, str);
}

// 判断节点是否在三角形单元内
bool CFemUnit::IsInUnit(double nX, double nY)
{
	double x1 = m_pNode[0]->GetX();
	double x2 = m_pNode[1]->GetX();
	double x3 = m_pNode[2]->GetX();
	double x4 = nX;
	double y1 = m_pNode[0]->GetY();
	double y2 = m_pNode[1]->GetY();
	double y3 = m_pNode[2]->GetY();
	double y4 = nY;

	// 如果是逆时针旋转的向量,向量的积分值为正,否则为负;
	// 点在三角行内的必要条件是, 第一象限, 点和三条边分别形成的向量组的积分都为正值.

	if( (y1+y2) * (x1-x2) + (y2+y4) * (x2-x4) + (y4+y1) * (x4-x1) < 0 )
		return false;

	if( (y2+y3) * (x2-x3) + (y3+y4) * (x3-x4) + (y4+y2) * (x4-x2) < 0 )
		return false;

	if( (y3+y1) * (x3-x1) + (y1+y4) * (x1-x4) + (y4+y3) * (x4-x3) < 0 )
		return false;

	return true;
}

// 取得三角形单元内的插值
double CFemUnit::GetPosValue(double nX, double nY)
{
	// 插值函数 V = a + b*x + c*y
	return m_nCoefA + m_nCoefB * nX + m_nCoefC * nY;
}

// 计算插值函数的参数
void CFemUnit::MakeInterposeFuncCoef(void)
{
	double x1 = m_pNode[0]->GetX();
	double x2 = m_pNode[1]->GetX();
	double x3 = m_pNode[2]->GetX();
	double y1 = m_pNode[0]->GetY();
	double y2 = m_pNode[1]->GetY();
	double y3 = m_pNode[2]->GetY();
	double v1 = m_pNode[0]->GetNodeValue();
	double v2 = m_pNode[1]->GetNodeValue();
	double v3 = m_pNode[2]->GetNodeValue();

	double A = ( (x2-x1)*(y3-y1) - (x3-x1)*(y2-y1) ) / 2.0;

	m_nCoefA = ( (x2*y3-x3*y2) * v1 + (x3*y1-x1*y3) * v2 + (x1*y2-x2*y1) * v3 ) / (2 * A);
	m_nCoefB = ( (y2-y3) * v1 + (y3-y1) * v2 + (y1-y2) * v3 ) / (2 * A);
	m_nCoefC = ( (x3-x2) * v1 + (x1-x3) * v2 + (x2-x1) * v3 ) / (2 * A);

}

// 画等位线
void CFemUnit::DrawEqualLine(CDC* pDC, double nZoom)
{
	//if( m_nUnitNo != 5 ) return;

	CList<CFemNode> lstVolPnt;

	double nMaxVal = 100;
	double nMinVal = 0;
	double nStep = 2;

	CFemNode *pNode1;
	CFemNode *pNode2;

	int nIdxS = 0;
	int nIdxE = 1;

	// 检查分布到三条边上的所有电位,找出符合需要输出的电位值,并计算该点的坐标.
	do
	{
		if( nIdxE > 2 ) nIdxE = 0;
		if( m_pNode[nIdxS]->GetNodeValue() < m_pNode[nIdxE]->GetNodeValue() )
		{
			pNode1 = m_pNode[nIdxS];
			pNode2 = m_pNode[nIdxE];
		}
		else
		{
			pNode1 = m_pNode[nIdxE];
			pNode2 = m_pNode[nIdxS];
		}

		nIdxS++;
		nIdxE++;

		for(double i=nMinVal; i<=nMaxVal; i+=nStep)
		{
			if(pNode1->GetNodeValue() <= i && pNode2->GetNodeValue() >= i)
			{
				double x1 = pNode1->GetX();
				double x2 = pNode2->GetX();
				double y1 = pNode1->GetY();
				double y2 = pNode2->GetY();

				// 起始点到终止点的距离
				double a = y2 - y1;
				double b = x2 - x1;
				double l = sqrt( a*a + b*b );

				// 等值点到起始点的距离
				double lx = ( i - pNode1->GetNodeValue() ) * l / ( pNode2->GetNodeValue() - pNode1->GetNodeValue() );

				// 等值点为起始点的判断
				if( l == 0 ){
					lstVolPnt.AddTail( CFemNode(x1, y1, i) );
					continue;
				}

				// 等值点为终止点的判断
				if( l - lx > -0.000000001 && l - lx < 0.000000001  ){
					lstVolPnt.AddTail( CFemNode(x2, y2, i) );
					continue;
				}

				// 斜角度
				double r = asin( a / l );

				// 角度调整
				if( a <= 0 && b < 0)
				{
					//第四象限
					r = -3.1415926 - r;
				}
				else if( a >= 0 && b < 0)
				{
					//第二象限
					r = 3.1415926 - r;
				}

				// 等值点的坐标
				double x = x1 + lx * cos( r );
				double y = y1 + lx * sin( r );

				lstVolPnt.AddTail( CFemNode(x, y, i) );

			}
		}
	}while(nIdxE > 1);

	// 画线
	POSITION pos = lstVolPnt.GetHeadPosition();
	while( pos )
	{
		pNode1 = &lstVolPnt.GetNext( pos );

		// 查找除自身之外的相同电位的点
		POSITION npos = pos;
		while( npos )
		{
			pNode2 = &lstVolPnt.GetNext( npos );

			double dif = pNode1->GetNodeValue() - pNode2->GetNodeValue();
			if( dif >= -0.00000001 && dif <= 0.00000001 )
			{
				pDC->MoveTo(pNode1->GetX() * nZoom, pNode1->GetY() * nZoom );
				pDC->LineTo(pNode2->GetX() * nZoom, pNode2->GetY() * nZoom );

				break;
			}
		}
	}

}