## 1 Introduction  
The accurate assessment of the ultimate bearing capacity of shallow foundations is critical for ensuring the safety, reliability, and cost-effectiveness of civil infrastructure. In geotechnical engineering, this problem has been extensively studied through various analytical and numerical methods [1]. Traditional approaches generally fall into four categories: limit equilibrium methods [2, 3, 4, 5], slip line methods [6, 7, 8], limit analysis methods [9, 10, 11, 12], and numerical methods, such as finite element or finite difference methods [13, 14, 15]. Classical bearing capacity theories developed by [PERSON] [2], [PERSON] [3], and [PERSON] [4] offer simple, closed-form solutions based on idealized assumptions. While these solutions remain widely used in practice, they are often inadequate for addressing complex conditions involving heterogeneous soils, irregular geometries, and nonlinear soil-structure interactions. The limit analysis offers a more rigorous theoretical framework by directly evaluating collapse loads without modeling the full-load deformation behavior. Based on the assumption of elastic-perfectly plastic material behavior, the upper and lower bounds of a system's load-bearing capacity depend on the chosen internal failure mechanisms and equilibrium constraints [16; 17; 18; 19]. The upper-bound theorem, in particular, estimates the maximum load a structure can sustain using kinematically admissible velocity fields and plastic dissipation. This approach has been extensively applied to geotechnical stability problems involving foundations, slopes, and retaining structures [20; 21]. In upper-bound limit analysis, it is essential to define a compatible failure mechanism whose geometry satisfies the principles of energy dissipation and external work. For cohesive soils, failure mechanisms typically comprise straight or circular slip surfaces, while for C-phi soils, mechanisms may include logarithmic spirals. The assumed failure geometry must closely reflect realistic collapse patterns to ensure that the estimated bearing capacity approximates actual behavior [22; 23; 24]. A key challenge in upper-bound limit analysis is the accurate representation of failure mechanisms. The assumed geometry must conform to realistic collapse patterns to ensure reliable capacity estimates, particularly in soils exhibiting cohesion and internal friction [22; 23; 24]. Although modern finite element limit analysis has enhanced the flexibility and accuracy of upper-bound methods, difficulties remain in capturing localized failures such as slip surfaces, cracks, and detachment at interfaces [22; 25; 26]. To overcome these limitations, zero-thickness interface elements have been introduced. These elements allow for displacement discontinuities and frictional sliding, making them well suited for modeling stress transfer across soil layers and soil-structure boundaries. When coupled with classical plasticity models, like Mohr-Coulomb, they offer a more realistic depiction of interfacial behavior under failure conditions [27; 28]. Recent advancements in automatic collapse mechanism identification have introduced robust optimization-based strategies. Discontinuity layout optimization (DLO), as proposed by [PERSON] and [PERSON] [29], offers a direct and efficient means to evaluate upper-bound solutions without assuming predefined failure paths. This method has been further advanced through open source implementations [30] and adaptive strategies [31], which allow for complex geometries and stress fields to be considered in a fully automated way. Finite element limit analysis (FELA) is a powerful numerical method used to estimate the ultimate load-bearing capacity of geotechnical systems without performing full stress-strain simulations. Based on the plasticity theory, FELA provides rigorous upper- or lower-bound solutions to collapse loads. In upper-bound FELA, kinematically admissible failure mechanisms are analyzed to compute collapse loads through optimization.