#include "DLMFD_reverse_Uzawa.h"
#include "navier-stokes/perfs.h"
#include "view.h"

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
double deltau;
scalar un[];

u.n[left] = dirichlet(1);
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);

  TOLERANCE = 1e-4;
  NITERMIN = 2;
  
  run();
}

event init (t = 0)
{

  /* Set dynamic viscosity */
  const face vector muc[] = {1/Re, 1/Re};
  mu = muc;

  /* set density of the flow */ 
  const scalar rhoc[] = 1.;
  rho = rhoc;

  /* initial condition: particles position */
  particle * pp = particles;

  GeomParameter gp = {0};
    
  pp[0].iscube = 0;
  pp[0].iswall = 0;
  /* write the position with respect to 0 */
  gp.center.x = L0/2;
  gp.center.y = L0/2;
  gp.radius   = 0.5;
      
  pp[0].g = gp;
   
  /* particle id */
  pp[0].pnum = 0;

  if (!restore ("restart")) {

    init_file_pointers(pdata, fdata, 0); // 0 - initial simulation / 1 - restarting simulation

    refine (sqrt (sq(x-gp.center.x) + sq(y-gp.center.y)) < 1.2*gp.radius && sqrt (sq(x-gp.center.x) + sq(y-gp.center.y)) > 0.8*gp.radius && level < maxlevel);

  }
  else { // restart
  
    init_file_pointers(pdata, fdata, 1);
 
  }

  foreach()
    un[] = u.x[];
}

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;
}

void theta_omega_dlmfd (const particle * pp, ThetaOmega ** to)
{
  scalar vortz[];
  vorticity (u, vortz);
   
  Array * a = array_new();
  coord center = (pp->g).center;
  for (int ii = 0; ii < pp->s.m; ii++) {

    ThetaOmega t;
    t.omega = interpolate (vortz, pp->s.x[ii], pp->s.y[ii]);
    t.theta = atan2 (pp->s.y[ii] - center.y, pp->s.x[ii] - center.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);
  *to =  p;
}

void omega_zero_dlmfd (FILE * fp, const particle * pp)
{
   /* fixme: this function will not work with MPI because of the
      interpolations which are thread-dependent */
  ThetaOmega * a = NULL;
  theta_omega_dlmfd (pp, &a);

  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++) {
  deltau = change (u.x, un);
  fprintf (stderr, "log output %d %g %d %d %g %g %g %ld\n", i, t, mgp.i, mgu.i, mgp.resa, mgu.resa, deltau, grid->tn);

}

void display_omega (int width, int height)
{
 view (fov = 2.16967, quat = {-0.00418115,0.0542593,-0.000125816,0.998518}, tx = -0.53112, ty = -0.499909, bg = {1,1,1}, width = width, height = height, samples = 1);

  squares ("omega", min = -12, max = 12, linear = 1, map = cool_warm);
}

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

  if (t == 1. || t == 3.) {
    // Fig. 4.
    view (fov = 2.16967, quat = {-0.00418115,0.0542593,-0.000125816,0.998518}, tx = -0.53112, ty = -0.499909, bg = {1,1,1}, width = 966, height = 775, samples = 1);
    
    double max = Re == 9500 ? 40 : 12;
    squares ("omega", min = -max, max = +max, linear = 1, map = cool_warm);
    char name[80];
    sprintf (name, "zoom-%g.png", t);
    save (name);

    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");
}


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

void cpout_dlmfd (FILE * fp, particle * pp)
{
  scalar omega[];
  vorticity(u, omega);
  
  for (int ii = 0; ii < pp->s.m; ii++) {
    GeomParameter gc = pp->g;
    coord c = gc.center;
    double w = atan2 (pp->s.y[ii] - c.y, pp->s.x[ii] - c.x);
    if (w < 0) w += 2*pi; 
    fprintf (fp, "%g %g %g %g \n",
	     pp->s.x[ii], // 1
	     pp->s.y[ii], // 2
	     w, // 3
	     interpolate (omega, pp->s.x[ii], pp->s.y[ii])
	     );
  }
}

event adapt (i++) {

  particle * pp = particles;
  double rhoval = 1.;
  
  DLMFD_subproblem (pp, i, rhoval);

  /* Save forces acting on particles before the adapting the mesh */
  sumLambda (pp, fdata, t, dt, flagfield, DLM_lambda, index_lambda, rhoval);

  static FILE * fp = fopen ("omega", "w");
  omega_zero_dlmfd (fp, pp);

  if (t == 0.5|| t == 0.9|| t ==1.5 || t == 2.5) {
    printf("inside cp block\n");
    char name[80];
    sprintf (name, "cp-%g", t);
    FILE * fp2 = fopen (name, "w");
    cpout_dlmfd (fp2, pp);
    fclose (fp2);
  }
  
  /* Free particle structures (we dont need them anymore) */
  free_particles (pp, NPARTICLES);

  /* Save particles trajectories */
  particle_data(pp, t, i, pdata);

  /* char name[80] = "test-init-starting"; */
  /* scalar * list = {p, index_lambda.x, flagfield}; */
  /* vector * vlist = {u}; */
  /* save_data(list, vlist, i, t, name); */
  /* return 1; */
  
  adapt_wavelet ({flagfield,u}, (double[]){1e-4,cmax,cmax}, maxlevel, 5);
}

event surface_profiles (t = {0.5,0.9,1.5,2.5})
{
  printf("forcing output data at t= {0.5,0.9,1.5,2.5}\n")
}