#include "embed.h"
#include "navier-stokes/centered.h"
#include "navier-stokes/double-projection.h"
#include "navier-stokes/perfs.h"
#include "view.h"

double D = 1.;
int maxlevel = 12;  // 15/16 for Re = 9500, 12 for Re = 1000
double Re = 1000;   // or 9500
double cmax = 3e-3; // 1e-3 for Re = 9500, 3e-3 for Re = 1000

face vector muv[];

double u0 = 1.;

u.n[left] = dirichlet(u0);

u.n[right] = neumann(0);
p[right] = dirichlet(0);
pf[right] = dirichlet(0);

int main (int argc, char * argv[])
{
  if (argc > 1)
    maxlevel = atoi (argv[1]);
  if (argc > 2)
    Re = atof (argv[2]);
  if (argc > 3)
    cmax = atof (argv[3]);

  size (18 [1]);
  origin (- L0/2.);

  mu = muv;
  TOLERANCE = 1e-4 [*];
  NITERMIN = 2;
  
  run();
}

event properties (i++) {
  foreach_face()
    muv.x[] = fm.x[]*u0*D/Re;
}

event init (t = 0)
{

  if (!restore ("restart")) {

    refine (level <= maxlevel*(1. - sqrt(fabs(sq(x) + sq(y) - sq(D/2.)))/2.));

    solid (cs, fs, sq(x) + sq(y) - sq(D/2.));
    
    foreach()
      u.x[] = u0*cs[]; // fixme: with 1 this results in sub-optimal adaptation
  }
  else { // restart

    solid (cs, fs, sq(x) + sq(y) - sq(D/2.));
    
  }

  u.n[embed] = dirichlet(0);
  u.t[embed] = dirichlet(0);  
}

typedef struct {
  double theta, omega;
} ThetaOmega;

int compar_theta (const void * a, const void * b)
{
  const ThetaOmega * p1 = a, * p2 = b;
  return p1->theta > p2->theta ? 1 : -1;
}

ThetaOmega * theta_omega()
{
  Array * a = array_new();
  foreach (serial)
    if (cs[] > 0. && cs[] < 1.) {
      coord n, b;
      embed_geometry (point, &b, &n);
      x += b.x*Delta, y += b.y*Delta;

      ThetaOmega t;
      t.omega = embed_vorticity (point, u, b, n);
      t.theta = atan2(y, x);
      array_append (a, &t, sizeof (ThetaOmega));
    }
  qsort (a->p, a->len/sizeof(ThetaOmega), sizeof(ThetaOmega), compar_theta);
  ThetaOmega t = {nodata, nodata};
  array_append (a, &t, sizeof (ThetaOmega));
  ThetaOmega * p = a->p;
  free (a);
  return p;
}

void omega_zero (FILE * fp)
{
  // fixme: this function will not work with MPI
  ThetaOmega * a = theta_omega();
  for (ThetaOmega * o = a; (o + 1)->theta != nodata; o++) {
    ThetaOmega * o1 = o + 1;
    if (o1->omega*o->omega < 0.)
      fprintf (fp, "%g %g\n", t,
	       o->theta + o->omega*(o1->theta - o->theta)/
	       (o->omega - o1->omega));
  }
  free (a);
  fflush (fp);
}

event logfile (i += 10)
{
  coord Fp, Fmu;
  embed_force (p, u, mu, &Fp, &Fmu);
  fprintf (stderr, "%d %g %g %.6f %.6f %.6f %.6f\n",
	   i, t, dt, Fp.x, Fp.y, Fmu.x, Fmu.y);

  static FILE * fp = fopen ("omega", "w");
  omega_zero (fp);
}

void display_omega (int width, int height)
{
  view (fov = 1.68, tx = -0.0251023,
	ty = 1e-12, // fixme: this is necessary to re-center the view
	width = width, height = height);
  draw_vof ("cs", filled = -1, fc = {1,1,1});
  draw_vof ("cs", lw = 6);
  squares ("omega", min = -12, max = 12, linear = true, map = blue_white_red);
  mirror ({0,1}) {
    draw_vof ("cs", lw = 6);
    squares (color = "level", min = 6, max = maxlevel);
  }
}

event snapshots (t = 1.; t <= 3.; t += 1.)
{
  scalar omega[];
  vorticity (u, omega);

  if (t == 1. || t == 3.) {
    // Fig. 4.
    view (fov = 0.470916, tx = -0.0197934, ty = -0.0258659,
	  width = 524, height = 480);
    draw_vof ("cs", filled = -1, fc = {1,1,1});
    draw_vof ("cs", filled = 0, lw = 1);
    double max = Re == 9500 ? 40 : 12;
    squares ("omega", min = -max, max = +max, linear = 1, map = blue_white_red);
    char name[80];
    sprintf (name, "zoom-%g.png", t);
    save (name);

    draw_vof ("cs", lw = 2);
    squares ("level");
    sprintf (name, "cells-%g.png", t);
    save (name);
  }
  
  if (t == 3.) {
    // Fig. 3.
    display_omega (640, 480);
    save ("omega-3.png");
  }
  
  p.nodump = false;
  char name [80];
  sprintf (name, "dump-%g", t);
  dump (name);
}

event movie (i += 10*(1 << (maxlevel - 12)))
{
  scalar omega[];
  vorticity (u, omega);
  display_omega (1280, 960);
  save ("omega.mp4");
}

void cpout (FILE * fp)
{
  foreach (serial)
    if (cs[] > 0. && cs[] < 1.) {
      coord b, n;
      double area = embed_geometry (point, &b, &n);
      x += b.x*Delta, y += b.y*Delta;

      fprintf (fp, "%g %g %g %g %g %g\n",
	       x, // 1
	       y, // 2
	       atan2(y, x), // 3
	       embed_interpolate (point, p, b), // 4
	       area*Delta, // 5
	       embed_vorticity (point, u, b, n)); // 6
    }
}

event surface_profiles (t = {0.5,0.9,1.5,2.5})
{
  char name[80];
  sprintf (name, "cp-%g", t);
  FILE * fp = fopen (name, "w");
  cpout (fp);
  fclose (fp);
}

#if 0
event snapshot (i += 10) {
  scalar omega[];
  vorticity (u, omega);
  p.nodump = false;
  dump();

  FILE * fp = fopen ("cp", "w");
  cpout (fp);
  fclose (fp);
}
#endif

event adapt (i++) {
  adapt_wavelet ({cs,u}, (double[]){1e-2,cmax,cmax}, maxlevel, 5);
}