1function S = plate_fd_setup(params)
2% PLATE_FD_SETUP Everything the explicit finite-difference time loop needs, computed once.
3%
4% S = plate_fd_setup(params) does the part of plate_fd that comes before the
5% time stepping: the derived constants (plate_setup), the grid and time step
6% (plate_fd_plan), the plate's footprint on the grid (plate_domain), the load
7% footprint and force history, the probe node, the frame schedule, the
8% listener geometry and, whenever the loop needs them, the sparse operators
9% and the rigid-body projector (plate_matrices, plate_rigid_modes).
10% plate_fd_loop advances the plate with S, and plate_fd_output turns what the
11% loop recorded into the result struct; plate_fd calls the three in turn.
12%
13% The split lets a host run the loop somewhere else while the setup and the
14% output stay this code: the lab has a WebGPU port of plate_fd_loop that reads
15% exactly the fields below and hands its record to plate_fd_output.
16%
17% Two forms of the spatial operators. A rectangle with supported or clamped
18% edges uses the dense stencils on the whole ny x nx array (plate_pad,
19% plate_laplacian). Everything else, free edges on any shape and every
20% condition on an ellipse, uses the sparse operators of plate_matrices on the
21% unknowns, embedded in vectors over all N = nx ny nodes so that the loop
22% keeps one array layout: rows and columns of nodes that are not unknowns are
23% zero, their node weight is one, and S.keep zeroes them after every step.
24%
25% Fields.
26% dom the footprint (plate_domain)
27% nx, ny, dx, dy, x, y, X, Y grid; w(i, j) lives at (x(j), y(i))
28% dt, nsteps, tt, F, fs time step, step count, tt = (0:nsteps) dt,
29% force F(tt), audio sample rate
30% s load footprint, ny x nx [1/m^2], zero off the plate
31% ip, jp row and column of the probe node
32% frame_steps, nf steps at which the field is kept
33% Lis, bins, wgt, na plate_listener geometry; delay bin and
34% weight of every node in column order; number
35% of audio samples to accumulate
36% free the edges are free (with or without the pin)
37% nunknowns nodes in the system (the pinned node included)
38% sparse the loop uses Lm, Km, mvec, keep, Rb, Qb
39% Lm, Km, mvec N x N Laplacian and stiffness, N node weights
40% keep ny x nx, 1 at the unknowns and 0 elsewhere
41% (off the plate, on a supported ring, at the pin)
42% nr, Rb, Qb rigid-body motions to remove (0, 2 or 3):
43% N x nr basis and projector, zero off the plate
44% fixed node index of the pinned node, or empty
45% c_old, c_new, k2dt2, s1dt2, constants of the update
46% fdt2, blowup wn = c_new (2 w - c_old wm - k2dt2 B
47% + s1dt2 (Lw - Lw_m) + fdt2 F(n) s)
48% with c_old = 1 - sigma0 dt, c_new =
49% 1 / (1 + sigma0 dt), k2dt2 = kappa^2 dt^2,
50% s1dt2 = 2 sigma1 dt, fdt2 = dt^2 / mu; blowup
51% is the deflection that declares the run
52% unstable
53% notes remarks about this particular run
54% t_start tic at entry, for the wall time
56S.t_start = tic;
57P = plate_setup(params);
58plan = plate_fd_plan(params);
59nx = plan.nx; ny = plan.ny; dx = plan.dx; dy = plan.dy; dt = plan.dt; nsteps = plan.nsteps;
60dom = plate_domain(P.plate, nx, ny);
61x = dom.x; y = dom.y; X = dom.X; Y = dom.Y; % w(i, j) lives at (x(j), y(i))
62N = nx * ny;
63fs = P.output.fs;
64notes = {};
65if P.fd.cfl > 1
66 notes{end + 1} = sprintf('cfl = %g exceeds the stability limit; expect the solution to blow up', P.fd.cfl);
67end
68rect = strcmp(dom.shape, 'rectangle');
69free = any(strcmp(P.bc, {'free', 'pinned-center'}));
70use_sparse = free || ~rect;
71if free
72 unknown = dom.inside;
73else
74 unknown = dom.interior;
75end
77% Load on the unknowns only: a supported ring and the outside carry none
78s = plate_force_footprint(P.excitation, X, Y, dx * dy, unknown);
79tt = (0:nsteps) * dt;
80F = plate_force_time(P.excitation, tt);
82% Probe: the nearest node whose eight neighbours belong to the plate, so the
83% curvatures there are central differences
84d2 = (X - P.probe.x).^2 + (Y - P.probe.y).^2;
85d2(~dom.full3) = Inf;
86[~, kp] = min(d2(:));
87ip = mod(kp - 1, ny) + 1;
88jp = floor((kp - 1) / ny) + 1;
89if abs(x(jp) - P.probe.x) > 0.5 * dx || abs(y(ip) - P.probe.y) > 0.5 * dy
90 notes{end + 1} = sprintf('probe moved to the nearest interior node (%.4g, %.4g) m', x(jp), y(ip));
91end
93% Frames at the steps nearest to nframes equally spaced times in [0, duration]
94frame_steps = unique(max(1, min(nsteps, round(linspace(0, P.fd.duration, P.output.nframes) / dt))));
96% Listener: every grid node is a patch; the acceleration of a step is added into
97% the audio sample that contains it, delayed by the node's travel time.
98Lis = plate_listener(P.listener, X, Y, dx * dy, fs);
100% Sparse operators on the unknowns, embedded over all N nodes, and the
101% projector onto the complement of the rigid-body motions
102fixed = zeros(0, 1);
103Lm = []; Km = []; mvec = []; Rb = []; Qb = []; nr = 0;
104keep = double(unknown);
105nunknowns = nnz(unknown); % the pinned node counts: it is in the system, held at zero
106if use_sparse
107 [Lu, Ku, Mu, idx, fixed_u] = plate_matrices(dom, P.bc, P.material.nu);
108 [ki, kj, kv] = find(Ku);
109 Km = sparse(idx(ki), idx(kj), kv, N, N);
110 [li, lj, lv] = find(Lu);
111 Lm = sparse(idx(li), idx(lj), lv, N, N);
112 mvec = ones(N, 1);
113 mvec(idx) = full(diag(Mu));
114 if ~isempty(fixed_u)
115 fixed = idx(fixed_u);
116 keep(fixed) = 0;
117 end
118 if free
119 if isempty(fixed)
120 [Ru, Qu] = plate_rigid_modes(X(idx), Y(idx), mvec(idx));
121 notes{end + 1} = 'free edges: the rigid-body translation and tilts are removed from the displacement after every step';
122 else
123 [Ru, Qu] = plate_rigid_modes(X(idx), Y(idx), mvec(idx), [X(fixed) Y(fixed)]);
124 notes{end + 1} = sprintf('pinned at the node (%.4g, %.4g) m nearest the centre; the two rigid rotations about it are removed after every step', X(fixed), Y(fixed));
125 end
126 nr = size(Ru, 2);
127 Rb = zeros(N, nr); Rb(idx, :) = Ru;
128 Qb = zeros(N, nr); Qb(idx, :) = Qu;
129 end
130end
131if ~rect
132 notes{end + 1} = sprintf('%s on a staircase of %d grid nodes; the edge is first order in the grid spacing', dom.shape, nnz(dom.inside));
133end
135S.P = P; S.plan = plan; S.dom = dom;
136S.nx = nx; S.ny = ny; S.dx = dx; S.dy = dy; S.x = x; S.y = y; S.X = X; S.Y = Y;
137S.dt = dt; S.nsteps = nsteps; S.tt = tt; S.F = F; S.fs = fs;
138S.s = s;
139S.ip = ip; S.jp = jp;
140S.frame_steps = frame_steps; S.nf = numel(frame_steps);
141S.Lis = Lis; S.bins = Lis.delay(:) + 1; S.wgt = Lis.weight(:);
142S.na = ceil(nsteps * dt * fs) + Lis.nbins + 1;
143S.free = free; S.sparse = use_sparse; S.fixed = fixed; S.nunknowns = nunknowns;
144S.Lm = Lm; S.Km = Km; S.mvec = mvec; S.keep = keep; S.nr = nr; S.Rb = Rb; S.Qb = Qb;
145S.c_old = 1 - P.damping.sigma0 * dt;
146S.c_new = 1 / (1 + P.damping.sigma0 * dt);
147S.k2dt2 = P.kappa^2 * dt^2;
148S.s1dt2 = 2 * P.damping.sigma1 * dt;
149S.fdt2 = dt^2 / P.mu;
150S.blowup = 1e3 * P.plate.h; % deflection beyond which the run is declared unstable
151S.notes = notes;
152end