2004981Plate dynamics lab 2: a hammer on a hung plate by Rayleigh-Ritz and by Mindlin-Reissner finite elementsJeremy Magland 1function A = plate_fem_assemble(params, P)
2% PLATE_FEM_ASSEMBLE The global stiffness and mass of the meshed plate on its thread.
3%
4% A = plate_fem_assemble(params, P) meshes the plate into
5% params.fem.mesh_x by params.fem.mesh_y rectangular elements, sums the
6% element matrices of plate_fem_element into sparse K and M (three
7% unknowns per node, nodes numbered with x fastest: node = ey (nx + 1) + ex),
8% and adds the thread's pendulum spring m g / L to the transverse
9% displacement at the hole, interpolated within the element that holds it
10% (plate_fem_interp). The hole itself is not cut from the mesh: its
11% 0.8 mm diameter is below the element size at every mesh offered here,
12% and its mass is a part in a thousand of the plate's.
13%
14% Returns K, M, support (the interpolation vector of the hole as a sparse
15% column), ndof, nn (nodes), nx, ny, dx, dy, grid_x, grid_y (node
16% coordinates) and wdofs (the indices of the transverse unknowns in node
17% order).
19nx = params.fem.mesh_x;
20ny = params.fem.mesh_y;
21if nx < 4 || ny < 4
22 error('plate_fem_assemble: the mesh needs at least 4 elements per side');
23end
24W = params.plate.width;
25H = params.plate.height;
26dx = W / nx;
27dy = H / ny;
28[ke, me] = plate_fem_element(dx, dy, params, P);
30[EX, EY] = ndgrid(0:nx - 1, 0:ny - 1);
31ex = EX(:);
32ey = EY(:);
33nel = numel(ex);
34n1 = ey * (nx + 1) + ex;
35nodes = [n1, n1 + 1, n1 + nx + 2, n1 + nx + 1]; % anticlockwise, zero-based
36dofs = zeros(nel, 12);
37for a = 1:4
38 dofs(:, 3 * a - 2) = 3 * nodes(:, a) + 1;
39 dofs(:, 3 * a - 1) = 3 * nodes(:, a) + 2;
40 dofs(:, 3 * a) = 3 * nodes(:, a) + 3;
41end
42rows = kron(dofs, ones(1, 12));
43cols = repmat(dofs, 1, 12);
44kv = repmat(reshape(ke', 1, 144), nel, 1);
45mv = repmat(reshape(me', 1, 144), nel, 1);
46nn = (nx + 1) * (ny + 1);
47ndof = 3 * nn;
48K = sparse(rows(:), cols(:), kv(:), ndof, ndof);
49M = sparse(rows(:), cols(:), mv(:), ndof, ndof);
51% the thread at the hole
52[idx, wts] = plate_fem_interp(P.hole_x, P.hole_y, params);
53nyn = ny + 1;
54% node grid index (i + (j - 1) nyn, rows y) -> node number ey (nx + 1) + ex -> w unknown
55i_ = mod(idx - 1, nyn);
56j_ = floor((idx - 1) / nyn);
57node = i_ * (nx + 1) + j_;
58wd = 3 * node + 1;
59support = sparse(wd(:), ones(4, 1), wts(:), ndof, 1);
60K = K + P.thread_k * (support * support');
62A.K = K;
63A.M = M;
64A.support = support;
65A.ndof = ndof;
66A.nn = nn;
67A.nx = nx;
68A.ny = ny;
69A.dx = dx;
70A.dy = dy;
71A.grid_x = linspace(-W / 2, W / 2, nx + 1);
72A.grid_y = linspace(-H / 2, H / 2, ny + 1);
73A.wdofs = (1:3:ndof)';
74end