Contact analysis is performed within each iteration until the equilibrium condition is satisfied. The residual vector that results from the global equilibrium is computed before the contact analysis starts. This vector is used to compute the contact forces resulting from the contact inequality constraints applied on the system. The constraint matrix is constructed and assembled to the left-hand-side coefficient matrix in each iteration. This would require the entire coefficient matrix to be reformulated and decomposed after each iteration, which is computationally expensive. Alternatively, using Crout factorization technique and isolating the Lagrange multiplier degrees of freedom at the end of the left-hand-side coefficient matrix, makes this action unnecessary. Only the part of the coefficient matrix corresponding to the Lagrange multiplier degrees of freedom is assembled and decomposed in each iteration without having to reconstruct and decompose the entire coefficient matrix. This technique saves a tremendous amount of computation.

The overall contact analysis procedure is performed, for each node on the contactor surface, as follows:

- 1.
- Read the contact forces from the residual vector
- 2.
- For the current contactor node
, perform the geometric analysis*c* - 3.
- Make the contact decision to determine the contact condition and compute the constraint matrix
- 4.
- Compute the contact force vector
- 5.
- Transfer the contact forces to the target nodes as follows
(122) (123) (124) - 6.
- Assemble the constraint matrix and force vectors to the global left-hand-side coefficient matrix and the right-hand-side load vector, respectively.