A finite element analysis system is developed to investigate the viscoelastic deformation of polymer during hot embossing process, employing a nonlinear viscoelastic model as a polymer constitutive model. The interface between air and polymer is captured by a level set method. The viscoelastic flow problem is solved using a discrete elastic-viscous stress splitting (DEVSS) formulation along with a matrix logarithm of the conformation tensor. A discontinuous Galerkin formulation is employed to treat convective problems. The developed method is applied to the filling stage of hot embossing for a two-dimensional cavity. The details of the polymer viscoelastic deformation are investigated in terms of the conformation tensor distribution at the end of the embossing. The embossing speed and the cavity aspect ratio are found to have significant effects on the polymer conformation development in the molded part. As for the difference in filling pattern between the Newtonian viscous fluid and the viscoelastic fluid, the difference grows with time, but not significant in the particular geometry and processing conditions chosen.