1function [L, K, M, idx, fixed] = plate_matrices(dom, bc, nu)
2% PLATE_MATRICES Sparse Laplacian, stiffness and mass matrices for the plate.
3%
4% [L, K, M, idx, fixed] = plate_matrices(dom, bc, nu) assembles the spatial
5% operators of the finite-difference schemes on the footprint dom from
6% plate_domain, as sparse matrices acting on the vector w(idx) of unknowns.
7% idx lists the unknown nodes in column-major order of the ny x nx field
8% (y running fastest, as w(:) does). fixed is the unknown held at w = 0
9% anyway (empty except for the pinned centre); the solvers zero it after
10% every step. The semi-discrete equation is
11%
12% M w_tt + 2 sigma0 M w_t - 2 sigma1 L w_t + kappa^2 K w = F(t) M s / mu
13%
14% Simply supported and clamped edges (w = 0 on the boundary ring). The
15% unknowns are dom.interior, M is the identity, L the five-point Laplacian
16% with zero Dirichlet data and K = L L, the thirteen-point biharmonic. On
17% the rectangle this is exactly the ghost-node scheme of plate_fd and
18% plate_pad: L L encodes the odd reflection (w_nn = 0). For clamped edges
19% the even reflection makes the Laplacian nonzero on the ring, and the term
20%
21% K = L L + sum_q m_q v_q v_q'
22%
23% adds it: v_q is the Laplacian at ring node q with each outside neighbour
24% mirrored onto the inside one, and m_q the trapezoidal weight of q. On the
25% rectangle this is the 7/dx^4 end coefficient in place of 5/dx^4; on a
26% staircase edge it treats each face of the staircase as clamped.
27%
28% Simply supported curved edges and the plate paradox. With w = 0 on the
29% boundary the bending energy is
30%
31% U = D/2 * [ integral (lap w)^2 dA - (1 - nu) * contour kappa w_n^2 ds ]
32%
33% where kappa is the curvature of the edge. L L is the first term alone.
34% On a polygon, and on a staircase, kappa vanishes, so the discrete plate
35% converges to one with lap w = 0 on the edge instead of zero bending
36% moment: the Babuska (Sapondzhyan) paradox, which puts the fundamental of a
37% simply supported disk 16 percent too high. For the ellipse the missing
38% term is added with the curvature and normal of the true edge nearest each
39% ring node (dom.kappa), w_n taken by one-sided differences into the plate
40% and the arc length shared out over the faces of the staircase. It is a
41% negative correction, small against the entries of L L but of order one
42% in the energy of a smooth mode, and it restores the simply supported
43% limit. Nothing is added for the rectangle, whose edges are straight.
44%
45% Free edges. Every node of the plate is an unknown (dom.inside). The
46% free-edge conditions, no bending moment and no Kirchhoff shear on an edge
47% and no twisting moment at a corner, involve third derivatives and do not
48% fit in a ghost layer. K is therefore assembled from the bending energy
49%
50% U = D/2 * integral( w_xx^2 + w_yy^2 + 2 nu w_xx w_yy + 2 (1-nu) w_xy^2 ) dA
51%
52% using second differences at the nodes that have both neighbours, the
53% cross difference at the cell centres and the trapezoidal node weights
54% dom.m (the finite-difference energy method; Bushnell, Computers &
55% Structures 1973). The free-edge conditions are then natural boundary
56% conditions of the discrete energy: the scheme is second order inside and
57% first order at the boundary. M holds the node weights and L is the
58% Laplacian energy operator with the natural zero-normal-derivative
59% condition, each grid edge weighted by its share of cells. K is symmetric
60% positive semidefinite and its null space is exactly the three rigid-body
61% motions 1, x and y, which plate_rigid_modes removes from the state. All
62% three matrices are per unit cell area, so inside the plate M = 1 and K
63% and L reduce to the usual stencils, and the largest eigenvalue of M \ K
64% stays near the interior bound (4/dx^2 + 4/dy^2)^2 that plate_fd_plan
65% uses for the time step.
66%
67% Pinned centre (a cymbal on a stand). The edges are free and the unknown
68% nearest the centre is held at w = 0: its row and column are zeroed in K
69% and L, so its equation decouples (M keeps its unit diagonal and A stays
70% positive definite) and the other nodes never see it. One pinned point
71% leaves the two rigid rotations about it at zero frequency;
72% plate_rigid_modes with the pin as pivot removes them, standing in for the
73% felt washers of a stand. Under a point support the curvature grows like
74% log r, so the stress right at the support depends on the grid, as it does
75% under a point load.
77if nargin < 3
78 nu = NaN;
79end
80nx = dom.nx; ny = dom.ny; dx = dom.dx; dy = dom.dy;
81dirichlet = any(strcmp(bc, {'simply-supported', 'clamped'}));
82energy = any(strcmp(bc, {'free', 'pinned-center'}));
83if ~dirichlet && ~energy
84 error('plate_matrices: unknown boundary condition "%s"', bc);
85end
86if energy && ~isfinite(nu)
87 error('plate_matrices: free edges need Poisson''s ratio nu');
88end
89if dirichlet
90 unk = dom.interior;
91else
92 unk = dom.inside;
93end
94idx = find(unk);
95n = numel(idx);
96me = (1:n)';
98% Unknown number of every node (0 for none), on a padded grid so that every
99% node has four neighbours to look up
100num_in = zeros(ny, nx);
101num_in(idx) = me;
102num = zeros(ny + 2, nx + 2);
103num(2:end - 1, 2:end - 1) = num_in;
104nbL = num(2:end - 1, 1:end - 2); nbR = num(2:end - 1, 3:end); % neighbours at x - dx, x + dx
105nbD = num(1:end - 2, 2:end - 1); nbU = num(3:end, 2:end - 1); % neighbours at y - dy, y + dy
106uL = nbL(idx); uR = nbR(idx); uD = nbD(idx); uU = nbU(idx); % of the unknowns, in unknown order
108if dirichlet
109 % ---- Dirichlet Laplacian: full diagonal, off-diagonals only to unknowns
110 ii = me; jj = me; vv = -(2 / dx^2 + 2 / dy^2) * ones(n, 1);
111 nbrs = {uL, uR, uD, uU};
112 coef = [1 / dx^2, 1 / dx^2, 1 / dy^2, 1 / dy^2];
113 for k = 1:4
114 nb = nbrs{k};
115 has = find(nb > 0);
116 ii = [ii; me(has)]; jj = [jj; nb(has)]; vv = [vv; coef(k) * ones(numel(has), 1)];
117 end
118 L = sparse(ii, jj, vv, n, n);
119 K = L * L;
120 M = speye(n);
121 fixed = zeros(0, 1);
123 % Ring nodes and their unknown neighbours
124 kb = find(dom.boundary);
125 nb = numel(kb);
126 qb = (1:nb)';
127 bL = nbL(kb); bR = nbR(kb); bD = nbD(kb); bU = nbU(kb);
128 if strcmp(bc, 'clamped')
129 % Laplacian at each ring node with the outside neighbour mirrored inside
130 ii = []; jj = []; vv = [];
131 for k = 1:2
132 if k == 1
133 p1 = bL; p2 = bR; c = 1 / dx^2;
134 else
135 p1 = bD; p2 = bU; c = 1 / dy^2;
136 end
137 both = find(p1 > 0 & p2 > 0);
138 one1 = find(p1 > 0 & p2 == 0);
139 one2 = find(p2 > 0 & p1 == 0);
140 ii = [ii; qb(both); qb(both); qb(one1); qb(one2)];
141 jj = [jj; p1(both); p2(both); p1(one1); p2(one2)];
142 vv = [vv; c * ones(2 * numel(both), 1); 2 * c * ones(numel(one1) + numel(one2), 1)];
143 end
144 V = sparse(ii, jj, vv, nb, n);
145 mq = dom.m(kb);
146 K = K + V' * spdiags(mq, 0, nb, nb) * V;
147 elseif any(dom.kappa(:) ~= 0)
148 % Simply supported curved edge: - (1 - nu) contour kappa w_n^2 ds, with
149 % w_n^2 = w_x^2 + w_y^2 on an edge where w = 0, the gradient by one-sided
150 % (or, with unknowns on both sides, central) differences at each ring
151 % node, and ds the faces of the staircase projected onto the true edge
152 if ~isfinite(nu)
153 error('plate_matrices: a simply supported curved edge needs Poisson''s ratio nu');
154 end
155 ds = dy * abs(dom.normal_x(kb)) .* ((bL > 0) + (bR > 0)) + dx * abs(dom.normal_y(kb)) .* ((bD > 0) + (bU > 0));
156 wq = (1 - nu) * dom.kappa(kb) .* ds / (dx * dy);
157 for k = 1:2
158 if k == 1
159 p1 = bL; p2 = bR; h = dx;
160 else
161 p1 = bD; p2 = bU; h = dy;
162 end
163 both = find(p1 > 0 & p2 > 0);
164 one1 = find(p1 > 0 & p2 == 0);
165 one2 = find(p2 > 0 & p1 == 0);
166 ii = [qb(both); qb(both); qb(one1); qb(one2)];
167 jj = [p1(both); p2(both); p1(one1); p2(one2)];
168 vv = [-ones(numel(both), 1) / (2 * h); ones(numel(both), 1) / (2 * h); ones(numel(one1) + numel(one2), 1) / h];
169 G = sparse(ii, jj, vv, nb, n);
170 K = K - G' * spdiags(wq, 0, nb, nb) * G;
171 end
172 end
173else
174 % ---- Energy method on every node of the plate
175 mu = dom.m(idx);
176 hasX = uL > 0 & uR > 0;
177 hasY = uD > 0 & uU > 0;
178 kx = find(hasX); nkx = numel(kx);
179 ky = find(hasY); nky = numel(ky);
180 Dxx = 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);
181 Dyy = 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);
182 Kxx = Dxx' * spdiags(mu(kx), 0, nkx, nkx) * Dxx; % w_xx^2 where both x-neighbours exist
183 Kyy = Dyy' * spdiags(mu(ky), 0, nky, nky) * Dyy; % w_yy^2 where both y-neighbours exist
184 kc = find(hasX & hasY); nkc = numel(kc);
185 Dxc = 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);
186 Dyc = 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);
187 Wc = spdiags(mu(kc), 0, nkc, nkc);
188 Kc = nu * (Dxc' * Wc * Dyc + Dyc' * Wc * Dxc); % 2 nu w_xx w_yy where both exist
189 % Cells with four unknown corners: the twist at the cell centre
190 A00 = num_in(1:end - 1, 1:end - 1); A10 = num_in(2:end, 1:end - 1);
191 A01 = num_in(1:end - 1, 2:end); A11 = num_in(2:end, 2:end);
192 cu = A00 > 0 & A10 > 0 & A01 > 0 & A11 > 0;
193 kcell = find(cu); ncell = numel(kcell);
194 Dxy = sparse([(1:ncell)'; (1:ncell)'; (1:ncell)'; (1:ncell)'], [A00(kcell); A10(kcell); A01(kcell); A11(kcell)], ...
195 [ones(ncell, 1); -ones(ncell, 1); -ones(ncell, 1); ones(ncell, 1)] / (dx * dy), ncell, n);
196 Kxy = 2 * (1 - nu) * (Dxy' * Dxy); % 2 (1-nu) w_xy^2 at the cell centres
197 K = Kxx + Kyy + Kc + Kxy;
198 % Laplacian energy: first differences along the grid edges, each weighted
199 % by its share of adjacent cells (1 inside, 1/2 along the boundary)
200 cpad = zeros(ny + 1, nx + 1);
201 cpad(2:end - 1, 2:end - 1) = cu;
202 WX = (cpad(2:end, 2:end) + cpad(1:end - 1, 2:end)) / 2; % x-edge from (i,j) to (i,j+1)
203 WY = (cpad(2:end, 2:end) + cpad(2:end, 1:end - 1)) / 2; % y-edge from (i,j) to (i+1,j)
204 ex = find(uR > 0); nex = numel(ex);
205 ey = find(uU > 0); ney = numel(ey);
206 Dx1 = sparse([(1:nex)'; (1:nex)'], [me(ex); uR(ex)], [-ones(nex, 1); ones(nex, 1)] / dx, nex, n);
207 Dy1 = sparse([(1:ney)'; (1:ney)'], [me(ey); uU(ey)], [-ones(ney, 1); ones(ney, 1)] / dy, ney, n);
208 wx = WX(idx(ex)); wy = WY(idx(ey));
209 L = -(Dx1' * spdiags(wx, 0, nex, nex) * Dx1 + Dy1' * spdiags(wy, 0, ney, ney) * Dy1);
210 M = spdiags(mu, 0, n, n);
211 fixed = zeros(0, 1);
212 if strcmp(bc, 'pinned-center')
213 xc = dom.x(end) / 2; yc = dom.y(end) / 2;
214 [~, fixed] = min((dom.X(idx) - xc).^2 + (dom.Y(idx) - yc).^2); % the unknown nearest the centre
215 sel = ones(n, 1); sel(fixed) = 0;
216 Pz = spdiags(sel, 0, n, n);
217 K = Pz * K * Pz; % its row and column drop out
218 L = Pz * L * Pz;
219 end
220end
221end