% Two ways to solve the Dirichlet problem on the L-shaped domain.
%
%   1. the five-point scheme on a uniform grid, solved by conjugate gradients
%      until the residual has fallen by 1e-9, which is far below the
%      discretisation error, so that what is measured is the discretisation
%      and not the solver;
%   2. an expansion in harmonic functions, fitted to the boundary data in the
%      least-squares sense: the method of particular solutions.
%
% The second method's basis is the interesting part. Harmonic polynomials
% Re z^m and Im z^m handle everything that is smooth up to the boundary, and
% they alone would converge slowly, because the solution is not smooth at the
% re-entrant corner. Adding the corner's own singular functions
%
%     r^(2k/3) sin(2k theta / 3),      k = 1, 2, 4, 5, 7, ...
%
% (the multiples of 3 are omitted because they are already harmonic
% polynomials) puts the singularity into the basis instead of trying to
% approximate it. This is the Fox-Henrici-Moler construction, and it is why a
% few dozen basis functions can beat a grid of a hundred thousand points.

function names = method_names()
names = {'5-point + CG', 'harmonic basis, least squares'};
end

function v = solve_l(mid, n, g, px, py)
if mid == 1
    v = fd_l(n, g, px, py);
else
    v = mps_l(n, g, px, py);
end
end

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

function v = fd_l(N, g, px, py)
% N intervals per unit length, so the grid is (-1 : 1/N : 1) in both
% directions and every corner of the domain is a grid point.
h = 1 / N;
m = 2*N + 1;
xs = (-N:N)' * h;
[X, Y] = meshgrid(xs, xs);

inL = ~(X > 1e-12 & Y < -1e-12);                 % the closed domain
% a node is an unknown when it lies in the open domain: inside the square, in
% the L, and not on either edge of the slit
onslit = (abs(Y) < 1e-12 & X > -1e-12) | (abs(X) < 1e-12 & Y < 1e-12);
unknown = inL & abs(X) < 1 - 1e-12 & abs(Y) < 1 - 1e-12 & ~onslit;
mask = double(unknown);

U = zeros(m, m);
U(~unknown) = g(X(~unknown), Y(~unknown));       % boundary values
U(~inL) = 0;                                      % outside: never read

S = cg_solve(U, mask, m, 1e-9, 4000);   % U holds the boundary data, S the solution

% 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) = S(i, j);
end
end

function U = cg_solve(U, mask, m, tol, maxit)
R = neg_residual(U, mask, m);
P = R;
rs = dot2(R, R);
stop = tol * tol * rs;
for it = 1:maxit
    AP = applyA(P, mask, m);
    al = rs / dot2(P, AP);
    U = U + al * P;
    R = R - al * AP;
    rsn = dot2(R, R);
    if rsn < stop
        break;
    end
    P = R + (rsn / rs) * P;
    rs = rsn;
end
end

function Y = applyA(V, mask, m)
Y = zeros(m, m);
for j = 2:m-1
    for i = 2:m-1
        Y(i,j) = mask(i,j) * (4*V(i,j) - V(i-1,j) - V(i+1,j) - V(i,j-1) - V(i,j+1));
    end
end
end

function Y = neg_residual(V, mask, m)
Y = zeros(m, m);
for j = 2:m-1
    for i = 2:m-1
        Y(i,j) = -mask(i,j) * (4*V(i,j) - V(i-1,j) - V(i+1,j) - V(i,j-1) - V(i,j+1));
    end
end
end

function s = dot2(A, B)
s = 0;
n = numel(A);
for k = 1:n
    s = s + A(k) * B(k);
end
end

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

function v = mps_l(n, g, px, py)
M = n;                       % harmonic polynomials up to degree M
K = n;                       % corner functions r^(2k/3) sin(2k theta/3)

nb = (2*M + 1) + K;
[bx, by] = boundary_points(max(8, ceil(nb / 2)));   % about three points per unknown
A = basis_matrix(bx, by, M, K);
c = A \ g(bx, by);           % least squares: many more points than unknowns

v = basis_matrix(px, py, M, K) * c;
end

function A = basis_matrix(x, y, M, K)
np = numel(x);
nb = (2*M + 1) + numel(corner_orders(K));
A = zeros(np, nb);
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);
ks = corner_orders(K);
for q = 1:numel(ks)
    k = ks(q);
    A(:, col) = r.^(2*k/3) .* sin(2*k*th/3);
    col = col + 1;
end
end

function ks = corner_orders(K)
% k = 1, 2, 4, 5, 7, ...: the multiples of 3 are harmonic polynomials already
% in the basis, and including them twice would make the least-squares matrix
% exactly rank deficient.
ks = [];
k = 1;
while numel(ks) < K
    if mod(k, 3) ~= 0
        ks = [ks, k];
    end
    k = k + 1;
end
end

function [bx, by] = boundary_points(nper)
% Points along the six edges of the boundary, walking round from the origin:
%
%   1: y = 0,  x from 0 to 1     the upper edge of the slit
%   2: x = 1,  y from 0 to 1
%   3: y = 1,  x from 1 to -1
%   4: x = -1, y from 1 to -1
%   5: y = -1, x from -1 to 0
%   6: x = 0,  y from -1 to 0    the left edge of the slit
%
% The two edges that meet at the re-entrant corner carry points clustered
% towards it, where the data varies fastest.
t = ((1:nper)' - 0.5) / nper;
c = t .^ 2;
z = zeros(nper, 1);
o = ones(nper, 1);
bx = [c;      o;         1 - 2*t;  -o;       -1 + t;  z];
by = [z;      t;         o;        1 - 2*t;  -o;      -c];
end