public class MeshRefine {

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

  // non-constant fields

  private static int single baselev, toplev;

  private static RectDomain<1> single levinds, levindsnew;

  // array of domains of the patches at each height

  private static RectDomain<2> [1d][1d] single oldGrids;

  private static Domain<2> [1d] single tagged;  // baselev:toplev

  single public static final void Refine(int single base, int single top) {

    // Call this routine from main program.
    // Replace the patches at heights base+1:top
    // and add new patches at height top+1.
    // (If base == top then no patches are replaced; a new level is added.)

    baselev = base;
    toplev = top;
    levinds = [base : top];
    levindsnew = [base+1 : top+1];

    // Get single list of patches and tagged cells at each height.

    // System.out.println("MeshRefine.GetOldGrids");
    GetOldGrids();   // in parallel

    int myProc = Ti.thisProc();
    RectDomain<1> single procTeam = [0 : Ti.numProcs() - 1];

    // Get new grids, for levels [base+1:top+1].
    
    RectDomain<2> [1d][1d][1d] myRecList =
      broadcast new RectDomain<2>[levindsnew][procTeam][1d] from 0;

    // System.out.println("defined myRecList");

    Ti.barrier();

    if (Ti.thisProc() == 0) {

      // System.out.println("setting oldGridsLoc");

      RectDomain<2> [1d][1d] oldGridsLoc = new RectDomain<2>[levinds][1d];
      foreach (h in levinds) oldGridsLoc[h] = oldGrids[h];

      Domain<2> [1d] taggedLoc = tagged;

      // Coarsen oldGridsLoc and taggedLoc, from [base:top] to [base-1:top-1].
      oldGridsLoc = CoarsenBoxes(oldGridsLoc);
      taggedLoc = CoarsenDomains(taggedLoc);

      // newGrids indexed by [base:top][1d]
      RectDomain<2> [1d][1d] newGrids =
	MeshRefineD.Refine(oldGridsLoc, taggedLoc);

      // Refine newGrids, from [base:top] to [base+1:top+1].
      newGrids = RefineBoxes(newGrids);

      // OLD AssignBoxes(newGrids, myRecList);
      foreach (h in levindsnew) {
	myRecList[h] = AssignBoxes(newGrids[h], h[1]);
      }
    }
    Ti.barrier();

    // System.out.println("define recListProc");

    RectDomain<2> [1d][1d] recListProc;
    recListProc = new RectDomain<2>[levindsnew][1d];

    // System.out.println("entering loop to assign recListProc");
    Ti.barrier();
    for (int single height = base + 1; height <= top + 1; height++) {
      recListProc[height] = myRecList[height][myProc];
    }
      
    for (int single height = base; height <= top; height++) {

      /*
	Assign a new Level at height+1 for each proc.
      */

      // Retrieve fields of Level Amr2.processes[myProc].levels[height].
      Level levelh = Amr2.processes[myProc].levels[height];
      RectDomain<2> D = levelh.domain;
      Point<2> extent = D.max() - D.min() + Point<2>.all(1);
      double [1d] dsh = levelh.ds;

      // For level height+1, expand the domain and shrink ds.
      extent = extent * Amr2.nRefine;
      D = [Point<2>.all(0) : extent - Point<2>.all(1)];
      double [1d] ds = new double[Ddims];
      foreach (pdim in Ddims) ds[pdim] = dsh[pdim] / Amr2.dRefine;

      // New Level for height+1.
      Level levelNew = new Level(recListProc[height+1].domain(), D, height+1, ds);
      if (Amr2.finestLevel < height + 1) Amr2.finestLevel = height + 1;
	
      Amr2.processes[myProc].levels[height+1] = levelNew;

      /*
	Assign patches to the new Level.
      */

      foreach (indpatch in recListProc[height+1].domain()) {
	levelNew.patches[indpatch] =
	  new Patch2(levelNew, recListProc[height+1][indpatch]);
	// Fill in other patch data later.
      }
    }
  }


