1% Create the advection constrain.
2%
3% c = divand_constr_advec(s,velocity)
4%
5% Create the advection constrain using the specified velocity.
6%
7% Input:
8%   s: structure created by divand_background
9%   velocity: cell array of velocity vectors
10%
11% Output:
12%   c: structure to be used by divand_addc with the following fields: R (a
13%     covariance matrix), H (extraction operator) and yo (specified value for
14%     the constrain).
15
16function c = divand_constr_advec(s,velocity)
17
18velocity = cat_cell_array(velocity);
19
20% check for NaNs
21nancount = sum(isnan(velocity(:)));
22if nancount > 0
23    error('%d velocity values are equal to NaN',nancount);
24end
25
26mask = s.mask;
27n  = s.n;
28iscyclic = s.iscyclic;
29
30sz = size(mask);
31
32A = sparse(0);
33
34vel = reshape(velocity,numel(mask),n);
35
36
37
38for i=1:n
39  S = sparse_stagger(sz,i,iscyclic(i));
40  m = (S * mask(:) == 1);
41
42  d = vel(:,i);
43
44  A = A + sparse_diag(d(mask==1)) * sparse_pack(mask) * S' * sparse_pack(m)' * s.Dx{i};
45end
46
47l = size(A,1);
48
49c.H = A;
50c.yo = sparse(l,1);
51c.R = speye(l);
52
53
54% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
55%
56% This program is free software; you can redistribute it and/or modify it under
57% the terms of the GNU General Public License as published by the Free Software
58% Foundation; either version 2 of the License, or (at your option) any later
59% version.
60%
61% This program is distributed in the hope that it will be useful, but WITHOUT
62% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
63% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
64% details.
65%
66% You should have received a copy of the GNU General Public License along with
67% this program; if not, see <http://www.gnu.org/licenses/>.
68