1function [wxx, wyy, wxy] = plate_curvature_domain(w, dom)
2% PLATE_CURVATURE_DOMAIN Second derivatives of the displacement on an embedded plate.
3%
4% [wxx, wyy, wxy] = plate_curvature_domain(w, dom) returns the curvatures of
5% the ny x nx field w on the footprint dom from plate_domain (rows index y,
6% columns index x). At every node whose eight neighbours belong to the plate
7% (dom.full3) they are the usual central differences. The remaining nodes of
8% the plate form the staircase ring along a curved edge, where no single
9% formula is reliable: an edge node of the staircase has neighbours outside
10% the plate in directions that depend on where the true edge happens to cut
11% the grid. Those nodes take the curvatures of the nearest node with a full
12% stencil (dom.near), a first-order extrapolation that matches the accuracy
13% of the staircase itself. Outside the plate the curvatures are NaN, so the
14% stress fields drawn from them are blank there.
15%
16% For a clamped edge the largest stress sits on the edge, and this puts it
17% about one grid spacing inside; refine the grid to sharpen it. The
18% rectangle keeps its own edge formulas (plate_curvature, plate_curvature_free).
20dx = dom.dx; dy = dom.dy;
21[ny, nx] = size(w);
22wxx = zeros(ny, nx); wyy = zeros(ny, nx); wxy = zeros(ny, nx);
23c = w(2:end - 1, 2:end - 1);
24wxx(2:end - 1, 2:end - 1) = (w(2:end - 1, 3:end) - 2 * c + w(2:end - 1, 1:end - 2)) / dx^2;
25wyy(2:end - 1, 2:end - 1) = (w(3:end, 2:end - 1) - 2 * c + w(1:end - 2, 2:end - 1)) / dy^2;
26wxy(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);
27ring = find(dom.inside & ~dom.full3);
28src = dom.near(ring);
29wxx(ring) = wxx(src);
30wyy(ring) = wyy(src);
31wxy(ring) = wxy(src);
32out = ~dom.inside;
33wxx(out) = NaN;
34wyy(out) = NaN;
35wxy(out) = NaN;
36end