#include "grid/octree.h"
#include "utils.h"
#include "distance.h"
#include "fractions.h"

#include "embed.h"
#include "navier-stokes/centered.h"

#include "view.h"

double U_IN = 0.01;

face vector muv[];

int main()
{
    init_grid(16);
    size(10 * 0.15);
    origin(-0.3, -L0 / 2, -L0 / 2);
    mu = muv;
    run();
}

u.n[left] = dirichlet(U_IN);
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.);
u.r[embed] = dirichlet(0.);

event properties(i++)
{
    foreach_face()
        muv.x[] = 0.001;
}

event init(t = 0)
{
    coord *point = input_stl(fopen("../six_spheres.stl", "r"));
    coord min, max;
    bounding_box(point, &min, &max);
    double maxl = -HUGE;
    foreach_dimension()
    {
        if (max.x - min.x > maxl)
            maxl = max.x - min.x;
    }

    scalar dist[];
    distance(dist, point);
    while (adapt_wavelet({dist}, {0.001 * maxl}, 8, 4).nf)
        ;

    solid(cs, fs, (dist[] + dist[-1] + dist[0, -1] + dist[-1, -1] + dist[0, 0, -1] + dist[-1, 0, -1] + dist[0, -1, -1] + dist[-1, -1, -1]) / 8.);

    // init vel
    foreach ()
        u.x[] = cs[] ? U_IN : 0.;

    view(camera = "iso", bg = {1, 1, 1},
         width = 1000, height = 1000);
    cells(n = {0, 0, 1});
    cells(n = {0, 1, 0});
    box();
    save("mesh_init.png");
}

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

event adapt(i++)
{
    adapt_wavelet({cs, u},
                  {1e-2, 0.02, 0.02, 0.02},
                  8, 4);
}

event movies(t = 0; t += 0.01; t <= 0.1)
{
    scalar u_mag[];
    foreach ()
        u_mag[] = norm(u);

    view(camera = "iso", bg = {1, 1, 1},
         width = 1000, height = 1000);
    squares("u_mag", n = {0, 0, 1}, map = cool_warm);
    squares("u_mag", n = {0, 1, 0}, map = cool_warm);
    box();
    save("movie.mp4");
}