The numerical simulation of complex flows has been a subject of intense research in the last years with important industrial applications in many fields. In this paper we present a finite element method to solve the two immiscible fluid flow problems using the level set method. When the interface between both fluids cuts an element, the discontinuity in the material properties leads to discontinuities in the gradients of the unknowns which cannot be captured using a standard finite element interpolation. The method presented in this work features a local enrichment for the pressure unknowns which allows one to capture pressure gradient discontinuities in fluids presenting different density values. The method is tested on two problems: the first example consists of a sloshing case that involves the interaction of a Giesekus and a Newtonian fluid. This example shows that the enriched pressure functions permit the exact resolution of the hydrostatic rest state. The second example is the classical jet buckling problem used to validate our method. To permit the use of equal interpolation between the variables, we use a variational multiscale formulation proposed recently by Castillo and Codina (2014) [21], that has shown very good stability properties, permitting also the resolution of the jet buckling flow problem in the the range of Weissenberg number 0 < We < 100, using the Oldroyd-B model without any sign of numerical instability. Additional features of the work are the inclusion of a discontinuity capturing technique for the constitutive equation and some comparisons between a monolithic resolution and a fractional step approach to solve the viscoelastic fluid flow problem from the point of view of computational requirements. (C) 2015 Elsevier B.V. All rights reserved.