This paper presents a series of smoothed finite element methods (S-FEM) to solve transcranial stimulation problems. The analysis domain is first discretized into a set of four-node tetrahedral elements, and linear shape functions are employed to interpolate the field variables. Then, the smoothing domains are further formed with respect to the nodes, edges or faces of elements. In order to improve the accuracy of lower order interpolation, the gradient field of the problem is smoothed using gradient smoothing operations over each smoothing domain. Based on the generalized smoothed Galerkin weakform, the discretized system equations are finally obtained. Numerical examples, including transcranial direct current stimulation (tDCS) and transcranial magnetic stimulation (TMS) problems, demonstrate that the S-FEM method possesses the following important properties: (1) temporal stability; (2) super accuracy and super convergence; (3) higher computational efficiency;