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