I know that spectral-methods which are weighted-residual methods can be applied to solve for the incompressible N-S equations. In particular, they are applied to Direct Numerical Simulations (DNS) for studying turbulence. However, I'm still unclear how the continuity condition is enforced in such cases, and what trial functions are chosen.
For incompressible spectral DNS, often you'll see reference to John Kim, Parviz Moin and Robert Moser (1987): Turbulence statistics in fully developed channel flow at low Reynolds number, Journal of Fluid Mechanics, 177, pp 133-166, which includes details on how continuity is enforced. You might also find the book Spectral Methods in Fluid Dynamics by Canuto et al., Springer 1988 of interest.
Frequently some sort of velocity-pressure splitting is employed which gives a step in the method where the velocity and pressure are coupled directly.