The Ericksen-Leslie model describes dynamics of low molar-mass nematic liquid crystals, where the spatiotemporal distribution of defects defining texture is represented by the director unit vector field ${\bf d}: \Omega_T \rightarrow \mathbb{S}^2$. It consists of the Navier-Stokes equations with an extra viscous stress tensor, and a convective harmonic map heat flow equation to govern the dynamics of the director field. Two fully discrete finite element methods, for a regularized system using the Ginzburg-Landau regularization and for the limiting Ericksen-Leslie model, are proposed, and well-posedness and related discrete energy laws are established. For the regularized model, unconditional convergence of finite element solutions towards weak solutions of the continuum model as well as convergence towards measure-valued solutions of the limiting Ericksen-Leslie model are verified when the mesh parameters and the regularization parameter successively tend to zero. Computational experiments are also presented to show the importance of balancing numerical and regularization parameters, to compare regularized and direct approaches, and to show numerically the finite-time formation, annihilation, and evolution of point defects.