function [edges, elem_ids] = get_edges(elements)
%GET_EDGES Extract the edges of elements
%   [EDGES, ELEM_IDS] = GET_EDGES(ELEMENTS) Extract all edges from the 
%       element definition matrix ELEMENTS of a finite element mesh.
%       The result matrix EDGES is an N x 2 matrix containing the IDs of
%       the nodes of all edges. ELEM_IDS is an N x 1 column vector
%       containing the ID of the parent element of each edge.
%       Note: The elements are traversed sorted by element type.
%   

% 2020-2021, Mihaly Ulveczki & Peter Rucz  

elem_types = [12, 23, 24, 34, 36, 38];
num_edges = [1, 3, 4, 6, 9, 12];
[~, type_idx] = ismember(elements(:,2), elem_types);
n_edges = sum(num_edges(type_idx(type_idx ~= 0)));

% Preallocate
edges = zeros(n_edges, 2);
elem_ids = zeros(n_edges, 1);

% Counter
c = 0;

% Main cycle through element types
for i = 1 : length(elem_types)
    % Elements of the current types
    act_elem_inds = find(elements(:,2) == elem_types(i));
    % Skip empty
    if isempty(act_elem_inds)
        continue
    end
    
    switch elem_types(i)
        case 12     % 1D Line element
            act_edges = [1 2];
        case 23     % 2D Tria element
            act_edges = [1 2; ...
                2 3; ...
                3 1];
        case 24     % 2D Quad element
            act_edges = [1 2; ...
                2 3; ...
                3 4; ...
                4 1];
        case 34     % 3D Tetra element
            act_edges = [1 2; ...
                1 3; ...
                1 4; ...
                2 3; ...
                2 4; ...
                3 4];
        case 36     % 3D Wedge element
            act_edges = [1 2; ...
                2 3; ...
                3 1; ...
                4 5; ...
                5 6; ...
                6 4; ...
                1 4; ...
                2 5; ...
                3 6];
        case 38     % 3D Hexa element
            act_edges = [1 2; ...
                2 3; ...
                3 4; ...
                4 1; ...
                5 6; ...
                6 7; ...
                7 8; ...
                8 5; ...
                1 5; ...
                2 6; ...
                3 7; ...
                4 8];
    end
    
    edges_part = elements(act_elem_inds, act_edges.'+4);
    edges_part = reshape(edges_part.',2,[]).';
    
    n_part = size(edges_part, 1);
    
    edges(c + (1 : n_part), :) = edges_part;
    elem_ids(c + (1 : n_part)) = ...
        reshape(repmat(elements(act_elem_inds, 1), 1, num_edges(i)).', [], 1);
    
    c = c + n_part;
end % of cycle through element types

end % of function get_edges
