DivergenceSpace.m

classdef DivergenceSpace < ILocalizedOperator & ISeparableOperator

  properties (Constant)
    id = 'L_divergence';

    ret_vars = {'pressure'};
    arg_vars = {'velocity'};
    arg_vars_short = {'U'};
  end

  methods
    function [INC, INCoff, J] = apply(this, descr, model_data, sim_data, NU_ind)
    %function INC = apply(model, model_data, sim_data, NU_ind)

      grid = model_data.grid;

      U = sim_data.U;

      if nargin <= 4 || isempty(NU_ind)
        NU_ind = 1:grid.nelements;
      end

      INC = sum(U(model_data.gEI(NU_ind,:)).*(grid.NX(NU_ind,:)+grid.NY(NU_ind,:)), 2) ...
        - model_data.opdata.S_upper + model_data.opdata.S_lower;

      INCoff = [];

      if nargout > 2
        [J] = full_matrix(this, descr, model_data, sim_data);
      end
    end

  end

  methods

    function [LL, bb] = full_matrix(this, descr, model_data, params)

      grid = model_data.grid;
      n = model_data.grid.nelements;
      if nargin <=3
        params = [];
      end
      if ~isfield(params, 'm')
        params.m = model_data.gn_edges;
      end
      if ~isfield(params, 'Uoff')
        params.Uoff = 0;
      end

      sp3i    = repmat(1:n, 1, 4);
      sp3j    = model_data.gEI + params.Uoff;
      sp3vals = grid.NX + grid.NY;
      LL      = sparse(sp3i,sp3j,sp3vals, n, params.m);

      bb = -model_data.opdata.S_upper + model_data.opdata.S_lower;
    end

    function [LL_comps, bb_comps] = components(this, descr, model_data)
      [LL,bb] = full_matrix(this, descr, model_data);
      LL_comps = {LL};
      bb_comps = {bb};
    end

    function [LL_coeffs, bb_coeffs] = coefficients(this, descr)
      LL_coeffs = 1;
      bb_coeffs = 1;
    end

    function ret_size = ret_size(this, model_data, NU_ind)
      if nargin == 2 || isempty(NU_ind)
        ret_size = model_data.grid.nelements;
      else
        ret_size = length(NU_ind);
      end
    end

    function arg_size = arg_size(this, model_data, NU_ind)
      arg_size = model_data.gn_edges;
    end

    function ipm = inner_product_matrix(this, model_data)
      ipm = speye(size(model_data.W));
    end

    function eind = compute_TM_global(this, grid, TM)
      eind = ILocalizedOperator.compute_TM_global_edge(grid, TM);
    end


    function [LL, bb] = gram_projection(this, model, detailed_data)

      dd_vel = get_by_description(detailed_data, 'velocity');
      dd_prs = get_by_description(detailed_data, 'pressure');
      RB_vel = dd_vel.RB;
      RB_prs = dd_prs.RB;
      Nvel = get_by_description(model.N, 'velocity');
      Nprs = get_by_description(model.N, 'pressure');

      model_data = dd_vel.model_data;
      jn = model_data.gn_edges;
      jm = model_data.grid.nelements;
      ii = model_data.gEI;
      jj = repmat((1:jm)', 1, 4);
      spvals = model_data.grid.NX + model_data.grid.NY;
      Jump = sparse(ii(:), jj(:), spvals(:), jn, jm);
      RB_prs_jump = Jump * RB_prs(:,1:Nprs);

      [~,bb] = full_matrix(this, model.descr, model_data);

%      diamondW = model_data.diamondW;
%      LL = RB_prs_jump' * diamondW * RB_vel;
      LL = RB_prs_jump' * RB_vel(:,1:Nvel);

%      bb = RB_prs' * model_data.W * bb;
      bb = RB_prs(:,1:Nprs)' * bb;
    end

    function [LL_comps, bb_comps] = gram_project_components(this, rmodel, detailed_data)
      [LL, bb] = gram_projection(this, rmodel, detailed_data);
      LL_comps = {LL};
      bb_comps = {bb};
    end
  end

end