diffuse-react/diffreactstencil.m at master · verbof/diffuse-react · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
function S = diffreactstencil(U, x_old, nx, ny, bnd, alpha, dxs)
% DIFFREACTSTENCIL computes stencil as inner + bounds + corners
%
% Fabio VERBOSIO, Universita` della Svizzera italiana, November 2018
%

bndN = bnd.bndN;
bndS = bnd.bndS;
bndE = bnd.bndE;
bndW = bnd.bndW;

S = zeros(size(U));

%% Interior grid points
for j = 2:ny-1          %(int j=1; j < jend; j++) {
    for i = 2:nx-1      %(int i=1; i < iend; i++) {
        S(i,j) = -(4.0 + alpha) * U(i,j) ...               % central point
                                + U(i-1,j) + U(i+1,j) ... % east and west
                                + U(i,j-1) + U(i,j+1) ... % north and south
                                + alpha * x_old(i,j)  ...
                                + dxs * U(i,j) * (1.0 - U(i,j));
    end
end
% end Interior grid points

%% East boundary
i = nx;                 %int i = nx - 1;
for j = 2:ny-1          %(int j = 1; j < jend; j++)
    S(i,j) = -(4.0 + alpha) * U(i,j)     ...
                + U(i-1,j) + U(i,j-1) + U(i,j+1) ...
                + alpha*x_old(i,j) + bndE(j)     ...
                + dxs * U(i,j) * (1.0 - U(i,j));
end
% end East boundary

%% West boundary
i = 1;                  %int i = 0;
for j = 2:ny-1 %(int j = 1; j < jend; j++)
    S(i,j) = -(4.0 + alpha) * U(i,j)         ...
                + U(i+1,j) + U(i,j-1) + U(i,j+1) ...
                + alpha * x_old(i,j) + bndW(j)   ...
                + dxs * U(i,j) * (1.0 - U(i,j));
end
% end West boundary

%% North boundary, plus NE and NW corners
j = ny;           %int j = ny - 1;

% NW corner
i = 1;            %int i = 0; // NW corner
S(i,j) = -(4.0 + alpha) * U(i,j)         ...
            + U(i+1,j) + U(i,j-1)       ...
            + alpha * x_old(i,j) + bndW(j) + bndN(i) ...
            + dxs * U(i,j) * (1.0 - U(i,j));
% end NW corner

% North boundary w/o corners
for i = 2:nx-1        %(int i = 1; i < iend; i++)
    S(i,j) = -(4.0 + alpha) * U(i,j)        ...
                + U(i-1,j) + U(i+1,j) + U(i,j-1) ...
                + alpha*x_old(i,j) + bndN(i)    ...
                + dxs * U(i,j) * (1.0 - U(i,j));
end
% end North boundary w/o corners

% NE corner
i = nx;        % int i = nx-1; // NE corner
S(i,j) = -(4.0 + alpha) * U(i,j)        ...
            + U(i-1,j) + U(i,j-1)       ...
            + alpha * x_old(i,j) + bndE(j) + bndN(i)    ...
            + dxs * U(i,j) * (1.0 - U(i,j));
% end NE corner

%% South boundary, plus SW and SE corners
j = 1;        % int j = 0;

% SW corner
i = 1;        % int i = 0; // SW corner
S(i,j) = -(4.0 + alpha) * U(i,j)    ...
            + U(i+1,j) + U(i,j+1)   ...
            + alpha * x_old(i,j) + bndW(j) + bndS(i) ...
            + dxs * U(i,j) * (1.0 - U(i,j));
% end SW corner

% South boundary w/o corners
for i = 2:nx-1      %(int i = 1; i < iend; i++)
    S(i,j) = -(4.0 + alpha) * U(i,j)                ...
            + U(i-1,j) + U(i+1,j) + U(i,j+1)    ...
            + alpha * x_old(i,j) + bndS(i)      ...
            + dxs * U(i,j) * (1.0 - U(i,j));
end
% end South boundary w/o corners

% SE corner
i = nx;     % int i = nx - 1; // SE corner
S(i,j) = -(4.0 + alpha) * U(i,j)        ...
            + U(i-1,j) + U(i,j+1)       ...
            + alpha * x_old(i,j) + bndE(j) + bndS(i)    ...
            + dxs * U(i,j) * (1.0 - U(i,j));
% end SW corner

end