/**Remember to add stride, the default value now is 1. May 06, 2002*/
package AMRTools;
import BoxTools.BoxLayout;
import BoxTools.BoxLayoutData;
import BoxTools.LayoutIterator;
import BoxTools.BoxedArray;
import BoxTools.DataIndex;
import BoxTools.Copier;
import BoxTools.Util;
/**
* A <code>QuadCFInterp</code> object is used to interpolate the values on the fine ghost cells
* at the coarse-fine interface. The accuracy is O(h^3) so that for example we can obtain at least 
* an O(h) approximation to the Laplacian. Both the data from the fine level and the coarse level
* are used. Hence, this class is defined on two levels of <code>BoxLayout</code>. The coarse-fine
* boundary is calculated once then cached. This algorithm works on <code>BoxLayoutData<\code>
* whose components are <code>BoxedArray<\code>s.
* Please refer to the Chombo design document for the details of this algorithm. 
* <p> Usage: a <code>QuadCFInterp</code> object must be declared local. 
* @see BoxTools.BoxLayout, BoxTools.BoxLayoutData, BoxTools.LayoutIterator, BoxTools.BoxedArray
*      and BoxTools.DataIndex.
* @see <a href="http://davis.lbl.gov/NASA/chombomanual/classQuadCFInterp.html">Chombo Specification</a>
* @version 1.1 
*		Modified on Jul 09, 2004. change of syntax of foreach, see BoxTools/foreachLoopTest.ti.
* @author Tong Wen, LBNL
* @since 1.0
*/ 
template<class T> public class QuadCFInterp{
	private static final int SPACE_DIM=BoxTools.Util.SPACE_DIM;
	private static final RectDomain<SPACE_DIM> NULL_BOX=BoxTools.Util.NULL_BOX;
	private final static int highBoundary=1;
	private final static int lowBoundary=0;

	private BoxLayout local single m_layoutf, m_layoutc;
	private LayoutIterator local m_Iteratorf;
	//private  double   single m_dx;
	private GridSpacing local m_DX;
	private int    single m_refRatio;
	//private int    single m_nComp;
	public RectDomain<SPACE_DIM>   single m_domf;
	private RectDomain<SPACE_DIM>   single m_domc;
	
	//RectDomain<SPACE_DIM> [] local [] local m_boundaryf,m_boundaryc;
	//Domain<SPACE_DIM> local [] local [] local m_boundaryf,m_boundaryc;//changed on 11/25/03
	CFBoundary  [] local [] local m_boundaryf,m_boundaryc;//changed on 09/15/04

	//template BoxLayoutData<template BoxedArray<T > > single local [] local m_phiBoundaryc;
	template BoxLayoutData<template BoxedArray<T > > single local m_phiBoundaryc;
	//template Copier<template BoxedArray<T > > local single m_copier;

	Timer tmr1=new Timer();public static int Tint1;
	Timer tmr2=new Timer();public static int Tint2;
	Timer tmrCopy=new Timer();public static int TintCopy;
	public static int IsDomainCount;

	public QuadCFInterp(){
		m_refRatio=1;//m_nComp=1;
	}
	public QuadCFInterp(BoxLayout local single a_layoutf, double   single a_dx, int single a_refRatio, 
				int single a_nComp, RectDomain<SPACE_DIM> single  a_domf, 
				BoxLayout local single a_layoutc, RectDomain<SPACE_DIM> single  a_domc){
		GridSpacing local a_DX=new GridSpacing(a_dx);
		define(a_layoutf, a_DX, a_refRatio, a_nComp, a_domf, a_layoutc, a_domc);
	}

	public QuadCFInterp(BoxLayout local single a_layoutf, GridSpacing local a_DX, int single a_refRatio, 
				int single a_nComp, RectDomain<SPACE_DIM> single  a_domf, 
				BoxLayout local single a_layoutc, RectDomain<SPACE_DIM> single  a_domc){
		define(a_layoutf, a_DX, a_refRatio, a_nComp, a_domf, a_layoutc, a_domc);
	}

	public final local void define(BoxLayout local single a_layoutf, GridSpacing local  a_DX, 
			int single  a_refRatio, int single  a_nComp, RectDomain<SPACE_DIM> single  a_domf, 
			BoxLayout local single a_layoutc, RectDomain<SPACE_DIM> single  a_domc){
		PrivateRegion local pr = new PrivateRegion(); 
	    	Region dr;
		m_layoutf=a_layoutf;
		//m_Iteratorf=new LayoutIterator(a_layoutf.getRoot());
		m_Iteratorf=new LayoutIterator(a_layoutf);
		m_DX=a_DX;
		m_refRatio=a_refRatio;
		//m_nComp=a_nComp;
		m_domf=a_domf;
		m_layoutc=a_layoutc;
		m_domc=a_domc;
		//m_phiBoundaryc=new template BoxLayoutData<template BoxedArray<T > > single local [2*SPACE_DIM] local;
		
		int single  orient,high_low,plus_minus,d;
		BoxLayout local single  boundaryf, boundaryc;
		//size=(int single)m_layoutf.size();
		int size=m_layoutf.numBoxesAt(Ti.thisProc());
		//m_boundaryf=new  RectDomain<SPACE_DIM> [2*SPACE_DIM] local [size] local;
		//m_boundaryc=new  RectDomain<SPACE_DIM> [2*SPACE_DIM] local [size] local;
		m_boundaryf=new  CFBoundary [2*SPACE_DIM] local [size] local;
		m_boundaryc=new  CFBoundary [2*SPACE_DIM] local [size] local;

		Domain<SPACE_DIM>   boxCut, boxCutc;
		RectDomain<SPACE_DIM > [] local   neighbors;
		DataIndex local DI;
		int single i,j,k;

		for (orient=2;orient<2*(SPACE_DIM+1);orient++){
			d=orient/2;
			high_low=orient%2;
			
			if (high_low==highBoundary){
				plus_minus=-1;
				boundaryf=m_layoutf.highBoundary(d,1,false);
				boundaryc=m_layoutf.highBoundary(d,1,true);
				//boundaryc=m_layoutf.highBoundary(d,1,false);
			    }
			else{
				plus_minus=1;
				boundaryf=m_layoutf.lowBoundary(d,1,false);
				boundaryc=m_layoutf.lowBoundary(d,1,true);
				//boundaryc=m_layoutf.highBoundary(d,1,false);
			}
			boundaryc=boundaryc.coarsen(m_refRatio);
			m_Iteratorf.reset();j=0;
			//for (k=0;k<size;k++){
			while ((boolean single)!m_Iteratorf.isEnded()){
			   DI=m_Iteratorf.getDataIndex();
			   boxCut=boundaryf[DI]; 
			   boxCutc=boundaryc[DI]; 
			   if ((boolean single)(boxCut<m_domf)){
				neighbors=m_layoutf.getNeighbors(DI);
	
				dr = Domain<SPACE_DIM>.setRegion(pr);
				for (i=0;i<(int single)neighbors.length-1;i++){
					boxCut=boxCut-neighbors[i];
					boxCutc=boxCutc-Util.coarsen(neighbors[i],m_refRatio);
				}
				Domain<SPACE_DIM>.setRegion(dr);
				
				boxCut= boxCut-neighbors[i];
				//boxCutc=(boxCutc-Util.coarsen(neighbors[i],m_refRatio));
				boxCutc=(boxCutc-Util.coarsen(neighbors[i],m_refRatio))*m_domc;
				
				//m_boundaryf[orient-2][j]=(RectDomain<SPACE_DIM>)boxCut;
				//m_boundaryc[orient-2][j]=(RectDomain<SPACE_DIM>)(boxCutc*m_domc);
				//if (!boxCut.isRectangular()) IsDomainCount++;
				if (boxCut.isRectangular())
					m_boundaryf[orient-2][j]=new CFBoundary((RectDomain<SPACE_DIM>)boxCut);
				else
					m_boundaryf[orient-2][j]=new CFBoundary(boxCut);
				if (boxCutc.isRectangular())
					m_boundaryc[orient-2][j]=new CFBoundary((RectDomain<SPACE_DIM>)boxCutc);
				else  
					m_boundaryc[orient-2][j]=new CFBoundary(boxCutc);
								
			   }
			   else{
				//m_boundaryf[orient-2][j]=NULL_BOX;
				m_boundaryf[orient-2][j]=new CFBoundary(NULL_BOX);
				//m_boundaryc[orient-2][j]=NULL_BOX;
				m_boundaryc[orient-2][j]=new CFBoundary(NULL_BOX);
			   }
			   //System.out.println("boxCutC before: "+boxCutc.toString()+"; boxCutc after: "+((RectDomain<SPACE_DIM>)boxCutc).toString());
			   m_Iteratorf.advance();j++;
			}

			//m_phiBoundaryc[orient-2]=new template BoxLayoutData<template BoxedArray<T > >((boundaryc.coarsen(m_refRatio),0);
			//m_phiBoundaryc[orient-2]=new template BoxLayoutData<template BoxedArray<T > >(boundaryc,0);
			//m_phiBoundaryc[orient-2].myprint();
		}
		//m_copier=new template Copier<template BoxedArray<T > > (m_phiBoundaryc.boxLayout(),m_layoutc);
		try {
		    pr.delete();
		} 
		catch (RegionInUse e){
			Util.printErrMsg("QuadCFInterp::define(): failed to delete PrivateRegion pr");
		}
		m_phiBoundaryc=new template BoxLayoutData<template BoxedArray<T > >(m_layoutf.coarsen(m_refRatio).accrete(1),0);
	}
	final local double quadPoly(double x, double y, double h,double bias){
		double z;
		double a=1;double b=1;double c=1;
		z=(x+y+1+bias)*h;
		return a*(z*z)+b*z+c;
	}
	/*//versin before Dec 10 
	public final local void coarseFineInterp(template BoxLayoutData<template BoxedArray<T > > local single a_phif){
		
		if (!a_phif.isDisjointed() || a_phif.boxLayout()!=m_layoutf){
			BoxTools.Util.printErrMsg("Invalid layout or boxes are not guaranteed to be disjointed!");
			return;
		}
		tmr1.reset();tmr1.start();
		int single  orient,high_low,plus_minus,d;
		//RectDomain<SPACE_DIM> boxf;
		Domain<SPACE_DIM> local boxf;
		Point<SPACE_DIM> p, pc,ed,e2d;
		DataIndex local DI;
		T [SPACE_DIM d] local TArrayf;
		T a,b,c;

		for (orient=2;orient<2*(SPACE_DIM+1);orient++){
			d=orient/2;
			plus_minus=1-(orient%2)*2;
			ed=Point<SPACE_DIM>.direction(d);
			e2d=Point<SPACE_DIM>.direction(d,2);
			m_Iteratorf.reset();
			while (!m_Iteratorf.isEnded()){
				DI=m_Iteratorf.getDataIndex();
				boxf=m_boundaryf[orient-2][DI.index()];
				TArrayf=a_phif[DI].getLocalArray();
				boxf=boxf*TArrayf.domain();

				foreach (p in boxf){				
					c=TArrayf[p+e2d*plus_minus];
					a=2*((m_refRatio+1)*c-(m_refRatio+3)*TArrayf[p+ed*plus_minus])/((m_refRatio+3)*(m_refRatio+1));
					b=TArrayf[p+ed*plus_minus]-a-c;
					TArrayf[p]=4*a+2*b+c;	
				}
						
				m_Iteratorf.advance();
			}
		}
		tmr1.stop();Tint1+=tmr1.millis();
	}*/
	
	public final local void coarseFineInterp(template BoxLayoutData<template BoxedArray<T > > local single a_phif){
		if (!a_phif.isDisjointed() || a_phif.boxLayout()!=m_layoutf){
			BoxTools.Util.printErrMsg("Invalid layout or boxes are not guaranteed to be disjointed!");
			return;
		}
		tmr1.reset();tmr1.start();
		
		PrivateRegion local region = new PrivateRegion();
		Region dr = Domain<SPACE_DIM>.setRegion(region);
		
		int single  orient,high_low,plus_minus,d;
		RectDomain<SPACE_DIM> boxf;
		Domain<SPACE_DIM>   domainf;
		RectDomain<SPACE_DIM> rectDomainf;
		//Point<SPACE_DIM> p, pc; /* for change in the syntax of foreach */
		Point<SPACE_DIM> pc;
		Point<SPACE_DIM> [1d] local ed =new(region) Point<SPACE_DIM> [[1:SPACE_DIM]];
		Point<SPACE_DIM> [1d] local e2d =new(region) Point<SPACE_DIM> [[1:SPACE_DIM]];
		int i;
		for (i=1;i<=SPACE_DIM;i++){ed[i]=Point<SPACE_DIM>.direction(i);e2d[i]=Point<SPACE_DIM>.direction(i,2);}
		DataIndex local DI;
		T [SPACE_DIM d] local TArrayf;
		T a,b,c;
		final int refRatio=m_refRatio;
		//int temp=0;tmr1.reset();tmr2.reset();tmrCopy.reset();
		boolean IsRectangular;

		m_Iteratorf.reset();i=0;
		while (!(boolean single)m_Iteratorf.isEnded()){
			//tmr1.start();
			DI=m_Iteratorf.getDataIndex();
			TArrayf=a_phif[DI].getLocalArray();
			//tmr1.stop();temp++;
			//rectDomainf=TArrayf.domain(); 
		 	for (orient=2;orient<2*(SPACE_DIM+1);orient++){
				
				d=orient/2;
				//ed=Point<SPACE_DIM>.direction(d);
				//e2d=Point<SPACE_DIM>.direction(d,2);
				plus_minus=1-(orient%2)*2;
				//domainf=m_boundaryf[orient-2][DI.index()]*TArrayf.domain();
				//domainf=m_boundaryf[orient-2][i]*TArrayf.domain();
				//tmrCopy.start();

				//domainf=m_boundaryf[orient-2][i]; //on sep 14, 2004
	
				//tmrCopy.stop();
				//domainf=domainf*rectDomainf;
				
				//tmr2.start();
				/* on sep 14, 2004
				foreach (p in domainf){				
					c=TArrayf[p+e2d[d]*plus_minus];
					a=2*((refRatio+1)*c-(refRatio+3)*TArrayf[p+ed[d]*plus_minus])/((refRatio+3)*(refRatio+1));
					b=TArrayf[p+ed[d]*plus_minus]-a-c;
					TArrayf[p]=4*a+2*b+c;	
				}
				*/
				//tmr2.stop();

				/* on sep 14, 2004 */
				if (m_boundaryf[orient-2][i].m_isRectangular){
					rectDomainf=m_boundaryf[orient-2][i].m_rectDomain;
					foreach (p in rectDomainf){				
					  c=TArrayf[p+e2d[d]*plus_minus];
					  a=2*((refRatio+1)*c-(refRatio+3)*TArrayf[p+ed[d]*plus_minus])/((refRatio+3)*(refRatio+1));
					  b=TArrayf[p+ed[d]*plus_minus]-a-c;
					  TArrayf[p]=4*a+2*b+c;	
					}
				}
				else{
					domainf=m_boundaryf[orient-2][i].m_domain;
					foreach (p in domainf){				
					  c=TArrayf[p+e2d[d]*plus_minus];
					  a=2*((refRatio+1)*c-(refRatio+3)*TArrayf[p+ed[d]*plus_minus])/((refRatio+3)*(refRatio+1));
					  b=TArrayf[p+ed[d]*plus_minus]-a-c;
					  TArrayf[p]=4*a+2*b+c;	
				        }
				}
				
			}
			m_Iteratorf.advance();i++;
		}
		//System.out.println("cost for CFInterp to get data on "+Ti.thisProc()+" is "+tmr1.secs());
		//System.out.println("cost for intersection "+Ti.thisProc()+" is "+tmrCopy.secs());
		//System.out.println("cost for CFInterp to interp on "+Ti.thisProc()+" is "+tmr2.secs()+" and # is "+temp);
		Domain<SPACE_DIM>.setRegion(dr);
		try{
		    region.delete();
		}
		catch(RegionInUse e){
			Util.printErrMsg("QuadCFInterp::coarseFineInterp1: failed to delete PrivateRegion region");
		}
		tmr1.stop();Tint1+=tmr1.millis();
	}
	
	
	public final local  void coarseFineInterp(template BoxLayoutData<template BoxedArray<T > > local single a_phif, template BoxLayoutData<template BoxedArray<T > > local single  a_phic){
		if (!(a_phif.isDisjointed() && a_phic.isDisjointed()) || !(a_phif.boxLayout()==m_layoutf && a_phic.boxLayout()==m_layoutc)){
			BoxTools.Util.printErrMsg("Invalid layout or boxes are not guaranteed to be disjointed!");
			return;
		}
		tmr2.reset();tmr2.start();
		
		PrivateRegion local region = new PrivateRegion();
		Region dr;
		
		int single  orient,high_low,plus_minus,d;
		int [] local   d1=new(region) int [SPACE_DIM-1];
		int   i,j,k;

		RectDomain<SPACE_DIM> boxf, boxc;
		Domain<SPACE_DIM>  domainf, domainc;
		
		//Point<SPACE_DIM> p, pc,ed,e2d,ed0,ed1,eD,e2D; /* for change in the syntax of foreach */
		Point<SPACE_DIM> pc,ed,e2d,ed0,ed1,eD,e2D;
		Point<SPACE_DIM> pc_old=Point<SPACE_DIM>.all(-1000000);
		
		DataIndex local DI;
		T [SPACE_DIM d] local TArrayf, TArrayc;
	
		j=(SPACE_DIM-1)*2+((SPACE_DIM<3)? 0:1);
		double [] local   dx=new(region) double[j];
		T [] local Diff=new(region) T[j];
		T phi,a,b,c;
		tmrCopy.reset();tmrCopy.start();m_phiBoundaryc.copy(a_phic,false);tmrCopy.stop(); TintCopy+=tmrCopy.millis();
		//tmrCopy.reset();tmrCopy.start();m_copier.copy(m_phiBoundaryc,a_phic);tmrCopy.stop(); TintCopy+=tmrCopy.millis();
		for (orient=2;orient<2*(SPACE_DIM+1);orient++){
			d=orient/2;
			high_low=orient%2;
			
			if (high_low==highBoundary) plus_minus=-1; else plus_minus=1;
			eD=Point<SPACE_DIM>.direction(d);
			e2D=Point<SPACE_DIM>.direction(d,2);

			//tmrCopy.reset();tmrCopy.start();m_phiBoundaryc[orient-2].copy(a_phic,false);tmrCopy.stop(); TintCopy+=tmrCopy.millis();//change to false on Oct 22, 03
			k=0;
			for (i=1;i<SPACE_DIM+1;i++) 
				if (i!=d){
					d1[k]=i;k++;
				}
			if (SPACE_DIM==3){
				ed0=Point<SPACE_DIM>.direction(d1[0]);
				ed1=Point<SPACE_DIM>.direction(d1[1]);
			}

			m_Iteratorf.reset();j=0;
			while (!m_Iteratorf.isEnded()){
				DI=m_Iteratorf.getDataIndex();
				//boxf=m_boundaryf[orient-2][DI.index()];
				//domainf=m_boundaryf[orient-2][j]; //changed on sep 15, 2004	
				TArrayf=a_phif[DI].getLocalArray();
				//TArrayc=m_phiBoundaryc[orient-2][DI].getLocalArray();
				TArrayc=m_phiBoundaryc[DI].getLocalArray();
				//System.out.println("box_"+j+": "+boxCut.toString());
				//boxc=TArrayc.domain()*m_domc*m_boundaryc[orient-2][DI.index()]; 
				//boxc=TArrayc.domain()*m_boundaryc[orient-2][DI.index()]; //change on Nov 24,03
				//boxc=m_boundaryc[orient-2][DI.index()];
				//domainc=m_boundaryc[orient-2][j];//changed on sep 15, 2004
				//boxf=boxf*TArrayf.domain();
				if (m_boundaryf[orient-2][j].m_isRectangular){
				   dr=Domain<SPACE_DIM>.setRegion(region);
				   boxf=m_boundaryf[orient-2][j].m_rectDomain;
				   foreach (p in boxf){
					pc=p/m_refRatio;
					phi=TArrayc[pc];	
										
					for (i=0;i<SPACE_DIM-1;i++) {
					   k=i*2;
					   dx[k]=(p[d1[i]]+0.5)/m_refRatio-(pc[d1[i]]+0.5);
					   dx[k+1]=dx[k]*dx[k]*0.5;
					   if (pc!=pc_old){
						ed=Point<SPACE_DIM>.direction(d1[i]);
						e2d=Point<SPACE_DIM>.direction(d1[i],2);
						Diff[k]=(TArrayc[pc+ed]-TArrayc[pc-ed])*0.5;
						Diff[k+1]=TArrayc[pc+ed]-2*TArrayc[pc]+TArrayc[pc-ed];   
					   }
					}
					if (SPACE_DIM==3){
				     		dx[4]=dx[0]*dx[2];						
				    		if (pc!=pc_old){
				        	 Diff[4]=TArrayc[pc]-TArrayc[pc];
					
						 //ed=Point<SPACE_DIM>.direction(d1[0]);
						 //ed1=Point<SPACE_DIM>.direction(d1[1]);
						 for (i=-1;i<2;i+=2)
							for (k=-1;k<2;k+=2){
						        	Diff[4]=Diff[4]+(i*k)*(TArrayc[pc+ed0*i+ed1*k]+TArrayc[pc]-TArrayc[pc+ed0*i]-TArrayc[pc+ed1*k]);
						    	}
						 Diff[4]=Diff[4]/4;
					   	}
					}
			
					for (i=0;i<dx.length;i++) phi=phi+Diff[i]*dx[i];
				
					c=TArrayf[p+e2D*plus_minus];	
					a=2*(2*phi-(m_refRatio+3)*TArrayf[p+eD*plus_minus]+(m_refRatio+1)*c)/((m_refRatio+3)*(m_refRatio+1));
					b=TArrayf[p+eD*plus_minus]-a-c;
					TArrayf[p]=4*a+2*b+c;
					pc_old=pc;

				   }
				  Domain<SPACE_DIM>.setRegion(dr);
				}
				else{
				dr=Domain<SPACE_DIM>.setRegion(region);
				domainf=m_boundaryf[orient-2][j].m_domain;
				
				if (m_boundaryc[orient-2][j].m_isRectangular){
				  boxc=m_boundaryc[orient-2][j].m_rectDomain;
				  foreach (p in domainf){
						pc=p/m_refRatio;
						phi=TArrayc[pc];
					
						for (i=0;i<SPACE_DIM-1;i++){
							k=i*2;
							dx[k]=(p[d1[i]]+0.5)/m_refRatio-(pc[d1[i]]+0.5);
							dx[k+1]=dx[k]*dx[k]*0.5;
						   if (pc!=pc_old){
							ed=Point<SPACE_DIM>.direction(d1[i]);
							e2d=Point<SPACE_DIM>.direction(d1[i],2);
							if (boxc.contains(pc+ed))
								if (boxc.contains(pc-ed)){
									Diff[k]=(TArrayc[pc+ed]-TArrayc[pc-ed])*0.5;
									Diff[k+1]=TArrayc[pc+ed]-2*TArrayc[pc]+TArrayc[pc-ed];
							   	}
								else 
								  if (boxc.contains(pc+e2d)){
									Diff[k]=(TArrayc[pc+ed]-TArrayc[pc])*1.5-(TArrayc[pc+e2d]-TArrayc[pc+ed])*0.5;
									Diff[k+1]=TArrayc[pc+e2d]-2*TArrayc[pc+ed]+TArrayc[pc];
								  }
								  else{
									Diff[k]=TArrayc[pc+ed]-TArrayc[pc];
									Diff[k+1]=0;
								  }							
							else
								if (boxc.contains(pc-ed))
								  if (boxc.contains(pc-e2d)){
									Diff[k]=(TArrayc[pc]-TArrayc[pc-ed])*1.5-(TArrayc[pc-ed]-TArrayc[pc-e2d])*0.5;
									Diff[k+1]=TArrayc[pc]-2*TArrayc[pc-ed]+TArrayc[pc-e2d];
							   	  }
								  else{
									Diff[k]=TArrayc[pc]-TArrayc[pc-ed];
									Diff[k+1]=0;	
								  }
								else{
									Diff[k]=0;
									Diff[k+1]=0;
								}
						   }
						}
					
					
					if (SPACE_DIM==3){
					     dx[4]=dx[0]*dx[2];						
					     if (pc!=pc_old){
					        Diff[4]=TArrayc[pc]-TArrayc[pc];
						int   n=0;
						//ed=Point<SPACE_DIM>.direction(d1[0]);
						//ed1=Point<SPACE_DIM>.direction(d1[1]);
						for (i=-1;i<2;i+=2)
							for (k=-1;k<2;k+=2)
							    if (boxc.contains(pc+ed0*i+ed1*k)&&boxc.contains(pc+ed0*i)&&boxc.contains(pc+ed1*k)){
							        n++;
							        Diff[4]=Diff[4]+(i*k)*(TArrayc[pc+ed0*i+ed1*k]+TArrayc[pc]-TArrayc[pc+ed0*i]-TArrayc[pc+ed1*k]);
							    }
						if (n>0) Diff[4]=Diff[4]/n;
					     }

					}
					
					for (i=0;i<dx.length;i++) phi=phi+Diff[i]*dx[i];
				
					c=TArrayf[p+e2D*plus_minus];	
					a=2*(2*phi-(m_refRatio+3)*TArrayf[p+eD*plus_minus]+(m_refRatio+1)*c)/((m_refRatio+3)*(m_refRatio+1));
					b=TArrayf[p+eD*plus_minus]-a-c;
					TArrayf[p]=4*a+2*b+c;
					pc_old=pc;
				  }
				}
				else{
				  domainc=m_boundaryc[orient-2][j].m_domain;
				  foreach (p in domainf){
						pc=p/m_refRatio;
						phi=TArrayc[pc];
					
						for (i=0;i<SPACE_DIM-1;i++){
							k=i*2;
							dx[k]=(p[d1[i]]+0.5)/m_refRatio-(pc[d1[i]]+0.5);
							dx[k+1]=dx[k]*dx[k]*0.5;
						   if (pc!=pc_old){
							ed=Point<SPACE_DIM>.direction(d1[i]);
							e2d=Point<SPACE_DIM>.direction(d1[i],2);
							if (domainc.contains(pc+ed))
								if (domainc.contains(pc-ed)){
									Diff[k]=(TArrayc[pc+ed]-TArrayc[pc-ed])*0.5;
									Diff[k+1]=TArrayc[pc+ed]-2*TArrayc[pc]+TArrayc[pc-ed];
							   	}
								else 
								  if (domainc.contains(pc+e2d)){
									Diff[k]=(TArrayc[pc+ed]-TArrayc[pc])*1.5-(TArrayc[pc+e2d]-TArrayc[pc+ed])*0.5;
									Diff[k+1]=TArrayc[pc+e2d]-2*TArrayc[pc+ed]+TArrayc[pc];
								  }
								  else{
									Diff[k]=TArrayc[pc+ed]-TArrayc[pc];
									Diff[k+1]=0;
								  }							
							else
								if (domainc.contains(pc-ed))
								  if (domainc.contains(pc-e2d)){
									Diff[k]=(TArrayc[pc]-TArrayc[pc-ed])*1.5-(TArrayc[pc-ed]-TArrayc[pc-e2d])*0.5;
									Diff[k+1]=TArrayc[pc]-2*TArrayc[pc-ed]+TArrayc[pc-e2d];
							   	  }
								  else{
									Diff[k]=TArrayc[pc]-TArrayc[pc-ed];
									Diff[k+1]=0;	
								  }
								else{
									Diff[k]=0;
									Diff[k+1]=0;
								}
						   }
						}
					
					
					if (SPACE_DIM==3){
					     dx[4]=dx[0]*dx[2];						
					     if (pc!=pc_old){
					        Diff[4]=TArrayc[pc]-TArrayc[pc];
						int   n=0;
						//ed=Point<SPACE_DIM>.direction(d1[0]);
						//ed1=Point<SPACE_DIM>.direction(d1[1]);
						for (i=-1;i<2;i+=2)
							for (k=-1;k<2;k+=2)
							    if (domainc.contains(pc+ed0*i+ed1*k)&&domainc.contains(pc+ed0*i)&&domainc.contains(pc+ed1*k)){
							        n++;
							        Diff[4]=Diff[4]+(i*k)*(TArrayc[pc+ed0*i+ed1*k]+TArrayc[pc]-TArrayc[pc+ed0*i]-TArrayc[pc+ed1*k]);
							    }
						if (n>0) Diff[4]=Diff[4]/n;
					     }

					}
					
					for (i=0;i<dx.length;i++) phi=phi+Diff[i]*dx[i];
				
					c=TArrayf[p+e2D*plus_minus];	
					a=2*(2*phi-(m_refRatio+3)*TArrayf[p+eD*plus_minus]+(m_refRatio+1)*c)/((m_refRatio+3)*(m_refRatio+1));
					b=TArrayf[p+eD*plus_minus]-a-c;
					TArrayf[p]=4*a+2*b+c;
					pc_old=pc;
				  }	
				}
				Domain<SPACE_DIM>.setRegion(dr);
				}
				m_Iteratorf.advance();j++;
			   
			}
		}
		
		try{
		    region.delete();
		}
		catch(RegionInUse e){
			Util.printErrMsg("QuadCFInterp::coarseFineInterp2: failed to delete PrivateRegion region");
		}
		
		tmr2.stop();Tint2+=tmr2.millis();
	}

	/*//version before Dec 10, 03
	public final local  void coarseFineInterp(template BoxLayoutData<template BoxedArray<T > > local single a_phif, template BoxLayoutData<template BoxedArray<T > > local single  a_phic){
		if (!(a_phif.isDisjointed() && a_phic.isDisjointed()) || !(a_phif.boxLayout()==m_layoutf && a_phic.boxLayout()==m_layoutc)){
			BoxTools.Util.printErrMsg("Invalid layout or boxes are not guaranteed to be disjointed!");
			return;
		}
		tmr2.reset();tmr2.start();
		int single  orient,high_low,plus_minus,d;
		int [] local   d1=new int [SPACE_DIM-1];
		int   i,j,k;

		RectDomain<SPACE_DIM> boxf, boxc;
		Domain<SPACE_DIM>  domainf, domainc;
		
		Point<SPACE_DIM> p, pc,ed,e2d,ed0,ed1,eD,e2D;
		Point<SPACE_DIM> pc_old=Point<SPACE_DIM>.all(-1000000);
		
		DataIndex local DI;
		T [SPACE_DIM d] local TArrayf, TArrayc;
	
		j=(SPACE_DIM-1)*2+((SPACE_DIM<3)? 0:1);
		double [] local   dx=new double[j];
		T [] local Diff=new T[j];
		T phi,a,b,c;
		tmrCopy.reset();tmrCopy.start();m_phiBoundaryc.copy(a_phic,false);tmrCopy.stop(); TintCopy+=tmrCopy.millis();
		//tmrCopy.reset();tmrCopy.start();m_copier.copy(m_phiBoundaryc,a_phic);tmrCopy.stop(); TintCopy+=tmrCopy.millis();
		m_Iteratorf.reset();
		while (!(boolean single)m_Iteratorf.isEnded()){
			DI=m_Iteratorf.getDataIndex();
			TArrayf=a_phif[DI].getLocalArray();
			TArrayc=m_phiBoundaryc[DI].getLocalArray();
	
			for (orient=2;orient<2*(SPACE_DIM+1);orient++){
				d=orient/2;
				plus_minus=1-(orient%2)*2;

				eD=Point<SPACE_DIM>.direction(d);
				e2D=Point<SPACE_DIM>.direction(d,2);
			
				k=0;
				for (i=1;i<SPACE_DIM+1;i++) 
					if (i!=d){
						d1[k]=i;k++;
					}
				if (SPACE_DIM==3){
					ed0=Point<SPACE_DIM>.direction(d1[0]);
					ed1=Point<SPACE_DIM>.direction(d1[1]);
				}			
				
				domainf=m_boundaryf[orient-2][DI.index()];	
				//domainc=TArrayc.domain()*m_boundaryc[orient-2][DI.index()]; //change on Nov 24,03
				domainc=m_boundaryc[orient-2][DI.index()];
				domainf=domainf*TArrayf.domain();
				if (domainf.isRectangular()){
				   boxf=(RectDomain<SPACE_DIM>)domainf;
				   //for (i=0;i<k;i++) {boxf=boxf.shrink(m_refRatio,d1[i]);boxf=boxf.shrink(m_refRatio,-d1[i]);}
				   //domainf=domainf-boxf;
				   foreach (p in boxf){
					pc=p/m_refRatio;
					phi=TArrayc[pc];	
										
					for (i=0;i<SPACE_DIM-1;i++) {
					   k=i*2;
					   dx[k]=(p[d1[i]]+0.5)/m_refRatio-(pc[d1[i]]+0.5);
					   dx[k+1]=dx[k]*dx[k]*0.5;
					   if (pc!=pc_old){
						ed=Point<SPACE_DIM>.direction(d1[i]);
						e2d=Point<SPACE_DIM>.direction(d1[i],2);
						Diff[k]=(TArrayc[pc+ed]-TArrayc[pc-ed])*0.5;
						Diff[k+1]=TArrayc[pc+ed]-2*TArrayc[pc]+TArrayc[pc-ed];   
					   }
					}
					if (SPACE_DIM==3){
				     		dx[4]=dx[0]*dx[2];						
				    		if (pc!=pc_old){
				        	 Diff[4]=TArrayc[pc]-TArrayc[pc];
					
						 //ed=Point<SPACE_DIM>.direction(d1[0]);
						 //ed1=Point<SPACE_DIM>.direction(d1[1]);
						 for (i=-1;i<2;i+=2)
							for (k=-1;k<2;k+=2){
						        	Diff[4]=Diff[4]+(i*k)*(TArrayc[pc+ed0*i+ed1*k]+TArrayc[pc]-TArrayc[pc+ed0*i]-TArrayc[pc+ed1*k]);
						    	}
						 Diff[4]=Diff[4]/4;
					   	}
					}
			
					for (i=0;i<dx.length;i++) phi=phi+Diff[i]*dx[i];
				
					c=TArrayf[p+e2D*plus_minus];	
					a=2*(2*phi-(m_refRatio+3)*TArrayf[p+eD*plus_minus]+(m_refRatio+1)*c)/((m_refRatio+3)*(m_refRatio+1));
					b=TArrayf[p+eD*plus_minus]-a-c;
					TArrayf[p]=4*a+2*b+c;
					pc_old=pc;

				   }
				}
				else{
				foreach (p in domainf){
					pc=p/m_refRatio;
					phi=TArrayc[pc];
					
						for (i=0;i<SPACE_DIM-1;i++){
							k=i*2;
							dx[k]=(p[d1[i]]+0.5)/m_refRatio-(pc[d1[i]]+0.5);
							dx[k+1]=dx[k]*dx[k]*0.5;
						   if (pc!=pc_old){
							ed=Point<SPACE_DIM>.direction(d1[i]);
							e2d=Point<SPACE_DIM>.direction(d1[i],2);
							if (domainc.contains(pc+ed))
								if (domainc.contains(pc-ed)){
									Diff[k]=(TArrayc[pc+ed]-TArrayc[pc-ed])*0.5;
									Diff[k+1]=TArrayc[pc+ed]-2*TArrayc[pc]+TArrayc[pc-ed];
							   	}
								else 
								  if (domainc.contains(pc+e2d)){
									Diff[k]=(TArrayc[pc+ed]-TArrayc[pc])*1.5-(TArrayc[pc+e2d]-TArrayc[pc+ed])*0.5;
									Diff[k+1]=TArrayc[pc+e2d]-2*TArrayc[pc+ed]+TArrayc[pc];
								  }
								  else{
									Diff[k]=TArrayc[pc+ed]-TArrayc[pc];
									Diff[k+1]=0;
								  }							
							else
								if (domainc.contains(pc-ed))
								  if (domainc.contains(pc-e2d)){
									Diff[k]=(TArrayc[pc]-TArrayc[pc-ed])*1.5-(TArrayc[pc-ed]-TArrayc[pc-e2d])*0.5;
									Diff[k+1]=TArrayc[pc]-2*TArrayc[pc-ed]+TArrayc[pc-e2d];
							   	  }
								  else{
									Diff[k]=TArrayc[pc]-TArrayc[pc-ed];
									Diff[k+1]=0;	
								  }
								else{
									Diff[k]=0;
									Diff[k+1]=0;
								}
						   }
						}
					
					
					if (SPACE_DIM==3){
					     dx[4]=dx[0]*dx[2];						
					     if (pc!=pc_old){
					        Diff[4]=TArrayc[pc]-TArrayc[pc];
						int   n=0;
						//ed=Point<SPACE_DIM>.direction(d1[0]);
						//ed1=Point<SPACE_DIM>.direction(d1[1]);
						for (i=-1;i<2;i+=2)
							for (k=-1;k<2;k+=2)
							    if (domainc.contains(pc+ed0*i+ed1*k)&&domainc.contains(pc+ed0*i)&&domainc.contains(pc+ed1*k)){
							        n++;
							        Diff[4]=Diff[4]+(i*k)*(TArrayc[pc+ed0*i+ed1*k]+TArrayc[pc]-TArrayc[pc+ed0*i]-TArrayc[pc+ed1*k]);
							    }
						if (n>0) Diff[4]=Diff[4]/n;
					     }

					}
					
					for (i=0;i<dx.length;i++) phi=phi+Diff[i]*dx[i];
				
					c=TArrayf[p+e2D*plus_minus];	
					a=2*(2*phi-(m_refRatio+3)*TArrayf[p+eD*plus_minus]+(m_refRatio+1)*c)/((m_refRatio+3)*(m_refRatio+1));
					b=TArrayf[p+eD*plus_minus]-a-c;
					TArrayf[p]=4*a+2*b+c;
					pc_old=pc;
				
			   	}
				}
			}
			m_Iteratorf.advance();
		}
		tmr2.stop();Tint2+=tmr2.millis();
	}
*/

}