[stokes_solver] Fix mass matrix for linear CVFE and allow direct solver for P

The mass matrix in the Stokes solver for linear CVFE was only an approximation and is now the correct mass matrix (which also has some off-diagonal entries). This reduces slightly the number of iterations.

For linear CVFE it's possible to implement more efficiently, so I suggest adding the code for higher order schemes in a different MR (if needed e.g. for P3P2).