1% Two ways to put the edge of a round drum on a square grid.
2%
3% The grid covers the square [-1, 1]^2 with spacing h = 2/N, and the drum is
4% the unit disk. Every node strictly inside the disk is an unknown; every
5% other node is held at zero. What differs is what the five-point Laplacian
6% does at an unknown whose neighbour lies outside.
7%
8% 'staircase' nothing: the outside neighbour is simply zero, as if the rim
9% ran through it. The drum the grid sees is the union of the
10% cells around the inside nodes, a disk with a jagged edge, and
11% its radius is wrong by up to h.
12% 'sw' Shortley and Weller (1938): if the rim crosses the grid line
13% at a fraction theta of the way to the outside neighbour, use
14% the three-point second difference for unequal spacing,
15%
16% u_xx ~ 2/h^2 [ u_W / (1 + theta) - u_P / theta ]
17% (with u = 0 at the crossing),
18%
19% which is consistent to first order at the rim and gives
20% second-order frequencies. The price is the 1/theta: a node
21% very close to the rim has a very large coefficient, and an
22% explicit time step has to respect it.
23%
24% The result is the stencil, times h^2, one number per node and direction:
25% C at the node itself and E, W, Nn, S for its neighbours at (i, j+1), (i, j-1),
26% (i+1, j), (i-1, j), where i counts y and j counts x. `in` marks the
27% unknowns. The driver runs leapfrog with it.
29function [C, E, W, Nn, S, in] = diskstencil(N, method)
30h = 2 / N;
31x = -1 + (0:N) * h;
32[X, Y] = meshgrid(x, x); % X(i, j) = x(j), Y(i, j) = x(i)
33in = X .^ 2 + Y .^ 2 < 1 - 1e-12;
34n1 = N + 1;
35C = zeros(n1, n1);
36E = zeros(n1, n1);
37W = zeros(n1, n1);
38Nn = zeros(n1, n1);
39S = zeros(n1, n1);
40sw = strcmp(method, 'sw');
41for i = 2:n1 - 1
42 for j = 2:n1 - 1
43 if ~in(i, j)
44 continue;
45 end
46 % x direction: distances to the east and west, as fractions of h
47 te = 1;
48 tw = 1;
49 if ~in(i, j + 1) && sw
50 te = (sqrt(1 - Y(i, j) ^ 2) - X(i, j)) / h;
51 end
52 if ~in(i, j - 1) && sw
53 tw = (X(i, j) + sqrt(1 - Y(i, j) ^ 2)) / h;
54 end
55 tn = 1;
56 ts = 1;
57 if ~in(i + 1, j) && sw
58 tn = (sqrt(1 - X(i, j) ^ 2) - Y(i, j)) / h;
59 end
60 if ~in(i - 1, j) && sw
61 ts = (Y(i, j) + sqrt(1 - X(i, j) ^ 2)) / h;
62 end
63 % the unequal-spacing second difference in each direction; with every
64 % theta equal to 1 this is the ordinary five-point stencil
65 E(i, j) = 2 / (te * (te + tw)) * in(i, j + 1);
66 W(i, j) = 2 / (tw * (te + tw)) * in(i, j - 1);
67 Nn(i, j) = 2 / (tn * (tn + ts)) * in(i + 1, j);
68 S(i, j) = 2 / (ts * (tn + ts)) * in(i - 1, j);
69 C(i, j) = -2 / (te * tw) - 2 / (tn * ts);
70 end
71end
72end