A detailed and systematic analysis is performed on the local and global properties of the recently developed Harmonic Polynomial Cell (HPC) method, a very accurate and efficient field solver for problems governed by the Laplace equation. At the local cell level, a simple rule is identified for the proper choice of harmonic polynomials in the local representation of the velocity potential in cells with symmetry properties. The local solution error, its convergence rate, its dependence on the cell topology, its distribution inside the cell and its features across cells with different dimensions, are carefully examined with relevant findings for HPC numerical implementations. At the global level, the error convergence rate is analytically estimated in terms of error contributions from the boundary conditions and from inside the liquid domain. In most cases, the error associated with boundary conditions dominates the global error. In order to minimize it, Quadtree grid strategies or high-order local expressions of the velocity potential are proposed for cells near critical boundary portions. To model accurately the boundary conditions on rigid or deformable surfaces with generic geometries, three different grid strategies are proposed by adopting concepts of immersed boundary method and overlapping grids. They are comparatively studied for a circular rigid cylinder in infinite fluid and for the propagation of a free-surface wave. Then, an immersed boundary strategy, using numerical choices suggested in this paper, is successfully compared against a fully nonlinear Boundary Element Method for the case of a surface-piercing circular cylinder heaving in otherwise calm water.