80c0b6dPlate dynamics lab: plates in realizable mounts under gravity, explicit finite differences on WebGPUJeremy Magland 1function [R, Q, G] = plate_rigid_basis(dom, mnt, Mvec, unk)
2% PLATE_RIGID_BASIS Rigid-body coordinates of the plate about the mount's pivot.
3%
4% [R, Q, G] = plate_rigid_basis(dom, mnt, Mvec, unk) returns the N x nr
5% basis R of the rigid-body motions the mount leaves free, evaluated at the
6% grid nodes and zero at the nodes that are not unknowns (unk false), the
7% projector Q = M R (R' M R)^-1 with M = diag(Mvec) the node masses, so that
8%
9% c = Q' w are the rigid-body coordinates of a field w and
10% w - R c its flexural remainder, mass-orthogonal to R,
11%
12% and the nr x nr rigid mass matrix G = R' M R. For nr = 3 the columns are
13%
14% 1 translation along the normal (heave, or the
15% sway of a hanging plate's hole, in metres)
16% (P - pivot) . ed rotation about the horizontal in-plane axis eh
17% through the pivot, in radians (the swing of a
18% hung plate, the tilt of one on pads)
19% (P - pivot) . eh rotation about the in-plane axis ed through the
20% pivot (the spin of a hung plate about its thread)
21%
22% with mnt.pivot, mnt.ed (in-plane down) and mnt.eh (in-plane horizontal)
23% from plate_mount. A mount that holds nodes still (the vise) has nr = 0.
24%
25% The time stepping carries the two parts separately (plate_fd_loop): the
26% rigid coordinates as three scalar equations with the exact reduced
27% stiffness R' A R, the remainder on the grid. Bending annihilates R, so
28% the split costs nothing in accuracy and it is what keeps the metre-scale
29% sway and the micrometre-scale flexure apart in single precision.
31N = dom.nx * dom.ny;
32nr = mnt.nr;
33if nr == 0
34 R = zeros(N, 0); Q = zeros(N, 0); G = zeros(0, 0);
35 return
36end
37X = dom.X(:); Y = dom.Y(:);
38rd = (X - mnt.pivot(1)) * mnt.ed(1) + (Y - mnt.pivot(2)) * mnt.ed(2);
39rh = (X - mnt.pivot(1)) * mnt.eh(1) + (Y - mnt.pivot(2)) * mnt.eh(2);
40R = [ones(N, 1), rd, rh];
41R(~unk(:), :) = 0;
42m = Mvec(:) .* double(unk(:));
43G = R' * (m .* R);
44Q = (m .* R) / G;
45end