EiReducedDataNode.m

classdef EiReducedDataNode < Greedy.User.IReducedDataNode
  % Implementation of a Greedy.User.IReducedDataNode storing reduced data
  % depending on collateral reduced basis space information only (e.g.
  % interpolation DOFs and a local grid).

  properties
    % @copybrief NonlinEvol.EiRbReducedDataNode.BM
    %
    % @copydetails NonlinEvol.EiRbReducedDataNode.BM
    BM;

    % @copybrief NonlinEvol.EiRbReducedDataNode.grid_local_ext
    %
    % @copydetails NonlinEvol.EiRbReducedDataNode.grid_local_ext
    grid_local_ext;

    % @copybrief NonlinEvol.EiRbReducedDataNode.TM_local
    %
    % @copydetails NonlinEvol.EiRbReducedDataNode.TM_local
    TM_local;

    % @copybrief NonlinEvol.EiRbReducedDataNode.TM_global
    %
    % @copydetails NonlinEvol.EiRbReducedDataNode.TM_global
    TM_global;

    % @copybrief NonlinEvol.EiRbReducedDataNode.Mmass
    %
    % @copydetails NonlinEvol.EiRbReducedDataNode.Mmass
    Mmass;

    % @copybrief NonlinEvol.EiRbReducedDataNode.Mstrich
    %
    % @copydetails NonlinEvol.EiRbReducedDataNode.Mstrich
    Mstrich = 0;

    factor = 0;
  end

  properties (SetAccess = private, Dependent)
    % @copybrief NonlinEvol.EiRbReducedDataNode.M
    %
    % @copydetails NonlinEvol.EiRbReducedDataNode.M
    M;

    % @copybrief NonlinEvol.EiRbReducedDataNode.grid
    %
    % @copydetails NonlinEvol.EiRbReducedDataNode.grid
    grid;
  end

  methods
    function rd = EiReducedDataNode(rmodel, detailed_data)
      % function rd = EiReducedDataNode(rmodel, detailed_data)
      % constructor for the generation of the reduced data.
      %
      % Parameters:
      %   rmodel: of type NonlinEvol.ReducedModel
      %   detailed_data: of type IDetailedData
      %
      % Alternative synopsis:
      % @code EiReducedDataNode(rmodel, ei_reduced_data_node) @endcode copying
      % a reduced data node and (optionally) extracting smaller matrices.
      if nargin == 0
        error('NonlinEvol.EiReducedDataNode constructor needs an argument');
      elseif nargin == 2 && isa(rmodel, 'IReducedModel') ...
          && isa(detailed_data, 'NonlinEvol.EiReducedDataNode')
        copy   = detailed_data;
        copy_extract(rd, copy, rmodel);
      elseif nargin == 2
        fill_fields(rd, rmodel, detailed_data);
      else
        error('Did not find constructor for your arguments');
      end
    end

    function conds = get_conds(this)
      if ~isempty(this.Mmass)
        conds.Mmass  = condest(this.Mmass);
      end
      if ~isempty(this.BM)
        conds.BM     = condest(this.BM);
      end
    end
    function grid = get.grid(this)
      grid = this.grid_local_ext;
    end

    function M = get.M(this)
      M = size(this.BM, 2) - this.Mstrich;
    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 NonlinEvol.ReducedModel
      yesno = this.M ~= rmodel.M || this.Mstrich ~= rmodel.Mstrich;
    end

  end
  methods(Access=private)

    function fill_fields(this, rmodel, detailed_data)
      model_data = detailed_data.model_data;
      A          = model_data.W;
      grid       = model_data.grid;

      this.BM = detailed_data.BM;
      if rmodel.verbose > 11
        disp(['stencil mode: ' , rmodel.stencil_mode]);
        disp(['stencil size: ' , num2str(rmodel.local_stencil_size)]);
      end

      eind = compute_TM_global(this, grid, detailed_data.TM, rmodel);

      this.TM_global = eind;

      if isa(grid,'rectgrid') || isa(grid,'triagrid') || isa(grid,'onedgrid')
        this.grid_local_ext = gridpart(grid,eind);
      %  else % e.g. struct, or not required local grid:
    %    disp('please gen nice grid class for 1d grid...')
      %    reduced_ei_data.grid_local_ext{oi} = onedgrid_gridpart(grid,eind);  
      %    keyboard;
      end
      glob2loc       = zeros(grid.nelements,1);
      glob2loc(eind) = 1:length(eind);
      this.TM_local  = glob2loc(detailed_data.TM);

      this.Mmass = detailed_data.QM' * A * detailed_data.QM;
    end

    function copy_extract(this, reduced_data, rmodel)
      MM = rmodel.M + rmodel.Mstrich;
      if MM > reduced_data.M + reduced_data.Mstrich
        error('M and M'' too large for current size of reduced basis!');
      end

      grid                = reduced_data.grid_local_ext;

      % new version with external routine index_ext:
      eind = compute_TM_global(this, grid, reduced_data.TM_local(1:MM), rmodel);
      this.TM_global      = eind;
      this.grid_local_ext = gridpart(grid,eind);

      glob2loc       = zeros(grid.nelements,1);
      glob2loc(eind) = 1:length(eind);
      this.TM_local = glob2loc(reduced_data.TM_local(1:MM));
      this.BM       = reduced_data.BM(1:MM,1:MM);
      this.Mmass    = reduced_data.Mmass(1:MM,1:MM);

      % get Q^t W Q for M to M+M'
      if rmodel.Mstrich > 0
        QWQ         = this.Mmass(rmodel.M+1:end,MM);
        this.factor = norm(QWQ);
      else
        this.factor = NaN;
      end
    end
  end

  methods (Access = public)
    function eind = compute_TM_global(this, grid, TM, rmodel)
      if isequal(rmodel.stencil_mode, 'vertex')
        mask = zeros(1, grid.nelements);
        nbi  = grid.NBI(TM,:);
        i    = find(nbi > 0);
        % get indices of TM's neighbour-neighbours
        nnbi            = grid.NBI(unique(nbi(i)),:);
        ni              = find(nnbi > 0);
        [nnbuniq,tally] = unique(sort(nnbi(ni)),'first');
        tally           = [tally(2:end);length(nnbi(ni))+1] - tally;
        mask(nnbuniq)   = tally;
        mask(nbi(i))    = 5;
        mask(TM)        = 6;
        mask(mask < 2)  = 0; % skip the elements which have no vertex with the
        % given elements in TM in common
        eind            = find(mask);
      else
        eind = index_ext(grid, TM, rmodel.local_stencil_size);
      end

    end
  end
end