|
| 1 | +% SPDX-License-Identifier: GPL-3.0-only |
| 2 | +% |
| 3 | +% Copyright 2008-2024 San Diego State University Research Foundation (SDSURF). |
| 4 | +% |
| 5 | +% This program is free software: you can redistribute it and/or modify |
| 6 | +% it under the terms of the GNU General Public License as published by |
| 7 | +% the Free Software Foundation, version 3. |
| 8 | +% |
| 9 | +% This program is distributed in the hope that it will be useful, |
| 10 | +% but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 11 | +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 12 | +% LICENSE file or on the web GNU General Public License |
| 13 | +% <https://www.gnu.org/licenses/> for more details. |
| 14 | +% |
| 15 | +% ------------------------------------------------------------------------ |
| 16 | + |
| 17 | +function G = gradPer(k, m, dx) |
| 18 | +% Returns a m by m+2 one-dimensional mimetic gradient operator |
| 19 | +% when the boundary condition is periodic |
| 20 | +% |
| 21 | +% Parameters: |
| 22 | +% k : Order of accuracy |
| 23 | +% m : Number of cells |
| 24 | +% dx : Step size |
| 25 | + |
| 26 | + % Assertions: |
| 27 | + assert(k >= 2, 'k >= 2'); |
| 28 | + assert(mod(k, 2) == 0, 'k % 2 = 0'); |
| 29 | + assert(m >= 2*k, ['m >= ' num2str(2*k) ' for k = ' num2str(k)]); |
| 30 | + |
| 31 | + % constructing circulant matrix |
| 32 | + V = sparse(1, m); % vector of values for circulant matrix |
| 33 | + idx = repmat(-1,m, m); % matrix of indices for circulant matrix |
| 34 | + idx(:,1) = 1:m; |
| 35 | + idx = cumsum(idx, 2); |
| 36 | + idx = rem(idx+m, m) + 1; |
| 37 | + |
| 38 | + switch k |
| 39 | + case 2 |
| 40 | + V(2:3) = [1, -1]; |
| 41 | + |
| 42 | + case 4 |
| 43 | + V(1:4) = [-1/24, 9/8, -9/8, 1/24]; |
| 44 | + |
| 45 | + case 6 |
| 46 | + V(1:5) = [-25/384, 75/64, -75/64, 25/384, -3/640]; V(end) = 3/640; |
| 47 | + |
| 48 | + case 8 |
| 49 | + V(1:6) = [-245/3072, 1225/1024, -1225/1024, 245/3072, -49/5120, 5/7168]; V(end-1) = -5/7168; V(end) = 49/5120; |
| 50 | + end |
| 51 | + |
| 52 | + % G constructed as a circulant matrix |
| 53 | + G = V(idx); |
| 54 | + G = (1/dx).*G; |
| 55 | +end |
0 commit comments