// Titanium 3-D AMR Poisson solver. // The Patch class. An element of this class contains all the grids // used in the computation of a single rectangular patch. It also // contains structures that specify relationships with the grids in // other patches. // // We use patch-level parallelism. Every processor owns a small // number of patches. class Patch { // The basic stuff. DOMAIN is the domain of this patch. GDOMAIN // is DOMAIN + ghost cells. The grids containing the field values // have domain DOMAIN or GDOMAIN. For instance, PHI has boundary // (ghost) cells, but RHO doesn't. RectDomain<3> domain; RectDomain<3> gdomain; RectDomain<3> coarseDomain; double [3d] local phi; double [3d] local savedPhi; double [3d] local residual; double [3d] local rho; double [3d] local error; double [3d] local laplacian; double [3d] local relaxFactor; int [3d] local tagged; double rhoNorm; // Even though they are associated with this patch, the flux // registers are one level coarser. double [3d] local [1d] local [1d] local fluxRegisters; int [1d] [1d] signatures; double delta; // The boundaries of DOMAIN. They are cell-centered and just // outside DOMAIN. They are indexed by [DIRECTION][SIDE]. // ISPHYSICAL is true when the boundary on that direction and side // is a physical boundary. The problem domain is a box, so a // boundary may not be only partially physical. RectDomain<3> [1d][1d] boundaries = new RectDomain<3>[[1:4]][[-1:2:2]]; boolean [1d][1d] isPhysical = new boolean[[1:4]][[-1:2:2]]; // Red-black domains. REDBLACK[0] = RED, REDBLACK[1] = BLACK, and // RED + BLACK = DOMAIN. Domain<3> [1d] redblack; // Interrelationships. // Patches partially covered by this one. RelatedPatch [1d] coarserPatches; // Patches which partially cover this one. RelatedPatch [1d] finerPatches; // Boundaries with coarser cells, indexed by [DIRECTION][SIDE]. CoarseFineBoundary [1d][1d] coarseBoundaries; // Boundary patches at the same level of this patch. RelatedPatch [1d] adjacentPatches; // Boundaries with finer level. Need these to access the flux // corrections during computation of the laplacian. RelatedPatch [1d] fineBoundaryPatches; private int [3d] local [1d] local scan; // see tagErrors() for explanation private local final /*inline*/ boolean isTagged(Point<3> p) { return tagged[p] > 0; } private static final /*inline*/ int tagsInRect2(int [2d] local f, Point<2> fmin, Point<2> p) { if (p[1] < fmin[1] || p[2] < fmin[2]) return 0; else return f[p]; } /* f is constructed to contain a 2d cumulative histogram, i.e., f[[q, r]] contains the number of tagged points in this patch with coordinate 1 < q and coordinate 2 < r. (The coords can be (y, z) or (x, z) or (x, y) depending on whether dim is 0, 1, or 2.) The last few lines are to compute the number of tagged points with q0 < coord 1 < q1 and r0 < coord2 < r1 given f. Return value: the number of tagged points in this patch in the box and in the plane dim=value. Note that the dim=value plane must intersect the box. */ private final /*inline*/ int tagsInRect(int dim, int value, RectDomain<3> box) { int [2d] local f = scan[dim].slice(dim + 1, value); Point<2> fmin = f.domain.min(), boxmin = box.slice(dim + 1).min(), boxmax = box.slice(dim + 1).max(); if (fmin == boxmin) return f[boxmax]; // speeds up common case else { boxmin -= [1, 1]; return tagsInRect2(f, fmin, boxmax) + tagsInRect2(f, fmin, boxmin) - tagsInRect2(f, fmin, [boxmax[1], boxmin[2]]) - tagsInRect2(f, fmin, [boxmin[1], boxmax[2]]); } } /* Return value, ret, should have ret[dim][value] contain the number of tagged cells in plane dim=value of box intersected with this. */ public local final int [1d] local [1d] tagsInPlanes(RectDomain<3> box) { int [1d] local [1d] ret = new int [[0 : 3]] [1d]; RectDomain<3> N = taggedRect(); box *= N; if (box.isNull()) { RectDomain<1> R = [0 : 0]; ret[0] = new int[R]; ret[1] = new int[R]; ret[2] = new int[R]; } else { for (int dim = 0; dim < 3; dim++) { RectDomain<1> R = [box.min()[dim] : box.max()[dim] + 1]; ret[dim] = new int[R]; ///foreach (i within R) ret[dim][i] = tagsInRect(dim, i, box); } } return ret; } private static final /*inline*/ void min3(int[] x, Point<3> p) { if (p[1] < x[0]) x[0] = p[1]; if (p[2] < x[1]) x[1] = p[2]; if (p[3] < x[2]) x[2] = p[3]; } private static final /*inline*/ void max3(int[] x, Point<3> p) { if (p[1] > x[0]) x[0] = p[1]; if (p[2] > x[1]) x[1] = p[2]; if (p[3] > x[2]) x[2] = p[3]; } private static final void rescan(int [3d] local a, int dim) { RectDomain<3> R = a.domain; foreach (p2 within R.slice(dim + 1)) { int s = 0, q = R.min()[dim], qmax = R.max()[dim]; Point<3> p = (dim == 0) ? [q, p2[1], p2[2]] : ((dim == 1) ? [p2[1], q, p2[2]] : [p2[1], p2[2], q]); Point<3> pinc = Point<3>.direction(dim); for (; q <= qmax; q++, p += pinc) { s += a[p]; a[p] = s; } } } /* tagErrors() finds the cells with high error and sets up data structures to quickly query the spatial structure of the errors later. In particular, we need to answer questions of the form: Given a box b and a plane p, how many tagged cells are in the intersection of this patch, b, and p? */ public local final void tagErrors() { RectDomain<3> R = error.domain; RectDomain<3> N; // subset of R that contains tagged points int [3d] local [1d] local scan = new int [[0 : 3]][3d]; this.scan = scan; /* Define f_{x=i}(j,k) to be the number of tagged cells in the intersection of the x = i plane and the y <= j and z <= k halfspaces. Define f_{y=j}(i,k) and f_{z=k}(i,j) analogously. The goal of the code below is to produce the array scan, such that: scan[0][i,j,k] = f_{x=i}(j,k), scan[1][i,j,k] = f_{y=j}(i,k), and scan[2][i,j,k] = f_{z=k}(i,j). */ // X scan then Y, and compute N { int[] lo = new int[3], hi = new int[3]; scan[2] = new int[R]; boolean seenAny = false; foreach (yz within R.slice(1)) { for (int s = 0, x = R.min()[1], xmax = R.max()[1]; x <= xmax; x++) { Point<3> p = [x, yz[1], yz[2]]; if (isTagged(p)) { s++; if (seenAny) { min3(lo, p); max3(hi, p); } else { seenAny = true; lo[0] = p[1]; lo[1] = p[2]; lo[2] = p[3]; hi[0] = p[1]; hi[1] = p[2]; hi[2] = p[3]; } } scan[2][p] = s; } } N = [[lo[0], lo[1], lo[2]] : [hi[0], hi[1], hi[2]] + [1, 1, 1]]; scan[2] = scan[2].restrict(N); rescan(scan[2], 1); } // Y scan then Z scan[0] = new int[N]; foreach (xz within N.slice(2)) { for (int s = 0, y = R.min()[1], ymax = R.max()[1]; y <= ymax; y++) { Point<3> p = [xz[1], y, xz[2]]; if (isTagged(p)) s++; scan[0][p] = s; } } rescan(scan[0], 2); // Z scan then X scan[1] = new int[N]; foreach (xy within N.slice(3)) { for (int s = 0, z = R.min()[2], zmax = R.max()[2]; z <= zmax; z++) { Point<3> p = [xy[1], xy[2], z]; if (isTagged(p)) s++; scan[1][p] = s; } } rescan(scan[1], 0); } /* Returns a rectangle that just covers the tagged cells in the patch. */ public local final RectDomain<3> taggedRect() { return scan[0].domain; } // Allocate and initialize grids for a patch. DOMAIN and GDOMAIN // must be set before calling this. void allocateGrids() { phi = new double[gdomain]; savedPhi = new double[gdomain]; error = new double[gdomain]; residual = new double[domain]; rho = new double[domain]; laplacian = new double[domain]; relaxFactor = new double[domain]; RectDomain<3> sdomain = domain.shrink(1); foreach (p within sdomain) { relaxFactor[p] = 1./4.; } foreach (p within domain - sdomain) { relaxFactor[p] = 3./16.; } } void coarseToFine(boolean doPhi, boolean reflux) { double [3d] fine = doPhi ? phi : error; for (int direction = 1; direction <= 3; direction++) { for (int side = -1; side <= 1; side += 2) { // Useful numbers and vectors. // Normal directions. Point<3> normDir = Point<3>.direction(direction); Point<3> orient = normDir * side; Point<3> offset = (normDir - orient) / 2; // Tangential directions. int k1 = (direction + 1) % 3 + 1; int k2 = (direction + 2) % 3 + 1; Point<3> t1 = Point<3>.direction(k1); Point<3> t2 = Point<3>.direction(k2); Point<3> tt = t1 + t2; CoarseFineBoundary cb = coarseBoundaries[direction][side]; updateCoarseBoundary(cb, doPhi); // Interpolate the ghost cell values. double [3d] cv = cb.coarseValues; foreach (p within cb.coarseBoundary) { double [3d] fr = fluxRegisters[direction][side]; if (reflux) { // Initialize flux register. Sign convention // and centering is done in a way that the // correction is done by adding the contents // of the flux register to the value of the // operator at the coarse cell centered in the // same location: foreach (direction, side) // foreach (p in coarseBoundary) // laplacian[kCoarse][p] += // fluxRegister[direction][side][p]. fr[p] = (cv[p] - cv[p - orient]) / (AMR.nRefine * delta); } double x1 = (cv[p + t1] - cv[p - t1]) / 2; double x2 = (cv[p + t2] - cv[p - t2]) / 2; double x1x1 = (cv[p + t1] + cv[p - t1] - 2 * cv[p]) / 2; double x2x2 = (cv[p + t2] + cv[p - t2] - 2 * cv[p]) / 2; double x1x2 = (cv[p + t1 + t2] + cv[p - t1 - t2] - cv[p - t1 + t2] + cv[p + t1 - t2]) / 4; // Compute interpolant. Point<3> pf0 = (p + offset) * AMR.nRefine - offset; foreach (q within AMR.refinedCell) { double xx1 = (q[k1] + .5) / AMR.nRefine - .5; double xx2 = (q[k2] + .5) / AMR.nRefine - .5; fine[pf0 + q] = cv[p] + xx1 * x1 + xx2 * x2 + xx1 * xx1 * x1x1 + xx2 * xx2 * x2x2 + xx1 * xx2 * x1x2; fine[pf0 + q] = AMR.c1 * fine[pf0 + q] + AMR.c2 * fine[pf0 + q - orient] + AMR.c3 * fine[pf0 + q - [2, 2, 2] * orient]; if (reflux) { // Increment values in flux register with // properly normalized fine fluxes. fr[p] -= (cv[pf0 + q] - fine[pf0 + q - orient]) / delta / Math.pow(AMR.nRefine, 2); } } } } } } void updateCoarseBoundary(CoarseFineBoundary cb, boolean doPhi) { // Update with coarse values where available. cb.coarseValues.set(0); foreach (i within cb.coarsePatches.domain) { cb.coarseValues.copy(doPhi ? cb.coarsePatches[i].phi : cb.coarsePatches[i].error); } // Collect fine values. foreach (i within cb.finePatches.domain) { cb.fineValues.copy(doPhi ? cb.finePatches[i].phi : cb.finePatches[i].error); } // Average and correct where covered by fine patches. foreach (p within cb.fineValues.domain) { cb.coarseValues[p / AMR.nRefine] += cb.fineValues[p] / Math.pow(AMR.nRefine, 3); } foreach (p within cb.insideLap) { double L = -6 * cb.fineValues[p]; for (int k = 1; k <= 3; k++) /*unroll*/ { L += cb.fineValues[p + Point<3>.direction(k)] + cb.fineValues[p - Point<3>.direction(k)]; } cb.coarseValues[p / AMR.nRefine] += L / 2 / cb.numLaps[p]; } } // Before calling this, make sure that the flux registers at L+1 // are up to date, by calling coarseToFine(true, true) over all // L+1 grids. void computePhiLaplacian(boolean reflux) { // Update phi boundaries from adjacent grids. foreach (p within adjacentPatches.domain) { RelatedPatch rp = adjacentPatches[p]; phi.copy(rp.patch.phi); } // Set ghost cells on boundaries with coarser grids, // unless this is the coarsest level. if (! coarserPatches.domain.isNull()) { coarseToFine(true, false); } computeLaplacian(true); if (reflux) { foreach (i within fineBoundaryPatches.domain) { RelatedPatch rp = fineBoundaryPatches[i]; for (int direction = 1; direction <= 3; direction++) { for (int side = -1; side <= 1; side += 2) { double [3d] fr = rp.patch.fluxRegisters[direction][side]; // This domain intersection is one of the few // things I don't bother to precompute. foreach (p within laplacian.domain * fr.domain) { laplacian[p] += fr[p]; } } } } } } // Set ghost cells to impose physical BCs. void physicalBC(boolean doPhi) { int direction, side; double [3d] field = doPhi ? phi : error; for (direction = 1; direction <= 3; direction++) { for (side = -1; side <= 1; side += 2) { if (isPhysical[direction][side]) { foreach (p within boundaries[direction][side]) { Point<3> v = Point<3>.direction(direction); field[p + v] = field[p - v] / 3 - 2 * field[p]; } } } } } void computeLaplacian(boolean doPhi) { double [3d] field = doPhi ? phi : error; foreach (p within laplacian.domain) { laplacian[p] = ( field[p + [ 0, 0, 1]] + field[p + [ 0, 0,-1]] + field[p + [ 0, 1, 0]] + field[p + [ 0,-1, 0]] + field[p + [ 1, 0, 0]] + field[p + [-1, 0, 0]] - 6 * field[p]) / Math.pow(delta, 2); } } void setResidualErrorLapDifference(boolean doLap) { double [3d] dest = doLap ? laplacian : residual; physicalBC(false); coarseToFine(false, false); computeLaplacian(false); foreach (p within dest.domain) { dest[p] = residual[p] - laplacian[p]; } } // Compute the difference between the laplacian of the coarsened // solution, and the coarsened right-hand side. Compare with norm // of fine rho, and tag points where too large. boolean estimateError() { double tolerance = rhoNorm * 0.01; boolean anyTagged = false; foreach (p within domain / AMR.nRefineP) { double error = 0; foreach (q within AMR.refinedCell) /*unroll*/ { double coarseLaplacian = phi[p * AMR.nRefineP + q + [ 0, 0, 1]] + phi[p * AMR.nRefineP + q + [ 0, 1, 0]] + phi[p * AMR.nRefineP + q + [ 1, 0, 0]] + phi[p * AMR.nRefineP + q + [ 0, 0, -1]] + phi[p * AMR.nRefineP + q + [ 0, -1, 0]] + phi[p * AMR.nRefineP + q + [-1, 0, 0]] - phi[p * AMR.nRefineP + q] * 6; error += coarseLaplacian - rho[p * AMR.nRefineP + q]; } error *= 1 / (Math.pow(delta, 2) * Math.pow(AMR.nRefine, 4)); tagged[p] = Math.abs(error) > tolerance ? 1 : 0; anyTagged = anyTagged || (tagged[p] == 1); } return anyTagged; } // Go through all other patches at level L (this level), // above, and below, and compute relationships. void computeRelationships(int l) { BoxedList_RelatedPatch adjacentList = new BoxedList_RelatedPatch(); BoxedList_RelatedPatch coarserList = new BoxedList_RelatedPatch(); BoxedList_RelatedPatch finerList = new BoxedList_RelatedPatch(); foreach (p within AMR.processes.domain) { AMRProcess ap = AMR.processes[p]; foreach (ii within ap.levels[l].patches.domain) { Patch otherPatch = ap.levels[l].patches[ii]; if (this != otherPatch) { RectDomain<3> b = gdomain * otherPatch.domain; if (! b.isNull()) { RelatedPatch rp = new RelatedPatch(); rp.domain = b; rp.patch = otherPatch; adjacentList.push(rp); } } } if (l > 0) { foreach (ii within ap.levels[l-1].patches.domain) { Patch otherPatch = ap.levels[l-1].patches[ii]; RectDomain<3> x = otherPatch.domain * coarseDomain; if (! x.isNull()) { RelatedPatch rp = new RelatedPatch(); rp.patch = otherPatch; rp.domain = x; coarserList.push(rp); } } // Compute coarse-fine boundaries. Must // find all coarse patches that contribute // to the boundary computation. foreach (dir within [1:4]) { foreach (side within [-1:2:2]) { BoxedList_Patch patches = new BoxedList_Patch(); RectDomain<3> cb = ((boundaries[dir][side] / AMR.nRefineP) - Point<3>.direction(dir[1], side[1])); // Accrete cb tangentially. Point<3> tangent = [1, 1, 1] - Point<3>.direction(dir[1]); RectDomain<3> acb = [cb.lwb() - tangent : cb.upb() + tangent]; foreach (ii within ap.levels[l-1].patches.domain) { Patch otherPatch = ap.levels[l-1].patches[ii]; if (! (acb * otherPatch.domain).isNull()) { patches.push(otherPatch); } } coarseBoundaries[dir][side] = patches.length() == 0 ? null : new CoarseFineBoundary(cb, acb, patches, l, dir, side); } } } if (l < AMR.nLevels - 1) { foreach (ii within ap.levels[l+1].patches.domain) { Patch otherPatch = ap.levels[l+1].patches[ii]; RectDomain<3> x = this.domain * (otherPatch.domain / AMR.nRefineP); if (! x.isNull()) { RelatedPatch rp = new RelatedPatch(); rp.patch = otherPatch; rp.domain = x; finerList.push(rp); } } } } this.adjacentPatches = adjacentList.toArray(); this.coarserPatches = coarserList.toArray(); this.finerPatches = finerList.toArray(); } } // Auxiliary classes. class RelatedPatch { Patch patch; RectDomain<3> local domain; // The coarser domain. } // Coarse-fine boundary patches require additional information // used by the quadratic interpolation in the coarse-fine BC // calculation. class CoarseFineBoundary { RectDomain<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; // OK to compute the laplacian of FINEVALUES in INSIDELAP. Domain<3> insideLap; // NUMLAPS is the number of fine laplacians available // per coarse cell. int [3d] numLaps; // CB is the coarse side of the boundary. ACB is the // tangentially accreted CB. PATCHES is a list of // coarse patches that intersect ACB. L is the index of // this level. DIR and SIDE have the usual geometrical // meanings. public CoarseFineBoundary(RectDomain<3> cb, RectDomain<3> acb, BoxedList_Patch patches, int l, Point<1> dir, Point<1> side) { coarseBoundary = cb; coarsePatches = patches.toArray(); // Accrete normally. RectDomain<3> nacb = acb.accrete(1, - dir[1] * side[1]); coarseValues = new double[acb]; RectDomain<3> fineBoundary = [cb.lwb() * AMR.nRefineP : cb.upb() * AMR.nRefineP]; // Find fine patches that cover cells of CB (if // any). Domain<3> coveringFine = [0:0, 0:0, 0:0]; foreach (ii within AMR.levels[l].patches.domain) { Patch finePatch = AMR.levels[l].patches[ii]; coveringFine += finePatch.domain * fineBoundary; } if (! coveringFine.isNull()) { // LAPDOMAIN extends COOVERINGFINE with the // boundary cells needed by the laplacian. Domain<3> lapDomain = coveringFine; foreach (dir2 within [1:4]) { foreach (side2 within [-1:2:2]) { lapDomain += coveringFine + Point<3>.direction(dir2[1]) * side2[1]; } } // Do another pass over all fine patches to find // the ones we need. BoxedList_Patch fineList = new BoxedList_Patch(); Domain<3> coveredLap = [0:0, 0:0, 0:0]; foreach (ii within AMR.levels[l].patches.domain) { Patch finePatch = AMR.levels[l].patches[ii]; Domain<3> x = finePatch.domain * lapDomain; if (! x.isNull()) { fineList.push(finePatch); coveredLap += x; } } fineValues = new double [coveredLap.boundingBox()]; finePatches = fineList.toArray(); // 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 within [1:4]) { foreach (side2 within [-1:2:2]) { insideLap *= insideLap + Point<3>.direction(dir2[1]) * side2[1]; } } numLaps = new int[coveringFine.boundingBox() / AMR.nRefineP]; numLaps.set(0); } } } // Per-process data structure. class AMRProcess { Level local [1d] local levels; double localResidualNorm; boolean anyTagged; // Iterate on a fixed grid hierarchy. // // Sometimes UP and DOWN stand for, respectively, FINE and COARSE. // Finest is top, coarsest is bottom. It's better to use FINE and // COARSE explicitly. // // 0 is the coarsest level, 1 is finer, etc. single boolean single mainLoop(int single maxIterations, boolean checkConvergence) { int single iterations; for (iterations = 0; iterations < maxIterations; iterations++) { relaxFinest(); relaxFineToCoarse(); relaxCoarsest(); relaxCoarseToFine(); // Check for convergence. boolean single converged; for (int single l = 0; l < AMR.nLevels; l++) { localResidualNorm = 0; foreach (i within levels[l].patches.domain) { localResidualNorm += AMR.norm(levels[l].patches[i].residual, levels[l].delta); } // There are other ways of doing this that may be // faster, but this is the most concise. // double [1d] norms = new double[Process.myTeam().domain]; // Process.exchange(localResidualNorm, norms); // double residualNorm = 0; // foreach (i within norms.domain) { // residualNorm += norms[i]; // } boolean cv = true; if (Ti.thisProc() == 0) { double residualNorm = 0; foreach (i within AMR.processes.domain) { residualNorm += AMR.processes[i].localResidualNorm; } if (residualNorm > AMR.epsilon * AMR.rhoNorms[l]) { cv = false; break; } } converged = broadcast cv from 0; } if (converged) break; } return iterations < maxIterations; } void relaxFinest() { Level local l = levels[AMR.finest]; foreach (i within l.patches.domain) { Patch local patch = l.patches[i]; patch.computePhiLaplacian(false); foreach (p within patch.residual.domain) { patch.residual[p] = patch.rho[p] - patch.laplacian[p]; patch.error[p] = 0; } } } single void relaxFineToCoarse() { for (int single l = AMR.finest; l > 0; l--) { // Clear correction at l-1. foreach (i within levels[l-1].patches.domain) { Patch patch = levels[l-1].patches[i]; patch.error.set(0); } foreach (i within levels[l].patches.domain) { Patch patch = levels[l].patches[i]; // Relax at this level. AMR.gsrbRelax(patch.error, patch.residual, patch.relaxFactor, patch.boundaries, patch.delta); // Save solution. patch.savedPhi.copy(patch.phi); // Apply correction. foreach (p within patch.error.domain) { patch.phi[p] += patch.error[p]; } // Update flux registers for forthcoming computation // of L(phi) one level down. patch.coarseToFine(true, true); } Ti.barrier(); foreach (i within levels[l-1].patches.domain) { Patch patch = levels[l-1].patches[i]; patch.computePhiLaplacian(true); // Set the residual at l-1, over all the domain, then // later set again the part under the fine region. // (The fine region is small, so we don't bother to // set only the exposed part.) foreach (p within patch.domain) { patch.residual[p] = patch.rho[p] - patch.laplacian[p]; } } Ti.barrier(); foreach (i within levels[l].patches.domain) { Patch patch = levels[l].patches[i]; // Compute residual - LNF(error) and temporarily store // in laplacian. patch.setResidualErrorLapDifference(true); // Set the residual at l-1 under the fine region, by // simple average. foreach (j within patch.coarserPatches.domain) { Patch rp = patch.coarserPatches[j].patch; Domain<3> rd = patch.coarserPatches[j].domain; foreach (p within rd) { rp.residual[p] = 0; foreach (q within AMR.refinedCell) /*unroll*/ { rp.residual[p] += patch.laplacian[p * AMR.refineVector + q]; } rp.residual[p] /= 8; } } } Ti.barrier(); } } void relaxCoarsest() { foreach (i within levels[0].patches.domain) { Patch patch = levels[0].patches[i]; // Clear correction. patch.error.set(0); // Relax/solve at level 0. Multigrid.relax(patch.error, patch.rho, patch.residual, patch.relaxFactor, patch.delta, 0); // Apply correction. foreach (p within patch.phi.domain) { patch.phi[p] += patch.error[p]; } } } single void relaxCoarseToFine() { for (int single l = 1; l < AMR.nLevels; l++) { foreach (j within levels[l].patches.domain) { Patch patch = levels[l].patches[j]; // Correct the correction by interpolation from // coarser levels. foreach (i within patch.coarserPatches.domain) { RelatedPatch rp = patch.coarserPatches[i]; foreach (p within rp.patch.domain * patch.coarseDomain) { foreach (q within AMR.refinedCell) /*unroll*/ { patch.error[p * AMR.nRefineP + q] += rp.patch.error[p]; } } } // Decrement the residual by the laplacian of the // error. patch.setResidualErrorLapDifference(false); // Compute correction to apply to error. Use phi for // temporary space. patch.phi.set(0); AMR.gsrbRelax(patch.phi, patch.residual, patch.relaxFactor, patch.boundaries, patch.delta); foreach (p within patch.error.domain) { // Correct correction. patch.error[p] += patch.phi[p]; // Compute new value of solution. patch.phi[p] = patch.savedPhi[p] + patch.error[p]; } } Ti.barrier(); } } // Create new grids for hierarchy. // // Construct new grid, from finest to coarsest level. A new level // may be added here, or several levels may disappear. Need one // extra level above level HIGHEST-BAD-LEVEL. The regridding goes // from that extra level down to 1 (don't regrid level 0, that's // the problem domain). single void regridAsNeeded() { int single highestBadLevel; int single newNLevels; int single newMaxLevel; for (int single l = AMR.nLevels - 1; l >= 0; l--) { // Compute tagged cells on each patch. Level level = levels[l]; anyTagged = false; foreach (i within level.patches.domain) { anyTagged = anyTagged || level.patches[i].estimateError(); } int hbl = -1; if (Ti.thisProc() == 0) { foreach (i within AMR.processes.domain) { if (AMR.processes[i].anyTagged) { hbl = l; break; } } } highestBadLevel = broadcast hbl from 0; if (highestBadLevel > -1) break; } if (highestBadLevel == -1) return; RectDomain<3> [1d] [1d] boxes; newMaxLevel = highestBadLevel + 1; newNLevels = newMaxLevel + 1; // Need to check for too many levels (l >= maxLevels). for (int single l = newMaxLevel; l > 0; l--) levels[l].boxes = Refiner.refinement(levels[l-1], l > 1 ? levels[l-2] : null); allocatePatches(newNLevels); } single void allocatePatches(int single newNLevels) { // No need to free the old patches. AMR.nLevels = newNLevels; AMR.finest = newNLevels - 1; for (int l = 0; l < newNLevels; l++) { levels[l].patches = new Patch[levels[l].boxes.domain]; foreach (i within levels[l].patches.domain) { Patch patch = levels[l].patches[i]; patch.domain = levels[l].boxes[i]; patch.coarseDomain = patch.domain / AMR.nRefineP; patch.gdomain = patch.domain.accrete(1); for (int direction = 1; direction <= 3; direction++) { for (int side = -1; side <= 1; side +=2) { RectDomain<3> b = (RectDomain<3>) (patch.domain - (patch.domain + Point<3>.direction(-direction))); patch.boundaries[direction][side] = b; RectDomain<3> physical = levels[0].patches[0].boundaries[direction][side]; patch.isPhysical[direction][side] = ! ((b / Point<3>.all((int) Math.pow(AMR.nRefine, l))) * physical).isNull(); } } } } Ti.barrier(); for (int l = 0; l < newNLevels; l++) { levels[l].patches = new Patch[levels[l].boxes.domain]; foreach (i within levels[l].patches.domain) { Patch patch = levels[l].patches[i]; // Allocate various grids. patch.allocateGrids(); // Compute relationships. patch.computeRelationships(l); } } } } // The Level class. All information pertaining to one level. class Level { double delta; RectDomain<3> [1d] boxes; Patch [1d] patches; } // The Multigrid class. Used at the coarsest level. All methods are // sequential. class Multigrid { final int levels = 10; static double [3d][1d] rhos; static double [3d][1d] errors; static double [3d][1d] residuals; static double [3d][1d] factors; static RectDomain<3> [1d][1d][1d] boundaries; static void relax(double [3d] field, double [3d] rhs, double [3d] residual, double [3d] factor, double h, int level) { // This is backward from AMR. Coarser levels have higher // numbers. if (level < levels) { int coarse = level + 1; double [3d] coarseField = errors[coarse]; double [3d] coarseRho = rhos[coarse]; RectDomain<3> coarseInterior = coarseRho.domain; // Relax on this level. AMR.gsrbRelax(field, rhs, factor, Multigrid.boundaries[level], h); // Initialize solution at next level. coarseField.set(0); // Compute fine residual. First, set dirichlet boundary // conditions like in Patch.physicalBC(). for (int direction = 1; direction <= 3; direction++) { for (int side = -1; side <= 1; side += 2) { foreach (p within boundaries[level][direction][side]) { Point<3> v = Point<3>.direction(direction); field[p + v] = field[p - v] / 3 - 2 * field[p]; } } } foreach (p within residual.domain) { residual[p] = rhs[p] - ( field[p + [ 0, 0, 1]] + field[p + [ 0, 0, -1]] + field[p + [ 0, 1, 0]] + field[p + [ 0, -1, 0]] + field[p + [ 1, 0, 0]] + field[p + [-1, 0, 0]] - 6 * field[p]) / (h * h); } // Average fine residual onto coarse rhs. foreach (p within coarseField.domain) { coarseField[p]= 0; foreach (q within AMR.refinedCell) /*unroll*/ { coarseField[p] += residual[p * AMR.nRefineP + q]; } coarseField[p] /= 8; } // Recursive call. relax(coarseField, coarseRho, residuals[coarse],factors[coarse], h * 2, coarse); // Correct fine solution by interpolating coarse // correction. foreach (p within coarseInterior) { foreach (q within AMR.refinedCell) /*unroll*/ { field[p * AMR.nRefineP + q] += coarseField[p]; } } // Relax again. AMR.gsrbRelax(field, rhs, factor, Multigrid.boundaries[level], h); } else { // Coarsest level. relax_1x1(field, rhs, h); } } static void relax_1x1(double [3d] field, double [3d] rhs, double h) { field[[0, 0, 0]] += 0.5 * Math.pow(h, 2) * rhs[[0, 0, 0]]; } } // The AMR class, representing the full hierarchy of grids and all // else. class AMR { static AMRProcess [1d] single processes; static int single maxLevels; static int single nLevels; static int single finest; // nLevels - 1 static final int maxRegriddings = 10; static final int coarsest = 0; static final int nRefine = 2; static final Point<3> nRefineP = [2, 2, 2]; static final double epsilon = 0.01; static final Point<3> refineVector = [nRefine, nRefine, nRefine]; static final Domain<3> refinedCell = [[0, 0, 0] : [nRefine, nRefine, nRefine]]; static double [1d] rhoNorms = new double [[0 : maxLevels]]; static final double c1 = 8. / (nRefine + 1) * (nRefine + 3); static final double c2 = (nRefine - 1) * 2. * (nRefine + 3) / (nRefine + 1) * (nRefine + 3); static final double c3 = - (nRefine + 1.) * (nRefine - 1) / (nRefine + 1) * (nRefine + 3); static Multigrid mg = new Multigrid(); // A global list of all patches at all levels. There is // a per-process list as well, but sometimes it's more // convenient to have this one. static Level [1d] levels; // Compute the square of the norm. static double norm(double [3d] field, double h) { double n = 0; foreach (p within field.domain) { n += Math.pow(field[p] * h, 2); } return n; } static void gsrbRelax(double [3d] field, double [3d] rhs, double [3d] factor, RectDomain<3> [1d][1d] boundaries, double h) { // Set boundary conditions. for (int direction = 1; direction <= 3; direction++) { for (int side = -1; side <= 1; side += 2) { foreach (p within boundaries[direction][side]) { Point<3> p1 = Point<3>.direction(direction, side); field[p - p1] = field[p + p1] / 3 - field[p] * 2; } } } // Relax. for (int i = 0; i < 2; i++) { foreach (p within (rhs.domain / [2, 2, 2]) + [i, i, i]) { field[p] += factor[p] * ( field[p + [ 0, 0, 1]] + field[p + [ 0, 0,-1]] + field[p + [ 0, 1, 0]] + field[p + [ 0,-1, 0]] + field[p + [ 1, 0, 0]] + field[p + [-1, 0, 0]] - 6 * field[p] - Math.pow(h, 2) * rhs[p]); } } } single static void main() { // Initialization of single AMR static variables goes here. AMR.maxLevels = 20; AMRProcess single p = new AMRProcess(); // Make every process know about every other process. processes = new AMRProcess[[0 : Ti.numProcs() - 1]]; processes.exchange(p); if (Ti.thisProc() == 0) { // Sequential initialization code. // Construct initial hierarchy and initialize right-hand // side. // // The problem domain, in real-world coordinates, starts // at (0,0) and ends at (1,1). N can change, but the // problem domain doesn't. } Ti.barrier(); int single regriddings = 0; for (;;) { if (p.mainLoop(10, true)) break; regriddings++; if (regriddings >= maxRegriddings) break; p.regridAsNeeded(); } } } //_____________________________________________________________________________ class TagProfileContribution { public BoxedList_int_1d__1d_local local Q; public local /*final*/ TagProfileContribution() { Q = new BoxedList_int_1d__1d_local(); } public local /*final*/ void add(Patch local patch, RectDomain<3> box) { Q.push(patch.tagsInPlanes(box)); } } // end of class TagProfileContribution class TagProfile { Level level; RectDomain<3> box; // a[0][k] tells how many tagged in x = k plane // a[1][k] tells how many tagged in y = k plane // a[2][k] tells how many tagged in z = k plane int [1d] local [1d] local a; TagProfileContribution [1d] local tpc; TagProfileContribution local myContrib; single private final void collectProfiles() { Ti.barrier(); if (Ti.thisProc() == 0) { a = new int [[0 : 3]][1d]; a[0] = new int[box.slice(2).slice(2)]; a[1] = new int[box.slice(1).slice(2)]; a[2] = new int[box.slice(1).slice(1)]; foreach (i within tpc.domain) { BoxedList_int_1d__1d_local Q = tpc[i].Q; while (!Q.isEmpty()) { int [1d] local [1d] b = Q.pop(); for (int dim = 0; dim < 3; dim++) { int [1d] local c = new int[b[dim].domain]; c.copy(b[dim]); ///if (dim == 0) foreach (p within c.domain) patchTaggedTotal += c[p]; foreach (p within c.domain) a[dim][p] += c[p]; } } } } } single public TagProfile(Level level, RectDomain<3> box) { this.level = level; this.box = box; if (Ti.thisProc() == 0) tpc = new TagProfileContribution[Ti.myTeam().domain]; myContrib = new TagProfileContribution(); foreach (i within level.patches.domain) myContrib.add(level.patches[i], box); // Largest left-hand side in the program! (broadcast tpc from 0)[Ti.thisProc()] = myContrib; collectProfiles(); } // Return the list of all maximal planar cuts of thickness >= 1 that don't // touch any tagged cells. By maximal I mean the cut cannot be made bigger // without touching a tagged cell or going outside the box. public BoxedList_Cut blanks() { BoxedList_Cut output = new BoxedList_Cut(); for (int d = 0; d < 3; d++) for (int v = a[d].domain.min()[1], vmax = a[d].domain.max()[1]; v <= vmax; v++) if (a[d][v] == 0) { int thickness = 1, v0 = v; while (v < vmax && a[d][++v] == 0) thickness++; output.push(new Cut(d, v0, thickness)); } return output; } // Return the list of all planar cuts that are at points of high // derivative in a. This is intended to cut (2D example) // xxxxxxxxx // xxxxxxxxx // xxxxx // xxxxx // ^ // // at the spot marked with a caret (and/or between the top half and bottom). public BoxedList_Cut highDeriv() { return new BoxedList_Cut(); /* for (int d = 0; d < 3; d++) for (int v = a[d].min()[0], vmax = a[d].max()[0]; v <= vmax; v++) if (highDerivAt(d, v)) output.push(new Cut(d, v, 0)); return output; */ } // Return (the singleton list of) the cut that cuts // box's longest dimension in half. public BoxedList_Cut cutInHalf() { Point<3> boxgirth = box.max() - box.min(); int biggestDim = 0; int biggest = boxgirth[1]; for (int d = 1; d < 3; d++) if (biggest < boxgirth[d]) { biggestDim = d; biggest = boxgirth[d]; } return new BoxedList_Cut(new Cut(biggestDim, box.min()[biggestDim] + biggest / 2, 0)); } } // end of class TagProfile class Cut { int dim, value, thickness; /* sample: dim = 1, value = 3, thickness = k applied to [[0, 0, 0] : [N, N, N]] yields [[0, 0, 0] : [N, 3, N]] and [[0, 3 + k, 0] : [N, N, N]] */ public Cut(int dim, int value, int thickness) { this.dim = dim; this.value = value; this.thickness = thickness; } // Apply the given cut to the given box, and cons the result(s) onto out. // If the given cut doesn't pass through the given box, should be // equivalent to out.push(box). private static final void apply1(Cut cut, RectDomain<3> box, BoxedList_RectDomain_3 out) { int loU[] = new int[3]; int hiL[] = new int[3]; Point<3> boxU = box.max() + [1, 1, 1], boxL = box.min(); RectDomain<3> c; for (int i = 0; i < 3; i++) if (i == cut.dim) { loU[i] = cut.value; hiL[i] = cut.value + cut.thickness; } else { loU[i] = boxU[i]; hiL[i] = boxL[i]; } c = [boxL : [loU[0], loU[1], loU[2]]] * box; if (!c.isNull()) out.push(c); c = [[hiL[0], hiL[1], hiL[2]] : boxU] * box; if (!c.isNull()) out.push(c); } // Apply the given cuts to the given box, resulting in a list of boxes. public static final BoxedList_RectDomain_3 apply(BoxedList_Cut cuts, RectDomain<3> box) { BoxedList_RectDomain_3 boxes = new BoxedList_RectDomain_3(box); while (!cuts.isEmpty()) { Cut cut = cuts.pop(); BoxedList_RectDomain_3 newboxes = new BoxedList_RectDomain_3(); while (!boxes.isEmpty()) apply1(cut, boxes.pop(), newboxes); boxes = newboxes; } return boxes; } } // end of class Cut class Refiner { static final RectDomain<3> emptyRect = [[0, 0, 0] : [0, 0, 0]]; /* Figure out how to chop up the given box into pieces, add each piece either to finished, if it is satisfactory, or to toDo. */ static single final void chopBox(RectDomain<3> box, BoxedList_RectDomain_3 finished, BoxedList_RectDomain_3 toDo, Level level) { TagProfile tp = new TagProfile(level, box); BoxedList_Cut cuts = tp.blanks(); if (cuts.isEmpty()) cuts = tp.highDeriv(); if (cuts.isEmpty()) cuts = tp.cutInHalf(); { BoxedList_RectDomain_3 newBoxes = Cut.apply(cuts, box); while (!newBoxes.isEmpty()) { RectDomain<3> newBox = newBoxes.pop(); ///if (level.satisfactoryPatch(newBox)) finished.push(newBox); ///else toDo.push(newBox); } } } /* Given a list, l, of RectDomains, return a different 1d array on each processor such that the concatenation of all the arrays returned contains exactly what l contains. */ static final RectDomain<3> [1d] local distributeRects(BoxedList_RectDomain_3 l) { int mylength = 0, length = l.length(); while (length > Ti.thisProc()) { mylength++; length -= Ti.numProcs(); } { RectDomain<3> [1d] local ret = new RectDomain<3>[[0 : mylength]]; while (mylength > 0) if (length-- % Ti.numProcs() == Ti.thisProc()) ret[--mylength] = l.pop(); else l.pop(); return ret; } } /* The idea for refinement: RectDomain<3> [1d] local refinement(Level level, Level coarserLevel) { in parallel: compute scans in each direction on each patch; boxes = new Set(bounding box of union of all patches); // outputRects = new Set(); // need not be shared while (boxes is not empty) { take a box from boxes; in parallel: look for gaps, and points of high derivative; sequentially: decide where to chop; choppedBoxes = chop up the box; for all b in choppedBoxes { if satisfactory(b) outputRects.adjoin(b); else boxes.adjoin(b); } } in parallel: for all patches p in outputRects { if there is a patch in coarserLevel that needs to grow, make it grow; } return outputRects; } To improve parallelism, one could do multiple iterations of the while loop in parallel. I think an asynchronous model (task queue/RPC) would be even better. */ single static final RectDomain<3> [1d] local refinement(Level level, Level coarserLevel) { BoxedList_RectDomain_3 outputRects = new BoxedList_RectDomain_3(); BoxedList_RectDomain_3 boxes = new BoxedList_RectDomain_3(); RectDomain<3> localBBox, globalBBox; foreach (i within level.patches.domain) { Patch patch = level.patches[i]; patch.tagErrors(); localBBox = boundingBox(localBBox, patch.taggedRect()); } globalBBox = allReduceBBox(localBBox); if (Ti.thisProc() == 0) boxes.push(globalBBox); while (true) { RectDomain<3> single box; RectDomain<3> b; if (Ti.thisProc() == 0) b = boxes.isEmpty() ? emptyRect : boxes.pop(); /// else assert(boxes.isEmpty()); box = broadcast b from 0; if (box.isNull()) break; else chopBox(box, outputRects, boxes, level); } // if (coarserLevel != null) // ensureProperNesting(outputRects, coarserLevel); return distributeRects(outputRects); } // Return the smallest unit stride RectDomain that contains both a and b. static final RectDomain<3> boundingBox(RectDomain<3> a, RectDomain<3> b) { if (a.isNull()) return b; if (b.isNull()) return a; return [ [a.lwb()[1] < b.lwb()[1] ? a.lwb()[1] : b.lwb()[1], a.lwb()[2] < b.lwb()[2] ? a.lwb()[2] : b.lwb()[2], a.lwb()[3] < b.lwb()[3] ? a.lwb()[3] : b.lwb()[3]] : [a.upb()[1] > b.upb()[1] ? a.upb()[1] : b.upb()[1], a.upb()[2] > b.upb()[2] ? a.upb()[2] : b.upb()[2], a.upb()[3] > b.upb()[3] ? a.upb()[3] : b.upb()[3]] ]; } static final RectDomain<3> allReduceBBox(RectDomain<3> r) { Point<3> lo = r.min(), hi = r.max() + [1, 1, 1]; // return [[allReduceMin(lo[0]), allReduceMin(lo[1]), allReduceMin(lo[2])] : // [allReduceMax(hi[0]), allReduceMax(hi[1]), allReduceMax(hi[2])]]; return r; } } // end of class Refiner class List_RelatedPatch { public RelatedPatch car; public List_RelatedPatch cdr; // Main constructor. public List_RelatedPatch(RelatedPatch x, List_RelatedPatch rest) { car = x; cdr = rest; } public List_RelatedPatch cons(RelatedPatch x) { return new List_RelatedPatch(x, this); } public /*inline*/ RelatedPatch first() { return car; } public /*inline*/ List_RelatedPatch rest() { return cdr; } public int length() { int i = 0; List_RelatedPatch r = this; while (r != null) { i++; r = r.rest(); } return i; } } // The Boxed List class. Used to construct a list // element-by-element and convert it into an array. class BoxedList_RelatedPatch { List_RelatedPatch l = null; public BoxedList_RelatedPatch() { } public BoxedList_RelatedPatch(RelatedPatch x) { push(x); } public void push(RelatedPatch item) { l = new List_RelatedPatch(item, l); } public RelatedPatch pop() { RelatedPatch t = l.first(); l = l.rest(); return t; } public RelatedPatch [1d] toArray() { RelatedPatch [1d] a = new RelatedPatch [[0 : l.length() - 1]]; int i = 0; for (List_RelatedPatch r = l; r != null; r = r.rest()) { a[i] = r.first(); i++; } return a; } public boolean isEmpty() { return l == null; } public int length() { return l.length(); } } class List_CoarseFineBoundary { public CoarseFineBoundary car; public List_CoarseFineBoundary cdr; // Main constructor. public List_CoarseFineBoundary(CoarseFineBoundary x, List_CoarseFineBoundary rest) { car = x; cdr = rest; } public List_CoarseFineBoundary cons(CoarseFineBoundary x) { return new List_CoarseFineBoundary(x, this); } public /*inline*/ CoarseFineBoundary first() { return car; } public /*inline*/ List_CoarseFineBoundary rest() { return cdr; } public int length() { int i = 0; List_CoarseFineBoundary r = this; while (r != null) { i++; r = r.rest(); } return i; } } // The Boxed List class. Used to construct a list // element-by-element and convert it into an array. class BoxedList_CoarseFineBoundary { List_CoarseFineBoundary l = null; public BoxedList_CoarseFineBoundary() { } public BoxedList_CoarseFineBoundary(CoarseFineBoundary x) { push(x); } public void push(CoarseFineBoundary item) { l = new List_CoarseFineBoundary(item, l); } public CoarseFineBoundary pop() { CoarseFineBoundary t = l.first(); l = l.rest(); return t; } public CoarseFineBoundary [1d] toArray() { CoarseFineBoundary [1d] a = new CoarseFineBoundary [[0 : l.length() - 1]]; int i = 0; for (List_CoarseFineBoundary r = l; r != null; r = r.rest()) { a[i] = r.first(); i++; } return a; } public boolean isEmpty() { return l == null; } public int length() { return l.length(); } } class List_Patch { public Patch car; public List_Patch cdr; // Main constructor. public List_Patch(Patch x, List_Patch rest) { car = x; cdr = rest; } public List_Patch cons(Patch x) { return new List_Patch(x, this); } public /*inline*/ Patch first() { return car; } public /*inline*/ List_Patch rest() { return cdr; } public int length() { int i = 0; List_Patch r = this; while (r != null) { i++; r = r.rest(); } return i; } } // The Boxed List class. Used to construct a list // element-by-element and convert it into an array. class BoxedList_Patch { List_Patch l = null; public BoxedList_Patch() { } public BoxedList_Patch(Patch x) { push(x); } public void push(Patch item) { l = new List_Patch(item, l); } public Patch pop() { Patch t = l.first(); l = l.rest(); return t; } public Patch [1d] toArray() { Patch [1d] a = new Patch [[0 : l.length() - 1]]; int i = 0; for (List_Patch r = l; r != null; r = r.rest()) { a[i] = r.first(); i++; } return a; } public boolean isEmpty() { return l == null; } public int length() { return l.length(); } } class List_int_1d__1d_local { public int [1d] [1d] local car; public List_int_1d__1d_local cdr; // Main constructor. public List_int_1d__1d_local(int [1d] [1d] local x, List_int_1d__1d_local rest) { car = x; cdr = rest; } public List_int_1d__1d_local cons(int [1d] [1d] local x) { return new List_int_1d__1d_local(x, this); } public /*inline*/ int [1d] [1d] local first() { return car; } public /*inline*/ List_int_1d__1d_local rest() { return cdr; } public int length() { int i = 0; List_int_1d__1d_local r = this; while (r != null) { i++; r = r.rest(); } return i; } } // The Boxed List class. Used to construct a list // element-by-element and convert it into an array. class BoxedList_int_1d__1d_local { List_int_1d__1d_local l = null; public BoxedList_int_1d__1d_local() { } public BoxedList_int_1d__1d_local(int [1d] [1d] local x) { push(x); } public void push(int [1d] [1d] local item) { l = new List_int_1d__1d_local(item, l); } public int [1d] [1d] local pop() { int [1d] [1d] local t = l.first(); l = l.rest(); return t; } public int [1d] [1d] local [1d] toArray() { int [1d] [1d] local [1d] a = new int [[0 : l.length() - 1]] [1d] local [1d]; int i = 0; for (List_int_1d__1d_local r = l; r != null; r = r.rest()) { a[i] = r.first(); i++; } return a; } public boolean isEmpty() { return l == null; } public int length() { return l.length(); } } class List_Cut { public Cut car; public List_Cut cdr; // Main constructor. public List_Cut(Cut x, List_Cut rest) { car = x; cdr = rest; } public List_Cut cons(Cut x) { return new List_Cut(x, this); } public /*inline*/ Cut first() { return car; } public /*inline*/ List_Cut rest() { return cdr; } public int length() { int i = 0; List_Cut r = this; while (r != null) { i++; r = r.rest(); } return i; } } // The Boxed List class. Used to construct a list // element-by-element and convert it into an array. class BoxedList_Cut { List_Cut l = null; public BoxedList_Cut() { } public BoxedList_Cut(Cut x) { push(x); } public void push(Cut item) { l = new List_Cut(item, l); } public Cut pop() { Cut t = l.first(); l = l.rest(); return t; } public Cut [1d] toArray() { Cut [1d] a = new Cut [[0 : l.length() - 1]]; int i = 0; for (List_Cut r = l; r != null; r = r.rest()) { a[i] = r.first(); i++; } return a; } public boolean isEmpty() { return l == null; } public int length() { return l.length(); } } class List_RectDomain_3 { public RectDomain<3> car; public List_RectDomain_3 cdr; // Main constructor. public List_RectDomain_3(RectDomain<3> x, List_RectDomain_3 rest) { car = x; cdr = rest; } public List_RectDomain_3 cons(RectDomain<3> x) { return new List_RectDomain_3(x, this); } public /*inline*/ RectDomain<3> first() { return car; } public /*inline*/ List_RectDomain_3 rest() { return cdr; } public int length() { int i = 0; List_RectDomain_3 r = this; while (r != null) { i++; r = r.rest(); } return i; } } // The Boxed List class. Used to construct a list // element-by-element and convert it into an array. class BoxedList_RectDomain_3 { List_RectDomain_3 l = null; public BoxedList_RectDomain_3() { } public BoxedList_RectDomain_3(RectDomain<3> x) { push(x); } public void push(RectDomain<3> item) { l = new List_RectDomain_3(item, l); } public RectDomain<3> pop() { RectDomain<3> t = l.first(); l = l.rest(); return t; } public RectDomain<3> [1d] toArray() { RectDomain<3> [1d] a = new RectDomain<3> [[0 : l.length() - 1]]; int i = 0; for (List_RectDomain_3 r = l; r != null; r = r.rest()) { a[i] = r.first(); i++; } return a; } public boolean isEmpty() { return l == null; } public int length() { return l.length(); } }