/* Mapping items of a dyadic multigrid into some range space. */
/* Last edited on 2005-02-21 11:27:53 by stolfi */

#ifndef dg_itemmap_H
#define dg_itemmap_H

#include <dg_grid.h>
#include <dg_tree.h>
#include <bz_patch.h>
#include <affirm.h>
#include <stdlib.h>
   
/*
   ITEM MAPS
   
   An /item map/ is a function from some item of a dyadic multigrid
   into some range space {R^n}, together with some derived
   information. If {mp} is a map for some item {E}, we denote the
   map's image of a point {p} of {E} by {mp(p)}
   
   An item map {mp} for an item {E} is typically defined by a Bézier
   patch {mp.b}. The axes of the domain of the patch correspond to the
   spanning axes {Spn(E)} of the item, in increasing order; i.e.,
   {mp(p) = mp.b(x)} where argument {x[j]} of the patch is the
   coordinate of {p} relative to {E} along axis {Spn(E)[j]}. Note that
   {x[j]} varies from 0 to 1 as {p} traverses the item {E} in the
   positive direction.
   
   Item maps are pasted together to build smooth functions.  We can see 
   them as 
   
   Item maps may be used, in particular, to define the true geometry
   of the dyadic multigrid. In that case, some or all of the
   coordinates of the range space {R^n} are interpreted as a point in
   the true domain of the grid.
   
   SHARING CONTROL POINTS
   
   For maps 
   
   To implement these procedures without wasting too much work
   and storage, we store in each item map {mp} only the 
   /interior/ control points of the batch {mp.b}, namely those
   that are not associated with any face.  In other words,
   
   FACE SNAPPING
   
   Many applications, such as plotting or integration, demand that
   complicated (high-degree) item maps be approximated by simpler
   (low-degree) ones, with subdivision if necessary. 
   
   The approximation must be done with care in order to maintain a
   desired level of continuity between the maps of adjacent items;
   failure to do so may result in ``cracks'' or ``foldovers'' in the
   approximating mesh.  Basically, (1) whenever we decide to replace
   the map of an item {E} by an approximation, we must modify
   the maps of all items which are faces of {E} or which have
   {E} as a face.  Moreover, (2) once an item is approximated, that
   approximation must be used for all its sub-items
   in lower levels of the multigrid.
   
   On the other hand, in a recursive algorithm, (3) the decision to
   approximate each item must be made independently for each item map,
   without looking at the maps of neighboring items. In order to
   satisfy requirements (1), (2), and (3), the decision to approximate an
   item map {mp} by a lower-degree one must depend exclusively on the
   item's rank and the patch {mp.b} itself. In particular, the
   decision must not depend on properties of its ancestors in the
   recurson tree.
   
   
   a /semi-flattened/ version 
   {b'} of the original patch {b}.  Whenever a face of {b'}
   is found to be flat enough, 
   
   Therefore, an
   item map {mp} also carries an approximation {b'} of the Bézier
   patch {b}, that is as simple as possible.
   
   This process is called /snapping/, and consists in detecting 
   faces of {b} that are almost flat, and forcing them to be
   perfectly flat. Here /flat/ means that {b(x) = b#(x)}
   for all points {x} on that face, where {b#} is the 
   multiaffine patch with same corners as {b}; and /almost
   flat/ means that the distance between {b(x)} and {b#(x)}
   is smaller than a given threshhold. 
   
   
   
   {f} that was found to be
   sufficiently close to the corresponding face of {bFlat}
   was forced to coincide with it.  Note that snapping 
   preserves the corners of the patch; and when a face is snapped.
   
 */

typedef struct dg_item_map_t   /* Map and location of an item {E}. */   
  { dg_Item E;            /* The multigrid item. */ 
    bz_Patch bTrue;       /* The original geometry of {E}. */
    bz_Patch bSnap;       /* The partially flattened geometry of {E}. */
    double_vec_t dev;     /* Non-flatness of each sub-face of {E}. */
    dg_node_ref_vec_t ta;  /* Tree nodes associated to adjacent cells. */
  } dg_item_map_t;
  /* A {dg_item_map_t} describes the geometry of an item {E} 
    in a {d}-dimensional dyadic multigrid.  
    
    The original geometry of item {E} is given by the Bézier patch
    {bTrue}. The vector {dev} has one element for each
    face {f} of {E}, and records the /deviation/ of that face -- that is, an
    upper bound for the distance between {bTrue(x)} and {bFlat(x)},
    for any point {x} in the face's domain; where {bFlat(x)}
    is the multiaffine Bézier patch with the same corners as {bTrue}.
    
    The patch {bSnap} contains a version of {bTrue} that has
    been /snapped/ --
    where {f} and
    the corresponding point of a multiaffine patch with same dimension
    and same corners as {f}. The vector {dev} satisfies the invariant:
    {dev[h] \leq dev[f]}, for any face {f} of {E}, and any sub-face
    {h} of {f}.
    
    In particular, a face {f} with {dev[f] == 0} is flat -- that is,
    coincides with that multiaffine patch; and ist sub-faces are all
    flat, too.
    
    The field {ta} is related to the set {Adj(E)} of all the 
    cells of the grid whose closure contains {E}.  If {E} has
    {p} normal axes, then there are {2^p} such cells, {Adj(E)[0..2^p-1]}.
    These cells are listed in such an order that cells {[i]}
    and {[i+2^j]} are separated by a {d-1}-dimensional plane
    orthogonal to axis {Nrm(E)[j]}, with the former on the {-oo} side
    of that plane. This property holds for each {j = 0..p-1} and each
    {i} in {0..2^p-1} such that {i & 2^j == 0}. In particular, if the
    face {E} is a full-dimension cell (i.e., {m == d}), then there is
    only one cell in the list {ka.el[0]} is the index of that cell.
    
    Fields {ta} is a list {ta[0..2^p-1]} of tree nodes, such that
    {ta[j]} , if not NULL, is the lowest node in some finite cell tree
    that contains the corresponding cell {Adj(E)[j]}. Note that the
    rank of node {ta[j]} may be less than the rank {r} of cell
    {Adj(E)[j]}. */

