public class MeshRefineD {

  // constants

  private static final int dimensions = 2;
  private static final RectDomain<1> Ddims = [1:dimensions], Ddirs = [-1:1:2];

  // non-constant fields

  private static Domain<2> [1d] ProperNestingDomain; // [base:top]
  private static RectDomain<2> [1d] PhysicalDomain;  // [base:top]

  // Methods in this class also use the following fields from Amr2:
  // domain, nRefine, nRefineP

  public static final RectDomain<2> [1d][1d]  // [base+1:top+1][1d] returned
    Refine(RectDomain<2> [1d][1d] oldGrids,   // [base:top][1d]
	   Domain<2> [1d] tagged) {        // [base:top]

    // Given grids at levels base:top
    // and domains containing tagged cells at levels base:top,
    // return arrays of rectangles that cover them at levels base+1:top+1.

    RectDomain<1> levels = oldGrids.domain();
    int base = levels.min()[1], top = levels.max()[1];
    // System.out.println(base + " to " + top);

    // (1) Compute proper nesting domains for levels base:top
    // private static Domain<2> [1d] ProperNestingDomain;
    // System.out.println("Get physical domains");
    GetPhysicalDomains(base, top);

    // System.out.println("Get PNDs");
    GetPNDs(oldGrids);

    // (2) Turn off tags at points outside PND.
    // System.out.println("Clip tagged");
    ClipTagged(tagged);

    RectDomain<2> [1d][1d] returnGrids = new RectDomain<2>[base+1:top+1][1d];

    for (int h = top; h >= base; h--) {

      // (3) Find grids covering tagged points in patches at level h.

      // System.out.println("Make grids");
      RectDomain<2> [1d] newGrids =
	GridGenerator.MakeGrids(tagged[h], ProperNestingDomain[h]);

      // (4) Refine grids and store in returnGrids[h+1].
      
      // System.out.println("Refine grids");
      RefineGrids(newGrids);
      returnGrids[h+1] = new RectDomain<2>[newGrids.domain()];
      foreach (indgrid in newGrids.domain())
	returnGrids[h+1][indgrid] = newGrids[indgrid];

      // (5) Expand grids by one cell in each direction (clipping at
      // physical boundaries) and then coarsen and union with tags
      // at level h-1.
      // Omit this step for the base level.

      if (h > base) {
	// System.out.println("Coarsen etc.");
	CoarsenGrids(newGrids);  // back to level h

	// Extend each grid one cell in each direction.
	foreach (indgrid in newGrids.domain())
	  foreach (pdim in Ddims) foreach (pdir in Ddirs) {
	  Point<2> shift = Point<2>.direction(pdim[1] * pdir[1], 1);
	  newGrids[indgrid] = (RectDomain<2>)
	    (newGrids[indgrid] +
	     ((newGrids[indgrid] + shift) * PhysicalDomain[h]));
	}

	CoarsenGrids(newGrids);  // to level h-1

	foreach (indgrid in newGrids.domain())
	  tagged[h-1] = tagged[h-1] + newGrids[indgrid];
	// are these tagged at the right level?
      }
    }
    return returnGrids;
  }


  public static final void GetPhysicalDomains(int base, int top) {
    // Find physical domain for each level.

    PhysicalDomain = new RectDomain<2>[base:top];
    Point<2> zero = Point<2>.all(0);
    Point<2> one = Point<2>.all(1);
    Point<2> extent = Amr2.domain.max() + one;
    for (int h = 0; h > base; h--) extent = extent / Amr2.nRefineP;
    for (int h = 0; h < base; h++) extent = extent * Amr2.nRefineP;
    for (int h = base; h <= top; h++) {
      PhysicalDomain[h] = [zero : extent - one];
      // System.out.println("PhysicalDomain[" + h + "] = " +
      // PhysicalDomain[h].toString());
      extent = extent * Amr2.nRefineP;
    }
  }


  public static final void GetPNDs(RectDomain<2> [1d][1d] oldGrids) {

    // Compute proper nesting domains.

    ProperNestingDomain = new Domain<2>[oldGrids.domain()];

    foreach (h in oldGrids.domain()) {

      // First find the domain covered by grids at level h.

      // System.out.println("GetPND for " + h[1]);
      Domain<2> levDomain = [Point<2>.all(0) : Point<2>.all(-1)];
      foreach (indgrid in oldGrids[h].domain())
        levDomain = levDomain + oldGrids[h][indgrid];

      Domain<2> PND = levDomain;
      foreach (pdim in Ddims) foreach (pdir in Ddirs) {
        // extend one cell in each direction
        Point<2> shift = Point<2>.direction(pdim[1] * pdir[1], 1);

	Domain<2> remove = (((PND + shift) - PND) * PhysicalDomain[h]) - shift;
	// System.out.println("removing " + remove);
        PND = PND - remove;
      }
    
      ProperNestingDomain[h] = PND;
      // System.out.println("ProperNestingDomain[" + h[1] + "] = " + PND);
    }
  }


  public static final void ClipTagged(Domain<2> [1d] tagged) {

    // At each level, remove points outside proper nesting domain
    // from the tagged domain,

    foreach (h in tagged.domain())
      tagged[h] = tagged[h] * ProperNestingDomain[h];
  }


  private static final void RefineGrids(RectDomain<2> [1d] grids) {
    Point<2> one = Point<2>.all(1);
    foreach (ind in grids.domain())
      grids[ind] = 
        [grids[ind].min() * Amr2.nRefineP :
	 (grids[ind].max() + one) * Amr2.nRefineP - one];
  }


  private static final void CoarsenGrids(RectDomain<2> [1d] grids) {
    foreach (ind in grids.domain()) grids[ind] = grids[ind] / Amr2.nRefineP;
  }
}