  In the last two decades, increased need for high-fidelity simulations of the
time evolution and propagation of forces in granular media has spurred renewed
interest in discrete element method (DEM) modeling of frictional contact. Force
penalty methods, while economic and accessible, introduce artificial stiffness,
requiring small time steps to retain numerical stability. Optimization-based
methods, which enforce contacts geometrically through complementarity
constraints, allow the use of larger time steps at the expense of solving a
nonlinear complementarity problem (NCP) each time step. We review the latest
efforts to produce solvers for this NCP, focusing on its relaxation to a cone
complementarity problem (CCP) and solution via an equivalent quadratic
optimization problem with conic constraints. We distinguish between linearly
convergent first order methods and second order methods, which gain quadratic
convergence and more robust performance at the expense of the solution of large
sparse linear systems.
  We propose a novel acceleration for the solution of Newton step linear
systems in second order methods using low-rank compression based fast direct
solvers. We use the Quantized Tensor Train (QTT) decomposition to produce
efficient approximate representations of the system matrix and its inverse.
This provides a robust framework to accelerate its solution in a direct or a
preconditioned iterative method. In a number of numerical tests, we demonstrate
that this approach displays sublinear scaling of precomputation costs, may be
efficiently updated across Newton iterations as well as across time steps, and
leads to a fast, optimal complexity solution of the Newton step. This allows
our method to gain an order of magnitude speedups over state-of-the-art
preconditioning techniques for moderate to large-scale systems, mitigating the
computational bottleneck of second order methods.
