/* See dgitemshape.h */
/* Last edited on 2005-02-21 11:28:21 by stolfi */

#include <dg_grid.h>
#include <dg_tree.h>
#include <bz_basic.h>
#include <dgitemshape.h>
#include <affirm.h>
#include <vec.h>
#include <math.h>

/* INTERNAL PROTOTYPES */

void dg_shape_dev_check(double_vec_t dev, box_face_index_t f, bz_DDim m);
  /* Checks whether the non-flatness vector {dev} satisfies the invariant,
    namely that the facets of face {f} must be at least as flat as {f}. */

/* IMPLEMENTATIONS */

dg_ItemShape dg_ItemShape_new
  ( bz_DDim m,
    bz_RDim n,
    const bz_Degree *g
  )
  { dg_ItemShape sh;
    sh.b = bz_Patch_new(m, n, g);
    sh.dev = double_vec_new(ipow(3,m));
    return sh;
  }

dg_ItemShape dg_ItemShape_facet_new
  ( bz_DDim m,
    bz_RDim n,
    const bz_Degree *g,
    dg_axis_t ax
  )
  { dg_ItemShape sh;
    bz_Degrees gF;
    int iC, iF;
    for (iC = 0, iF = 0; iC < m; iC++)
      { if (iC != ax) { gF[iF] = g[iC]; iF++; } }
    sh.b = bz_Patch_new(m-1, n, gF);
    sh.dev = double_vec_new(ipow(3,m-1));
    return sh;
  }
  
void dg_ItemShape_free(dg_ItemShape sh)
  { bz_Patch_free(sh.b);
    free(sh.dev.e); 
    free(sh.ta.e);
  }

dg_ItemShape *dg_shape_for_cell
  ( bz_Patch *b,
    dg_rank_t r,
    dg_cell_index_t k, 
    dg_node_t *t,
    double tol
  )
  { int d = b->m;
    dg_ItemShape *sh = (dg_ItemShape *)malloc(sizeof(dg_ItemShape));

    fprintf(stderr, "+ dg_shape_from_bezier\n");

    sh->E = dg_item(k,0,r);
    sh->b = (*b);
    sh->dev = dg_shape_snap(b);
    sh->ta = dg_node_ref_vec_new(1); sh->ta.e[0] = t;
    
    fprintf(stderr, "- dg_shape_from_bezier\n");
    return sh;
  }

