After doing some mathematics related to the stability of elements in 3D Stokes problem I was slightly shocked to realize that $P_2-P_1$ is not stable for an arbitrary tetrahedral mesh. More precisely, in case you have an element where all nodes and three out of four facets lie on the boundary of the domain with a Dirichlet condition, you end up getting a singular matrix. This is in fact fairly trivial to conclude from the weak form of the Stokes system.
I tested the only commercial Stokes code I have an access to (COMSOL) and it allowed me to create such a mesh. Upon clicking solve I get 'Error: Singular matrix' as expected. (I am under the impression that COMSOL uses $P_2-P_1$ for its creeping flow module.)
To further test that the problem was not related to other configurations I tried the following mesh and everything works as expected.
Questions: Is this kind of constraint taken into account in (adaptive or non-adaptive) mesh generators? I see from various research papers that this element seems to be quite popular. Are these kind of boundary instabilities generally disregarded as insignificant when choosing a method to use? More importantly, what does it really mean to have a stable finite element, i.e., what kind of mesh-dependent instabilities are too much so that we conclude that the method is bad?