The Stokes system with constant viscosity can be cast into different formulations by exploiting the incompressibility constraint. For instance, the rate of strain tensor in the weak formulation can be replaced by the velocity-gradient yielding a decoupling of the velocity components in the different coordinate directions. Consequently, the discretization of this partly decoupled formulation leads to fewer nonzero entries in the stiffness matrix. This is of particular interest in large scale simulations where a reduced memory bandwidth requirement can help to significantly accelerate the computations. In the case of a piecewise constant viscosity, as it typically arises in multi-phase flows, or when the boundary conditions involve traction, the situation is more complex, and one has to treat the cross derivatives in the original Stokes system with care. A naive application of the standard vectorial Laplacian results in a physically incorrect solution, while formulations based on the rate of strain tensor increase the computational effort globally. Here, we propose a new approach that is consistent with the stress-divergence formulation and preserves the decoupling advantages of the velocity-gradient-divergence formulation in isoviscous subdomains. The modification is equivalent to locally changing the discretization stencils at interfaces or boundaries. Hence, the more expensive discretization is of lower complexity, making the additional computational cost in large scale simulations negligible. We establish consistency and convergence properties and show that in a massively parallel setup, the multigrid solution of the resulting discrete systems is faster than for the classical stress-divergence formulation.