/* See dgplotextra.h */
/* Last edited on 2004-08-17 23:59:46 by stolfi */

#include <dgplot.h>
#include <dgplotextra.h>
#include <dgtree.h>
#include <dgbezier.h>
#include <dgfaceshape.h>
#include <js.h>
#include <vec.h>
#include <ia.h>
#include <math.h>

/* IMPLEMENTATIONS */

void dg_enum_flattened
  ( dg_Dim d,            /* Dimension of grid. */
    dg_FaceShape *sh,    /* Shape of face {F}. */
    dg_Axis ax[],        /* Axes which are perpendicular and parallel to {F}. */
    dg_Rank r,           /* Rank of some cell of which {F} is face. */
    dg_Index_vec_t K,    /* Cells of rank {r} that are incident to {F}. */
    dg_NodeRef_vec_t N,  /* Nodes containing those cells. */
    double tol,          /* Flattening tolerance. */
    dg_Rank rMax,        /* Never split below this rank. */
    dg_FaceProc process  /* Client face-processing procedure. */  
  )
  { dg_Dim m = sh->b.m;        /* Dimension of domain of face {F}. */
    dg_Dim s = d-m;            /* Co-dimension of domain of face {F}. */
    int nFaces = ipow(3, m);   /* Number of faces of {F}. */
    int nCells = ipow(2, s);   /* Number of adjacent cells in {*k,*n}. */
    bool split;

    
    fprintf(stderr, "+ dg_enum_flattened(dim = %d r = %d)\n", sh->b.m, r);
    
    /* Check whether splitting-and-recursion is needed. */ 
    /* The face must be split if it is not flat, or if any adjacent cell
      has a split along an axis parallel to the face: */
    if (m == 0)
      { /* A vertex never needs splitting. */
        split = FALSE;
      }
    else if (r >= rMax)
      { /* Too deep to split; force the entire face to be flat: */
        dg_FaceShape_flatten(sh);
        split = FALSE;
      }
    else
      { int i;
        split = (sh->dev.el[0] > 0.0); 
        /* Could be improved --- should skip levels with parallel splits. */
        for (i = 0; (!split) && (i < nCells); i++)
          { split |= ((N.el[i] != NULL) && (N.el[i]->ch[BLO] != NULL)); }
      }

    if (! split) 
      { /* No need to split, process it: */
        dg_Bezier a = dg_Bezier_uniform_new(sh->b.m, sh->b.n, 1);
        dg_Bezier_multiaffine_approx(&(sh->b), &a);
        process(d, &a, ax, r, K, N);
        free(a.c);
      }
    else 
      { /* Split the adjacent cells, recurse: */
        dg_Axis splax = r % d;  /* Split axis. */
        int splj;               /* Position of split axis in {ax} list. */

        /* Find position {splj} of {splax} in the axis sequence: */
        for (splj = 0; (splj < d) && (ax[splj] != splax); splj++) { }
        affirm(splj < d, "inconsistent ax[]");
        fprintf(stderr, "split on axis ax[%d] = %d\n", splj, splax);

        if (splj < s)
          { /* The split axis is perpendicular to the face. */
            int smask = ipow(2, splj);
            dg_Index_vec_t Kch = dg_Index_vec_new(nCells);
            dg_NodeRef_vec_t Nch = dg_NodeRef_vec_new(nCells);
            int i;
            /* Split each adj cell, discard halves not incident to face: */
            for (i = 0; i < nCells; i++)
              { int Ki = K.el[i]; dg_Node *Ni = N.el[i];
                /* Decide which child of {K[i]} is to be kept: */
                int keep = ((i & smask) == 0 ? 1 : 0);
                /* Keep the child of {K[i]} that is adjacent to {F}: */
                if (Ki == DG_NULL_INDEX) 
                  { Kch.el[i] = DG_NULL_INDEX; }
                else
                  { Kch.el[i] = 2*Ki + keep; }
                /* If the node {N[i]} is level with {K[i]}, split it too: */
                if ((Ni == NULL) || (Ki == DG_NULL_INDEX))
                  { Nch.el[i] = NULL; }
                else if (Ni->index < Ki)
                  { Nch.el[i] = Ni; }
                else
                  { /* If {N[i]} is a leaf, keep it: */
                    dg_Node *Nchi = Ni->ch[keep];
                    Nch.el[i] = (Nchi == NULL ? Ni : Nchi);
                  }
              }
            /* Recurse on same face, with same shape. */
            dg_enum_flattened(d, sh, ax, r+1, Kch, Nch, tol, rMax, process);
            free(Kch.el); free(Nch.el);
          }
        else
          { /* The split axis is parallel to the face. */
            dg_FaceShape shMD, shLO, shHI;
            /* Split face and each adj cell: */
            /* Data for subfaces: */
            ?? dg_Index_vec_t Kch = dg_Index_vec_new(2*nCells);
            ?? dg_NodeRef_vec_t Nch = dg_NodeRef_vec_new(2*nCells);
            ?? /* Data for splitting facet: */
            ?? dg_FaceShape shMD = dg_FaceShape_facet_new(m, sh->b.n, sh->b.g, splj-s);
            ?? dg_Index_vec_t KMD = (dg_Index_vec_t){ 2*nCells, &(Kch.el[0]) };
            ?? dg_NodeRef_vec_t NMD = (dg_NodeRef_vec_t){ 2*nCells, &(Nch.el[0]) };
            ?? /* Data for low half-face: */
            ?? dg_FaceShape shLO = dg_FaceShape_new(m, sh->b.n, sh->b.g);
            ?? dg_Index_vec_t KLO = (dg_Index_vec_t){ nCells, &(Kch.el[0]) };
            dg_NodeRef_vec_t NLO = (dg_NodeRef_vec_t){ nCells, &(Nch.el[0]) };
            ?? /* Data for high half-face: */
            ?? dg_FaceShape shHI = dg_FaceShape_new(m, sh->b.n, sh->b.g);
            ?? dg_Index_vec_t KHI = (dg_Index_vec_t){ nCells, &(Kch.el[nCells]) };
            ?? dg_NodeRef_vec_t NHI = (dg_NodeRef_vec_t){ nCells, &(Nch.el[nCells]) };
              
            int i;

            /* Split the face shape: */
            dg_FaceShape_split(sh, splj - s, tol, &shLO, &shHI, &shMD);

            /* Recurse for each half of face {F}: */
            dg_enum_flattened(&shLO, tol, rMax, process);
            dg_enum_flattened(&shHI, tol, rMax, process);

            /* Recurse for the separating face {S} of {F}. */
            dg_enum_flattened(&shMD, tol, rMax, process);

            dg_FaceShape_free(* shMD);
            dg_FaceShape_free(* shLO);
            dg_FaceShape_free(* shHI);
            free(Kch.el); 
            free(Nch.el);
          }
      }

    fprintf(stderr, "- dg_enum_flattened(dim = %d r = %d)\n", sh->b.m, r);
    
  }