A generic projection maps one vector to another such that their difference is a gradient field and the projected vector does not have to be solenoidal. By studying the commutator of Laplacian and the generic projection, the incompressible Navier-Stokes equations with no-slip conditions are reformulated as the sole evolution of a divergent velocity with the incompressibility constraint enforced by a pressure Poisson equation. This GePUP formulation yields numerical methods with distinguishing features as follows. (1) The velocity divergence is governed by a heat equation, and thus bounded by the maximum principle. (2) The method is fourth-order accurate both in time and in space. (3) The linear systems are solved with geometric multigrid, yielding optimal complexity. (4) The time integrator is treated as a black box in the framework of method-of-lines so that changing from explicit time integration to implicit-explicit time integration is trivial, making the formulation applicable to flows with both low and high Reynolds numbers. This finite-volume solver is further augmented with parallel computing and adaptive mesh refinement. We have generalized the solver to domains with irregular boundaries and are working on domains with moving boundaries.