docs: def faces_per_edge(self):
          """Returns an Ex2 array of adjacencies between faces, where
          each element in the array is a face index. Each edge is included
          only once. Edges that are not shared by 2 faces are not included."""
          import numpy as np
          import scipy.sparse as sp
          from blmath.numerics.matlab import col
          IS = np.repeat(np.arange(len(self.f)), 3)
          JS = self.f.ravel()
          data = np.ones(IS.size)
          f2v = sp.csc_matrix((data, (IS, JS)), shape=(len(self.f), np.max(self.f.ravel())+1))
          f2f = f2v.dot(f2v.T)
          f2f = f2f.tocoo()
          f2f = np.hstack((col(f2f.row), col(f2f.col), col(f2f.data)))
          which = (f2f[:, 0] < f2f[:, 1]) & (f2f[:, 2] >= 2)
          return np.asarray(f2f[which, :2], np.uint32)