1function [K, L, ops] = plate_bending(dom, held, nu, D)
2% PLATE_BENDING Sparse bending and Laplacian operators on the plate by the energy method.
3%
4% [K, L, ops] = plate_bending(dom, held, nu, D) assembles, on every node of
5% the plate (dom.inside from plate_domain), the stiffness matrix K of the
6% bending energy
7%
8% U = D/2 * integral( w_xx^2 + w_yy^2 + 2 nu w_xx w_yy + 2 (1-nu) w_xy^2 ) dA
9%
10% and the operator L of the Laplacian energy integral |grad w|^2 dA, both as
11% sparse N x N matrices over all nx ny grid nodes (rows and columns of nodes
12% off the plate are empty). Second differences are taken at the nodes that
13% have both neighbours on the plate, the cross difference at the centre of
14% each cell whose four corners are on the plate, each weighted by its area
15% (the finite-difference energy method; Bushnell, Computers & Structures
16% 1973). Inside the plate K is D times the 13-point biharmonic stencil
17% times the cell area; on the edges the free-edge conditions, no bending
18% moment and no Kirchhoff shear, are the natural conditions of the discrete
19% energy, met to first order in the spacing. K is symmetric positive
20% semidefinite with the three rigid-body motions 1, x, y as its null space.
21%
22% Held nodes (the logical ny x nx array held, the jaws of a vise) are part
23% of the energy but not of the unknowns: their rows and columns are dropped
24% after assembly. A second difference centred on a held node that reaches a
25% free node then penalizes the slope at the clamp, which is exactly the
26% clamped condition w = w_n = 0; it needs a held node behind every held edge
27% node, which plate_mount guarantees.
28%
29% L is the Laplacian energy operator with the natural condition of zero
30% normal slope, each grid edge weighted by its share of cells; the loss
31% term of the plate equation uses it divided by the node areas.
32%
33% ops returns the difference operators the prestress term shares
34% (plate_prestress): Bx, By first differences along the x and y edges of the
35% plate (edges x areas), Bxc, Byc the cell-centred gradients, and the edge and
36% cell weights, all on the numbering ops.idx of the plate's nodes.
38nx = dom.nx; ny = dom.ny; dx = dom.dx; dy = dom.dy;
39N = nx * ny;
40inside = dom.inside;
41idx = find(inside);
42n = numel(idx);
43me = (1:n)';
44unk = inside & ~held;
46% Node number of every grid node (0 off the plate), padded so every node has four neighbours to look up
47num_in = zeros(ny, nx);
48num_in(idx) = me;
49num = zeros(ny + 2, nx + 2);
50num(2:end - 1, 2:end - 1) = num_in;
51nbL = num(2:end - 1, 1:end - 2); nbR = num(2:end - 1, 3:end);
52nbD = num(1:end - 2, 2:end - 1); nbU = num(3:end, 2:end - 1);
53uL = nbL(idx); uR = nbR(idx); uD = nbD(idx); uU = nbU(idx);
55area = dx * dy;
56mu = dom.m(idx) * area; % node areas
57hasX = uL > 0 & uR > 0;
58hasY = uD > 0 & uU > 0;
59kx = find(hasX); nkx = numel(kx);
60ky = find(hasY); nky = numel(ky);
61Dxx = sparse([(1:nkx)'; (1:nkx)'; (1:nkx)'], [uL(kx); me(kx); uR(kx)], [ones(nkx, 1); -2 * ones(nkx, 1); ones(nkx, 1)] / dx^2, nkx, n);
62Dyy = sparse([(1:nky)'; (1:nky)'; (1:nky)'], [uD(ky); me(ky); uU(ky)], [ones(nky, 1); -2 * ones(nky, 1); ones(nky, 1)] / dy^2, nky, n);
63Kxx = Dxx' * spdiags(mu(kx), 0, nkx, nkx) * Dxx;
64Kyy = Dyy' * spdiags(mu(ky), 0, nky, nky) * Dyy;
65kc = find(hasX & hasY); nkc = numel(kc);
66Dxc = sparse([(1:nkc)'; (1:nkc)'; (1:nkc)'], [uL(kc); me(kc); uR(kc)], [ones(nkc, 1); -2 * ones(nkc, 1); ones(nkc, 1)] / dx^2, nkc, n);
67Dyc = sparse([(1:nkc)'; (1:nkc)'; (1:nkc)'], [uD(kc); me(kc); uU(kc)], [ones(nkc, 1); -2 * ones(nkc, 1); ones(nkc, 1)] / dy^2, nkc, n);
68Wc = spdiags(mu(kc), 0, nkc, nkc);
69Kc = nu * (Dxc' * Wc * Dyc + Dyc' * Wc * Dxc);
71% Cells with four corners on the plate: the twist at the centre
72A00 = num_in(1:end - 1, 1:end - 1); A10 = num_in(2:end, 1:end - 1);
73A01 = num_in(1:end - 1, 2:end); A11 = num_in(2:end, 2:end);
74cu = A00 > 0 & A10 > 0 & A01 > 0 & A11 > 0;
75kcell = find(cu); ncell = numel(kcell);
76c00 = A00(kcell); c10 = A10(kcell); c01 = A01(kcell); c11 = A11(kcell);
77ic = (1:ncell)';
78Dxy = sparse([ic; ic; ic; ic], [c00; c10; c01; c11], [ones(ncell, 1); -ones(ncell, 1); -ones(ncell, 1); ones(ncell, 1)] / (dx * dy), ncell, n);
79Kxy = 2 * (1 - nu) * area * (Dxy' * Dxy);
80Ku = D * (Kxx + Kyy + Kc + Kxy);
82% Laplacian energy: first differences along the grid edges, weighted by the
83% cell area each edge shares (a whole cell inside, half a cell on the boundary)
84cpad = zeros(ny + 1, nx + 1);
85cpad(2:end - 1, 2:end - 1) = cu;
86WX = (cpad(2:end, 2:end) + cpad(1:end - 1, 2:end)) / 2; % x-edge from (i, j) to (i, j+1)
87WY = (cpad(2:end, 2:end) + cpad(2:end, 1:end - 1)) / 2; % y-edge from (i, j) to (i+1, j)
88ex = find(uR > 0); nex = numel(ex);
89ey = find(uU > 0); ney = numel(ey);
90Bx = sparse([(1:nex)'; (1:nex)'], [me(ex); uR(ex)], [-ones(nex, 1); ones(nex, 1)] / dx, nex, n);
91By = sparse([(1:ney)'; (1:ney)'], [me(ey); uU(ey)], [-ones(ney, 1); ones(ney, 1)] / dy, ney, n);
92wx = WX(idx(ex)) * area; wy = WY(idx(ey)) * area;
93Lu = -(Bx' * spdiags(wx, 0, nex, nex) * Bx + By' * spdiags(wy, 0, ney, ney) * By);
95% Cell-centred gradients: the average of the two x-edges (bottom and top) and the two y-edges (left and right)
96xe_of = zeros(n, 1); xe_of(me(ex)) = (1:nex)';
97ye_of = zeros(n, 1); ye_of(me(ey)) = (1:ney)';
98Cx = sparse([ic; ic], [xe_of(c00); xe_of(c10)], 0.5 * ones(2 * ncell, 1), ncell, nex);
99Cy = sparse([ic; ic], [ye_of(c00); ye_of(c01)], 0.5 * ones(2 * ncell, 1), ncell, ney);
100Bxc = Cx * Bx;
101Byc = Cy * By;
103% Drop the held nodes, then embed over all N grid nodes
104sel = double(unk(idx));
105Pz = spdiags(sel, 0, n, n);
106Ku = Pz * Ku * Pz;
107Lu = Pz * Lu * Pz;
108[ki, kj, kv] = find(Ku);
109K = sparse(idx(ki), idx(kj), kv, N, N);
110[li, lj, lv] = find(Lu);
111L = sparse(idx(li), idx(lj), lv, N, N);
113ops.idx = idx; ops.n = n; ops.mu = mu;
114ops.Bx = Bx; ops.By = By; ops.Bxc = Bxc; ops.Byc = Byc; ops.Cx = Cx; ops.Cy = Cy;
115ops.wx = wx; ops.wy = wy; ops.ncell = ncell; ops.area = area;
116ops.unk = sel;
117ops.corners = [c00 c10 c01 c11];
118% The factors of K: strain operators and the weights between them, so that
119% K = D (Dxx' Mx Dxx + Dyy' My Dyy + nu (Dxc' Wc Dyc + Dyc' Wc Dxc) + 2 (1 - nu) A Dxy' Dxy).
120% The GPU port applies K in this form (a second difference, weights, a second
121% difference transposed) because in single precision the assembled 13-point
122% stencil loses the smooth modes to rounding; see plate_fd_setup.
123ops.Dxx = Dxx; ops.Dyy = Dyy; ops.Dxc = Dxc; ops.Dyc = Dyc; ops.Dxy = Dxy;
124ops.mkx = mu(kx); ops.mky = mu(ky); ops.mkc = mu(kc);
125ops.D = D; ops.nu = nu;
126end