function A = plate_fem_assemble(params, P)
% PLATE_FEM_ASSEMBLE  The global stiffness and mass of the meshed plate on its thread.
%
%   A = plate_fem_assemble(params, P) meshes the plate into
%   params.fem.mesh_x by params.fem.mesh_y rectangular elements, sums the
%   element matrices of plate_fem_element into sparse K and M (three
%   unknowns per node, nodes numbered with x fastest: node = ey (nx + 1) + ex),
%   and adds the thread's pendulum spring m g / L to the transverse
%   displacement at the hole, interpolated within the element that holds it
%   (plate_fem_interp).  The hole itself is not cut from the mesh: its
%   0.8 mm diameter is below the element size at every mesh offered here,
%   and its mass is a part in a thousand of the plate's.
%
%   Returns K, M, support (the interpolation vector of the hole as a sparse
%   column), ndof, nn (nodes), nx, ny, dx, dy, grid_x, grid_y (node
%   coordinates) and wdofs (the indices of the transverse unknowns in node
%   order).

nx = params.fem.mesh_x;
ny = params.fem.mesh_y;
if nx < 4 || ny < 4
    error('plate_fem_assemble: the mesh needs at least 4 elements per side');
end
W = params.plate.width;
H = params.plate.height;
dx = W / nx;
dy = H / ny;
[ke, me] = plate_fem_element(dx, dy, params, P);

[EX, EY] = ndgrid(0:nx - 1, 0:ny - 1);
ex = EX(:);
ey = EY(:);
nel = numel(ex);
n1 = ey * (nx + 1) + ex;
nodes = [n1, n1 + 1, n1 + nx + 2, n1 + nx + 1];           % anticlockwise, zero-based
dofs = zeros(nel, 12);
for a = 1:4
    dofs(:, 3 * a - 2) = 3 * nodes(:, a) + 1;
    dofs(:, 3 * a - 1) = 3 * nodes(:, a) + 2;
    dofs(:, 3 * a) = 3 * nodes(:, a) + 3;
end
rows = kron(dofs, ones(1, 12));
cols = repmat(dofs, 1, 12);
kv = repmat(reshape(ke', 1, 144), nel, 1);
mv = repmat(reshape(me', 1, 144), nel, 1);
nn = (nx + 1) * (ny + 1);
ndof = 3 * nn;
K = sparse(rows(:), cols(:), kv(:), ndof, ndof);
M = sparse(rows(:), cols(:), mv(:), ndof, ndof);

% the thread at the hole
[idx, wts] = plate_fem_interp(P.hole_x, P.hole_y, params);
nyn = ny + 1;
% node grid index (i + (j - 1) nyn, rows y) -> node number ey (nx + 1) + ex -> w unknown
i_ = mod(idx - 1, nyn);
j_ = floor((idx - 1) / nyn);
node = i_ * (nx + 1) + j_;
wd = 3 * node + 1;
support = sparse(wd(:), ones(4, 1), wts(:), ndof, 1);
K = K + P.thread_k * (support * support');

A.K = K;
A.M = M;
A.support = support;
A.ndof = ndof;
A.nn = nn;
A.nx = nx;
A.ny = ny;
A.dx = dx;
A.dy = dy;
A.grid_x = linspace(-W / 2, W / 2, nx + 1);
A.grid_y = linspace(-H / 2, H / 2, ny + 1);
A.wdofs = (1:3:ndof)';
end