2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 1function result = plate_implicit(params)
2% PLATE_IMPLICIT Implicit finite differences (theta scheme) for the damped Kirchhoff-Love plate.
3%
4% result = plate_implicit(params) advances
5%
6% w_tt + 2 sigma0 w_t - 2 sigma1 lap(w_t) + kappa^2 lap(lap(w)) = F(t) s(x,y) / mu
7%
8% on the same grid and with the same spatial stencils as plate_fd, but the
9% stiffness term is averaged over three time levels,
10%
11% (w+ - 2w + w-)/dt^2 + sigma0 (w+ - w-)/dt - sigma1 L (w+ - w-)/dt
12% + kappa^2 B (theta w+ + (1 - 2 theta) w + theta w-) = F s / mu,
13%
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 14% with L and B the sparse Laplacian and biharmonic matrices of plate_matrices
15% on the plate's footprint (plate_domain), which also supplies the diagonal
16% mass matrix M: the identity on the interior nodes for simply supported or
17% clamped edges, the trapezoidal node weights on every node for free edges.
18% Multiplying by dt^2 and collecting terms gives one linear system per step
19% with a constant matrix,
78d40e9Free edges for both finite-difference methodsJeremy Magland 21% A w+ = C w - A_old w- + dt^2 F M s / mu
22% A = (1 + sigma0 dt) M - sigma1 dt L + theta kappa^2 dt^2 B
23% C = 2 M - (1 - 2 theta) kappa^2 dt^2 B
24% A_old = (1 - sigma0 dt) M + sigma1 dt L + theta kappa^2 dt^2 B.
25%
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 26% For free edges the rigid-body motions are projected out of w+ after every
27% solve (plate_rigid_modes); for the centre-pinned plate the pinned node is
28% zeroed first.
30% A is symmetric positive definite, so it is factored once by sparse Cholesky
31% and each step costs two sparse triangular solves plus two sparse products.
32% (R' is formed once; forming it every step would cost more than the solves.)
33% For theta >= 1/4 the scheme is unconditionally stable (von Neumann analysis
34% of the amplification polynomial); theta = 1/4 is the trapezoidal rule. The
35% step is therefore chosen for accuracy, dt = 1 / (fs steps_per_sample).
36%
37% Approximations. Second order in space and time. Implicit averaging
38% lengthens the period of a mode by about (omega dt)^2 / 12, so the highest
39% modes are detuned rather than unstable. Probe, frames, stress envelope and
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 40% audio follow the conventions of plate_fd (see plate_result_schema); on an
41% ellipse the edge is a staircase and first order, and the nodes outside the
42% plate are NaN in the fields.
44t_start = tic;
45P = plate_setup(params);
46plan = plate_implicit_plan(params);
47nx = plan.nx; ny = plan.ny; dx = plan.dx; dy = plan.dy; dt = plan.dt; nsteps = plan.nsteps;
48theta = P.implicit.theta;
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 49dom = plate_domain(P.plate, nx, ny);
50x = dom.x; y = dom.y; X = dom.X; Y = dom.Y;
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 51fs = P.output.fs;
52notes = {};
53if theta < 0.25
54 notes{end + 1} = sprintf('theta = %g is below 1/4, so the scheme is only conditionally stable', theta);
55end
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 56rect = strcmp(dom.shape, 'rectangle');
78d40e9Free edges for both finite-difference methodsJeremy Magland 58% ---- Matrices on the unknowns, factored once
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 59free = any(strcmp(P.bc, {'free', 'pinned-center'}));
60[L, B, M, idx, fixed] = plate_matrices(dom, P.bc, P.material.nu);
62mvec = full(diag(M));
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 63unknown = false(ny, nx); unknown(idx) = true;
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 64k2dt2 = P.kappa^2 * dt^2;
78d40e9Free edges for both finite-difference methodsJeremy Magland 65A = (1 + P.damping.sigma0 * dt) * M - P.damping.sigma1 * dt * L + theta * k2dt2 * B;
66C = 2 * M - (1 - 2 * theta) * k2dt2 * B;
67A_old = (1 - P.damping.sigma0 * dt) * M + P.damping.sigma1 * dt * L + theta * k2dt2 * B;
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 68t_factor = tic;
69R = chol(A); % A = R' * R
70Rt = R'; % transposed once, not once per step
71factor_s = toc(t_factor);
73% ---- Load, probe, frames, listener (as in plate_fd)
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 74s = plate_force_footprint(P.excitation, X, Y, dx * dy, unknown);
75load = (dt^2 / P.mu) * (mvec .* s(idx)); % M s on the unknowns
76if free && isempty(fixed)
77 [Rb, Qb] = plate_rigid_modes(X(idx), Y(idx), mvec);
78d40e9Free edges for both finite-difference methodsJeremy Magland 78 notes{end + 1} = 'free edges: the rigid-body translation and tilts are removed from the displacement after every step';
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 79elseif free
80 [Rb, Qb] = plate_rigid_modes(X(idx), Y(idx), mvec, [X(idx(fixed)) Y(idx(fixed))]);
81 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(idx(fixed)), Y(idx(fixed)));
82end
83if ~rect
84 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));
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 86tt = (0:nsteps) * dt;
87F = plate_force_time(P.excitation, tt);
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 89% Probe: the nearest node whose eight neighbours belong to the plate
90d2 = (X - P.probe.x).^2 + (Y - P.probe.y).^2;
91d2(~dom.full3) = Inf;
92[~, kp] = min(d2(:));
93ip = mod(kp - 1, ny) + 1;
94jp = floor((kp - 1) / ny) + 1;
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 95if abs(x(jp) - P.probe.x) > 0.5 * dx || abs(y(ip) - P.probe.y) > 0.5 * dy
96 notes{end + 1} = sprintf('probe moved to the nearest interior node (%.4g, %.4g) m', x(jp), y(ip));
97end
99frame_steps = unique(max(1, min(nsteps, round(linspace(0, P.implicit.duration, P.output.nframes) / dt))));
100nf = numel(frame_steps);
101Wf = zeros(ny, nx, nf); VMf = zeros(ny, nx, nf); SPf = zeros(ny, nx, nf);
102tf = zeros(1, nf);
104Lis = plate_listener(P.listener, X, Y, dx * dy, fs);
105bins = Lis.delay(:) + 1;
106wgt = Lis.weight(:);
107na = ceil(nsteps * dt * fs) + Lis.nbins + 1;
108Pbins = zeros(Lis.nbins, na);
78d40e9Free edges for both finite-difference methodsJeremy Magland 110% ---- Time stepping on the vector of unknowns
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 111w = zeros(m, 1); wm = zeros(m, 1);
112wfull = zeros(ny, nx);
113w_p = zeros(1, nsteps + 1); wxx_p = zeros(1, nsteps + 1); wyy_p = zeros(1, nsteps + 1); wxy_p = zeros(1, nsteps + 1);
114blowup = 1e3 * P.plate.h;
115diverged_step = 0;
116kf = 1;
118for nstep = 1:nsteps
119 rhs = C * w - A_old * wm + F(nstep) * load;
120 wn = R \ (Rt \ rhs); % forward then back substitution
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 122 if ~isempty(fixed)
123 wn(fixed) = 0;
124 end
126 end
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 128 wfull(idx) = wn;
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 129 w_p(nstep + 1) = wfull(ip, jp);
130 wxx_p(nstep + 1) = (wfull(ip, jp + 1) - 2 * wfull(ip, jp) + wfull(ip, jp - 1)) / dx^2;
131 wyy_p(nstep + 1) = (wfull(ip + 1, jp) - 2 * wfull(ip, jp) + wfull(ip - 1, jp)) / dy^2;
132 wxy_p(nstep + 1) = (wfull(ip + 1, jp + 1) - wfull(ip + 1, jp - 1) - wfull(ip - 1, jp + 1) + wfull(ip - 1, jp - 1)) / (4 * dx * dy);
134 acc = zeros(ny, nx);
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 135 acc(idx) = (wn - 2 * w + wm) / dt^2;
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 136 ia = floor(nstep * dt * fs) + 1;
137 Pbins(:, ia) = Pbins(:, ia) + accumarray(bins, wgt .* acc(:), [Lis.nbins 1]) * (dt * fs);
139 if kf <= nf && nstep == frame_steps(kf)
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 140 if ~rect
141 [wxx, wyy, wxy] = plate_curvature_domain(wfull, dom);
142 elseif free
78d40e9Free edges for both finite-difference methodsJeremy Magland 143 [wxx, wyy, wxy] = plate_curvature_free(wfull, dx, dy, P.material.nu);
144 else
145 [wxx, wyy, wxy] = plate_curvature(plate_pad(wfull, P.bc), dx, dy);
146 end
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 147 S = plate_stress(P, wxx, wyy, wxy);
148 Wf(:, :, kf) = wfull; VMf(:, :, kf) = S.vm; SPf(:, :, kf) = S.sp; tf(kf) = nstep * dt;
149 kf = kf + 1;
150 end
152 peak = max(abs(wn));
153 if ~(peak < blowup)
154 diverged_step = nstep;
155 notes{end + 1} = sprintf('solution diverged at step %d (t = %.3g s); results are truncated there', nstep, nstep * dt);
156 break
157 end
158 wm = w; w = wn;
159end
161if diverged_step > 0
162 nkeep = diverged_step;
163 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);
164 keepf = frame_steps < diverged_step;
165 Wf = Wf(:, :, keepf); VMf = VMf(:, :, keepf); SPf = SPf(:, :, keepf); tf = tf(keepf);
166end
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 167nf = numel(tf);
168if isempty(tf)
169 env_vm = zeros(ny, nx); env_sp = zeros(ny, nx);
170else
171 env_vm = max(VMf, [], 3); env_sp = max(abs(SPf), [], 3);
172end
173if ~rect
174 out = ~dom.inside;
175 for kf = 1:nf
176 tmp = Wf(:, :, kf); tmp(out) = NaN; Wf(:, :, kf) = tmp;
177 tmp = VMf(:, :, kf); tmp(out) = NaN; VMf(:, :, kf) = tmp;
178 tmp = SPf(:, :, kf); tmp(out) = NaN; SPf(:, :, kf) = tmp;
179 end
180 env_vm(out) = NaN; env_sp(out) = NaN;
181end
183p_audio = zeros(1, na);
184for d = 1:Lis.nbins
185 p_audio(d:end) = p_audio(d:end) + Pbins(d, 1:na - d + 1);
186end
187n_audio = min(na, ceil(numel(tt) * dt * fs));
188p_audio = p_audio(1:n_audio);
190Sp = plate_stress(P, wxx_p, wyy_p, wxy_p);
191[fw, Aw] = plate_spectrum(tt, w_p);
192[~, Avm] = plate_spectrum(tt, Sp.vm);
194result.method = 'implicit';
195result.grid = struct('x', x, 'y', y);
196result.t = tt;
197result.probe = struct('x', x(jp), 'y', y(ip), 'w', w_p, 'sxx', Sp.sxx, 'syy', Sp.syy, 'txy', Sp.txy, 'vm', Sp.vm, 'sp', Sp.sp);
198result.frames = struct('t', tf, 'w', Wf, 'vm', VMf, 'sp', SPf);
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 199result.envelope = struct('vm', env_vm, 'sp', env_sp);
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 200result.spectrum = struct('f', fw, 'w', Aw, 'vm', Avm);
201result.audio = struct('fs', fs, 'p', p_audio, 't0', Lis.t0);
202result.info = struct('nx', nx, 'ny', ny, 'dx', dx, 'dy', dy, 'dt', dt, 'theta', theta, ...
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 203 'steps_per_sample', P.implicit.steps_per_sample, 'nunknowns', m, 'shape', dom.shape, 'nnz_factor', nnz(R), 'factor_s', factor_s, ...
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 204 'nsteps', nsteps, 'ncells', nx * ny, 'cell_steps', plan.cell_steps, 'fs', fs, ...
205 'stable', diverged_step == 0, 'diverged_step', diverged_step, 'wall_s', toc(t_start), 'notes', {notes});
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 206if ~isempty(fixed)
207 result.info.support_x = X(idx(fixed));
208 result.info.support_y = Y(idx(fixed));
209end