compute_edge_indices.m

Go to the documentation of this file.

function [gEI, in_edges, b_edges, i_ints_bigger, b_ints] = compute_edge_indices(grid)
% edge index matrix. This matrix maps each edge specified by the tuple
% '(element_id, local_edge_id)' to a unique and continuous edge index
% enumeration.

gEI  = zeros(grid.nelements, 4);
% I1 = inner intersections with left element id 'K' smaller than right id 'L'
i_ints_smaller  = grid.NBI > 0 ...
  & grid.NBI > repmat((1:grid.nelements)', 1, grid.nneigh);
% number of innter intersections
in_edges = sum(i_ints_smaller(:));

% set half of the edge index matrix for all edges 'KL'
gEI(i_ints_smaller) = 1:in_edges;

% for all inner edges, 'itmp' has the smaller element id 'K' of this edge ...
[itmp, jtmp]                   = find(i_ints_smaller);
% and 'neighbours' stores the larger element id 'L'.
neighbours = grid.NBI(i_ints_smaller);

% jtmp stores the 'local edge id' for all edges 'LK', for this we create the
% logical matrix 'local_edges'.
local_edges                    = grid.NBI(neighbours, :) ...
    == repmat(itmp, 1, grid.nneigh);
[itmp, jtmp]                   = find(local_edges);

% I2 = inner intersections with left element id 'L' bigger than right id 'K'
i_ints_bigger                  = sub2ind(size(grid.NBI), neighbours, jtmp);

% set other half of edge index matrix for all inner edges 'LK'
gEI(i_ints_bigger)  = 1:in_edges;

% boundary intersections
b_ints = grid.NBI < 0 & grid.NBI > -10;
b_edges = sum(b_ints(:));
gn_edges = in_edges + b_edges;
% set edge index matrix entries for all boundary edges
gEI(b_ints)       = in_edges + 1:gn_edges;