1function [idx, wts, ex, ey] = plate_fem_interp(x, y, params)
2% PLATE_FEM_INTERP Where a point of the plate falls in the mesh, and its shape weights.
3%
4% [idx, wts, ex, ey] = plate_fem_interp(x, y, params) finds the element
5% (ex, ey, zero-based) containing the point (x, y) (metres from the
6% plate's centre; a point on the far edge belongs to the last element) and
7% returns the bilinear weights wts (1 x 4) of its four nodes, with idx
8% (1 x 4) the linear index of each node in the node grid stored as an
9% (ny + 1) x (nx + 1) array with rows indexing y. A field sampled at the
10% nodes is interpolated as wts * field(idx).
12nx = params.fem.mesh_x;
13ny = params.fem.mesh_y;
14W = params.plate.width;
15H = params.plate.height;
16dx = W / nx;
17dy = H / ny;
18ux = min(max((x + W / 2) / dx, 0), nx - 1e-12);
19uy = min(max((y + H / 2) / dy, 0), ny - 1e-12);
20ex = min(floor(ux), nx - 1);
21ey = min(floor(uy), ny - 1);
22xi = 2 * (ux - ex) - 1;
23eta = 2 * (uy - ey) - 1;
24wts = plate_fem_shape(xi, eta);
25nyn = ny + 1;
26% nodes anticlockwise from (ex, ey): rows are y, so the linear index is (ey + 1) + ex * nyn
27idx = [(ey + 1) + ex * nyn, (ey + 1) + (ex + 1) * nyn, (ey + 2) + (ex + 1) * nyn, (ey + 2) + ex * nyn];
28end