/* Abstract dyadic multigrids in arbitrary dimension */ /* Last edited on 2004-08-24 00:16:14 by stolfi */ #ifndef dgbasic_H #define dgbasic_H #include #include #include #include /* Integer types: */ typedef int8_t int8; typedef uint8_t nat8; typedef int16_t int16; typedef uint16_t nat16; typedef int32_t int32; typedef uint32_t nat32; typedef int64_t int64; typedef uint64_t nat64; /* SPACE DIMENSIONS AND COORDINATE AXES */ typedef int32 dg_Dim; /* The dimension of a space, set, basis, etc.. */ typedef int32 dg_Axis; /* Identifies a coordinate axis. We assume that the axes of {R^d} are numbered from {0} to {d-1}. */ #define DG_MAX_GRID_DIM (4) #define DG_MAX_GRID_AXIS (DG_MAX_GRID_DIM-1) /* Just to be safe. High grid 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}. */ #define DG_NO_AXIS (DG_MAX_GRID_DIM) /* A {dg_Axis} value which means ``no such axis''. */ /* INTERVALS A {dg_Interval} {I} is a pair of {double}s {lo(I) = I.end[BLO]} and {hi(I) = I.end{BHI}}. It usually represents all real numbers {x} such that {lo(I) < x < hi(I)}. Depending on the context, the endpoints {lo(I)} and/or {hi(I)} may be included too. An open or half-open interval {I} is empty iff {lo(I) \geq hi(I)}. A closed interval {I} is empty iff {lo(I) > hi(I)}. */ typedef struct dg_Interval { double end[2]; } dg_Interval; /* A real interval, open or closed. If {end[0] > end[1]}, it represents the empty set. */ typedef enum { BLO = 0, BHI = 1 } dg_BinaryDir; /* Binary direction along an axis : {BLO} towards {-oo}, {BHI} towards {+oo}. Used to identify the two halves of a cell, the endpoints of an interval, etc. */ dg_Interval dg_interval_split(dg_Interval v, dg_BinaryDir dir); /* Returns the lower or upper half of {v}, depending on {dir}. */ dg_Interval dg_interval_join(dg_Interval *u, dg_Interval *v); /* Returns the smallest interval enclosing both {u} and {v}. */ dg_Interval dg_interval_meet(dg_Interval *u, dg_Interval *v); /* Returns the intersection of {u} and {v}. */ /* BOXES An /{m}-dimensional box of {R^d}/ is the Cartesian product {B} of {d} real sets {B_i}, {i = 0,..d-1}; of which {m} are open intervals, and {d-m} are singleton sets. A box {B} of {R^d} is here represented by an array of {dg_Intervals} {B[0..d-1]}, where {B[i]} is interpreted as open if {B[i].lo < B[i].hi}, and as closed (i.e., singleton) if {B[i].lo == B[i].hi}. */ void dg_box_lo_corner(dg_Dim d, dg_Interval B[], double p[]); /* Stores in {p[0..d-1]} the inferior corner of the box {B[0..d-1]}, i.e. sets {p[i] = lo(B[i])} for all {i}. */ void dg_box_hi_corner(dg_Dim d, dg_Interval B[], double p[]); /* Stores in {p[0..d-1] the superior corner of box {B[0..d-1]}. */ void dg_box_corner(dg_Dim d, dg_Interval B[], dg_BinaryDir dir[], double p[]); /* Stores in {p[0..d-1]} the corner of the box {B[0..d-1]} which, along each axis {i}, lies in the direction {dir[i]} -- namely, {p[i] = B[i].end[dir[i]]} for all {i}. */ double dg_box_size(dg_Dim d, dg_Interval B[]); /* Returns the maximum extent of box {B[0..d-1]} along any axis. */ void dg_box_join(dg_Dim d, dg_Interval A[], dg_Interval B[], dg_Interval C[]); /* Sets box {C[0..d-1]} to the smallest box enclosing both boxes {A} and {B}. */ void dg_box_meet(dg_Dim d, dg_Interval A[], dg_Interval B[], dg_Interval C[]); /* Sets box {C[0..d-1]} to the intersection of boxes {A} and {B}. */ /* BOX MAPPING AND UNMAPPING An {m}-dimensional box {B[0..d-1]} with spanning axes {s[0..m-1]} and stabbing axes {t[0..d-m-1]} defines a mapping from {R^m} to {R^d}, whereby each coordinate {z[j]} of the argument is scaled and shifted so that the interval [0_1] becomes the interval {x[s[j]] = B[s[j]]}, and {x[t[j]]} is set to the singleton {B[t[j]}. We will denote this transformation by {x = S_B(z)}. The numbers {z[0..m-1]} are the /coordinates of {x} relative to {B}/. */ void dg_point_box_map(dg_Dim d, double z[], dg_Interval B[], double x[]); /* Computes the point of of {R^d} which has coordinates {z[0..m-1]} relative to the {m}-dimensional box {B[0..d-1]}, and stores it in {x[0..d-1]}. */ void dg_point_box_unmap(dg_Dim d, double x[], dg_Interval B[], double z[]); /* Computes the cordinates of point {x[0..d-1]} relative to the {m}-dimensional box {B}, and stores them in {z[0..m-1]}. The point {x} is implicitly projected onto the affine span of {B}. */ void dg_box_box_map(dg_Dim d, dg_Interval Z[], dg_Interval B[], dg_Interval X[]); /* Computes the box of {R^d} which has coordinate ranges {Z[0..m-1]} relative to the {m}-dimensional box {B[0..d-1]}, and stores it in {X[0..d-1]}. */ void dg_box_box_unmap(dg_Dim d, dg_Interval X[], dg_Interval B[], dg_Interval Z[]); /* Computes the coordinate ranges of box {X[0..d-1]} relative to the {m}-dimensional box {B[0..d-1]}, and stores them in {Z[0..m-1]}. The box {X} is implicitly projected onto the affine span of {B}. */ /* BOX ORIENTATION A box {B} is/parallel/ to a coordinate axis {i} if its projection on that axis is a non-trivial interval; otherwise {B} is /perpendicular/ to axis {i}. We also say that axis {i} is /spanning/ or /normal/ to the box, respectively. We denote by {Spn(B)} the spanning axes of box {B}, and by {Nrm(B)} the normal ones. The set {Nrm(B)} is the /orientation/ of the box. */ typedef nat32 dg_Axes; /* A {dg_Axes} value {A} represents a set of axes, namely all axes {i} such that {A & 2^i != 0}. It can be used to represent the orientation of a box, etc.. */ #define dg_axis_belongs(ax,A) \ ((A >> ax) & 1) /* Returns TRUE iff axis {ax} belongs to axis set {A}. */ #define dg_axis_include(ax,A) \ ((A) | (1 << (ax))) /* Returns the axis set {A} augmented with axis {ax} . */ #define dg_axis_exclude(ax,A) \ ((A) & (~(1 << (ax)))) /* Returns the axis set {A} minus axis {ax} . */ int dg_axes_count(dg_Axes A); /* The cardinality of the axis set {A}. */ dg_Axes dg_axes_complement(dg_Dim d, dg_Axes A); /* The complement of the axis set {A}, namely the axes in {0..d-1} which are not present in {A}. */ typedef nat8 dg_AxisIndex; /* Relative index of an axis among a set {A} of axes: 0 = the smallest axis in {X}, 1 = second smallest, etc.. */ #define DG_AXIS_NOT_FOUND (DG_MAX_GRID_DIM) /* A {dg_AxisIndex} value which means ``no such axis''. */ dg_Axis dg_axis(dg_Axes A, dg_AxisIndex j); /* Returns the {j}th smallest axis in the set {A}, counting from {j=0} to {dg_axes_count(A)-1}. Returns {DG_NO_AXIS} if the set {A} has {j} or less elements. */ dg_AxisIndex dg_find_axis(dg_Axes A, dg_Axis ax); /* Returns the position of axis {ax} in the set {A}, sorted by increasing order. Returns {DG_AXIS_NOT_FOUND} if {ax} is not in the set {A}. */ /* FACETS AND FACES A /facet/ of an {m}-dimensional box {B} is a maximal box with dimension {m-1} contained in the boundary of {B}. There are {2*m} facets, two for each axis parallel to {B}: the /inferior/ facet, in the {-oo} direction from {B}, and the /superior/ one, in the {+oo} direction. The /faces/ of {B} are the box {B} itself, and the facets of all its faces. The /inferior/ (resp. /superior/) faces are the inferior (resp. superior) facets of all faces of {B}. An {m}-dimensional box has {3^m} faces, of which {2^d-1} are inferior and {2^d-1} are superior. RELATIVE FACE POSITIONS An {m}-dimensional box {B} has {3^m} faces. Let {Spn(B) = ax[0..m-1]} be the spanning axes of {B}. Each face {F} of {B} can be identified by its /signature/, a vector {dir[0..m-1]} where {dir[j]} specifies the projection {F_j} of {F} onto axis {ax[j]}, in terms of the projection {B_j} of {B} on that axis. Namely, {F_j = lo(B_j)} if {dir[j] = -1} {F_j = B_j} if {dir[j] = 0}, and, {F_j = hi(B_j)} if {dir[j] = +1}. In particular, if {dir = (0,.. 0)}, then {F} is {B} itself. An inferior face of a box has a signature consisting entirely of {0}s and {-1}s; the corresponding superior face has the {-1}s replaced by {+1}s. The signatures {(-1,.. -1)} and {(+1,.. +1)} specify the inferior and superior corners of {B}; respectively. */ typedef enum { SLO = -1, SMD = 0, SHI = +1 } dg_SignedDir; /* Three-valued direction along an axis: {SLO} towards {-oo}, {SHI} towards {+oo}, {SMD} in the middle/interior. Used to specify the position of a face relative to a box, etc. */ /* RELATIVE FACE INDICES The signature {dir[0..m-1]} of a face {F} relative to a box {B} can be packed into a single /relative face index/ {fi}, by the formula {fi = SUM{ 3^j * (dir[j] \bmod 3) : j = 0..m-1}}. Note that that the relative face index {fi} ranges from 0 to {3^m-1}. The face index {fi = 0} identifies the box {B} itself. The inferior corner has index {fi = (3^{m}-1)/2}; and the superior corner has index {fi = 3^{m}-1}. */ typedef nat dg_FaceIndex; /* Identifies one of the {3^m} faces of an {m}-dimensional box. */ dg_Dim dg_face_dimension(dg_Dim m, dg_FaceIndex fi); /* Returns the dimension the face of an {m}-dimensional box {B} that is identified by the relative face index {fi}. */ dg_SignedDir dg_face_position(dg_FaceIndex fi, dg_AxisIndex j); /* Given the index {fi} of a face {F} of some box {B}, returns the position of {F} relative to {B} along axis {Spn(B)[j]}. */ void dg_face_signature(dg_Dim m, dg_FaceIndex fi, dg_SignedDir dir[]); /* Given the index {fi} of a face {F} of an {m}-dimensional box {B}, stores in {dir[j]} the position of {F} relative to {B} along axis {Spn(B)[j]}. */ void dg_face_rel_box(dg_Dim m, dg_FaceIndex fi, dg_Interval F[]); /* Given the relative index {fi} of a face {F} of an {m}-dimensional box {F}, stores in each {F[j]} the coordinate range of face {fi} along axis {Spn(B)[j]}, for {j = 0..m-1}. The coordinates are relative to the box {B}, i.e. 0 means the low side of the box, 1 means the high side. Thus, {F[j]} will be set to either {0,0}, {1,1}, or {0,1}. Note that the first two mean singleton sets {\set{0}} and {\set{1}}, while the latter means the open interval {(0_1)}. */ void dg_face_print(FILE *wr, dg_Dim m, dg_FaceIndex fi); /* Prints the signature of face number {fi} of an {m}-dimensional box. */ /* GRIDS A /{d}-dimensional grid/ is a partition {G} of {R^d} into boxes of various dimensions, the /items/ of the grid, such that (1) the intersection of the closure of any two items of {G} is the union of items of {G}; (2) for any axis {i}, the projections of all items parallel to {i} are open intervals with the same length {\delta_i} -- the /step/ of the grid along that axis; and (3) one of the items consists of a single point, the origin of {R^d}. It follows that all items with the same spanning axes are congruent translates of each other, and that every item with dimension {m < d} lies on the boundary of {2^{d-m}} items with dimension {m+1}. A /cell/ of a grid is an item {C} with maximum dimension, i.e. {dim(C) == d}. Grid items with dimensions 0, 1, 2, 3 are called /vertices/, /edges/, /faces/, and /blocks/ of the grid, respectively. GRID POSITIONS A point {p} on a fixed grid can be identified by a /position vector/ {pos[0..d-1]}, where {pos[i]} is the displacement from the origin to {p} along the coordinate axis {i}, in multiples of the corresponding grid step {\delta_i}. In particular, the cell whose inferior corner is the origin (the /base cell/ of the grid) has position vector {(0,..0)}. A grid position vector can also be used to identify a cell (or tile) of the grid, or, more generally, any item with a fixed set of normal and spanning vectors: it suffices to give the position of the item's inferior corner. */ typedef int64 dg_GridPos; /* Displacement of a grid item along an axis, in multiples of the corresponding grid step. Also relative displacement between two items. */ /* TOROIDAL GRIDS A /toroidal grid/ is the quotient of a grid by an equivalence relation that equates any two points of {R^d} which, along every coordinate axis {i}, differ by an integer multiple of {W_i \delta_i}, where {W_i} is some positive integer. Note that {R^d} is then reduced to a {d}-dimensional torus {T^d} where {T} is the unit interval {[0_1)} taken with circular topology. In other words, a toroidal grid consists of a finite collection {G} of items of a {d}-dimensional grid, whose union {\U G} is some {d}-dimensional box plus all the inferior faces of that box; with the assumption that every superior face of that box coincides with the corresponding inferior one. The set {\U G} is called the /domain/ of the grid {G}. Note that, in a toroidal grid {G}, there is a one-to-one correspondence between the the cells of {G} and the items of {G} with any fixed orientation. */ typedef int64 dg_GridSize; /* Size of a toroidal grid along an axis, in multiples of the corresponding grid step. */ /* THE DYADIC MULTIGRID An (/infinite dyadic/) /{d}-dimensional multigrid/ {G*} is the union of an infinite collection of {d}-dimensional toroidal grids -- the /levels/ or /layers/ of {G*} -- each identified by a natural number, its /rank/. By definition, level 0 has a single cell, the /root cell/, whose inferior corner is at the origin. Every cell in level {r} contains exactly two cells of rank {r+1}, its /children/; and, except for the root, every cell is contained in exactly one cell of rank {r-1}, its /parent/. The two children of a cell are /sisters/ to each other, and the root cell is its own sister by definition. The children of a cell are obtained by splitting {C} in half, perpendicularly to axis {r \bmod d}. The half-cell that lies towards {-oo} along that axis is the /low child/ or /0-/child, and the other half is the /high child/ or /1-child/. Note that all levels of a multigrid {G*} cover the same domain {D = \U G*}, namely the union of the root cell plus all its inferior faces, with toroidal topology. */ typedef int32 dg_Rank; /* Identifies a layer in a multigrid, or the depth of a node in a binary tree, etc. It may also be a relative rank (i.e. difference of two ranks). */ dg_Axis dg_split_axis(dg_Dim d, dg_Rank r); /* The axis perpendicular to the cut that divides a {d}-dimensional cell of level {r} into its children cells; namely, axis {r \bmod d}. */ void dg_grid_size(dg_Dim d, dg_Rank r, dg_GridSize sz[]); /* Stores in {sz[0..d-1]} the number of grid cells along axis {i} in level {r} of the {d}-dimensional multigrid. */ dg_GridPos dg_max_grid_pos(dg_Dim d, dg_Rank r, dg_Axis i); /* Maximum grid position along axis {i} of level {r} in a {d}-dimensional multigrid. */ #define DG_MAX_RANK (64 - 1) /* Maximum rank of any level in a multigrid. The limit is such that a {dg_CellIndex} (below) will fit in one 64-bit word. In a 4-dimensional domain (for example, 3D space + time), the finest level has at least {2^15 = 32768} cells along each axis. That is enough for 1m resolution in a field 30 km wide, together with 5 min resolution in a 100-day simulation period. In a 3-dimensional domain (3D space, or 2D space + time), the finest level has {2^21 = 2097152} cells along each axis, which means 5 cm resolution in a field 100 km wide, together with 1 min resolution in a 1400-day simulation period. In a 2-dimensional domain (2D space, or 1D space + time), the finest level has at least {2^31 = 2147483648} cells along each axis, which means 1 mm resolution in a field 2000 km wide, together with 1 sec resolution in a 68-year simulation period. In a 1-dimensional domain, the finest level (rank 62) has {2^62 = 4.6×10^18} cells, which means 10 picometer (sub-atomic) resolution over the Earth's circumference, together with 1 millisecond resolution in a 150,000-year simulation period. That should be enough for most practical purposes. */ /* CELL INDICES We can uniquely identify a cell of a fixed multigrid {G*}, at any level, by a a single /cell index/ {k}. By definition, the root cell in layer 0 has index 1; and the children of cell {k} are the cells with indices {2*k} (low) and {2*k + 1} (high). The parent of cell {k} has therefore index {floor(k/2)}. Note that these formulas are independent of the grid's dimension {d}. These rules assign consecutive indices to all cells of each layer of {G*}; in layer {r} the indices range from {2^r} to {2^{r+1}-1}, inclusive. Therefore, the position of the leading 1 bit in the index of a cell is its rank. The remaining bits are the bits of the cell's position vector, interleaved: that is, bit {2^j} of {pos[i]} is bit {2^{d*j + i}} of the cell index {k}, for each axis {i} and each bit position {j} such that {d*j + i < r}. */ typedef nat64 dg_CellIndex; /* Identifies a cell in some multigrid. */ #define DG_NO_CELL (0) /* A cell index that means /no such cell/. */ #define DG_ROOT_CELL (1) /* The index of the root cell of any multigrid. */ #define dg_cell_min_index(r) (1 << (r)) #define dg_cell_max_index(r) ((1<<(r+1))-1) /* Cells in level {r} have conscutive indices in the range {[dg_cell_min_index(r) .. dg_cell_max_index(r)]}. The minimum index is the base cell of layer {r}, the maximum index is the cell adjacent to the superior corner of the domain. */ dg_Rank dg_cell_rank(dg_CellIndex k); /* The rank of the cell whose index is {k}. */ #define dg_check_cell_rank(k,r) \ (((k)>>(r)) == 1) /* TRUE iff cell {k} belongs to level {r}, */ void dg_cell_box_root_relative(dg_Dim d, dg_CellIndex k, dg_Interval B[]); /* Stores in {B[0..d-1]} the lower and upper coordinates of the cell {k}, relative to the root cell, in a {d}-dimensional grid. */ void dg_cell_position(dg_Dim d, dg_CellIndex k, dg_GridPos pos[]); /* For each axis {i} in {0..d-1}, stores in {pos[i]} the position of cell {k} within its layer. */ dg_CellIndex dg_cell_shift(dg_Dim d, dg_CellIndex k, dg_Axis ax, dg_GridPos dp); /* The index of the cell that is displaced from cell {k} by {dp} grid steps along axis {ax}. */ dg_CellIndex dg_cell_translate(dg_Dim d, dg_CellIndex k, dg_GridPos dp[]); /* Returns the index of the cell in a {d}-dimensional multigrid whose position vector is {dp[0..d-1]} relative to the inferior corner of cell {k}, within the layer {r} that contains {k}. The resulting position is implicitly reduced modulo {dg_max_grid_pos(d,r,i)}. */ dg_CellIndex dg_cell_index_add(dg_Dim d, dg_CellIndex ka, dg_CellIndex kb); /* The index of the cell whose position vector is the sum of the position vectors corresponding to indices {ka} and {kb}. The two cells should have the same rank. */ dg_CellIndex dg_cell_neighbor(dg_Dim d, dg_CellIndex k, dg_FaceIndex fi); /* The cell {k'} of rank as {k}, such that the highest-dimensional face shared by {k} and {k'} is face {fi} of {k}. In particular, if {fi == 0}, returns {k} itself. */ /* LOCUS IN A MULTIGRID Note that, in any {d}-dimensional multigrid {G*}, an item {B} with dimension {m < d} may occur in two or more consecutive levels. This situation happens whenever the splitting axis {r \bmod d} is in the set {Nrm{B}}. In particular, any vertex of layer {r} is also a vertex of layer {r+t}, for all {t > 0}. A /locus/ in a {d}-dimensional multigrid is an occurrence of a specific item at a specific level. We can designate a grid locus uniquely by giving its normal axes {Nrm(B)}, and its /owner cell/ {Cll(B,r)}, the {d}-dimensional cell of level {r} of which {B} is an inferior face along all those axes. (In particular, if {B} itself is a cell, then {Cll(B,r) = B}.) Note that the toroidal (periodic) topology of each level means that {Cll(B,r)} always exists, even for items contained in the superior faces of the domain {D}. In these comments, will use the notation {} for the locus consisting of the inferior face of the cell with index {k} with normal axes encoded by the integer {A}, at the level that contains cell {k}. Thus, for example, the root cell is locus {<0:1>}, and the vertex at the origin of layer {r} is locus {<2^d-1:2^r>}. */ typedef struct dg_Locus { dg_Axes norm; dg_CellIndex cell; } dg_Locus; /* Identifies a locus in a multigrid, namely locus {}. */ #define dg_locus(A,k) \ ((dg_Locus){(A),(k)}) /* The locus with orientation {A} which is the inferior face of cell {k} along those axes. */ #define dg_locus_rank(E) \ (dg_cell_rank(E.cell)) /* The rank of locus {E}. */ dg_Dim dg_locus_dimension(dg_Dim d, dg_Locus E); /* The dimension of locus {E} in a {d}-dimensional dyadic multigrid. */ dg_Axis dg_locus_normal_axis(dg_Dim d, dg_Locus E, dg_AxisIndex j); /* Returns the {j}th normal axis of locus {E}. */ dg_Axis dg_locus_spanning_axis(dg_Dim d, dg_Locus E, dg_AxisIndex j); /* Returns the {j}th spanning axis of locus {E}. */ dg_Locus dg_locus_minimize_rank(dg_Dim d, dg_Locus E); /* Returns the locus of minimum rank which consists of the same face of the grid as {E}. */ /* FACES AND SUB-FACES OF A LOCUS The inferior facet of locus {} along its spanning axis {i} is {}. The superior facet is {} where {k'} is the cell adjacent to cell {k} in the {+oo} direction along axis {i}. */ dg_Locus dg_locus_facet(dg_Dim d, dg_Locus E, dg_Axis ax, dg_SignedDir dir); /* Returns the facet of locus {E} that lies in the direction {dir} along axis {ax}. If {E} is orthogonal to axis {ax}, returns {E} itself. */ dg_Locus dg_locus_face(dg_Dim d, dg_Locus E, dg_FaceIndex fi); /* The face of locus {E} that is identified by the relative face index {fi}. */ dg_Locus dg_locus_neighbor(dg_Dim d, dg_Locus E, dg_FaceIndex fi); /* The locus {E'} of same dimension, rank and orientation as {E}, such that the largest common face of {E} and {E'} has index {fi} relative to {E}. In particular, if {fi == 0}, returns {E} itself. */ dg_Locus dg_locus_extend(dg_Dim d, dg_Locus E, dg_Axis ax, dg_SignedDir dir); /* Returns the result of sweeping locus {E} by one grid step along the axis {ax}, in the direction {dir}. If {dir == SMD}, or if {E} is parallel to axis {ax}, returns {E} itself. */ dg_Locus dg_locus_face_from_signature(dg_Dim d, dg_Locus E, dg_SignedDir dir[]); /* The face {F} of locus {E} that lies in the direction {dir[j]} along axis {ax[j]}, where {ax[0..m-1] = Spn(E)}. */ void dg_locus_box_root_relative(dg_Dim d, dg_Locus E, dg_Interval B[]); /* Stores in {B[0..d-1]} the lower and upper root-relative coordinates of the box described by locus {E} in a {d}-dimensional grid. */ /* STARS For each {m}-dimensional locus {E} at level {r} of the infinite dyadic grid, we define the /star of {E}/, {K(E)}, as being the list of {2^{d-m}} cells of level {r} whose boundary contains {E}. These cells are the /wings/ of {E}. WING INDICES The entries of {K(E)} are identified by a /wing index/ in {[0..2^m-1]} where {m = d - dim(E)} is the co-dimension of {E}. Let {C0,C1} be two cells of {K(E)} that are adjacent across a facet stabbed by axis {Nrm(E)[j]}, with {C0} on the low side of that facet; if the wing index of {C0} is {ix}, that of {C1} is {ix + 2^j}. Thus the locus {E} is the superior face of cell {K(E)[0]}, and the inferior face of cell {K(E)[2^m-1]}, along the axes {Nrm(E)}. DUPLICATED WINGS Thanks to the toroidal topology, every locus has a complete star. Note however that the star of a locus {E} of rank {r < d-1} will include the same cell more than once. More precisely, suppose that {dim(E) = m} and {Nrm(E)} includes {k} coordinate axes in the range {0 .. r}. Then the star {K(E)} actually has only {2^k} distinct cells, whose wing indices are {0..2^k-1}; these cells are then repeated {2^{m-k}} times in the list {K(E)}. */ typedef int32 dg_WingIndex; /* Identifies a cell slot in the star of some locus {E}. */ #define DG_NOT_WING ((dg_WingIndex)-1) /* The {dg_WingIndex} of a cell that is not a wing. */ void dg_locus_star(dg_Dim d, dg_Locus E, dg_CellIndex k[]); /* Returns in {k[0..2^m-1]} the indices of the cells of {K(E)}, in order of their wing indices, where {m} is the co-dimension of {E}. */ dg_WingIndex dg_cell_index_in_star(dg_Dim d, dg_CellIndex k, dg_Locus E); /* Returns the wing index of cell {k} among the cells of {K(E)}. If {k} occurs multiple times in {K(E)}, returns the first occurrence. If {k} is not part of {K(E)}, returns {2^d}. */ bool dg_cell_intersects_star(dg_Dim d, dg_CellIndex k, dg_Locus E); /* Returns TRUE if cell {k} intersects some cell of {K(E)}. */ void dg_star_box_root_relative(dg_Dim d, dg_Locus E, dg_Interval B[]); /* Stores in {B[0..d-1]} the lower and upper root-relative coordinates of the star of locus {E} in a {d}-dimensional grid. BEWARE: due to the toroidal geometry, the interval {B[i]} will extend into negative coordinates whenever {cell(E)} has position 0 along axis {i}. In particular, if {rank(E)