Classical Finite Volume methods for multi-dimensional problems include stabilization (e.g. via a Riemann solver), that is derived by considering several one-dimensional problems in different directions. Such methods therefore ignore a possibly existing balance of contributions coming from different directions, such as the one characterizing multi-dimensional stationary states. Instead of being preserved, they are usually diffused away by such methods. Stationarity preserving methods use a better suited stabilization term that vanishes at the stationary state, allowing the method to preserve it. This work presents a general approach to stationarity preserving Finite Volume methods for nonlinear conservation/balance laws. It is based on a multi-dimensional stationarity preserving quadrature strategy that allows to naturally introduce genuinely multi-dimensional numerical fluxes. The new methods are shown to significantly outperform existing ones even if the latter are of higher order of accuracy and even on non-stationary solutions. These methods are by construction defined on Cartesian structured grids, hence, to extend the method to more complex geometries, we introduce a stable immersed boundaries techniques, that allows the application on any geometry, without losing the stationarity preserving property.