% Driver for the L-shaped domain.
%
% The editable script supplies
%
%     names = method_names()
%     v = solve_l(mid, n, g, px, py)
%
% and must return the computed solution at the points (px, py), whatever grid
% or basis it used internally. Comparing two methods at a common set of points,
% rather than each on its own grid, is what makes a finite difference scheme
% and a basis expansion comparable at all.
%
% g(x, y) is the Dirichlet data, which here is the exact solution restricted to
% the boundary. n is the method's own resolution parameter: intervals per unit
% length for the difference scheme, number of basis terms for the expansion.
p = lp_params();
kind = p.kind;
names = method_names();
nm = numel(names);
g = @(xx, yy) lp_lexact(kind, xx, yy);

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

if strcmp(p.mode, 'field')
    % a picture on a fixed grid of spacing 1/32, masked outside the domain
    m = 32;
    xs = (-m:m)' / m;
    [X, Y] = meshgrid(xs, xs);
    inside = ~(X > 1e-12 & Y < -1e-12);
    px = X(inside);
    py = Y(inside);
    t0 = tic;
    v = solve_l(p.mid, p.n, g, px, py);
    out.ms = toc(t0) * 1000;
    ve = lp_lexact(kind, px, py);
    out.maxerr = max(abs(v - ve));

    U = NaN(2*m+1, 2*m+1);
    E = NaN(2*m+1, 2*m+1);
    U(inside) = v;
    E(inside) = abs(v - ve);
    out.x = xs';
    out.y = xs';
    out.U = lp_rows(U);
    out.Err = lp_rows(E);
else
    % errors at the nodes of a fixed grid of spacing 1/8, well away from the
    % plotting grid and coarse enough that every difference grid contains it
    m = 8;
    xs = (-m:m)' / m;
    [X, Y] = meshgrid(xs, xs);
    inside = ~(X > 1e-12 & Y < -1e-12) & abs(X) < 1 - 1e-12 & abs(Y) < 1 - 1e-12;
    px = X(inside);
    py = Y(inside);
    ve = lp_lexact(kind, px, py);

    % each method gets its own list of resolutions: intervals per unit length
    % for the difference scheme, basis terms for the expansion
    % numbl's jsondecode turns a JSON array of arrays into a cell of column
    % vectors, where MATLAB would give a matrix; unpack it by hand
    NSc = p.nsm;
    nn = numel(NSc{1});
    NS = zeros(nm, nn);
    for k = 1:nm
        row = NSc{k};
        NS(k, :) = row(:)';
    end
    E = zeros(nm, nn);
    T = zeros(nm, nn);
    Q = zeros(nm, nn);
    ETall = drv_methodsl(1, NaN(2*nm, nn), NS, g, px, py, ve, nm);
    for k = 1:nm
        E(k, :) = ETall(2*k-1, :);
        T(k, :) = ETall(2*k, :);
    end
    for k = 1:nm
        Q(k, :) = lp_orders(E(k, :), NS(k, :));
    end
    out.ns = NS(1, :);
    out.NS = lp_rows(NS);
    out.E = lp_rows(E);
    out.T = lp_rows(T);
    out.orders = lp_rows(Q);
end

lp_write(out);

function ET = drv_methodsl(k, ET, NS, g, px, py, ve, nm)
% Recursion rather than a loop: see the note in the one-dimensional driver.
if k > nm
    return;
end
vw = solve_l(k, NS(k, 1), g, px(1), py(1));     % warm-up
R = drv_resl(1, NaN(2, size(NS, 2)), k, NS(k, :), g, px, py, ve);
ET(2*k-1, :) = R(1, :);
ET(2*k, :) = R(2, :);
ET = drv_methodsl(k+1, ET, NS, g, px, py, ve, nm);
end

function R = drv_resl(j, R, mid, ns, g, px, py, ve)
if j > numel(ns)
    return;
end
n = ns(j);
if isfinite(n) && n > 0
    t0 = tic;
    v = solve_l(mid, n, g, px, py);
    R(2, j) = toc(t0) * 1000;
    R(1, j) = max(abs(v - ve));
end
R = drv_resl(j+1, R, mid, ns, g, px, py, ve);
end

function A = drv_store(A, ET, row, k)
for c = 1:size(A, 2)
    A(k, c) = ET(row, c);
end
end