/* See somegrids.h */
/* Last edited on 2011-09-15 17:59:51 by stolfilocal */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

// #include <dg_grid.h>
// #include <dg_tree.h>
// #include <bz_patch.h>
// #include <dg_makegrid.h>
// #include <dg_test_tools.h>

#include <affirm.h>
#include <bool.h>

#include <somegrids.h>
#include <figtools.h>

/* INTERNAL PROTOTYPES */

dg_tree_node_t *somegrid_adaptive(dg_dim_t d, dg_rank_t maxDp, bz_patch_t *b)
  {
    interval_t bbox[MAX_RDIM]; /* BBox of the cell or Bezier patch {b}. */
    double ctr[MAX_RDIM];
    int n = (b == NULL ? 2 : b->n);

    auto bool_t omit(dg_cell_index_t k, dg_rank_t r, bz_patch_t *b, dg_tree_node_t *p);
      /* Tells {dg_make_grid} when to discard a cell. */ 

    auto bool_t split(dg_cell_index_t k, dg_rank_t r, bz_patch_t *b, dg_tree_node_t *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_t omit(dg_cell_index_t k, dg_rank_t r, bz_patch_t *b, dg_tree_node_t *p)
      { 
        /* fprintf(stderr, "  omit k = %lld\n", k); */
        return FALSE;
      }

    bool_t split(dg_cell_index_t k, dg_rank_t r, bz_patch_t *b, dg_tree_node_t *p)
      { int i; double r2 = 0.0;
        double rad, radMax;
        if (r >= maxDp) return FALSE;
        if (b != NULL)
          { bz_patch_compute_bbox(b, bbox); }
        else
          { dg_cell_box_canonical(d, k, bbox); } 
        fprintf(stderr, "  split k = %lld = ", k);
        print_2d_cell(stderr, bbox);
        for (i = 0; i < n; i++) 
          { interval_t bi = bbox[i];
            double ri = (HI(bi) - LO(bi))/2; /* Half-width of interval. */
            ctr[i] = (LO(bi) + HI(bi))/2;
            r2 += ri*ri;
          }
        rad = sqrt(r2);
        radMax = dg_test_local_tol(n, ctr, tol_eps);
        fprintf(stderr, "  rad = %7.4f radMax = %7.4f\n", rad, radMax);
        return rad >= radMax;
      }

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

dg_tree_node_t *somegrid_varied(dg_dim_t d, dg_rank_t splitDp, dg_rank_t maxDp, bz_patch_t *b)
  {
    auto bool_t omit(dg_cell_index_t k, dg_rank_t r, bz_patch_t *b, dg_tree_node_t *p);
      /* Tells {dg_make_grid} when to discard a cell. */ 

    auto bool_t split(dg_cell_index_t k, dg_rank_t r, bz_patch_t *b, dg_tree_node_t *p);
      /* Tells {dg_make_grid} when to split a cell. Basically, 
        a random function of the parent so that about 3/4 of the cells
        split at each level. */ 

    bool_t omit(dg_cell_index_t k, dg_rank_t r, bz_patch_t *b, dg_tree_node_t *p)
      { 
        /* fprintf(stderr, "  omit k = %lld\n", k); */
        return FALSE;
      }

    dg_cell_index_t firstCell[DG_MAX_RANK];
    dg_rank_t r;
    for (r = 0; r <= DG_MAX_RANK; r++) { firstCell[r] = DG_NO_CELL; }
    
    bool_t split(dg_cell_index_t k, dg_rank_t r, bz_patch_t *b, dg_tree_node_t *p)
      { /* Always split cells up to {splitDp}: */
        if (r < splitDp) { return TRUE; }
        /* Stop at depth {maxDp}: */
        if (r >= maxDp) { return FALSE; }
        
        if (firstCell[r] == DG_NO_CELL)
          { /* First cell in this level, split it: */
            firstCell[r] = k; return TRUE;
          }
        else if (firstCell[r]+1 == k)
          { /* Second cell, split it only in level {splitDp+1}: */
            return (r == splitDp + 1);
          }
        /* Split randomly 6/7 of the cells. */
        dg_cell_index_t parent = (k >> 1);
        double burp = ((double)parent) * ((sqrt(5) + 1.0)/2.0);
        burp -= floor(burp);
        dg_cell_index_t scramble = (int)floor(8*burp);
        return ((k & 7L) != (scramble & 7L));
      }

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