function fe = completeFEconnections( fe )
%fe = completeFEconnections( fe )
%   Given a standard FE with just vertexes, edgeends and faceedges defined,
%   construct all the other connectivity data.

    % Map each edge to its incident faces.
    fe.edgefaces = zeros( size(fe.edgeends,1), 2 );
    for fi=1:size(fe.faceedges,1)
        for fj = 1:size(fe.faceedges,2)
            ei = fe.faceedges(fi,fj);
            if ei ~= 0
                if fe.edgefaces( ei, 1 ) == 0
                    fe.edgefaces( ei, 1 ) = fi;
                else
                    fe.edgefaces( ei, 2 ) = fi;
                end
            end
        end
    end

    % Map each face to its neighbouring faces.
    fe.facefaces = zeros( size(fe.faceedges,1), 2 );
    for fi=1:size(fe.faceedges,1)
        for fj = 1:size(fe.faceedges,2)
            ei = fe.faceedges(fi,fj);
            if ei ~= 0
                efaces = fe.edgefaces(ei,:);
                if efaces(1)==fi
                    otherface = efaces(2);
                else
                    otherface = efaces(1);
                end
                fe.facefaces(fi,fj) = otherface;
            end
        end
    end
    
    % Map each vertex to its incident edges and faces.
    
    % Map each face to its vertexes.
    fe.facevertexes = zeros( size(fe.faceedges) );
    for fi=1:size(fe.faceedges,1)
        fes = fe.faceedges(fi,:);
        fes = fes(fes>0);
        faceedgeends = fe.edgeends( fes, : );
        if (faceedgeends(1,1)==faceedgeends(end,1)) || (faceedgeends(1,1)==faceedgeends(end,2))
            lastend = 1;
        else
            lastend = 2;
        end
        fe.facevertexes( fi, 1 ) = faceedgeends( 1, lastend );
        nextvertex = faceedgeends( 1, 3-lastend );
        for fj = 2:length(fes)
            if nextvertex==faceedgeends( fj, 1 );
                nextend = 2;
            else
                nextend = 1;
            end
            fe.facevertexes( fi, fj ) = nextvertex;
            nextvertex = faceedgeends( fj, nextend );
        end
    end
    
    numvxs = size(fe.vertexes,1);
    fe.vertexedges = cell( 1, numvxs );
    fe.vertexfaces = cell( 1, numvxs );
    for fi=1:size(fe.facevertexes,1)
        for fj = 1:size(fe.facevertexes,2)
            % Pick an vertex.
            vi = fe.facevertexes(fi,fj);
            % Check we haven't already processed this vertex.
            if (vi ~= 0) && isempty( fe.vertexedges{vi} )
                currentface = fi;
                fei = fj;
                % Get the corresponding edge.
                currentedge = fe.faceedges(currentface,fei);
                ves = currentedge;
                vfs = currentface;
                while true
                    % Find the previous edge.
                    if fei==1
                        fei1 = find( fe.faceedges(currentface,:) ~= 0, 1, 'last' );
                        currentedge = fe.faceedges(currentface,fei1);
                    else
                        currentedge = fe.faceedges(currentface,fei-1);
                    end
                    % Get the face on the other side of the current edge.
                    nextface = fe.edgefaces(currentedge,1);
                    if nextface==currentface
                        nextface = fe.edgefaces(currentedge,2);
                    end
                    % If we're back to the face we started with, stop.
                    if nextface==fi
                        break;
                    end
                    currentface = nextface;
                    % Find the current edge in that face.
                    fei = find( fe.faceedges( currentface, : )==currentedge, 1 );
                    ves(end+1) = currentedge;
                    vfs(end+1) = currentface;
                end
                fe.vertexedges{vi} = ves;
                fe.vertexfaces{vi} = vfs;
            end
        end
    end
end
