Let $E$ be a metric space, $\lambda$ be a measure on $\mathcal B(E)$ (in my application $E=[0,1)^2$ and $\lambda$ is the Lebesgue measure on $\mathcal B(\mathbb R^2)$ restricted to $E$), $p:E\to[0,\infty)$ with $$c:=\int p\:{\rm d}\lambda\in(0,\infty)$$ and $\varepsilon>0$.
Are we able to sample a sequence $X_1,\ldots,X_N$ of random variables such that $$d(X_m,X_n)\ge\varepsilon\;\;\;\text{for all }m,n\in\{1,\ldots,N\}\tag1$$ and $$X_i\sim\frac pc\lambda\tag2\;\;\;\text{for all }i\in\{1,\ldots,N\}$$ or at least the distribution of $X_n$ converges (in a suitable sense) to $\frac pc\lambda$?
An obvious approach would be to sample $(X_i)_{i\in\mathbb N}$ independent and identically distributed and discard all samples which do not satisfy the constraint $(1)$. However, is the resulting process still distributed according to $(2)$? (This question is related.)
EDIT: As far as I can tell from this question, it should at least be possible to sample from some distributions (in the question it is a normal distribution) with constraints. So, the generation I'm asking for might be possible as well ...
EDIT 2: $N$ does not need to be arbitrary here. I would to find the largest $N$ possible satisfying $(1)$ and either $(2)$ or $(3)$.
This is clearly not possible for arbitrary $N$, as when $N$ is large enough, there is no such configuration at all. The choice of $p$ will also affect the the maximum possible $N$, as the more peaked the distribution, the fewer points will "fit." For instance, taking $p(x) = \mathbf{1}_{0 \leq x < \epsilon}$ will only allow $N = 4$.
In the special case of the uniform distribution $p \equiv 1$, you can drop a hexagonal lattice with edge length $\epsilon$ uniformly at random over $[0, 1)^2$ (this can be done by sampling a base point uniformly at random from $[0, 1)^2$ and then choosing an edge direction uniformly at random, and this will determine the whole lattice).
Now, $X_n$ can be chosen by sampling vertices from the lattice without replacement.