femarix.cpp

有限单元

#include "StdAfx.h"
#include "FeMarix.h"
#include ".\femarix.h"
CFeMarix::CFeMarix()
{
}
CFeMarix::CFeMarix(long m_lNIn,CArray<long,long>& m_lMarixBandIn)
{
	m_lN = m_lNIn;
	Init(m_lN, m_lMarixBand);
}
CFeMarix::CFeMarix(CFeMarix& a)
{
	m_lN = a.m_lN;
	m_lMarixBand.RemoveAll();
	for(long i = 0; i < a.m_lMarixBand.GetCount(); i++)
	{
		m_lMarixBand.Add(a.m_lMarixBand[i]);
	}
	data = new double*[m_lMarixBand.GetCount()];
	ic_data = new double*[m_lMarixBand.GetCount()];
	for(long i = 0; i < m_lMarixBand.GetCount(); i++)
	{
		data[i] = new double[m_lN - m_lMarixBand[i]];
		ic_data[i] = new double[m_lN - m_lMarixBand[i]];
	}
	//zero
	for(long i = 0; i < m_lMarixBand.GetCount(); i++)
	{
		for(long j = 0; j < m_lN - m_lMarixBand[i]; j++)
		{
			data[i][j] = a.data[i][j];
			ic_data[i][j] = 0;
		}
	}
}
CFeMarix::CFeMarix(long m_lNIn, CArray<long,long>& m_lMarixBandIn,double** dataIn)
{
	m_lN = m_lNIn;
	for(long i = 0; i < m_lMarixBandIn.GetCount(); i++)
	{
		m_lMarixBand.Add(m_lMarixBandIn[i]);
	}
	data = new double*[m_lMarixBand.GetCount()];
	ic_data = new double*[m_lMarixBand.GetCount()];
	for(long i = 0; i < m_lMarixBand.GetCount(); i++)
	{
		data[i] = new double[m_lN - m_lMarixBand[i]];
		ic_data[i] = new double[m_lN - m_lMarixBand[i]];
	}
	
	//zero
	for(long i = 0; i < m_lMarixBand.GetCount(); i++)
	{
		for(long j = 0; j < m_lN - m_lMarixBand[i]; j++)
		{
			data[i][j] = dataIn[i][j];
			ic_data[i][j] = 0;
		}
	}
}
void CFeMarix::Init(long m_lNIn,CArray<long,long>& m_lMarixBandIn)
{
	m_lN = m_lNIn;
	for(long i = 0; i < m_lMarixBandIn.GetCount(); i++)
	{
		m_lMarixBand.Add(m_lMarixBandIn[i]);
	}
	data = new double*[m_lMarixBand.GetCount()];
	ic_data = new double*[m_lMarixBand.GetCount()];
	for(long i = 0; i < m_lMarixBand.GetCount(); i++)
	{
		data[i] = new double[m_lN - m_lMarixBand[i]];
		ic_data[i] = new double[m_lN - m_lMarixBand[i]];
	}
	
	//zero
	for(long i = 0; i < m_lMarixBand.GetCount(); i++)
	{
		for(long j = 0; j < m_lN - m_lMarixBand[i]; j++)
		{
			data[i][j] = 0;
			ic_data[i][j] = 0;
		}
	}
}
CFeMarix::~CFeMarix(void)
{	
	if(m_lMarixBand.GetCount())
	{
		for(long i = 0; i < m_lMarixBand.GetCount(); i++)
		{
			delete[] data[i];
			delete[] ic_data[i];
		}
		delete[] data;
		delete[] ic_data;
		m_lMarixBand.RemoveAll();
	}
}
CFeMarix& CFeMarix::operator =(CFeMarix& a )
{
	m_lN = a.m_lN;
	m_lMarixBand.RemoveAll();
	for(long i = 0; i < a.m_lMarixBand.GetCount(); i++)
	{
		m_lMarixBand.Add(a.m_lMarixBand[i]);
	}
	data = new double*[m_lMarixBand.GetCount()];
	ic_data = new double*[m_lMarixBand.GetCount()];
	for(long i = 0; i < m_lMarixBand.GetCount(); i++)
	{
		data[i] = new double[m_lN - m_lMarixBand[i]];
		ic_data[i] = new double[m_lN - m_lMarixBand[i]];
	}
	//zero
	for(long i = 0; i < m_lMarixBand.GetCount(); i++)
	{
		for(long j = 0; j < m_lN - m_lMarixBand[i]; j++)
		{
			data[i][j] = a.data[i][j];
			ic_data[i][j] = a.ic_data[i][j];
		}
	}	
	return* this;
}

