1% Create the laplacian operator.
2%
3% Lap = divand_laplacian(mask,pmn,nu,iscyclic)
4%
5% Form a Laplacian using finite differences
6% assumes that gradient is zero at "coastline"
7%
8% Input:
9%   mask: binary mask delimiting the domain. 1 is inside and 0 outside.
10%         For oceanographic application, this is the land-sea mask.
11%   pmn: scale factor of the grid.
12%   nu: diffusion coefficient of the Laplacian
13%      field of the size mask or cell arrays of fields
14%
15% Output:
16%   Lap: sparce matrix represeting a Laplaciant
17%
18%
19function Lap = divand_laplacian(mask,pmn,nu,iscyclic)
20
21pmn = cat_cell_array(pmn);
22
23% number of dimensions
24n = size(pmn,ndims(pmn));
25
26sz = size(mask);
27
28if nargin == 2
29  nu = ones(sz);
30end
31
32if isequal(sz,size(nu))
33  % assume diffusion is the same along all dimensions
34  nu = repmat(nu,[ones(1,n) n]);
35end
36
37% nu: 1st index represents the dimensions
38nu = permute(nu,[n+1 1:n]);
39
40% extraction operator of sea points
41H = sparse_pack(mask==1);
42
43
44sz = size(mask);
45
46%DD = sparse(0);
47DD = sparse(prod(sz),prod(sz));
48
49Pmn = reshape(pmn,[prod(sz) n]);
50for i=1:n
51  % operator for staggering in dimension i
52  S = sparse_stagger(sz,i,iscyclic(i));
53
54  % d = 1 for interfaces surounded by water, zero otherwise
55  d = S * mask(:) == 1;
56
57  % metric
58  for j = 1:n
59    tmp = S * Pmn(:,j);
60
61    if j==i
62      d = d .* tmp;
63    else
64      d = d ./ tmp;
65    end
66  end
67
68  % tmp: "diffusion coefficient"  along dimension i
69  tmp = nu(i,:);
70
71  %whos tmp d S nu
72  %keyboard
73  d = d .* (S * tmp(:));
74
75  % Flux operators D
76  % zero at "coastline"
77
78  D = sparse_diag(d) * sparse_diff(sz,i,iscyclic(i));
79  szt = sz;
80
81  if ~iscyclic(i)
82    % extx: extension operator
83    szt(i) = szt(i)+1;
84    extx = sparse_trim(szt,i)';
85
86    D = extx * D;
87  else
88    % shift back one grid cell along dimension i
89    D = sparse_shift(sz,i,iscyclic(i))' * D;
90  end
91
92  % add laplacian along dimension i
93  DD = DD + sparse_diff(szt,i,iscyclic(i)) * D;
94end
95
96ivol = prod(pmn,n+1);
97
98% Laplacian on regualar grid DD
99DD = sparse_diag(ivol) * DD;
100
101% Laplacian on grid with on sea points Lap
102Lap = H * DD * H';
103
104
105
106% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
107%
108% This program is free software; you can redistribute it and/or modify it under
109% the terms of the GNU General Public License as published by the Free Software
110% Foundation; either version 2 of the License, or (at your option) any later
111% version.
112%
113% This program is distributed in the hope that it will be useful, but WITHOUT
114% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
115% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
116% details.
117%
118% You should have received a copy of the GNU General Public License along with
119% this program; if not, see <http://www.gnu.org/licenses/>.
120