docs: def compute_curl(self, vector_field):
          """Computes the curl of a vector field over the mesh. While the vector
          field is point-based, the curl will be cell-based. The approximation is
          based on
  
          .. math::
              n\\cdot curl(F) = \\lim_{A\\to 0} |A|^{-1} <\\int_{dGamma}, F> dr;
  
          see <https://en.wikipedia.org/wiki/Curl_(mathematics)>. Actually, to
          approximate the integral, one would only need the projection of the
          vector field onto the edges at the midpoint of the edges.
          """
          # Compute the projection of A on the edge at each edge midpoint.
          # Take the average of `vector_field` at the endpoints to get the
          # approximate value at the edge midpoint.
          A = 0.5 * numpy.sum(vector_field[self.idx_hierarchy], axis=0)
          # sum of <edge, A> for all three edges
          sum_edge_dot_A = numpy.einsum("ijk, ijk->j", self.half_edge_coords, A)
  
          # Get normalized vector orthogonal to triangle
          z = numpy.cross(self.half_edge_coords[0], self.half_edge_coords[1])
  
          # Now compute
          #
          #    curl = z / ||z|| * sum_edge_dot_A / |A|.
          #
          # Since ||z|| = 2*|A|, one can save a sqrt and do
          #
          #    curl = z * sum_edge_dot_A * 0.5 / |A|^2.
          #
          curl = z * (0.5 * sum_edge_dot_A / self.cell_volumes ** 2)[..., None]
          return curl