1function N = plate_fem_nodal(E, nyn, nxn)
2% PLATE_FEM_NODAL Element-centre values averaged to the nodes.
3%
4% N = plate_fem_nodal(E, nyn, nxn) takes values on the ny x nx elements
5% (any number of layers along the third dimension) and returns the mean of
6% the elements around each node of the (nyn = ny + 1) x (nxn = nx + 1) grid:
7% four inside, two along an edge, one at a corner.
9nl = size(E, 3);
10acc = zeros(nyn, nxn, nl);
11cnt = zeros(nyn, nxn);
12acc(1:end - 1, 1:end - 1, :) = acc(1:end - 1, 1:end - 1, :) + E;
13acc(2:end, 1:end - 1, :) = acc(2:end, 1:end - 1, :) + E;
14acc(1:end - 1, 2:end, :) = acc(1:end - 1, 2:end, :) + E;
15acc(2:end, 2:end, :) = acc(2:end, 2:end, :) + E;
16cnt(1:end - 1, 1:end - 1) = cnt(1:end - 1, 1:end - 1) + 1;
17cnt(2:end, 1:end - 1) = cnt(2:end, 1:end - 1) + 1;
18cnt(1:end - 1, 2:end) = cnt(1:end - 1, 2:end) + 1;
19cnt(2:end, 2:end) = cnt(2:end, 2:end) + 1;
20N = acc ./ repmat(cnt, [1, 1, nl]);
21end