mg_utils.cpp

pic 模拟程序！面向对象

#include "ovector.h"
#include "grid.h"

void boundary_Dirichlet(Scalar **, const intVector2&);
void boundary_Dirichlet_Zero(Scalar **, const intVector2&);
void boundary_Neumann_Zero(Scalar **, const intVector2&);
void GSRB_Boundary(Scalar**, Scalar**, const intVector2&, Scalar***, 
									 const int&, const int&, const int&);
void Residual_Boundary(Scalar**, Scalar**, Scalar**, const intVector2&, 
											 Scalar***, const int&, const int&);

#define NORTH 0
#define SOUTH 1
#define WEST 2
#define EAST 3
#define SOURCE 4

void GSRB(Scalar** p, Scalar** rhs, const intVector2& nhigh, Scalar ***Coeff, 
					const int& PeriodicFlagX1, const int& PeriodicFlagX2)
{
	int ioff, j, i, iinc;
	int nhigh1, nhigh2;
	nhigh1 = nhigh.e1();
	nhigh2 = nhigh.e2();
	Scalar* Cij;

	
	// Gauss - Seidel w/ Red - Black ordering, optimal parameter.
	for (ioff=0; ioff<2; ioff++){

		GSRB_Boundary(p, rhs, nhigh, Coeff, ioff, PeriodicFlagX1, PeriodicFlagX2);

		for (j=1; j<nhigh2; j++){
			iinc = (j+ioff)%2;
			for(i = 1 + iinc; i < nhigh1; i=i+2){
				Cij = Coeff[i][j];
				if (Cij[SOURCE]!=0.0)
					p[i][j] = p[i][j+1]*Cij[NORTH]+ Cij[SOUTH]*p[i][j-1] 
						+Cij[WEST]*p[i-1][j]+Cij[EAST]*p[i+1][j] + Cij[SOURCE]*rhs[i][j];
			}
		}
	}
}

void Interpolate(Scalar** fine,Scalar** coarse,const intVector2& nfhigh,
								 const intVector2& nchigh, Scalar ***FineCoeff)
{
	int nchigh1, nchigh2, i, j;
	
	nchigh1 = nchigh.e1();
	nchigh2 = nchigh.e2();
	if ((((Scalar)nfhigh.e1()/2)==nchigh.e1())&&(((Scalar)nfhigh.e2()/2)==nchigh.e2())){
		for(j = 1; j<=nchigh2; j++)
			for(i = 1; i <= nchigh1; i++){
				fine[2*i-1][2*j  ] +=  (FineCoeff[2*i-1][2*j  ][SOURCE]==0) ? 0 :  .5*(coarse[i-1][j]+coarse[i][j]); 
				fine[2*i  ][2*j-1] +=  (FineCoeff[2*i  ][2*j-1][SOURCE]==0) ? 0 :  .5*(coarse[i][j]+coarse[i][j-1]);
				fine[2*i-1][2*j-1] +=  (FineCoeff[2*i-1][2*j-1][SOURCE]==0) ? 0 : .25*(coarse[i][j]+coarse[i][j-1]+coarse[i-1][j-1]+coarse[i-1][j]); 
				fine[2*i  ][2*j  ] +=  (FineCoeff[2*i  ][2*j  ][SOURCE]==0) ? 0 :  coarse[i][j];
			}
		j=0;
		for(i = 1; i <= nchigh1; i++){
			fine[2*i-1][2*j  ] +=  (FineCoeff[2*i-1][2*j  ][SOURCE]==0) ? 0 : .5*(coarse[i-1][j]+coarse[i][j]); 
			fine[2*i  ][2*j  ] +=  (FineCoeff[2*i  ][2*j  ][SOURCE]==0) ? 0 :  coarse[i][j];
		}
		i=0;
		for(j = 1; j<=nchigh2; j++){
			fine[2*i  ][2*j-1] += (FineCoeff[2*i  ][2*j-1][SOURCE]==0) ? 0 : .5*(coarse[i][j]+coarse[i][j-1]);
			fine[2*i  ][2*j  ] += (FineCoeff[2*i  ][2*j  ][SOURCE]==0) ? 0 :  coarse[i][j];
		}
		j=0;
		fine[2*i  ][2*j  ] += (FineCoeff[2*i  ][2*j  ][SOURCE]==0) ? 0 : coarse[i][j];
	}
	else 	if (((nfhigh.e1()/2)==nchigh.e1())&&(nfhigh.e2()/2!=nchigh.e2()))	{
		for(j = 0; j<nchigh2; j++)
			for(i = 1; i < nchigh1; i++){
				fine[2*i-1][j] += (FineCoeff[2*i-1][j][SOURCE]==0) ? 0 : .5*(coarse[i][j]+coarse[i-1][j]);
				fine[2*i  ][j] += (FineCoeff[2*i  ][j][SOURCE]==0) ? 0 : coarse[i][j];
			}
		i=0;
		for(j = 0; j<nchigh2; j++)
			fine[2*i  ][j] += (FineCoeff[2*i  ][j][SOURCE]==0) ? 0 : coarse[i][j];
	}
	else if ((nfhigh.e1()/2!=nchigh.e1())&&(nfhigh.e2()/2==nchigh.e2())) {
		for(j = 1; j<nchigh2; j++)
			for(i = 0; i < nchigh1; i++)
				{
					fine[i][2*j-1] += (FineCoeff[i][2*j-1][SOURCE]==0) ? 0 :  .5*(coarse[i][j]+coarse[i][j-1]);
					fine[i][2*j  ] += (FineCoeff[i][2*j  ][SOURCE]==0) ? 0 :   coarse[i][j];
				}
		j=0;
		for(i = 0; i < nchigh1; i++)
			fine[i][2*j  ] += (FineCoeff[i  ][2*j  ][SOURCE]==0) ? 0 : coarse[i][j];
	}
}