dg_item_map_t dg_item_map_t_new
  ( bz_DDim m,
    bz_RDim n,
    const bz_Degree *g
  );
  /* Allocates a new {dg_item_map_t} record {sh}, with given dimensions 
    and degrees. The Bézier patch {sh->b} and the adjoining cell lists
    {sh->ka} and {sh->ta} are allocated but left undefined. */

dg_item_map_t dg_item_map_t_facet_new
  ( bz_DDim m,
    bz_RDim n,
    const bz_Degree *g,
    dg_axis_t jax
  );
  /* Allocates a new {dg_item_map_t} record {sh}, with dimensions 
    and degrees suitable to represent a facet {E} of a cell {C},
    orthogonal to axis {jax} among the domain axes of {C}.
    
    The procedure assumes that the parent cell {C} has domain dimension {m},
    range dimension {n}, and its map is a Bézier patch
    with degrees {g[0..m-1]}.  The facet {E} has therefore dimension
    {m-1}, and its degree vector is {g} with {g[jax]} excluded. */

void dg_item_map_t_free(dg_item_map_t sh);
  /* De-allocates the internal storage nodes of {sh} (but not {sh} itself). */

dg_item_map_t *dg_map_for_cell
  ( bz_Patch *b, 
    dg_rank_t r,
    dg_cell_index_t k, 
    dg_node_t *t,
    double tol
  );
  /* Returns a {dg_item_map_t} record {sh} for a full-dimensional cell
    {C} of rank {r} and index {k}, whose geometry is given by the
    Bézier patch {b}. Node {t}, if not NULL, must be the
    lowest node in some tree whose cell contains {C}.
    
    Also computes the deviation {sh->dev[f]}, for each face {f} of
    {C}. Each face that is flat to within tolerance {tol} is forced to
    be perfectly flat. */

void dg_map_split
  ( dg_item_map_t *sh,
    dg_axis_t ax,
    double tol, 
    dg_item_map_t *shLO,
    dg_item_map_t *shHI,
    dg_item_map_t *shMD
  );
  /* Splits a {dg_item_map_t} in half, by a plane perpendicular to axis
    number {ax} of its domain. Allocates and returns the maps
    {shLO,shHI} of the two half-faces (with the same dimensions as
    {sh}) and the map {shMD} of the separating cut (with axis {ax}
    excluded from its domain). */

void dg_map_get_facet
  ( dg_item_map_t *sh, 
    dg_axis_t jax, 
    dg_BinaryDir dir, 
    dg_item_map_t *shF 
  );
  /* Extracts from the face map {sh} the facet that lies in
    direction {dir} along axis {jax}, and stores it in {shF}. The
    deviation {shF->dev[h]} of each sub-face {h} of {shF} is copied
    from {sh->dev[f]}, where {f} is the index of that same subface
    among the faces of {sh}. */
    
void dg_map_print_dev(FILE *wr, double_vec_t dev, bz_DDim m, char *fmt);
  /* Given the deviation vector {dev} of an {m}-dimensional face
    {f}, writes it to file {wr}, using the {fmt} string for each element. */

void dg_map_flatten(dg_item_map_t *sh);
  /* Replaces the control points of {sh->b} by the control points of a
    spline with the same corners an degrees. which coincides with a
    multiaffine (degree 1) spline. Also clears out the {sh->dev}
    vector. */ 

void dg_map_snap_half
  ( dg_item_map_t *sh, 
    dg_axis_t ax, 
    dg_BinaryDir dir, 
    double_vec_t dev,
    double tol
  );
  /* This procedure assumes that {sh} is the map of one half of some
    parent face {E}, whose deviation vector is {dev}; specifically,
    the half of {E} that lies in direction {dir} along axis {ax} of
    the domain. The procedure recomputes the nonflatness {sh->dev[f]}
    of every face {f} of the half-face {sh}. 
    
    Let {F} be the facet of {E} that lies in the direction {dir} along
    axis {ax}. If the face {f} is contained in {F}, the nonflatness is
    simply inherited from the parent cell, i.e. {sh->dev[f] = dev[f]}.
    Otherwise the procedure attempts to flatten {f}, if that can be
    done with error not exceeding {tol}. */
    
void dg_map_propagate_dev(double_vec_t *dev, box_face_index_t f, bz_DDim m); 
  /* Given the face deviation vector {dev} of some {m}-dimensional 
    cell or face {C}, and the index {f} of some face of {C},
    sets {dev.el[f1] = 0} for all facets {f1} of face {f}. */ 

/* TESTING */

dg_item_map_t *dg_map_make_test(bz_DDim d, bz_DDim n, bz_Degree g, double tol);
  /* Creates a map for the root cell of a {d}-dimensional
    finite dyadic grid, calling {dg_make_test_bezier}
    and {dg_map_from_bezier}. */

#endif