just to develop what Fedja said above:
say you have a map $\phi(x,y)=(u(x,y),v(x,y),w(x,y))$ from $D = [a;b]^2$ to $\mathbb{R}^3$ that represents a surface (Klein bottle in your initial example). The infinitesimal surface $(x;x+dx) \times (y;y+dy)$ is mapped to a surface of area $|\partial_x \Phi \wedge \partial_y \Phi| dx dy = A(x,y) dx dy$. So it would suffice to sample points inside $D$ distributed according to a probability distribution $p(x,y) \propto A(x,y) dx dy$.
1: as Fedja mentioned it, you can use a rejection sampling approach - you just need to find a bound on $\|A\|_{\infty}$, which should be very easy since $D$ is simple.
2: you can use a MCMC approach, which is equally easily implemented (e.g. independence sampler) - might be a little bit faster.
