/* Finite elements for multidimensional dyadic spline spaces. */
/* Last edited on 2005-02-04 16:28:36 by stolfi */

#ifndef dgtent_H
#define dgtent_H

#include <dgbasic.h>
#include <dgbezier.h>
#include <dgpulse.h>
#include <dgbezpatch.h>
#include <vec.h>
#include <js.h>

/*
  STANDARD DYADIC TENTS
  
  For multidimensional grids, the finite-element bases that we
  consider here consist of /standard dyadic finite tents/, or /dyadic
  tents/ for short, which are characterized by the following
  parameters:

    * a /pulse kind vector/ {pkind[0..d-1]} 
    * a /continuity vector/ {c[0..d-1]}
    * a /degree vector/ {g[0..d-1]}
    * a /pulse index vector/ {pix \geq 0}
    * a /reference cell/ {C} 
    
  Any dyadic tent {t} is a polynomial spline, defined on level
  {r=rank(C)} of the unidimensional dyadic multigrid {G*}, with
  continuity {c[i]} and maximum degree {g[i]} along cordinate axis
  {i}, for each {i} in {0..d-1}. It is the product of {d} univariate
  dyadic pulses, each depending on one coordinate:
  
    { t(z) = PROD { b[i](z[i]) : i \in 0..d-1 } }
  
  Each pulse factor {b[i]} has pulse kind {pkind[i]}, continuity
  {c[i]}, degree {g[i]}, and index {pix[i]}.
  
  TENT SUPPORT
  
  It folows that the support of a standard dyadic tent is some
  fully-dimensional box {B} of {T^d}, which is the Cartesian product
  of the supports of its factors. The box {B} is the union of a
  compact block of cells and faces of rank {r = rank(C)}.
  
  By definition, the reference cell {C} of the tent is the cell in {B}
  that is incident to its the upper corner.  The number of cells {wd[i]}
  in {B} along axis {i} can be determined by 
  
    {wd[i] = dg_pulse_supp_count(pkind[i],c[i],g[i],pix[i],sz[i])}
    
  where {sz[i]} is the grid size of level {r} of {G*} along direction {i}.
  
  TENT FACTORS
  
  Recall that, at present, all standard dyadic pulses are supported on
  two consecutive intervals of some level of the uni-dimensional grid.
  It follows that a standard dyadic tent is supported on the star
  {K(E)} of a vertex, the lower corner of the reference cell {C}. 
  The vertex {E} is If
  the rank of {C} is {d} or more, the star {K(E)} comprises {2^d} cells
  and their separating faces.

  Let {g} be at least {2c+1}, and {C} be any cell in an infinite
  regular grid {L}. Associated with {C} we have {(g-c)^d} /tent
  functions/ which are {d}-dimensional splines of degree {g} and
  continuity {c}. These splines are supported on {C} and on the cells
  opposite to {C} across its inferior faces.
  
  Specifically, each tent function is identified by a tuple
  of indices {pix[0..d-1]} where each {pix_i} ranges in {0,..g-c-1}.
  The corresponding tent funtion is the product of univariate
  pulses (B-pulse or Hermite),
  
    { W^g_c[pix](z) = PROD { w^g_c[pix[i]](z[i]) : i in {0,..d-1} } }
    
  where {z[0..d-1]} are the coordinates of the point relative to
  the cell {C}.  
  
  TENT CENTER
  
  Let {t} be a dyadic tent function which is given by {W^g_c[pix]}
  relative to a dyadic cell {C}. Note that the support of the basis
  element {w^g_c[pix[i]]} is {[0_+1]} if {i > c}, and {[-1_+1]} if {i
  \leq c}. It follows that the support {D(t)} of {t} is the star
  {K(E)}, where {E} is the inferior face of {C} whose normal axes are
  all axes {i} such that {pix[i] \leq c}. The locus {E} is the /center/
  of the tent.
  
  TENTS WITH FOLDED-OVER SUPPORTS
  
  The definitions above hold also for tents on periodic grids.
  However, if a cell {C} spans the whole extent of the root cell along
  some axis {i}, the coordinate {z[i]} relative to {C} is
  indeterminate up to an integer shift; and, for the tent value
  {W^g_c[e](z)}, it makes a difference whether we use {z[i]} (which
  lies in {[0_+1]}) or {z[i]-1} (which lies in {[-1_0]}.
  
  In such cases, we must evaluate the tent function {W^g_c[e](z)} with
  both {z[i]} and {z[i]-1}, and add the results together.
  
  TENT INDICES
  
  The tent index {tix} is a single integer that conveniently packs
  these {d} pulse indices into a single integer, in a way that depends
  on the tent kind.
  
  For "H" and "B" tents, the index tuple {pix[0..d-1]} can be encoded into a single 
  /tent index/ { tix = SUM { e[i] (g-c)^i : i in {0,..d-1} } }. */
  
