We propose a novel method to address the computational modeling of contact between rigid and flexible bodies in the large-deformation setting. The rigid bodies can be expressed with analytical formulations or discrete meshes, where different contact detection schemes are adopted to accelerate the detection process, respectively. The flexible bodies are modeled with finite element (FE) meshes, where different constitutive models can be assigned, e.g., elastic, and hyperelastic. Following the convention, the rigid surface is selected as the target surface, as it only translates or rotates but does not deform. Instead of using acceleration as a variable in multibody dynamics (MBF), we use displacement as the primal variable just as that in traditional static FEM structural modeling. Contact formulation with augmented Lagrangian multipliers is used in an updated Lagrange framework to account for both the large deformation and rotation. An unconditionally stable time integration scheme is then adopted so that the motions of rigid parts and deformations and motions of flexible parts can be evaluated seamlessly. Besides, this time integration scheme is readily used in both quasi-static and dynamic analysis. Several examples are made to demonstrate that the results obtained with the above method are accurate without the need for contact stabilization. In addition, the second-order convergence is retained with a Newton-Raphson iterative approach.