This thesis considers the numerical simulation of rigid and deformable bodies, as well as compressible fluids. We consider each of these types of simulations independently, and in particular we focus on what it takes to make these simulations both efficient and scalable. First, we develop a robust parallelized method for simulating cloth and we demonstrate simulations consisting of up to 2 million triangles. This added level of detail allows us to achieve high detailed folds and wrinkles. We propose a robust history-based repulsion/collision framework where repulsions are treated accurately and efficiently on a per time step basis. Distributed memory parallelism is used for both time evolution and collisions and we specifically address Gauss-Seidel ordering of repulsion/collision response. Next, we propose a method for alleviating the stringent CFL condition imposed by the sound speed in simulating inviscid compressible flow with shocks, contacts and rarefactions. Our method is based on the pressure evolution equation, so it works for arbitrary equations of state, chemical species, etc. The relaxed CFL condition allows us to simulate shocks, contacts and rarefactions accurately while taking much larger time steps than before. Then, we turn to the simulation of rigid bodies, where we present an algorithm for conserving energy and momentum when advancing rigid body orientations. Furthermore, we develop a technique for clamping energy gain during contact and collisions. Together, these methods allow us to prevent energy increase during rigid body simulations, regardless of the time step size. This allows us to reduce the computation needed while still producing stable and physically plausible simulations. We also introduce a technique for fast and realistic fracture of rigid bodies using a novel collision-centered prescoring algorithm. Finally, we extend the use of energy preservation techniques to the simulation of deformable bodies, again with the goal of reducing the cost of these simulations. We propose a new spring that, in one spatial dimension, gives the exact solution regardless of the size of the time step chosen. In multiple spatial dimensions, the problem becomes nonlinear because the direction of the spring changes over time, and thus we propose an iterative approach. Then, we consider the simulation of more complicated elements such as triangles, tetrahedra, and finally full meshes and propose a novel technique that allows us to cut the iterative approach short and instead apply a final correction globally to the mesh.