1function dom = plate_domain(plate, nx, ny)
2% PLATE_DOMAIN The footprint of the plate on a uniform nx x ny grid.
3%
4% dom = plate_domain(plate, nx, ny) lays the grid over the bounding box
5% [0, Lx] x [0, Ly] of plate and marks which nodes belong to the plate.
6% For plate.shape = 'rectangle' every node does. For 'ellipse' the plate is
7% the ellipse inscribed in the box, (x - Lx/2)^2 / (Lx/2)^2 +
8% (y - Ly/2)^2 / (Ly/2)^2 <= 1, a disk when Lx = Ly, and its nodes are those
9% within half a grid spacing of the inside, so that the staircase boundary
10% they form straddles the true edge rather than lying to one side of it.
11% Nodes that do not belong to a grid cell with all four corners inside are
12% dropped, so every node carries some area. Both finite-difference methods
13% assemble their operators on this footprint (plate_matrices); the modal
14% method has closed-form modes only for the rectangle.
15%
16% Fields.
17% shape, nx, ny, dx, dy, x, y, X, Y the grid; w(i, j) lives at (x(j), y(i))
18% inside ny x nx logical, the plate's nodes
19% cell (ny-1) x (nx-1) logical, cells with all four corners inside
20% m ny x nx trapezoidal node weights: adjacent plate cells / 4
21% (1 inside, 1/2 on a straight edge, 1/4 at a convex corner)
22% boundary inside nodes with a 4-neighbour outside: the ring where
23% w = 0 is imposed for supported and clamped edges
24% interior inside and not boundary
25% full3 inside nodes whose eight neighbours are inside, where the
26% curvatures can be taken by central differences
27% near ny x nx index of the nearest full3 node, for every inside
28% node (0 outside); plate_curvature_domain copies from it
29% nunknowns_dirichlet, nunknowns_free sizes of the two unknown sets
30% kappa, normal_x, normal_y ellipse only: curvature and outward unit
31% normal of the true edge nearest each boundary node (zero
32% elsewhere); the simply supported ellipse needs them, see
33% plate_matrices
34%
35% Accuracy. A staircase approximates a curved edge to first order in the
36% grid spacing: the geometry, the edge conditions and the stress near the
37% edge all carry O(dx) errors there, against O(dx^2) inside. Refining the
38% grid is the check.
40if ~isfield(plate, 'shape')
41 shape = 'rectangle';
42else
43 shape = plate.shape;
44end
45if ~any(strcmp(shape, {'rectangle', 'ellipse'}))
46 error('plate_domain: unknown plate shape "%s" (expected rectangle or ellipse)', shape);
47end
48Lx = plate.Lx; Ly = plate.Ly;
49dx = Lx / (nx - 1);
50dy = Ly / (ny - 1);
51x = linspace(0, Lx, nx);
52y = linspace(0, Ly, ny);
53[X, Y] = meshgrid(x, y);
54N = nx * ny;
56kappa = zeros(ny, nx);
57normal_x = zeros(ny, nx);
58normal_y = zeros(ny, nx);
59if strcmp(shape, 'rectangle')
60 in0 = true(ny, nx);
61else
62 xc = Lx / 2; yc = Ly / 2; a = Lx / 2; b = Ly / 2;
63 phi = ((X - xc) / a).^2 + ((Y - yc) / b).^2; % 1 on the edge
64 gx = 2 * (X - xc) / a^2; gy = 2 * (Y - yc) / b^2; % gradient of phi
65 g = sqrt(gx.^2 + gy.^2);
66 g(g == 0) = eps; % the centre
67 d = (phi - 1) ./ g; % signed distance to the edge, first order
68 hn = 1 ./ sqrt((gx ./ g / dx).^2 + (gy ./ g / dy).^2); % grid spacing along the normal
69 in0 = d <= 0.5 * hn; % round the edge to the nearest node
70end
71cell = in0(1:end - 1, 1:end - 1) & in0(2:end, 1:end - 1) & in0(1:end - 1, 2:end) & in0(2:end, 2:end);
72cnt = zeros(ny, nx);
73cnt(1:end - 1, 1:end - 1) = cnt(1:end - 1, 1:end - 1) + cell;
74cnt(2:end, 1:end - 1) = cnt(2:end, 1:end - 1) + cell;
75cnt(1:end - 1, 2:end) = cnt(1:end - 1, 2:end) + cell;
76cnt(2:end, 2:end) = cnt(2:end, 2:end) + cell;
77inside = cnt > 0;
78m = cnt / 4;
80% Neighbourhoods, with everything off the grid counting as outside
81pad = false(ny + 2, nx + 2);
82pad(2:end - 1, 2:end - 1) = inside;
83nb4 = pad(1:end - 2, 2:end - 1) & pad(3:end, 2:end - 1) & pad(2:end - 1, 1:end - 2) & pad(2:end - 1, 3:end);
84nb8 = nb4 & pad(1:end - 2, 1:end - 2) & pad(1:end - 2, 3:end) & pad(3:end, 1:end - 2) & pad(3:end, 3:end);
85interior = inside & nb4;
86boundary = inside & ~nb4;
87full3 = inside & nb8;
88if ~any(full3(:))
89 error('plate_domain: the %d x %d grid is too coarse for this %s; no node has a full stencil inside the plate', nx, ny, shape);
90end
92% Nearest node with a full stencil, by flooding outwards from those nodes
93near = zeros(ny, nx);
94kf = find(full3);
95near(kf) = kf;
96todo = inside & ~full3;
97for pass = 1:max(nx, ny)
98 if ~any(todo(:))
99 break
100 end
101 np = zeros(ny + 2, nx + 2);
102 np(2:end - 1, 2:end - 1) = near;
103 shifts = [-1 0; 1 0; 0 -1; 0 1];
104 for k = 1:4
105 sh = np(2 + shifts(k, 1):end - 1 + shifts(k, 1), 2 + shifts(k, 2):end - 1 + shifts(k, 2));
106 take = todo & sh > 0;
107 near(take) = sh(take);
108 todo(take) = false;
109 end
110end
112% Ellipse: curvature and normal of the true edge at the point nearest each
113% boundary node (the radial projection of the node onto the ellipse)
114if strcmp(shape, 'ellipse')
115 kb = find(boundary);
116 t = 1 ./ sqrt(max(phi(kb), eps));
117 xe = xc + (X(kb) - xc) .* t;
118 ye = yc + (Y(kb) - yc) .* t;
119 px = 2 * (xe - xc) / a^2; py = 2 * (ye - yc) / b^2; % gradient of phi on the edge
120 pg = sqrt(px.^2 + py.^2);
121 normal_x(kb) = px ./ pg;
122 normal_y(kb) = py ./ pg;
123 kappa(kb) = ((2 / a^2) * py.^2 + (2 / b^2) * px.^2) ./ pg.^3; % curvature of the level set
124end
126dom.shape = shape;
127dom.nx = nx; dom.ny = ny; dom.dx = dx; dom.dy = dy;
128dom.x = x; dom.y = y; dom.X = X; dom.Y = Y;
129dom.inside = inside; dom.cell = cell; dom.m = m;
130dom.boundary = boundary; dom.interior = interior; dom.full3 = full3; dom.near = near;
131dom.nunknowns_dirichlet = nnz(interior);
132dom.nunknowns_free = nnz(inside);
133dom.kappa = kappa; dom.normal_x = normal_x; dom.normal_y = normal_y;
134dom.N = N;
135end