In the quest for a numerical method for surface waves and wave‐induced effects applicable when linear or weakly nonlinear methods are insufficient, a three‐dimensional numerical wave tank assuming fully‐nonlinear potential‐flow theory is proposed. When viscous‐flow effects, breaking waves or other violent flow‐phenomena are not of primary importance, potential‐flow methods may have similar capability in capturing the involved physics as Navier–Stokes solvers while being potentially more accurate in handling wave‐propagation mechanism and more computationally efficient. If made sufficiently accurate, efficient and numerically robust, fully‐nonlinear potential flow models can therefore represent a powerful tool in the study of ocean waves and their interaction with marine structures, which is the main motivation behind the present work. The governing Laplace equation for the velocity potential is solved using the harmonic polynomial cell method, which is a field method giving high‐order accuracy provided that the cells used to describe the water domain have no stretching or distortion. This can only be achieved in a grid with cubic cells, which leads to poor numerical efficiency unless measures are introduced to refine the grid locally. Here, to improve the efficiency using strictly cubic cells, an adaptive grid refinement technique is introduced. It is shown that this has the ability to improve the computational speed with a factor of up to 20 without sacrificing accuracy. Numerical results are shown to be in good agreement with highly accurate nonlinear reference solutions for regular and irregular waves of various steepness up to the limit of theoretical wave breaking. For long‐crested irregular waves, significant discrepancies with a second‐order theory for the crest‐height distribution are identified, while the second‐order theory appears to provide a better description of the crest height for the single short‐crested irregular sea state simulated. Having demonstrated that the proposed numerical method accurately models nonlinear wave phenomena up to the limit of wave breaking, future work should seek to implement wave‐body interaction capabilities. The adaptive grid refinement technique, which refines the grid dynamically depending on the position of boundaries of interest, is developed with this application in mind. Except from providing a robust way of dealing with wave‐body intersection points, extending the method to account for wave‐body interactions should therefore involve limited difficulty. [ABSTRACT FROM AUTHOR]