femath.h

线性方程组求解。ICCG法；LU分解法； 东华理工 张志勇

#pragma once
#include "math.h"
#include "FeMarix.h"
/*有限单元常用函数*/
//贝塞尔函数
#ifndef MYMATH
#define MYMATH
double bessi0(double x)
{
	double ax,ans;
	double y;
	if((ax = fabs(x)) <3.75)
	{
		y = x/3.75;
		y*=y;
		ans = 1.0 + y*(3.5156229 + y*(3.0899424 + y*(1.2067492 + 
			y*(0.2659732 + y*(0.360768e-1 + y*0.45813e-2)))));
	}
	else
	{
		y = 3.75/ax;
		ans = (exp(ax)/sqrt(ax)*(0.39894228 + y*(0.1328592e-1 
			+ y*(0.225319e-2 + y*(-0.157565e-2 + y*(0.916281e-2
			+ y*(-0.2057706e-1 + y*(0.2635537e-1 + y*(-0.1647633e-1
			+ y*(0.392377e-2))))))))));
	}
	return ans;
}
double bessk0(double x)
{
	double y,ans;
	if(x <= 2.0)
	{
		y = x*x/4.0;
		ans = (-log(x/2.0)*bessi0(x)) + (-0.57721566 + y*(0.42278420
			+ y*(0.23069756 + y*(0.3488590e-1 + y*(0.262698e-2
			+ y*(0.10750e-3 + y*0.74e-5))))));
	}
	else
	{
		y = 2.0/x;
		ans = (exp(-x)/sqrt(x)*(1.25331414 + y*(-0.7832358e-1
			+ y*(0.2189568e-1 + y*(-0.1062446e-1 + y*(0.587872e-2
			+ y*(-0.251540e-2 + y*0.53208e-3)))))));
	}
	return ans;
}

double bessi1(double x)
{
	double ax,ans;
	double y;
	if((ax = fabs(x)) < 3.75)
	{
		y = x/3.75;
		y *= y;
		ans = ax*(0.5 + y*(0.87890594 + y*(0.51498869 + y*(0.15084934
			+ y*(0.2658733e-1 + y*(0.301532e-2 + y*0.32411e-3))))));
	}
	else
	{
		y = 3.75/ax;
		ans =0.2282967e-1 + y*(-0.2895312e-1 + y*(0.1787654e-1
			- y*0.420059e-2));
		ans = 0.39894228 + y*(-0.3988024e-1 + y*(-0.362018e-2 
			+ y*(0.163801e-2 + y*(-0.1031555e-1 + y*ans))));
		ans *= (exp(ax)/sqrt(ax));
	}
	return x < 0.0 ? -ans : ans;
}
double bessk1(double x)
{
	double y,ans;
	if( x <= 2.0)
	{
		y = x*x/4.0;
		ans = (log(x/2.0)*bessi1(x) + (1.0/x)*(1.0 + y*(0.15443144
			+ y*(-0.67278579 + y*(-0.18156897 + y*(-0.1919402e-1
			+ y*(-0.110404e-2 + y*(-0.4686e-4))))))));
	}
	else
	{
		y = 2.0/x;
		ans = (exp(-x)/sqrt(x))*(1.25331414 + y*(0.234988619
			+ y*(-0.3655620e-1 + y*(0.1504268e-1 + y*(-0.780353e-2
			+ y*(0.325614e-2 + y*(-0.68245e-3)))))));
	}
	return ans;
}
//线性方程组列主消元
void SolveEquateLU(float** m_pfDGA, float* m_pfX, float* m_pfB, long m_lN)
{
	//列主消元解方程
	for(long i = 0; i < m_lN-1; i++)
	{
		//找主元并交换
		//EquateME( m_pfDGA,  m_pfB,  m_lN, i);
		//消元
		for(long j = i+1; j < m_lN; j++)
		{
			m_pfDGA[j][i] /= m_pfDGA[i][i];
			for(long k = i+1; k < m_lN; k++)
			{
				m_pfDGA[j][k] -= m_pfDGA[i][k]*m_pfDGA[j][i];
			}
			m_pfB[j] -= m_pfB[i]*m_pfDGA[j][i];
			m_pfDGA[j][i] = 0;
		}
	}
	//回带求解
	for(long i = m_lN-1; i >= 0; i--)
	{
		for(long j = i+1; j < m_lN; j++)
		{
			m_pfB[i] -= m_pfX[j]*m_pfDGA[i][j];
		}
		m_pfX[i] = m_pfB[i]/m_pfDGA[i][i];
	}
}

