#include "grid/octree.h"
#include "embed.h"
#include "poisson.h"
#include "run.h"
#include "colonization.h"
#include "../ABL/perlin.h"
#include "../tracer-particles.h"
#include "view.h"
#include "../scatter2.h"

scalar pot[], laplace[]; 
Particles P;

pot[front] = dirichlet (1);
pot[back] = dirichlet (0);
pot[embed] = neumann (0.);
vector u[];
u.n[front] = neumann(0);
u.n[back] = neumann(0);
int ja = 1000; double R1 = 3, H = 1.5, Hg = 1;

double sy = 1.4;
int main () {
  L0 = 10;
  X0 = Z0 = -L0/2;
  N = 64;
  run();
}

face vector fp[];
scalar cp[];
event init (t = 0) {
  //Seed atraction points in a prolate semi spheriod
  srand (0);
  coord atr[ja];
  int j = 0;
  while (j < ja) {
    double R = R1*pow(fabs(noise()), 1./3.1);
    coord vec = {noise(), noise(), noise()};
    normalize (&vec);
    if (vec.y > 0) {
      foreach_dimension()
	atr[j].x = R*vec.x;
      atr[j].y *= sy;
      atr[j].y += H + Hg;
      j++;
    }
  }

  //Get tree nodes
  int tn = 1;
  Tnode * nodes = calloc (tn, sizeof(Tnode));
  nodes[0].x = 0;
  nodes[0].y = Hg;
  nodes[0].z = 0;
  tn = colonize (atr, ja, &nodes, 1, 2, 0.2, 0.3);

 // Convert nodes to branches
  Branch * trees = calloc (499, sizeof(Branch)); //?
  printf ("# nodes: %d\n", tn);
  nb = nodestotree (nodes, tn, &trees);
  printf ("# branches: %d\n", nb);
  fflush (stdout);

  //Convert branches to volume fractions
  int it = 0;
  do {
    tree_interface (trees, cs, fs);
    boundary({cs});
    printf ("# Tree Reconstruction Iteration %d\n", ++it);
  } while (adapt_wavelet ({cs}, (double[]){0.1}, 9, 4).nf > (grid->tn/100));
  tree_interface (trees, cs, fs);
  boundary({cs, fs});
  fractions_cleanup (cs, fs);

  //Obtain surface volume and face fractions
  cp.prolongation = cs.prolongation;
  vertex scalar phi[];
  int nx = 5;
  int ny = 3;
  init_perlin (nx, ny);
  it = 0;
  do {
    foreach_vertex()
      phi[] =  y - Hg - (perlin (x, z, nx, ny));
    boundary ({phi});
    fractions (phi, cp, fp);
    boundary ({cp});
    printf ("# Surface Reconstruction Iteration %d\n", ++it);
  }  while (adapt_wavelet ({cs, cp}, (double[]){0.1, 0.01}, 8, 4).nf > (grid->tn/1000));
  foreach_vertex()
    phi[] =  y - Hg - (perlin (x, z, nx, ny));
  boundary ({phi});
  fractions (phi, cp, fp);
 
  //Combine tree and surface
  foreach()
    cs[] = min (cs[], cp[]);
  foreach_face()
    fs.x[] = min(fs.x[], fp.x[]);
  boundary ({cs, fs});
  fractions_cleanup (cs, fs);
  
  //Obain potential
  printf ("# Start solving for the Potential\n");
  NITERMIN = 10;
  mgstats mg = poisson (pot, laplace, alpha = fs,
			embed_flux = embed_flux, tolerance = HUGE);
  printf ("# Solved with a max. residual %g\n", mg.resa);

  //Get potential flow field
  foreach() {
    foreach_dimension()
      u.x[] = center_gradient(pot) * cs[];
  }
  boundary ((scalar*){u});

  //Place tracer particles
  P = init_tp_square (n = 30, ym = H + Hg, l = 2*R1, zp = -L0/4);
  DT = 0.1;

  //free stuff!
  free (trees);
  free (nodes);
  free (gradp);
}

event set_dtmax (i++) {
  dt = dtnext(DT);
}

double th = 0;
event movie (t += 1.) {
  view (bg = {0.3, 0.3, 0.9}, theta = th, phi = 0.2);
  view (ty = -.3);
  draw_vof ("cs", "fs", fc = {0.73,0.4,0.2});
  translate (y = 0.01)
    draw_vof ("cp", "fp", fc = {0.2,0.9,0.2});
  scatter (P);
  box();
  save ("tree.mp4");
  th += 0.006;
}

event stop (t = 500);