The capability of calculation of non-Newtonian flows: finite elements and stochastic simulation technique (CONNFFESSIT) to compute transient free surface flows is extended to deal with transient viscoelastic flows for domains in which the topology (connectivity) can change in time and also to deal with multivalued free surfaces, i.e. those that cannot be described by a single-valued height function nor be mapped onto an equivalent problem. Such flows are of great practical relevance as they are prevalent in thick-part mold filling, compression molding, multilayer extrusion and thermoforming, peeling of a pressure sensitive adhesive and numerous other processes in which a viscoelastic fluid stream can fold back onto and reconnect with itself. The method is based on an Eulerian approach for the solution of the conservation equations on a fixed grid, while molecular models are used both as stress calculators and to track the position of the free surface. The finite element scheme uses continuous piecewise elements for both velocity and pressure on an unstructured mesh. The scheme is stabilized by a Galerkin/least-squares method and is convergent. The method is quantitatively tested by comparing with previous numerical results. (C) 2003 Published by Elsevier B.V.