package AMRElliptic;
import BoxTools.BoxLayout;
import BoxTools.BoxLayoutData;
import BoxTools.BoxedArray;
import BoxTools.Util;
import BoxTools.LayoutIterator;
import BoxTools.DataIndex;
import BoxTools.PAPIWrapper;
import AMRTools.QuadCFInterp;
import AMRTools.CoarseAverage;
import AMRTools.ConstantFineInterp;
import java.lang.Math;
/**
* <code>LevelOp</code> defines operations involving the Possion operator at one level. It also
* enforces the physical boundary condition. The physical boundary is cached for better performance.
* The definition of <code>LevelOp</code> depends on the current level layout and the layout at the
* coarser level. Among the operations it defines, included are two multigrid solvers for one level of
* data.	One for the case where there is a coarser level data, and another one for the bottom level.
* Intermediate data structures and operators are reused to save resources.
* <code>LevelOp</code> should be an interface or abstract class, which defines the interface between
* <code>AMRSolver</code> and specific definition of an elliptic solver. However, since Titanium 
* currently does not support inheritance involving a template class, it is impossible for a class 
* for example <code>PossionOp</code> to extends or implements <code>LevelOp</code>. The current
* implementation just combines these two classes into one.
* The details of the multigrid algorithms are available in Chombo Design Document.
* <p> Usage: a <code>LevelOp </code> object must be declared local. 
* @see BoxLayout,BoxLayoutData,BoxedArray,LayoutIterator,QuadCFInterp,CoarseAverage,ConstantFineInterp.
* @see <a href="http://davis.lbl.gov/NASA/chombomanual/classLevelOp.html">Chombo Specification</a>
* @version 1.1
*		Modified on Jul 09, 2004. change of syntax of foreach, see BoxTools/foreachLoopTest.ti.
*			 on Aug 26, 2004. BoxLayoutData.<<=constant covers ghost cells.
* @author Tong Wen, LBNL
* @since 1.0
*/ 
public abstract class LevelOp{
	protected static final int  SPACE_DIM  =Util.SPACE_DIM;
	protected static final int  SORT_DIM  =Util.SORT_DIM;
	protected final static int highBoundary=Util.HIGH_SIDE;
	protected final static int lowBoundary=Util.LOW_SIDE;
	protected final static RectDomain<SPACE_DIM> NULL_BOX=Util.NULL_BOX;
	protected final static int numVCycle=Util.NUMVCYCLE;
	protected final static int numBVCycle=Util.NUMBVCYCLE;

	final static double R=0.25;
	final static double a=2/(25*R*R*R);
	final static double b=3/(16*R*R);

	/*the above info should be obtained only from Util class.*/
	protected  static Point<SPACE_DIM> DISP1, DISP2, DISP3;
	
	BoxLayout local single m_layout,m_layoutC;
  	double single m_Dx;
	int single m_refRatio; //must be a power of 2
	int single m_mgLevels;
	boolean single m_homogeneousOnly;
	RectDomain<SPACE_DIM> single m_domain;
	RectDomain<SPACE_DIM> single m_domainc;
	int single m_numSmoothDown, m_numSmoothUp;
	boolean single m_isBottom;
	//add public for test
	public LayoutIterator local m_layoutItr, m_layoutItrC;
	BoxLayout local [] local m_boundaries=new BoxLayout local [2*SPACE_DIM];
	boolean [] local [] local m_isPhyBoundary=new boolean [2*SPACE_DIM] local [] local;
	template QuadCFInterp<double > local single m_CFinterp;
	template CoarseAverage<double > local single m_Average;
        template ConstantFineInterp<double > local single m_Ipwc;	

	/** data structures for mgRelax */
	template BoxLayoutData<template BoxedArray<double> > local single [] local mg_phi, mg_rhs, mg_R;
	template QuadCFInterp<double > local single [] local mg_CFinterp;
	template CoarseAverage<double > local single [] local mg_Average;
	template ConstantFineInterp<double > local single [] local  mg_Ipwc;
	double single [] local mg_dx;
	
	/** date structures for bottom smoother */	
	int single bs_levels;
	LayoutIterator local bs_layoutItr;
	RectDomain<SPACE_DIM> [] local bs_boundaries=new RectDomain<SPACE_DIM> [2*SPACE_DIM];
	boolean [] local bs_isPhyBoundary=new boolean [2*SPACE_DIM];
	template BoxLayoutData<template BoxedArray<double> > local single [] local bs_phi, bs_rhs, bs_R;
	template QuadCFInterp<double > local single [] local bs_CFinterp;
	template CoarseAverage<double > local single [] local bs_Average;
	template ConstantFineInterp<double > local single [] local  bs_Ipwc;
	double single [] local bs_dx;

	/** Timer **/
	public static int countAverage; Timer tmrAverage =new Timer();
	public static int countIpwc; Timer tmrIpwc=new Timer();
	public static int countSmooth; Timer tmrSmooth=new Timer();
	public static int countResidual; Timer tmrResidual=new Timer();
	public static int countGSRB; Timer tmrGSRB=new Timer();
	public static int countUpdate; Timer tmrUpdate=new Timer();
	public static int countSExchange; Timer tmrSExchange=new Timer();
	public static int countSHPhyBC; Timer tmrSHPhyBC=new Timer();
	public static int countSCFInterp; Timer tmrSCFInterp=new Timer();
	public static Timer local tmrBSsmooth=new Timer(), tmrBS=new Timer();
	public static double flopsSmooth; 
	
