A direct numerical simulation of the three-dimensional elektrokinetic instability near a charge selective surface (electric membrane, electrode, or system of micro-/nanochannels) is carried out and analyzed. A special finite-difference method was used for the space discretization along with a semi-implicit $3frac{1}{3}$-step Runge-Kutta scheme for the integration in time. The calculations employed parallel computing. Three characteristic patterns, which correspond to the overlimiting currents, are observed: (a) two-dimensional electroconvective rolls, (b) three-dimensional regular hexagonal structures, and (c) three-dimensional structures of spatiotemporal chaos, which are a combination of unsteady hexagons, quadrangles and triangles. The transition from (b) to (c) is accompanied by the generation of interacting two-dimensional solitary pulses.