docs: def jac_uniform(X, cells):
      """The approximated Jacobian is
  
        partial_i E = 2/(d+1) (x_i int_{omega_i} rho(x) dx - int_{omega_i} x rho(x) dx)
                    = 2/(d+1) sum_{tau_j in omega_i} (x_i - b_{j, rho}) int_{tau_j} rho,
  
      see Chen-Holst. This method here assumes uniform density, rho(x) = 1, such that
  
        partial_i E = 2/(d+1) sum_{tau_j in omega_i} (x_i - b_j) |tau_j|
  
      with b_j being the ordinary barycenter.
      """
      dim = 2
      mesh = MeshTri(X, cells)
  
      jac = numpy.zeros(X.shape)
      for k in range(mesh.cells["nodes"].shape[1]):
          i = mesh.cells["nodes"][:, k]
          fastfunc.add.at(
              jac,
              i,
              ((mesh.node_coords[i] - mesh.cell_barycenters).T * mesh.cell_volumes).T,
          )
  
      return 2 / (dim + 1) * jac