	public static double GSRBL1DCM, GSRBL2DCM;
	public static double GSRBL1DCA, GSRBL2DCA;
	PAPIWrapper papiFPINS=null, papiFMAINS=null;// papiL1DCM=null, papiL2DCM=null, papiL1DCA=null, papiL2DCA=null;

	public LevelOp(){
		m_Dx=1;m_refRatio=1; 	
	}

	public local final void define(BoxLayout local single a_boxLayout, BoxLayout local single a_boxLayoutC,
	     double single a_DxLevel, 
	     int single a_refRatio, boolean single a_homogeneousOnly, RectDomain<SPACE_DIM> single a_domain,
	     RectDomain<SPACE_DIM> single a_domainc, int single a_numSmoothDown, int single a_numSmoothUp){
		//System.out.println("This is levelOp constructor");
		m_layout=a_boxLayout; 
		m_layoutC=a_boxLayoutC;
		m_Dx=a_DxLevel;
		m_refRatio=a_refRatio;
		//m_mgLevels=(int single)(Math.log(a_refRatio)/Math.log(2)+1.5);
		m_homogeneousOnly=a_homogeneousOnly;
		m_domain=a_domain;
		m_domainc=a_domainc;
		m_numSmoothDown=a_numSmoothDown;
		m_numSmoothUp=a_numSmoothUp;

		m_layoutItr=new LayoutIterator(m_layout);
		m_layoutItrC=new LayoutIterator(m_layoutC);
		int single i, dir, side;
		int j,k;
		k=m_layout.numBoxesAt(Ti.thisProc());
		for (i=0;i<2*SPACE_DIM;i++){
			dir=(i+2)/2;
			side=(i+2)%2;
			if (side==highBoundary)
				m_boundaries[i]=m_layout.highBoundary(dir,1,false);
			else
				m_boundaries[i]=m_layout.lowBoundary(dir,1,false);
			m_isPhyBoundary[i]=new boolean [k];
			m_layoutItr.reset();
			j=0;
			while (!m_layoutItr.isEnded()){
				if (m_boundaries[i][m_layoutItr.getDataIndex()]<m_domain)
					m_isPhyBoundary[i][j]=false;
				else
					m_isPhyBoundary[i][j]=true;
				j++;
				m_layoutItr.advance();
			}
		}
		//System.out.println("Start CFInterp0");
		m_CFinterp=new template QuadCFInterp<double >(m_layout,m_Dx,m_refRatio,1,m_domain, m_layoutC, m_domainc);
		//System.out.println("Finish CFInterp0");
		m_Average=new template CoarseAverage<double >(m_layout, m_refRatio);
		m_Ipwc=new template ConstantFineInterp<double >(m_layoutC, m_refRatio);
		
		/*Here, it is assumed that the base domain is square or cubic, and that it is decomposed evenly.*/
		int single numOfBoxes=(int single)m_layout.size();
		final int single adjustment=2;
		if (m_refRatio==1 || m_layout==m_layoutC){
			m_isBottom=true;
			if (numOfBoxes==1){
				m_mgLevels=1;
				bs_levels=(int single)(Math.log(Math.pow(m_domain.size(),1.0/3))/Math.log(2)+1.5);
			}
			else{	
				/*(int ) rounds off a double toward zero!.
				m_mgLevels=(int single)(Math.log(Math.pow(m_domain.size(),1.0/3))/Math.log(2)+1.5);
				bs_levels=(int single)(Math.log(Math.pow(numOfBoxes,1.0/3))/Math.log(2)+1.5)+adjustment;
				m_mgLevels-=(bs_levels-1);
				bs_levels+=1;*/
				m_mgLevels=3;
				bs_levels=7;
			}
		}
		else{
			m_isBottom=false;
			m_mgLevels=(int single)(Math.log(a_refRatio)/Math.log(2)+1.5);
			bs_levels=0;
		}
		Util.printMsg("Levelop: m_mgLevels "+m_mgLevels+"bs_levels "+bs_levels+"? "+m_isBottom);
		mg_phi= new template BoxLayoutData<template BoxedArray<double> > local single [m_mgLevels];
		mg_rhs= new template BoxLayoutData<template BoxedArray<double> > local single [m_mgLevels];
		mg_R= new template BoxLayoutData<template BoxedArray<double> > local single [m_mgLevels];
		//System.out.println("To construct AMRTools:CFInterp");
		mg_CFinterp= new template QuadCFInterp<double > local single [m_mgLevels];
		//System.out.println("To construct AMRTools:Average");
		mg_Average= new template CoarseAverage<double > local single [m_mgLevels];
		//System.out.println("To construct AMRTools:Ipwc");
		mg_Ipwc= new template ConstantFineInterp<double > local single [m_mgLevels];
		mg_dx=new double single [m_mgLevels];
		
		i=m_mgLevels-1;
		double single dx=m_Dx;
		int single refRatio=m_refRatio;
		RectDomain<SPACE_DIM> single domain=m_domain;
		RectDomain<SPACE_DIM> single domainc;
		BoxLayout local single bl=m_layout;
		BoxLayout local single blc;
		//System.out.println("To construct"+m_mgLevels+"level AMRTools.");
		while (i>0){ //i>=0
		  if (i!=m_mgLevels-1){
			mg_phi[i]=new template BoxLayoutData<template BoxedArray<double> >(bl,1);
			mg_rhs[i]=new template BoxLayoutData<template BoxedArray<double> >(bl,0);
		  }
		  mg_dx[i]=dx;
		  //domainc=(RectDomain<SPACE_DIM> single)Util.coarsen(domain,2);
		  //mg_CFinterp[i]=new template QuadCFInterp<double >(bl,dx,2,1,domain,bl.coarsen(2), domainc);
		  //System.out.println("This is mg level"+i);
		  if (i!=0){
			mg_R[i]=new template BoxLayoutData<template BoxedArray<double> >(bl,0);
			mg_Average[i]=new template CoarseAverage<double >(bl, 2);
			//System.out.println("1: "+domain.toString());
			domainc=(RectDomain<SPACE_DIM> single)Util.coarsen(domain,2);
			blc=bl.coarsen(2);
		  	mg_CFinterp[i]=new template QuadCFInterp<double >(bl,dx,refRatio,1,domain,blc,domainc);
			bl=blc;
			mg_Ipwc[i]=new template ConstantFineInterp<double >(bl,2);
			dx*=2;
			refRatio=(int single)refRatio/2;
			domain=domainc;
		  }
		  else
		  	mg_CFinterp[i]=new template QuadCFInterp<double >(bl,dx,1,1,domain,bl,domain);
		  i--;
		}
		
		if (DISP1==Point<SPACE_DIM>.all(0))
		  switch (SPACE_DIM){
			case 3: 
				DISP3=Point<SPACE_DIM>.direction(3);
			case 2:
				DISP2=Point<SPACE_DIM>.direction(2);
			case 1: 
				DISP1=Point<SPACE_DIM>.direction(1);	
		  }
		//for (i=0;i<m_mgLevels;i++) System.out.print(" "+mg_dx[i]); System.out.print("\n"); System.out.flush();
		/** Constructing data structures for bottom smoother. */
		if (m_isBottom){
			Domain<SPACE_DIM> domain1;
			RectDomain<SPACE_DIM> single rectdomain;
					 
			 if (m_mgLevels==1){  
				rectdomain=m_domain; dx=m_Dx;
			 }
			 else{ 
				rectdomain=(RectDomain<SPACE_DIM> single)Util.refine(domain,2);dx=dx/2;
			 }

			   for (i=0;i<2*SPACE_DIM;i++){
				dir=(i+2)/2;
				side=(i+2)%2;
				if (side==highBoundary)
					bs_boundaries[i]=rectdomain.border(1,dir,1);
				else
					bs_boundaries[i]=rectdomain.border(1,-dir,1);
				/*this if statement is redundant
				if (bs_boundaries[i]<m_domain)
					bs_isPhyBoundary[i]=false;
				else
					bs_isPhyBoundary[i]=true;
				*/
			   }	
			   RectDomain<SPACE_DIM> single [] single domains={rectdomain};	
		 	   int single [] single procIDs={0};
			   bl=new BoxLayout(domains,procIDs);
			   //bl.myprint();Ti.barrier();
			   bs_layoutItr=new LayoutIterator(bl);
			   
			   i=bs_levels-1;
			   //dx=dx/2;
			   Util.printMsg("define BS:"+rectdomain.toString()+"bs_levels"+bs_levels+" dx "+dx);
			   domain=rectdomain;
			   bs_phi= new template BoxLayoutData<template BoxedArray<double> > local single [bs_levels];
			   bs_rhs= new template BoxLayoutData<template BoxedArray<double> > local single [bs_levels];
			   bs_R= new template BoxLayoutData<template BoxedArray<double> > local single [bs_levels];
			   bs_CFinterp= new template QuadCFInterp<double > local single [bs_levels];
			   bs_Average= new template CoarseAverage<double > local single [bs_levels];
			   bs_Ipwc= new template ConstantFineInterp<double > local single [bs_levels];
			   bs_dx=new double single [bs_levels];
			   
			   while (i>=0){
			     bs_phi[i]=new template BoxLayoutData<template BoxedArray<double> >(bl,1);
			     bs_rhs[i]=new template BoxLayoutData<template BoxedArray<double> >(bl,0);
			     bs_dx[i]=dx;
			     if (i!=0){
				bs_R[i]=new template BoxLayoutData<template BoxedArray<double> >(bl,0);
				bs_Average[i]=new template CoarseAverage<double >(bl, 2);
				domainc=(RectDomain<SPACE_DIM> single)Util.coarsen(domain,2);
				blc=bl.coarsen(2);
			     	//bs_CFinterp[i]=new template QuadCFInterp<double >(bl,dx,2,1,domain,blc,domainc);
				bl=blc;
				bs_Ipwc[i]=new template ConstantFineInterp<double >(bl,2);
				dx*=2;
				domain=domainc;
		  	     }
			     //else
				//bs_CFinterp[i]=new template QuadCFInterp<double >(bl,dx,1,1,domain,bl,domain);
		 	     i--;
			   }
		//System.out.println("Finishing bottome construction");
		}
	}

