1% Generate the gradient and Laplacian operators.
2%
3% s = divand_operators(mask,pmn,nu,iscyclic)
4%
5% Form sparse matrixes representing the gradient and Laplacian using
6% finite differences
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%
14% Output:
15%   s: stucture containing
16%   s.Dx: cell array of the gradient
17%   s.D: Laplaciant
18%   s.sv: structure describing the state vector
19%   s.mask: land-sea mask
20%   s.WE: diagonal matrix where each element is the surface
21%     of a grid cell
22
23
24function s = divand_operators(mask,pmn,nu,iscyclic,mapindex)
25
26pmn = cat_cell_array(pmn);
27
28% number of dimensions
29n = size(pmn,ndims(pmn));
30
31sz = size(mask);
32
33sv = statevector_init(mask);
34
35if ~isempty(mapindex)
36  % mapindex is unpacked and referers to unpacked indices
37
38  % range of packed indices
39  % land point map to 1, but those points are remove by statevector_pack
40  i2 = statevector_unpack(sv,(1:sv.n)',1);
41  mapindex_packed = statevector_pack(sv,i2(mapindex));
42
43  % applybc*x applies the boundary conditions to x
44  i = [1:sv.n]';
45  applybc = sparse(i,mapindex_packed(i),ones(1,sv.n),sv.n,sv.n);
46
47  % a halo point is points which maps to a (different) interior point
48  % a interior point maps to itself
49  s.isinterior = (i == mapindex_packed(i));
50  s.isinterior_unpacked = statevector_unpack(sv,s.isinterior);
51
52  s.mapindex_packed = mapindex_packed;
53end
54
55D = divand_laplacian(mask,pmn,nu,iscyclic);
56
57% XXX remove this WE
58d = statevector_pack(sv,1./(prod(pmn,n+1)));
59WE = sparse_diag(sqrt(d));
60
61for i=1:n
62  S = sparse_stagger(sz,i,iscyclic(i));
63
64  % mask on staggered grid
65  ma = (S * mask(:) == 1);
66  s.mask_stag{i} = ma;
67
68  d = sparse_pack(ma) * prod(S * reshape(pmn,numel(mask),n),2);
69  d = 1./d;
70  s.WEs{i} = sparse_diag(sqrt(d));
71end
72
73s.Dx = cell(n,1);
74[s.Dx{:}] = sparse_gradient(mask,pmn,iscyclic);
75
76
77if ~isempty(mapindex)
78  D = applybc * D * applybc;
79  WE = sparse_diag(s.isinterior) * WE;
80
81  for i=1:n
82    S = sparse_stagger(sz,i,iscyclic(i));
83    s.isinterior_stag{i} =  sparse_pack(s.mask_stag{i}) * S * s.isinterior_unpacked(:);
84
85    % the results of 's.Dx{i} * field' satisfies the bc is field does
86    % there is need to reapply the bc on the result
87    s.Dx{i} = s.Dx{i} * applybc;
88  end
89
90  s.applybc = applybc;
91end
92
93s.D = D;
94s.sv = sv;
95s.mask = mask;
96s.WE = WE;
97s.n = n;
98
99% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
100%
101% This program is free software; you can redistribute it and/or modify it under
102% the terms of the GNU General Public License as published by the Free Software
103% Foundation; either version 2 of the License, or (at your option) any later
104% version.
105%
106% This program is distributed in the hope that it will be useful, but WITHOUT
107% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
108% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
109% details.
110%
111% You should have received a copy of the GNU General Public License along with
112% this program; if not, see <http://www.gnu.org/licenses/>.
113