mg_utils.cpp

pic 模拟程序！面向对象

//			resnorm += sqr(res[i][j]*grid->get_halfCellVolumes()[i][j]);

	resnorm = sqrt(resnorm);
	
	return resnorm;
}

void Residual(Scalar** p,Scalar** rhs,Scalar** res, 
							const intVector2 &nhigh, Scalar*** Coeff, const int& PeriodicFlagX1,
							const int& PeriodicFlagX2)
{
	int j,k;
	Scalar* Cjk;

	Residual_Boundary(p, rhs, res, nhigh, Coeff, PeriodicFlagX1, PeriodicFlagX2);	

	for(j = 1; j<nhigh.e1(); j++)
		for(k = 1; k<nhigh.e2(); k++){
			Cjk = Coeff[j][k];
			if (Cjk[SOURCE]!=0)
				res[j][k] = -p[j][k]*Cjk[SOURCE] - p[j-1][k]*Cjk[WEST] 
					- p[j+1][k]*Cjk[EAST] 
						- p[j][k-1]*Cjk[SOUTH] - p[j][k+1]*Cjk[NORTH] + rhs[j][k];

		}
}

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

	if (PeriodicFlagX2){
		j=0;
		for(i = 1; i < nhigh1; i++){
			Cij = Coeff[i][j];
			if (Cij[SOURCE]!=0){
				res[i][j] = -p[i][j+1]*Cij[NORTH] - Cij[WEST]*p[i-1][j] - 
					Cij[EAST]*p[i+1][j] - Cij[SOUTH]*p[i][nhigh2-1] 
						- Cij[SOURCE]*p[i][j] + rhs[i][j];
				res[i][nhigh2] = res[i][j];
			}
		}
	}
	else{
		j=0;
		for(i = 1; i < nhigh1; i++){
			Cij = Coeff[i][j];
			if (Cij[SOURCE]!=0)
				res[i][j] = -p[i][j+1]*Cij[NORTH] - Cij[WEST]*p[i-1][j] - Cij[EAST]*p[i+1][j] 
					- Cij[SOURCE]*p[i][j] + rhs[i][j];
		}
		
		j=nhigh2;
		for(i = 1; i < nhigh1; i++){
			Cij = Coeff[i][j];
			if (Cij[SOURCE]!=0)
				res[i][j] = -Cij[SOUTH]*p[i][j-1]-Cij[WEST]*p[i-1][j]-Cij[EAST]*p[i+1][j] 
					- Cij[SOURCE]*p[i][j] + rhs[i][j];
		}
	}
	
	if (PeriodicFlagX1){
		i = 0;
		for (j=1; j<nhigh2; j++){
			Cij = Coeff[i][j];
			if (Cij[SOURCE]!=0){
				res[i][j] = -p[i][j+1]*Cij[NORTH]- Cij[SOUTH]*p[i][j-1] 
					- Cij[EAST]*p[i+1][j] - Cij[WEST]*p[nhigh1-1][j]  
						- Cij[SOURCE]*p[i][j] + rhs[i][j];
				res[nhigh1][j] = res[i][j];
			}
		}
	}
	else{
		i = 0;
		for (j=1; j<nhigh2; j++){
			Cij = Coeff[i][j];
			if (Cij[SOURCE]!=0)
				res[i][j] = -p[i][j+1]*Cij[NORTH]- Cij[SOUTH]*p[i][j-1] 
					-Cij[EAST]*p[i+1][j] - Cij[SOURCE]*p[i][j] + rhs[i][j];
		}
		
		i = nhigh1;
		for (j=1; j<nhigh2; j++){
			Cij = Coeff[i][j];
			if (Cij[SOURCE]!=0)
				res[i][j] = -p[i][j+1]*Cij[NORTH]- Cij[SOUTH]*p[i][j-1] 
					-Cij[WEST]*p[i-1][j] - Cij[SOURCE]*p[i][j] +rhs[i][j];
		}
	}

	if (!PeriodicFlagX1&&!PeriodicFlagX2){
		i=0;
		j=0;
		Cij = Coeff[i][j];
		if (Cij[SOURCE]!=0)
			res[i][j] = -p[i][j+1]*Cij[NORTH]-Cij[EAST]*p[i+1][j] 
				- Cij[SOURCE]*p[i][j]+rhs[i][j];
		
		i=nhigh1;
		Cij = Coeff[i][j];
		if (Cij[SOURCE]!=0)
			res[i][j] = -p[i][j+1]*Cij[NORTH]-Cij[WEST]*p[i-1][j] 
				- Cij[SOURCE]*p[i][j] + rhs[i][j];
		
		j=nhigh2;
		Cij = Coeff[i][j];
		if (Cij[SOURCE]!=0)
			res[i][j] = -Cij[SOUTH]*p[i][j-1] -Cij[WEST]*p[i-1][j]
				- Cij[SOURCE]*p[i][j] + rhs[i][j];
		
		i=0;
		Cij = Coeff[i][j];
		if (Cij[SOURCE]!=0)
			res[i][j] = -Cij[SOUTH]*p[i][j-1] -Cij[EAST]*p[i+1][j] 
				- Cij[SOURCE]*p[i][j] + rhs[i][j];
	}
	else if (PeriodicFlagX1&&PeriodicFlagX2){
		i=0;
		j=0;
		Cij = Coeff[i][j];
		if (Cij[SOURCE]!=0){
			res[i][j] = -p[i][j+1]*Cij[NORTH]-Cij[EAST]*p[i+1][j] 
				-Cij[WEST]*p[nhigh1-1][j] - Cij[SOUTH]*p[i][nhigh2-1] - 
					Cij[SOURCE]*p[i][j]+rhs[i][j];
			res[nhigh1][j] = res[i][nhigh2] = res[nhigh1][nhigh2] = res[j][j];
		}				
	}
	else if (PeriodicFlagX1){
		i=0;
		j=0;
		Cij = Coeff[i][j];
		if (Cij[SOURCE]!=0){
			res[i][j] = -p[i][j+1]*Cij[NORTH]-Cij[EAST]*p[i+1][j]
				- Cij[WEST]*p[nhigh1-1][j]
					- Cij[SOURCE]*p[i][j]+rhs[i][j];
			res[nhigh1][j] = res[i][j];
		}
		
		j=nhigh2;
		Cij = Coeff[i][j];
		if (Cij[SOURCE]!=0){
			res[i][j] = -Cij[SOUTH]*p[i][j-1] -Cij[WEST]*p[nhigh1-1][j]
				- Cij[EAST]*p[i+1][j] 
					- Cij[SOURCE]*p[i][j] + rhs[i][j];
			res[nhigh1][j] = res[i][j];
		}
	}
	else{
		i=0;
		j=0;
		Cij = Coeff[i][j];
		if (Cij[SOURCE]!=0){
			res[i][j] = -p[i][j+1]*Cij[NORTH]-Cij[EAST]*p[i+1][j]
				- Cij[SOUTH]*p[i][nhigh2-1]
					- Cij[SOURCE]*p[i][j]+rhs[i][j];
			res[i][nhigh2] = res[i][j];
		}
		i=nhigh1;
		Cij = Coeff[i][j];
		if (Cij[SOURCE]!=0){
			res[i][j] = -p[i][j+1]*Cij[NORTH]-Cij[WEST]*p[i-1][j] 
				- Cij[SOUTH]*p[i][nhigh2-1]
					- Cij[SOURCE]*p[i][j] + rhs[i][j];
			res[i][nhigh2] =res[i][j];
		}
	}
}
void GSRB_Boundary(Scalar** p, Scalar** rhs, const intVector2& nhigh, 
									 Scalar*** Coeff, const int& ioff, const int& PeriodicFlagX1,
									 const int& PeriodicFlagX2)
{
	int i,j,iinc;
	Scalar* Cij;

	int nhigh1, nhigh2;
	
	nhigh1 = nhigh.e1();
	nhigh2 = nhigh.e2();

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

	if (!PeriodicFlagX1&&!PeriodicFlagX2){
		if (ioff){
			i=0;
			j=0;
			Cij = Coeff[i][j];
			if (Cij[SOURCE]!=0)
				p[i][j] = p[i][j+1]*Cij[NORTH]+Cij[EAST]*p[i+1][j] 
					+ Cij[SOURCE]*rhs[i][j];
		}
		
		if ((nhigh1%2)!=ioff){
			i=nhigh1;
			j=0;
			Cij = Coeff[i][j];
			if (Cij[SOURCE]!=0)
				p[i][j] = p[i][j+1]*Cij[NORTH]+Cij[WEST]*p[i-1][j] 
					+ Cij[SOURCE]*rhs[i][j];
		}
		
		if (((nhigh1+nhigh2)%2)!=ioff){
			i=nhigh1;
			j=nhigh2;
			Cij = Coeff[i][j];
			if (Cij[SOURCE]!=0)
				p[i][j] = Cij[SOUTH]*p[i][j-1] +Cij[WEST]*p[i-1][j]
					+ Cij[SOURCE]*rhs[i][j];
		}
		
		if ((nhigh2%2)!=ioff){
			i=0;
			j=nhigh2;
			Cij = Coeff[i][j];
			if (Cij[SOURCE]!=0)
				p[i][j] = Cij[SOUTH]*p[i][j-1] +Cij[EAST]*p[i+1][j] 
					+ Cij[SOURCE]*rhs[i][j];
		}
	}
	else if (PeriodicFlagX1&&PeriodicFlagX2){
		if (ioff){
			i=0;
			j=0;
			Cij = Coeff[i][j];
			if (Cij[SOURCE]!=0){
				p[i][j] = p[i][j+1]*Cij[NORTH]+Cij[EAST]*p[i+1][j]
					+ Cij[SOUTH]*p[i][nhigh2-1] + Cij[WEST]*p[nhigh1-1][j]
						+ Cij[SOURCE]*rhs[i][j];
				p[i][nhigh2]=p[nhigh1][j]=p[nhigh1][nhigh2]=p[i][j];
			}
		}
	}
	else if (PeriodicFlagX1){
		if (ioff){
			i=0;
			j=0;
			Cij = Coeff[i][j];
			if (Cij[SOURCE]!=0){
				p[i][j] = p[i][j+1]*Cij[NORTH]+Cij[EAST]*p[i+1][j] 
					+ Cij[WEST]*p[nhigh1-1][j]
						+ Cij[SOURCE]*rhs[i][j];
				p[nhigh1][j] = p[i][j];
			}
		}
		
		
		if ((nhigh2%2)!=ioff){
			j=nhigh2;
			Cij = Coeff[i][j];
			if (Cij[SOURCE]!=0){
				p[i][j] = Cij[SOUTH]*p[i][j-1] +Cij[EAST]*p[i+1][j] 
					+ Cij[WEST]*p[nhigh1-1][j]
						+ Cij[SOURCE]*rhs[i][j];
				p[nhigh1][j] = p[i][j];
			}	
		}
	}
	else{
		if (ioff){
			i=0;
			j=0;
			Cij = Coeff[i][j];
			if (Cij[SOURCE]!=0){
				p[i][j] = p[i][j+1]*Cij[NORTH]+Cij[EAST]*p[i+1][j] 
					+ Cij[SOUTH]*p[i][nhigh2]
						+ Cij[SOURCE]*rhs[i][j];
				p[i][nhigh2-1] = p[i][j];
			}
			
			if ((nhigh1%2)!=ioff){
				i=nhigh1;
				Cij = Coeff[i][j];
				if (Cij[SOURCE]!=0){
					p[i][j] = p[i][j+1]*Cij[NORTH]+Cij[WEST]*p[i-1][j]
						+ Cij[SOUTH]*p[i][nhigh2]
							+ Cij[SOURCE]*rhs[i][j];
					p[i][nhigh2-1] = p[i][j];
				}
			}
		}	
	}	
}