	public LevelOp(BoxLayout local single a_boxLayout,BoxLayout local single a_boxLayoutC,
	     double single a_DxLevel,  int single a_refRatio, 
	     boolean single a_homogeneousOnly, RectDomain<SPACE_DIM> single a_domain,
	     RectDomain<SPACE_DIM> single a_domainc, int single a_numSmoothDown, int single a_numSmoothUp){
		define(a_boxLayout, a_boxLayoutC, a_DxLevel, a_refRatio, a_homogeneousOnly, 
			a_domain, a_domainc, a_numSmoothDown, a_numSmoothUp);
		   papiFPINS=new PAPIWrapper(PAPICounter.PAPI_FP_INS);
		   papiFMAINS=new PAPIWrapper(PAPICounter.PAPI_FMA_INS);

		/*
		   papiL1DCM=new PAPIWrapper(PAPICounter.L1_DataCache_Misses);
		   papiL2DCM=new PAPIWrapper(PAPICounter.L2_DataCache_Misses);
		   papiL1DCA=new PAPIWrapper(PAPICounter.L1_DataCache_Accesses);
		   papiL2DCA=new PAPIWrapper(PAPICounter.L2_DataCache_Accesses);
		*/
	}

	/** Interfaces to m_CFinterp, m_Average, m_Ipwc. Operations are done without exchange.*/
	public local final void CFInterp(template BoxLayoutData<template BoxedArray<double > > local single a_phif, template BoxLayoutData<template BoxedArray<double > > local single  a_phic){
		m_CFinterp.coarseFineInterp(a_phif,a_phic);
	}
	public final local void CFInterp(template BoxLayoutData<template BoxedArray<double > > local single a_phif){
		m_CFinterp.coarseFineInterp(a_phif);
	}
	public final local void Average(template BoxLayoutData<template BoxedArray<double > > local single a_phic, template BoxLayoutData<template BoxedArray<double > > local single  a_phif){
		m_Average.averageToCoarse(a_phic, a_phif);
	}
	public final local inline void Ipwc(template BoxLayoutData<template BoxedArray<double > > local single a_phif, template BoxLayoutData<template BoxedArray<double > > local single  a_phic){
		m_Ipwc.interpToFine(a_phif, a_phic);
	}
	public final inline local double getDx(){
		return m_Dx;
	}
	public final local void applyPhyBC(template BoxLayoutData<template BoxedArray<double > > local single a_phi){
		if (a_phi.boxLayout()==m_layout)
		  if (m_homogeneousOnly)	
			applyHPhyBC(a_phi,m_mgLevels-1);
		  else
			applyImHPhyBC(a_phi);
		else
			Util.printErrMsg("PoissonOp::applyPhyBC: invalid layout!");
	}
	public final local void applyHPhyBC(template BoxLayoutData<template BoxedArray<double > > local single a_phi){
		if (a_phi.boxLayout()==m_layout)	
			applyHPhyBC(a_phi,m_mgLevels-1);
		else
			Util.printErrMsg("PoissonOp::applyPhyBC: invalid layout!");
	}
	//public local inline double gfunction(Point<SPACE_DIM> a_point, Point<SPACE_DIM> a_disp){return 0;}
	abstract public local double gfunction(Point<SPACE_DIM> a_point, Point<SPACE_DIM> a_disp);
	public final local void applyImHPhyBC(template BoxLayoutData<template BoxedArray<double > > local single a_phi){
		if (a_phi.boxLayout()==m_layout){
			PrivateRegion local region=new PrivateRegion();
			int i,j,dir,side;
			Point<SPACE_DIM> [1d] local disp =new(region) Point<SPACE_DIM> [[1:SPACE_DIM]];
			for (i=1;i<=SPACE_DIM;i++) disp[i]=Point<SPACE_DIM>.direction(i);
			RectDomain<SPACE_DIM> box;
			//Point<SPACE_DIM> point; /* for change in the syntax of foreach */
			DataIndex local DI;
			double [SPACE_DIM d] local TiArray;
			
			for (i=0;i<2*SPACE_DIM;i++){
				dir=(i+2)/2;
				side=(i+2)%2;
				m_layoutItr.reset();
				j=0;
				while (!m_layoutItr.isEnded()){
					if (m_isPhyBoundary[i][j]){
						DI=m_layoutItr.getDataIndex();
						box=m_boundaries[i][DI];	
						TiArray=a_phi[DI].getLocalArray();
						/*foreach (point in box ) {
						   if (side==highBoundary)
							TiArray[point]=2*gfunction(point, disp[dir]*(-1))-TiArray[point-disp[dir]];
						   else
							TiArray[point]=2*gfunction(point, disp[dir])-TiArray[point+disp[dir]];
						}*/
						/*** Note that in Titanium "if" statement is very expensive.***/
						if (side==highBoundary)
							foreach (point in box ) 
								TiArray[point]=2*gfunction(point, disp[dir]*(-1))-TiArray[point-disp[dir]];
						   else
							foreach (point in box ) 
								TiArray[point]=2*gfunction(point, disp[dir])-TiArray[point+disp[dir]];
					}				
					j++;
					m_layoutItr.advance();
				}
			}
			try {
		    		region.delete();
			} 
			catch (RegionInUse e){
				Util.printErrMsg("LevelOp::applyImHPhyBC(): failed to delete PrivateRegion pr");
			}
		}			
		else
			Util.printErrMsg("PoissonOp::applyImHPhyBC: invalid layout!");
	}
	final local void applyHPhyBC(template BoxLayoutData<template BoxedArray<double > > local single a_phi, int single a_mgLevel){
		int i,j,dir,side;
		PrivateRegion local region=new PrivateRegion();
		Point<SPACE_DIM> [1d] local disp =new(region) Point<SPACE_DIM> [[1:SPACE_DIM]];
		for (i=1;i<=SPACE_DIM;i++) disp[i]=Point<SPACE_DIM>.direction(i);
		RectDomain<SPACE_DIM> box;
		//Point<SPACE_DIM> point; /* for change in the syntax of foreach */
		DataIndex local DI;
		double [SPACE_DIM d] local TiArray;
		int factor=(int)Math.pow(2,m_mgLevels-1-a_mgLevel);

		for (i=0;i<2*SPACE_DIM;i++){
			dir=(i+2)/2;
			side=(i+2)%2;
			m_layoutItr.reset();
			j=0;
			while (!m_layoutItr.isEnded()){
				if (m_isPhyBoundary[i][j]){
					DI=m_layoutItr.getDataIndex();
					box=m_boundaries[i][DI];	
					TiArray=a_phi[DI].getLocalArray();
					if (factor !=1) box=Util.coarsen(box,factor);
					/*foreach (point in box ) {
					   if (side==highBoundary)
						TiArray[point]=-TiArray[point-disp[dir]];
					   else
						TiArray[point]=-TiArray[point+disp[dir]];
					}*/
					/*** Note that in Titanium "if" statement is very expensive.***/
					if (side==highBoundary)
						foreach (point in box )
							TiArray[point]=-TiArray[point-disp[dir]];
					else
						foreach (point in box )
							TiArray[point]=-TiArray[point+disp[dir]];
				}				
				j++;
				m_layoutItr.advance();
			}
		}
		try {
		    region.delete();
		} 
		catch (RegionInUse e){
			Util.printErrMsg("LevelOp::applyHPhyBC(): failed to delete PrivateRegion pr");
		}
	}
	final local void BSapplyHPhyBC(template BoxLayoutData<template BoxedArray<double > > local single a_phi, int single a_bsLevel){
		int i,dir,side;
		PrivateRegion local region=new PrivateRegion();
		Point<SPACE_DIM> [1d] local disp =new(region) Point<SPACE_DIM> [[1:SPACE_DIM]];
		for (i=1;i<=SPACE_DIM;i++) disp[i]=Point<SPACE_DIM>.direction(i);
		RectDomain<SPACE_DIM> box;
		//Point<SPACE_DIM> point; /* for change in the syntax of foreach */
		double [SPACE_DIM d] local TiArray;
		DataIndex local DI;
		int factor=(int)Math.pow(2,bs_levels-1-a_bsLevel);
		if (Ti.thisProc()==0){
		  bs_layoutItr.reset();
		  DI=bs_layoutItr.getDataIndex();
		  TiArray=a_phi[DI].getLocalArray();
		  for (i=0;i<2*SPACE_DIM;i++){
			dir=(i+2)/2;
			side=(i+2)%2;
			/*** The following "if" statment may be redundant because everyone is physical boundary ***/
			//if (bs_isPhyBoundary[i]){
				box=bs_boundaries[i];	
				if (factor !=1) box=Util.coarsen(box,factor);
				if (side==highBoundary)
					foreach (point in box ) {TiArray[point]=-TiArray[point-disp[dir]];}
				 else
					foreach (point in box ) {TiArray[point]=-TiArray[point+disp[dir]];}
		  	//}				
		  }
		}
		try {
		    region.delete();
		} 
		catch (RegionInUse e){
			Util.printErrMsg("LevelOp::BSapplyHPhyBC(): failed to delete PrivateRegion pr");
		}
	}

