The Scott-Vogelius \(P_k\) - \(P_{k-1}^{\hbox {\tiny disc}}\) mixed finite element is stable for \(k\ge 4\) if the underlying triangular mesh is nearly-singular vertex free. The \(P_2\) - \(P_1^{\hbox {\tiny disc}}\) and \(P_3\) - \(P_2^{\hbox {\tiny disc}}\) mixed finite elements are unstable on general triangular meshes. In this work, we enrich the \(P_k(\ge 2)\) finite element velocity space by Guzman-Neilan bubble functions, (2 \(P_2\) bubbles for \(k=2\) , 3 \(P_3\) bubbles for \(k=3\) and 0 or 1 bubble for \(k\ge 4\) on one triangle,) so that the divergence of the \(P_k\) velocity space is exactly the full discontinuous \(P_{k-1}\) space and the inf-sup condition holds on all quasiuniform triangular meshes. A Guzman-Neilan \(P_2\) bubble vector is a continuous, three-piece polynomial vector on a triangle, which vanishes on the three edges and its divergence is a one-piece \(P_1\) polynomial on the triangle. Numerical comparisons are presented showing all failing cases of the Scott-Vogelius \(P_k\) - \(P_{k-1}^{\hbox {\tiny disc}}\) finite element, while the Guzman-Neilan bubble enriched \(P_k\) - \(P_{k-1}^{\hbox {\tiny disc}}\) finite element is quasi-optimal, divergence-free and pressure-robust.
