femview.cpp

用vc写的有限元计算代码

			{
				CFemUnit femUnit = lstUnit.GetNext(pos);
				nCoefMatrix[i*nCol+j] += femUnit.GetNodeCoefficient( i+1, j+1 );
			}
			
		}
	}

#if _DEBUG
	CFile file;
	file.Open(_T("D:\\FemCoefMatrixTest.log"), CFile::modeCreate | CFile::modeWrite );
	CString str;
	char cMat[4096] = {0};
	char cTmp[20] = {0};
	int nIdx = 0;
	for(int x=0; x<nRow; x++)
	{
		nIdx = 0;
		for(int y = 0; y<nCol; y++)
		{
			sprintf( cTmp, "%0.6f ", nCoefMatrix[x*nCol + y] );
			sprintf( cMat + nIdx, "%s", cTmp );
			nIdx += strlen( cTmp );
			if( nIdx > sizeof( cMat ) - 20 )
			{
				file.Write( cMat, nIdx);
				nIdx = 0;
			}
		}

		file.Write( cMat, nIdx);
		file.Write("\r\n", 2);
	}
	file.Close();
#endif

	memset( nCoefMatrix, 0, sizeof(double) * nRow * nCol );
	
	// 生成全局系数矩阵
	for(int i=0; i<nRow; i++)
	{
		if( femNode[i].IsFitValue() )
		{
			// 加入定解条件
			nCoefMatrix[i*nCol+i] = 1.0;
			nCoefMatrix[i*nCol+nRow] = femNode[i].GetNodeValue();
		}
		else
		{
			// 生成叠加节点系数矩阵
			for(int j=0; j<nRow; j++)
			{
				pos = lstUnit.GetHeadPosition();
				while( pos )
				{
					CFemUnit femUnit = lstUnit.GetNext(pos);
					nCoefMatrix[i*nCol+j] += femUnit.GetNodeCoefficient( i+1, j+1 );
				}
				
			}
		}
	}

	double *pResult = GaussMatrixSolve(nCoefMatrix, nRow);

	// 赋节点值
	for(int i=0; i<nRow; i++)
	{
		femNode[i].SetNodeValue(pResult[i]);
	}

	// 计算插值函数参数
	pos = lstUnit.GetHeadPosition();
	while( pos )
	{
		CFemUnit &femUnit = lstUnit.GetNext(pos);
		femUnit.MakeInterposeFuncCoef();
	}

	ShowFemGrid(lstUnit);

	delete nCoefMatrix;
}

// 使用高斯消去法求解矩阵方程
double * CFEMView::GaussMatrixSolve(double * pCoefMatrix, int nNodeNum)
{
	int r, c, w ;
	int nRow = nNodeNum;
    int nCol = nRow + 1;

    double nDiv;
    double *arrResult = new double[nCol];
	double *arrMatrix = pCoefMatrix;

    //变为上三角矩阵
    for( w = 0; w<nRow; w++)
	{
        for( r = w; r<nRow; r++)
		{
            nDiv = arrMatrix[r*nCol + w];
            if( nDiv == 0 ) continue;

            // 使左端的非零值变为1
			for( c = w; c < nCol; c++ ){
                arrMatrix[r*nCol + c] /= nDiv;
			}
		}

        // 使下三角部分的值都为0
        for( r = w + 1; r < nRow; r++)
		{
            if( arrMatrix[r*nCol + w] == 0 ) continue;

			for( c = w; c < nCol; c++){
                arrMatrix[r*nCol + c] -= arrMatrix[w*nCol + c];
			}
		}
	}

	for( c = 0; c < nCol; c++ ){
        arrResult[c] = 1;
	}

    // 取得最后一行的解
    arrResult[nRow - 1] = arrMatrix[(nRow - 1) * nCol + (nCol - 1)];

	// 计算其他的解
	for( w = nRow - 2; w >= 0; w-- ){
		for( c = w + 1; c < nCol; c++ ){
            arrMatrix[w * nCol + c] *= arrResult[c];
		}

		for( c = w + 1; c < nCol-1; c++ ){
            arrMatrix[w * nCol + (nCol - 1)] -= arrMatrix[w * nCol + c];
		}

        arrResult[w] = arrMatrix[w * nCol + (nCol - 1)];
	}

	memcpy(pCoefMatrix, arrResult, sizeof(double) * nCol);
	delete arrResult;

    return pCoefMatrix;
}

