#include "grid/octree.h" 		// For 3D
#include "view.h"			// For bview
#include "navier-stokes/centered.h"     // Navier stokes 
// #include "tracer.h"			// Tracers
// #include "diffusion.h"			// Diffusion 

#include "embed.h"
#include "PointTriangle.h"

/** Global variables */
int minlevel, maxlevel;         	// Grid depths
double meps, eps;			// Maximum error and error in u fields

double TEND = 15;
double H = 5, W = 80, D = 2; //Height (y), width (z) and `depth`(x)
double xp = 200, zp = 200;

#define CP 1005.	// C_p for air 
#define gCONST 9.81	// Gravitational constant
#define TREF 273.	// Kelvin
#define INVERSION .3 	// Kelvin per meter
#define karman 0.4      // von Karman constant 

#define roughY0u 0.1    // roughness wind length 
#define roughY0h 0.1     // roughness heat length

#define WIND(s) (max(0.1*log(s/roughY0u),0.))   // log Wind profile 
// #define QFLX 0. 	// 0 (0.001 = 20wm2)
// #define BSURF (1.5*b[] - 0.5*b[0, 1])
// #define GFLX (-Lambda*(BSURF - bd))
// double Lambda = 0.005, bd = 0.;   // Grass coupling
// #define STRAT(s) gCONST/TREF*s*INVERSION

// scalar b[];
// scalar * tracers = {b};

double crho = 1.;	
// face vector av[]; 

// #include "physics.h"			// Physics of the simulation 
// #include "fan.h"			// Include a fan
#include "lambda2.h"

face vector muc[];
double nu;

int main() {	
    nu = 1./300.;
    minlevel = 4;
    maxlevel = 9;

    L0 = 400.;
    X0 = Y0 = Z0 = 0.;
    mu = muc;

    init_grid(1<<8);
    // a = av; 
    // meps = 10.;					// Maximum adaptivity criterion
    // DT = 0.5;					// For poisson solver 
    // TOLERANCE=10E-4;				// For poisson solver 
    // CFL = 0.8;					// CFL condition
    run();						// Start simulation 
}

event init(t=0) {
    eps = .05;
    
	  // u.n[bottom] = dirichlet(0.);
	  // u.t[bottom] = dirichlet(0.);
	  // u.n[top] = dirichlet(0.);
	  // u.t[top] = dirichlet(WIND(y));
    // u.n[embed] = dirichlet(0.);
    // u.t[embed] = dirichlet(0.);
    // periodic (left);
    // #if dimension == 3
    //     u.r[embed] = dirichlet(0.);
    //     periodic(front);
    // #endif 
        
    u.n[left]  = dirichlet(WIND(y));
    p[left]    = neumann(0.);
    pf[left]   = neumann(0.);

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

    u.n[embed] = dirichlet(0.);
    u.t[embed] = dirichlet(0.);
    #if dimension == 3
        u.r[embed] = dirichlet(0.);
        // periodic(back);
    #endif


    vertex scalar phw[];
    D /= 2.;
    H -= D/2;
    W -= W/2;
    foreach_vertex() {
        coord cc = {x, y, z};
        if ((z - zp) > W/2) { // distance to a side edge
        coord p1 = {xp, Y0, zp + W/2.}, p2 = {xp, Y0 + H, zp + W/2.};
        coord tmp1[1]; double tmp2[1];
        phw[] = sqrt (PointSegmentDistance (&cc, &p1, &p2, tmp1, tmp2));
        }
        else if ((z - zp) < -W/2) {
        coord p1 = {xp, Y0, zp - W/2.}, p2 = {xp, Y0 + H, zp - W/2.};
        coord tmp1[1]; double tmp2[1];
        phw[] = sqrt (PointSegmentDistance (&cc, &p1, &p2, tmp1, tmp2));
        }
        else if (y - Y0 > H) {
        coord p1 = {xp, Y0 + H, zp - W/2.}, p2 = {xp, Y0 + H, zp + W/2.};
        coord tmp1[1]; double tmp2[1];
        phw[] = sqrt (PointSegmentDistance (&cc, &p1, &p2, tmp1, tmp2));
        }
        else
        phw[] = fabs(cc.x - xp);
    }
    fractions (phw, cs, fs, D);
    boundary ({cs, fs});
    foreach()
      u.x[] = cs[] ? WIND(y) : 0.;
    // while(adapt_wavelet((scalar *){cs, u},(double []){0.01,eps,eps,eps},maxlevel,minlevel).nf) {
    //   foreach() {
    //       u.x[] = cs[]? WIND(y): 0.;
    //   }
    // }
}

