/* Abstract unidimensional dyadic grids. */ /* Last edited on 2009-08-21 18:24:13 by stolfilocal */ #ifndef udg_grid_H #define udg_grid_H #define _GNU_SOURCE #include #include #include "interval.h" #include "bool.h" #include "sign.h" #include "vec.h" /* UNIDIMENSIONAL GRIDS A /unidimensional grid/ is a partition {G} of the reals {R} into open intervals and single real numbers, the /items/ of the grid, such that (1) one of the items is the single number {0}, and (2) all intervals have the the same length {delta} -- the /step/ of the grid. The singleton items and the interval items are also called the /vertices/ and /cells/ of the grid, respectively. Note that every item has the same geometry and topology as any other item of the same dimension. In particular, every vertex is incident to two cells, and every cell is bounded by two vertices. GRID POSITIONS A vertex {v} on an unidimensional grid can be identified by an integer /position/ {pos}, the displacement from the origin to {p} divided by the corresponding grid step {delta[i]}. In particular, the vertex at 0 has position 0. Similarly, a cell can be identified by the position of its inferior endpoint. In particular, the cell whose inferior endpoint is 0 (the /base cell/ of the grid) has position 0. */ typedef int64_t udg_grid_pos_t; /* Displacement of a grid item, in multiples of the corresponding grid step. Also a relative displacement between two items. */ /* CIRCULAR TOPOLOGY A /flat circle/ is an interval {[0_wd)}, for some positive real number {wd} (the /period/ or /extent/ of the flat circle), endowed with the topology that results from assuming that the superior endpoint {wd} is identical to the inferior endpoint {0}. It can be understood as the quotient of the real line {R} by an equivalence relation that equates any two reals which differ by an integer multiple of {wd}. CIRCULAR GRIDS A /circular grid/ is a grid on a flat circle whose period {wd} is some integer multiple {gsz*delta} of the grid's step {delta}. The integer {gsz} is called the grid's /size/. Thus a circular grid consists of a finite set {G} of items, taken from some (infinite) unidimensional grid, whose union {D = \U G} is the interval {[0 _ wd)} The interval {D} is called the /domain/ of the grid. Note that in a circular grid {G}, as in an infinite grid, all items of the same dimension are topologically alike. Namely, every vertex (including the origin 0) separates two cells, and every cell (including the last cell {[wd-delta _ wd)}) is bounded by two vertices. (However, if the grid has a single cell, then the two ends of the cell are the same vertex, and that vertex has the same cell on both sides.) */ typedef uint64_t udg_grid_size_t; /* Size of a circular grid, in multiples of the grid step. */ /* THE UNIDIMENSIONAL DYADIC GRID An (/infinite/) dyadic uni-dimensional multigrid/ {G*} is the union of an infinite collection of unidimensional circular 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/. 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. The half-cell that lies towards {-oo} along that axis is the /low/, /inferior/, or /0-/child; and the other half is the /high/, /superior/, or /1-/child. Note that all levels of a dyadic multigrid {G*} cover the same domain {D = \U G*}, namely the union of the root cell plus its inferior vertex, with circular topology. */ typedef int32_t udg_rank_t; /* Identifies a layer in a dyadic 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). */ sign_t udg_rank_cmp(udg_rank_t *x, udg_rank_t *y); /* A comparison function for ranks. */ udg_grid_size_t udg_grid_size(udg_rank_t r); /* The number of grid cells in level {r} of a unidimensional dyadic multigrid, namely {2^r}. */ udg_grid_pos_t udg_max_grid_pos(udg_rank_t r); /* Maximum position of an item in level {r} of a unidimensional dyadic multigrid. */ #define udg_MAX_RANK (64 - 1) /* Maximum rank of any level in a multigrid. The limit is such that a {udg_cell_index_t} (below) will fit in one 64-bit word. With this definition, 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, or 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 dyadic 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 within that layer. */ typedef uint64_t udg_cell_index_t; /* Identifies a cell in some multigrid. */ #define udg_NO_CELL (0LL) /* A cell index that means /no such cell/. */ #define udg_ROOT_CELL (1LL) /* The index of the root cell of any multigrid. */ #define udg_cell_min_index(r) (1LL << (r)) #define udg_cell_max_index(r) ((1LL << (r+1))-1LL) /* Cells in level {r} have conscutive indices in the range {[udg_cell_min_index(r) .. udg_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. */ udg_rank_t udg_cell_rank(udg_cell_index_t k); /* The rank of the cell whose index is {k}. */ #define udg_check_cell_rank(k,r) (((k)>>(r)) == 1LL) /* TRUE iff cell {k} belongs to level {r}, */ udg_grid_pos_t udg_cell_position(udg_cell_index_t k); /* Returns the position {pos} of cell {k} within its layer. */ udg_cell_index_t udg_cell_from_position(udg_rank_t r, udg_grid_pos_t pos); /* Returns the index of the cell of level {r} whose position is {pos}. */ udg_cell_index_t udg_cell_shift(udg_cell_index_t k, udg_grid_pos_t dp); /* The index of the cell that is displaced from cell {k} by {dp} grid steps. Takes into account the circular topology of each layer. */ /* ROOT-RELATIVE ARGUMENT A point {x} of {D} can identified by its /root-relative position/, which is 0 at the low end and 1 at the high end of the root cell. Root-relative coordinates that differ by an integer are implicitly equivalent, and therefore they are often reduced modulo 1 to real fractions in the range {[0_1)}. */ void udg_cell_interval_root_relative(udg_cell_index_t k, interval_t *B); /* Stores in {B[0..d-1]} the inferior and superior root-relative coordinates of the cell {k}, in a {d}-dimensional grid. The inferior bounds will be in the range {[0_1)}; the superior bounds will be strictly greater than the inferior bounds, in the range {(0_1]}. */ /* CELL PACKS A /cell pack/ is a sequence {\X} of consecutive cells in a given layer {r} of {G*}, defined by the layer index {r}, an /initial position/ {pos}, the position of the lowest cell in that layer, and a positive /pack size/ {psz}, the number of cells in the pack. Alternatively (and in several procedures below), a pack is specified by its size {psz} and by the cell index {b} of the inferior cell. Let {gsz} be the size of layer {r}. A cell {C'} of layer {r} belongs to the pack {\X} if and only if its position {pos'} satisfies {pos' = pos + e (mod gsz)}, where {e} is some integer in the range {0..psz-1}. The integer {e} is the /relative position/ or /slot index/ of {C'} within the pack. The lowest (inferior) cell of the pack has slot index 0, and the uppermost (superior) has slot index {psz-1}. */ typedef uint64_t udg_pack_slot_index_t; /* Identifies a cell slot in a cell pack. */ #define udg_NOT_SLOT ((udg_pack_slot_index_t)-1) /* The {udg_pack_slot_index_t} of a cell that is not in the pack. */ void udg_pack_cells(udg_cell_index_t b, udg_grid_size_t psz, udg_cell_index_t k[]); /* Returns in {k[0..psz-1]} the indices of all cells in the pack whose lower cell is {b} and whose size is {psz}. */ udg_pack_slot_index_t udg_pack_find_cell ( udg_cell_index_t k, udg_cell_index_t b, udg_grid_size_t psz ); /* Returns the slot index {pkx} of cell {k} in the cell pack with lower cell {b} and size {psz}. If {k} occurs multiple times in the pack, returns the first occurrence. If {k} is not part of the pack, returns {udg_NOT_SLOT} */ bool_t udg_pack_intersects_cell ( udg_cell_index_t k, udg_cell_index_t b, udg_grid_size_t psz ); /* Returns TRUE if cell {k} intersects some cell of the pack with lower cell {b} and size {psz}. */ void udg_pack_interval_root_relative ( udg_cell_index_t b, udg_grid_size_t psz, interval_t *B ); /* Stores in {B} the interval that is the union of all cells in the cell pack with lower cell {b} and size {psz}. The lower endpoint will be in the range {[0_1)}; the upper endpoint will be strictly greater. Due to the toroidal geometry, the upper endpoint {HI(B[i])} will be 1 or more whenever the position of {b} is {gsz-psz} or more, where {gsz} is the grid's size. */ /* MISCELLANEOUS */ sign_t udg_rank_cmp(udg_rank_t *x, udg_rank_t *y); /* Returns {-1} if {*x < *y}, {0} if {*x == *y}, {+1} if {*x > *y}. Useful for sorting things by rank. */ /* PRINTOUT */ void udg_cell_print(FILE *wr, udg_cell_index_t k); /* Prints the cell {k} in the format "{RANK}:{POS}" */ vec_typedef(udg_cell_index_vec_t,udg_cell_index_vec,udg_cell_index_t); /* A vector of {udg_cell_index_t}. */ #endif