A backwards-tracking Lagrangian–Eulerian method is used to simulate planar viscoelastic squeeze flow. The momentum and continuity equations are discretized with the finite volume method and implicit immersed boundary conditions are used to describe objects in the domain. The viscoelastic squeeze flow, which involves moving solid geometry as well as free surface flow, is chosen for its relevance in industrial applications, such as adhesive parts assembly and hemming. The main objectives are to validate the numerical method for such flows and to outline the grid resolution dependence of important flow quantities. The main part of the study is performed with the Oldroyd-B model, for which the grid dependence is assessed over a wide range of Weissenberg numbers. An important conclusion is that the load exerted on the solids can be predicted with reasonable accuracy using a relatively coarse grid. Furthermore, the results are found to be in excellent agreement with theoretical predictions as well as in qualitative resemblance with numerical results from the literature. The effects of different viscoelastic properties are further investigated using the PTT model, revealing a strong influence of shear-thinning for moderate Weissenberg numbers. Finally, a reverse squeeze flow is simulated, highlighting important aspects in the context of adhesive joining applications.