/* Last edited on 2004-06-19 21:42:35 by stolfi */
/* Test of libdygrid routines. */

#include <stdio.h>
#include <math.h>
#include <dgbasic.h>
#include <dgtree.h>
#include <dgbezier.h>
#include <dgplot.h>
#include <dgmakegrid.h>
#include <js.h>
#include <pswr.h>
#include <pswr_color.h>
#include <pswr_iso.h>

/* Dimension of grid's domain: */
#define DDIM (2)  

/* Dimension of grid's range space {(X,Y,F(X,Y))}: */
#define MAX_RDIM (5)

/* Number of faces in a cell of dimension {DDIM}, namely {3^DDIM}: */
#define NF (9)

/* Flatness tolerance (mm): */
#define TOL_MM (0.5)

/* TRUE to debug {make_shape} instead of {dg_enum_faces}: */
#define DEBUG_SHAPE FALSE

/* INTERNAL PROTOTYPES */

int main(int argc, char **argv);
  
dg_Node *call_make_grid(dg_Dim d, dg_Rank maxDepth, bz_Patch *b);
  /* Creates a random finite dyadic grid of dimension {d} and
    maximum depth {maxDepth}, whose root cell has shape {b}. */
   
PSStream *init_plot
  ( char *name, 
    bool epsformat,
    char *paperSize,
    dg_Interval bbox[]
  );
  /* Opens a Postsript plotfile with given name, plus the 
    appropriate extension (".eps" or ".ps"), and initializes
    it for plotting.  The scales are set so that the 
    rectangle {bbox[0] × bbox[1]} fits snugly in the available 
    plotting area, with equal scales. */

dg_LocusProc plot2D;
  /* Plots a leaf item of a finte multigrid. */
     
bz_Patch *do_extract_face(bz_Patch *b, dg_FaceIndex f);
/* Extracts face {f} of the Bezier patch {b} and stores it into a newly
  allocated record. */

void finish_plot
  ( PSStream *ps,  
    bool epsformat
  );
  /* Finalizes and closes the Postscript plotfile {f}. */

void slow_plot(PSStream *ps, bz_Patch *b);
void slow_plot_vertex(PSStream *ps, bz_Patch *b, int steps);
void slow_plot_edge(PSStream *ps, bz_Patch *b, int steps);
void slow_plot_face(PSStream *ps, bz_Patch *b, int steps, double R, double G, double B);
  /* Slow and dumb plotting routines, for debugging. */

double tol(double x, double y);
  /* Maximum box radius allowed near {p}. */ 

double inside(double x, double y, double rad);
  /* TRUE if a box of radius {rad} centered at {p} intersects the 
    desired domain. */ 

double dsoft(double dx, double dy, double eps);
  /* Length of {(dx,dy)}, cooked to be at least {eps}. */

void test_plot(char *name, bz_Degree g, bz_RDim n, dg_Rank treeDp, dg_Rank plotDp);
  /* Tests of plotting routines, with a root patch of degree {g},
    a random tree of maximum depth {treeDp}, and maximum plotting split 
    depth {plotDp}. */

/* IMPLEMENTATIONS */

int main(int argc, char **argv)
  {   
    // test_plot("test", 1, 3, 1, 3);
    // test_plot("test", 1, 4, 5, 5);
    test_plot("test", 1, 3, 14, 14);
    // test_plot("test", 2, 3, 5, 5);
    // test_plot("test", 3, 3, 14, 12);
    fprintf(stderr, "OK so far...\n");
    return(0);
  }
  
static PSStream *ps;
static dg_Rank maxPlotDepth;