typedef int32 dg_TentIndex;
  /* Index of a tent among the {(g-c)^d} {d}-dimensional
    tents with given degree, continuity, and reference cell. */
  
/*
  TENT KINDS

  This module presently implements three kinds of finite tents:
  
    /H-tents/ ({kind = DG_TK_H}) are useful for Hermite-style
    interpolation of data given at grid corners, given all
    derivatives to a maximum order {c[i]} on each coordinate {i}.
    Each pulse index {pix[i]} ranges independently over {0..c}.
    
    /B-tents/ ({kind = DG_TK_B}) have the property that they are a
    partition of unity (non-negative and add to 1 everywhere).
    Each pulse index {pix[i]} ranges independently over {0..c}.
  
    /D-tents/ ({kind = DG_TK_D}) are a subset of the H-tents 
    that allow interpolation of data at cell corners, given
    all derivatives to some total order {c}.
    Each pulse index {pix[i]} ranges over {0..c},
    but their total sum is at most {c}.
  
  All three kinds are equivalent for {c=(0,..0), g=(1,..1)}. */
  
typedef enum { DG_TK_B, DG_TK_H, DG_TK_D } dg_TentKind;

typedef struct dg_Tent  /* A generalized tent function. */ 
  { dg_CellIndex cell;  /* Reference cell of the tent. */
    dg_TentIndex tix;   /* Packed tent index tuple. */
  } dg_Tent;
  /* 
    Identifies the tent whose tent index is {tix} and whose reference
    cell has cell index is {cell}. The other parameters (grid
    dimension, tent kind, continuity order, and degree) should be
    provided separately.  */

dg_TentIndex dg_tent_index
  ( dg_Dim d, 
    dg_TentKind kind,
    dg_Cont c, 
    dg_Degree g, 
    dg_PulseIndex pix[]
  );
  /* Computes the packed tent index {tix} from the given tuple of pulse
    indices {pix[0..d-1]}, assuming that the tent has dimension {d}, 
    tent kind {kind}, continuity {c}, degree {g}. */

dg_TentIndex dg_tent_index_max
  ( dg_Dim d, 
    dg_TentKind kind,
    dg_Cont c, 
    dg_Degree g
  );
  /* The maximum value of a packed tent index {tix} for tents
    of dimension {d}, tent kind {kind}, continuity {c}, and degree {g}. */

void dg_tent_get_pulse_indices
  ( dg_Dim d, 
    dg_TentKind kind,
    dg_Cont c, 
    dg_Degree g, 
    dg_TentIndex tix, 
    dg_PulseIndex pix[]
  );
  /* Stores in {pix[0..d-1]} the index tuple corresponding to the 
    packed tent index {tix}, assuming that the tent has dimension
    {d}, and continuity {c}, degree {g}. */

dg_Locus dg_tent_center
  ( dg_Dim d, 
    dg_TentKind kind,
    dg_Cont c,
    dg_Degree g, 
    dg_CellIndex k,
    dg_PulseIndex e[]
  );
  /* Returns the locus at the center of the support of the tent with
    dimension {d}, continuity {c}, degree {g}, and index tuple
    {e[0..d-1]}, associated to the cell with index {k}. */