  single public static final void GetOldGrids() {
    // Retrieve single array of grids at each height, and tagged domains.
    // You need to go into each process explicitly.
    oldGrids = new RectDomain<2>[levinds][1d] single;
    tagged = new Domain<2>[baselev:toplev] single;

    int myProc = Ti.thisProc();
    Amr2Process proc = Amr2.processes[myProc];
    RectDomain<1> single procTeam = [0 : Ti.numProcs() - 1];

    // DomainUnion unite = new DomainUnion(); // doesn't compile for NOW

    foreach (h in levinds) {

      Level lev = proc.levels[h];
      int mynumpatches = lev.indices.size();

      int single totpatches = Reduce.add(mynumpatches);
      // System.out.println("Level " + h[1] + ":  totpatches = " + totpatches);

      // oldGrids = new RectDomain<2>[levinds][1d] single;

      // oldGrids[h] =
      // (RectDomain<2> [1d] single) new RectDomain<2>[0:totpatches-1];
      oldGrids[h] = broadcast new RectDomain<2>[0:totpatches-1] from 0;

      // Scan.add(mynumpatches) gives total number of patches in processes
      // up to AND INCLUDING this one
      int indgrid = Scan.add(mynumpatches) - mynumpatches;

      // System.out.println("now setting old grids and tags");

      Domain<2> mytaggedh = [Point<2>.all(0) : Point<2>.all(-1)];

      foreach (indpatch in lev.indices) {
	Patch2 patch = lev.patches[indpatch];
	oldGrids[h][indgrid] = patch.domain;
	// System.out.println("oldGrids adding " + patch.domain + " at " + indgrid);
	mytaggedh = mytaggedh + patch.tagged;
	indgrid++;
      }
      // tagged[h] = (Domain<2>) Reduce.gen(unite, mytaggedh);

      Domain<2> [1d] single taggedall =
	(Domain<2> [1d] single) new Domain<2>[procTeam];
      taggedall.exchange(mytaggedh);
      Domain<2> taggedunion = [Point<2>.all(0) : Point<2>.all(-1)];
      if (myProc == 0) {
	foreach (proc in procTeam) taggedunion = taggedunion + taggedall[proc];
      }
      Ti.barrier();
      tagged[h] = broadcast taggedunion from 0;
    }
  }


  private static final RectDomain<2> [1d][1d]
    RefineBoxes(RectDomain<2> [1d][1d] boxes) {
    RectDomain<2> [1d][1d] newboxes = boxes.translate([+1]);
    foreach (h in newboxes.domain()) RefineBoxes(newboxes[h]);
    return newboxes;
  }


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


  private static final RectDomain<2> [1d][1d]
    CoarsenBoxes(RectDomain<2> [1d][1d] boxes) {
    RectDomain<2> [1d][1d] newboxes = boxes.translate([-1]);
    foreach (h in newboxes.domain()) CoarsenBoxes(newboxes[h]);
    return newboxes;
  }


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


  private static final Domain<2> [1d] CoarsenDomains(Domain<2> [1d] domains) {
    Domain<2> [1d] newdomains = domains.translate([-1]);
    foreach (h in newdomains.domain()) {
      newdomains[h] = newdomains[h] / Amr2.nRefineP;
      // The line above wouldn't give correct result earlier.  See Bug#146.
      // Domain<2> tempdomains = newdomains[h] / Amr2.nRefineP;
      // newdomains[h] = [Point<2>.all(0) : Point<2>.all(-1)];
      // foreach (p in tempdomains) newdomains[h] = newdomains[h] + [p:p];
    }
    return newdomains;
  }


