function [Kg, ps] = plate_prestress(P, dom, mnt, ops)
% PLATE_PRESTRESS  In-plane stress of the plate's own weight, and its effect on bending.
%
%   [Kg, ps] = plate_prestress(P, dom, mnt, ops) solves the plane-stress
%   problem of the plate loaded by the in-plane part of its weight,
%   rho h g (gx, gy) per unit area, with the in-plane constraints of the mount
%   (the holes a thread passes through, the jaws of a vise, the pinched patch),
%   and returns the sparse N x N operator Kg of the membrane term
%
%     U_N = 1/2 integral( N_xx w_x^2 + N_yy w_y^2 + 2 N_xy w_x w_y ) dA
%
%   that the in-plane stress resultants N add to the bending energy.  This
%   is how gravity acts on the vibration of a hanging plate: a tilted element
%   under tension is pulled back, and for a rigid rotation of the plate about
%   its hole the term integrates to the compound-pendulum stiffness m g d
%   exactly (d the distance from the hole to the centre of mass), so the
%   swing of a hung plate and the small stiffening of its lowest flexural
%   modes come out of the same operator.  Parts of the plate above the hole
%   are in compression and the term softens them slightly.  Without an
%   in-plane gravity component (a horizontal plate) Kg is empty.
%
%   Discretization.  The in-plane displacement (u, v) lives on the nodes.
%   The direct strains are first differences along the grid edges,
%   e_xx = (u_R - u_L) / dx on the x-edges and e_yy on the y-edges, the shear
%   strain gamma = u_y + v_x and the coupling e_xx e_yy are taken at the cell
%   centres from the averaged corner differences (ops from plate_bending).
%   The strain energy
%
%     U = h/2 sum_edges C11 w_e e_xx^2 + h/2 sum_edges C11 w_e e_yy^2
%       + h sum_cells A_c ( C12 e_xx e_yy + G/2 gamma^2 ),
%     C11 = E / (1 - nu^2), C12 = nu C11, G = E / (2 (1 + nu)),
%
%   with the edge weights w_e their share of cell area, has no spurious zero
%   modes: only the two translations and the rotation are strain-free.  The
%   stress resultants are the derivatives of U with respect to the strains,
%   n_xx = dU/de_xx on each x-edge, n_yy on each y-edge, n_xy = dU/dgamma on
%   each cell (each N times its area), and Kg is assembled from the same
%   difference operators applied to w:
%
%     Kg = Bx' diag(n_xx) Bx + By' diag(n_yy) By + Bxc' diag(n_xy) Byc + Byc' diag(n_xy) Bxc.
%
%   Because both U and U_N use the same strain operators, the identity
%   r' Kg r = m g d for the rigid tilt r holds to rounding, not just to
%   truncation error.
%
%   Constraints (mnt.fix) are imposed by a change of variables: a node with
%   one fixed direction keeps its displacement along the other, a node with
%   two keeps none.  A plate on one thread can still rotate about the hole;
%   the load has no moment about it (the plate hangs with its centre of mass
%   below the hole, plate_mount), so fixing the tangential displacement of
%   the node farthest from the hole removes the rotation at zero reaction.
%
%   ps returns the solution for inspection: u, v (ny x nx), the reactions at
%   the fixed nodes and the tension of each thread, and the node-averaged
%   stress resultants Nxx, Nyy, Nxy [N/m].

nx = dom.nx; ny = dom.ny;
N = nx * ny;
ps = struct('solved', false);
if P.g == 0 || (mnt.gx == 0 && mnt.gy == 0)
    Kg = sparse(N, N);
    return
end
idx = ops.idx; n = ops.n;
Bx = ops.Bx; By = ops.By; Bxc = ops.Bxc; Byc = ops.Byc; Cx = ops.Cx; Cy = ops.Cy;
wx = ops.wx; wy = ops.wy; area = ops.area; ncell = ops.ncell;
nex = numel(wx); ney = numel(wy);
E = P.material.E; nu = P.material.nu; h = P.plate.h;
C11 = E / (1 - nu^2); C12 = nu * C11; G = E / (2 * (1 + nu));

