## 2 Methodology  
Finite Element discretization of a structure with linear elastic material undergoing large displacement/rotation and frictional contact leads to a second-order ordinary differential equation (see, e.g., Refs. [22; 28; 29; 30; 31] for more details):  
\[\mathbf{M}\bar{\mathbf{u}}(t)+\mathbf{C}\dot{\mathbf{u}}(t)+\mathbf{f}_{int}( \mathbf{u}(t))=\mathbf{f}_{c}(t)+\mathbf{f}_{c}(\dot{\mathbf{u}},\mathbf{u},t) \tag{1}\]  
in which \(\mathbf{M}\in\mathbb{R}^{N\times N}\) and \(\mathbf{C}\in\mathbb{R}^{N\times N}\) represent the mass and viscous damping matrices of the structure and \(\mathbf{u}(t)\in\mathbb{R}^{N}\) is the vector of nodal displacement. The force vectors are the elastic internal forces \(\mathbf{f}_{int}:\mathbb{R}^{N}\mapsto\mathbb{R}^{N}\), the linear excitation \(\mathbf{f}_{c}\in\mathbb{R}^{N}\), and the nonlinear contact forces \(\mathbf{f}_{c}:\mathbb{R}^{N}\mapsto\mathbb{R}^{N}\). The vector of internal elastic force \(\mathbf{f}_{int}\) accounts for both the linear and the nonlinear components. In the case of a geometrically linear structure, it is a linear function of the nodal displacement vector \(\mathbf{f}_{int}(t)=\mathbf{K}\mathbf{u}(t)\), where \(\mathbf{K}\) is the stiffness matrix. This study considers the mass proportional damping as \(\mathbf{C}=\mathbf{\alpha}\mathbf{M}\).  
The solution of Equation (1) for a full-order model provides accurate results, but its computing cost is significant if the size of the system \(N\) is large. A model order reduction technique is employed to mitigate this cost by reducing the equation size. The linear projection of Equation (1) on a reduced space is performed by introducing the coordinate transformation:  
\[\mathbf{u}(t)=\mathbf{R}\mathbf{q}(t) \tag{2}\]where \(\mathbf{R}\in\mathbb{R}^{N\times n}\), the so-called projection matrix is constructed by a set of orthogonal vectors that represent the reduced space basis. Time-dependent vector \(\mathbf{q}\in\mathbb{R}^{n}(n\ll N)\) consists of the reduced order coordinates. Upon substitution of Equation (2) in Equation (1) and premultiplying \(\mathbf{R}^{\mathsf{T}}\) in the resulting equation, the following reduced-order equation is obtained using Galerkin projection, which requires the residual to be orthogonal to the subspace span:  
\[\mathbf{\tilde{M}}\mathbf{\dot{q}}(t)+\mathbf{\tilde{C}}\mathbf{\dot{q}}(t)+ \mathbf{\tilde{f}}_{\text{int}}(\mathbf{q}(t))=\mathbf{\tilde{f}}_{\text{e}}( t)+\mathbf{\tilde{f}}_{\text{c}}(\mathbf{\dot{q}},\mathbf{q},t) \tag{3}\]  
in which the matrices \(\mathbf{\tilde{M}}\in\mathbb{R}^{n\times n}\) and \(\mathbf{\tilde{C}}\in\mathbb{R}^{n\times n}\) represent the reduced mass and damping matrices of the structure, respectively, where \(\mathbf{\tilde{A}}=\mathbf{R}^{\mathsf{T}}\mathbf{A}\mathbf{R}\). The generalized force vectors, \(\mathbf{\tilde{f}}\), are also computed as \(\mathbf{\tilde{f}}=\mathbf{R}^{\mathsf{T}}\mathbf{f}\).  
The linear structural matrices \((\mathbf{M},\mathbf{C})\) and the reduced contact force \(\mathbf{f}_{\text{c}}\) are projected before the solution of the equations in the offline stage when the reduction basis is also computed. On the contrary, the projection of the elastic forces \(\mathbf{f}_{\text{int}}\) has to be performed in the online stage at each iteration of the solution of Equation (3) based on the predicted full displacement field.