void Average(Scalar** fine,Scalar** coarse,const intVector2& nfhigh, 
						 const intVector2& nchigh, const int& PeriodicFlagX1, const int& PeriodicFlagX2)
{
	int i,j;

	if ((((Scalar)nfhigh.e1()/2)==nchigh.e1())&&(((Scalar)nfhigh.e2()/2)==nchigh.e2())){
		for(j = 1; j < nchigh.e2(); j++)
			for(i = 1; i< nchigh.e1(); i++)
				coarse[i][j] = (fine[2*i-1][2*j-1] + fine[2*i-1][2*j+1] + 
												fine[2*i+1][2*j-1]+ fine[2*i+1][2*j+1] + 
												2*(fine[2*i-1][2*j] + 
													 fine[2*i][2*j-1] + fine[2*i][2*j+1]+ fine[2*i+1][2*j]) +
												4*fine[2*i][2*j])/16;
		if (PeriodicFlagX1){
			i=nchigh.e1();
			for(j = 1; j < nchigh.e2(); j++){
				coarse[i][j] = (fine[2*i-1][2*j-1] + fine[2*i-1][2*j+1] + 
												fine[1][2*j-1]+ fine[1][2*j+1] + 
												2*(fine[2*i-1][2*j] + 
													 fine[2*i][2*j-1] + fine[2*i][2*j+1]+ fine[1][2*j]) +
												4*fine[2*i][2*j])/16;
				coarse[0][j] = coarse[i][j];
			}
		}
		else {
			for(j = 1; j < nchigh.e2(); j++){
				i=0;
				coarse[i][j] = (fine[2*i+1][2*j-1]+ fine[2*i+1][2*j+1] + 
												2*(fine[2*i][2*j-1] + fine[2*i][2*j+1]+ fine[2*i+1][2*j]) +
												4*fine[2*i][2*j])/12;
				i=nchigh.e1();
				coarse[i][j] = (fine[2*i-1][2*j-1] + fine[2*i-1][2*j+1] + 
												2*(fine[2*i-1][2*j] + 
													 fine[2*i][2*j-1] + fine[2*i][2*j+1]) +
												4*fine[2*i][2*j])/12;
			}
		}
		if (PeriodicFlagX2){
			for(i = 1; i< nchigh.e1(); i++){
				j=nchigh.e2();
				coarse[i][j] = (fine[2*i-1][2*j-1] + fine[2*i-1][1] + 
												fine[2*i+1][2*j-1]+ fine[2*i+1][1] + 
												2*(fine[2*i-1][2*j] + 
													 fine[2*i][2*j-1] + fine[2*i][1]+ fine[2*i+1][2*j]) +
												4*fine[2*i][2*j])/16;
				coarse[i][0] = coarse[i][j];
			}
		}
		else{
			for(i = 1; i< nchigh.e1(); i++){
				j=0;
				coarse[i][j] = (fine[2*i-1][2*j+1] + 
												fine[2*i+1][2*j+1] + 
												2*(fine[2*i-1][2*j] + 
													 fine[2*i][2*j+1]+ fine[2*i+1][2*j]) +
												4*fine[2*i][2*j])/12;
				j=nchigh.e2();
				coarse[i][j] = (fine[2*i-1][2*j-1] + fine[2*i+1][2*j-1]+ 
												2*(fine[2*i-1][2*j] + 
													 fine[2*i][2*j-1] + fine[2*i+1][2*j]) +
												4*fine[2*i][2*j])/12;
			}
		}
		
		if (PeriodicFlagX1&&PeriodicFlagX2){
			i = nchigh.e1();
			j = nchigh.e2();
			coarse[i][j] = (fine[2*i-1][2*j-1] + fine[2*i-1][1] + 
												fine[1][2*j-1]+ fine[1][1] + 
												2*(fine[2*i-1][2*j] + 
													 fine[2*i][2*j-1] + fine[2*i][1]+ fine[1][2*j]) +
												4*fine[2*i][2*j])/16;
			coarse[0][j] = coarse[i][0] = coarse[0][0] = coarse[i][j];
		}  
		else if (PeriodicFlagX1){
			j = 0;
			i=nchigh.e1();
			coarse[i][j] = (fine[2*i-1][2*j+1] + 
											fine[1][2*j+1] + 
											2*(fine[2*i-1][2*j] + 
												 fine[2*i][2*j+1]+ fine[1][2*j]) +
											4*fine[2*i][2*j])/12;
			coarse[0][j] = coarse[i][j];
			j = nchigh.e2();
			coarse[i][j] = (fine[2*i-1][2*j-1] + 
											fine[1][2*j-1] + 
											2*(fine[2*i-1][2*j] + 
												 fine[2*i][2*j-1] + fine[1][2*j]) +
											4*fine[2*i][2*j])/12;
			coarse[0][j] = coarse[i][j];
		}
		else if (PeriodicFlagX2){
			i = nchigh.e1();
			j = nchigh.e2();
			coarse[i][j] = (fine[2*i-1][2*j-1] + fine[2*i-1][1] + 
											2*(fine[2*i-1][2*j] + 
												 fine[2*i][2*j-1] + fine[2*i][1]) +
											4*fine[2*i][2*j])/12;
			coarse[i][0] = coarse[i][j];
			i = 0;
			coarse[i][j] = (fine[2*i+1][2*j-1]+ fine[2*i+1][1] + 
												2*(fine[2*i][2*j-1] + fine[2*i][1]+ fine[2*i+1][2*j]) +
												4*fine[2*i][2*j])/12;
			coarse[i][0] = coarse[i][j];
		}
		else{
			i=0;j=0;
			coarse[i][j] = (fine[2*i+1][2*j+1] + 
											2*(fine[2*i][2*j+1]+ fine[2*i+1][2*j]) +
											4*fine[2*i][2*j])/9;
			i=nchigh.e1();
			coarse[i][j] = (fine[2*i-1][2*j+1] + 
											2*(fine[2*i-1][2*j] + fine[2*i][2*j+1])+
											4*fine[2*i][2*j])/9;	

			j=nchigh.e2();
			coarse[i][j] = (fine[2*i-1][2*j-1] + 
											2*(fine[2*i-1][2*j] +	fine[2*i][2*j-1]) + 
											4*fine[2*i][2*j])/9;

			i=0;
			coarse[i][j] = (fine[2*i+1][2*j-1]+
											2*(fine[2*i][2*j-1] + fine[2*i+1][2*j]) +
											4*fine[2*i][2*j])/9;
		}
	}
	
	else 	if (((nfhigh.e1()/2)==nchigh.e1())&&(nfhigh.e2()/2!=nchigh.e2())){
		for(j = 0; j <= nchigh.e2(); j++)
			for(i = 1; i< nchigh.e1(); i++)
				coarse[i][j] = .25*(fine[2*i+1][j]+2*fine[2*i][j] + fine[2*i-1][j]);

		if	(PeriodicFlagX1){
			for (j = 0; j <= nchigh.e2(); j++){
				i=nchigh.e1(); 
				coarse[i][j] = .25*(fine[1][j]+2*fine[2*i][j] + fine[2*i-1][j]);
				coarse[0][j] = coarse[i][j];
			}
		}
		else {
			for(j = 0; j <= nchigh.e2(); j++){
				i=0;
				coarse[i][j] = (fine[2*i+1][j]+2*fine[2*i][j])/3;
				i=nchigh.e1();
				coarse[i][j] = (2*fine[2*i][j] + fine[2*i-1][j])/3;
			}
		}
	}

	
	else if ((nfhigh.e1()/2!=nchigh.e1())&&(nfhigh.e2()/2==nchigh.e2())){
		for(j = 1; j < nchigh.e2(); j++)
			for(i = 0; i<= nchigh.e1(); i++)
				coarse[i][j] = .25*(fine[i][2*j+1] + 2*fine[i][2*j] + fine[i][2*j-1]);
		
	
		if (PeriodicFlagX2){
			for(i = 0; i<= nchigh.e1(); i++){
				j = nchigh.e2();
				coarse[i][j] = .25*(fine[i][1] + 2*fine[i][2*j] + fine[i][2*j-1]);
				coarse[i][0] = coarse[i][j];
			}
		}
		else{
			for(i = 0; i<= nchigh.e1(); i++){
				j=0;
				coarse[i][j] = (fine[i][2*j+1] + 2*fine[i][2*j])/3;
				j=nchigh.e2();
				coarse[i][j] = (2*fine[i][2*j] + fine[i][2*j-1])/3;
			}
		}
	}
	
}

