docs: 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