% Four methods for the Dirichlet problem on the Pac-Man domain: the unit disk
% with the quadrant x > 0, y < 0 removed. The boundary is a curved arc and two
% straight edges; the edges meet the arc in two right-angle corners, and each
% other in a re-entrant corner of interior angle 3*pi/2 at the origin. Every
% difficulty this tutorial has met is present at once: a boundary that cuts
% through grid cells (Chapter 10's curve), and a corner singularity r^(2/3)
% (Chapter 9's corner).
%
%   1. SHORTLEY-WELLER + SOR. The five-point scheme, with the legs of the
%      stencil shortened wherever a neighbour falls outside the arc, so the
%      boundary condition is imposed on the arc itself and not on a staircase
%      approximation of it. The shortened stencil is no longer symmetric, so
%      conjugate gradients is not available; the solver is SOR with the
%      overrelaxation parameter near its optimum.
%
%   2. THE HARMONIC BASIS. The method of particular solutions from Chapter 9,
%      unchanged in spirit: harmonic polynomials Re z^m, Im z^m for the smooth
%      part, plus the corner's own functions r^(2k/3) sin(2k theta/3), fitted
%      to the boundary data by least squares.
%
%   3. DOUBLE LAYER, GRADED PANELS. The Nystrom method of Chapter 10 cannot
%      use the trapezoid rule here, because the boundary is not smooth: the
%      density has its own singularity at each corner. The quadrature is
%      composite Gauss-Legendre on panels refined dyadically towards the three
%      corners, which resolves the density algebraically panel by panel and
%      geometrically in the refinement depth.
%
%   4. CHEBYSHEV IN POLAR COORDINATES. In polar coordinates the Pac-Man is a
%      rectangle, (r, theta) in (0,1) x (0, 3*pi/2), and Laplace's equation
%      becomes r^2 u_rr + r u_r + u_tt = 0, which is separable: the matrix
%      decomposition of Chapter 7 applies verbatim. The whole side r = 0 of
%      the rectangle maps to the corner, where u is known, so the coordinate
%      degeneracy costs nothing to impose; what it costs is accuracy, because
%      the r^(2/3) singularity now sits on the edge of the expansion domain.

function names = method_names()
names = {'Shortley-Weller + SOR', 'harmonic basis, least squares', ...
    'double layer, graded panels', 'Chebyshev, polar rectangle'};
end

function v = solve_pac(mid, n, g, px, py)
if mid == 1
    v = fd_pac(n, g, px, py);
elseif mid == 2
    v = mps_pac(n, g, px, py);
elseif mid == 3
    v = bie_pac(n, g, px, py);
else
    v = spec_pac(n, g, px, py);
end
end

function d = dof_pac(mid, n)
if mid == 1
    d = fd_count(n);
elseif mid == 2
    d = 3*n + 1;
elseif mid == 3
    [qx, qy, qnx, qny, qw, qkap] = pac_panels(n);
    d = numel(qx);
else
    d = (n - 1)^2;
end
end

% --------------------------------------------------- finite differences -----

function v = fd_pac(N, g, px, py)
% N intervals per unit length on (-1:1/N:1)^2. Nodes on the two straight
% edges are grid points and carry the data directly; nodes adjacent to the
% arc get Shortley-Weller legs, shortened to where the arc actually is.
h = 1 / N;
m = 2*N + 1;
xs = (-N:N)' * h;
[X, Y] = meshgrid(xs, xs);
R2 = X.^2 + Y.^2;

% start from the data everywhere: boundary nodes keep it, and for the
% unknowns it is a continuous extension, which is a decent initial guess
U = g(X, Y);
inquad = X > 1e-12 & Y < -1e-12;
U(inquad) = 0;
onedge = (abs(Y) < 1e-12 & X > -1e-12) | (abs(X) < 1e-12 & Y < 1e-12);
unknown = R2 < 1 - 1e-12 & ~inquad & ~onedge;
mask = double(unknown);

% Shortley-Weller coefficients. Every leg is h except where the neighbour
% lies beyond the arc; there the leg is shortened to the crossing and the
% boundary value at the crossing moves to the right-hand side.
CE = zeros(m, m);
CW = zeros(m, m);
CN = zeros(m, m);
CS = zeros(m, m);
CP = ones(m, m);
B = zeros(m, m);
for j = 2:m-1
    for i = 2:m-1
        if mask(i, j) == 0
            continue;
        end
        x = xs(j);
        y = xs(i);
        hE = h; hW = h; hN = h; hS = h;
        gE = 0; gW = 0; gN = 0; gS = 0;
        if (x + h)^2 + y^2 >= 1 - 1e-12
            hE = max(sqrt(1 - y^2) - x, 1e-8 * h);
            gE = 1;
        end
        if (x - h)^2 + y^2 >= 1 - 1e-12
            hW = max(x + sqrt(1 - y^2), 1e-8 * h);
            gW = 1;
        end
        if x^2 + (y + h)^2 >= 1 - 1e-12
            hN = max(sqrt(1 - x^2) - y, 1e-8 * h);
            gN = 1;
        end
        if x^2 + (y - h)^2 >= 1 - 1e-12
            hS = max(y + sqrt(1 - x^2), 1e-8 * h);
            gS = 1;
        end
        cE = 2 / (hE * (hE + hW));
        cW = 2 / (hW * (hE + hW));
        cN = 2 / (hN * (hN + hS));
        cS = 2 / (hS * (hN + hS));
        CP(i, j) = cE + cW + cN + cS;
        % a shortened leg reaches a point of the arc, not a grid node: its
        % (known) value multiplies the coefficient and moves to the rhs
        if gE == 1
            B(i, j) = B(i, j) + cE * g(x + hE, y);
        else
            CE(i, j) = cE;
        end
        if gW == 1
            B(i, j) = B(i, j) + cW * g(x - hW, y);
        else
            CW(i, j) = cW;
        end
        if gN == 1
            B(i, j) = B(i, j) + cN * g(x, y + hN);
        else
            CN(i, j) = cN;
        end
        if gS == 1
            B(i, j) = B(i, j) + cS * g(x, y - hS);
        else
            CS(i, j) = cS;
        end
    end
end

% SOR, swept until the largest Gauss-Seidel correction is far below the
% discretisation error. The shortened stencils make the matrix
% nonsymmetric, which is what rules out the CG solver of Chapter 9.
om = 2 / (1 + sin(pi * h / 2));
maxit = 60 * N;
for it = 1:maxit
    du = 0;
    for j = 2:m-1
        for i = 2:m-1
            if mask(i, j) > 0
                unew = (CE(i,j) * U(i,j+1) + CW(i,j) * U(i,j-1) ...
                    + CN(i,j) * U(i+1,j) + CS(i,j) * U(i-1,j) + B(i,j)) / CP(i,j);
                d = unew - U(i, j);
                U(i, j) = U(i, j) + om * d;
                if abs(d) > du
                    du = abs(d);
                end
            end
        end
    end
    if du < 1e-9
        break;
    end
end

% the requested points are grid points, so sampling is a lookup
v = zeros(numel(px), 1);
for k = 1:numel(px)
    i = round(py(k) / h) + N + 1;
    j = round(px(k) / h) + N + 1;
    v(k) = U(i, j);
end
end

function d = fd_count(N)
h = 1 / N;
xs = (-N:N)' * h;
[X, Y] = meshgrid(xs, xs);
R2 = X.^2 + Y.^2;
inquad = X > 1e-12 & Y < -1e-12;
onedge = (abs(Y) < 1e-12 & X > -1e-12) | (abs(X) < 1e-12 & Y < 1e-12);
d = sum(sum(double(R2 < 1 - 1e-12 & ~inquad & ~onedge)));
end

% ------------------------------------------- the method of particular sols --

function v = mps_pac(n, g, px, py)
% Harmonic polynomials up to degree n plus n corner functions, fitted by
% least squares at about three boundary points per unknown.
nb = 3*n + 1;
[bx, by] = pac_bpoints(3 * nb);
A = pac_basis(bx, by, n, n);
c = A \ g(bx, by);
v = pac_basis(px, py, n, n) * c;
end

function A = pac_basis(x, y, M, K)
% Re z^m, Im z^m for m = 0..M, then r^(2k/3) sin(2k theta/3) for K values of
% k that are not multiples of 3 (those already are harmonic polynomials, and
% listing them twice would make the matrix exactly rank deficient).
np = numel(x);
A = zeros(np, 2*M + 1 + K);
z = x(:) + 1i * y(:);
col = 1;
p = ones(np, 1);
for mdeg = 0:M
    A(:, col) = real(p);
    col = col + 1;
    if mdeg > 0
        A(:, col) = imag(p);
        col = col + 1;
    end
    p = p .* z;
end
th = atan2(y(:), x(:));
th(th < -1e-13) = th(th < -1e-13) + 2*pi;
r = abs(z);
k = 0;
q = 0;
while q < K
    k = k + 1;
    if mod(k, 3) ~= 0
        q = q + 1;
        A(:, col) = r.^(2*k/3) .* sin(2*k*th/3);
        col = col + 1;
    end
end
end

function [bx, by] = pac_bpoints(ntot)
% Boundary points: on each straight edge, clustered quadratically towards
% the re-entrant corner, where the data varies fastest; on the arc, spaced
% uniformly in angle. Points per piece in proportion to its length.
ne = max(6, round(ntot / 6.7));
na = max(12, ntot - 2*ne);
t = ((1:ne)' - 0.5) / ne;
c = t.^2;
s = 3*pi/2 * ((1:na)' - 0.5) / na;
bx = [c; cos(s); zeros(ne, 1)];
by = [zeros(ne, 1); sin(s); -c];
end

% --------------------------------------------- boundary integral, panels ----

function v = bie_pac(d, g, px, py)
% The double-layer representation and jump relation of Chapter 10, with the
% trapezoid rule replaced by 12-point Gauss-Legendre panels, refined
% dyadically to depth d towards each of the three corners. On the straight
% edges the kernel between two points of the same edge is identically zero,
% and on the arc it is the constant -1/(4*pi); everything nontrivial in this
% matrix happens where two pieces meet at a corner.
[qx, qy, qnx, qny, qw, qkap] = pac_panels(d);
n = numel(qx);

DX = qx * ones(1, n) - ones(n, 1) * qx';
DY = qy * ones(1, n) - ones(n, 1) * qy';
R2 = DX.^2 + DY.^2;
R2(1:n+1:end) = 1;                          % keep 0/0 off the diagonal
K = ((ones(n, 1) * qnx') .* DX + (ones(n, 1) * qny') .* DY) ./ (2*pi*R2);
K(1:n+1:end) = -qkap / (4*pi);              % the smooth-limit diagonal
A = -0.5 * eye(n) + K .* (ones(n, 1) * qw');

sigma = A \ g(qx, qy);

np = numel(px);
EX = px(:) * ones(1, n) - ones(np, 1) * qx';
EY = py(:) * ones(1, n) - ones(np, 1) * qy';
KE = ((ones(np, 1) * qnx') .* EX + (ones(np, 1) * qny') .* EY) ...
    ./ (2*pi*(EX.^2 + EY.^2));
v = (KE .* (ones(np, 1) * qw')) * sigma;
end

function [qx, qy, qnx, qny, qw, qkap] = pac_panels(d)
% Panel endpoints on each piece: dyadic towards both ends, with no panel
% longer than 0.35, then 12-point Gauss-Legendre on every panel. All three
% pieces are parametrised by arclength, so the speed is 1 throughout.
[gx, gw] = lp_gauss(12);
be = graded_breaks(1, d, 0.35);
ba = graded_breaks(3*pi/2, d, 0.35);

[t1, w1] = panel_nodes(be, gx, gw);
[t2, w2] = panel_nodes(ba, gx, gw);
[t3, w3] = panel_nodes(be, gx, gw);

% edge theta = 0 from the corner to (1,0); the arc; the edge theta = 3*pi/2
% from (0,-1) back towards the corner. Outward normals for a domain kept on
% the left of this traversal.
qx = [t1; cos(t2); zeros(numel(t3), 1)];
qy = [zeros(numel(t1), 1); sin(t2); t3 - 1];
qnx = [zeros(numel(t1), 1); cos(t2); ones(numel(t3), 1)];
qny = [-ones(numel(t1), 1); sin(t2); zeros(numel(t3), 1)];
qw = [w1; w2; w3];
qkap = [zeros(numel(t1), 1); ones(numel(t2), 1); zeros(numel(t3), 1)];
end

function bp = graded_breaks(T, d, base)
% Breakpoints on [0, T]: halve towards each end d times, then split anything
% still longer than base into equal pieces.
half = T / 2;
left = [0, half * 2.^(-d:0)];
bp0 = [left, T - left(end-1:-1:1)];
bp = bp0(1);
for k = 2:numel(bp0)
    len = bp0(k) - bp0(k-1);
    ns = max(1, ceil(len / base));
    for q = 1:ns
        bp = [bp, bp0(k-1) + len * q / ns];
    end
end
end

function [t, w] = panel_nodes(bp, gx, gw)
np = numel(bp) - 1;
p = numel(gx);
t = zeros(np * p, 1);
w = zeros(np * p, 1);
for k = 1:np
    a = bp(k);
    b = bp(k+1);
    t((k-1)*p+1 : k*p) = (a + b)/2 + (b - a)/2 * gx;
    w((k-1)*p+1 : k*p) = (b - a)/2 * gw;
end
end

% --------------------------------------- Chebyshev on the polar rectangle ---

function v = spec_pac(N, g, px, py)
% Collocation on the tensor Chebyshev grid of the rectangle
% (r, theta) in [0,1] x [0, 3*pi/2], where Laplace's equation reads
% r^2 u_rr + r u_r + u_tt = 0. The side r = 1 carries the arc data, the
% sides theta = 0 and theta = 3*pi/2 the edge data, and the side r = 0 is
% the corner itself, where u equals the corner value of the data.
[D, s] = lp_cheb(N);
r = (s + 1) / 2;                 % from 1 down to 0
th = (s + 1) * (3*pi/4);         % from 3*pi/2 down to 0
Dr = 2 * D;
Dt = (4 / (3*pi)) * D;
Ar = diag(r.^2) * (Dr * Dr) + diag(r) * Dr;
Bt = Dt * Dt;

G = zeros(N+1, N+1);
G(1, :) = g(cos(th'), sin(th'));            % r = 1: the arc
G(N+1, :) = g(0, 0) * ones(1, N+1);         % r = 0: the corner
G(:, 1) = g(zeros(N+1, 1), -r);             % theta = 3*pi/2
G(:, N+1) = g(r, zeros(N+1, 1));            % theta = 0

ii = 2:N;
nn = N - 1;
A = Ar(ii, ii);
Bi = Bt(ii, ii);
F = -(Ar(ii, [1, N+1]) * G([1, N+1], ii)) - (G(ii, [1, N+1]) * Bt(ii, [1, N+1])');

% the matrix decomposition of Chapter 7, with two different one-dimensional
% operators instead of one
[V, LA] = eig(A);
[W, LB] = eig(Bi);
la = diag(LA);
lb = diag(LB);
Ft = V \ (F / W.');
Ut = Ft ./ (la * ones(1, nn) + ones(nn, 1) * lb.');
U = G;
U(ii, ii) = real(V * Ut * W.');

% evaluate off the grid by barycentric interpolation in each direction
rr = sqrt(px.^2 + py.^2);
tt = atan2(py, px);
tt(tt < -1e-13) = tt(tt < -1e-13) + 2*pi;
v = bary2(U, s, 2*rr - 1, tt / (3*pi/4) - 1);
end

function v = bary2(U, s, sr, st)
% Tensor barycentric interpolation on the Chebyshev grid s (both
% directions), rows of U indexed by the first coordinate.
n1 = numel(s);
wb = [0.5; ones(n1-2, 1); 0.5] .* (-1).^((0:n1-1)');
np = numel(sr);
v = zeros(np, 1);
for k = 1:np
    dr = sr(k) - s;
    dt = st(k) - s;
    hr = abs(dr) < 1e-14;
    ht = abs(dt) < 1e-14;
    if any(hr)
        cr = double(hr);
    else
        cr = wb ./ dr;
        cr = cr / sum(cr);
    end
    if any(ht)
        ct = double(ht);
    else
        ct = wb ./ dt;
        ct = ct / sum(ct);
    end
    v(k) = cr' * U * ct;
end
end