void Sample(Scalar** fine,Scalar** coarse,const intVector2& nfhigh, 
						const intVector2& nchigh)
{
	int i,j;

	if ((((Scalar)nfhigh.e1()/2)==nchigh.e1())&&(((Scalar)nfhigh.e2()/2)==nchigh.e2())){
		for(j = 0; j <= nchigh.e2(); j++)
			for(i = 0; i <= nchigh.e1(); i++)
				coarse[i][j] = fine[2*i][2*j];
	}

	else 	if (((nfhigh.e1()/2)==nchigh.e1())&&(nfhigh.e2()/2!=nchigh.e2())){
		for(j = 0; j <= nchigh.e2(); j++)
			for(i = 1; i <= nchigh.e1(); i++)
				coarse[i][j] = fine[2*i][j];
	}
	
	else if ((nfhigh.e1()/2!=nchigh.e1())&&(nfhigh.e2()/2==nchigh.e2())){
		for(j = 0; j <= nchigh.e2(); j++)
			for(i = 0; i<= nchigh.e1(); i++)
				coarse[i][j] = fine[i][2*j];
	}
	
}
void Insert(Scalar** fine,Scalar** coarse,const intVector2& nfhigh,
						const intVector2& nchigh, Scalar ***FineCoeff)
{
	int nchigh1, nchigh2, i, j;
	
	nchigh1 = nchigh.e1();
	nchigh2 = nchigh.e2();
	if ((((Scalar)nfhigh.e1()/2)==nchigh.e1())&&(((Scalar)nfhigh.e2()/2)==nchigh.e2())){
		for(j = 0; j<=nchigh2; j++)
			for(i = 0; i <= nchigh1; i++)
				fine[2*i  ][2*j  ] +=  (FineCoeff[2*i  ][2*j  ][SOURCE]==0) ? 0 :  coarse[i][j];
	}
	
	else 	if (((nfhigh.e1()/2)==nchigh.e1())&&(nfhigh.e2()/2!=nchigh.e2()))	{
		for(j = 0; j<nchigh2; j++)
			for(i = 0; i < nchigh1; i++)
				fine[2*i  ][j] += (FineCoeff[2*i  ][j][SOURCE]==0) ? 0 : coarse[i][j];
	}

	else if ((nfhigh.e1()/2!=nchigh.e1())&&(nfhigh.e2()/2==nchigh.e2())) {
		for(j = 0; j<nchigh2; j++)
			for(i = 0; i < nchigh1; i++)
				fine[i][2*j  ] += (FineCoeff[i][2*j  ][SOURCE]==0) ? 0 : coarse[i][j];
	}
}


