/* Last edited on 2004-08-17 23:11:15 by stolfi */ ---------------------------------------------------------------------- /* B-PULSES The dyadic B-pulses are defined for any {c \geq 0} and any degree {g = 2*c+1}. For a given {g} and {c}, there are {c+1} distinct mother pulses, with index {pix} ranging from 0 to {c}. The mother B-pulse {w^g_c[pix]} is supported on the interval {(-1_+1)}. In the positive half {[0_1)}, the mother pulse can be expressed as a combination of Bernstein-Bézier polynomials: { w^g_c[pix](z) = SUM { choose(c-k,pix-k)/2^{c-k} BB^g_{k}(z) : k = 0..pix } } For {z} in {[-1_0]}, the following relation holds { w^g_c[pix](z) = w^g_c[c-pix](-z) } It follows that the dyadic B-pulses of rank {r} and shift {j} are supported by the interval {[(j-1)/2^r_ (j+1)/2^r)}, reduced modulo 1. */ ---------------------------------------------------------------------- double dg_mother_B_pulse ( dg_Degree g, dg_Cont c, dg_PulseIndex i, double z, dg_GridSize sz ) { affirm(c >= -1, "bad contin for B-pulse"); affirm(g >= 0, "negative degree for B-pulse"); affirm(g >= 2*c + 1, "broad B-pulses not implemented yet"); affirm(g <= DG_MAX_B_PULSE_DEGREE, "degree too high for B-pulse"); affirm((i >= 0) && (i <= g - c - 1), "invalid index"); if (i > c) { /* A Bernstein-Bézier element with single-cell support: */ if ((z <= 0.0) || (z >= 1.0)) { return 0.0; } else { return dg_bezier_poly(g, i, z); } } else { if (sz == 1) { /* Take care of wrap-around: */ return dg_mother_bsp_poly(g, c, i, z) + dg_mother_bsp_poly(g, c, c - i, 1.0 - z); } else { /* Try to reduce to {[-1 _ +1]} modulo {sz}: */ if (z > sz-1) { z = z - sz; } /* Evaluate pulse: */ if ((z <= -1.0) || (z >= +1.0)) { return 0.0; } else { if (z >= 0.0) { return dg_mother_bsp_poly(g, c, i, z); } else { return dg_mother_bsp_poly(g, c, c - i, -z); } } } } } ---------------------------------------------------------------------- double dg_mother_poly ( dg_PulseKind kind, dg_Degree g, dg_Cont c, dg_PulseIndex i, double z ); /* Either {dg_mother_H_pulse_poly} or {dg_mother_bsp_poly}, depending on {kind}. */ double dg_mother_poly ( dg_PulseKind kind, dg_Degree g, dg_Cont c, dg_PulseIndex i, double z ) { switch(kind) { case DG_PK_H: affirm(g == 2*c + 1, "bad Hermite degree"); return dg_mother_H_pulse_poly(c, i, z); case DG_PK_B: return dg_mother_bsp_poly(g, c, i, z); default: affirm(FALSE, "unknown pulse kind"); return 0.0; /* To keep the compiler happy. */ } } ---------------------------------------------------------------------- /* GENERAL GRIDS An /orthogonal set/ is a subset of {R^d} which is the union of finitely many {d}-dimensional boxes and all their faces. The /membership function/ of an orthogonal set {X}, also denoted by {X}, is a function from {R^d} to {[0 _ 1]} such that {X(p)} is zero if {p} lies in the exterior of {X}, and {k*2^{-d}} otherwise, where {k} is the maximum number of disjoint {d}-dimensional boxes that are contained in {X} and have {p} as a corner. We define the /restriction of a function {f} to a box {B}/ as {f*B}, the pointwise product of {f} and the membership function of {B}. */ ---------------------------------------------------------------------- bool dg_node_star_is_whole(dg_NodeStar *N); /* Returns TRUE if the star of the locus {E} represented by the node star {N} is contained in the grid's domain, i.e. if there are no missing leaves in any subtree {N[i]}. */ bool dg_node_star_is_whole(dg_NodeStar *N) { int i; for (i = 0; i < N->nel; i++) { dg_Node *Ni = N->el[i]; if (Ni == NULL) { return FALSE; } } return TRUE; } ---------------------------------------------------------------------- if (sz == 1) { /* Take care of wrap-around: */ if (kind == DG_HER) { double s = (i % 2 == 0 ? 1.0 : -1.0); return dg_mother_her_poly(g, c, i, z) + s * dg_mother_her_poly(g, c, i, 1.0 - z); } else if (kind == DG_BSP) { return dg_mother_bsp_poly(g, c, i, z) + dg_mother_bsp_poly(g, c, c - i, 1.0 - z); } else { affirm(FALSE, "invalid spline kind"); return 0.0; } } else { /* Try to reduce to {[-1 _ +1]} modulo {sz}: */ if (z > sz-1) { z = z - sz; } /* Evaluate spline: */ if ((z <= -1.0) || (z >= +1.0)) { return 0.0; } else { if (z >= 0.0) { return dg_mother_poly(kind, g, c, i, z); } else { if (kind == DG_HER) { double s = (i % 2 == 0 ? 1.0 : -1.0); return s * dg_mother_her_poly(g, c, i, -z); } else if (kind == DG_BSP) { return dg_mother_bsp_poly(g, c, c - i, -z); } else { affirm(FALSE, "invalid spline kind"); return 0.0; } } } } ---------------------------------------------------------------------- ---------------------------------------------------------------------- /* 1/2*3^(1/2) */ #define HER_0_0_NORM ( 0.86602540378443864675 ) /* 1/26*455^(1/2) */ #define HER_1_0_NORM ( 0.82041265414236699041 ) /* 1/2*105^(1/2) */ #define HER_1_1_NORM ( 5.1234753829797991915 ) /* 1/362*83622^(1/2) */ #define HER_2_0_NORM ( 0.79882510411311308942 ) /* 3/52*5005^(1/2) */ #define HER_2_1_NORM ( 4.0815014205747676248 ) /* 1/2*2310^(1/2) */ #define HER_2_2_NORM ( 24.031229681395831560 ) ---------------------------------------------------------------------- dg_StarIndex dg_cell_index_in_star(dg_Dim d, dg_CellIndex k, dg_Locus E) { dg_Rank kr = dg_cell_rank(k); dg_Rank Er = dg_cell_rank(E.cell); if (kr != Er) { return (1 << d); } dg_GridPos kp[DG_MAX_GRID_DIM]; dg_cell_position(d, k, kp); dg_GridPos Ep[DG_MAX_GRID_DIM]; dg_cell_position(d, E.cell, Ep); int i, j; dg_StarIndex ix = 0; for (i = 0, j = 0; i < d; i++) { if (dg_axis_belongs(i, E.norm)) { if (kp[i] == Ep[i]) { ix += (1 << j); } j++; } } return ix; } ---------------------------------------------------------------------- /* !!! Should redefine the star to avoid repeating cells when there is wraparound. Namely define the /double axes/ of {E} as the axes in {0..d-1} along which the star has two distinct cells. An axis {ax} is double if {ax \in Nrm(E)} and the level of {E} has more than one cell along {ax}. */ ---------------------------------------------------------------------- dg_Locus_vec_t dg_find_critical_loci(dg_Dim d, dg_Node *root, bool vtxOnly) { affirm(root->index = 1, "must start at root"); dg_Locus_vec_t crit = dg_Locus_vec_new(100); int ncrit = 0; auto void dg_find_critical_sub_loci(dg_Locus E, dg_Rank r, dg_NodeStar *NE); void dg_find_critical_sub_loci(dg_Locus E, dg_Rank r, dg_NodeStar *NE) { fprintf(stderr, " enter rank = %d locus = ", r); dg_locus_print(stderr, d, E); fprintf(stderr, "\n"); if (! dg_node_star_is_valid(r, NE)) { /* Locus is not in grid; backtrack. */ fprintf(stderr, " [invalid locus: star = "); dg_node_star_print(stderr, d, NE); fprintf(stderr, "]\n"); } else if (dg_node_star_is_leaf(NE)) { /* Leaf locus; save unless {vtxOnly} and not vertex. */ fprintf(stderr, " [leaf locus: "); affirm(dg_node_star_is_whole(NE), "incomplete splines not supported yet"); if ((! vtxOnly) || (dg_locus_dimension(d, E) == 0)) { if (ncrit >= crit.nel) { dg_Locus_vec_expand(&crit, ncrit); } crit.el[ncrit] = E; ncrit++; fprintf(stderr, " collected"); } else { fprintf(stderr, " not vertex, ignored"); } fprintf(stderr, "]\n"); } else { /* Locus is not a leaf locus; split and recurse. */ fprintf(stderr, " [nonleaf locus: "); dg_Axis splax = dg_split_axis(d, r); /* Splitting axis. */ /* Generate subtrees: */ if (dg_axis_belongs(splax, E.norm)) { /* The split axis is perpendicular to the item. */ fprintf(stderr, " shrink star]\n"); dg_AxisIndex splj = dg_find_axis(E.norm, splax); dg_Locus E1 = dg_locus(E.norm, 2*E.cell); /* Same item in next level. */ /* Shrink node star around element: */ dg_NodeRef_vec_t NE1; dg_node_star_shrink(NE, splj, &NE1); /* Recurse on same item in the next level, with same shape. */ dg_find_critical_sub_loci(E1, r+1, &NE1); free(NE1.el); } else { /* The split axis is parallel to item {E}, split it: */ fprintf(stderr, " bisect locus:"); /* Compute normal axes of splitting item: */ dg_Axes NrmMD = dg_axis_include(splax, E.norm); /* Find relative position {splj} of {splax} among the {Nrm(EMD)}: */ dg_AxisIndex splj = dg_find_axis(NrmMD, splax); affirm(splj < d, "inconsistent axis sets"); fprintf(stderr, " split on axis %d = Nrm(EMD)[%d]\n", splax, splj); fprintf(stderr, "]\n"); /* Pieces of this locus: */ dg_Locus ELO = dg_locus(E.norm, 2*E.cell); dg_Locus EHI = dg_locus(E.norm, 2*E.cell + 1); dg_Locus EMD = dg_locus(NrmMD, 2*E.cell + 1); /* Node stars of the pieces: */ dg_NodeRef_vec_t NLO, NHI, NMD; dg_node_star_split(NE, splj, &NLO, &NHI, &NMD); /* Recurse for each half of item {F}: */ dg_find_critical_sub_loci(ELO, r+1, &NLO); dg_find_critical_sub_loci(EHI, r+1, &NHI); /* Recurse for the separating item {S} of {F}. */ dg_find_critical_sub_loci(EMD, r+1, &NMD); /* Reclaim temporary storage. */ free(NMD.el); free(NLO.el); free(NHI.el); } } fprintf(stderr, " exit rank = %d locus = ", r); dg_locus_print(stderr, d, E); fprintf(stderr, "\n"); } /* Start recursion with root cell and its faces: */ { /* Allocate a root {dg_NodeStar} of sufficient size: */ int max_ncells = (1 << d); dg_StarIndex ix; dg_NodeRef_vec_t NE = dg_NodeRef_vec_new(max_ncells); for (ix = 0; ix < max_ncells; ix++) { NE.el[ix] = root; } /* Enumerate the faces of the root cell: */ dg_Axes nrm; for (nrm = 0; nrm < (1<= 2, "invalid point coordinates"); /* Compute mean function values normalized to [-1 _ +1]: */ for (j = 0; j < n-2; j++) { /* Compute mean function value: */ f[j] = (a[j+2]+b[j+2]+c[j+2])/3.0; /* Clip to {[fMin _ fMax]}: */ if (f[j] < fMin) { f[j] = fMin; } if (f[j] > fMax) { f[j] = fMax; } /* Normalize to [-1 _ +1]: */ if (f[j] > 0.0) { f[j] /= fMax; } else if (f[j] < 0.0) { f[j] /= -fMin; } } /* Compute color {R,G,B} for white-background plotting: */ double R0 = 1.0, G0 = 0.342, B0 = 0.000; if (n == 2) { R = 1.0; G = 0.9; B = 0.8; } else if (n == 3) { pswr_color_scale_1(f[0], R0, G0, B0, &R, &G, &B); } else if (n == 4) { pswr_color_scale_2(f[0], f[1], R0, G0, B0, &R, &G, &B); } else if (n >= 5) { pswr_color_scale_3(f[0], f[1], f[2], &R, &G, &B); } pswr_set_fill_color(ps, R,G,B); pswr_triangle(ps, a[0], a[1], b[0], b[1], c[0], c[1], TRUE, FALSE); } ---------------------------------------------------------------------- /* For color shading, at most three function values {z[0..m-1]} are considered, i.e {p[2..n-1]} where {m = min(3,n-2)}. Each function component {z[j]} is clipped to the range {[fMin[j] _ fMax[j]}, then normalized to a value {s[j]} in {[-1 _ +1]}. Specifically, if {z[j]} (after clipping) is negative, then {s[j] = -z[j]/fMin}, else {s[j] = z[j]/fMax}. The vector {s} is converted to a color value with {ps_map_to_color}. */ /* If {m = 1} or {m = 2}, the Euclidean length {sm = |s|} of {s} determines the brightness, by the formula {log(y) = sm * log(y0)} where {y0} is the brightness of {r,g,b}. Note that {s = (0,..)} implies {y = 1} (white), and {s != (0,..)} implies {y < 1}. If {m = 3}, the The color is then {(y,y,y) + m*(s[0]*A + s[1]*B)}, where {A} is the zero-brightness vector {(r,g,b)-(y0,y0,y0)} and {B} is the zero-brightness vector orthogonal to {A} and ahead of {A} spectrumwise, with {|B| = |A|} in the YUV model. The factor {m} is a number no greater than 1, depending only on {s[1]/s[0]}, such that the combination { If {m = 1}, let {fu = s[0]} Positive values of {fu} are mapped to the line through {1,1,1} and {r,g,b}, negative values to the line through {1,1,1} and {r',g',b'} which is the complement of {r,g,b} relative to the gray of same brightness. If {m = 2}, let {fu = s[0],fv = s[1]} at each point. The modulus {f = hypot(fu,fv)} determines the brightness, as above. The ratios {fu/f,fv/f} are mapped to a point on the plane of colors with that brightness, in such a way that {fv = 0} implies a color on the line through {1,1,1} and {r,g,b}, and {fu = 0} implies a color 90 degrees ahead of that, spectrumwise. If {m = 3}, let {fu = s[0],fv = s[1],fw = s[2]}. The In the last three cases, the color scales are such that the function value {(0,...)} gets mapped to WHITE, not black, and higher values mean deeper colors. If {n = 3}, there is a single function value {f} at each point. The value is clipped to the range {[-fMax _ +fMax]} and converted to a color scale where the given color is {fMax}, white is 0, and the complementary color of the same brightness is {-fMax}. For {n = 4}, there are two function values {fu,fv} at each point. These values are clipped to the unit disk, then mapped to a two-dimensional color scale. The isolines show the modulus {hypot(fu,fv)}. If {n >= 5}, the first three function values {fu,fv,fw} at each point are considered, and the rest is ignored. The function values are clipped to the range {[0 _ fMax]}, then turned into color values by an {RGB}-like convention. The isolines show the average value {(fu+fv+fw)/3}. In the last three cases, */ dg_Rank rank = 0; for (i = 0; i < d; i++) { z[i] = x[i]; } { dg_CellIndex cell = t.cell; i = 0; while (cell > 1) { if ((cell & 1LL) == 1) { z[i] -= 0.5; } z[i] *= 2.0; i = (i + 1) % d; rank++; } } /* If the rank of {t.cell} is less than {d}, the star {K(t)} may include the same cell {C} several times, with different star indices. If the point {x} belongs to one such cell, we evaluate the tent for each star index {six}, and /* Reduce {z[i]} to 0-centered range modulo the period: */ for (i = 0; i < d; i++) { double wdi = (double)dg_max_grid_pos(d, rank, i); if (z[i] > wdi/2) { z[i] -= wdi; } } double dg_tent_general_eval ( dg_Dim d, dg_Degree g, dg_Cont c, dg_Tent t, double *z ) { int ax; /* Extract Bézier index tuple {tg} from tent indices {ix0,ix1}: */ int M; if ((M = (g - 2*c - 1)) > 0) { dg_BezIndex ix0 = t.ix0; dg_Axes SpnE = dg_axes_complement(d, t.E.norm); dg_Axis ax = 0; while (SpnE > 0) { if ((SpnE & 1) == 1) { tg[ax] = c+1 + (ix0 % M); } ax++; ix0 /= M; } } else { affirm(t.ix0 == 0, "bad ix0"); } if ((M = (c + 1)) > 0) { dg_BezIndex ix1 = t.ix1; dg_Axes NrmE = t.E.norm; dg_Axis ax = 0; while (NrmE > 0) { if ((NrmE & 1) == 1) { tg[ax] = ix1 % M; } ax++; ix1 /= M; } } else { affirm(t.ix1 == 0, "bad ix1"); } /* Evaluate the Bernstein-Bézier polynomial on {z}: */ double fx = 1.0; for (ax = 0; ax < d; ax++) { affirm(FALSE, "not implemented yet"); } return fx; } /* If point is outside the tent's support, return 0: */ for (ax = 0; ax < d; ax++) { if (dg_axis_belongs(ax, t.E.norm)) { if (fabs(z[ax]) > 1.0) { return 0.0; } } else { if ((z[ax] > 1.0) || (z[ax] < 0.0)) { return 0.0; } } } affirm(g >= 2*c + 1, "overconstrained tents not implemented yet"); if ((c == -1) && (g == 0)) { affirm(t.E.norm == 0, "C-1 constant tent on non-cell locus"); affirm(t.ix0 == 0, "bad ix0 in C-1 constant tent"); affirm(t.ix1 == 0, "bad ix1 in C-1 constant tent"); return 1.0; } else if ((c == 0) && (g == 1)) { double fx = 1.0; affirm(t.E.norm == ((1 << d) - 1), "C0 multiaffine tent on non-vert locus"); affirm(t.ix0 == 0, "bad ix0 in C0 multiaffine tent"); affirm(t.ix1 == 0, "bad ix1 in C0 multiaffine tent"); for (ax = 0; ax < d; ax++) { fx *= 1.0 - fabs(z[ax]); } return fx; } else { return dg_tent_general_eval(d,g,c,t,z); } if ((g == 1) && (c == 0)) { affirm(dg_locus_dimension(d, t.E) == 0, "wrong dim"); double val = 1.0; int i; for (i = 0; i < d; i++) { val *= (ix && (1 << i) == 0 ? z[i] : 1.0 - z[i]); } return val; } else { /* For now, return something nonzero: */ double dx = z[0]-0.5, dy = z[1]-0.5, r = hypot(dx,dy); return cos(6*r); double compute_plot_depth ( dg_Interval rootCell[], /* The root cell (client units). */ double scale, /* Plot scale (mm per client unit). */ double meshSize /* Nominal plot step (mm). */ dg_Rank *minDepth, /* (OUT) Always split until this level. */ dg_Rank *extraDepth, /* (OUT) Split leaves this many levels. */ dg_Rank *maxDepth /* (OUT) Never go below this level. */ ) { double dx = rootCell[0].end[1] - rootCell[0].end[0]; double dy = rootCell[1].end[1] - rootCell[1].end[0]; double diam = scale*hypot(dx, dy); int dp = (int)ceil(log(diam/meshSize)/log(2.0)); if (dp < 0) { dp = 0; } fprintf(stderr, "compute_plot_depth: dp = %d\n", dp); *minDepth = dp - 4; *extraDepth = 4; *maxDepth = DG_MAX_RANK; /* Always plot the whole tree. */ } /* Compute zero-brightness vector orthogonal to {uR,uG,uB}: */ double vR = (0.587*uB - 0.114*uG)/0.668555; double vG = (0.114*uR - 0.299*uB)/0.668555; double vB = (0.299*uG - 0.587*uR)/0.668555; dg_Tent_vec_t dg_minimal_tent_basis(dg_Tree *t) { nat ntv = 0; dg_Tent_vec_t tv = dg_Tent_vec_new(100); int_vec_t perp_axis = int_vec_new(t->d); auto void dg_find_tents(dg_Rank r, nat s, dg_NodeStar *NE); void dg_find_tents(dg_Rank r, nat s, dg_NodeStar *NE) { int i, mask; int nCells = NE->nel; dg_Node **CX = (dg_Node **)malloc(nCells*sizeof(dg_Node *)); dg_Axis long_axis = r % t->d; bool none_is_leaf; for( i=0; ich[0] && NE[i]->ch[1]) ) { none_is_leaf = FALSE; } if( none_is_leaf ) // There are no leaves { bool parallel_split = TRUE; // Check parallelism for( i=0; ich[0]; dg_find_tents ( r+1, s, CX ); for( i=0; ich[1]; dg_find_tents ( r+1, s, CX ); perp_axis.el[s] = long_axis; for( i=0; ich[0]; CX[2*i+1] = NE[i]->ch[1]; } dg_find_tents ( r+1, s+1, CX ); } else { for( i=0; ich[0]; else CX[i] = NE[i]->ch[1]; } dg_find_tents ( r+1, s, CX); } } else if( s == t->d ) { dg_Tent_vec_expand(&tv, ntv); tv.el[ntv] = (dg_Tent){NE[0]->index, r}; ntv++; } } dg_find_tents ( 0, 0, &(t->root) ); dg_Tent_vec_trim(&tv, ntv); return tv; } /* MINIMAL DYADIC TENTS IN A FINITE DYADIC GRID */ /* A dyadic tent {t} is said to be `compatible' with a grid {G} if each of its bricks is the union of one or more whole cells of {G}. A compatible tent {t} is `minimal' if there is no other dyadic tent function compatible {G} whose support is properly contained in that of {t}. If {G} is a finite dyadic grid, it is easy to see that a dyadic tent {t} is minimal if, and only if, all its bricks are cells (leaf or non-leaf) of {G}, of the same rank {r}, and at least one of them is a leaf cell. It then follows that the support {B} of {t} is congruent to a dyadic cell of rank {r-d}, translated so that it is centered at the common vertex {p} of the bricks. (Note that {B} itself may not be a cell of the dyadic cell tree; it may it may straddle two or more such cells.) It is conjectured that, for any dyadic grid {G}, the set of all the unit minimal tent functions of {G} constitute a basis for the space {S^1_0(G)}. */ /* MULTIAFFINE INTERPOLATION An important special case of the basis construction above is {c = 0}, {g = 1} -- the /continuous multilinear splines/. A tensor polynomial of degree 1 is an affine function of each coordinate separtely (and degree {d} along a general direction). This is a case where the only stars that give rise to basis elements are the vertex stars. Each vertex star contributes a single basis element, a /tent function/, which is 1 at the central vertex of the star, and decays linearly to 0 on the periphery of the star. BICUBIC INTERPOLATION Another important special case is {c = 1}, {g = 3} --- the /smooth bicubic splines/. Here too, the only stars that matter are vertex stars. Each star contributes 2^d basis elements, which can be labeled with {d}-bit strings. The element with bit vector {e_0,.. e_{d-1}} controls the derivatives {\pder_1^{e_1} \pder_2^{e_2} .. f} of {f} at the vertex {v}. Note that if such a function is zero at two points of an axis-aligned line segment, it is zero on the whole segment. Thus, in particular, if a spline {f} from {S^1_0(L)} is zero at the two endpoints of some edge {e} of the grid, it is zero along the entire edge {e}. If {f} is zero on two opposite edges of some face {H} of the grid, it is zero on the whole face {H}. And so on. /* */ /* TENT FUNCTIONS */ /* A `tent function' is a spline {t} in {S^1_0(Q)}, where {Q} is a grid consisting of {2^d} axis-aligned boxes {b[0..2^d-1]}, of the same shape and size, that share a common corner {p}. By definition, within each of these boxes (the `bricks' of {t}), the tent function is a single non-zero multilinear function, which decreases from some central value {t(p)} to zero on the boundary of the tent's support {B} (the union of its bricks). A `unit tent function' is a tent function that has value 1 at the central point. Observe that all the tent functions with the same set of bricks {b[i]} are multiples of the corresponding unit tent function. A `dyadic tent' is a tent funtion whose bricks are cells of the dyadic cell tree. */ /* TENT FUNCTION EVALUATION */ /* The evaluation of a unit tent function {t} at a point {p} of {R^d} is best divided three steps. First, we compute the `support-relative' coordinates {x} of {p}, by the formula {x[i] = (p[r+i]-c[r+i])/h[r+i]}, where {r} is the rank of the bricks, {c} is the center of the tent's support, and {h[j]} is the extent along axis {j} of the bricks of {t}. (All indices are taken modulo {d}.) Note that the support-relative coordinate {x[0]} is the relative position of {p} along the longest axis of the bricks. Then, if any of the {x[i]} lies outside the range {[-1 _ +1]}, the value of {t(p)} is zero. Otherwise, the value is simply the product of {1 - |x[i]|}, for all {i}. In practice, the grid may be deformed by a complicated shape function {f}, which needs to be inverted before computing {t(p)}. So, part of the computation of {x} is probbaly best merged with the computation of {f^{-1}(p)}. */ /* It follows that the If a cell {C} in {before applying these rules to a set of adjacent pieces whose tiles have different rank, the larger piece must be subdivided, by DeCasteljau's algorithm, until the sizes agree. BOUNDARY CONDITIONS Note that since the functions of {S(F,L)} are defined to be zero outside the diagram's domain {\U L}, a function from {S_c(F,L)} will be zero to order {c} at any point of the boundary bertween {\U L} and {C0 \setminus \U L}. To model functions which have continuity {c >= 0} within {\U L}, but are unconstrained on the boundary of {\U L}, one must choose a complete dyadic diagram {L'} that contains {L}, and work with the space {S(F,L')}, ignoring any portions that lie outside {L}. In that case, one should chose the grid sizes so that the extent of {\U L} along any axis is atmost 1/2 of the extent of {C0}; otherwise the periodic topology of {C0} will be felt inside {L}. In favourable cases, a basis for the set {S_c(P^g,L)} can be formed of the /primitive elements/ of {G(L)}, which in turn are derived from a suitable set of /stars/. PRIMITIVE ELEMENTS We say that the star {K(E)} of an item {E} /respects/ {L} if the intersection of {\U K(E}} with {\U L} is the union of tiles of {L}. In that case, we say that {K(E}} /fits in/ {L} if {\U K(r,E}} is contained in {\U L}. Let {c} and {g} be such that {g \geq 2 c + 1}. For each {E} whose star {K(E)} fits in {L} we will now define a set {H(E)} of splines whose support is {K(E)}, the /primitive elements/ associated with {E}. These splines are linearly independent among themselves and from all primitive elements {H(E´)} where {E´} ???? ??? {N(E) = (c+1)^{d-m}*(g-2*c-1)^m} splines of {S_c(P^g,L)} {K(E)} elements for each {m}-dimensional item {E} which has a star that hits {L}. Note that if {g = 2*c+1} then {N(E)} is zero for any {m > 0}; that is, we only need to consider the vertices of {L}. Let {K = K(r,E)} be such a star (If there is more than one star for the same item {E}, choose one of them arbitrarily.) Let {t} be the inferior tile of {K}. Among the array of Bézier coefficients of {t}, let {\Gamma_{K}} be the set of those which are within {c} steps of {E}, inclusive, but at least {c+1} steps away from any face of {E}. Let {\Delta_{K}} be the set of coeffs which satisfy these same conditions in the other tiles of {K}. For each coeff in {\Gamma_{K}}, construct a spline where that coeff is 1, and the rest is 0; coeffs not in {\Gamma_{K} \U \Delta_{K}} are 0; the coeffs in {\Delta_{K}} are computed by the continuity constraints. MINIMAL-SUPPORT BASES FOR DYADIC POLYNOMIAL SPLINES A minimal-support basis for the set {S_c(P^g,L)} can be constructed as follows. For each {m}-dimensional item {E} of level {r} of the infinite dyadic grid, we define the /star of rank {r} at {E}/, {K(r,E)}, as being the {2^{d-m}} tiles of level {r} whose boundary contains {E}. (This definition assumes that layer {r} has toroidal topology.) We are now ready to describe the construction of the basis for {S_c(P^g,L)}: * For {c = -1} and any degree {g \geq 0}, the basis has {(g+1)^d} elements for each tile {t} of {L}. Each of these elements has one of the Bézier coefficients set to 1, all the rest set to 0. * For any {c \geq 0} and {g \leq c}, the space {S_c(P^g,L)} is trivial, and the basis is empty. * For any {c \geq 0}, {g \geq 2*c + 1}, * For other combinations of {c} and {g}, the problem is harder; must think harder. Conjecture: the basis described above is a basis for the space {S_c(P^g,L)}. */ /* A diyadic tent function is completely determined by the center point {c} and the rank {r} of its bricks. The center {c}, in turn, is uniquely identified by the index of the the cell of which {cx} is the center [[??? Or the index of lowest-numbred brick ???]]. Thus, for example, the largest possible unit dyadic tent, whose support is the entire root cell, spans the {2^d} bricks of rank {d}, whose indices are {2^d .. 2^{d+1}-1}; thus it can be represented by the pair { cx = 1, r = d } [[??? or {cx = 2^d, r = d} ???]]. */ dg_Dim d; /* Dimension of grid. */ dg_Axis gax[DG_MAX_DIM]; /* Axes which are perpendicular and parallel to {S}. */ dg_Rank r; /* Rank of grid cells adjacent to {S}. */ dg_Index_vec_t ka; /* Indices of those cells. */ dg_GridIndex dg_cell_grid_index(dg_CellIndex k) { return k - (1 << dg_cell_rank(k)); } typedef struct dg_Tree { dg_Dim d; /* Dimension of the grid's domain (positive). */ dg_Node *root; /* Root node of the tree. */ } dg_Tree; /* Header for a finite tree of {dg_Node}s. */ typedef nat8 dg_Height; /* The `height' of a finite, non-empty binary tree is one plus the maximum rank of any node relative to the root node. In particular, a tree with just one node has height 1. The empty tree, by definition, has height 0. */ dg_Tree *dg_Tree_new(dg_Dim d); /* Creates a new {d}-dimensional tree with a single node (the root). */ void dg_Tree_free(dg_Tree* t); /* Recursively reclaims the whole tree {t}, including the header. */ dg_Tree *dg_Tree_new(dg_Dim d) { void *v = notnull(malloc(sizeof(dg_Tree)), "no mem for dg_Tree"); dg_Tree *h = (dg_Tree*)v; h->d = d; h->root = dg_Node_new(NULL, 1); return h; } void dg_Tree_free(dg_Tree* t) { if (t != NULL) { dg_free_subtree(t->root); free(t); } } /* The procedures in this section are concerned with the sub-space {S^1_0(L)} of {S^1(L)}, consisting of the *continuous* splines of degree 1 for some finite *dyadic* grid {L}. Note that since the union of all cells of {L} is the root cell {C_0}, any such spline is by definition zero outside {L}, and (by continuity) also on the boundary of {L}. */ typedef struct dg_Star /* A star of some item of a dyadic multigrid. */ { dg_Item E; /* Central item. */ } dg_Star; /* A star of an element {E}, consisting of the point {c} and the rank {r} of its bricks. The center {c}, in turn, is uniquely identified by the index of the the cell of which {cx} is the center [[??? Or the index of lowest-numbred brick ???]]. Thus, for example, the largest possible unit dyadic tent, whose support is the entire root cell, spans the {2^d} bricks of rank {d}, whose indices are {2^d .. 2^{d+1}-1}; thus it can be represented by the pair { cx = 1, r = d } [[??? or {cx = 2^d, r = d} ???]]. */ /* ITEM EXTENSION For the following procedures, {F-} is some face of an {m}-dimensional index {E-}, both orthogonal to a given coordinate axis {ax}. Also, {E+} is the {ax}-extension of {E-}, that is, the item which results from sweeping {E-} forward along axis {ax} by a whole grid step; and {F+} is {ax}-extension of {F-}. The following procedures relate the face index of {F-} relative to {E-} to the face index of {F+} relative to {E+}. */ dg_FaceIndex dg_face_project(dg_FaceIndex fi, dg_AxisIndex j); /* Given the face index {fi} of {F+} relative to {E+}, returns the index of {F-} relative to {E-}, where the extension axis is {ax = Spn(E)[j]}. */ dg_FaceIndex dg_project_face(dg_FaceIndex fc, dg_Axis ax) { dg_FaceIndex h = ipow(3,ax); return ((fc/h)/3)*h + (fc%h); } /* RELATIVE FACE INDICES */ dg_Dim dg_face_dimension(dg_FaceIndex f, dg_Dim m) { dg_Dim d = 0; while (m > 0) { if ((f % 3) == 0) { d++; } f /= 3; m--; } return d; } dg_SignedDir dg_face_position(dg_FaceIndex f, dg_Dim m, dg_Axis ax) { dg_FaceIndex h = ipow(3,ax); return ((f/h + 1) % 3) - 1; } void dg_face_signature(dg_FaceIndex f, dg_Dim m, dg_SignedDir loc[]) { int i; for (i = m-1; i >= 0; i--) { loc[i] = ((f + 1) % 3) - 1; f /= 3; } } /* A {dg_Node} may also be interpreted the point that is the center of that cell. (Note that there is at most one cell of the dyadic cell tree which is centered at a given a point {p} of {R^d}.) Then the leaf nodes, in this interpretation, correspond to the centers of undivided cells of the former interpretation; and internal nodes nodes correspond to centers of divided cells. If we interpret each node as a dyadic cell, then the tree represents a hierarchical partition of the root cell into dyadic cells of various sizes. */ #define DG_MAX_INDEX (1152921504606846975L) #define DG_MAX_TREE_HEIGHT (60) #define DG_MAX_RANK (DG_MAX_TREE_HEIGHT-1) /* This package is designed to handle finite binary trees whose node indices are in the range {0 .. 2^60-1}. Thus the maximum height of an actual tree is 60, and the maximum node rank is 59. */ #define DG_MAX_AXIS_HEIGHT(d) (DG_MAX_TREE_HEIGHT/(d)) #define DG_MAX_AXIS_RANK(d) (DG_MAX_AXIS_HEIGHT(d)-1) /* If the tree is used to represent a grid in 4-dimensional space (for instance, 3D space + time), the finest achievable grid has {2^15 = 32768} cells along each axis. That is enough for 1m resolution in a field 30 km wide, together with 10 min resolution in a 200-day simulation period. In three dimensions (3D space, or 2D space + time), the finest achievable grid has {2^20 = 1048576} cells along each axis, which means 10 cm resolution in a field 100 km wide, together with 1 min resolution in a 700-day simulation period. */ #define DG_MAX_DIM (4) /* Just to be safe. High dimensions {d} are dangerous because too many things have size {c^d} for some {c>1}. In particular, a {dg_FaceIndex} ranges from {0..3^d-1}. */ dg_Rank dg_rank(dg_Index k) { int r = 1; while (k > 1) { k = (k/2); r++; } return r; } /* Vectors of {dg_Axis}: */ typedef struct dg_Axis_vec_t { nat nel; dg_Axis *el; } dg_Axis_vec_t; dg_Axis_vec_t dg_Axis_vec_new(nat nel); #define dg_Axis_vec_expand(nv,index) \ vec_expand(vec_cast_ref(nv), index, sizeof(dg_Axis)) #define dg_Axis_vec_trim(nv,nel) \ vec_trim(vec_cast_ref(nv), nel, sizeof(dg_Axis)) dg_Axis_vec_t dg_Axis_vec_new(nat nel) { /* This is not a macro only because gcc does not allow cast of struct: */ vec_t v = vec_new(nel, sizeof(dg_Axis)); dg_Axis_vec_t r; r.nel = v.nel; r.el = (dg_Axis *)v.el; return r; } /* A {dg_FaceIndex} {f} is an integer that identifies a particular face {F} of some dyadic cell {C}. It encodes the sign vector {loc[0..d-1]} that defines the face, according to the formula {f = SUM {(loc[i]%3) * 3^i }}. Thus, for example, the face index {f = 0} encodes the sign sequence {(0,..0)}, which means the interior of the cell, i.e. {(0 _ 1)^d}. The lower corner of the cell, with sign sequence {(-1,.. -1)}, has face index {f = 3^d-1}; the upper corner, with sign sequence {(+1,.. +1)}, has face index {f = (3^{d+1}-1)/2}. The two facets perpendicular to coordinate axis axis {i} have face indices {f = 3^i} ({SHI} direction) and {f = 2*(3^i)} ({SLO} direction). A {dg_FaceIndex} {s} can also be used to specify a sub-face {S} of some {m}-dimensional face {F}, with {m < d}, provided that one has decided on a specific ordering {ax[0..m-1]} of the {m} coordinate axes which are parallel to {F}. In that case the face index {s} ranges over {0..3^m-1}, and element {loc[i]} of the the signature encoded by {s} is the position of {S} relative to {F} along axis {ax[i]}. */ /* DYADIC CELL FACES */ /* A /face/ of a {d}-dimensional cell {C} is a subset of the closure of {C} which, in cell-relative coordinates, can be written as {I[0] \x I[1] \x ... \x I[d-1]}, where each {I[i]} is either the single value {0}, or the single value {1}, or the open interval {(0 _ 1)}. Each cell has therefore {3^d} faces. Each face can be identified by its /signature/ {dir[0..d-1]}, a string of {d} signs {-1,0,+1}, where {dir[i]} tells whether {I[i]} is the single value {0} if {dir[i] = -1}, the single value {1} if {dir[i] = +1}, and the open interval {(0 _ 1)} if {dir[i] = 0} --- all in coordinates relative to the cell {C}. The only {d}-dimensional face of a {d}-dimensional cell {C} is the interior of {C}, i.e. the set {(0 _ 1)^d}, which has signature {(0,.. 0)}. There are {2*d} faces of dimension {d-1}, the /facets/ of {C}; and {2^d}} faces of dimension 0, the /corners/ of {C}. In general, there are {2^k * \choose(d,k)} faces of dimension {d-k}. If {E} is an {m}-dimensional dyadic face of some cell {C}, then the /sub-faces/ of {E} are the {3^m} dyadic faces of {C} contained in the closure of {E}. Given a specific ordering {ax[0..m-1]} of the coordinate axes parallel to {E}, then any sub-face {S} of {E} can be identified by its /relative signature/, a string of {m} signs {dir[0..m-1]}, where {dir[i]} specifies the projection of {S} on axis {ax[i]}, under the above convention. */ /* Compute non-flatness of sub-faces of {FLO,FHI}: */ for (f = 0; f < nFaces; f++) { shLO->dev.el[f] = 1.0e200; shHI->dev.el[f] = 1.0e200; } for (f = 0; f < nFaces; f++) { dg_SignedDir fsign = dg_face_position(f, m, ax); double devf; /* Inherit or compute the non-flatness of sub-face {f} of {FLO}: */ if (shLO->dev.el[f] == 0.0) { /* {f} is a subface of a face which was previously found to be flat: */ devf = 0.0; } else if (fsign == SLO) { /* {f} lies in the low face of {FLO}, inherit from {F}: */ devf = sh->dev.el[f]; } else if (fsign == SHI) { /* {f} lies in the high face of {FLO}, inherit from {FMD}: */ devf = shMD->dev.el[dg_project_face(f, ax)]; } else if (sh->dev.el[f] == 0.0) { /* {f} is part of a flat face of {F}, it is flat too: */ devf = 0.0; } else { devf = dg_Bezier_try_flatten_face(&(shLO->b), f, tol); } shLO->dev.el[f] = devf; if (devf == 0.0) { dg_propagate_flatness(&(shLO->dev), f, m); } /* Inherit or compute the non-flatness of sub-face {f} of {FHI}: */ if (shHI->dev.el[f] == 0.0) { /* {f} is a subface of a face which was previously found to be flat: */ devf = 0.0; } else if (fsign == SHI) { /* {f} is on the high face of {FHI}, inherit from {F}: */ devf = sh->dev.el[f]; } else if (SLO) { /* {f} is on the low face of {FHI}, inherit from {FMD}: */ devf = shMD->dev.el[dg_project_face(f, ax)]; } else if (sh->dev.el[f] == 0.0) { /* {f} is part of a flat face of {F}, it is flat too: */ devf = 0.0; } else { devf = dg_Bezier_try_flatten_face(&(shHI->b), f, tol); } shHI->dev.el[f] = devf; if (devf == 0.0) { dg_propagate_flatness(&(shHI->dev), f, m); } } dg_Degree *dg_uniform_degrees(dg_Degree g, dg_Dim d); /* Returns a constant vector of at least {d} degrees, all set to {g}. This vector may be shared so it should not be altered. */ #define DG_MAX_DEGREE 6 static dg_Degree **dg_unif_cache = NULL; /* Cache of uniform degree vectors. */ dg_Degree *dg_uniform_degrees(dg_Degree g, dg_Dim d) { affirm(d <= DG_MAX_DIM, "invalid domain dimension"); affirm(g <= DG_MAX_DEGREE, "invalid degree"); if (dg_unif_cache == NULL) { int i; dg_unif_cache = (dg_Degree **)notnull(malloc(DG_MAX_DEGREE*sizeof(dg_Degree*)), "no mem"); for (i = 0; i < DG_MAX_DEGREE; i++) { int j; dg_unif_cache[i] = (dg_Degree *)notnull(malloc(DG_MAX_DIM*sizeof(dg_Degree)), "no mem"); for (j = 0; j < DG_MAX_DIM; j++) { dg_unif_cache[i][j] = i; } } } return dg_unif_cache[g]; } The domain of {F} is the face of the unit cube which is perpendicular to axes {ax[0..s-1]}, which must be an ordered subset of {0..m-1}, where {m = b.m}. Along axis {ax[i]}, {F} lies in the direction {dir[i]}. Axes {ax[s..m-1]} should be the axes parallel to {F}, in some order. (I.e., {ax[0..m-1]} should be a permutation of {0..m-1}.) Axis {i} of the domain of {t} corresponds to axis {ax[s+i]} of the domain of {b}. In particular, a call with {s = 0} will permute the axes of {b}'s domain in the order given by {ax}. /* The dimension of the grid should not be compiled-in! The same library should be usable for grids of different dimensions! */ #define DIM 3 /* Changed false and true to FALSE and TRUE for compatibility with other libs. */ typedef enum{false, true} bool; /* It's not a good idea to tie the dimension indices to specific axes. */ /* Dimensions should be just integers [0..d-1] where d is the dimension of the grid. */ typedef enum{X, Y, Z, T} dimensions; #define NUM_AXIS 2 typedef enum{x_axis, y_axis} axis; #define next_axis(P) ((P+1)%NUM_AXIS) /* You sould never compare two floats for equality. Eps-equality is bogus --- there is no "right" value for eps. */ #define PRECISION 1e-10 #define EQUAL(A,B) (fabs(A - B) < PRECISION) #define GT(A,B) ((A - B) > PRECISION) #define LT(A,B) ((B - A) > PRECISION) /* Bad name! Use an array of intervals, of varaible length. */ typedef struct dyg_dimensions { double min[DIM]; double max[DIM]; }dyg_dimensions; /* Bad name: DYG_MAX_DEPTH should be 60. The name should be DYG_MAX_AXIS_DEPTH #define DYG_MAX_DEPTH(d) (60/d) /* "initial" is spelled wrong. */ /* By definition the initial dimension is always 1 × s × s^2 × ... × s^{d-1} */ dyg_dimensions* inicial_dimension; /* By definition the initial axis is always 0 */ int inicial_axis; typedef bool dg_EdgePred ( dg_Rank r, dg_Axis a, dg_Bezier *b ); /* A predicate used when plotting edges of a 2D grid. */ void dg_plot_2D_grid_edges ( dg_Rank r, /* Rank of initial cell. */ dg_Index k, /* Index of initial cell. */ dg_Node *p, /* Node of initial cell, or NULL. */ dg_Bezier *b, /* Shape of initial cell, or NULL. */ bool border, /* TRUE plots outer border of initial cell. */ dg_EdgePred plotEdge, /* Plots edge, or tells to split further. */ ); /* Enumerates the edges of all the cells in the subtree rooted at cell {k}, calling {plotEdge} on them. If {border} is {TRUE}, the procedure will consider all the edges of all these cells, taking care not to visit the same edge twice. if {border} is {FALSE}, the procedure will consider only edges which are interior to the cell {C}, and omit those which lie on the boundary of {C}. The predicate {plotEdge(r1,a1,kLO,kHI,pLO,pHI,b1)} is called once for each edge {e1} considered in the recursion. Here {b1} is a Bézier curve, i.e. a patch with {b1-> = 1}, which defines the shape of that edge, and {a1} is the coordinate axis which is transverse to edge in the nominal coordinate system. The argument {r1} is the rank of the two cells of which this edge is a side; {kLO,kHI} are their indices (in increasing order along the axis {a}), and {pLO,pHI} are the corresponding tree nodes (if they exist). If {e1} is a border edge, either {kLO} or {kHI} will be 0. */ void dg_plot_2D_grid_edges ( dg_Rank r, /* Rank of cell. */ dg_Index k, /* Index of cell. */ dg_Node *p, /* Node of tree. */ dg_Bezier *b, /* Shape of cell. */ bool border, /* TRUE to draw the outer border of the cell. */ dg_EdgePred plotEdge, /* Plots face, or tells to split further. */ ) { affirm((b == NULL) || (b->m == 2), "wrong shape dimension"); affirm(k > 0, "bad cell index"); /* Plot outer edges, if so requested: */ if (border) { dg_plot_2D_cell_side(r, k, 0, LO, NULL, p, b, plotEdge); dg_plot_2D_cell_side(r, k, 0, HI, p, NULL, b, plotEdge); dg_plot_2D_cell_side(r, k, 1, LO, NULL, p, b, plotEdge); dg_plot_2D_cell_side(r, k, 1, HI, p, NULL, b, plotEdge); } /* Plot dividing edge, then plot insides of children: */ !!! Cannot work -- must have a plotFace here too !!! if (p->c[LO] != NULL) { dg_Axis a = dg_longest_axis(2, r); affirm(p->c[HI] != NULL, "broken tree"); dg_Node *pLO = p->c[LO]; dg_Node *pHI = p->c[HI]; dg_Bezier *bLO = dg_Bezier_new(b->d, b->d, b->g); dg_Bezier *bHI = dg_Bezier_new(b->d, b->d, b->g); dg_Bezier_split(b, 0, bLO, bHI); dg_plot_2D_cell_side(f, bLO, r+1, 0, HI, pLO, pHI, tol); dg_plot_2D_tree(f, pLO, r+1, bLO, FALSE); dg_plot_2D_tree(f, pHI, r+1, bHI, FALSE); } } void dg_plot_2D_cell_side ( FILE *f, dg_Bezier *b, /* Shape of cell in question. */ dg_Rank r, /* Rank of cell in question. */ dg_Axis a, /* Draw side of {b} that is transverse to axis {a}. */ dg_Direction dir, /* Which side to draw, {LO} or {HI}. */ dg_Node *pLO, /* Tree node on {LO} side of edge to plot. */ dg_Node *pHI, /* Tree node on {HI} side of edge to plot. */ double tol /* Maximum error in plotted coords. */ ) { dg_Bezier *bside = dg_Bezier_new(1, b->n, b->g); /* Extract Bézier description of edge: */ dg_Bezier_extract_slice(b, a, (dir == LO ? 0 : b->g), bside); dg_plot_2D_edge(f, bside, r, a, NULL, p, tol); } void dg_plot_2D_edge ( FILE *f, dg_Bezier *b, /* Shape of edge segment. */ dg_Rank r, /* Rank of cells (even if sub-leaf). */ dg_Axis a, /* Edge is transverse to axis {a}. */ dg_Node *pLO, /* Tree node on {LO} side, or {NULL}. */ dg_Node *pHI, /* Tree node on {HI} side of edge, or {NULL}. */ double tol /* Maximum error in plotted coords. */ ) { /* When drawing an edge, cells must be split recursively until small enough to plot, even below the fringe of the tree. To ensure that cells join properly, even when the tolerance {tol} is large, we must artificially subdivide any actual leaf cell of the grid which shares some edge with two or more grid cells. For that purpose, we pass along the two tree nodes of rank {r}, {pLO} and {pHI}, which are separated by the edge being plotted, if there are any. As long as either of these nodes exists and is not a leaf, we must continue subdividing. Eventually the nodes {pLO} and/or {pHI} become {NULL}, and we lose information about the identity of the adjacent cells. We know only the shape of the edge {b}, the rank {r} of those sub-cells, and that the edge is transverse to axis {a}. */ if ( ((p == NULL) || (p->c[LO] == NULL)) && ((q == NULL) || (q->c[LO] == NULL)) && ) { Compute straight line approx: if (err <= tol) { draw segment; return } } /* Must split and recurse: */ { dg_Axis asplit = dg_longest_axis(2, r); affirm(r < MAX_RANK, "recursion too deep"); if (asplit == a) { /* Split cells, but not edge: */ dg_plot_2D_edge(f, b, r+1, a, pLO->c[HI], pHI->c[LO], tol); } else { /* Split cells and edge: */ dg_Bezier *bLO = dg_Bezier_new(1, b->d, b->g); dg_Bezier *bHI = dg_Bezier_new(1, b->d, b->g); dg_Bezier_split(b, asplit, 0.5, bLO, bHI); dg_plot_2D_edge(f, bLO, r+1, a, pLO->c[LO], pHI->c[LO], tol); dg_plot_2D_edge(f, bHI, r+1, a, pLO->c[HI], pHI->c[HI], tol); } } } void dg_enum_faces ( dg_Dim d, /* Dimension of grid. */ dg_Rank r, /* Rank of initial cell. */ dg_Index k, /* Index of initial cell. */ bool *fIncl, /* Tells which faces of the cell {k} to include. */ dg_Bezier *b, /* Shape of initial cell. */ double tol, /* Maximum allowed deviation from flatness. */ dg_Node *n, /* Node of initial cell, or NULL. */ dg_FaceProc process /* Client procedure that uses a face. */ ) { int_vec_t ax = int_vec_new(d); /* A permutation of the axes {0..d-1}. */ { int nFaces = ipow(3,d); /* Number of faces in a {d}-dimensional cell. */ bool_vec_t flat = bool_vec_new(nFaces); int f; /* Flatten any faces of {b} that can be flattened: */ for (f = 0; f < nFaces; f++) { flat.el[f] = FALSE; } flatten(b, flat); /* Now process faces that were requested by client: */ for (f = 0; f < nFaces; f++) { if (fIncl[f]) { /* Extract its Bézier patch and its {flat} vector, build the adjacent cell array {k}, {n}, } double_vec_t corner = double_vec_new(ncorners*n); /* Choose sensible values for corners: */ for (i = 0; i < ncorners; i++) { double p[m]; /* Get coordinates {p[0..m-1]} of corner {i} of the unit {m}-cube: */ int t = i; for (j = 0; j < m; j++) { p[j] = (double)(t % 2); t /= 2; } /* Map corner by some map: */ for (j = 0; j < m; j++) { int mask = ipow(2,j); p[j] = (double)(i & mask); } /* Compute