function [wxx, wyy, wxy] = plate_curvature_domain(w, dom)
% PLATE_CURVATURE_DOMAIN  Second derivatives of the displacement on an embedded plate.
%
%   [wxx, wyy, wxy] = plate_curvature_domain(w, dom) returns the curvatures of
%   the ny x nx field w on the footprint dom from plate_domain (rows index y,
%   columns index x).  At every node whose eight neighbours belong to the plate
%   (dom.full3) they are the usual central differences.  The remaining nodes of
%   the plate form the staircase ring along a curved edge, where no single
%   formula is reliable: an edge node of the staircase has neighbours outside
%   the plate in directions that depend on where the true edge happens to cut
%   the grid.  Those nodes take the curvatures of the nearest node with a full
%   stencil (dom.near), a first-order extrapolation that matches the accuracy
%   of the staircase itself.  Outside the plate the curvatures are NaN, so the
%   stress fields drawn from them are blank there.
%
%   For a clamped edge the largest stress sits on the edge, and this puts it
%   about one grid spacing inside; refine the grid to sharpen it.  The
%   rectangle keeps its own edge formulas (plate_curvature, plate_curvature_free).

dx = dom.dx; dy = dom.dy;
[ny, nx] = size(w);
wxx = zeros(ny, nx); wyy = zeros(ny, nx); wxy = zeros(ny, nx);
c = w(2:end - 1, 2:end - 1);
wxx(2:end - 1, 2:end - 1) = (w(2:end - 1, 3:end) - 2 * c + w(2:end - 1, 1:end - 2)) / dx^2;
wyy(2:end - 1, 2:end - 1) = (w(3:end, 2:end - 1) - 2 * c + w(1:end - 2, 2:end - 1)) / dy^2;
wxy(2:end - 1, 2:end - 1) = (w(3:end, 3:end) - w(3:end, 1:end - 2) - w(1:end - 2, 3:end) + w(1:end - 2, 1:end - 2)) / (4 * dx * dy);
ring = find(dom.inside & ~dom.full3);
src = dom.near(ring);
wxx(ring) = wxx(src);
wyy(ring) = wyy(src);
wxy(ring) = wxy(src);
out = ~dom.inside;
wxx(out) = NaN;
wyy(out) = NaN;
wxy(out) = NaN;
end