1function M = plate_fem_modes(params)
2% PLATE_FEM_MODES Modes of the hung plate by the Q4 Mindlin-Reissner finite-element method.
3%
4% M = plate_fem_modes(params) assembles the mesh (plate_fem_assemble) and
5% solves the generalized eigenproblem K v = omega^2 M v for its lowest
6% modes (plate_fem_eigs), enough of them to pass 0.47 fs. Three modes lie
7% below 20 Hz: the two rigid rotations about axes through the hole, which
8% nothing restrains, and the pendulum mode of the thread's spring. They
9% are kept, so that the collision can put momentum into the freely hanging
10% plate, and marked by their frequency; plate_synthesize leaves them out of
11% the sound. The audible modes are those between 20 Hz and 0.47 fs, at
12% most params.render.modes of them.
13%
14% Returns M with method 'fem', freqs (1 x nm, the low ones first), the
15% mass-normalized transverse shapes on the node grid, shapes ((ny + 1) x
16% (nx + 1) x nm, rows index y), the curvatures kxx, kyy, kxy of each shape
17% (from the rotations: w_xx = -d theta_x / dx and so on, element-centre
18% values averaged to the nodes), support_part (the shape at the hole),
19% volume_part (the trapezoidal integral of the shape over the plate),
20% grid_x, grid_y, the mesh sizes and the eigensolver's report.
22t_start = tic;
23P = plate_setup(params);
24A = plate_fem_assemble(params, P);
25nx = A.nx;
26ny = A.ny;
27nyn = ny + 1;
28nxn = nx + 1;
29notes = {};
31% how many modes to ask for: the modal-density estimate plus a margin, then
32% enlarge until the wanted set reaches past f_max or the cap
33plan = plate_plan_common(params, P);
34n_rigid = 3;
35lam_max = (2 * pi * P.f_max)^2;
36nwant = min(params.render.modes, plan.modes_expected + 12) + n_rigid;
37nwant = max(nwant, n_rigid + 4);
38tries = 0;
39while true
40 tries = tries + 1;
41 nb = min(A.ndof, nwant + max(nwant, 8));
42 [lam, V, einfo] = plate_fem_eigs(A.K, A.M, nwant, nb, 60, 1e-10);
43 n_below = nnz(lam < lam_max);
44 reached = lam(end) >= lam_max;
45 capped = n_below - n_rigid >= params.render.modes;
46 if reached || capped || nb >= A.ndof || tries >= 4
47 break
48 end
49 nwant = ceil(1.5 * nwant) + 4;
50end
51freqs = sqrt(max(lam, 0))' / (2 * pi);
52low = freqs < 20;
53audible = (lam' > 0) & (freqs >= 20) & (freqs < P.f_max);
54n_audible = nnz(audible);
55count = min(params.render.modes, n_audible);
56sel_a = find(audible);
57chosen = [find(low), sel_a(1:count)];
58freqs = freqs(chosen);
59V = V(:, chosen);
60nm = numel(freqs);
61% explicit M-normalization, as the reference simulator does
62for k = 1:nm
63 V(:, k) = V(:, k) / sqrt(V(:, k)' * (A.M * V(:, k)));
64end
65if ~einfo.converged
66 notes{end + 1} = sprintf('the eigensolver stopped after %d iterations without meeting its tolerance; largest residual %.1e', einfo.iterations, einfo.residual);
67end
68if ~reached && ~capped
69 notes{end + 1} = sprintf('the eigensolver found %d modes but did not reach 0.47 fs; some audible modes may be missing', n_audible);
70end
71if count < n_audible
72 notes{end + 1} = sprintf('%d modes lie below 0.47 fs; the first %d are kept (params.render.modes)', n_audible, count);
73end
74if nnz(low) ~= n_rigid
75 notes{end + 1} = sprintf('%d modes below 20 Hz where 3 (two free rotations and the pendulum) were expected', nnz(low));
76end
78% shapes on the node grid: node ey (nx + 1) + ex -> array (ey + 1, ex + 1)
79Wn = V(A.wdofs, :);
80shapes = permute(reshape(Wn, nxn, nyn, nm), [2, 1, 3]);
81RX = permute(reshape(V(2:3:end, :), nxn, nyn, nm), [2, 1, 3]);
82RY = permute(reshape(V(3:3:end, :), nxn, nyn, nm), [2, 1, 3]);
83dx = A.dx;
84dy = A.dy;
85% element-centre derivatives of the rotations, then the average of the elements around each node
86drx_dx = ((RX(1:end - 1, 2:end, :) + RX(2:end, 2:end, :)) - (RX(1:end - 1, 1:end - 1, :) + RX(2:end, 1:end - 1, :))) / (2 * dx);
87dry_dy = ((RY(2:end, 1:end - 1, :) + RY(2:end, 2:end, :)) - (RY(1:end - 1, 1:end - 1, :) + RY(1:end - 1, 2:end, :))) / (2 * dy);
88drx_dy = ((RX(2:end, 1:end - 1, :) + RX(2:end, 2:end, :)) - (RX(1:end - 1, 1:end - 1, :) + RX(1:end - 1, 2:end, :))) / (2 * dy);
89dry_dx = ((RY(1:end - 1, 2:end, :) + RY(2:end, 2:end, :)) - (RY(1:end - 1, 1:end - 1, :) + RY(2:end, 1:end - 1, :))) / (2 * dx);
90kxx = -plate_fem_nodal(drx_dx, nyn, nxn);
91kyy = -plate_fem_nodal(dry_dy, nyn, nxn);
92kxy = -plate_fem_nodal(drx_dy + dry_dx, nyn, nxn) / 2;
94% participations
95support_part = (V' * A.support)';
96area = ones(nyn, nxn) * dx * dy;
97area([1, end], :) = area([1, end], :) * 0.5;
98area(:, [1, end]) = area(:, [1, end]) * 0.5;
99volume_part = (reshape(shapes, nyn * nxn, nm)' * area(:))';
101M.method = 'fem';
102M.freqs = freqs;
103M.support_part = support_part;
104M.volume_part = volume_part;
105M.shapes = shapes;
106M.kxx = kxx;
107M.kyy = kyy;
108M.kxy = kxy;
109M.grid_x = A.grid_x;
110M.grid_y = A.grid_y;
111M.ndof = A.ndof;
112M.nn = A.nn;
113M.nx = nx;
114M.ny = ny;
115M.dx = dx;
116M.dy = dy;
117M.n_rigid = nnz(low);
118M.n_found = numel(lam);
119M.plate_mass = P.plate_mass;
120M.eig_iterations = einfo.iterations;
121M.eig_residual = einfo.residual;
122M.eig_block = einfo.block;
123M.size_note = sprintf('%d x %d Q4 elements, %d unknowns; %d Ritz vectors, %d subspace iterations', nx, ny, A.ndof, einfo.block, einfo.iterations);
124M.wall_s = toc(t_start);
125M.notes = notes;
126end