In this paper, we introduce an \(H({\text {div}})\) finite element method on polygonal and polyhedral meshes for solving the Stokes equations in the primary velocity–pressure formulation. An \(H({\text {div}})\) finite element on polygons or polyhedra is introduced to approximate the velocity so that the method is pressure robust and produces exact divergence free solutions, when the discontinuous \(P_k\) finite element is adopted for the pressure approximation. In addition, this method is stabilizer-free, in handling the \(H^1\) non-conformity. Optimal order error estimates are established for the method. Numerical tests are conducted to confirm the theory.