  // OLD private static final void AssignBoxes(RectDomain<2> [1d][1d] newGrids,
  // RectDomain<2> [1d][1d][1d] myRecList) {
  private static final RectDomain<2> [1d][1d]
    AssignBoxes(RectDomain<2> [1d] boxes, int myh) {

    // Assign boxes to procs.

    int numProcs = Ti.numProcs();
    RectDomain<1> procTeam = [0 : numProcs - 1];
    RectDomain<2> [1d][1d] myRecList = new RectDomain<2>[procTeam][1d];

    /*
      procLoad[proc] is total load assigned to process proc

      The load associated with a grid is taken to be 
      grid.accrete(2) * levDomain
                   ^
                   2 ghost cells in each direction
    */
    int [1d] procLoad = new int[procTeam];
    foreach (proc in procTeam) procLoad[proc] = 0;

    Point<2> extent = Amr2.domain.max() - Amr2.domain.min() + Point<2>.all(1);
    for (int h = 0; h < myh; h++) extent = extent * Amr2.nRefine;
    RectDomain<2> levDomain = [Point<2>.all(0) : extent - Point<2>.all(1)];
    System.out.println("levDomain = " + levDomain);

    // Compute estimated load of each new patch.
    RectDomain<1> gridindices = boxes.domain();
    int [1d] load = new int[gridindices];

    foreach (indgrid in gridindices) {
      load[indgrid] = (boxes[indgrid].accrete(2) * levDomain).size();
      System.out.println("Load of " + boxes[indgrid] + " is " + load[indgrid]);
    }

    /*
      Figure out how to allocate to processes.
      This algorithm comes from
      William Y. Crutchfield, ``Load Balancing Irregular Algorithms''
    */

    Domain<1> unassigned = gridindices;
    Domain<1> [1d] procGrids = new Domain<1>[procTeam];
    // procGrids[proc] containing ind
    // means that boxes[ind] is assigned to process proc.
    foreach (proc in procTeam) procGrids[proc] = [0 : -1];

    // Find heaviest grid and assign it to the lightest process.
    while (! unassigned.isNull()) {
      // find heaviest grid among those that are as yet unassigned
      int maxloadGrid = 0;
      Point<1> heaviestGrid;
      foreach (gridp in unassigned)
	if (load[gridp] > maxloadGrid) {
	  heaviestGrid = gridp; maxloadGrid = load[gridp];
	}

      // find lightest process
      int lightestProc = 0;
      foreach (proc in procTeam)
	if (procLoad[proc] < procLoad[lightestProc]) lightestProc = proc[1];

      // assign heaviest grid to lightest process
      Domain<1> heaviestGridD = [heaviestGrid : heaviestGrid];
      procGrids[lightestProc] = procGrids[lightestProc] + heaviestGridD;
      procLoad[lightestProc] += load[heaviestGrid];
      unassigned = unassigned - heaviestGridD;
    }

    // System.out.println("now finding heaviest");

    boolean unbalanced = true;
    findswaps:
    while (unbalanced) {
      unbalanced = false;

      // procGridsPrint(procGrids);

      // find heaviest process
      int heaviestProc = 0;
      foreach (proc in procTeam)
	if (procLoad[proc] > procLoad[heaviestProc]) heaviestProc = proc[1];

      // for each grid in heaviest process
      foreach (gridph in procGrids[heaviestProc]) {
	int wth = load[gridph];

	// for each grid NOT in heaviest process
	foreach (otherProc in procTeam - [heaviestProc:heaviestProc]) {
	  Domain<1> procGridsOther = procGrids[otherProc];
	  foreach (gridpo in procGridsOther) {

	    int wto = load[gridpo];  // weight of the other grid
	    int loadhnew = procLoad[heaviestProc] - wth + wto;
	    int loadonew = procLoad[otherProc] - wto + wth;

	    // if interchanging gridph and gridpo improves load balance
	    // (i.e., difference in process loads goes down)
	    if (Math.abs(loadhnew - loadonew) <
		procLoad[heaviestProc] - procLoad[otherProc]) {

	      // System.out.println("interchanging " + gridph[1] +
	      // " in [" + heaviestProc + "] with " +
	      //		 gridpo[1] + " in " + otherProc);

	      // perform interchange
	      procGrids[heaviestProc] =
		procGrids[heaviestProc] - [gridph:gridph] + [gridpo:gridpo];
	      procLoad[heaviestProc] = loadhnew;
	      procGrids[otherProc] =
		procGrids[otherProc] - [gridpo:gridpo] + [gridph:gridph];
	      procLoad[otherProc] = loadonew;
	      // find heaviest process again...
	      unbalanced = true;
	      continue findswaps;
	    }
	  }
	}
      }
    } // while (unbalanced)

    // print out loads
    String s = "Loads:";
    for (int proc = 0; proc < Ti.numProcs(); proc++) {
      s = s + "  " + procLoad[proc] +
	" (" + procGrids[proc].size() + " grids)";
    }
    System.out.println(s);
    
    // Now assign grids to processes.
    foreach (proc in procTeam) {
      Domain<1> procGridsProc = procGrids[proc];
      int numgrids = procGridsProc.size();
      // assign domains
      myRecList[proc] = new RectDomain<2>[0 : numgrids - 1];
      int igrid = 0;
      foreach (gridp in procGridsProc) {
	myRecList[proc][igrid] = boxes[gridp];
	igrid++;
      }
      // should have igrid == numgrids now
    }
    return myRecList;
  }


  public static void procGridsPrint(Domain<1> [1d] procGrids) {
    String s = "procGrids: ";
    foreach (p in procGrids.domain()) {
      Domain<1> grids = procGrids[p];
      s = s + ' ' + p;
      foreach (g in grids) s = s + ' ' + g[1];
    }
    System.out.println(s);
  }
}