def element_neighbors(self):
        """Return a list with element numbers (zero indexed) of neighboring
        elements. Note that the elements are not sorted. No spacial orientation
        can be inferred from the order of neighbors.

        WARNING: This function is slow due to a nested loop. This would be a
        good starting point for further optimizations.

        In order to speed things up, we could search using the raw data, i.e.,
        with CutMcK enabled sorting, and then restrict the loops to 2x the
        bandwidth (before - after).

        While not being returned, this function also sets the variable
        self.element_neighbors_edges, in which the common nodes with each
        neighbor are stored.

        Returns
        -------
        neighbors : list
            a list (length equal to nr of elements) with neighboring elements

        Examples
        --------


        """
        if self.element_neighbors_data is not None:
            return self.element_neighbors_data

        max_nr_edges = self.header['element_infos'][0, 2]

        # initialize the neighbor array
        self.element_neighbors_data = []
        self.element_neighbors_edges = []

        # determine neighbors
        print('Looking for neighbors')
        for nr, element_nodes in enumerate(self.elements):
            # print('element {0}/{1}'.format(nr + 1, self.nr_of_elements))
            # print(element_nodes)
            neighbors = []
            neighbors_edges = []  # store the edges to this neighbor
            for nr1, el in enumerate(self.elements):
                # we look for elements that have two nodes in common with this
                # element
                intersection = np.intersect1d(element_nodes, el)
                if intersection.size == 2:
                    neighbors.append(nr1)
                    neighbors_edges.append(intersection)
                    # stop if we reached the maximum number of possible edges
                    # this saves us quite some loop iterations
                    if len(neighbors) == max_nr_edges:
                        break
            self.element_neighbors_data.append(neighbors)
            self.element_neighbors_edges.append(neighbors_edges)
        return self.element_neighbors_data