For an unconstrained linear SPD problem, BFGS convergence matches unpreconditioned conjugate gradients. Differential operators appearing in structural mechanics problems are usually ill-conditioned with many large eigenvalues, making unpreconditioned iterative methods hopelessly slow. To deal with the troublesome spectrum, structural mechanics solvers depends on either sparse direct methods or preconditioned iterative methods (often using algebraic multigrid or domain decomposition).
In contrast, BFGS was developed in an optimization setting where most problems have a spectrum that looks like an integral operator under grid refinement, e.g., $I + K$ where $K$ is compact. For linear versions of these problems, Krylov methods would converge in a constant number of iterations, independent of resolution. The Hessian is often dense and practically inaccessible, so it is applied in unassembled form and only the identity or a diagonal scaling can be used as the "seed" Hessian $H_0$. BFGS is thus a way to improve a crude approximation to the inverse Hessian.
When solving differential problems, we can use BFGS as a lagging method, starting with a good approximation that is merely expensive to construct. This is useful when nonlinearity represents a part of the difficulty in solving a problem, but ill-conditioning is significant enough that preconditioning is needed. This turns out to work pretty well for some problems like elasticity.
Peter Brune and I wrote a short paper on this recently:
See Table II on page 9 for a large-deformation elasticity example. The standard Newton method needs 10 Jacobian evaluations, 11 residual evaluations, and 77 preconditioner applications (a linear V-cycle of algebraic multigrid; note that we can use fewer total preconditioner applications using inexact solves, but this will increase the number of Jacobian evaluations). LBFGS-6 with the same V-cycle of AMG used for $H_0$ requires 4 Jacobian evaluations, 49 residual evaluations, and 24 preconditioner applications. LBFGS lets us trade Jacobian work for more residual evaluations while retaining faster convergence than normal Jacobian lagging.
You can try this out in PETSc using only command-line options, e.g.
-snes_type qn -snes_qn_restart_type periodic -snes_qn_scale_type jacobian
The utility of this method for constrained elasticity will depend on how quickly the constraints are changing, but I expect it to be useful. We are planning to extend PETSc's support for solving constrained systems so that these methods can be used for elastic contact problems. Please let us know if you're interested in this.
