In the spirit of the “Principle of Equipresence” introduced by Truesdell & Toupin, The Classical Field Theories (1960), we use the full version of the viscous stress tensor ν∇u+∇uT−23∇⋅uI which was originally derived for compressible flows, instead of the classical incompressible stress tensor ν∇u. (Note that, here ν is the dynamic viscosity coefficient, and u is the velocity field.) In our approach, the divergence-free constraint for the viscous stress term is not enforced ahead of discretization. Instead, our formulation allows the scheme itself to “choose” a consistent way to interpret the divergence-free constraint: i.e., the divergence-free constraint is interpreted (or enforced) in a consistent fashion in both the mass conservation equation and the stress tensor term (in the momentum equation). Furthermore, our approach preserves the original symmetrical properties of the stress tensor, e.g. its rotational invariance, and it remains physically correct in the context of compressible flows. As a result, our approach facilitates versatility and code reuse. In this paper, we introduce our approach and establish some important mathematical properties for the resulting class of finite element schemes. More precisely, for general mixed methods, which are not necessarily pointwise divergence-free, we establish the existence of a new norm induced by the full, viscous bilinear form. Thereafter, we prove the coercivity of the viscous bilinear form and the semi-coercivity of a convective trilinear form. In addition, we demonstrate L2-stability of the discrete velocity fields for the general class of methods and (by deduction) the H(div)-conforming methods. Finally, we run some numerical experiments to illustrate the behavior of the versatile mixed methods, and we make careful comparisons with a conventional H(div)-conforming scheme.