function faces = get_faces(elements, opt)
%GET_FACES Extract all faces of a FE mesh
%   FACES = GET_FACES(ELEMENTS) extracts all the faces of a fe mesh given
%   by the matrix model.ELEMENTS.
%   The structure of the matrix ELEMENTS is
%   [ElemID ElemType MatID PropID NodID1 NODID2 ... NodIDn]
%   where ElemType can be 23, 24, 34, 36 or 38 for standard elements and
%   122, 133 or 134 for infinite elements.
%   The structure of the matrix FACES is
%   [ElemID FaceType NodID1 ... NodIDn]
%   where ElemID refers to the parent element of the face.
%   FaceType can be 12, 23 or 24

%   Copyright 2008-2012 P. Fiala
%   Budapest University of Technology and Economics
%   Dept. of Telecommunications

if nargin < 2
    % Preserve normal
    opt = 'normal';
end

% zero padding
elements = [elements zeros(size(elements,1), 12-size(elements,2))];

tria = elements(elements(:,2) == 23,[1 (5:7)]);
quad = elements((elements(:,2) == 24) | (elements(:,2) == 122),[1 (5:8)]);
tetra = elements(elements(:,2) == 34,[1 (5:8)]);
pyramid = elements(elements(:,2) == 35,[1 (5:9)]);
penta = elements((elements(:,2) == 36) | (elements(:,2) == 133),[1 (5:10)]);
hexa = elements((elements(:,2) == 38) | (elements(:,2) == 134),[1 (5:12)]);

triafaces = [
    0 1 2 
    0 2 3 
    0 3 1 
    ].';
quadfaces = [
    0 1 2 
    0 2 3 
    0 3 4 
    0 4 1 
    ].';
tetrafaces = [
    0 1 3 2
    0 2 3 4
    0 1 4 3
    0 1 2 4
    ].';
switch opt
    case {'n', 'norm', 'normal'}
        pyramidtriafaces = [
            0 1 5 2
            0 2 5 3
            0 3 5 4
            0 4 5 1
        ].';
    
    
        pyramidquadfaces = [
            0 1 2 3 4
        ].';
        
        pentatriafaces = [
            0 1 2 3
            0 4 6 5
            ].';

        pentaquadfaces = [
         %   0 1 2 5 4
         %   0 2 3 6 5
         %   0 3 1 4 6
            0 1 4 5 2
            0 2 5 6 3
            0 3 6 4 1
            ].';
        
        % NOTE: This face ordering is correct for OpenFOAM cyclicity
        hexafaces = [
        %    0 1 4 3 2
        %    0 5 6 7 8
        %    0 1 2 6 5
        %    0 2 3 7 6
        %    0 3 4 8 7
        %    0 4 1 5 8
            0 1 2 3 4
            0 5 8 7 6
            0 1 5 6 2
            0 4 3 7 8
            0 2 6 7 3
            0 1 4 8 5
            ].';
    case {'o', 'ord', 'order'}
        pentatriafaces = [
            0 1 2 3
            0 4 5 6
            ].';

        pentaquadfaces = [
         %   0 1 2 5 4
         %   0 2 3 6 5
         %   0 3 1 4 6
            0 1 4 5 2
            0 2 5 6 3
            0 3 6 4 1
            ].';
        
        hexafaces = [
        %    0 1 4 3 2
        %    0 5 6 7 8
        %    0 1 2 6 5
        %    0 2 3 7 6
        %    0 3 4 8 7
        %    0 4 1 5 8
            0 1 2 3 4
            0 5 6 7 8
            0 1 5 6 2
            0 4 8 7 3
            0 2 6 7 3
            0 1 5 8 4 
            ].';
    otherwise
        error('get_faces:invalid_option', 'Invalid option');
end

triafac = reshape(tria(:,triafaces(:)+1).',2+1,3*size(tria,1)).';
quadfac = reshape(quad(:,quadfaces(:)+1).',2+1,4*size(quad,1)).';
tetrafac = reshape(tetra(:,tetrafaces(:)+1).',3+1,4*size(tetra,1)).';
pyramidtriafac = reshape(pyramid(:,pyramidtriafaces(:)+1).',3+1,4*size(pyramid,1)).';
pyramidquadfac = reshape(pyramid(:,pyramidquadfaces(:)+1).',4+1,1*size(pyramid,1)).';
pentatriafac = reshape(penta(:,pentatriafaces(:)+1).',3+1,2*size(penta,1)).';
pentaquadfac = reshape(penta(:,pentaquadfaces(:)+1).',4+1,3*size(penta,1)).';
hexafac = reshape(hexa(:,hexafaces(:)+1).',4+1,6*size(hexa,1)).';

ntr = size(triafac,1);
nq = size(quadfac,1);
nte = size(tetrafac,1);
npy3 = size(pyramidtriafac,1);
npy4 = size(pyramidquadfac,1);
np3 = size(pentatriafac,1);
np4 = size(pentaquadfac,1);
nh = size(hexafac,1);

faces = zeros(ntr+nq+nte+np3+np4+nh,6);
faces(1:ntr,[1 (3:4)]) = triafac;
faces(ntr+(1:nq),[1 (3:4)]) = quadfac;
faces(ntr+nq+(1:nte),[1 (3:5)]) = tetrafac;
faces(ntr+nq+nte+(1:npy3),[1 (3:5)]) = pyramidtriafac;
faces(ntr+nq+nte+npy3+(1:npy4),[1 (3:6)]) = pyramidquadfac;
faces(ntr+nq+nte+npy3+npy4+(1:np3),[1 (3:5)]) = pentatriafac;
faces(ntr+nq+nte+npy3+npy4+np3+(1:np4),[1 (3:6)]) = pentaquadfac;
faces(ntr+nq+nte+npy3+npy4+np3+np4+(1:nh),[1 (3:6)]) = hexafac;
faces(faces(:,6) ~= 0, 2) = 24;
faces(faces(:,6) == 0, 2) = 23;
faces(faces(:,5) == 0, 2) = 12;

faces = faces(:,any(faces, 1));

end
