quicwind/mat_wind_flux_div.m

   1 function DM=mat_wind_flux_div(X,type)
   2 % get matrix of flux divergence by probing its action
   3 % input:
   4 %     X       mesh
   5 %     type    'M': return just the matrix M, 'D': return just the matrix D,
   6 %             'DM': return D*M, default is 'DM'
   7 % output:
   8 %     DM  matrix
   9 if ~exist('type','var')
  10     type = 'DM';
  11 end
  12
  13 check_mesh(X);
  14 wind_template=grad3z(rand(size(X{1})-1),[1 1 1]);  % cell matrix with entries size of u,v,w
  15 n = sum(cellfun(@numel,wind_template));   % size of vector this will act on
  16
  17 Mfun=@(u)cell2vec(wind2flux(vec2cell(u,wind_template),X));
  18 Mmat=fun2mat(Mfun,[n,1,1]);
  19 Dfun=@(u)div3(vec2cell(u,wind_template));
  20 Dmat=fun2mat(Dfun,[n,1,1]);
  21
  22 if (type == 'M')
  23     DM = Mmat;
  24 elseif (type == 'D')
  25     DM = Dmat;
  26 else
  27     DM = Dmat * Mmat;
  28 end
  29
  30 end