// 显示有限元网格
void CFEMView::ShowFemGrid(CList<CFemUnit> &lstFemUnit)
{
	CDC *pDC = this->GetDC();
	double nZoom = 1;//50;

	POSITION pos = lstFemUnit.GetHeadPosition();
	while( pos )
	{
		CFemUnit &femUnit = lstFemUnit.GetNext(pos);

		pDC->MoveTo(femUnit.GetNode(0)->GetX()*nZoom, femUnit.GetNode(0)->GetY()*nZoom);
		pDC->LineTo(femUnit.GetNode(1)->GetX()*nZoom, femUnit.GetNode(1)->GetY()*nZoom);
		pDC->LineTo(femUnit.GetNode(2)->GetX()*nZoom, femUnit.GetNode(2)->GetY()*nZoom);
		pDC->LineTo(femUnit.GetNode(0)->GetX()*nZoom, femUnit.GetNode(0)->GetY()*nZoom);

		femUnit.Show( pDC, nZoom );
		femUnit.DrawEqualLine( pDC, nZoom );
		femUnit.ShowNodeValue( pDC, nZoom );
	}

	this->ReleaseDC( pDC );
}

BOOL CFEMView::CalculateCurveParameter(CDoubleArray *X,CDoubleArray *Y,long M,long N,CDoubleArray *A)
{
	//X,Y --  X,Y两轴的坐标
	//M   --  结果变量组数
	//N   --  采样数目
	//A   --  结果参数

	register long i,j,k;
	double Z,D1,D2,C,P,G,Q;
	CDoubleArray B,T,S;
	B.SetSize(N);
	T.SetSize(N);
	S.SetSize(N);
	A->SetSize(M);

	if(M>N)M=N;
	for(i=0;i<M;i++)
		(*A)[i]=0;
	Z=0;
	B[0]=1;
	D1=N;
	P=0;
	C=0;
	for(i=0;i<N;i++)
	{
		P=P+(*X)[i]-Z;
		C=C+(*Y)[i];
	}
	C=C/D1;
	P=P /D1;
	(*A)[0]=C*B[0];
	if(M>1)
	{
		T[1]=1;
		T[0]=-P;
		D2=0;
		C=0;
		G=0;
		for(i=0;i<N;i++)
		{
			Q=(*X)[i]-Z-P;
			D2=D2+Q*Q;
			C=(*Y)[i]*Q+C;
			G=((*X)[i]-Z)*Q*Q+G;
		}
		C=C/D2;
		P=G/D2;
		Q=D2/D1;
		D1=D2;
		(*A)[1]=C*T[1];
		(*A)[0]=C*T[0]+(*A)[0];
	}
	for(j=2;j<M;j++)
	{
		S[j]=T[j-1];
		S[j-1]=-P*T[j-1]+T[j-2];
		if(j>=3)
		{
			for(k=j-2;k>=1;k--)
				S[k]=-P*T[k]+T[k-1]-Q*B[k];
		}
		S[0]=-P*T[0]-Q*B[0];
		D2=0;
		C=0;
		G=0;
		for(i=0;i<N ;i++)
		{
			Q=S[j];
			for(k=j-1;k>=0;k--)
				Q=Q*((*X)[i]-Z)+S[k];
			D2=D2+Q*Q;
			C=(*Y)[i]*Q+C;
			G=((*X)[i]-Z)*Q*Q+G;
		}
		C=C/D2;
		P=G/D2;
		Q=D2/D1;
		D1=D2;
		(*A)[j]=C*S[j];
		T[j]=S[j];
		for(k=j-1;k>=0;k--)
		{
			(*A)[k]=C*S[k]+(*A)[k];
			B[k]=T[k];
			T[k]=S[k];
		}
	}
	return TRUE;
}





void spline(double x[],double y[],int n,double yp1,double ypn,double y2[]) 
{ 
	double u[100],aaa,sig,p,bbb,ccc,qn,un; 
	int i,k; 
	if(yp1 > 9.9e+29) 
	{ 
		y2[1] = 0; 
		u[1] = 0; 
	} 
	else 
	{ 
		y2[1] = - 0.5; 
		aaa = (y[2] - y[1])/(x[2] - x[1]); 
		u[1] = (3.0/(x[2] - x[1]))*(aaa- yp1);
	} 
	for(i = 2;i<=n-1;i++) 
	{ 
		sig = (x[i] - x[i-1])/(x[i+1] - x[i-1]); 
		p = sig*y2[i-1] + 2.0; 
		y2[i] = (sig-1.0)/p; 
		aaa = (y[i+1] - y[i])/(x[i+1]-x[i-1]); 
		bbb = (y[i] - y[i-1])/(x[i]-x[i-1]); 
		ccc = x[i+1] - x[i-1]; 
		u[i] = (6.0*(aaa-bbb)/ccc-sig*u[i-1])/p; 
	} 
	if(ypn>9.9e+29) 
	{ 
		qn = 0.0; 
		un = 0.0; 
	} 
	else 
	{ 
		qn = 0.5; 
		aaa = ypn - (y[n] - y[n-1])/(x[n]-x[n-1]); 
		un = (3.0/(x[n] - x[n-1]))*aaa; 
	} 
	y2[n] =(un-qn*u[n-1])/(qn*y2[n-1] + 1.0); 
	for(k = n-1;k>=1;k--) 
		y2[k] = y2[k]*y2[k+1] + u[k]; 
} 