void EquateME(float** m_pfGA, float* m_pfB, long m_lN, long m_lNCur)
{
	float maxcur = fabs(m_pfGA[m_lNCur][m_lNCur]);
	long maxcuri = m_lNCur;
	for(long i = m_lNCur; i < m_lN; i++)
	{
		if(fabs(m_pfGA[i][m_lNCur]) > maxcur)
		{
			maxcur = fabs(m_pfGA[i][m_lNCur]);
			maxcuri = i;
		}
	}
	if(maxcuri != m_lNCur)
	{
		for(long i = m_lNCur; i < m_lN; i++)
		{
			maxcur = m_pfGA[m_lNCur][i];
			m_pfGA[m_lNCur][i] = m_pfGA[maxcuri][i];
			m_pfGA[maxcuri][i] = maxcur;
		}
		maxcur = m_pfB[m_lNCur];
		m_pfB[m_lNCur] = m_pfB[maxcuri];
		m_pfB[maxcuri] = maxcur;
	}
}
//ICCG
void SolveEqutionICCG(CFeMarix& A, CFeVector& B, CFeVector& X,long m_lTotal)
{
	CFeVector  R(m_lTotal);
	CFeVector  P(m_lTotal);
    CFeVector  T(m_lTotal);
	CFeVector  AP(m_lTotal);
    CFeVector  alfAP(m_lTotal);
	CFeVector  ICCR(m_lTotal);
	double alf = 0;
	double bt = 0;
	double eps = 0;
	long	k = 0;

	//ICCG
	R = B - A*X;
	A.ICC();
	P =  A.ICCINV(R);
	do
	{
		AP = A*P;		
		bt = P*AP;
		alf = R*A.ICCINV(R);
		alf /= bt;

		T   = P*alf;		
		X   = X + T;
		alfAP = AP*alf;		
		
		alf *= bt;
		R = R - alfAP;

		T  = A.ICCINV(R);
		bt = R* T;
		bt /= alf;

		P = T+ P*bt;
		eps = sqrt(R*R);
		k++;
	}while(eps > 1.0e-5);
	k++;
}
//计算阶乘
long Factorial(long n)
{
	if(0 == n) 
	{
		return 1;
	}
	if(1 == n)
	{
		return 1;
	}
	else
	{
		return n*Factorial(n-1);
	}
}
//三角单元面积分
long IntegralLij(long a, long b)
{
	return Factorial(a)*Factorial(b)/Factorial(a+b+1);
}
long IntegralLijm(long a, long b, long c)
{
	return 2*Factorial(a)*Factorial(b)*Factorial(c)/Factorial(a+b+c+2);
}
long IntegralLpijm(long a, long b, long c,long d)
{
	return 6*Factorial(a)*Factorial(b)*Factorial(c)*Factorial(d)/Factorial(a+b+c+d+3);
}

//计算行列式
double Determinant(long n, double** m_pdfData)
{
	if(2 == n)
	{
		return (m_pdfData[0][0]*m_pdfData[1][1] - m_pdfData[0][1]*m_pdfData[1][0]);
	}
	else
	{
		//create new data
		double** m_pdfDataNew;
		double sum = 0;
		m_pdfDataNew = new double*[n-1];
		for(long i = 0; i < n-1; i++)
		{
			m_pdfDataNew[i] = new double[n];
		}
		//		
		for(long l = 0; l < n; l++)
		{
			for(long i = 0,curi = 0; i < n-1; i++,curi++)
			{
				if(i == l) curi++;
				for(long j = 0; j < n-1; j++)
				{
					m_pdfDataNew[i][j] = m_pdfData[i+1][curi];
				}
			}
			sum += pow(-1,(l+2))*Determinant(n-1, m_pdfDataNew);
		}		
		//		
		for(long i = 0; i < n-1; i++)
		{
			delete[] m_pdfDataNew[i];;
		}
		delete[] m_pdfDataNew;
		return sum;
	}
}
#endif