1function [R, Q] = plate_rigid_modes(X, Y, mvec)
2% PLATE_RIGID_MODES Basis of the rigid-body motions of a free plate and its M-orthogonal projector.
3%
4% [R, Q] = plate_rigid_modes(X, Y, mvec) takes the node coordinates X, Y and
5% the diagonal mvec of the mass matrix (the trapezoidal weights of
6% plate_matrices) and returns R = [1, x - xc, y - yc] evaluated at the nodes
7% (n x 3) and Q = M R (R' M R)^-1 (n x 3), so that
8%
9% w - R * (Q' * w)
10%
11% is the part of w orthogonal to the rigid-body motions in the mass inner
12% product: what remains when the mean displacement and the two mean tilts are
13% taken out.
14%
15% Why remove them. A completely free plate has three zero-frequency modes,
16% and a strike gives it net momentum and angular momentum: with only the loss
17% sigma0 to stop it, the plate drifts to a displacement of order
18% (impulse / mass) / (2 sigma0), millimetres against the micrometres of the
19% bending motion. A real free-free measurement rests the plate on foam or
20% hangs it from threads, which absorb exactly this motion and nothing else,
21% so the solvers project it out after every step. The bending stiffness K
22% annihilates R exactly, so this changes nothing about the flexural dynamics;
23% the sigma1 loss operator does not, and the projection also keeps that term
24% from feeding a rigid tilt back into flexure.
26x = X(:);
27y = Y(:);
28m = mvec(:);
29xc = sum(m .* x) / sum(m);
30yc = sum(m .* y) / sum(m);
31R = [ones(numel(x), 1), x - xc, y - yc];
32G = R' * (m .* R); % 3 x 3 Gram matrix in the mass inner product
33Q = (m .* R) / G;
34end