function N = plate_fem_nodal(E, nyn, nxn)
% PLATE_FEM_NODAL  Element-centre values averaged to the nodes.
%
%   N = plate_fem_nodal(E, nyn, nxn) takes values on the ny x nx elements
%   (any number of layers along the third dimension) and returns the mean of
%   the elements around each node of the (nyn = ny + 1) x (nxn = nx + 1) grid:
%   four inside, two along an edge, one at a corner.

nl = size(E, 3);
acc = zeros(nyn, nxn, nl);
cnt = zeros(nyn, nxn);
acc(1:end - 1, 1:end - 1, :) = acc(1:end - 1, 1:end - 1, :) + E;
acc(2:end, 1:end - 1, :) = acc(2:end, 1:end - 1, :) + E;
acc(1:end - 1, 2:end, :) = acc(1:end - 1, 2:end, :) + E;
acc(2:end, 2:end, :) = acc(2:end, 2:end, :) + E;
cnt(1:end - 1, 1:end - 1) = cnt(1:end - 1, 1:end - 1) + 1;
cnt(2:end, 1:end - 1) = cnt(2:end, 1:end - 1) + 1;
cnt(1:end - 1, 2:end) = cnt(1:end - 1, 2:end) + 1;
cnt(2:end, 2:end) = cnt(2:end, 2:end) + 1;
N = acc ./ repmat(cnt, [1, 1, nl]);
end