CFeMarix& CFeMarix::operator-= (CFeMarix& a)
{
	for(long i = 0; i < m_lMarixBand.GetCount(); i++)
	{
		for(long j = 0; j < m_lN - m_lMarixBand[i]; j++)
		{
			data[i][j] -= a.data[i][j];
		}
	}
	return *this;
}
CFeMarix& CFeMarix::operator+= (CFeMarix& a)
{
	for(long i = 0; i < m_lMarixBand.GetCount(); i++)
	{
		for(long j = 0; j < m_lN - m_lMarixBand[i]; j++)
		{
			data[i][j] += a.data[i][j];
		}
	}
	return *this;
}
CFeMarix& CFeMarix::operator*= (double a)
{
	for(long i = 0; i < m_lMarixBand.GetCount(); i++)
	{
		for(long j = 0; j < m_lN - m_lMarixBand[i]; j++)
		{
			data[i][j] *= a;
		}
	}
	return *this;
}
CFeVector CFeMarix::operator* (CFeVector& a)
{
	double* v  = new double[m_lN];
	for(long i = 0; i < m_lN; i++)
	{
		v[i] = 0.0;
		for(long j = Judge(i); j > 0; j--)
		{
			v[i] += data[j][i - m_lMarixBand[j]]*a.v[i - m_lMarixBand[j]];
		}
		for(long j = 0; j < m_lMarixBand.GetCount() && m_lMarixBand[j]+i <m_lN ; j++)
		{
			v[i] += data[j][i]*a.v[i+m_lMarixBand[j]];
		}
	}
	return CFeVector(m_lN,v);
}
long CFeMarix::Judge(long i)
{
	for(long j = m_lMarixBand.GetCount() - 1; j >= 0; j --)
	{
		if( i >= m_lMarixBand[j])
			return j;
	}
	return -1;
}
void CFeMarix::Set(CFeMarix& a)
{
	for(long i = 0; i < m_lMarixBand.GetCount(); i++)
	{
		for(long j = 0; j < m_lN - m_lMarixBand[i]; j++)
		{
			data[i][j] = a.data[i][j];
		}
	}
}

void CFeMarix::Set_E(long i, long j, double a)
{
	if(j>=i)
	{
		for(long k = 0; k < m_lMarixBand.GetCount(); k++)
		{
			if(i + m_lMarixBand[k] == j)
			{
				data[k][i] = a;
			}
		}
	}	
}
double CFeMarix::Get_E(long i, long j)
{
	if(j>=i)
	{
		for(long k = 0; k < m_lMarixBand.GetCount(); k++)
		{
			if(i + m_lMarixBand[k] == j)
			{
				return data[k][i];
			}
			else
			{
				return 0;
			}
		}
	}
	else
	{
		return CFeMarix::Get_E(j, i);
	}
}
void CFeMarix::Add_E(long i, long j, double a)
{
	if(j>=i)
	{
		for(long k = 0; k < m_lMarixBand.GetCount(); k++)
		{
			if(i + m_lMarixBand[k] == j)
			{
				data[k][i] += a;
			}
		}
	}	
}

