Usually the divergence-free constraint is enforced in a discrete sense, i.e., with respect to some finite dimensional set of pressure test functions:
where \mathcal{Q}^h is the discrete pressure test space. A velocity satisfying such a condition can be solved for alongside the corresponding pressure by using a mixed element, as in the Stokes demo here:
which can be modified in a straightforward way to include the nonlinear convection term for the Navier–Stokes problem. Some care must be taken to choose the velocity and pressure spaces in a stable way, though. Alternatively, one can include stabilization in the formulation, with the most popular choice being the so-called pressure-stabilizing Petrov–Galerkin (PSPG) method (typically in conjunction with streamline-upwind Petrov–Galerkin (SUPG) when there is convection). With PSPG, you can use whatever pressure and velocity spaces you want, and the resulting equation system tends to be friendlier to iterative linear solvers, but this modifies the approximate div-free condition satisfied by \mathbf{u}^h (in a way that still converges toward a div-free exact solution under refinement).
It is possible to get an exactly (pointwise) div-free \mathbf{u}^h, if \mathcal{V}^h and \mathcal{Q}^h are chosen such that the divergence of every function in \mathcal{V}^h is in \mathcal{Q}^h. This can be done using standard finite elements in FEniCS as outlined here
but is subject to some special requirements on the mesh. Exactly div-free solutions can also be obtained using isogeometric spline spaces, which are accessible in FEniCS using a library that I wrote. In isogeometric analysis, you can also solve directly for a vector potential whose curl is the velocity, such that the velocity is div-free by construction. However, exactly div-free approaches are not currently widely-used.