	//local void applyOp(template BoxLayoutData<template BoxedArray<double > > local single 
	//	a_phi, template BoxLayoutData<template BoxedArray<double > > local single a_LOfPhi){}
	abstract local void applyOp(template BoxLayoutData<template BoxedArray<double > > local single 
		a_phi, template BoxLayoutData<template BoxedArray<double > > local single a_LOfPhi);

	public final local void applyOpH(template BoxLayoutData<template BoxedArray<double > > local single 
		a_phi, template BoxLayoutData<template BoxedArray<double > > local single a_LOfPhi,
		boolean single a_CFBoundaryCheckOn, boolean single a_PhysBoundaryCheckOn){
		/***
		 There are three kinds of ghost cells:
		 Type 1 should be filled by exchange;
		 Type 2 should be filled by CFInterp;
		 Type 3 should be filled by Physical Boundary Condition.
		 For a_phi, it is assumed that the exchange operation has been done, and that
		 all the ghost cells have zero values.
		***/
		if (a_phi.boxLayout()==m_layout && a_LOfPhi.boxLayout()==m_layout)	
			applyOp(a_phi,a_LOfPhi);
		else
			Util.printErrMsg("PoissonOp::applyOpH: invalid layout!");
	}

	public final local void applyOpIcfHphys(template BoxLayoutData<template BoxedArray<double > > local single 
		a_phi, template BoxLayoutData<template BoxedArray<double > > local single a_LOfPhi,
		boolean single a_CFInterpOn, boolean single a_PhysBoundaryCheckOn){
		/***
		 For a_phi, it is assumed that the exchange operation has been done.
		***/
		if (a_phi.boxLayout()==m_layout && a_LOfPhi.boxLayout()==m_layout){
			if (a_CFInterpOn) CFInterp(a_phi);
			if (a_PhysBoundaryCheckOn) applyHPhyBC(a_phi);	
			applyOp(a_phi,a_LOfPhi);
		}
		else
			Util.printErrMsg("PoissonOp::applyOpIcfHphys: invalid layout!");
	}

