1% Sparse operator for staggering. 2% 3% S = sparse_stagger(sz1,m,cyclic) 4% 5% Create a sparse operator for staggering a field in dimension m. 6% The field is a "collapsed" matrix of the size sz1. 7% 8% Input: 9% sz1: size of rhs 10% m: dimension to stagger 11% cyclic: true if domain is cyclic along dimension m. False is the 12% default value 13 14function S = sparse_stagger(sz1,m,cyclic) 15 16 17if nargin == 2 18 cyclic = false; 19end 20 21n1 = prod(sz1); 22 23sz2 = sz1; 24 25if ~cyclic 26 sz2(m) = sz2(m)-1; 27end 28n2 = prod(sz2); 29 30n = length(sz1); 31 32for i=1:n 33 vi{i} = 1:sz2(i); 34end 35 36IJ = cell(n,1); 37 38[IJ{:}] = ndgrid(vi{:}); 39 40for i=1:n 41 IJ{i} = IJ{i}(:); 42end 43 44L1 = [1:n2]'; 45 46L2 = sub2ind(sz1,IJ{:}); 47v = ones(size(L1))/2; 48 49IJ{m} = IJ{m} + 1; 50 51if cyclic 52 IJ{m} = mod(IJ{m}-1,sz1(m))+1; 53end 54 55L2o = sub2ind(sz1,IJ{:}); 56 57S = sparse( ... 58 [L1; L1; ], ... 59 [L2; L2o; ], ... 60 [v; v ], n2 , n1 ); 61 62 63 64% Copyright (C) 2009 Alexander Barth <a.barth@ulg.ac.be> 65% 66% This program is free software; you can redistribute it and/or modify 67% it under the terms of the GNU General Public License as published by 68% the Free Software Foundation; either version 2 of the License, or 69% (at your option) any later version. 70% 71% This program is distributed in the hope that it will be useful, 72% but WITHOUT ANY WARRANTY; without even the implied warranty of 73% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 74% GNU General Public License for more details. 75% 76% You should have received a copy of the GNU General Public License 77% along with this program; If not, see <http://www.gnu.org/licenses/>. 78 79