double dg_tent_eval_cell_relative
  ( dg_Dim d, 
    dg_TentKind kind,
    dg_Cont c, 
    dg_Degree g, 
    dg_TentIndex tix, 
    double  *z,
    dg_GridSize sz[]
  );
  /* Evaluates the tent function of dimension {d}, continuity {c},
    degree {g} and packed index tuple {tix} at a point {z[0..d-1]} in
    the star {K(t)}.
    
    Each coordinate {z[i]} must be relative to the reference cell
    {t.cell}, and must be reduced modulo the corresponding grid size
    {sz[i]}, in such a way that, for a point in the tent's support
    {D(t)}, {z[i]} lies in the range {[-1_+1]}. */

double dg_tent_eval_root_relative
  ( dg_Dim d, 
    dg_TentKind kind,
    dg_Cont c, 
    dg_Degree g, 
    dg_Tent t, 
    double  *x
  );
  /* Evaluates the tent function {t}, assumed to have dimension {d},
    continuity {c}, and degree {g}, at a point whose coordinates are
    {x[0..d-1]} relative to the root cell. This procedure
    properly accounts for any wrap-around of the support {K(t)}. */

void dg_tent_eval_diff_root_relative
  ( dg_Dim d, 
    dg_TentKind kind,
    dg_Cont c, 
    dg_Degree g, 
    dg_Tent t, 
    double *x,
    double *f
  );
  /* Stores in {f[0...(c+1)^d]} the value and all derivatives up to
    max order {c} of the tent function {t} at the point whose
    coordinates are {x[0..d-1]} relative to the root cell. Assumes
    that {t} has dimension {d}, continuity {c}, and degree {g}.
    
    The elements of {f} are sorted by total order and then
    lexicographically by axes. That is, {f[0]} is the tent's value;
    {f[1..d]} are the first derivatives along each axis;
    {f[d+1..d+choose(d+1,d-1)]} are the second derivatives, and so on.
    
    Note that the derivatives are computed in root-relative
    coordinates, namely assuming that the root cell has unit size
    along all axes.  This procedure
    properly accounts for any wrap-around of the support {K(t)}. */
/* 
  BÉZIER REPRESENTATION */

void dg_mother_tent_to_bezier
  ( dg_Dim d, 
    dg_TentKind kind,
    dg_Cont c, 
    dg_Degree g,
    dg_TentIndex tix, 
    double **bp
  );
  /* 
    Splits the {d}-dimensional mother tent {t} of the given {kind},
    continuity order {c}, degree {d}, and tent index {tix}
    into its {NP} constituent polynomial patches {tp[0..NP-1]},
    computes the array of {NB} Bézier coeffs of each patch {tp[k]},
    and adds that array to the given array {bp[k][0..NB-1]}.
    If {bp[k]} is NULL, an appropriate coeff array is automatically
    allocated and initialized to zeros. */

/* PRINTOUT: */

void dg_tent_kind_print(FILE *wr, dg_TentKind kind);
  /* Prints the tent {kind} ("H" for {DG_TK_H}, "B" for {DG_TK_B}, etc.) */
  
void dg_tent_print(FILE *wr, dg_Dim d, dg_TentKind kind, dg_Cont c, dg_Degree g, dg_Tent t);
  /* Prints the tent {t} in the format 
    "{KIND}:{RANK}:({POS[0]},..{POS[d-1]})[{IX[0]},..{IX[d-1]}" */

/* Vectors of {dg_Tent}: */

typedef struct dg_Tent_vec_t { nat nel; dg_Tent *el; } dg_Tent_vec_t;
  /* A list of unit tent functions. */

dg_Tent_vec_t dg_Tent_vec_new(nat nel);

#define dg_Tent_vec_expand(nv,index) \
  vec_expand(vec_cast_ref(nv), index, sizeof(dg_Tent))
  /* Makes sure that the vector {nv} has element {nv.el[index]}.
    Reallocates {nv.el} if necessary, to have at least {index+1}
    elements -- but possibly more. */

#define dg_Tent_vec_trim(nv,nel) \
  vec_trim(vec_cast_ref(nv), nel, sizeof(dg_Tent))
  /* Makes sure that the vector {nv} has exactly {nel} elements.
    Reallocates {nv.el} if necessary. */

#endif