void splint(double xa[],double ya[],double y2a[],int n,double& x,double& y) 
{ 
	int klo,khi,k; 
	double h,a,b,aaa,bbb; 
	klo = 1; 
	khi = n; 
loop:
	if(khi-klo>1) 
	{ 
	 k = (khi+klo)/2; 
	 if(xa[k]>x) 
		 khi = k; 
	 else 
	 { 
		 klo = k; 
	 } 
	 goto loop; 
	} 
	h = xa[khi] - xa[klo]; 
	if(h==0) 
	{ 
	 //cout<<"pause 'bad xa input'"<<endl; 
	 return; 
	} 
	a = (xa[khi] - x)/h; 
	b = (x-xa[klo])/h; 
	aaa = a*ya[klo] + b*ya[khi]; 
	bbb = (a*a*a - a)*y2a[klo] + (b*b*b - b)*y2a[khi]; 
	y = aaa + bbb*h*h/6.0; 
} 


void CFEMView::OnFemBezier()
{
	CClientDC dc( this );
	/*CPoint *point = new CPoint[10];
	point[0] = CPoint(10, 10);
	point[1] = CPoint(20, 12);
	point[2] = CPoint(25, 18);
	point[3] = CPoint(35, 22);
	point[4] = CPoint(40, 38);
	point[5] = CPoint(50, 42);
	point[6] = CPoint(55, 48);
	point[7] = CPoint(75, 60);
	point[8] = CPoint(80, 55);
	point[9] = CPoint(90, 40);

	

	for(int i=0; i<10; i++)
	{
		point[i].x *= 5;
		point[i].y *= 5;

		dc.Ellipse(point[i].x - 3, point[i].y - 3, point[i].x + 3, point[i].y + 3);
	}

	
	dc.PolyBezier(point, 10);*/

	int   px[12]={50,80,120,140,180,230,270,380,430,460,500,550};//x坐标,10个点！   
	int   py[12]={100,230,230,160,50,50,120,230,150,170,150,190};//y坐标   

	for(int i=0; i<12; i++)
	{
		dc.Ellipse(px[i] - 3, py[i] - 3, px[i] + 3, py[i] + 3);
		Sleep(50);
	}

	int   *x = px;
	int   *y = py;
	int   n = 10;
	unsigned   int   k = 20;    
      
	unsigned   int   i,j;     
	float   t1,t2,t3,t,a,b,c,d,tx,ty;  

	/**x = *(x+1);
	*y = *(y+1);     
	*(x+n+1)=*(x+n);
	*(y+n+1)=*(y+n); */

	t=0.5/k;     

	CPen pen;
	pen.CreatePen(PS_SOLID, 1, RGB(255,0,0));
	dc.SelectObject(&pen);
  
	dc.MoveTo(*(x+1),*(y+1));
	for(i=0;i<n-1;i++)     
	{     
		for(j=1;j<k;j++)     
		{     
			t1=j*t;     
			t2=t1*t1;     
			t3=t2*t1;     

			a = 4 * t2 - t1 - 4*t3;     
			b = 1- 10*t2 + 12*t3;     
			c = t1 + 8*t2 - 12*t3;     
			d = 4*t3 - 2*t2;    

			tx=a*(*(x+i))+b*(*(x+i+1))+c*(*(x+i+2))+d*(*(x+i+3));     
			ty=a*(*(y+i))+b*(*(y+i+1))+c*(*(y+i+2))+d*(*(y+i+3));

			dc.LineTo(tx,ty);     
			Sleep(20);
		}
	}     
	dc.LineTo(*(x+i+2),*(y+i+2));  


	//unsigned   int   i,j;     
	//float   t,t1,t2,a,b,c,tx,ty;     
	/**x=*(x+1);
	*y=*(y+1);     
	*(x+n+1)=*(x+n);
	*(y+n+1)=*(y+n);  */

	t=1.0/k;      

	pen.DeleteObject();
	pen.CreatePen(PS_SOLID, 1, RGB(0,0,255));
	dc.SelectObject(&pen);

	dc.MoveTo((*x+(*(x+1)))/2.0,(*y+(*(y+1)))/2.0);     
	for(i=0;i<n;i++)     
	{     
		for(j=1;j<k;j++)     
		{     
			t1=j*t;     
			t2=t1*t1;     
			a=(t2-2*t1+1)/2.0;     
			b=t1-t2+1/2.0;     
			c=t2/2.0;     
			tx=a*(*(x+i))+b*(*(x+i+1))+c*(*(x+i+2));     
			ty=a*(*(y+i))+b*(*(y+i+1))+c*(*(y+i+2));     
			dc.LineTo(tx,ty);   
			Sleep(20);
		}     
	}

}