double_vec_t dg_shape_snap(bz_Patch *b)
  { int d = b->m;
    int nFaces = ipow(3, d);
    int i, f;
    double_vec_t dev = double_vec_new(nFaces); 
    
    /* Flatten near-flat faces, starting from lower-dimensional ones: */
    fprintf(stderr, "original patch b =\n");
    bz_print(stderr, b, "%6.2f");
    fprintf(stderr, "\n");

    for (f = nFaces-1; f >= 0; f--)
      { dev.e[f] = bz_try_flatten_face(b, f, tol);
        if (dev.e[f] == 0.0) { bz_print(stderr, b, "%6.2f"); }
      }
    if (
    dg_shape_print_dev(stderr, sh->dev, d, "%.2f");
    fprintf(stderr, "\n");
      sh->dev = dev;
    }

void dg_shape_split
  ( dg_ItemShape *sh,
    dg_axis_t ax,
    double tol, 
    dg_ItemShape *shLO,
    dg_ItemShape *shHI,
    dg_ItemShape *shMD
  )
  { bz_DDim m = sh->b.m;
    /* Allocate shape faces: */
    (*shMD) = dg_ItemShape_facet_new(m, sh->b.n, sh->b.g, ax);
    (*shLO) = dg_ItemShape_new(m, sh->b.n, sh->b.g);
    (*shHI) = dg_ItemShape_new(m, sh->b.n, sh->b.g);
    
    box_face_index_t nFaces = ipow(3, m); /* Number of sub-faces of shape. */
    
    /* Extract shapes of the two halves {SLO, SHI} of face {S}: */
    bz_split(&(sh->b), ax, 0.5, &(shLO->b), &(shHI->b));
    
    /* Copy global dimension and axes: */
    shLO->d = shHI->d = sh->d;
    for (i = 0; i < sh->d; i++) { shLO->gax[i] = shHI->gax[i] = sh->gax[i]; }
    
    /* Get the adjacent cells and nodes for the two half-faces: */
    shLO->E = dg_item_split...
    shHI->E = dg_item_split...;
    for (i = 0; i < nCells; i++)
      { int kai = sh->ka.e[i]; dg_node_t *tai = sh->ta.e[i];
        if (kai == dg_NULL_INDEX)
          { shLO->ka.e[i] = shHI->ka.e[i] = dg_NULL_INDEX; }
        else
          { shLO->ka.e[i] = 2*kai; shHI->ka.e[i] = 2*kai+1; }
        if ((tai == NULL) || (kai == dg_NULL_INDEX))
          { shLO->ta.e[i] = shHI->ta.e[i] = NULL; }
        else if (tai->index < kai)
          { shLO->ta.e[i] = shHI->ta.e[i] = tai; }
        else
          { /* If {ta[i]} is a leaf, keep it: */
            dg_node_t *taiLO = tai->ch[BLO];
            dg_node_t *taiHI = tai->ch[BHI];
            shLO->ta.e[i] = (taiLO == NULL ? tai : taiLO); 
            shHI->ta.e[i] = (taiHI == NULL ? tai : taiHI);
          }
      }

    /* Inherit and/or recompute non-flatness vectors: */
    dg_ItemShape_normalize_half(shLO, ax, BLO, sh->dev, tol);
    dg_ItemShape_normalize_half(shHI, ax, BHI, sh->dev, tol);
    
    /* Extract shape of separating face {SMD}: */
    dg_ItemShape_get_facet(shHI, ax, BLO, shMD);

    /* Checking for consistent split/flattening: */
    { dg_ItemShape shMD1 = dg_ItemShape_facet_new(m, sh->b.n, sh->b.g, ax);
      int h;
      dg_ItemShape_get_facet(shLO, ax, BHI, &shMD1);
      affirm(bz_dist(&(shMD->b), &(SHMD1.b)) < 1.0e-6*tol, "split with crack");
      for (h = 0; h < nFaces/3; h++) 
        { double devh = shMD->dev.e[h];
          double devh1 = shMD1.dev.e[h];
          affirm(fabs(devh - devh1) < 1.0e-6*tol, "inconsistent flatness");
        }
      dg_ItemShape_free(shMD1);
    }
  }

void dg_shape_get_facet
  ( dg_ItemShape *sh, 
    dg_axis_t ax, 
    dg_BinaryDir dir, 
    dg_ItemShape *shF 
  )
  {
    int m = sh->b.m;    /* Dimension of face. */
    int p = sh->d - m;  /* Number of global axes perpendicular to face. */
    int h, nFaces = sh->dev.ne;
    box_signed_dir_t sdir = 2*dir - 1;
    
    shF->b = bz_...
    bz_get_facet(&(sh->b), ax, dir, &(shF->b));
    
    shF->r = sh->r;
    
    shF->d = sh->d;
    
    /* Axis {splax} is perpendicular to facet, must permute {shF->ax}. */
    { int splj = p + ax, i; 
      for (i = 0; i < d; i++) { shF->gax[i] = sh->gax[i]; }
      { dg_axis_t tax = shF->gax[p]; shF->gax[p] = shF->gax[splj]; shF->gax[splj] = tax; }
    }
    
    shF->ka = dg_cell_index_vec_new(2*nCells);
    shF->ta = dg_cell_index_vec_new(2*nCells);
    
    
    /* Inherit the non-flatness of faces of {shF} from those of {sh}: */
    for (h = nFaces/3 - 1; h >= 0; h++) 
      { double devh;
        int f = dg_parent_face(h, m, ax, sdir);
        shF->dev.e[h] = sh->dev.e[f];
        /* Consistency check: */
        dg_shape_dev_check(shMD->dev, h, m-1); }
      }

void dg_shape_snap_half
  ( dg_ItemShape *sh, 
    dg_axis_t ax, 
    dg_BinaryDir dir, 
    double_vec_t dev,
    double tol
  )
  { int m = sh->b.m;
    int nFaces = dev.ne;
    box_signed_dir_t sdir = 2*dir - 1;
    int f;
    affirm(sh->dev.ne == dev.ne, "inconsistent face counts");
    for (f = 0; f < nFaces; f++) { sh->dev.e[f] = 1.0e200; }
    for (f = nFaces-1; f >= 0; f--) 
      { /* Inherit or compute the non-flatness of sub-face {f}: */
        box_signed_dir_t fsign = box_face_position(f, m, ax);
        double devf;
        if (fsign == sdir) 
          { /* {f} lies in the intact facet of the parent, inhert: */
            devf = dev.e[f];
          }
        else if (dev.e[f] == 0.0)
          { /* {f} is part of a flat face of {F}, it is flat too: */
            devf = 0.0;
          }
        else
          { devf = bz_try_flatten_face(&(sh->b), f, tol); }
        sh->dev.e[f] = devf;
      }
  }
  
void dg_shape_dev_check(double_vec_t dev, box_face_index_t f, bz_DDim m)
  { int f1 = f, df = 1, i;
    double devf = dev.e[f];
    for (i = m-1; i >= 0; i--)
      { if ((f1 % 3) == 0)
          { box_face_index_t fHI = f + df, fLO = fHI + df;
            affirm (fHI < dev.ne, "face indexing bug");
            assert (dev.e[fHI] <= devf;
            assert (dev.e[fLO] <= devf;
          }
        f1 /= 3; df *= 3;
      }
  }

void dg_shape_flatten(dg_ItemShape *sh)
  { /* Data for Bézier patch {b}: */
    bz_Patch *b = &(sh->b);
    int bm = b->m;
    bz_Degree *bg = &(b->g[0]);
    bool_t trivial = TRUE;
    { int ib;
      for (ib = 0; ib < bm; ib++)
        { if (bg[ib] != 1) { trivial = FALSE; } }
    }
    if (! trivial)
      { 
        bz_Patch a = bz_Patch_uniform_new(bm, b->n, 1);
        bz_multiaffine_approx(b, &a);
        bz_raise_degree(&a, b);
        free(a.c);
      }
    { int nFaces = ipow(3, bm), f;
      for (f = 0; f < nFaces; f++) { sh->dev.e[f] = 0.0; } 
    }
  }

void dg_shape_print_dev(FILE *wr, double_vec_t dev, bz_DDim m, char *fmt)
    { int nFaces = dev.ne, f;
      for (f = 0; f < nFaces; f++) 
        { if (f != 0) { fprintf(wr, " "); } 
          box_face_print(wr, f, m);
          fprintf(wr, ": ");
          fprintf(wr, fmt, dev.e[f]);
        } 
    }


dg_ItemShape *dg_shape_make_test(bz_DDim d, bz_DDim n, bz_Degree g, double tol)
  { 
    bz_Patch b = dg_make_test_bezier(d, n, g);
    return dg_shape_from_bezier(b, tol);
  }