	public final local void applyOpIcfIphys(template BoxLayoutData<template BoxedArray<double > > local single 
		a_phi, template BoxLayoutData<template BoxedArray<double > > local single a_phic,
		template BoxLayoutData<template BoxedArray<double > > local single a_LOfPhi,
		boolean single a_CFInterpOn, boolean single a_PhysBoundaryCheckOn){
		/* 
		 need to decide how to enforce the Homogeneous physical boundary condition. 
		 For a_phi and a_phic, it is assumed that the exchange operation has been done.
		*/
		if (a_phi.boxLayout()==m_layout && a_LOfPhi.boxLayout()==m_layout){
			if (a_CFInterpOn) CFInterp(a_phi,a_phic);
			if (a_PhysBoundaryCheckOn) applyImHPhyBC(a_phi);	
			applyOp(a_phi,a_LOfPhi);
		}
		else
			Util.printErrMsg("PoissonOp::applyOpIcfHphys: invalid layout!");
	}

	abstract public local void GetFlux(template BoxLayoutData<template BoxedArray<double > > local single 
		a_flux, template BoxLayoutData<template BoxedArray<double > > local single a_phi, 
		int single a_dir, int single a_side);

	abstract public local void getFlux(template BoxedArray<double > local a_flux, template BoxedArray<double > local a_phi, int a_dir, int  a_side);
	//add public for test purpose Dec 09,03
	abstract public local void residual(template BoxLayoutData<template BoxedArray<double > > local single 
		a_phi, template BoxLayoutData<template BoxedArray<double > > local single a_rhs,
		double single a_dx, 
		template BoxLayoutData<template BoxedArray<double > > local single a_residual);

