ILocalizedOperator.m

classdef ILocalizedOperator < handle
  % classdef ILocalizedOperator
  % Interface for a <em>localized operator</em> that fulfills the so-called
  % `H`-independent DOF dependence.
  %
  % A discrete operator `{\cal L}_h:{\cal V}_h \to {\cal W}_h` fulfills the
  % <em>`H`-independent DOF dependence</em> if any operator evaluation
  % restricted to a single DOF `\tau[{\cal L}_h[u_h]]` can be computed with a
  % constant complexity `J` independently of the dimension `H` for all DOFs
  % `\tau \in \Sigma({\cal W}_h)` and functions `u_h \in {\cal V}_h`.
  %
  % Example:
  %   - Grid based operators with local stencil for operator evaluations like
  %     most finite volume and finite element operators.
  %
  % Definitions:
  %   - Let `\Sigma({\cal V}_h):=\{\sigma_i\}_{i=1}^H` and
  %     `\Sigma({\cal W}_h):=\{\tau_i\}_{i=1}^H` be the DOFs of the discrete
  %     functions spaces `{\cal V}_h` and `{\cal W}_h`, respectively.
  %   - We denote a map of <em>interpolation DOF indices</em>
  %    `{\cal J}:\{1,\ldots,M\} \to \{1,\ldots,H\}` and
  %   - a map of <em>restriction DOF indices</em>
  %    `{\cal I}:\{1,\ldots,M\} \to \text{Pot}(\{1,\ldots,H\})`, such that
  %    ``\tau_{{\cal J}(m)}\left[
  %        {\cal L}_h \left[
  %           \sum_{m' \in {\cal I}(m)} \sigma_{m'}[u_h] \psi_{m'}
  %                  \right]
  %                       \right] =
  %       \tau_{{\cal J}(m)}[{\cal L}_h[u_h]]``
  %     for all `m=1,\ldots,M` and `u_h \in {\cal V}_h`.
  %     Here `\psi_i \in {\cal V}_h, i=1,\ldots,H` are basis functions of the
  %     discrete argument function space.
  %   - Sometimes it is desirable, to reduce the complexity of the discrete
  %     function spaces from `H` to `M` and `M_{\text{ext}}`. In this case, we
  %     obtain restricted DOF sets `\hat{\Sigma}_{M_{\text{ext}}}({\cal V}_h)` and
  %     `\hat{\Sigma}_{M}({\cal W}_h)`, and indices map `\hat{\cal J}_M` and
  %     `\hat{\cal I}_{M_{\text{ext}}}`.
  %
  properties (Constant, Abstract)
    id;      % a string name specifying the implemented operator
  end

  methods (Abstract)

    % function [INC, INCoff, J] = apply(this, model, model_data, arg, NU_ind);
    % compute operator evaluation on argument function for given
    % <em>interpolation DOF indices</em> specified by 'NU_ind'
    %
    % Parameters:
    %   arg:       the argument function `u_h \in {\cal V}_h`
    %   NU_ind:    interpolation DOF indices
    %              `\nu_{\text{ind}} \subset \text{range}(\hat{\cal J}_M)`.
    %              ( Default: [] meaning `\{1,\ldots,H\}` )
    %
    % Return values:
    %   INC:       vector of DOFs with indices `\nu_{\text{ind}}` of an
    %              operator evaluation
    %   INCoff:    optionally stores the DOFs of a constant part of an
    %              operator. This has been added for compatibility with
    %              separable affine operators.
    %   J:         Jacobian `D[{\cal L}_h]|_{u_h}` restricted to the vector of
    %              interpolation DOFs. This ought to be a matrix of size
    %              `M_{\text{ext}}\times M` with `M_{\text{ext}}:= \text{card}(
    %              \cup_{m\in\nu_{\text{ind}}} {\cal J}(m))`.
    [INC, INCoff, J] = apply(this, model, model_data, arg, NU_ind);


    % function ret_size = ret_size(this, model_data, NU_ind);
    % returns the length of the return 'INC' value of an apply() method call.
    %
    % This should usually return `\text{card}(\nu_{\text{ind}})` or `H`.
    %
    % Parameters:
    %   NU_ind:    interpolation DOF indices
    %              `\nu_{\text{ind}} \subset \text{range}(\hat{\cal J}_M)`.
    %              ( Default: [] meaning `\{1,\ldots,H\}` )
    %
    % Return values:
    %   ret_size: length of the return 'INC' value of an apply() method call.
    ret_size = ret_size(this, model_data, NU_ind);

    % function arg_size = arg_size(this, model_data, NU_ind);
    % returns the length of the required argument DOF vector for an apply()
    % method call.
    %
    % This should usually return `\text{card}( \cup_{m\in\nu_{\text{ind}}}
    % {\cal J}(m))` or `H`.
    %
    % Parameters:
    %   NU_ind:    interpolation DOF indices
    %              `\nu_{\text{ind}} \subset \text{range}(\hat{\cal J}_M)`.
    %              ( Default: [] meaning `\{1,\ldots,H\}` )
    %
    % Return values:
    %   arg_size: length of the required argument DOF vector for an apply()
    %             method call.
    arg_size = arg_size(this, model_data, NU_ind);

    % function [nmd,eind,eind_local] = compute_TM_global(this, model_data, TM_local);
    % computes a subset of <em>restriction DOF indices</em>
    % `\cup_{m \in T_{M_{\text{local}}}}
    %   \hat{{\cal I}}_{M_{\text{ext}}}(m)`
    %
    % The most common use-case is the computation of indices for grid cells in
    % the reduced grid 'model_data.grid_local_ext'.
    % For most operators the neighbors of the interpolation DOFs need to be
    % computed somehow. Default implementations for this exist as
    %  - compute_TM_global_vertex()
    %  - compute_TM_global_edge().
    %
    % Parameters:
    %   TM_local: a set of local DOF indices at which the operator shall be
    %             interpolated `T_{M_{\text{local}}}`.
    %
    % Return values:
    %    nmd:        updated 'model_data' structure with new index sets
    %                `{\cal I}_{M_{\text{ext}}}` and `{\cal J}_M`
    %    eind:       a subset of <em>restriction DOF indices</em>
    %                `\cup_{m \in T_{M_{\text{local}}}} {\cal I}(m)`
    %    eind_local: the new localized set of restriciton DOF indices (relative
    %                to the new 'model_data'
    %                `\cup_{m \in T_{M_{\text{local}}}} {\cal I}_{M_{\text{ext}}}(m)`
    [nmd,eind,eind_local] = compute_TM_global(this, model_data, TM_local);

    % function ipm = inner_product_matrix(this, model_data)
    % return inner product matrix for range function space `{\cal Wh}`
    %
    % Return values:
    %   ipm: the inner product matrix of the range function space
    ipm = inner_product_matrix(this, model_data);
  end

  methods (Static)
    function eind = compute_TM_global_vertex(model_data, TM_local)
      % function eind = compute_TM_global_vertex(model_data, TM_local)
      % Specialization of compute_TM_global() for grid based finite volume operators.
      %
      % If each DOF is connected to a grid cell of an underlying grid, this
      % method selects to a set of \selected DOFs all DOFs connected to grid cells
      % neighboring the vertices of all \em selected grid cells.
      %
      % Parameters:
      %   TM_local: a set of local DOF indices at which the operator shall be
      %             interpolated `T_{M_{\text{local}}}`.
      %
      % Return values:
      %    eind:    a subset of <em>restriction DOF indices</em>
      %             `\cup_{m \in T_{M_{\text{local}}}}
      %               \hat{{\cal I}}_{M_{\text{ext}}}(m)`
      grid = model_data.grid;
      mask = zeros(1, grid.nelements);
      nbi  = grid.NBI(TM_local,:);
      i    = find(nbi > 0);
      % get indices of TM_local'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_local)  = 6;
      mask(mask < 2)  = 0; % skip the elements which have no vertex with the
      % given elements in TM_local in common
      eind            = find(mask);
    end

    function eind = compute_TM_global_edge(model_data, TM_local, local_stencil_size)
      % function eind = compute_TM_global_edge(model_data, TM_local, local_stencil_size)
      % Specialization of compute_TM_global() for grid based finite volume operators.
      %
      % If each DOF is connected to a grid cell of an underlying grid, this
      % method selects to a set of \selected DOFs all DOFs connected to grid cells
      % neighboring the edges of all \em selected grid cells. The inclusion of
      % edge-neighbors is repeated as often as indicated by
      % 'local_stencil_size'
      %
      % Parameters:
      %   TM_local: a set of local DOF indices at which the operator shall be
      %             interpolated `T_{M_{\text{local}}}`.
      %   local_stencil_size: indicates how often the inclusion of
      %             edge-neighors shall be repeated. (Default: 1)
      %
      % Return values:
      %    eind:    a subset of <em>restriction DOF indices</em>
      %             `\cup_{m \in T_{M_{\text{local}}}}
      %               \hat{{\cal I}}_{M_{\text{ext}}}(m)`
      %
      if nargin < 3
        local_stencil_size = 1;
      end
      grid = model_data.grid;
      eind = index_ext(grid, TM_local, local_stencil_size);
    end

    function opdata = fill_opdata(rmodel, detailed_data)
    % function opdata = fill_opdata(rmodel, detailed_data)
    % This method can be called during empirical interpolation in order to
    % collect arbitrary data (for example parameter independent or constant
    % data) in a detailed data node.
    %
    % Parameters:
    %   detailed_data: of type @IDetailedData
    %
    % Return values:
    %   opdata:        struct storing arbitrary data on the operator
      opdata = [];
    end

    function opdata = copy_extract_opdata(opdata_copy, rmodel, Mid)
    % function opdata = fill_opdata(rmodel, detailed_data)
    % Used by a copy constructor of IDetailedData objects which collect
    % operator data with fill_opdata() methods.
    %
    % Parameters:
    %   opdata_copy: opdata structure to be copied
    %   Mid:         Id of the 'detailed_data' structure, i.e.\ the
    %                interpolation DOFs, this operator is connected to. This
    %                can differ from the Id id.
    %
    % Return values:
    %   opdata:        copied struct storing arbitrary data on the operator
      opdata = [];
    end

  end
end