A spectral element method (SEM) is introduced for accurate calculation of band structures of three-dimensional anisotropic photonic crystals. The method is based on the finite-element framework with curvilinear hexahedral elements. Gauss-Lobatto-Legendre polynomials are used to construct the basis functions. In order to suppress spurious modes, mixed-order vector basis functions are employed and the Bloch periodic boundary condition is imposed into the basis functions with tangential components at the boundary by multiplying a Bloch phase factor. The fields and coordinates in the curvilinear hexahedral elements are mapped to the reference domain by covariant mapping, which preserves the continuity of tangential components of the field. Numerical results show that the SEM has exponential convergence for both square-lattice and triangular-lattice photonic crystals. The sampling density as small as 3.4 points per wavelength can achieve accuracy as high as 99.9%. The band structures of several modified woodpile photonic crystals are calculated by using the SEM.