	final local void smooth(template BoxLayoutData<template BoxedArray<double > > local single a_phi, 
				template BoxLayoutData<template BoxedArray<double > > local single a_rhs,
		 		double single a_dx, template QuadCFInterp<double > local single a_CFInterp,
				int single a_mgLevel){
		tmrSCFInterp.start();a_CFInterp.coarseFineInterp(a_phi);tmrSCFInterp.stop();
		smooth(a_phi, a_rhs, a_dx,a_CFInterp);
		tmrSExchange.start();a_phi.exchange();tmrSExchange.stop();		
		tmrSHPhyBC.start();applyHPhyBC(a_phi,a_mgLevel);tmrSHPhyBC.stop();
		countSExchange+=tmrSExchange.millis();tmrSExchange.reset();
		countSHPhyBC+=tmrSHPhyBC.millis();tmrSHPhyBC.reset();
		countSCFInterp+=tmrSCFInterp.millis();tmrSCFInterp.reset();
	}

	final local void smoothNoUpdate(template BoxLayoutData<template BoxedArray<double > > local single a_phi, 
				template BoxLayoutData<template BoxedArray<double > > local single a_rhs,
		 		double single a_dx, template QuadCFInterp<double > local single a_CFInterp,
				int single a_mgLevel){
		tmrSCFInterp.start();a_CFInterp.coarseFineInterp(a_phi);tmrSCFInterp.stop();
		smooth(a_phi, a_rhs, a_dx,a_CFInterp);
		countSCFInterp+=tmrSCFInterp.millis();tmrSCFInterp.reset();
	}
	final local void BSsmooth(template BoxLayoutData<template BoxedArray<double > > local single a_phi, 
				template BoxLayoutData<template BoxedArray<double > > local single a_rhs,
		 		double single a_dx, template QuadCFInterp<double > local single a_CFInterp,
				int single a_bsLevel){
		//a_CFInterp.coarseFineInterp(a_phi); //must be redundant.
		//smooth(a_phi, a_rhs, a_dx,a_CFInterp);
		BSsmooth(a_phi, a_rhs, a_dx,a_CFInterp);
		BSapplyHPhyBC(a_phi,a_bsLevel);
	}
	/*** 
	 Note that the relaxation scheme fails at the 1 by 1 case. 
	 So we need an explicity solve. But the it seems that
	 it doesn't matter to use it or not.	
	***/
	final local void BS1by1(template BoxLayoutData<template BoxedArray<double > > local single a_phi, 
				template BoxLayoutData<template BoxedArray<double > > local single a_rhs,
		 		double single a_dx){
		
		double [SPACE_DIM d] local TiArray1, TiArray2;
		
		//Point<SPACE_DIM> point; /* for change in the syntax of foreach */
		DataIndex local DI;
		bs_layoutItr.reset();
		while(!bs_layoutItr.isEnded()){
			DI=bs_layoutItr.getDataIndex();
			TiArray1=a_phi[DI].getLocalArray();TiArray2=a_rhs[DI].getLocalArray();
			foreach (point in TiArray2.domain()) TiArray1[point]=-(a_dx*a_dx)/12*TiArray2[point];
			bs_layoutItr.advance();
		}	
	}

