GlobalSpace.m

classdef GlobalSpace < ILocalizedOperator

  properties(Constant)
    id = 'L_global';

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

  properties(Access=private)
    prs_op = Fv.TwoPhase.PressureMean;
    div_op = Fv.TwoPhase.DivergenceSpace;
    sat_op = Fv.TwoPhase.SaturationSpace;
    vel_op = Fv.TwoPhase.VelocitySpace;
  end

  methods
    function [INC, INCoff, J] = apply(this, model, model_data, sim_data, NU_ind)
      Uoff = model_data.grid.nelements;
      Poff = Uoff + model_data.gn_edges;
      NU_indS = NU_ind(NU_ind<=Uoff);
      NU_indU = NU_ind(NU_ind>Uoff & NU_ind<=Poff) - Uoff;
      NU_indP = NU_ind(NU_ind>Poff) - Poff;

      if nargout == 2
        [INCS, INCSoff] = sat_op.apply(this, model, model_data, sim_data, NU_indS);
        [INCU, INCUoff] = vel_op.apply(this, model, model_data, sim_data, NU_indU);
        [INCP, INCPoff] = div_op.apply(this, model, model_data, sim_data, NU_indP);
        [INCM, INCMoff] = prs_op.apply(this, model, model_data, sim_data, []);
      elseif nargout > 2
        [INCS, INCSoff, JS] = sat_op.apply(this, model, model_data, sim_data, NU_indS);
        [INCU, INCUoff, JU] = vel_op.apply(this, model, model_data, sim_data, NU_indU);
        [INCP, INCPoff, JP] = div_op.apply(this, model, model_data, sim_data, NU_indP);
        [INCM, INCMoff, JM] = prs_op.apply(this, model, model_data, sim_data, []);
      end

      INC = [INCS; INCU; INCP; INCM];
      INCoff = [INCSoff; INCUoff; INCPoff; INCMoff];

      if nargout > 2
        J = [JS;
             JU;
             JP;
             JM];
      end
    end


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

    function arg_size = arg_size(this, model_data)
      arg_size = 2 * model_data.grid.nelements + model_data.gn_edges;
    end

    function ipm = inner_product_matrix(this, model_data)
      nT = size(model_data.WT);
      nE = size(model_data.WE);
      ipm = [model_data.WT; sparse(nT(1),nE(2)+nT(2)+1);
             sparse(nE(1),nT(2)), model_data.WE, sparse(nE(1),nT(2)+1);
             sparse(nT(1),nT(2)+nE(2)), model_data.WT, sparse(nT(1),1);
             sparse(1, nT(2)*2+nE(2)), 1];
    end

    function [nmd, eind, eind_local] = compute_TM_global(this, model_data, TM)
      Uoff = model_data.grid.nelements;
      Poff = Uoff + model_data.gn_edges;

      TMS = TM(TM<=Uoff);
      TMU = TM(TM>Uoff & TM<Poff)-Uoff;
      TMP = TM(TM>Poff) - Poff;

      [nmdS, eindS, eindS_local] = ILocalizedOperator.compute_TM_global_edge(model_data, TMS);
      eindU = Fv.TwoPhase.VelocitySpace.compute_TM
      eindP = ILocalizedOperator.compute_TM_global_edge(model_data, TMP)+Poff;
    end

    function opdata = fill_opdata(this, rmodel, detailed_data)
      MM = get_by_description(rmodel.M, 'L_global');
      opdata.S_lower = detailed_data.model_data.opdata.S_lower(detailed_data.TM(1:MM));
      opdata.S_upper = detailed_data.model_data.opdata.S_upper(detailed_data.TM(1:MM));
    end

    function opdata = copy_extract_opdata(this, opdata_copy, rmodel, id)
      MM = get_by_description(rmodel.M, id);
      opdata.S_lower = opdata_copy.S_lower(1:MM);
      opdata.S_upper = opdata_copy.S_upper(1:MM);
    end

  end

end