$\newcommand\R{\mathbb R}\newcommand\si{\sigma}\newcommand\ep{\varepsilon}\newcommand\de{\delta}$Any probability distribution on $\R$ can be approximated by a discrete probability distribution on $\R$. Any discrete probability distribution on $\R$ is a mixture of Dirac probability distributions on $\R$. The Dirac probability distribution supported at a point $a\in\R$ can be approximated by the normal distribution with a small variance centered at $a$.
Thus indeed, if the set $D$ is large enough, then any probability distribution on $\R$ can be approximated by a mixture of normal distributions $N(m,\sigma^2)$ with $(m,\si^2)\in D$.
Another way: any probability distribution $\nu$ on $\R$ can be approximated by the convolution $\nu*N(0,\si^2)$ of $\nu$ with the centered normal distribution $N(0,\si^2)$ with a small variance $\si^2$, and such a convolution is a mixture of normal distributions:
\begin{equation*}
\nu*N(0,\si^2)=\int_\R N(y,\si^2)\nu(dy). \tag{0}
\end{equation*}
If now the set $D$ is large enough to, say, contain a set of the form $C_\de\times(0,\ep)$, where $\ep\in(0,\infty)$, $\de$ is a small positive number, and the set $C_\de\subseteq\R$ is such that $\nu(C_\de)>1-\de$, then for small $\si^2\in(0,\ep)$
\begin{equation*}
\nu\approx\nu*N(0,\si^2)\approx\int_{C_\de} N(y,\si^2)\nu(dy)
=\int_D N(y,s^2)\mu(dy\times d(s^2)), \tag{1}
\end{equation*}
where $\mu(dy\times d(s^2)):=\nu(dy)1(y\in C_\de,s^2=\si^2)$ --
so that $\nu$ will be approximated by the mixture $\int_D N(y,s^2)\mu(dy\times d(s^2))$ of normal distributions $N(m,s^2)$ with $(m,s^2)\in D$.
Details: The approximate equalities in (1) can be understood as follows. Suppose that $\si_n^2\downarrow0$ and $\de_n\downarrow0$ (as $n\to\infty$). Then, by Slutsky's theorem, $\nu*N(0,\si_n^2)\to\nu$ weakly; that is, for any bounded continuous function $g\colon\R\to\R$ we have
\begin{equation*}
(\nu*N(0,\si_n^2))(g)\to\nu(g),
\end{equation*}
where $\nu(g):=\int_\R g\,d\nu$; this is how the first approximate equality in (1) can be understood.
Next, by (0), for any bounded continuous function $g\colon\R\to\R$ with $M:=\sup_{x\in\R}|g(x)|<\infty$ we have
\begin{equation*}
\begin{aligned}
&\Big|(\nu*N(0,\si_n^2))(g)-\int_{C_{\de_n}} N(y,\si_n^2)(g)\nu(dy)\Big| \\
=&\Big|\int_{\R} N(y,\si_n^2)(g)\nu(dy)-\int_{C_{\de_n}} N(y,\si_n^2)(g)\nu(dy)\Big| \\
\le&\int_{\R\setminus C_{\de_n}} N(y,\si^2)(|g|)\nu(dy) \\
\le& M\nu(\R\setminus C_{\de_n})\le M\de_n\to0;
\end{aligned}
\end{equation*}
this is how the second approximate equality in (1) can be understood.
So, for any bounded continuous function $g\colon\R\to\R$,
\begin{equation*}
\Big(\int_{C_{\de_n}} N(y,\si_n^2)\nu(dy)\Big)(g)
=\int_{C_{\de_n}} N(y,\si_n^2)(g)\nu(dy)\to\nu(g).
\end{equation*}
That is, we have the weak convergence of the normal mixtures $\int_{C_{\de_n}} N(y,\si_n^2)\nu(dy)$ to $\nu$. (The equality in the latter display is an instance of the Fubini theorem.)
