Blame packages/core/src/m/plate_fem_strike.m at main - proof-of-concept/plate-dynamics-lab-2

proof-of-concept / plate-dynamics-lab-2

2004981 Plate dynamics lab 2: a hammer on a hung plate by Rayleigh-Ritz and by Mindlin-Reissner finite elementsJeremy Magland 1 month ago1function result = plate_fem_strike(params, M)

2% PLATE_FEM_STRIKE Strike the finite-element modes and render the result.

4% result = plate_fem_strike(params, M) takes the modes of plate_fem_modes,

5% interpolates each shape at the strike point and at the probe within the

6% element that holds them (plate_fem_interp), sets the modal damping

7% (plate_damping; the rigid-body modes get 0.01), integrates the coupled

8% hammer/plate contact (plate_fem_contact) and synthesizes the sound, the

9% series and the frames (plate_synthesize). Only this function needs to

10% run again when the hammer, the strike point, the microphone or the

11% render settings change.

13t_start = tic;

14if ~strcmp(M.method, 'fem')

15 error('plate_fem_strike: the modes were computed by another method (%s)', M.method);

16end

17P = plate_setup(params);

18pl = params.plate;

19if abs(params.strike.x) > pl.width / 2 || abs(params.strike.y) > pl.height / 2

20 error('plate_fem_strike: the strike point is outside the plate');

21end

22if abs(params.probe.x) > pl.width / 2 || abs(params.probe.y) > pl.height / 2

23 error('plate_fem_strike: the probe is outside the plate');

24end

25if params.fem.mesh_x ~= M.nx || params.fem.mesh_y ~= M.ny

26 error('plate_fem_strike: the modes were computed on a %d x %d mesh, not %d x %d', M.nx, M.ny, params.fem.mesh_x, params.fem.mesh_y);

27end

29nm = numel(M.freqs);

30nn = numel(M.shapes) / nm;

31Sh = reshape(M.shapes, nn, nm);

32[idx, wts] = plate_fem_interp(params.strike.x, params.strike.y, params);

33phi_hit = wts * Sh(idx, :);

34[zeta, Q] = plate_damping(M.freqs, M.support_part, params, P);

36fs = params.render.fs;

37na = round(params.render.duration * fs);

38t = (0:na - 1) / fs;

39C = plate_fem_contact(params, P, M.freqs, zeta, phi_hit, t);

41[ip, wp] = plate_fem_interp(params.probe.x, params.probe.y, params);

42D.zeta = zeta;

43D.Q = Q;

44D.phi_hit = phi_hit;

45D.phi_probe = wp * Sh(ip, :);

46Kxx = reshape(M.kxx, nn, nm);

47Kyy = reshape(M.kyy, nn, nm);

48Kxy = reshape(M.kxy, nn, nm);

49D.kxx_probe = wp * Kxx(ip, :);

50D.kyy_probe = wp * Kyy(ip, :);

51D.kxy_probe = wp * Kxy(ip, :);

53result = plate_synthesize(params, P, M, C, D);

54result.info.strike_wall_s = toc(t_start);

55end