Numerical simulations of galaxy formation provide a powerful technique for calculating the non-linear evolution of cosmic structure formation. In fact, they have played an instrumental role in establishing the current standard cosmological model known as LCDM. However, unlocking the predictive power of current petaflop and future exaflop computing platforms requires a targeted effort in developing new numerical methods that excel in accuracy, parallel scalability, and in physical fidelity to the processes relevant in galaxy formation. A new moving-mesh technique for hydrodynamics recently developed by us provides a significant opportunity for a paradigm shift in cosmological simulations of structure formation, replacing the established smoothed particle hydrodynamics technique with a much more accurate and flexible approach. Building on the first successes with this method, we here propose a comprehensive research program to apply this novel numerical framework in a new generation of hydrodynamical simulations of galaxy formation that aim to greatly expand the physical complexity and dynamic range of theoretical galaxy formation models. We will perform the first simulations of individual galaxies with several tens of billion hydrodynamical resolution elements and full adaptivity, allowing us to resolve the interstellar medium in global models of galaxies with an unprecedented combination of spatial resolution and volume. We will simultaneously and self-consistently follow the radiation field in galaxies down to very small scales, something that has never been attempted before. Through cosmological simulations of galaxy formation in representative regions of the Universe, we will shed light on the connection between galaxy formation and the large-scale distribution of gas in the Universe, and on the many facets of feedback processes that regulate galactic star formation, such as energy input from evolving and dying stars or from accreting supermassive black holes.