1function result = plate_fd_output(S, R)
2% PLATE_FD_OUTPUT The result struct of plate_fd from the record of its time loop.
3%
4% result = plate_fd_output(S, R), with S from plate_fd_setup and R 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 displacement
10% frames the loop kept, whichever precision the loop ran in. On an ellipse
11% the nodes outside the plate are set to NaN in every field, so they are
12% drawn blank.
14P = S.P;
15dom = S.dom;
16nx = S.nx; ny = S.ny; dx = S.dx; dy = S.dy; dt = S.dt; fs = S.fs;
17tt = S.tt;
18frame_steps = S.frame_steps;
19w_p = R.w_p; wxx_p = R.wxx_p; wyy_p = R.wyy_p; wxy_p = R.wxy_p;
20Wf = R.Wf;
21notes = [S.notes, R.notes];
22rect = strcmp(dom.shape, 'rectangle');
24% Truncate after a blow-up so the outputs stay finite
25if R.diverged_step > 0
26 nkeep = R.diverged_step; % keep the last finite level
27 tt = tt(1:nkeep); w_p = w_p(1:nkeep); wxx_p = wxx_p(1:nkeep); wyy_p = wyy_p(1:nkeep); wxy_p = wxy_p(1:nkeep);
28 keepf = frame_steps < R.diverged_step;
29 Wf = Wf(:, :, keepf); frame_steps = frame_steps(keepf);
30end
31nf = numel(frame_steps);
32tf = frame_steps * dt;
34% Stress at the saved frames, from the curvatures of each frame
35VMf = zeros(ny, nx, nf); SPf = zeros(ny, nx, nf);
36for kf = 1:nf
37 if ~rect
38 [wxx, wyy, wxy] = plate_curvature_domain(Wf(:, :, kf), dom);
39 elseif S.free
40 [wxx, wyy, wxy] = plate_curvature_free(Wf(:, :, kf), dx, dy, P.material.nu);
41 else
42 [wxx, wyy, wxy] = plate_curvature(plate_pad(Wf(:, :, kf), P.bc), dx, dy);
43 end
44 Sf = plate_stress(P, wxx, wyy, wxy);
45 VMf(:, :, kf) = Sf.vm; SPf(:, :, kf) = Sf.sp;
46end
47if isempty(tf)
48 env_vm = zeros(ny, nx); env_sp = zeros(ny, nx);
49else
50 env_vm = max(VMf, [], 3); env_sp = max(abs(SPf), [], 3);
51end
52if ~rect
53 % Off the plate there is nothing to show
54 out = ~dom.inside;
55 for kf = 1:nf
56 tmp = Wf(:, :, kf); tmp(out) = NaN; Wf(:, :, kf) = tmp;
57 tmp = VMf(:, :, kf); tmp(out) = NaN; VMf(:, :, kf) = tmp;
58 tmp = SPf(:, :, kf); tmp(out) = NaN; SPf(:, :, kf) = tmp;
59 end
60 env_vm(out) = NaN; env_sp(out) = NaN;
61end
63% Audio: sum the delay bins, then trim to the simulated duration
64na = S.na;
65p_audio = zeros(1, na);
66for d = 1:S.Lis.nbins
67 p_audio(d:end) = p_audio(d:end) + R.Pbins(d, 1:na - d + 1);
68end
69n_audio = min(na, ceil(numel(tt) * dt * fs));
70p_audio = p_audio(1:n_audio);
72Sp = plate_stress(P, wxx_p, wyy_p, wxy_p);
73[fw, Aw] = plate_spectrum(tt, w_p);
74[~, Avm] = plate_spectrum(tt, Sp.vm);
76result.method = 'fd';
77result.grid = struct('x', S.x, 'y', S.y);
78result.t = tt;
79result.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);
80result.frames = struct('t', tf, 'w', Wf, 'vm', VMf, 'sp', SPf);
81result.envelope = struct('vm', env_vm, 'sp', env_sp);
82result.spectrum = struct('f', fw, 'w', Aw, 'vm', Avm);
83result.audio = struct('fs', fs, 'p', p_audio, 't0', S.Lis.t0);
84result.info = struct('nx', nx, 'ny', ny, 'dx', dx, 'dy', dy, 'dt', dt, 'dt_max', S.plan.dt_max, 'cfl', P.fd.cfl, ...
85 'nsteps', S.nsteps, 'ncells', nx * ny, 'nunknowns', S.nunknowns, 'shape', dom.shape, 'cell_steps', S.plan.cell_steps, 'fs', fs, ...
86 'stable', R.diverged_step == 0, 'diverged_step', R.diverged_step, ...
87 'loop', R.backend, 'loop_s', R.loop_s, 'wall_s', toc(S.t_start), 'notes', {notes});
88if ~isempty(S.fixed)
89 result.info.support_x = S.X(S.fixed);
90 result.info.support_y = S.Y(S.fixed);
91end
92end