function [K, L, ops] = plate_bending(dom, held, nu, D)
% PLATE_BENDING  Sparse bending and Laplacian operators on the plate by the energy method.
%
%   [K, L, ops] = plate_bending(dom, held, nu, D) assembles, on every node of
%   the plate (dom.inside from plate_domain), the stiffness matrix K of the
%   bending energy
%
%     U = D/2 * integral( w_xx^2 + w_yy^2 + 2 nu w_xx w_yy + 2 (1-nu) w_xy^2 ) dA
%
%   and the operator L of the Laplacian energy integral |grad w|^2 dA, both as
%   sparse N x N matrices over all nx ny grid nodes (rows and columns of nodes
%   off the plate are empty).  Second differences are taken at the nodes that
%   have both neighbours on the plate, the cross difference at the centre of
%   each cell whose four corners are on the plate, each weighted by its area
%   (the finite-difference energy method; Bushnell, Computers & Structures
%   1973).  Inside the plate K is D times the 13-point biharmonic stencil
%   times the cell area; on the edges the free-edge conditions, no bending
%   moment and no Kirchhoff shear, are the natural conditions of the discrete
%   energy, met to first order in the spacing.  K is symmetric positive
%   semidefinite with the three rigid-body motions 1, x, y as its null space.
%
%   Held nodes (the logical ny x nx array held, the jaws of a vise) are part
%   of the energy but not of the unknowns: their rows and columns are dropped
%   after assembly.  A second difference centred on a held node that reaches a
%   free node then penalizes the slope at the clamp, which is exactly the
%   clamped condition w = w_n = 0; it needs a held node behind every held edge
%   node, which plate_mount guarantees.
%
%   L is the Laplacian energy operator with the natural condition of zero
%   normal slope, each grid edge weighted by its share of cells; the loss
%   term of the plate equation uses it divided by the node areas.
%
%   ops returns the difference operators the prestress term shares
%   (plate_prestress): Bx, By first differences along the x and y edges of the
%   plate (edges x areas), Bxc, Byc the cell-centred gradients, and the edge and
%   cell weights, all on the numbering ops.idx of the plate's nodes.

nx = dom.nx; ny = dom.ny; dx = dom.dx; dy = dom.dy;
N = nx * ny;
inside = dom.inside;
idx = find(inside);
n = numel(idx);
me = (1:n)';
unk = inside & ~held;

% Node number of every grid node (0 off the plate), padded so every node has four neighbours to look up
num_in = zeros(ny, nx);
num_in(idx) = me;
num = zeros(ny + 2, nx + 2);
num(2:end - 1, 2:end - 1) = num_in;
nbL = num(2:end - 1, 1:end - 2); nbR = num(2:end - 1, 3:end);
nbD = num(1:end - 2, 2:end - 1); nbU = num(3:end, 2:end - 1);
uL = nbL(idx); uR = nbR(idx); uD = nbD(idx); uU = nbU(idx);

