/* See somegrids.h */
/* Last edited on 2007-12-27 00:40:54 by stolfi */

#include "somegrids.h"


/* INTERNAL PROTOTYPES */

dg_tree_node_t *somegrid_adaptive(mdg_dim_t d, mdg_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(mdg_cell_index_t k, mdg_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(mdg_cell_index_t k, mdg_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(mdg_cell_index_t k, mdg_rank_t r, bz_patch_t *b, dg_tree_node_t *p)
      { 
        /* fprintf(stderr, "  omit k = %lld\n", k); */
        return FALSE;
      }

    bool_t split(mdg_cell_index_t k, mdg_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
          { mdg_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);
        double tol_eps = 0.001;
        radMax = dg_test_local_tol(n, ctr, tol_eps);
#warning eu... comeco
        /* tentando mudar a malha */
        
        auto double gaussiana(mdg_dim_t d, double x[]);
        double gaussiana(mdg_dim_t d, double x[])
        {
          double xc[d];/* = 0.500, yc = 0.500*M_SQRT1_2;*/
          //xc = 0.0, yc = 0.0; // so pq o dominio vai ser [-1,1]x[-1,1]
          double r = 0.2;
          /* Gaussiana centrada em {xc,yc} com raio medio {r}: */
          int i;
          double dx[d], dist = 0.0;
          for(i=0; i<d; i++)
          {
            xc[i] = 0.500*pow(2.0,-i/d);
            dx[i] = x[i] - xc[i];
            dist += dx[i]*dx[i];
          }

          if( dist< r*r)
          { return exp(-dist/r); }
          else
          { return 0.0; }
        }
        rad = gaussiana( d, ctr); /* valor pontual da gaussiana no centro da celula */
        Funcao *f = &gaussiana;
        
        double xmin[d], xmax[d], cvol=1.0;
        for (i=0; i<d; i++)
        {
          interval_t bi = bbox[i];
          xmin[i] = LO(bi);
          xmax[i] = HI(bi);
          cvol *= fabs(xmax[i]-xmin[i]);
        }
        radMax = integra( d, xmin, xmax, f)/cvol;
        /* fim malha */
        double eps = 0.001;
        mdg_rank_t tree_rank = 10;
        tol_eps = pow(2.0,-d*(tree_rank - (r+1)))*eps;
#warning eu... fim
        
        fprintf(stderr, "  rad = %7.4f radMax = %7.4f\n", rad, radMax);
//         return rad >= radMax;
        return fabs(rad - radMax)  >= tol_eps; //eu...
      }

//     return dg_make_grid(d, mdg_ROOT_CELL, 0, b, NULL, &omit, &split);
     return dg_make_grid ( d, 1, 0, b, NULL, &split);

  }

// dg_tree_node_t *somegrid_varied(mdg_dim_t d, mdg_rank_t splitDp, mdg_rank_t maxDp, bz_patch_t *b)
//   {
//     auto bool_t omit(mdg_cell_index_t k, mdg_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(mdg_cell_index_t k, mdg_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(mdg_cell_index_t k, mdg_rank_t r, bz_patch_t *b, dg_tree_node_t *p)
//       { 
//         /* fprintf(stderr, "  omit k = %lld\n", k); */
//         return FALSE;
//       }
// 
//     mdg_cell_index_t firstCell[DG_MAX_RANK];
//     mdg_rank_t r;
//     for (r = 0; r <= DG_MAX_RANK; r++) { firstCell[r] = DG_NO_CELL; }
//     
//     bool_t split(mdg_cell_index_t k, mdg_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. */
//         mdg_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);
//   }