% 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 = NaN(nm, nn);
    T = NaN(nm, nn);
    Q = NaN(nm, nn);
    nmax = p.nmaxs(:)';
    timed = p.timed;
    for k = 1:nm
        [xw, yw, Uw] = solve2d(k, 8, g);            % warm-up
        for j = 1:nn
            N = ns(j);
            if N > nmax(k)
                continue;
            end
            t0 = tic;
            [x, y, U] = solve2d(k, N, g);
            el = toc(t0);
            rep = 0;
            while timed > 0 && el < 0.005 && rep < 4
                rep = rep + 1;
                t0 = tic;
                [x, y, U] = solve2d(k, N, g);
                el = min(el, toc(t0));
            end
            [X, Y] = meshgrid(x(:), y(:));
            E(k, j) = max(max(abs(U - lp_exact2d(kind, X, Y, par))));
            T(k, j) = el * 1000;
        end
        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);