1function S = plate_fd_setup(params)
2% PLATE_FD_SETUP Everything the explicit time loop needs, computed once.
3%
4% S = plate_fd_setup(params) takes the model (plate_model: operators,
5% statics, rigid-body basis, time step) and adds what a run needs: the load
6% footprint and force history, the probe node, the frame schedule, the
7% listener geometry, the initial state and the constants of the update.
8% plate_fd_loop advances the plate with S, and plate_fd_output turns what
9% the loop recorded into the result struct; plate_fd calls the three in
10% turn. The split lets a host run the loop somewhere else while the setup
11% and the output stay this code: the lab has a WebGPU port of plate_fd_loop
12% that reads exactly the fields below and hands its record to plate_fd_output.
13%
14% State. The loop carries the departure from the static equilibrium w0 as
15% rigid-body coordinates c (nr of them, plate_rigid_basis) and a flexural
16% field phi on the grid; the displacement is w = w0 + R c + phi. A release
17% (excitation.type = 'release') starts from the static deflection under a
18% force at the pull point, scaled to the asked displacement there and split
19% into c0 and phi0; a strike starts from rest.
20%
21% Fields the loop reads (indices are 1-based, arrays over all N = nx ny nodes
22% in column-major order unless said otherwise).
23% nx, ny, dx, dy, dt, nsteps, fs, na grid, step, count, sample rate and
24% number of output samples
25% keep N x 1: 1 at the unknowns, 0 elsewhere
26% Mvec N x 1 node masses (1 where keep is 0)
27% KB, KG, Lm sparse N x N: bending, prestress plus
28% linear springs, Laplacian per node area
29% Gs, Gscale, WB, WG, WL, ng, inva the same operators in factored form
30% (see below), for the port
31% nr, R, Q, AR, Krr, Ginv rigid-body basis (N x nr), projector,
32% coupling forces, reduced stiffness and
33% inverse mass (nr x nr)
34% kb, cb, gb, kt, ct, gt bottom and top contacts: stiffness
35% [N/m], dashpot [N s/m], compression at
36% rest [m], per node
37% cs, ks linear dashpots and springs (the springs
38% for the force report; they are in KG)
39% fstat_s, fstat_b, fstat_t static force of the springs, the bottom
40% and the top contacts per node
41% sup_ptr, sup_nodes, nsup, sup_kind supports as a CSR list of node indices,
42% each reporting one kind of force
43% (0 spring, 1 bottom, 2 top)
44% sa force per node per unit F(t) [N/N]
45% F force history at the nsteps + 1 levels
46% ip, jp probe row and column
47% frame_steps, nf steps at which the field is kept
48% lis_node, lis_bin, lis_wgt, nbins listener entries (plate_listener)
49% c0, phi0 initial state
50% wsamp na x 1 weight of a step in its sample
51% c_old, c_new, s1f, LR, dt2 constants of the update: the sigma0 loss,
52% the sigma1 loss as a nodal force on the
53% total velocity (with Lm R for the rigid part)
54% blowup deflection that declares the run unstable
55% notes, t_start
57S.t_start = tic;
58mdl = plate_model(params);
59P = mdl.P; dom = mdl.dom; mnt = mdl.mnt;
60nx = mdl.nx; ny = mdl.ny; dx = mdl.dx; dy = mdl.dy; N = mdl.N;
61dt = mdl.dt; nsteps = mdl.nsteps;
62X = dom.X; Y = dom.Y; x = dom.x; y = dom.y;
63unk = mdl.unk;
64fs = P.output.fs;
65notes = mdl.notes;
67% Load on the unknowns: force per node per unit of F(t)
68s = plate_force_footprint(P.excitation, X, Y, dx * dy, unk);
69sa = s(:) .* mdl.a(:);
70tt = (0:nsteps) * dt;
71F = plate_force_time(P.excitation, tt);
73% Probe: the nearest node whose eight neighbours belong to the plate
74d2 = (X - P.probe.x).^2 + (Y - P.probe.y).^2;
75d2(~dom.full3 | mnt.held) = Inf;
76[~, kp] = min(d2(:));
77ip = mod(kp - 1, ny) + 1;
78jp = floor((kp - 1) / ny) + 1;
79if abs(x(jp) - P.probe.x) > 0.5 * dx || abs(y(ip) - P.probe.y) > 0.5 * dy
80 notes{end + 1} = sprintf('probe moved to the nearest interior node (%.4g, %.4g) m', x(jp), y(ip));
81end
83% Frames at the steps nearest to nframes equally spaced times in (0, duration]
84frame_steps = unique(max(1, min(nsteps, round(linspace(0, P.fd.duration, P.output.nframes) / dt))));
86% Listener, and the output samples: each series value is the mean over the
87% steps whose time falls in the sample, so the weight of a step is one over
88% the number of steps its sample contains
89Lis = plate_listener(P.listener, dom, dx * dy, fs);
90na = floor(nsteps * dt * fs) + 1;
91ia_of = floor((0:nsteps - 1) * dt * fs) + 1;
92cnt = accumarray(ia_of', 1, [na 1]);
93wsamp = 1 ./ max(cnt, 1);
95% Initial state
96nr = mdl.nr;
97c0 = zeros(nr, 1);
98phi0 = zeros(N, 1);
99if strcmp(P.excitation.type, 'release')
100 e = P.excitation;
101 sp = plate_force_footprint(e, X, Y, dx * dy, unk);
102 f = sp(:) .* mdl.a(:);
103 A = mdl.st.Alin;
104 if mnt.free_rotation && nr > 0
105 % the spin about the thread has no stiffness to speak of: take the
106 % spin moment out of the load and pin the spin at one node
107 rh = mdl.R(:, 3);
108 Mr = mdl.Mvec .* rh;
109 f = f - (rh' * f) / (rh' * Mr) * Mr;
110 [~, kf] = max(abs(rh));
111 A = A + sparse(kf, kf, mdl.Mvec(kf) * 1e6, N, N);
112 notes{end + 1} = 'the pull is applied without twisting the plate about its thread';
113 end
114 iu = find(unk(:));
115 wr = zeros(N, 1);
116 [Rc, p] = chol(A(iu, iu));
117 if p > 0
118 error('plate_fd_setup: no static deflection to release from; the plate is not held stably');
119 end
120 wr(iu) = Rc \ (Rc' \ f(iu));
121 wP = (sp(:) .* mdl.a(:))' * wr; % displacement at the pull point
122 if ~(abs(wP) > 0)
123 error('plate_fd_setup: the pull point does not move; is it on a held part of the plate?');
124 end
125 wr = wr * (e.displacement / wP);
126 if nr > 0
127 c0 = mdl.Q' * wr;
128 if mnt.free_rotation
129 c0(3) = 0;
130 end
131 phi0 = wr - mdl.R * c0;
132 else
133 phi0 = wr;
134 end
135 phi0 = phi0 .* double(unk(:));
136 notes{end + 1} = sprintf('released from a static pull of %.3g mm at (%.4g, %.4g) m, which needed %.3g N', e.displacement * 1e3, e.x, e.y, e.displacement / wP);
137end
139% The operators in factored form, for the WebGPU port. Every operator of the
140% loop is a sum of terms G_q' W G_q with G_q a difference stencil, and the
141% port applies them that way: strains e = G phi first (each an exact integer
142% stencil times a scale), then the weights, then G' back to the nodes. The
143% assembled matrices lose their smooth modes in single precision (the 13-point
144% biharmonic of a smooth field is the small residual of large cancelling
145% terms, and the rounding of the field alone is a few percent of it at low
146% wavenumbers); the factored form keeps them to about 1e-5. Gs stacks
147% Dxx, Dyy, Dxc, Dyc, Dxy, Bx, By, Bxc, Byc; Gscale is the scale of each row
148% (the stencil divided by it is integer); WB, WG, WL are the symmetric weight
149% matrices with KB = Gs' WB Gs, KG - Ks = Gs' WG Gs and Lm = diag(1/a) Gs' WL Gs
150% on the unknowns.
151ops = mdl.ops;
152ni = ops.n; idx = ops.idx;
153Gu = [ops.Dxx; ops.Dyy; ops.Dxc; ops.Dyc; ops.Dxy; ops.Bx; ops.By; ops.Bxc; ops.Byc];
154nrow = [size(ops.Dxx, 1), size(ops.Dyy, 1), size(ops.Dxc, 1), size(ops.Dyc, 1), size(ops.Dxy, 1), ...
155 size(ops.Bx, 1), size(ops.By, 1), size(ops.Bxc, 1), size(ops.Byc, 1)];
156scl = [1 / dx^2, 1 / dy^2, 1 / dx^2, 1 / dy^2, 1 / (dx * dy), 1 / dx, 1 / dy, 1 / (2 * dx), 1 / (2 * dy)];
157Gscale = zeros(sum(nrow), 1);
158off = cumsum([0 nrow]);
159for q = 1:9
160 Gscale(off(q) + 1:off(q + 1)) = scl(q);
161end
162ng = sum(nrow);
163[gi, gj, gv] = find(Gu);
164Gs = sparse(gi, idx(gj), gv, ng, N);
165D = ops.D; nu = ops.nu; area = ops.area;
166o = off;
167WB = sparse(ng, ng);
168WB = WB + sparse(o(1) + (1:nrow(1)), o(1) + (1:nrow(1)), D * ops.mkx, ng, ng);
169WB = WB + sparse(o(2) + (1:nrow(2)), o(2) + (1:nrow(2)), D * ops.mky, ng, ng);
170WB = WB + sparse(o(3) + (1:nrow(3)), o(4) + (1:nrow(4)), D * nu * ops.mkc, ng, ng);
171WB = WB + sparse(o(4) + (1:nrow(4)), o(3) + (1:nrow(3)), D * nu * ops.mkc, ng, ng);
172WB = WB + sparse(o(5) + (1:nrow(5)), o(5) + (1:nrow(5)), D * 2 * (1 - nu) * area * ones(nrow(5), 1), ng, ng);
173WG = sparse(ng, ng);
174if mdl.ps.solved
175 WG = WG + sparse(o(6) + (1:nrow(6)), o(6) + (1:nrow(6)), mdl.ps.nxx, ng, ng);
176 WG = WG + sparse(o(7) + (1:nrow(7)), o(7) + (1:nrow(7)), mdl.ps.nyy, ng, ng);
177 WG = WG + sparse(o(8) + (1:nrow(8)), o(9) + (1:nrow(9)), mdl.ps.nxy, ng, ng);
178 WG = WG + sparse(o(9) + (1:nrow(9)), o(8) + (1:nrow(8)), mdl.ps.nxy, ng, ng);
179end
180WL = sparse(o(6) + (1:nrow(6)), o(6) + (1:nrow(6)), -ops.wx, ng, ng) + sparse(o(7) + (1:nrow(7)), o(7) + (1:nrow(7)), -ops.wy, ng, ng);
182% Supports as a CSR list
183sup = mnt.supports;
184nsup = numel(sup);
185sup_ptr = ones(nsup + 1, 1);
186sup_nodes = zeros(0, 1);
187for k = 1:nsup
188 nd = sup(k).nodes(:);
189 sup_nodes = [sup_nodes; nd];
190 sup_ptr(k + 1) = sup_ptr(k) + numel(nd);
191end
192sup_names = cell(1, nsup); sup_units = cell(1, nsup); sup_kinds = cell(1, nsup);
193sup_kind = zeros(nsup, 1); % 0 spring, 1 bottom contact, 2 top contact
194for k = 1:nsup
195 sup_names{k} = sup(k).name; sup_units{k} = sup(k).unit; sup_kinds{k} = sup(k).kind;
196 switch sup(k).kind
197 case 'spring'
198 sup_kind(k) = 0;
199 case 'bottom'
200 sup_kind(k) = 1;
201 case 'top'
202 sup_kind(k) = 2;
203 otherwise
204 error('plate_fd_setup: unknown support kind "%s"', sup(k).kind);
205 end
206end
208S.mdl = mdl; S.P = P; S.dom = dom; S.mnt = mnt;
209S.nx = nx; S.ny = ny; S.dx = dx; S.dy = dy; S.x = x; S.y = y; S.X = X; S.Y = Y; S.N = N;
210S.dt = dt; S.nsteps = nsteps; S.tt = tt; S.F = F; S.fs = fs; S.na = na; S.wsamp = wsamp;
211S.keep = double(unk(:));
212S.Mvec = mdl.Mvec; S.a = mdl.a(:);
213S.KB = mdl.Kb; S.KG = mdl.Kg + mdl.Ks; S.Lm = mdl.Lm;
214S.Gs = Gs; S.Gscale = Gscale; S.WB = WB; S.WG = WG; S.WL = WL; S.ng = ng;
215S.inva = double(unk(:)) ./ max(mdl.a(:), eps);
216S.nr = nr; S.R = mdl.R; S.Q = mdl.Q; S.AR = mdl.AR; S.Krr = mdl.Krr; S.Ginv = mdl.Ginv;
217S.kb = mnt.kb(:); S.cb = mnt.cb(:); S.gb = mdl.st.gb;
218S.kt = mnt.kt(:); S.ct = mnt.ct(:); S.gt = mdl.st.gt;
219S.cs = mnt.cs(:); S.ks = mnt.ks(:);
220S.fstat_s = mdl.st.fstat_s; S.fstat_b = mdl.st.fstat_b; S.fstat_t = mdl.st.fstat_t;
221S.w0 = mdl.st.w0;
222S.sup_ptr = sup_ptr; S.sup_nodes = sup_nodes; S.nsup = nsup; S.sup_kind = sup_kind;
223S.sup_names = sup_names; S.sup_units = sup_units; S.sup_kinds = sup_kinds;
224S.sa = sa; S.s = s;
225S.ip = ip; S.jp = jp;
226S.frame_steps = frame_steps; S.nf = numel(frame_steps);
227S.Lis = Lis; S.lis_node = Lis.node; S.lis_bin = Lis.delay + 1; S.lis_wgt = Lis.weight; S.nbins = Lis.nbins;
228S.c0 = c0; S.phi0 = phi0;
229S.c_old = 1 - P.damping.sigma0 * dt;
230S.c_new = 1 / (1 + P.damping.sigma0 * dt);
231S.s1f = 2 * P.damping.sigma1 * mdl.Mvec / dt; % the sigma1 loss as a nodal force per unit Lm(w - wm)
232S.LR = mdl.Lm * mdl.R; % Laplacian per node of the rigid basis, N x nr
233S.dt2 = dt^2;
234S.blowup = 1e3 * P.plate.h;
235S.notes = notes;
236end