load "delaunay";
load "airflow-window";
load "conjugate-gradient";
load "airflow-utils";

% ************************************* %
% THE MAIN PROGRAM                      %
% ************************************* %

function airflow_demo(ignore) =
let
    % Generate the window %
    display = get_environment_variable("DISPLAY");
    win = w_make_window(((50,0),(735,685)),"Fluid Flow",w_cyan,display);

    graph_files =    nesl_path ++ "examples/airfoil";
    (points,edges) = read_graph(graph_files);

    bounding_box = w_bounding_box(points);
    win,box = w_add_box(((110,120),(600,540)),bounding_box,"mesh",w_white,win);
    win = w_add_button(((25,120),(60,30)),"quit",w_red,win);
    win = w_add_button(((25,170),(60,30)),"zoom up",w_red,win);
    win = w_add_button(((25,220),(60,30)),"step",w_red,win);
    win,tbox = w_add_text_box(((25,32),(150,56)),"help",w_cyan,win);
    win,tbox = w_add_text_box(((190,22),(520,76)),"note",w_white,win);

    % INTRO %
    window_state = empty_window_state(win,bounding_box);
    ignore = write_note("
This is an algorithm animation of using the control volume method to simulate
fluid flow over a wing.  It first shows how to set up the linear equations
and then solves them using Conjugate Gradient.@
Any comments are welcome.",
window_state);
    window_state = add_segments(zip(points,points),1,w_dark_gray,window_state);
    fl,window_state = wait_for_step(window_state);
    in if not(fl) then w_kill_window(win) else let

    % DELAUNAY TRIANGULATION %
    window_state = empty_window_state(win,bounding_box);
    ignore = write_note("
First we generate the Delaunay triangulation of the set of points.",
window_state);
    delaunay = delaunay_from_edgelist(points,edges);
    segs = get_from_edges(edges,points);
    window_state = add_segments(segs,1,w_dark_gray,window_state);
    fl,window_state = wait_for_step(window_state);
    in if not(fl) then w_kill_window(win) else let

    % VORONOI DIAGRAM %
    ignore = write_note("Based on the Delaunay triangulation, we generate
the Voronoi diagram of the points.  The Voronoi edges are clipped
so they don't extend outside the boundaries.
This clipping is important for the next step.",window_state);
    voronoi = voronoi_edges_from_delaunay(points,delaunay);
    segs = flatten({{seg: x,seg in cell}: cell in voronoi});
    segs = clean_segments(segs);
    window_state = add_segments(segs,1,w_blue,window_state);
    fl,window_state = wait_for_step(window_state);
    in if not(fl) then w_kill_window(win) else let

    % SET UP EQUATIONS %
    ignore = write_note("We generate a set of linear equations in
terms of the Flow Potential at each node.  The equation for a node
sets the net flow into its voronoi cell to zero.  
This is done by estimating the flow through
each Voronoi edge (green) to be the length of the edge times
the gradient of the flow potential along its Delaunay
edge (red).",window_state);
    a = {cell_equation(i,x0,x1s,cell)
	 : cell in voronoi
	 ; x1s in nested_get(points,delaunay)
	 ; i in [0:#voronoi]
	 ; x0 in points};
    example_point = 1349;
    segs = get_from_edges({example_point,n: 
			   n in delaunay[example_point]},points);
    window_state = add_segments(segs,4,w_red,window_state);
    segs = {seg: x,seg in voronoi[example_point]};
    window_state = add_segments(segs,4,w_green,window_state);
    fl,window_state = wait_for_step(window_state);
    in if not(fl) then w_kill_window(win) else let

    % SET UP BOUNDARY CONDITIONS %
    ignore = write_note("For the boundary cells we set 
the flow into the cell to be proportional to the component of 
its boundary edges perpendicular to the X axis.@
We also set the potential at one boundary node to 0 so that we have
a reference value (the black point on the upper right).",window_state);
    segs = get_from_edges({i,j: i in outside_boundary; 
			   j in rotate(outside_boundary,1)},points);
    b,arrows = generate_bs(outside_boundary,points);
    sink_point = max_val(outside_boundary);
    a = rep(a,[(sink_point,1.0)],sink_point);
    b = rep(b,0.0,sink_point);
    arrows = {center,(.06,0.): center in arrows};
    window_state = add_segments(segs,2,w_green,window_state);
    window_state = add_arrows(arrows,w_red,window_state);
    window_state = add_points([(0.0,points[sink_point])],window_state);
    fl,window_state = wait_for_step(window_state);
    in if not(fl) then w_kill_window(win) else let

    % SOLVE THE EQUATIONS USING CONJUGATE GRADIENT %
    ignore = write_note("We now solve the linear equations using the
Conjugate Gradient iterative method.  The colors represent the Flow
Potential and the flow velocity is the gradient of this potential.
Each pass represents 150 iterations of the Conjugate Gradient method.",
window_state);
    segs = get_from_edges(edges,points);
    window_state = empty_window_state(win,bounding_box);
    window_state = add_segments(segs,1,w_dark_gray,window_state);
    xx = cgsolinit(a,b);
    window_state,x = multisolve(a,b,xx,0,3,150,points,window_state);
    fl,Window_state = wait_for_step(window_state);
    in if not(fl) then w_kill_window(win) else let

    % SHOW FLOWS %
    ignore = write_note("We now generate the flow amplitude and direction
for each Delaunay triangle.  This is done by approximating the gradient
of the flow potential at the center of the triangle by looking 
at the value of the potential at the three corners.  The length of the
arrows represent the amplitude.",window_state);
    segs = get_from_edges(edges,points);
    window_state = empty_window_state(win,bounding_box);
    window_state = add_segments(segs,1,w_dark_gray,window_state);
    triangles = delaunay_triangles_from_edges(delaunay);
    arrows = {gradient_from_triangle(v,p): 
	      v in get_from_triangles(triangles,x);
	      p in get_from_triangles(triangles,points)};
    window_state = add_arrows(arrows,w_purple,window_state);
    fl,Window_state = wait_for_step(window_state);
    in if not(fl) then w_kill_window(win) else let

    % END OF DEMO %
    ignore = write_note("The demo is done: the next step will quit.",
			window_state);
    fl,Window_state = wait_for_step(window_state);

in w_kill_window(win) $

%
set arg_check off;
dump vcode "/afs/cs/project/scandal/fun/src/airflow.vcode" airflow_demo(0);
%