area = dx * dy;
mu = dom.m(idx) * area;                                   % node areas
hasX = uL > 0 & uR > 0;
hasY = uD > 0 & uU > 0;
kx = find(hasX); nkx = numel(kx);
ky = find(hasY); nky = numel(ky);
Dxx = sparse([(1:nkx)'; (1:nkx)'; (1:nkx)'], [uL(kx); me(kx); uR(kx)], [ones(nkx, 1); -2 * ones(nkx, 1); ones(nkx, 1)] / dx^2, nkx, n);
Dyy = sparse([(1:nky)'; (1:nky)'; (1:nky)'], [uD(ky); me(ky); uU(ky)], [ones(nky, 1); -2 * ones(nky, 1); ones(nky, 1)] / dy^2, nky, n);
Kxx = Dxx' * spdiags(mu(kx), 0, nkx, nkx) * Dxx;
Kyy = Dyy' * spdiags(mu(ky), 0, nky, nky) * Dyy;
kc = find(hasX & hasY); nkc = numel(kc);
Dxc = sparse([(1:nkc)'; (1:nkc)'; (1:nkc)'], [uL(kc); me(kc); uR(kc)], [ones(nkc, 1); -2 * ones(nkc, 1); ones(nkc, 1)] / dx^2, nkc, n);
Dyc = sparse([(1:nkc)'; (1:nkc)'; (1:nkc)'], [uD(kc); me(kc); uU(kc)], [ones(nkc, 1); -2 * ones(nkc, 1); ones(nkc, 1)] / dy^2, nkc, n);
Wc = spdiags(mu(kc), 0, nkc, nkc);
Kc = nu * (Dxc' * Wc * Dyc + Dyc' * Wc * Dxc);

% Cells with four corners on the plate: the twist at the centre
A00 = num_in(1:end - 1, 1:end - 1); A10 = num_in(2:end, 1:end - 1);
A01 = num_in(1:end - 1, 2:end); A11 = num_in(2:end, 2:end);
cu = A00 > 0 & A10 > 0 & A01 > 0 & A11 > 0;
kcell = find(cu); ncell = numel(kcell);
c00 = A00(kcell); c10 = A10(kcell); c01 = A01(kcell); c11 = A11(kcell);
ic = (1:ncell)';
Dxy = sparse([ic; ic; ic; ic], [c00; c10; c01; c11], [ones(ncell, 1); -ones(ncell, 1); -ones(ncell, 1); ones(ncell, 1)] / (dx * dy), ncell, n);
Kxy = 2 * (1 - nu) * area * (Dxy' * Dxy);
Ku = D * (Kxx + Kyy + Kc + Kxy);

% Laplacian energy: first differences along the grid edges, weighted by the
% cell area each edge shares (a whole cell inside, half a cell on the boundary)
cpad = zeros(ny + 1, nx + 1);
cpad(2:end - 1, 2:end - 1) = cu;
WX = (cpad(2:end, 2:end) + cpad(1:end - 1, 2:end)) / 2;   % x-edge from (i, j) to (i, j+1)
WY = (cpad(2:end, 2:end) + cpad(2:end, 1:end - 1)) / 2;   % y-edge from (i, j) to (i+1, j)
ex = find(uR > 0); nex = numel(ex);
ey = find(uU > 0); ney = numel(ey);
Bx = sparse([(1:nex)'; (1:nex)'], [me(ex); uR(ex)], [-ones(nex, 1); ones(nex, 1)] / dx, nex, n);
By = sparse([(1:ney)'; (1:ney)'], [me(ey); uU(ey)], [-ones(ney, 1); ones(ney, 1)] / dy, ney, n);
wx = WX(idx(ex)) * area; wy = WY(idx(ey)) * area;
Lu = -(Bx' * spdiags(wx, 0, nex, nex) * Bx + By' * spdiags(wy, 0, ney, ney) * By);

% Cell-centred gradients: the average of the two x-edges (bottom and top) and the two y-edges (left and right)
xe_of = zeros(n, 1); xe_of(me(ex)) = (1:nex)';
ye_of = zeros(n, 1); ye_of(me(ey)) = (1:ney)';
Cx = sparse([ic; ic], [xe_of(c00); xe_of(c10)], 0.5 * ones(2 * ncell, 1), ncell, nex);
Cy = sparse([ic; ic], [ye_of(c00); ye_of(c01)], 0.5 * ones(2 * ncell, 1), ncell, ney);
Bxc = Cx * Bx;
Byc = Cy * By;

% Drop the held nodes, then embed over all N grid nodes
sel = double(unk(idx));
Pz = spdiags(sel, 0, n, n);
Ku = Pz * Ku * Pz;
Lu = Pz * Lu * Pz;
[ki, kj, kv] = find(Ku);
K = sparse(idx(ki), idx(kj), kv, N, N);
[li, lj, lv] = find(Lu);
L = sparse(idx(li), idx(lj), lv, N, N);

ops.idx = idx; ops.n = n; ops.mu = mu;
ops.Bx = Bx; ops.By = By; ops.Bxc = Bxc; ops.Byc = Byc; ops.Cx = Cx; ops.Cy = Cy;
ops.wx = wx; ops.wy = wy; ops.ncell = ncell; ops.area = area;
ops.unk = sel;
ops.corners = [c00 c10 c01 c11];
% The factors of K: strain operators and the weights between them, so that
% K = D (Dxx' Mx Dxx + Dyy' My Dyy + nu (Dxc' Wc Dyc + Dyc' Wc Dxc) + 2 (1 - nu) A Dxy' Dxy).
% The GPU port applies K in this form (a second difference, weights, a second
% difference transposed) because in single precision the assembled 13-point
% stencil loses the smooth modes to rounding; see plate_fd_setup.
ops.Dxx = Dxx; ops.Dyy = Dyy; ops.Dxc = Dxc; ops.Dyc = Dyc; ops.Dxy = Dxy;
ops.mkx = mu(kx); ops.mky = mu(ky); ops.mkc = mu(kc);
ops.D = D; ops.nu = nu;
end