The use of the penalty method for the finite element solution of incompressible flow problems has attracted the attention of many researchers. Its popularity is mainly due to the fact that it reduces the number of nodal unknowns of the problem. On the other hand, if the pure incompressible mixed velocity-pressure formulation is adopted, the element stiffness matrix of the discrete finite element problem has zero diagonal terms. Most of the standard finite element direct solvers use the elements of the diagonal as pivots |31]. In this case, renumbering algorithms have to be modified in order to prevent the appearence of zeroes in the diagonal of the assembled matrix when reducing the equations in the solution process. This results in the increase of the bandwidth of the global matria. This problem is circumvented if the incompressibility constraint is penalised.