We describe a generic infrastructure for time evolution simulations in numerical relativity using multiple grid patches. After a motivation of this approach, we discuss the relative advantages of global and patch-local tensor bases. We describe both our multi-patch infrastructure and our time evolution scheme, and comment on adaptive time integrators and parallelization. We also describe various patch system topologies that provide spherical outer and/or multiple inner boundaries. We employ penalty inter-patch boundary conditions, and we demonstrate the stability and accuracy of our three-dimensional implementation. We solve both a scalar wave equation on a stationary rotating black hole background and the full Einstein equations. For the scalar wave equation, we compare the effects of global and patch-local tensor bases, different finite differencing operators and the effect of artificial dissipation onto stability and accuracy. We show that multi-patch systems can directly compete with the so-called fixed mesh refinement approach; however, one can also combine both. For the Einstein equations, we show that using multiple grid patches with penalty boundary conditions leads to a robustly stable system. We also show long-term stable and accurate evolutions of a one-dimensional nonlinear gauge wave. Finally, we evolve weak gravitational waves in three dimensions and extract accurate waveforms, taking advantage of the spherical shape of our grid lines.