80c0b6dPlate dynamics lab: plates in realizable mounts under gravity, explicit finite differences on WebGPUJeremy Magland 1function result = plate_fd_output(S, Rec)
2% PLATE_FD_OUTPUT The result struct of plate_fd from the record of its time loop.
3%
4% result = plate_fd_output(S, Rec), with S from plate_fd_setup and Rec from
5% plate_fd_loop (or its WebGPU port), truncates the record after a blow-up,
6% evaluates the surface stresses at the saved frames and at the probe, sums
7% the listener's delay bins into the radiated pressure, takes the probe
8% spectra and assembles the fields listed in plate_result_schema. The
9% stresses are computed here, in double precision, from the flexural frames
10% the loop kept, whichever precision the loop ran in; the rigid-body part
11% has no curvature. Nodes outside the plate are NaN in every field, so
12% they are drawn blank.
14P = S.P;
15dom = S.dom; mnt = S.mnt;
16nx = S.nx; ny = S.ny; dx = S.dx; dy = S.dy; dt = S.dt; fs = S.fs;
17nr = S.nr; N = S.N;
18frame_steps = S.frame_steps;
19notes = [S.notes, Rec.notes];
20w_p = Rec.w_p; wxx_p = Rec.wxx_p; wyy_p = Rec.wyy_p; wxy_p = Rec.wxy_p;
21C = Rec.C; Fsup = Rec.Fsup; Wf = Rec.Wf; Cf = Rec.Cf; Pbins = Rec.Pbins;
23% Complete output samples only: the last one is partly filled, and after a
24% blow-up everything past the last finite step is dropped
25nsteps_done = S.nsteps;
26if Rec.diverged_step > 0
27 nsteps_done = Rec.diverged_step - 1;
28 keepf = frame_steps < Rec.diverged_step;
29 Wf = Wf(:, :, keepf); Cf = Cf(:, keepf); frame_steps = frame_steps(keepf);
30end
31n = floor(nsteps_done * dt * fs);
32n = max(1, min(n, S.na));
33t = ((0:n - 1) + 0.5) / fs;
34w_p = w_p(1:n); wxx_p = wxx_p(1:n); wyy_p = wyy_p(1:n); wxy_p = wxy_p(1:n);
35C = C(:, 1:n); Fsup = Fsup(:, 1:n);
36nf = numel(frame_steps);
37tf = frame_steps * dt;
39% Total departure at the frames, and the stresses from the flexural frames
40Wtot = zeros(ny, nx, nf);
41VMf = zeros(ny, nx, nf); SPf = zeros(ny, nx, nf);
42rect = strcmp(dom.shape, 'rectangle');
43for kf = 1:nf
44 phi = Wf(:, :, kf);
45 if nr > 0
46 Wtot(:, :, kf) = phi + reshape(S.R * Cf(:, kf), ny, nx);
47 else
48 Wtot(:, :, kf) = phi;
49 end
50 if rect
51 [wxx, wyy, wxy] = plate_curvature_free(phi, dx, dy, P.material.nu);
52 else
53 [wxx, wyy, wxy] = plate_curvature_domain(phi, dom);
54 end
55 Sf = plate_stress(P, wxx, wyy, wxy);
56 VMf(:, :, kf) = Sf.vm; SPf(:, :, kf) = Sf.sp;
57end
58if isempty(tf)
59 env_vm = zeros(ny, nx); env_sp = zeros(ny, nx);
60else
61 env_vm = max(VMf, [], 3); env_sp = max(abs(SPf), [], 3);
62end
63out = ~dom.inside;
64w0 = reshape(S.w0, ny, nx);
65for kf = 1:nf
66 tmp = Wf(:, :, kf); tmp(out) = NaN; Wf(:, :, kf) = tmp;
67 tmp = Wtot(:, :, kf); tmp(out) = NaN; Wtot(:, :, kf) = tmp;
68 tmp = VMf(:, :, kf); tmp(out) = NaN; VMf(:, :, kf) = tmp;
69 tmp = SPf(:, :, kf); tmp(out) = NaN; SPf(:, :, kf) = tmp;
70end
71env_vm(out) = NaN; env_sp(out) = NaN; w0(out) = NaN;
73% Audio: sum the delay bins, then trim to the complete samples
74na = S.na;
75p_audio = zeros(1, na);
76for d = 1:S.nbins
77 p_audio(d:end) = p_audio(d:end) + Pbins(d, 1:na - d + 1);
78end
79p_audio = p_audio(1:n);
81Sp = plate_stress(P, wxx_p, wyy_p, wxy_p);
82[fw, Aw] = plate_spectrum(t, w_p);
83[~, Avm] = plate_spectrum(t, Sp.vm);
85% Static support forces per support
86fstat_sup = zeros(1, S.nsup);
87for q = 1:S.nsup
88 nd = S.sup_nodes(S.sup_ptr(q):S.sup_ptr(q + 1) - 1);
89 if S.sup_kind(q) == 0
90 fstat_sup(q) = sum(S.fstat_s(nd));
91 elseif S.sup_kind(q) == 1
92 fstat_sup(q) = sum(S.fstat_b(nd));
93 else
94 fstat_sup(q) = sum(S.fstat_t(nd));
95 end
96end
98result.method = 'fd';
99result.grid = struct('x', S.x, 'y', S.y);
100result.t = t;
101result.probe = struct('x', S.x(S.jp), 'y', S.y(S.ip), 'w', w_p, 'sxx', Sp.sxx, 'syy', Sp.syy, 'txy', Sp.txy, 'vm', Sp.vm, 'sp', Sp.sp);
102result.frames = struct('t', tf, 'w', Wf, 'w_total', Wtot, 'rigid', Cf, 'vm', VMf, 'sp', SPf);
103result.envelope = struct('vm', env_vm, 'sp', env_sp);
104result.spectrum = struct('f', fw, 'w', Aw, 'vm', Avm);
105result.audio = struct('fs', fs, 'p', p_audio, 't0', S.Lis.t0, 'model', S.Lis.model);
106result.rigid = struct('nr', nr, 'names', {mnt.rigid_names}, 'series', C, ...
107 'pivot_x', mnt.pivot(1), 'pivot_y', mnt.pivot(2), 'ed', mnt.ed, 'eh', mnt.eh, ...
108 'hz', S.mdl.rigid_hz(:)', 'mode_names', {S.mdl.rigid_mode_names});
109result.supports = struct('n', S.nsup, 'names', {S.sup_names}, 'units', {S.sup_units}, 'kinds', {S.sup_kinds}, ...
110 'series', Fsup, 'static', fstat_sup);
111result.static = struct('w0', w0, 'sag', S.mdl.st.sag, 'gx', mnt.gx, 'gy', mnt.gy, 'gz', mnt.gz, 'g', P.g, ...
112 'mass', mnt.mass, 'xc', mnt.xc, 'yc', mnt.yc, 'held', double(mnt.held));
113if S.mdl.ps.solved
114 result.static.Nxx = S.mdl.ps.Nxx; result.static.Nyy = S.mdl.ps.Nyy; result.static.Nxy = S.mdl.ps.Nxy;
115end
116result.info = struct('nx', nx, 'ny', ny, 'dx', dx, 'dy', dy, 'dt', dt, 'dt_max', S.mdl.dt_max, 'cfl', P.fd.cfl, ...
117 'nsteps', S.nsteps, 'ncells', N, 'nunknowns', S.mdl.nunknowns, 'shape', dom.shape, 'mount', mnt.type, ...
118 'cell_steps', N * S.nsteps, 'fs', fs, ...
119 'stable', Rec.diverged_step == 0, 'diverged_step', Rec.diverged_step, ...
120 'loop', Rec.backend, 'loop_s', Rec.loop_s, 'wall_s', toc(S.t_start), 'notes', {notes});
121end