void plot2D(
    dg_Dim d,
    dg_Locus E,
    dg_Rank r, 
    bz_Patch *b,
    dg_NodeStar *na
  )
  { int m = b->m; /* Dimension of item. */
    int n = b->n; /* Dimension of range (coords plus func values). */
    fprintf(stderr, "+ plot2D(E = <%d:%lld> r = %d)\n", E.norm, E.cell, r);
    fprintf(stderr, "item shape\n");
    bz_print(stderr, b, "%6.2f");
    affirm(m <= 2, "bad domain dimension");
    double fMax = 1.0;
    double fStep = 0.2;
    double fSync = fStep/2;
    int kMin = pswr_inf_isoline(fSync, fStep, -fMax);
    int kMax = pswr_sup_isoline(fSync, fStep, +fMax);
    int N;
    double *R, *G, *B;
    pswr_make_color_table(
      fSync, fStep, kMin, kMax,
      0.000, 0.500, 1.000,
      1.000, 1.000, 1.000,
      1.000, 0.167, 0.000,
      &N, &R, &G, &B);
    switch(m)
    { case 0: 
        { double *c = &(b->c[0]);
          dg_plot_2D_vertex(ps, c, n);
        }
        break;
      case 1: 
        affirm(b->g[0] == 1, "bad patch degree");
        { double *c0 = &(b->c[0]);
          double *c1 = &(b->c[n]);
          dg_plot_2D_edge(ps, c0, c1, n); 
        }
        break;
      case 2: 
        affirm((b->g[0] == 1) && (b->g[1] == 1), "bad patch degree");
        { double *c00 = &(b->c[0]);
          double *c01 = &(b->c[n]);
          double *c10 = &(b->c[2*n]);
          double *c11 = &(b->c[3*n]);
          int subDp = maxPlotDepth - r + 1;
          if (subDp < 0) { subDp = 0; }
          if (subDp > 3) { subDp = 3; }
          dg_plot_2D_face
            ( ps, c00, c01, c10, c11, n, subDp, fSync, fStep, kMin, kMax, R,G,B );
        }
        break;
      default: affirm(FALSE, "bad item dimension"); 
    }
    fprintf(stderr, "- plot2D\n");
  }

void test_plot(char *name, bz_Degree g, bz_RDim n, dg_Rank treeDp, dg_Rank plotDp)
  {
    bool epsformat = FALSE;
    bz_Patch b = bz_make_test_patch(DDIM, n, g);
    dg_Node *grid_root = call_make_grid(DDIM, treeDp, &b);
    dg_Locus E0 = dg_locus(0, DG_ROOT_CELL);
    dg_NodeStar na = dg_NodeRef_vec_new(1);
    dg_Interval bbox[MAX_RDIM];

    affirm(b.m == 2, "wrong dimension");
    fprintf(stderr, "grid cells = %d\n", grid_root->nodes);
    fprintf(stderr, "root shape =\n");
    bz_print(stderr, &b, "%6.2f");
    /* bz_Patch_dev_print(stderr, sh->dev, b.m, "%.2f"); */
    /* fprintf(stderr, "\n"); */
    maxPlotDepth = plotDp;

    bz_compute_bbox(&b, bbox);
    ps = init_plot(name, epsformat, "a3", bbox);
    
    if (DEBUG_SHAPE)
      { slow_plot(ps, &b); }
    else
      { na.el[0] = grid_root;
        dg_enum_locus_leaves(DDIM, E0, 0, &b, &na, plotDp, plot2D);
      }
    finish_plot(ps, epsformat);
  }

PSStream *init_plot
  ( char *name, 
    bool epsformat,
    char *paperSize,
    dg_Interval bbox[]
  )
  { double hpt, vpt;
    if (epsformat)
      { hpt = 5.0*72.0; vpt = hpt * 4.0/3.0; }
    else
      { pswr_get_paper_dimensions(paperSize, &hpt, &vpt); }
    PSStream *ps = pswr_new_stream(epsformat, name, NULL, paperSize, hpt+6, vpt+8);
    pswr_set_page_layout(ps, 3.0, 4.0, TRUE, 0.0, 0.0, (epsformat ? 0 : 1), 1,1);
    double dx = 0.1*(bbox[0].end[BHI] - bbox[0].end[BLO])/2;
    double dy = 0.1*(bbox[1].end[BHI] - bbox[1].end[BLO])/2;
    pswr_new_picture
      ( ps, 
        bbox[0].end[BLO]-dx, bbox[0].end[BHI]+dx, 
        bbox[1].end[BLO]-dy, bbox[1].end[BHI]+dy
      );
    pswr_set_pen(ps, 0,0,0, 0.15, 0,0);
    if (! epsformat)
      { pswr_add_caption(ps, "Bezier-distorted grid", 0.0); }
    pswr_set_pen(ps, 1,0,0, 0.25, 0,0);
    return ps;
  }

void finish_plot
  ( PSStream *ps,  
    bool epsformat
  )
  { pswr_close_stream(ps); }

