A state-of-the-art simulation technique that solves the equations of motion together with the set-valued contact and impulse laws by the time-stepping scheme of Moreau is introduced to the legged robotics community. An analysis is given that shows which of the many variations of the method fits best to legged robots. Two different methods to solve the discretized normal cone inclusions are compared: the projected over-relaxed Jacobi and Gauss-Seidel iteration. The methods are evaluated for an electrically-driven quadrupedal robot in terms of robustness, accuracy, speed and ease of use. Furthermore, the dependence of the simulation speed on the choice of the generalized coordinates is examined. The proposed technique is implemented in C++ and compared to a fast and simple approach based on compliant contact models. In conclusion, the introduced method with hard contacts is very beneficial for the simulation of legged robots.