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
|
% Sparse operator for a gradient.
%
% [Dx1,Dx2,...,Dxn] = sparse_gradient(mask,pmn)
%
% Form the gradient using finite differences in all n-dimensions
%
% Input:
% mask: binary mask delimiting the domain. 1 is inside and 0 outside.
% For oceanographic application, this is the land-sea mask.
%
% pmn: scale factor of the grid.
%
% Output:
% Dx1,Dx2,...,Dxn: sparce matrix represeting a gradient along
% different dimensions
function varargout = sparse_gradient(mask,pmn,iscyclic)
H = sparse_pack(mask);
sz = size(mask);
n = size(pmn,ndims(pmn));
Pmn = reshape(pmn,[prod(sz) n]);
if nargin == 2
iscyclic = zeros(n,1);
end
for i=1:n
% staggering operator
S = sparse_stagger(sz,i,iscyclic(i));
% mask for staggered variable
m = (S * mask(:) == 1);
d = m .* (S * Pmn(:,i));
varargout{i} = sparse_pack(m) * sparse_diag(d) * sparse_diff(sz,i,iscyclic(i)) * H';
% varargout{i} = sparse_diag(d) * sparse_diff(sz,i) * H';
% size(varargout{i})
end
% Copyright (C) 2009 Alexander Barth <a.barth@ulg.ac.be>
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; If not, see <http://www.gnu.org/licenses/>.
|
