EiRbReducedDataNode.m

classdef EiRbReducedDataNode < Greedy.User.IReducedDataNode
  % Reduced data implementation for non-linear evolution problems with finite
  % volume discretizations.
  %
  % See @ref DHO11 for details on the implementations of the reduced matrices
  % and vectors.

  properties(Dependent, SetAccess=private)
    % empirical interpolation matrix `{\bf B}`
    %
    % The matrix entries are:
    % ``({\bf B})_{ij} = q_i(x_j)``
    % for `1 \leq i \leq j \leq N` and vanish elsewhere.
    BM;

    % local grid of type .gridbase with added neighbours such that a sparse
    % evaluation of the empirically interpolated operators in the interpolation
    % points is possible on this grid.
    grid_local_ext;

    % synonym for #grid_local_ext of type .gridbase
    grid;

    % indices of grid entities in the #grid_local_ext structure where the
    % interpolation points are situated.
    TM_local;

    % indices of grid entities in the the global grid structure where the
    % interpolation points are situated.
    TM_global;

    gEI;

    gn_edges;

    gn_inner_edges;

    gn_boundary_edges;

    % empirical interpolation mass matrix.
    %
    % The matrix entries are:
    % ``({ \bf M })_{ij} = \int q_i q_j``
    % for `i,j=1,\ldots,M`
    Mmass;

  end

  properties(Dependent, SetAccess=private)
    % number of reduced basis vectors stored in this data node.
    N;

    % number of collateral reduced basis vectors stored in this data node.
    M;

    % number of collateral reduced basis vectors used for error estimation.
    Mstrich;
  end

  properties
    % handle of object of type TwoPhaseFlow.RbReducedDataNode storing all reduced
    % magnitudes based only on the reduced basis.
    rb_rd;

    % handle of object of type TwoPhaseFlow.EiReducedDataNode storing all reduced
    % magnitudes based only on the collateral reduced basis.
    ei_rd;

    opdata;

    % "transition" mass matrix
    %
    % The matrix entries are:
    % ``({ \bf D })_{ij} = \int q_i \varphi_j``
    % for `i=1,\ldots,M, j=1,\ldots,N`
    DE;

    % empirical interpolation matrix for operator with projection to reduced
    % basis space with `M+M'` basis vectors.
    %
    % This matrix is the solution `{\bf C}`of the equation system
    % ``{ \bf B C } = { \bf D }``
    % with @ref #BM "empirical interpolation matrix" `{ \bf B }` and @ref #DE
    % "transition matrix" `{\bf D}`.
    %
    % In case you have a vector of interpolation evaluations `\bf l` at
    % interpolation points, the interpolation can the be calculated as
    % ``{\bf C l}``.
    CE;

    % empirical interpolation matrix for operator with projection to reduced
    % basis space for the empirical interpolation with only `M` basis vectors.
    %
    % See also: #CE
    CE_red;

    % restriction of the reduced basis vectors to interpolation DOFs of the
    % local grid #grid_local_ext.
    RB_local_ext;
  end

  methods
    function rd = EiRbReducedDataNode(rmodel, ...
                                      rb_reduced_data, rb_detailed_data, ...
                                      ei_reduced_data, ei_detailed_data)
      % Constructor for the generation of the reduced matrices and vectors.
      %
      % Parameters:
      %  rb_reduced_data: pure reduced basis reduced data structure of type TwoPhaseFlow.RbReducedDataNode
      %  rb_detailed_data: detailed data structure of type Greedy.DataTree.Detailed.RbLeafNode
      %                    storing the reduced basis.
      %  ei_reduced_data: pure EI reduced data structure of type TwoPhaseFlow.EiReducedDataNode
      %  ei_detailed_data: detailed data structure of type Greedy.DataTree.Detailed.EiLeafNode
      %                    storing the collateral reduced basis.
      if nargin == 0
        error('TwoPhaseFlow.EiRbReducedDataNode constructor needs an argument');
      elseif nargin == 2 && isa(rmodel, 'IReducedModel') ...
          && isa(rb_reduced_data, 'TwoPhaseFlow.EiRbReducedDataNode')
        copy   = rb_reduced_data;
        copy_extract(rd, copy, rmodel);
      elseif nargin == 5
        fill_fields(rd, rmodel, rb_reduced_data, rb_detailed_data, ei_reduced_data, ei_detailed_data);
      else
        error('Did not find constructor for your arguments');
      end
    end

    function conds = get_conds(this)
      conds_rb = create_tree(this.rb_rd, ...
                             DataTree.ScalarCreator(@(X) get_conds(X)),...
                             [], [], []);
      conds_ei = get_conds(this.ei_rd);
      if ~isempty(this.CE)
        conds.CE     = condest(this.CE);
      end
      if ~isempty(this.CE_red)
        conds.CE_red = condest(this.CE_red);
      end
      conds = structcpy(conds_rb, conds_ei);
    end

  end

  methods(Access=private)

    function fill_fields(this, rmodel, rb_red, rb_det, ei_red, ei_det)
      this.rb_rd = rb_red;
      this.ei_rd = ei_red;

      ret_vars = ei_det.ophandle.ret_vars;
      MM = get_by_description(rmodel.M, ei_det.ophandle.id);
      A        = ei_det.ophandle.inner_product_matrix(ei_det.model_data);

      RB = cell(1, length(ret_vars));
      for i = 1:length(ret_vars)
        rb_det_i = get_by_description(rb_det, ret_vars{i});
        N_i = get_by_description(rmodel.N, ret_vars{i});
        RB{i} = rb_det_i.RB(:,1:N_i);
      end
      RB = cat(1,RB{:});
      this.DE = RB' * A * ei_det.QM(:,1:MM);

      arg_vars = ei_det.ophandle.arg_vars;
      arg_vars_short = ei_det.ophandle.arg_vars_short;

      this.RB_local_ext = [];
      for i = 1:length(arg_vars);
        rb_det_i = get_by_description(rb_det, arg_vars{i});
        TM_global_i = this.ei_rd.TM_global_args{i};
        avs = arg_vars_short{i};
        N_i = get_by_description(rmodel.N, arg_vars{i});
        this.RB_local_ext.(avs) = rb_det_i.RB(TM_global_i, 1:N_i);
      end

      this.opdata = ei_det.ophandle.fill_opdata(rmodel, ei_det);
    end

    function copy_extract(this, copy, rmodel)
      MM = get_by_description(rmodel.M, copy.ei_rd.ophandle.id);
      ophandle = copy.ei_rd.ophandle;
      arg_vars = ophandle.arg_vars;
      arg_vars_short = ophandle.arg_vars_short;
      ret_var = ophandle.ret_vars;
      N_ret = get_by_description(rmodel.N, ret_var{1});

      this.ei_rd        = copy.ei_rd;
      this.rb_rd        = copy.rb_rd;
      this.DE           = copy.DE(1:N_ret, 1:MM);
      for i = 1:length(arg_vars_short)
        av = arg_vars{i};
        avs = arg_vars_short{i};
        N_arg = get_by_description(rmodel.N, av);
        if MM == copy.M
          this.RB_local_ext.(avs) = copy.RB_local_ext.(avs)(:, 1:N_arg);
        else
          if ~isempty(this.ei_rd.TM_reduction_args)
            TM_reduction_i = this.ei_rd.TM_reduction_args{i};
            this.RB_local_ext.(avs) = copy.RB_local_ext(TM_reduction_i, 1:N_arg);
          else
            error('the empirical interpolation points are _not_ cropped!');
          end
        end
      end
      this.opdata = ophandle.copy_extract_opdata(copy.opdata, rmodel, this.ei_rd.ophandle.id);
    end
  end

  methods

    function N  = get.N(this)
      N = create_tree(this.rb_rd, ...
                      DataTree.ScalarCreator(@(X) X.N), ...
                      [], [], []);
    end

    function M  = get.M(this)
      M = this.ei_rd.M;
    end

    function Mstrich  = get.Mstrich(this)
      Mstrich = this.ei_rd.Mstrich;
    end

    function BM = get.BM(this)
      BM = this.ei_rd.BM;
    end

    function Mmass = get.Mmass(this)
      Mmass = this.ei_rd.Mmass;
    end

    function grid_local_ext = get.grid_local_ext(this)
      grid_local_ext = this.ei_rd.grid_local_ext;
    end

    function grid_local_ext = get.grid(this)
      grid_local_ext = this.ei_rd.grid_local_ext;
    end

    function TM_global = get.TM_global(this)
      TM_global = this.ei_rd.TM_global;
    end

    function gEI = get.gEI(this)
      gEI = this.ei_rd.gEI;
    end

    function gn_edges = get.gn_edges(this)
      gn_edges = this.ei_rd.gn_edges;
    end

    function gn_inner_edges = get.gn_inner_edges(this)
      gn_inner_edges = this.ei_rd.gn_inner_edges;
    end

    function gn_boundary_edges = get.gn_boundary_edges(this)
      gn_boundary_edges = this.ei_rd.gn_boundary_edges;
    end

    function TM_local = get.TM_local(this)
      TM_local = this.ei_rd.TM_local;
    end

    function yesno = needs_subset_copy(this, rmodel)
      % function yesno = needs_subset_copy(this, rmodel)
      % @copybrief .Greedy.User.IReducedDataNode.needs_subset_copy()
      %
      % @copydetails .Greedy.User.IReducedDataNode.needs_subset_copy()
      %
      % Parameters:
      %   rmodel: of type TwoPhaseFlow.ReducedModel

      yesno = needs_subset_copy(this.rb_rd, rmodel) || needs_subset_copy(this.ei_rd, rmodel);
    end


  end
end