#define x_well_1 (2.5)
#define y_well_1 (4.5)

#define x_well_2 (4.0)
#define y_well_2 (5.0)

#define y_fault (2.5)
#define r_fault (0.5)
#define w_fault (1.5)

double tol(double x, double y)
  { double d1 = dsoft(x - x_well_1, y - y_well_1, 0.01);
    // double d1 = 100000.00;
    double d2 = dsoft(x - x_well_2, y - y_well_2, 0.01);
    // double d2 = 100000.00;
    double yfx = y_fault + r_fault*sin(w_fault*x);
    double d3 = dsoft(0.0, y - yfx, 0.01);
    // double d3 = 100000.00;
    double maxr = 0.25 * 3.0/(1.0/d1 + 1.0/d2 + 1.0/d3);
    return maxr;
  }

double inside(double x, double y, double rad)
  { double dctr = hypot(x - 0.5, y - 0.5);
    return dctr - rad < 0.4; 
  }

double dsoft(double dx, double dy, double eps)
  { return sqrt(dx*dx + dy*dy + eps*eps); }

dg_Node *call_make_grid(dg_Dim d, dg_Rank treeDp, bz_Patch *b)
  {
    dg_Interval bbox[MAX_RDIM];
    double ctr[MAX_RDIM];
    int n = b->n;

    auto bool omit(dg_CellIndex k, dg_Rank r, bz_Patch *b, dg_Node *p);
      /* Tells {dg_make_grid} when to discard a cell. Basically, when it
        does not intersect a specified circle. */ 

    auto bool split(dg_CellIndex k, dg_Rank r, bz_Patch *b, dg_Node *p);
      /* Tells {dg_make_grid} when to split a cell. Basically, when the variation in
        the coordinates and function values within the cell is larger
        than some tolerance, which depends on the cell's center. */ 

    bool omit(dg_CellIndex k, dg_Rank r, bz_Patch *b, dg_Node *p)
      { // int i; double r2 = 0.0;
        // if (r >= treeDp) return FALSE;
        // bz_compute_bbox(b, bbox);
        // for (i = 0; i < d; i++) 
        //   { dg_Interval bi = bbox[i];
        //     double ri = (bi.hi - bi.lo)/2;
        //     ctr[i] = (bi.lo + bi.hi)/2;
        //     r2 += ri*ri;
        //   }
        // return ! inside(ctr[0], ctr[1], sqrt(r2));
        fprintf(stderr, "  omit k = %lld\n", k);
        return FALSE;
      }

    bool split(dg_CellIndex k, dg_Rank r, bz_Patch *b, dg_Node *p)
      { int i; double r2 = 0.0;
        double rad, radMax;
        if (r >= treeDp) return FALSE;
        bz_compute_bbox(b, bbox);
        for (i = 0; i < n; i++) 
          { dg_Interval bi = bbox[i];
            double ri = (bi.end[BHI] - bi.end[BLO])/2; /* Half-width of interval. */
            ctr[i] = (bi.end[BLO] + bi.end[BHI])/2;
            r2 += ri*ri;
          }
        rad = sqrt(r2);
        radMax = tol(ctr[0], ctr[1]);
        fprintf(stderr, "  split k = %lld  rad = %7.4f radMax = %7.4f\n", k, rad, radMax);
        return rad >= radMax;
      }

    return dg_make_grid(d, DG_ROOT_CELL, 0, b, NULL, &omit, &split);
  }

void slow_plot(PSStream *ps, bz_Patch *b)
  { 
    int BSTEPS = 10;
    double R = 1.0, G = 0.9, B = 0.6;
    int f;
    fprintf(stderr, "+slow_plot dim = %d\n", b->m);
    /* Plot faces */
    for (f = 0; f < NF; f++)
      { if (dg_face_dimension(f, b->m) == 2)
        { bz_Patch *t = do_extract_face(b, f);
          slow_plot_face(ps, t, BSTEPS, R, G, B);
          free(t->c); free(t);
        }
      }
    /* Plot edges */
    pswr_set_pen(ps, 0.0, 0.0, 0.0,   0.15, 0.0, 0.0);
    for (f = 0; f < NF; f++)
      { if (dg_face_dimension(f, b->m) == 1)
        { bz_Patch *t = do_extract_face(b, f);
          slow_plot_edge(ps, t, BSTEPS);
          free(t->c); free(t);
        }
      }
    /* Plot vertices */
    pswr_set_pen(ps, 0.0, 0.0, 0.0,   0.15, 0.0, 0.0);
    for (f = 0; f < NF; f++)
      if (dg_face_dimension(f, b->m) == 0) 
        { bz_Patch *t = do_extract_face(b, f);
          slow_plot_vertex(ps, t, BSTEPS);
          free(t->c); free(t);
        }
    fprintf(stderr, "-slow_plot\n");
  }

