In principle the inverse transform method can be used for discrete distributions as well. This boils down to first subdividing the interval $[0,1]$, where the sub-intervals are formed according to the distribution function $F_X$ of the target random variable $X$. Then a standard uniform random variable $U$ is drawn and one assigns a value to $X$ depending on the sub-interval in which $U$ falls. A more technical explanation can for example be found in < (or most introductory probability / statistics text books). As noted in the pdf this is not always the most convenient way to go about (for example in the Poisson case).
Once a sampling algorithm based on a uniform random variable is available a multivariate observation can be generated via copulas as usual. First a pseudo-realization $\mathbf{U} = (U_1,\ldots,U_d)$ is drawn from the copula. Then you generate $X_i$ from $U_i$ (as described above) to get $\mathbf{X} = (X_1,\ldots,X_d)$.