VelocitySpace.m

classdef VelocitySpace < ILocalizedOperator & ISeparableOperator

  properties (Constant)
    id = 'L_velocity';

    ret_vars = {'velocity'};
    arg_vars = {'saturation', 'pressure'};
    arg_vars_short = {'S', 'P'};
  end

  methods (Static)
    function [cell_eind, cell_glob2loc, lind_cands, nind_cands] = compute_TM_global_pre(model_data, TM)
      lind_cands = zeros(1, model_data.gn_edges);

      % valid boundary indices
      vis = (model_data.gEI ~= 0);
      % valid row indices
      rowinds = repmat((1:size(model_data.gEI,1))', 1, 4);
      vrowinds = rowinds(vis);
      % valid col indices
      colinds = repmat((1:4), size(model_data.gEI,1), 1);
      vcolinds = colinds(vis);

      lind_cands(model_data.gEI(vis)) = vrowinds;
      rind_cands = zeros(1, model_data.gn_edges);
      rind_cands(model_data.gEI(vis)) = model_data.grid.NBI(vis);

      cell_eind = unique([lind_cands(TM), rind_cands(TM(TM <= model_data.gn_inner_edges))]);

      nind_cands = zeros(1, model_data.gn_edges);
      nind_cands(model_data.gEI(vis)) = vcolinds;

      cell_glob2loc = zeros(max(cell_eind),1);
      cell_glob2loc(cell_eind) = 1:length(cell_eind);

    end

    function [nmd] = compute_TM_global_grid(model_data, TM, cell_eind, cell_glob2loc, lind_cands, nind_cands)

      grid = model_data.grid;

      llind_cands = cell_glob2loc(lind_cands(TM));
      lnind_cands = nind_cands(TM)';
      if isa(grid,'rectgrid') || isa(grid,'triagrid') || isa(grid,'onedgrid')
        nmd.grid_local_ext = gridpart(grid,cell_eind);

        gle = nmd.grid_local_ext;
        mask = zeros(size(gle.NBI));
        lr_inds = sub2ind(size(mask), llind_cands, lnind_cands);
        mask(lr_inds) = 1;
        lr_inds_pos = find(gle.NBI(lr_inds) > 0);
        lr_inds_pos = lr_inds(lr_inds_pos);
        mask(sub2ind(size(mask), gle.NBI(lr_inds_pos), gle.INB(lr_inds_pos))) = 1;
        gle.NBI(~mask) = -20;

%        [tmp_i, tmp_j] = find(nmd.grid_local_ext.NBI(lrind_cands, :) < 0);
%        nmd.grid_local_ext.NBI(sub2ind(size(nmd.grid_local_ext.NBI), lrind_cands, tmp_j)) = -10;

        [nmd.gEI, nmd.gn_inner_edges, nmd.gn_boundary_edges] ...
          = compute_edge_indices(nmd.grid_local_ext);
        nmd.gn_edges = nmd.gn_inner_edges + nmd.gn_boundary_edges;
      else
        error('unsupported grid');
      end

      nmd.TM_global_args{1} = cell_eind;
      nmd.TM_global_args{2} = cell_eind;

    end

    function eind_local = compute_TM_global_post(nmd, TM, cell_glob2loc, lind_cands, nind_cands)
      subset = sub2ind(size(nmd.grid_local_ext.NBI), cell_glob2loc(lind_cands(TM)), nind_cands(TM)');
      eind_local = nmd.gEI(subset);
    end

  end

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

      S = sim_data.S;
      P = sim_data.P;
      if ~isfield(sim_data, 'n')
        sim_data.n = model_data.gn_edges;
      end

      if nargin <= 4 || isempty(NU_ind)
        NU_ind = (1:sim_data.n);
      end

      if ~isfield(sim_data, 'm')
        sim_data.m = length(S) + length(P);
      end

      num_diff_flux_u = fv_num_diff_flux_pressure_gradient(model, model_data, P, S);

      UKL = -num_diff_flux_u.G;

      INC = UKL(NU_ind);

      INCoff = [];

      if nargout > 2

        u_diff_p_s = fv_frechet_operators_diff_flux_pressure_gradient(model, model_data, P, S, NU_ind);

        n = length(NU_ind);
        UN_ind = zeros(sim_data.n,1);
        UN_ind(NU_ind) = 1:n;
        if ~isfield(sim_data, 'Poff')
          sim_data.Poff = length(S);
        end
        m = sim_data.m;

        sp21 = sparse(UN_ind(u_diff_p_s.si), u_diff_p_s.sj, u_diff_p_s.svals, n, m);
        sp22 = sparse(UN_ind(u_diff_p_s.pi), u_diff_p_s.pj+sim_data.Poff, u_diff_p_s.pvals, n, m);
        J = -(sp21+sp22);
      end

    end

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

      grid = model_data.grid;
      n = model_data.gn_edges;
      if nargin <= 3
        params = [];
      end
      if ~isfield(params, 'm')
        params.m = grid.nelements;
      end
      u_diff_p_s = fv_frechet_operators_diff_flux_pressure_gradient(descr, model_data, ones(params.m, 1), ones(params.m, 1));

      LL = -sparse(u_diff_p_s.pi, u_diff_p_s.pj, u_diff_p_s.pvals, n, params.m);
      bb = [];
    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 = [];
    end


    function ret_size = ret_size(this, model_data, NU_ind)
      if nargin == 2 || isempty(NU_ind)
        NU_ind = 1:model_data.gn_edges;
      end
      ret_size = length(NU_ind);
    end

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

    end

    function [LL_comps, bb_comps] = gram_project_components(this, rmodel, detailed_data)

      dd_prs = get_by_description(detailed_data, 'pressure');
      RB_prs = dd_prs.RB;

      dd_vel = get_by_description(detailed_data, 'velocity');
      RB_vel = dd_vel.RB;

      L_comps  = components(this, rmodel.descr, dd_prs.model_data);
      LL_comps = cellfun(@(X) X * RB_prs, L_comps, 'UniformOutput', false);

      A = inner_product_matrix(this, dd_prs.model_data);

      LL_comps = cellfun(@(X) RB_vel' * A * X, ...
                      LL_comps, 'UniformOutput', false);
      bb_comps = {};
    end

    function ipm = inner_product_matrix(this, model_data)
%      ipm = 1;
      ipm = model_data.diamondWinv;
    end

    function [nmd, eind, eind_local] = compute_TM_global(this, model_data, TM)

      [cell_eind, cell_glob2loc, lind_cands, nind_cands] = ...
        Fv.TwoPhase.VelocitySpace.compute_TM_global_pre...
           (model_data, TM);

      [nmd] = ...
        Fv.TwoPhase.VelocitySpace.compute_TM_global_grid...
           (model_data, TM, cell_eind, cell_glob2loc, lind_cands, nind_cands);

      [eind_local] = ...
        Fv.TwoPhase.VelocitySpace.compute_TM_global_post...
           (nmd, TM, cell_glob2loc, lind_cands, nind_cands);

      eind = TM;
    end
  end
end