One should not confuse the problem of a particle in the real line \(\mathbb{R}\) subjected to an infinite square well potential to that of a particle confined in a circle \(S^1\). In the latter problem (particle in the circle), the Fourier modes \(\langle\theta|n\rangle=e^{inx/\hbar}/(2\pi\hbar)^{1/2}\) for \(n \in \mathbb{Z}\) are simultaneous eigenstates of the Hamiltonian and momentum operators, and they form a complete orthonormal set.
On the other hand, in the problem of a particle in the real line within an infinite square well potential, the momentum states do not need to satisfy yours boundary condition imposed on \(D_\boldsymbol{p}\), namely, that \(\psi(0)=e^{i\theta}\psi(L)\) (I have translated the boundaries of the box from \([-1/2,+1/2]\) to \([0,L]\)). The energy eigenfunctions need to have a support \(\subset [0,L]\). But the momentum states wave functions are still the ones from the whole line, namely, \(\langle x|p\rangle=e^{ipx/\hbar}/(2\pi \hbar)^{1/2} \forall p\in\mathbb{R}\), while the Hamiltonian eigenstates \(Hf_n=Ef_n\) are now given (in position representation) by the functions \(\langle x|f_n\rangle=(2/L)^{1/2}sin(2\pi n x/L)\forall x \in [0,L]\) and are zero otherwise, for \(n \in \mathbb{Z}\). In the momentum representation, these energy eigenstates reads
\[\tilde{f}_n(p)=\langle p|f_n\rangle = \int_{\mathbb{R}} dx\langle p|x\rangle\langle x|f_n\rangle\propto\int_0^L dx e^{-ipx/\hbar}sin(2\pi nx/L),\]
agreeing with Wikipedia's very much well-known result that these are just the Fourier transform of the position representation energy eigenfunctions. Just observe that now, generally, \(\tilde{f}_n(p) \neq 0\) for \(p \in \mathbb{R}\), which shows that the stationary states for yours Hamiltonian are now a superposition of uncountable many momentum states of the line. (This is needed to make the support of the energy eigenstates compact.)
From a physical viewpoint, you can understand qualitatively what is going on from the uncertainty principle. Since the particle is within \(x\in[0,L]\), the uncertainty in position is bounded from above, \(\Delta x \le L\). So, since \(\Delta p \Delta x \ge\hbar/2\) we cannot have \(\Delta p=0\), that is, we always must have an uncertainty in the momentum. Thus the energy eigenstates are indeed expected to be formed by a superposition of momentum waves.
Now you ask which one is realized in Nature. This sounds preposterous, since all this is very much idealized. The best one can do here is to confine a particle in a finite but very deep square well. You will have an amplitude for the particle leaking the well, described by an exponential decreasing term etc., but the above remarks still holds, namely, the available momentum are still a continuum. (I once again remark that this is not the same problem as that of a particle in the circle!)
Note. As an exercise, I recommend you to do the following: let \(k=2\pi n/L\) and, by taking the limit \(L \longrightarrow\infty\), \(k\) becomes a continuous variable. Then, using the exponential representation of the sine function prove that
\(\lim_{L \longrightarrow \infty}\tilde{f}_n(p)\propto \delta(k-p/\hbar)+\delta(k+p/\hbar).\)
So we recover the result that \(|p|=\hbar |k|\).