A numerical method for computing the ground state solution of Bose–Einstein condensates modeled by the Gross–Pitaevskii equation is presented. In this method, the three-dimensional computational domain is divided into hexahedral elements in which the solution is approximated by a sum of basis functions. Both polynomial and plane wave bases are considered for this purpose, and Lagrange multipliers are introduced to weakly enforce the interelement continuity of the solution. The ground state is computed by an iterative procedure for minimizing the energy. The performance results obtained for several numerical experiments demonstrate that the proposed method is more computationally efficient than similar solution approaches based on the standard higher-order finite element method.