1% Sparse operator for a gradient. 2% 3% [Dx1,Dx2,...,Dxn] = sparse_gradient(mask,pmn) 4% 5% Form the gradient using finite differences in all n-dimensions 6% 7% Input: 8% mask: binary mask delimiting the domain. 1 is inside and 0 outside. 9% For oceanographic application, this is the land-sea mask. 10% 11% pmn: scale factor of the grid. 12% 13% Output: 14% Dx1,Dx2,...,Dxn: sparce matrix represeting a gradient along 15% different dimensions 16 17function varargout = sparse_gradient(mask,pmn,iscyclic) 18 19H = sparse_pack(mask); 20 21sz = size(mask); 22n = size(pmn,ndims(pmn)); 23Pmn = reshape(pmn,[prod(sz) n]); 24 25if nargin == 2 26 iscyclic = zeros(n,1); 27end 28 29for i=1:n 30 % staggering operator 31 S = sparse_stagger(sz,i,iscyclic(i)); 32 33 % mask for staggered variable 34 m = (S * mask(:) == 1); 35 36 d = m .* (S * Pmn(:,i)); 37 38 varargout{i} = sparse_pack(m) * sparse_diag(d) * sparse_diff(sz,i,iscyclic(i)) * H'; 39% varargout{i} = sparse_diag(d) * sparse_diff(sz,i) * H'; 40% size(varargout{i}) 41end 42 43 44 45% Copyright (C) 2009 Alexander Barth <a.barth@ulg.ac.be> 46% 47% This program is free software; you can redistribute it and/or modify 48% it under the terms of the GNU General Public License as published by 49% the Free Software Foundation; either version 2 of the License, or 50% (at your option) any later version. 51% 52% This program is distributed in the hope that it will be useful, 53% but WITHOUT ANY WARRANTY; without even the implied warranty of 54% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 55% GNU General Public License for more details. 56% 57% You should have received a copy of the GNU General Public License 58% along with this program; If not, see <http://www.gnu.org/licenses/>. 59 60