docs: def _compute_edges_cells(self):
          """This creates interior edge->cells relations. While it's not
          necessary for many applications, it sometimes does come in handy.
          """
          if self.edges is None:
              self.create_edges()
  
          num_edges = len(self.edges["nodes"])
  
          counts = numpy.zeros(num_edges, dtype=int)
          fastfunc.add.at(
              counts,
              self.cells["edges"],
              numpy.ones(self.cells["edges"].shape, dtype=int),
          )
  
          # <https://stackoverflow.com/a/50395231/353337>
          edges_flat = self.cells["edges"].flat
          idx_sort = numpy.argsort(edges_flat)
          idx_start, count = grp_start_len(edges_flat[idx_sort])
          res1 = idx_sort[idx_start[count == 1]][:, numpy.newaxis]
          idx = idx_start[count == 2]
          res2 = numpy.column_stack([idx_sort[idx], idx_sort[idx + 1]])
          self._edges_cells = [
              [],  # no edges with zero adjacent cells
              res1 // 3,
              res2 // 3,
          ]
          # self._edges_local = [
          #     [],  # no edges with zero adjacent cells
          #     res1 % 3,
          #     res2 % 3,
          #     ]
  
          # For each edge, store the number of adjacent cells plus the index into
          # the respective edge array.
          self._edge_gid_to_edge_list = numpy.empty((num_edges, 2), dtype=int)
          self._edge_gid_to_edge_list[:, 0] = count
          c1 = count == 1
          l1 = numpy.sum(c1)
          self._edge_gid_to_edge_list[c1, 1] = numpy.arange(l1)
          c2 = count == 2
          l2 = numpy.sum(c2)
          self._edge_gid_to_edge_list[c2, 1] = numpy.arange(l2)
          assert l1 + l2 == len(count)
  
          return