void slow_plot_vertex(PSStream *ps, bz_Patch *b, int steps)
  {
    fprintf(stderr, "+slow_plot_vertex dim = %d\n", b->m);
    affirm(b->m == 0, "invalid vertex patch");
    pswr_set_fill_color(ps, 0.0,0.0,0.0);
    pswr_dot(ps, b->c[0], b->c[1], 0.25, TRUE, FALSE);
    fprintf(stderr, "-slow_plot_vertex\n");
  }

void slow_plot_edge(PSStream *ps, bz_Patch *b, int steps)
  {
    fprintf(stderr, "+slow_plot_edge dim = %d\n", b->m);
    affirm(b->m == 1, "invalid edge patch");
    double p[2], q[2];
    int i;
    fprintf(stderr, "shape =\n");
    bz_print(stderr, b, "%6.2f");
    for (i = 0; i <= steps; i++)
      { double ui = ((double)i)/((double)steps);
        bz_eval(b, &ui, p);
        if (i > 0) { pswr_segment(ps, q[0], q[1], p[0], p[1]); }
        q[0] = p[0]; q[1] = p[1];
      }
    fprintf(stderr, "-slow_plot_edge\n");
  }

void slow_plot_face(PSStream *ps, bz_Patch *b, int steps, double R, double G, double B)
  {
    int i, j;
    fprintf(stderr, "+slow_plot_face dim = %d\n", b->m);
    affirm(b->m == 2, "invalid face patch");
    fprintf(stderr, "shape =\n");
    bz_print(stderr, b, "%6.2f");
    double u[2], v[2], p[MAX_RDIM], q[MAX_RDIM], r[MAX_RDIM], s[MAX_RDIM];
    for (i = 1; i <= steps; i++)
      { for (j = 0; j <= steps; j++)
          { u[0] = ((double)i)/((double)steps);
            u[1] = ((double)j)/((double)steps);
            bz_eval(b, u, p);
            v[0] = ((double)i-1)/((double)steps);
            v[1] = ((double)j)/((double)steps);
            bz_eval(b, v, r);
            if (j > 0) 
              { double xm = (p[0]+q[0]+r[0]+s[0])/4.0;
                double ym = (p[1]+q[1]+r[1]+s[1])/4.0;
                pswr_set_fill_color(ps, R,G,B);
                pswr_triangle(ps, p[0], p[1], q[0], q[1], xm, ym, TRUE, FALSE);
                pswr_triangle(ps, p[0], p[1], r[0], r[1], xm, ym, TRUE, FALSE);
                pswr_triangle(ps, s[0], s[1], q[0], q[1], xm, ym, TRUE, FALSE);
                pswr_triangle(ps, s[0], s[1], r[0], r[1], xm, ym, TRUE, FALSE);
              }
            q[0] = p[0]; q[1] = p[1];
            s[0] = r[0]; s[1] = r[1];
          }
      }
    fprintf(stderr, "-slow_plot_face\n");
  }

bz_Patch *do_extract_face(bz_Patch *b, dg_FaceIndex f)
  {
    dg_Dim bm = b->m, tm = 0;
    dg_SignedDir loc[DDIM];
    dg_face_signature(f, bm, loc);
    bz_Patch *t = (bz_Patch *)notnull(malloc(sizeof(bz_Patch)), "no mem");
    bz_Degree tg[DDIM];
    int i;
    /* Extract degree sequence: */
    for(i = 0; i < bm; i++)
      { if (loc[i] == SMD) { tg[tm] = b->g[i]; tm++; } }
    (*t) = bz_Patch_new(tm, b->n, tg);
    bz_get_face(b, loc, t);
    return t;
  }