function dom = plate_domain(plate, nx, ny)
% PLATE_DOMAIN  The footprint of the plate on a uniform nx x ny grid.
%
%   dom = plate_domain(plate, nx, ny) lays the grid over the bounding box
%   [0, Lx] x [0, Ly] of plate and marks which nodes belong to the plate.
%   For plate.shape = 'rectangle' every node does.  For 'ellipse' the plate is
%   the ellipse inscribed in the box, (x - Lx/2)^2 / (Lx/2)^2 +
%   (y - Ly/2)^2 / (Ly/2)^2 <= 1, a disk when Lx = Ly, and its nodes are those
%   within half a grid spacing of the inside, so that the staircase boundary
%   they form straddles the true edge rather than lying to one side of it.
%   Nodes that do not belong to a grid cell with all four corners inside are
%   dropped, so every node carries some area.  Both finite-difference methods
%   assemble their operators on this footprint (plate_matrices); the modal
%   method has closed-form modes only for the rectangle.
%
%   Fields.
%     shape, nx, ny, dx, dy, x, y, X, Y   the grid; w(i, j) lives at (x(j), y(i))
%     inside      ny x nx logical, the plate's nodes
%     cell        (ny-1) x (nx-1) logical, cells with all four corners inside
%     m           ny x nx trapezoidal node weights: adjacent plate cells / 4
%                 (1 inside, 1/2 on a straight edge, 1/4 at a convex corner)
%     boundary    inside nodes with a 4-neighbour outside: the ring where
%                 w = 0 is imposed for supported and clamped edges
%     interior    inside and not boundary
%     full3       inside nodes whose eight neighbours are inside, where the
%                 curvatures can be taken by central differences
%     near        ny x nx index of the nearest full3 node, for every inside
%                 node (0 outside); plate_curvature_domain copies from it
%     nunknowns_dirichlet, nunknowns_free   sizes of the two unknown sets
%     kappa, normal_x, normal_y   ellipse only: curvature and outward unit
%                 normal of the true edge nearest each boundary node (zero
%                 elsewhere); the simply supported ellipse needs them, see
%                 plate_matrices
%
%   Accuracy.  A staircase approximates a curved edge to first order in the
%   grid spacing: the geometry, the edge conditions and the stress near the
%   edge all carry O(dx) errors there, against O(dx^2) inside.  Refining the
%   grid is the check.

if ~isfield(plate, 'shape')
    shape = 'rectangle';
else
    shape = plate.shape;
end
if ~any(strcmp(shape, {'rectangle', 'ellipse'}))
    error('plate_domain: unknown plate shape "%s" (expected rectangle or ellipse)', shape);
end
Lx = plate.Lx; Ly = plate.Ly;
dx = Lx / (nx - 1);
dy = Ly / (ny - 1);
x = linspace(0, Lx, nx);
y = linspace(0, Ly, ny);
[X, Y] = meshgrid(x, y);
N = nx * ny;

kappa = zeros(ny, nx);
normal_x = zeros(ny, nx);
normal_y = zeros(ny, nx);
if strcmp(shape, 'rectangle')
    in0 = true(ny, nx);
else
    xc = Lx / 2; yc = Ly / 2; a = Lx / 2; b = Ly / 2;
    phi = ((X - xc) / a).^2 + ((Y - yc) / b).^2;          % 1 on the edge
    gx = 2 * (X - xc) / a^2; gy = 2 * (Y - yc) / b^2;      % gradient of phi
    g = sqrt(gx.^2 + gy.^2);
    g(g == 0) = eps;                                        % the centre
    d = (phi - 1) ./ g;                                     % signed distance to the edge, first order
    hn = 1 ./ sqrt((gx ./ g / dx).^2 + (gy ./ g / dy).^2); % grid spacing along the normal
    in0 = d <= 0.5 * hn;                                    % round the edge to the nearest node
end
cell = in0(1:end - 1, 1:end - 1) & in0(2:end, 1:end - 1) & in0(1:end - 1, 2:end) & in0(2:end, 2:end);
cnt = zeros(ny, nx);
cnt(1:end - 1, 1:end - 1) = cnt(1:end - 1, 1:end - 1) + cell;
cnt(2:end, 1:end - 1) = cnt(2:end, 1:end - 1) + cell;
cnt(1:end - 1, 2:end) = cnt(1:end - 1, 2:end) + cell;
cnt(2:end, 2:end) = cnt(2:end, 2:end) + cell;
inside = cnt > 0;
m = cnt / 4;

% Neighbourhoods, with everything off the grid counting as outside
pad = false(ny + 2, nx + 2);
pad(2:end - 1, 2:end - 1) = inside;
nb4 = pad(1:end - 2, 2:end - 1) & pad(3:end, 2:end - 1) & pad(2:end - 1, 1:end - 2) & pad(2:end - 1, 3:end);
nb8 = nb4 & pad(1:end - 2, 1:end - 2) & pad(1:end - 2, 3:end) & pad(3:end, 1:end - 2) & pad(3:end, 3:end);
interior = inside & nb4;
boundary = inside & ~nb4;
full3 = inside & nb8;
if ~any(full3(:))
    error('plate_domain: the %d x %d grid is too coarse for this %s; no node has a full stencil inside the plate', nx, ny, shape);
end

% Nearest node with a full stencil, by flooding outwards from those nodes
near = zeros(ny, nx);
kf = find(full3);
near(kf) = kf;
todo = inside & ~full3;
for pass = 1:max(nx, ny)
    if ~any(todo(:))
        break
    end
    np = zeros(ny + 2, nx + 2);
    np(2:end - 1, 2:end - 1) = near;
    shifts = [-1 0; 1 0; 0 -1; 0 1];
    for k = 1:4
        sh = np(2 + shifts(k, 1):end - 1 + shifts(k, 1), 2 + shifts(k, 2):end - 1 + shifts(k, 2));
        take = todo & sh > 0;
        near(take) = sh(take);
        todo(take) = false;
    end
end

% Ellipse: curvature and normal of the true edge at the point nearest each
% boundary node (the radial projection of the node onto the ellipse)
if strcmp(shape, 'ellipse')
    kb = find(boundary);
    t = 1 ./ sqrt(max(phi(kb), eps));
    xe = xc + (X(kb) - xc) .* t;
    ye = yc + (Y(kb) - yc) .* t;
    px = 2 * (xe - xc) / a^2; py = 2 * (ye - yc) / b^2;   % gradient of phi on the edge
    pg = sqrt(px.^2 + py.^2);
    normal_x(kb) = px ./ pg;
    normal_y(kb) = py ./ pg;
    kappa(kb) = ((2 / a^2) * py.^2 + (2 / b^2) * px.^2) ./ pg.^3;   % curvature of the level set
end

dom.shape = shape;
dom.nx = nx; dom.ny = ny; dom.dx = dx; dom.dy = dy;
dom.x = x; dom.y = y; dom.X = X; dom.Y = Y;
dom.inside = inside; dom.cell = cell; dom.m = m;
dom.boundary = boundary; dom.interior = interior; dom.full3 = full3; dom.near = near;
dom.nunknowns_dirichlet = nnz(interior);
dom.nunknowns_free = nnz(inside);
dom.kappa = kappa; dom.normal_x = normal_x; dom.normal_y = normal_y;
dom.N = N;
end