void CFeMarix::ICC(void)
{
	long i = 0;
	long j = 0;
	long k = 0;
	long l = 0;
	long m = 0;
	long q = 0;
	long n = 0;
	double sum = 0;
	double ca = 0;
	//data -> ic_data
	for( j = 0; j < m_lMarixBand.GetCount() && i + m_lMarixBand[j] < m_lN; j++)
	{
		::CopyMemory(ic_data[j],data[j],(m_lN - m_lMarixBand[j])*sizeof(double));
	}

	//
	for(i = 0; i < m_lN; i++)
	{
		for(j = Judge(i); j > 0; j--)
		{
			if( fabs(ic_data[0][i-m_lMarixBand[j]]) < 1.0e-20)
			{
				return;
			}
			ca = ic_data[j][i-m_lMarixBand[j]]/ic_data[0][i-m_lMarixBand[j]];
			for( k = i-m_lMarixBand[j]+1; k < i; k++)
			{
				for( l = Judge(k); l >0; l--)
				{
					if(k - m_lMarixBand[l] == i - m_lMarixBand[j])
					{
						ic_data[j][i-m_lMarixBand[j]] = ic_data[j][i-m_lMarixBand[j]] - ca*ic_data[l][k - m_lMarixBand[l]];
						break;
					}
				}
			}
		}
		ca = fabs(ic_data[0][i]);
		if( ca < 1.0e-20 )
		{
			return;
		}
		for(l = Judge(i); l > 0; l--)
		{
			 ca +=  fabs(ic_data[l][i- m_lMarixBand[l]]);
		}
		ic_data[0][i] = ca;
	}
	
}
CFeVector CFeMarix::ICCINV(CFeVector& R)
{

	CFeVector reV(R.m_lN);
	
	//解方程
	long i = 0;
	long j = 0;
	double sum = 0;
	for( i = 0; i < m_lN; i++)
	{
		for( sum = R.v[i], j = Judge(i); j > 0; j--)
		{
			sum -= ic_data[j][i-m_lMarixBand[j]]*reV.v[i-m_lMarixBand[j]]/ic_data[0][i-m_lMarixBand[j]];
		}
		reV.v[i] = sum;//ic_data[0][i];
	}
	//
	for( i = m_lN-1; i >= 0 ; i--)
	{
		for(j = 1; j < m_lMarixBand.GetCount() && i+m_lMarixBand[j] < m_lN; j++)
		{
			reV.v[i] -= reV.v[i+m_lMarixBand[j]]*ic_data[j][i];
		}
		reV.v[i] /= ic_data[0][i];
	}
	return reV;
}

void CFeMarix::Save(void)
{
	FILE* fp = fopen("E:\\data.dat","w+");
	for(long i = 0; i < m_lN; i++)
	{
		for(long j = 0; j < m_lMarixBand.GetCount() && i+m_lMarixBand[j] < m_lN ; j++)
		{
			fprintf(fp,"%f\t",data[j][i]);
		}
		fprintf(fp,"\n");
	}
	fclose(fp);
}

void CFeMarix::CreateData(long m_lNin, CArray<long,long>& m_lMarixBandIn, double** datades)
{
	datades = new double*[m_lMarixBandIn.GetCount()];
	for(long i = 0; i < m_lMarixBandIn.GetCount(); i++)
	{
		datades[i] = new double[m_lN - m_lMarixBandIn[i]];
	}
}
void CFeMarix::CholeskyEquation(CFeVector& X, CFeVector& B)
{
	//乔莱斯基分解解方程组
	double* m_dfB = B.v;
	double* m_dfX = X.v;
	long i,j,k;
	double sum;	

	double** m_dfA;
	m_dfA = new double*[m_lN];
	for( i = 0; i < m_lN; i++)
	{
		m_dfA[i] = new double[m_lN];
	}

	for( i = 0; i < m_lN; i++)
	{
		for( j = 0; j < m_lN; j++)
		{
			m_dfA[i][j] = 0.0;
		}
	}
	for(i = 0; i < m_lN;i++)
	{
		for(j = 0; j < m_lMarixBand.GetCount() && i+m_lMarixBand[j] < m_lN; j++)
		{
			m_dfA[i][i+m_lMarixBand[j]] = data[j][i];
		}
	}
	for(i = 0; i < m_lN;i++)
	{
		for(j = i; j < m_lN; j++)
		{
			for(sum = m_dfA[i][j], k = i-1; k >= 0; k--)
			{
				sum -= m_dfA[i][k]*m_dfA[j][k];
			}
			if(i == j)
			{
				if( sum <= 0.0)
				{
					//error
					sum = 1;					
				}
				m_dfA[i][i] = sqrt(sum);
            }
			else
			{
				m_dfA[j][i] = sum/m_dfA[i][i];
			}
		}
	}
	//
	for( i = 0; i < m_lN; i++)
	{
		for(sum = m_dfB[i],k = i - 1;k >= 0;k--)
		{
			sum -= m_dfA[i][k]*m_dfX[k];			
		}
		m_dfX[i] = sum/m_dfA[i][i];
	}
	for( i = m_lN-1; i>= 0; i--)
	{
		for(sum = m_dfX[i],k = i+1; k < m_lN; k++)
		{
			sum -= m_dfA[k][i]*m_dfX[k];
		}
		m_dfX[i] = sum/m_dfA[i][i];
	}
	//
	
	for( i = 0; i < m_lN; i++)
	{
		delete[] m_dfA[i];
	}
	delete[] m_dfA;

}

void CFeMarix::Zero(void)
{
	for(long i = 0; i < m_lMarixBand.GetCount(); i++)
	{
		for(long j = 0; j < m_lN - m_lMarixBand[i]; j++)
		{
			data[i][j] = 0;
		}
	}
}