docs: def _compute_surface_areas(self, cell_ids):
          """For each edge, one half of the the edge goes to each of the end
          points. Used for Neumann boundary conditions if on the boundary of the
          mesh and transition conditions if in the interior.
          """
          # Each of the three edges may contribute to the surface areas of all
          # three vertices. Here, only the two adjacent nodes receive a
          # contribution, but other approaches (e.g., the flat cell corrector),
          # may contribute to all three nodes.
          cn = self.cells["nodes"][cell_ids]
          ids = numpy.stack([cn, cn, cn], axis=1)
  
          half_el = 0.5 * self.edge_lengths[..., cell_ids]
          zero = numpy.zeros([half_el.shape[1]])
          vals = numpy.stack(
              [
                  numpy.column_stack([zero, half_el[0], half_el[0]]),
                  numpy.column_stack([half_el[1], zero, half_el[1]]),
                  numpy.column_stack([half_el[2], half_el[2], zero]),
              ],
              axis=1,
          )
  
          return ids, vals