We report the first simulations of the Faraday instability using the full three-dimensional Navier–Stokes equations in domains much larger than the characteristic wavelength of the pattern. We use a massively parallel code based on a hybrid front-tracking/level-set algorithm for Lagrangian tracking of arbitrarily deformable phase interfaces. Simulations performed in square and cylindrical domains yield complex patterns. In particular, a superlattice-like pattern similar to those of Douady & Fauve (Europhys. Lett., vol. 6, 1988, pp. 221–226) and Douady (J. Fluid Mech., vol. 221, 1990, pp. 383–409) is observed. The pattern consists of the superposition of two square superlattices. We conjecture that such patterns are widespread if the square container is large compared with the critical wavelength. In the cylinder, pentagonal cells near the outer wall allow a square-wave pattern to be accommodated in the centre.