Obtaining a satisfactory numerical solution of the classical three-dimensional drift-diffusion (DD) model, widely used in semiconductor device simulations, is still challenging nowadays, especially when the convection dominates the diffusion. In this work, we propose a series of finite element schemes with different types of averaging techniques to discretize the three-dimensional continuity equations. Our methods are based on the classical finite element framework, quite different from those mixed finite element/volume methods that also employ inverse averaging techniques. At first, the Slotboom variables are employed to transform the continuity equations into self-adjoint second-order elliptic equations with exponentially behaved coefficients. Then four averaging techniques, denoted with A1-A4, are introduced to approximate the exponential coefficient with its average on every tetrahedral element of the grid. The first scheme calculates the harmonic average of the exponential coefficient on a whole tetrahedron, and the other three schemes calculate the average of the exponential coefficient on each edge of a tetrahedral element. Our methods can avoid the spurious non-physical numerical oscillations and guarantee the conservation of the computed terminal currents with a terminal current evaluation approach. In fact, these methods not only maintain numerical stability but also overcome the disadvantages of some stabilization methods that cannot guarantee the conservation of the terminal currents, such as the streamline-upwind Petrov-Galerkin (SUPG) method. Moreover, the derivation of these discretization methods does not need the dual Voronoi grid as that of the finite volume Scharfetter-Gummel (FVSG) method or other mixed finite element/volume methods with inverse averaging techniques, which greatly reduces the complexity of parallel implementations of our methods. Simulations on two realistic three-dimensional semiconductor devices are carried out to test the accuracy and stability of our methods. According to numerical results, we conclude that scheme A4 can produce more accurate numerical solutions than the other three schemes, especially when the bias applied on the electrode is high. Numerical results also show that the scheme A4 is more robust than the Zlámal finite element method [1] in high-bias cases, and it also performs better than the FVSG method and a tetrahedral mixed finite element method [2] on poor-quality grids. Scheme A4 is also employed to study rich physical properties of the n-channel MOSFET.