A novel framework for simulation of transient viscoelastic fluid flow is proposed. The viscoelastic stresses are calculated at Lagrangian nodes which are distributed in the computational domain and convected by the fluid. The coupling between the constitutive equation and the fluid momentum equations is established through robust interpolation with radial basis functions. The framework is implemented in a finite volume based flow solver that combines an octree background grid with immersed boundary techniques. Since the distribution of the Lagrangian node set is performed entirely based on spatial information from the fluid solver, the ability to simulate flows in complex geometries is therefore as general as for the fluid solver itself. In the Lagrangian formulation the discretization of the convective terms in the constitutive equations is avoided. No re-formulation of the constitutive equation is required for stable solutions. Numerical experiments are performed of UCM and Oldroyd-B fluids in a channel flow and of a four mode PTT fluid in a confined cylinder flow. The computed flow quantities consistently converge and agree excellently with analytical and numerical data for fully developed and transient flow.