def extract_boundary_edges(elements):
    """
    Some meshes do not specify edge information. Extract it from the
    connectivity matrix.

    Required Arguments
    ------------------
    * elements : element connectivity matrix. Assumes Lagrange elements
                 and GMSH-style ordering.

    Output
    ------
    The output is a list of edge tuples. For example:

        print(extract_boundary_edges(mesh.elements))

        => [(1, 5, 7, -1), (2, 10, 12, -1), (3, 15, 17, -1),
            (4, 20, 22, -1), (5, 6, 8, -1), (6, 2, 9, -1),
            (10, 11, 13, -1), (11, 3, 14, -1), (15, 16, 18, -1),
            (16, 4, 19, -1), (20, 21, 23, -1), (21, 1, 24, -1)]

        print(mesh.edges)

        => [(1, 5, 7, 1), (5, 6, 8, 1), (6, 2, 9, 1), (2, 10, 12, 2),
            (10, 11, 13, 2), (11, 3, 14, 2), (3, 15, 17, 3),
            (15, 16, 18, 3), (16, 4, 19, 3), (4, 20, 22, 4),
            (20, 21, 23, 4), (21, 1, 24, 4)]

    where the last index, rather than being an edge label (based on
    GMSHs line attribute), is -1.
    """
    # Keep track of the original edges. At the end return the non-sorted edges.
    original_order = dict()
    sorted_edges = set()

    side_nodes = int((len(elements[0])-3)/3)
    for element in elements:
        # Guarantee uniqueness of edges by sorting the nodes. At the end map
        # the sorted versions back to the original versions.
        local_edges = [(element[i], ) + (element[j], ) +
                       tuple(element[3+i*side_nodes:3+(i+1)*side_nodes])
                       + (-1, ) for (i, j) in [(0, 1), (1, 2), (2, 0)]]

        local_sorted_edges = [tuple(sorted(t)) for t in local_edges]
        original_order.update(zip(local_sorted_edges, local_edges))

        for edge in local_sorted_edges:
            if edge in sorted_edges:
                sorted_edges.remove(edge)
            else:
                sorted_edges.add(edge)

    return [original_order[t] for t in sorted_edges]