2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 1function cmp = plate_compare(A, B)
2% PLATE_COMPARE Compare two solutions of the same plate problem.
3%
4% cmp = plate_compare(A, B) takes two result structs (for example from
5% plate_modal and plate_fd) and compares them over their common time window.
6% B is interpolated onto the time axis of A, so A should be the solution
7% with the coarser sampling in time if they differ.
8%
9% cmp.t, cmp.w_a, cmp.w_b, cmp.vm_a, cmp.vm_b aligned probe series
10% cmp.w_rel_l2, cmp.vm_rel_l2 ||a - b|| / ||a||
11% cmp.w_peak_a, cmp.w_peak_b, cmp.vm_peak_a, cmp.vm_peak_b
12% cmp.f_peak_a, cmp.f_peak_b frequency of the largest spectral peak
13% cmp.frames_rel_l2 per-frame relative field error of w, when
14% the two results share a grid and frame times
15% cmp.frames_t_mismatch largest difference between the frame times [s]
16%
17% Relative errors in stress are expected to exceed those in displacement,
18% since stress is a second derivative of the displacement.
20T = min(A.t(end), B.t(end));
21keep = A.t <= T;
22t = A.t(keep);
23w_a = A.probe.w(keep);
24vm_a = A.probe.vm(keep);
25w_b = interp1(B.t, B.probe.w, t);
26vm_b = interp1(B.t, B.probe.vm, t);
28cmp.t = t;
29cmp.w_a = w_a; cmp.w_b = w_b;
30cmp.vm_a = vm_a; cmp.vm_b = vm_b;
31cmp.w_rel_l2 = norm(w_a - w_b) / max(norm(w_a), eps);
32cmp.vm_rel_l2 = norm(vm_a - vm_b) / max(norm(vm_a), eps);
33cmp.w_peak_a = max(abs(w_a)); cmp.w_peak_b = max(abs(w_b));
34cmp.vm_peak_a = max(vm_a); cmp.vm_peak_b = max(vm_b);
36fmax = min(A.spectrum.f(end), B.spectrum.f(end));
37ia = A.spectrum.f <= fmax & A.spectrum.f > 0; ib = B.spectrum.f <= fmax & B.spectrum.f > 0; % skip the DC bin
38[~, ka] = max(A.spectrum.w(ia)); fa = A.spectrum.f(ia); cmp.f_peak_a = fa(ka);
39[~, kb] = max(B.spectrum.w(ib)); fb = B.spectrum.f(ib); cmp.f_peak_b = fb(kb);
41% Frames are compared directly when the grids match and the frame times agree to
42% within half a sample of each solver (each snaps its frames to its own time
43% samples, so they cannot coincide exactly).
44tol = 0.5 * ((A.t(2) - A.t(1)) + (B.t(2) - B.t(1))) + 1e-12;
45same_grid = isequal(size(A.frames.w), size(B.frames.w)) && ...
46 max(abs(A.grid.x - B.grid.x)) < 1e-9 && max(abs(A.grid.y - B.grid.y)) < 1e-9 && ...
47 max(abs(A.frames.t - B.frames.t)) <= tol;
48cmp.frames_t_mismatch = max(abs(A.frames.t(1:min(end, numel(B.frames.t))) - B.frames.t(1:min(end, numel(A.frames.t)))));
49if same_grid
50 nf = numel(A.frames.t);
51 cmp.frames_t = A.frames.t;
52 cmp.frames_rel_l2 = zeros(1, nf);
53 for k = 1:nf
54 da = A.frames.w(:, :, k); db = B.frames.w(:, :, k);
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 55 fin = isfinite(da) & isfinite(db); % nodes off an elliptical plate are NaN
56 cmp.frames_rel_l2(k) = norm(da(fin) - db(fin)) / max(norm(da(fin)), eps);
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 57 end
58 cmp.frames_diff = A.frames.w - B.frames.w;
59else
60 cmp.frames_t = [];
61 cmp.frames_rel_l2 = [];
62 cmp.frames_diff = [];
63end
64end