	//add public for test purpose Dec 12,03
	abstract local public single void smooth(template BoxLayoutData<template BoxedArray<double > > local single a_phi, 
				template BoxLayoutData<template BoxedArray<double > > local single a_rhs,
		 		double single a_dx, template QuadCFInterp<double > local single a_CFInterp);

	abstract local  void BSsmooth(template BoxLayoutData<template BoxedArray<double > > local single a_phi, 
				template BoxLayoutData<template BoxedArray<double > > local single a_rhs,
		 		double single a_dx, template QuadCFInterp<double > local single a_CFInterp);

	public final local void levelGSRB(template BoxLayoutData<template BoxedArray<double > > local single 
		a_phi, template BoxLayoutData<template BoxedArray<double > > local single a_rhs){
		if ( a_phi.boxLayout()!=m_layout || a_rhs.boxLayout()!=m_layout)
			Util.printErrMsg("LevelOp::LevelGSRB(): invalid layout info!");
		else{
			smooth(a_phi,a_rhs,m_Dx, m_CFinterp, m_mgLevels-1);
		}
	}	

	final local void mgRelax(int single a_level){
		//System.out.println("   This is mgRelax start at level:"+a_level);
		int single l=a_level;
		int single i,j;
		if (m_isBottom && l==1){
			tmrBS.start();
			bs_phi[bs_levels-1].copy(mg_phi[1],false);//true is changed to false Oct 22,03
			BSapplyHPhyBC(bs_phi[bs_levels-1],bs_levels-1);
			bs_rhs[bs_levels-1].copy(mg_rhs[1], false);
			//Util.printMsg("m_domf: "+mg_CFinterp[1].m_domf.toString()); 
			tmrBSsmooth.start();bottomSmoother((int single)(bs_levels-1));tmrBSsmooth.stop();
			mg_phi[1].copy(bs_phi[bs_levels-1], false);
			//mg_phi[1].exchange();applyHPhyBC(mg_phi[1],1);//redundant
			tmrBS.stop();	
		}
		else{
			tmrSmooth.start();
			for (i=0;i<m_numSmoothDown;i++) smooth(mg_phi[l],mg_rhs[l], mg_dx[l], mg_CFinterp[l], l);
			tmrSmooth.stop();
		}
		
		if (l>1){
			mg_phi[l-1]<<=0;
			
			//tmrUpdate.start();
				//mg_phi[l-1].exchange();
				//applyHPhyBC(mg_phi[l-1],l-1);
			//tmrUpdate.stop();
			//mg_CFinterp[l-1].coarseFineInterp(mg_phi[l-1]);
			mg_CFinterp[l].coarseFineInterp(mg_phi[l]);
			tmrResidual.start();
				residual(mg_phi[l],mg_rhs[l],mg_dx[l],mg_R[l]);
			tmrResidual.stop();
			tmrAverage.start();
				mg_Average[l].averageToCoarse(mg_rhs[l-1], mg_R[l]);
			tmrAverage.stop();
			mgRelax(l-1);
			tmrIpwc.start();mg_Ipwc[l].interpToFine(mg_R[l],mg_phi[l-1]);tmrIpwc.stop();
			mg_phi[l]+=mg_R[l];
			tmrUpdate.start();
				mg_phi[l].exchange();
				applyHPhyBC(mg_phi[l],l);
			tmrUpdate.stop();
			tmrSmooth.start();
				for (i=0;i<m_numSmoothUp-1;i++) smooth(mg_phi[l],mg_rhs[l], mg_dx[l], mg_CFinterp[l],l);
				if (l==(m_mgLevels-1))
					smooth(mg_phi[l],mg_rhs[l], mg_dx[l], mg_CFinterp[l],l);	
				else	
					smoothNoUpdate(mg_phi[l],mg_rhs[l], mg_dx[l], mg_CFinterp[l],l);
			tmrSmooth.stop();
		}

		//System.out.println("   This is mgRelax end at level:"+a_level);
		countAverage+=tmrAverage.millis();
		countIpwc+=tmrIpwc.millis();
		countSmooth+=tmrSmooth.millis();
		countResidual+=tmrResidual.millis();
		countUpdate+=tmrUpdate.millis();
		tmrAverage.reset();
		tmrIpwc.reset();
		tmrSmooth.reset();
		tmrResidual.reset();
		tmrUpdate.reset();
	}

	public final local  void mgRelax(template BoxLayoutData<template BoxedArray<double > > local single a_phi, 
				  template BoxLayoutData<template BoxedArray<double > > local single a_rhs){
		if ( a_phi.boxLayout()!=m_layout || a_rhs.boxLayout()!=m_layout)
			Util.printErrMsg("LevelOp::mgRelax(.): invalid layout info!");
		else{
			mg_phi[m_mgLevels-1]=a_phi;
			mg_rhs[m_mgLevels-1]=a_rhs;
			//for (int single i=0;i<numVCycle;i++) 
				mgRelax((int single)(m_mgLevels-1));
		}

	}
	
