Two primary sources of stability and accuracy problems in time-stepping schemes for simulating multibody systems are the use of polyhedral representations of smooth bodies and the approximation of the distance function (arising from the decoupling of collision detection from the solution of the dynamic time-stepping subproblem). The effects of these common approximations can be seen using a simple example of a disc rolling without slip on a flat surface. The exact solution for this problem is: the disc will roll at constant speed ad infinitum. When the disc is approximated by a uniform regular polygon the figures above show: (a) For a given step size, the energy loss decreases as the number of edges increases, to a limit (left figure) (b) For a given number of edges, the energy loss decreases with decreasing step size, to a limit (right figure). Thus, using a discretized representation of the disc always results in a loss of energy, regardless of how small you set the step size or the number of edges used in the representation. In this work, we present an implicit time-stepping scheme for simulating multibody systems with intermittent contact by incorporating the contact constraints as a set of complementarity and algebraic equations within the dynamics model. We assume each object to be a convex object described by an intersection of convex inequalities. We write the contact constraints as complementarity constraints between the contact force and a distance function dependent on the closest points on the objects. The closest points satisfy a set of algebraic constraints obtained from the KKT conditions of the minimum distance problem. These algebraic equations and the complementarity constraints taken together ensure satisfaction of the contact constraints. This enables us to formulate a geometrically implicit time-stepping scheme (i.e., we do not need to approximate the distance function) as a nonlinear complementarity problem (NCP). The resulting time-stepper is therefore more accurate; further it is the first geometrically implicit time-stepper that does not rely on a closed form expression for the distance function. We first present our approach assuming the bodies to be rigid and then extend it to locally compliant or quasi-rigid bodies. We demonstrate through example simulations the fidelity of this approach to analytical solutions and previously described simulation results. In the left figure above the horizontal line (labeled computed value) is the result obtained using our proposed scheme (which clearly shows energy conservation within a numerical tolerance).