1function [wxx, wyy, wxy] = plate_curvature_free(w, dx, dy, nu)
2% PLATE_CURVATURE_FREE Second derivatives of the displacement of a plate with free edges.
3%
4% [wxx, wyy, wxy] = plate_curvature_free(w, dx, dy, nu) returns the curvatures
5% at every node of the ny x nx field w (rows index y, columns index x).
6% Inside they are the usual central differences. On the edges there is no
7% ghost layer to lean on, so the free-edge conditions themselves supply the
8% missing values:
9%
10% no bending moment on an edge w_nn = -nu w_tt (n normal, t tangent)
11% no twisting moment at a corner w_xy = 0, hence also w_xx = w_yy = 0 there
12%
13% The tangential curvature along an edge is a central difference along the
14% edge, and the twist on an edge is one-sided in the normal direction. Both
15% direct stresses therefore vanish at a corner, and on an edge the normal
16% stress vanishes while the tangential stress is -E z w_tt.
18[ny, nx] = size(w);
19wxx = zeros(ny, nx);
20wyy = zeros(ny, nx);
21wxy = zeros(ny, nx);
22wxx(:, 2:end - 1) = (w(:, 3:end) - 2 * w(:, 2:end - 1) + w(:, 1:end - 2)) / dx^2;
23wyy(2:end - 1, :) = (w(3:end, :) - 2 * w(2:end - 1, :) + w(1:end - 2, :)) / dy^2;
24% Edge columns and rows: the moment condition, using the tangential curvature just computed
25wxx(2:end - 1, [1 end]) = -nu * wyy(2:end - 1, [1 end]);
26wyy([1 end], 2:end - 1) = -nu * wxx([1 end], 2:end - 1);
27% Twist: central inside, one-sided in the normal direction on the edges, zero at the corners
28wxy(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);
29wxy(2:end - 1, 1) = (w(3:end, 2) - w(3:end, 1) - w(1:end - 2, 2) + w(1:end - 2, 1)) / (2 * dx * dy);
30wxy(2:end - 1, end) = (w(3:end, end) - w(3:end, end - 1) - w(1:end - 2, end) + w(1:end - 2, end - 1)) / (2 * dx * dy);
31wxy(1, 2:end - 1) = (w(2, 3:end) - w(2, 1:end - 2) - w(1, 3:end) + w(1, 1:end - 2)) / (2 * dx * dy);
32wxy(end, 2:end - 1) = (w(end, 3:end) - w(end, 1:end - 2) - w(end - 1, 3:end) + w(end - 1, 1:end - 2)) / (2 * dx * dy);
33end