Free boundary problems arise in various applications like coating, polymer processing, and biomedical engineering. In a problem with free boundaries, the physical domain is unknown a priori. Mesh generation equations must be added to the conservation equations to locate the free boundaries and map the unknown physical domain into a convenient reference one. Several mesh generation methods for free boundary problems have been developed, chiefly elliptic mesh generation and domain deformation. Whereas these methods are well-established for 2-D free boundary problems, 3-D free boundary problems still present challenges because of the sheer size of the computational problems, the complexity of 3-D free surface manipulation, and the need of using iterative solvers like GMRES [1]. This work applies domain deformation to three dimensional free boundary flows of viscoelastic liquids. The domain deformation method treats the mesh as a compressible solid. The momentum conservation equation for the elastic solid is coupled with the problem equations through appropriate boundary conditions to locate the free boundaries of the flow domain. The general conformation tensor model is used to represent the viscoelastic liquid [2]. The DEVSS-G/SUPG Finite Element method is applied to translate the differential equations of the problem into nonlinear algebraic equations. Fully coupled Newton's method with analytical Jacobian is employed to solve simultaneously the nonlinear equations for position, pressure, velocity, velocity gradient, and conformation. The restarted GMRES method with ILU preconditioner is parallelized by OpenMP to solve the large scale linear algebraic equations. The method is applied to study a model viscoelastic flow in a 3-D channel with a free surface section.