% Stiffness on [u; v]
Kuu = h * C11 * (Bx' * spdiags(wx, 0, nex, nex) * Bx);
Kvv = h * C11 * (By' * spdiags(wy, 0, ney, ney) * By);
Kuv = h * C12 * area * (Bxc' * Byc);
Gam = [Byc, Bxc];                                          % gamma = u_y + v_x at the cells
Kps = [Kuu, Kuv; Kuv', Kvv] + h * G * area * (Gam' * Gam);

% Load: the in-plane weight of each node
Mi = P.mu * ops.mu;
f = [Mi * P.g * mnt.gx; Mi * P.g * mnt.gy];

% Constraints, in the plate's node numbering
fix = mnt.fix;
loc = zeros(N, 1); loc(idx) = (1:n)';
fixn = loc(fix(:, 1));
if any(fixn == 0)
    error('plate_prestress: a constraint sits on a node off the plate');
end
fix_mount = unique(fixn);                                  % the mount's own constraints, for the reactions
if mnt.free_rotation
    % the tangential displacement of the node farthest from the pivot is fixed
    pk = fixn(1);
    X = dom.X(idx); Y = dom.Y(idx);
    d2 = (X - X(pk)).^2 + (Y - Y(pk)).^2;
    d2(ops.unk == 0) = -1;
    [~, far] = max(d2);
    t = [-(Y(far) - Y(pk)), X(far) - X(pk)];
    t = t / norm(t);
    fixn = [fixn; far];
    fix = [fix; idx(far) t];
end
% Rotation matrix T from the reduced unknowns to [u; v]
one = zeros(n, 2); cnt = zeros(n, 1);
for q = 1:numel(fixn)
    k = fixn(q);
    cnt(k) = cnt(k) + 1;
    one(k, :) = fix(q, 2:3) / norm(fix(q, 2:3));
end
if any(cnt > 2)
    error('plate_prestress: more than two constraints on one node');
end
both = cnt == 2;
single = cnt == 1;
none = cnt == 0;
kn = find(none); ks = find(single);
nq = 2 * numel(kn) + numel(ks);
ti = [kn; n + kn; ks; n + ks];
tj = [(1:numel(kn))'; numel(kn) + (1:numel(kn))'; 2 * numel(kn) + (1:numel(ks))'; 2 * numel(kn) + (1:numel(ks))'];
tv = [ones(2 * numel(kn), 1); -one(ks, 2); one(ks, 1)];    % the direction perpendicular to the fixed one
T = sparse(ti, tj, tv, 2 * n, nq);
Kred = T' * Kps * T;
fred = T' * f;
q = Kred \ fred;
U = T * q;
u = U(1:n); v = U(n + 1:end);
if any(~isfinite(U))
    error('plate_prestress: the in-plane solve failed; the mount does not hold the plate in its plane');
end
react = Kps * U - f;                                       % force the constraints apply, [x; y] per node

% Stress resultants times area at the edges and cells
exx = Bx * u; eyy = By * v;
exxc = Bxc * u; eyyc = Byc * v; gam = Gam * U;
nxx = h * (C11 * (wx .* exx) + C12 * area * (Cx' * eyyc));
nyy = h * (C11 * (wy .* eyy) + C12 * area * (Cy' * exxc));
nxy = h * G * area * gam;
Kgu = Bx' * spdiags(nxx, 0, nex, nex) * Bx + By' * spdiags(nyy, 0, ney, ney) * By ...
    + Bxc' * spdiags(nxy, 0, ncell, ncell) * Byc + Byc' * spdiags(nxy, 0, ncell, ncell) * Bxc;
Pz = spdiags(ops.unk, 0, n, n);
Kgu = Pz * Kgu * Pz;
[ki, kj, kv] = find(Kgu);
Kg = sparse(idx(ki), idx(kj), kv, N, N);

% Node-averaged resultants for display: each edge's value shared by its two
% nodes, each cell's by its four corners, divided by the areas they cover
ps.solved = true;
uf = zeros(ny, nx); vf = zeros(ny, nx);
uf(idx) = u; vf(idx) = v;
ps.u = uf; ps.v = vf;
Ex = abs(Bx) * (dom.dx / 2);                               % 1/2 at each node of an x-edge
Ey = abs(By) * (dom.dy / 2);
sx = Ex' * nxx; ax = Ex' * wx;
sy = Ey' * nyy; ay = Ey' * wy;
ic = (1:ncell)';
Qc = sparse([ic; ic; ic; ic], ops.corners(:), 0.25 * ones(4 * ncell, 1), ncell, n);
sxy = Qc' * nxy; ac = Qc' * (area * ones(ncell, 1));
Nxx = zeros(ny, nx); Nyy = zeros(ny, nx); Nxy = zeros(ny, nx);
Nxx(idx) = sx ./ max(ax, eps); Nyy(idx) = sy ./ max(ay, eps); Nxy(idx) = sxy ./ max(ac, eps);
Nxx(~dom.inside) = NaN; Nyy(~dom.inside) = NaN; Nxy(~dom.inside) = NaN;
ps.Nxx = Nxx; ps.Nyy = Nyy; ps.Nxy = Nxy;
ps.nxx = nxx; ps.nyy = nyy; ps.nxy = nxy;                  % the resultants times area, per edge and cell: the factors of Kg
% Reactions and thread tensions: the force along -ed (up) at each fixed node
rx = react(1:n); ry = react(n + 1:end);
ps.react_x = rx; ps.react_y = ry;
ps.up = -(rx * mnt.ed(1) + ry * mnt.ed(2));                % upward force each constrained node carries
ps.fixed_nodes = idx(fix_mount);
ps.mass_check = sum(ps.up(fix_mount)) / (mnt.mass * P.g);   % 1 when the constraints carry the weight
ps.energy = 0.5 * (U' * Kps * U);
end