	final local void bottomSmoother(int single a_level){
		int single l=a_level;
		int single i,j;
		if (l>0)
			for (i=0;i<m_numSmoothDown;i++) BSsmooth(bs_phi[l],bs_rhs[l], bs_dx[l], bs_CFinterp[l],l);
		else
			BS1by1(bs_phi[l],bs_rhs[l], bs_dx[l]);
		//System.out.println("This is first smooth at bs_level:"+l+" Down "+m_numSmoothDown);
		if (l>0){
			bs_phi[l-1]<<=0;
			//bs_phi[l-1].exchange();//must be redundant.
			//BSapplyHPhyBC(bs_phi[l-1],l-1);
			//bs_CFinterp[l-1].coarseFineInterp(bs_phi[l-1]);//must be redundant.
			residual(bs_phi[l],bs_rhs[l],bs_dx[l],bs_R[l]);
			bs_Average[l].averageToCoarse(bs_rhs[l-1], bs_R[l]);
			bottomSmoother(l-1);
			bs_Ipwc[l].interpToFine(bs_R[l],bs_phi[l-1]);
			bs_phi[l]+=bs_R[l];
			//bs_phi[l].exchange(); //must be redundant.
			BSapplyHPhyBC(bs_phi[l],l);
			for (i=0;i<m_numSmoothUp;i++) BSsmooth(bs_phi[l],bs_rhs[l], bs_dx[l], bs_CFinterp[l],l);
			//System.out.println("This is second smooth at bs_level:"+l+" Up "+m_numSmoothUp);	
		}
	}
	//for test purpose
	public final local double maxLocalNorm(template BoxLayoutData<template BoxedArray<double > > local single a_phi){
		LayoutIterator local LItr=new LayoutIterator(a_phi.boxLayout());
		double [SPACE_DIM d] local TiArray;
		RectDomain<SPACE_DIM> box;
		//Point<SPACE_DIM> point, pointMin, pointMax; /* for change in the syntax of foreach */
		Point<SPACE_DIM> pointMin, pointMax;
		DataIndex local DI;
		double maxNorm=0;
		double temp;
		double Max=-10000000000;
		double Min= 10000000000;
		while(!LItr.isEnded()){
			box=LItr.access();
			DI=LItr.getDataIndex();
			TiArray=a_phi[DI].getLocalArray();
			foreach (point in box){
				temp=TiArray[point];
				if (temp>Max){Max=temp;pointMax=point;}
				if (temp<Min){Min=temp;pointMin=point;}
			}
			LItr.advance();	
		}	
		System.out.println("maxLocalNorm: the range is { "+Min+" at "+pointMin.toString()+"; "+Max+" at "+pointMax.toString()+" }");
		if (maxNorm<Math.abs(Max)) maxNorm=Math.abs(Max);if (maxNorm<Math.abs(Min)) maxNorm=Math.abs(Min);		
		return maxNorm;	
	}
	//for test purpose
	public final local double LocalNorm1(template BoxLayoutData<template BoxedArray<double > > local single a_phi){
		LayoutIterator local LItr=new LayoutIterator(a_phi.boxLayout());
		double [SPACE_DIM d] local TiArray;
		RectDomain<SPACE_DIM> box;
		//Point<SPACE_DIM> point; /* for change in the syntax of foreach */
		DataIndex local DI;
		double Norm=0;
		while(!LItr.isEnded()){
			box=LItr.access();
			DI=LItr.getDataIndex();
			TiArray=a_phi[DI].getLocalArray();
			foreach (point in box) Norm+=Math.abs(TiArray[point]);
			LItr.advance();	
		}	
		return Norm;	
	}
	//for test purpose
	public final local void BUpDownTest(template BoxLayoutData<template BoxedArray<double > > local single a_phi){
		int single i,j,l;
		l=bs_levels-1;
		bs_phi[l].copy(a_phi,false);
		Util.printErrMsg("Up&Down: at level "+l+" 1norm "+LocalNorm1(bs_phi[l]));
		for (i=l;i>0;i--){
			bs_phi[i-1]<<=0;
			bs_Average[i].averageToCoarse(bs_phi[i-1], bs_phi[i]);
			Util.printErrMsg("Up&Down: at level "+i+" 1norm "+LocalNorm1(bs_phi[i-1]));	
		}
		for (i=0;i<l;i++){
			bs_phi[i+1]<<=0;
			bs_Ipwc[i+1].interpToFine(bs_phi[i+1],bs_phi[i]);
			Util.printErrMsg("Up&Down: at level "+i+" 1norm "+LocalNorm1(bs_phi[i+1]));	
		}
	}
	
	public final local  void bottomSmoother(template BoxLayoutData<template BoxedArray<double > > local single a_phi, 
				  template BoxLayoutData<template BoxedArray<double > > local single a_rhs){
		if ( !m_isBottom){
			Util.printErrMsg("LevelOp::bottomSmoother(.): invalid calling!");
			return;
		}
			
		if (m_mgLevels==1){
			tmrBS.start();
			//bs_phi[bs_levels-1]<<=0;
			bs_phi[bs_levels-1].copy(a_phi,false);//true is changed to false Oct 22,03
			BSapplyHPhyBC(bs_phi[bs_levels-1],bs_levels-1);
			bs_rhs[bs_levels-1].copy(a_rhs, false);
			tmrBSsmooth.start();
			for (int single i=0;i<numBVCycle;i++) 
				bottomSmoother((int single)(bs_levels-1));
			tmrBSsmooth.stop();
			a_phi.copy(bs_phi[bs_levels-1], false);
			//a_phi.exchange();applyHPhyBC(a_phi);//e[l-1] is used, but maybe redundant;
			tmrBS.stop();
		}
		else{
			 for (int single i=0;i<numBVCycle;i++) 
				mgRelax(a_phi,a_rhs);

		}

	}
}