Solid-state high-order harmonic generation (HHG) continues to attract a lot of interest. From the theory and simulation standpoint, two issues are still open; The first is the so-called transition-dipole phase problem. It has been recognized that the dipoles must be treated as complex-valued quantities and that their corresponding Berry connections must be included to ensure phase-gauge invariance. However, while this has been successfully implemented for lower-dimensional systems, fully vectorial and three-dimensional simulations remain challenging. The second issue concerns the symmetry of the high-order harmonic response, when simulations sometimes fail to honor the symmetry of the crystalline material. This work addresses both of these problems with the help of a HHG-simulation approach which (a) is manifestly free of the transition-dipole phase problem, (b) does not require calculation of dipole moments, (c) can account for the contributions from the entire Brillouin zone, and (d) faithfully preserves the symmetry of the simulated crystalline material. We use the method to show that high-order harmonic sources are distributed throughout the Brillouin zone with various phase shifts giving rise to significant cancellations. As a consequence, for the simulated response to correctly capture the material symmetry, contributions from the entire Brillouin zone must be included. Our results have important implications for a number of HHG applications, including all-optical band and dipole reconstruction.