function st = plate_static(P, dom, mnt, A0, Mvec, unk)
% PLATE_STATIC  Static equilibrium of the plate in its mount under gravity.
%
%   st = plate_static(P, dom, mnt, A0, Mvec, unk) solves
%
%     (A0 + K_contact) w0 = M g_n + contact offsets
%
%   for the deflection w0 of the plate at rest, with A0 the linear operator
%   (bending, self-weight prestress and the linear springs of the mount,
%   plate_model), M the node masses, g_n the normal component of gravity, and
%   K_contact the springs of the one-sided contacts that are compressed at
%   rest.  Which contacts those are is found by an active-set iteration: all
%   contacts start engaged, any found in tension is released, and the solve is
%   repeated until the set no longer changes.  A plate that nothing holds in
%   equilibrium (too few pads to rest on, or one standing on its edge that
%   buckles under its own weight) fails the Cholesky factorization, and the
%   error says so.
%
%   The dynamics are then computed for the departure from w0 (plate_fd_loop):
%   the static forces cancel exactly, the contact forces are linear in the
%   departure until a contact opens, and single precision has only the motion
%   to represent.
%
%   Fields.
%     w0          N x 1 static deflection (zero at nodes that are not unknowns)
%     actb, actt  N x 1 logical: bottom and top contacts engaged at rest
%     gb, gt      N x 1 compression of each contact at rest [m]: the bottom
%                 contact pushes while gb - w > 0, the top while gt + w > 0
%     fstat_s, fstat_b, fstat_t   N x 1 force at rest of the linear springs, the
%                 bottom contacts and the top contacts [N]
%     Alin        the operator linearized about rest, for the release solve
%     sag         largest |w0| [m]
%     notes

N = dom.nx * dom.ny;
u = unk(:);
kb = mnt.kb(:); kt = mnt.kt(:); ob = mnt.ob(:); ot = mnt.ot(:);
ks = mnt.ks(:);
fg = Mvec(:) .* u * P.g * mnt.gz;
notes = {};
actb = kb > 0; actt = kt > 0;
w0 = zeros(N, 1);
iu = find(u);
if ~any(fg) && ~any(actb) && ~any(actt)
    Alin = A0;
else
    for iter = 1:20
        kc = kb .* actb + kt .* actt;
        Alin = A0 + spdiags(kc, 0, N, N);
        rhs = fg + kb .* actb .* ob - kt .* actt .* ot;
        w0 = zeros(N, 1);
        w0(iu) = plate_static_solve(Alin(iu, iu), rhs(iu));
        nb = (kb > 0) & (ob - w0 > 0);
        nt = (kt > 0) & (ot + w0 > 0);
        if isequal(nb, actb) && isequal(nt, actt)
            break
        end
        actb = nb; actt = nt;
        if iter == 20
            notes{end + 1} = 'the set of engaged contacts did not settle in 20 iterations; the last one is used';
        end
    end
    if any(kb > 0) && ~all(actb(kb > 0))
        notes{end + 1} = sprintf('%d support node(s) carry no load at rest: the plate does not touch them', nnz(kb > 0 & ~actb));
    end
    if any(kt > 0) && ~all(actt(kt > 0))
        notes{end + 1} = sprintf('%d upper-felt node(s) are clear of the plate at rest', nnz(kt > 0 & ~actt));
    end
end
st.w0 = w0;
st.actb = actb; st.actt = actt;
st.gb = ob - w0; st.gt = ot + w0;
st.fstat_s = -ks .* w0;                        % linear springs
st.fstat_b = kb .* max(st.gb, 0);              % bottom contacts, pushing up
st.fstat_t = -kt .* max(st.gt, 0);             % top contacts, pushing down
st.Alin = Alin;
st.sag = max(abs(w0));
st.notes = notes;
end

function x = plate_static_solve(A, b)
[Rc, p] = chol(A);
if p > 0
    error('plate_static: the plate has no stable equilibrium in this mount (it tips, swings over or buckles under its own weight); check the pads, the holes or the tilt');
end
x = Rc \ (Rc' \ b);
end