This paper describes a numerical method for solving the (incompressible) Navier-Stokes equations that is based on the idea that the motivation for discretizing differential operators should be to mimic their fundamental conservation and dissipation properties. Therefore, the symmetry of the underlying differential operators is preserved. The resulting discretization is stable on any grid. Its accuracy is tested for a turbulent channel flow at Re-tau =180 by comparing the results to those of physical experiments and previous numerical studies. The method is generalized to compute flows in domains with arbitrarily-shaped boundaries, where the boundary is represented using the Cartesian grid approach. To that end, a novel cut-cell discretization has been developed. The boundary treatment is successfully tested for flow around a circular cylinder.