// Coarse-fine boundary patches require additional information
// used by the quadratic interpolation in the coarse-fine BC
// calculation.

import java.io.FileOutputStream;
import java.io.FileInputStream;
import java.io.PrintStream;
import java.io.IOException;
import java.io.StreamTokenizer;

// Temporary.
import ti.domains.*;
// Temporary end.

class CoarseFineBoundary {
    Domain<3> coarseBoundary;		// Boundary on coarse side.
    Patch [1d] coarsePatches;		// Coarse patches at boundary.

    // Coarse phi or error.  The domain is defined on a
    // tangentially extended boundary (to allow for quadratic
    // interpolation) and also accreted in the inward normal
    // direction, to optionally store additional values for the
    // flux computation.
    double [3d] coarseValues;

    // Some coarse cells in the extended boundary may be
    // covered by the fine grid.  The values on those cells
    // are computed from a corrected fine grid average
    // (where possible).  FINEPATCHES contains the
    // contributing patches.
    Patch [1d] finePatches;

    // FINEVALUES collects fine values for the averaging on
    // covered coarse cells.
    double [3d] fineValues;

    // COVERINGFINE is the extent of fine patches covering
    // COARSEVALUES.
    Domain<3> coveringFine;

    // OK to compute the laplacian of FINEVALUES in
    // INSIDELAP.
    Domain<3> insideLap;

    // NUMLAPS is the number of fine laplacians available
    // per coarse cell.
    double [3d] numLaps;

    static Domain<3> expand(Domain<3> x, Point<3> factor) {

	// This is the ``bad'' version.
	// Domain<3> r = (Domain<3>) ([[0, 0, 0] : [-1, -1, -1]]);
	// Domain<3> e = x * factor;
	// foreach (i in [[0, 0, 0] : factor - [1, 1, 1]]) {
	// r = r + (e + i);
	// }

	// RectDomain<3> [1d] dlist = x.RectDomainList();
	RectDomain<3> [1d] dlist = x.RectDomainList();
	Domain<3> r = (Domain<3>) ([[0, 0, 0] : [-1, -1, -1]]);
	foreach (i in dlist.domain()) {
	    Point<3> l = dlist[i].min() * factor;
	    // My pet peeve.
	    Point<3> u = (dlist[i].max() + [1, 1, 1]) * factor - [1, 1, 1];
	    // dlist[i] = [l : u];
	    r = r + [l : u];
	}

	// return Domain<3>.toDomain(dlist);
	return r;

	////	// Temporary fix for RectDomainList.  This will work
	////	// only with rectangular domains, so better check.
	////	RectDomain<3> box = [ x.min() : x.max() ];
	////	if (! (box - x).isEmpty()) {
	////		System.out.println("IRREGULAR DOMAIN IN EXPAND!");
	////	}
	////
	////	return (Domain<3>) [ x.min() * factor : x.max() * factor ];
    }

    // CFB is the coarse side of the boundary in coarse
    // coordinates.  ACFB is the tangentially accreted CFB.
    // PATCHES is a list of coarse patches that intersect
    // ACFB.  L is the index of this level.  DIR and SIDE
    // have the usual geometrical meanings.
    public CoarseFineBoundary(Domain<3> cfb, Domain<3> acfb, Domain<3> nacfb,
			      BoxedList_Patch patches,
			      int l, Point<1> dir, Point<1> side) {
	coarseBoundary = cfb;
	coarsePatches = patches.toArray();
	RectDomain<3> bbnacfb = nacfb.boundingBox();
	coarseValues = new double[bbnacfb];
	// unsafeBoundary is the part of the boundary which
	// may be covered by other fine patches.
	Domain<3> unsafeBoundary = expand(acfb - cfb, AMR.nRefineP);
	Domain<3> emptyDomain = (Domain<3>) ([0 : -1, 0 : -1, 0 : -1]);

	// Find if any fine patch covers cells of unsafeBoundary.
	coveringFine = emptyDomain;
	foreach (iproc in AMR.processes.domain()) {
	    AMRProcess ap = AMR.processes[iproc];
	    foreach (ii in ap.levels[l].patches.domain()) {
		Patch finePatch = ap.levels[l].patches[ii];
		coveringFine = coveringFine +
		    (Domain<3>)finePatch.domain * unsafeBoundary;
	    }
	}
	if (! coveringFine.isEmpty()) {
	    if (AMR.verbose) {
		System.out.println("cfb:");
		Util.printlnD3(cfb);
		System.out.println("acfb:");
		Util.printlnD3(acfb);
		System.out.println("unsafeBoundary:");
		Util.printlnD3(unsafeBoundary);
		System.out.println("coveringfine:");
		Util.printlnD3(coveringFine);
	    }
	    // LAPDOMAIN extends COVERINGFINE with the
	    // boundary cells needed by the laplacian.
	    Domain<3> lapDomain = coveringFine;
	    foreach (dir2 in [1:3]) {
		foreach (side2 in [-1:1:2]) {
		    lapDomain = lapDomain +
		      (coveringFine + Point<3>.direction(dir2[1]) * side2[1]);
		}
	    }
	    if (AMR.verbose) {
		System.out.println("lapdomain:");
		Util.printlnD3(lapDomain);
	    }
	    // Do another pass over all fine patches to find
	    // the ones we need.
	    BoxedList_Patch fineList = new BoxedList_Patch();
	    Domain<3> coveredLap = emptyDomain;
	    foreach (p in AMR.processes.domain()) {
		AMRProcess ap = AMR.processes[p];
		foreach (ii in ap.levels[l].patches.domain()) {
		    Patch finePatch = ap.levels[l].patches[ii];
		    Domain<3> x = ((Domain<3>)finePatch.domain) * lapDomain;
		    if (! x.isEmpty()) {
			fineList.push(finePatch);
			coveredLap = coveredLap + x;
		    }
		}
	    }
	    fineValues = new double [coveredLap.boundingBox()];
	    finePatches = fineList.toArray();
	    if (AMR.verbose) {
		System.out.println("coveredlap:");
		Util.printlnD3(coveredLap);
	    }

	    // COVEREDLAP is the cells where values are
	    // available for the laplacian computation.  The
	    // actual laplacian can only be evaluated in a
	    // subdomain of this, INSIDELAP.
	    insideLap = coveredLap;
		
	    foreach (dir2 in [1:3]) {
		foreach (side2 in [-1:1:2]) {
		    insideLap = insideLap *
		      (coveredLap + Point<3>.direction(dir2[1]) * side2[1]);
		}
	    }
	    if (AMR.verbose) {
		System.out.println("insidelap:");
		Util.printlnD3(insideLap);
	    }
	    numLaps = new double[coveringFine.boundingBox() / AMR.nRefineP];
	    numLaps.set(0.);
	    foreach (p in coveringFine.boundingBox()) {
		numLaps[p / AMR.nRefineP] += 1.;
	    }
	} else {
	    RectDomain<3> empty = [[0, 0, 0] : [-1, -1, -1]];
	    insideLap = (Domain<3>) empty;
	    fineValues = new double[empty];
	    finePatches = new Patch[[[0] : [-1]]];
	}
    }
}