% Driver for the two-dimensional labs on the unit square.
%
% The editable script supplies
%
%     names = method_names()
%     [x, y, U] = solve2d(mid, N, g)
%
% where g(x, y) is the Dirichlet data. Because every exact solution here is
% harmonic on the whole plane, g happens to be defined everywhere, but a method
% is only entitled to evaluate it on the boundary of the square; anything else
% is looking up the answer.
%
% U is the solution on the tensor grid x by y, with U(i, j) at (x(j), y(i)), and
% the driver measures the error against the exact solution on exactly that grid.
p = lp_params();
kind = p.kind;
par = p.par;
names = method_names();
nm = numel(names);
g = @(xx, yy) lp_exact2d(kind, xx, yy, par);

out = struct();
out.names = names;

if strcmp(p.mode, 'single')
    t0 = tic;
    [x, y, U] = solve2d(p.mid, p.N, g);
    out.ms = toc(t0) * 1000;
    x = x(:);
    y = y(:);
    [X, Y] = meshgrid(x, y);
    UE = lp_exact2d(kind, X, Y, par);
    Err = U - UE;
    out.maxerr = max(max(abs(Err)));
    out.unknowns = (numel(x) - 2) * (numel(y) - 2);
    [Us, xs, ys] = lp_thin(U, x, y, p.maxpix);
    Es = lp_thin(abs(Err), x, y, p.maxpix);
    out.x = xs(:)';
    out.y = ys(:)';
    out.U = lp_rows(Us);
    out.Err = lp_rows(Es);
else
    ns = p.ns(:)';
    nn = numel(ns);
    E = zeros(nm, nn);
    T = zeros(nm, nn);
    Q = zeros(nm, nn);
    nmax = p.nmaxs(:)';
    timed = p.timed;
    ET = drv_methods2(1, NaN(2*nm, nn), ns, g, kind, par, nmax, timed, nm);
    for k = 1:nm
        E(k, :) = ET(2*k-1, :);
        T(k, :) = ET(2*k, :);
        Q(k, :) = lp_orders(E(k, :), ns);
    end
    out.ns = ns;
    out.E = lp_rows(E);
    out.T = lp_rows(T);
    out.orders = lp_rows(Q);
end

lp_write(out);

function ET = drv_methods2(k, ET, ns, g, kind, par, nmax, timed, nm)
% One method per level of the recursion. See the note in the one-dimensional
% driver: a for loop here would put every timed solve on numbl's slow path.
if k > nm
    return;
end
[xw, yw, Uw] = solve2d(k, 8, g);      % warm-up, so the first timing is honest
R = drv_res2(1, NaN(2, numel(ns)), k, ns, g, kind, par, nmax(k), timed);
ET(2*k-1, :) = R(1, :);
ET(2*k, :) = R(2, :);
ET = drv_methods2(k+1, ET, ns, g, kind, par, nmax, timed, nm);
end

function R = drv_res2(j, R, mid, ns, g, kind, par, nmax, timed)
if j > numel(ns)
    return;
end
N = ns(j);
if N <= nmax
    t0 = tic;
    [x, y, U] = solve2d(mid, N, g);
    el = toc(t0);
    if timed > 0 && el < 0.005
        el = drv_best2(4, el, mid, N, g);
    end
    [X, Y] = meshgrid(x(:), y(:));
    R(1, j) = max(max(abs(U - lp_exact2d(kind, X, Y, par))));
    R(2, j) = el * 1000;
end
R = drv_res2(j+1, R, mid, ns, g, kind, par, nmax, timed);
end

function el = drv_best2(k, el, mid, N, g)
if k == 0
    return;
end
t0 = tic;
[x, y, U] = solve2d(mid, N, g);
e2 = toc(t0);
if e2 < el
    el = e2;
end
el = drv_best2(k-1, el, mid, N, g);
end