void CellAve(Scalar** fine,Scalar** coarse,const intVector2& nfhigh, 
						 const intVector2& nchigh)
{
	int i,j;

	if ((((Scalar)nfhigh.e1()/2)==nchigh.e1())&&(((Scalar)nfhigh.e2()/2)==nchigh.e2())){
		for(j = 0; j < nchigh.e2(); j++)
			for(i = 0; i < nchigh.e1(); i++)
				coarse[i][j] = (fine[2*i][2*j]+fine[2*i+1][2*j]+
												fine[2*i][2*j+1]+fine[2*i+1][2*j+1])/4;
	}

	else 	if (((nfhigh.e1()/2)==nchigh.e1())&&(nfhigh.e2()/2!=nchigh.e2())){
		for(j = 0; j < nchigh.e2(); j++)
			for(i = 0; i < nchigh.e1(); i++)
				coarse[i][j] = (fine[2*i][j]+fine[2*i+1][j])/2;
	}
	
	else if ((nfhigh.e1()/2!=nchigh.e1())&&(nfhigh.e2()/2==nchigh.e2())){
		for(j = 0; j < nchigh.e2(); j++)
			for(i = 0; i< nchigh.e1(); i++)
				coarse[i][j] = (fine[i][2*j]+fine[i][2*j+1])/2;
	}
	
}

Scalar Norm(Scalar** res, const intVector2& nhigh, Grid *grid)
{
	Scalar resnorm = 0.0;
	for (int j=0; j<= nhigh.e2(); j++)
		for (int i=0; i<= nhigh.e1(); i++)
			resnorm += sqr(res[i][j]);