event properties (i++) {
  foreach_face()
    muc.x[] = fm.x[]*nu; 
  boundary ((scalar*){muc});
}

// event inflow(i++){
//     double sides = 50;
//     double relaxtime = dt/50;
//     foreach(){
// 	if((x < sides || x > L0-sides) ||
// 	   (z < sides || z > L0-sides) ||
// 	   (y > L0-2*sides )) {
// 	    double a = (x < sides) ? x : fabs(x-L0);
// 	    a = 1.; 
// 	    u.x[] = u.x[] + a*(WIND(y)-u.x[])*relaxtime;
//  	    // b[] = b[] + a*(STRAT(y) - b[])*relaxtime;
// 	    u.y[] = u.y[] - a*u.y[]*relaxtime;
// 		  u.z[] = u.z[] - a*u.z[]*relaxtime;
// 	}
//     }
// }

// event tracer_diffusion(i++){
//     scalar r[];
//     foreach() {
//         r[] = 0;
//         if (y < Delta)
//             r[] = (QFLX + GFLX)/sq(Delta); // div needed as normalization 
//     }
    
//     double flx = 0, bt = 0;
//     double fctr = CP*TREF/gCONST;
//     foreach_boundary(bottom reduction(+:flx) reduction(+:bt)) {
//         flx = flx + (QFLX + GFLX) * sq(Delta);
//          bt = bt + BSURF * sq(Delta);
//     }
//     bt = bt/sq(L0);
//     flx = flx/sq(L0);
//     fprintf(stderr, "soil=%g %g %g %d\n", t, fctr*flx, fctr*bt/CP, i);  
    
//     mgb = diffusion(b, dt, mu, r = r);
// }

event adapt(i++) {
    adapt_wavelet((scalar *){cs,u},(double []){0.01,eps,eps,eps},maxlevel,minlevel);
}

event progress(i++) {
    fprintf (stderr, "%d %g %d %d\n", i, t, mgp.i, mgu.i);
}

event mov (t += 0.5) {
#if (dimension == 2)
  scalar omega[];
  vorticity (u, omega);
  boundary ({omega});
  // draw_vof ("cs", "fs", filled = -1, fc = {0,0,0});
  // draw_vof ("cs", "fs");
  squares ("omega", linear = true, map = cool_warm);
  mirror ({0,-1})
  cells();
#elif (dimension == 3)
  scalar l2[], vyz[];
  foreach()
    vyz[] = ((u.y[0,0,1] - u.y[0,0,-1]) - (u.z[0,1] - u.z[0,-1]))/(2.*Delta);
  boundary ({vyz});
  lambda2 (u, l2);
  view (fov = 30, theta = 0.5, phi = 0.4, 
	tx = -0.1, ty = 0.1, bg = {65./256,157./256,217./256},
	width = 1080, height = 1080);
  box();
  draw_vof ("cs", "fs", fc = {0.5,0.1,0.2});
  draw_vof ("cs", "fs");
  isosurface ("l2", -0.01);
  translate (y = -5){
    squares ("u.x", n = {0,1,0}, alpha = 5.,
	     min = -1.2, max = 2, map = cool_warm);}
  // translate (z = -L0/2.){
  //   squares ("b", n = {0,0,1}, alpha = L0/2.,
	//      min = -0.1, max = 0.6, map = cool_warm);}
  // translate (x = -L0/2.){
  //   squares ("b", n = {1,0,0}, alpha = L0/2.,
	//      min = -0.1, max = 0.6, map = cool_warm);}
#endif
  save ("lam_3D.mp4");
}

event end(t=TEND) {
}