Accurately modeling the interaction between landslides and flexible barriers in complex terrain remains challenging due to the limitations of conventional solvers in geometric fidelity, contact accuracy, and computational efficiency. This study presents an enhanced Material Point Method (MPM) framework coupled with a Discrete Element Method (DEM) representation of flexible barriers to resolve the large-deformation dynamics of landslides and the large-strain responses of flexible barriers. Flexible barrier components are discretized as DEM particles with predefined topological constraints, while a pre-contact list is generated to eliminate redundant contact searches that hinder GPU parallelization. Additionally, a grid-based Hertz-Mindlin contact formulation is implemented to robustly capture the interactions between MPM particles and DEM elements, significantly improving computational efficiency. The framework also leverages substratum virtualization technology and a digital terrain model (DTM) to efficiently represent complex topographies. Validation against benchmark tests, including rockfall impacts on flexible ring nets and dry debris flows impacting flexible barriers, demonstrates strong agreement with experimental observations. Its applicability to large-scale engineering problems is further verified through the Ganluo landslide along the ChengKun railway. In this case study, flexible barriers placed within the narrow gully at the landslide source dissipate approximately 50 % of the kinetic energy of landslide. Moreover, the global GPU acceleration achieved approximately a twofold improvement in computational efficiency compared with existing CPU-GPU frameworks. Overall, this framework provides a reliable and efficient computational tool for the design, optimization, and performance assessment of flexible protective structures subjected to complex landslide impact scenarios.