This contribution outlines a computational framework for the analysis of flexible multibody dynamics contact problems. The framework combines a flexible body formalism, specifically, the absolute nodal coordinate formulation (ANCF), with a discrete continuous contact force model to address many-body dynamics problems, i.e., problems with hundreds of thousands of rigid and deformable bodies. Since the computational effort associated with these problems is significant, the analytical framework is implemented to leverage the computational power available on today's commodity graphical processing unit (GPU) cards. The framework developed is validated against commercial and research finite element software. The robustness and efficiency of this approach is demonstrated through numerical simulations. The resulting simulation capability is shown to result in 2 orders of magnitude shorter simulation times for systems with a large number of flexible beams that might typically be encountered in hair or polymer simulations.