Discontinuity Layout Optimization for Lower-Bound Limit Analysis: A Statically Admissible Framework for Plane Plasticity

1. Introduction

Collapse loads are often needed without the cost and path dependence of an incremental plasticity analysis, particularly when several slip mechanisms or singular zones compete. Limit analysis replaces the loading history by two bounding problems. A stress field that satisfies equilibrium, prescribed tractions, and the yield condition throughout the body gives a lower bound. A kinematically admissible velocity field for which external work equals plastic dissipation gives an upper bound [1]. Their difference measures the remaining numerical uncertainty in the collapse load.

Numerical lower-bound methods have principally developed through finite-element stress discretization. Lysmer introduced nodal stress variables for plane soil-mechanics problems [2], and Sloan formulated a linear program using three-node stress triangles, element equilibrium, interelement traction continuity, and a conservative linearization of the yield surface [3]. Subsequent nonlinear and conic formulations, improved stress elements, and adaptive remeshing increased the accuracy and robustness of this approach [4,5,6]. Their lower-bound status follows from static and yield admissibility in the continuum. Accuracy then depends on the stress approximation, the permitted locations of stress jumps, and the resolution near critical zones.

Discontinuity layout optimization (DLO) was developed from the complementary kinematic viewpoint. Smith and Gilbert connected a set of points by numerous candidate discontinuities and used linear optimization to select the critical subset [7]. Because candidate lines may cross without an additional node at every intersection, fans, singularities, and competing mechanisms can emerge from a broad search space. A distinct extension, virtual displacement-based DLO (VDLO), connects continuum stress fields and line layouts. It integrates instantaneous finite-element stresses into normal and tangential line resultants [8]. Gilbert and co-workers demonstrated the method’s geotechnical scope and extended it to rotational and three-dimensional failure mechanisms [9,10,11,12]. Within this research programme, He and collaborators developed DLO formulations and tools for slabs, practical yield-line patterns, masonry panels and walls, and coupled in-plane/out-of-plane shell failure [13,14,15,16,17,18,19]. Zhang introduced a multi-slicing strategy for three-dimensional DLO and, with Zhuang and Sun, applied DLO to shallow and shotcrete-supported tunnels, fire-damaged and pre-grouted tunnel stability, and adaptive tailings dams [20,21,22,23,24,25]. Conventional DLO remains an upper-bound method whose unknowns are velocity jumps and plastic multipliers. Related lower-bound developments include Hawksbee’s relaxed pseudo bounds based on translational stress functions [26] and the stress-function interpretation of equilibrium-form DLO by Smith and Gilbert [27]. Neither these lower-bound constructions nor the VDLO projection makes a complete-line strength constraint equivalent to pointwise yield admissibility, which is the central difficulty in a lower-bound DLO.

We develop a statically admissible DLO lower-bound formulation (DLO-LB) for plane Mohr–Coulomb plasticity. The stress discretization follows Sloan [3]; the intersecting line network and companion kinematic problem follow Smith and Gilbert [7]; and line resultants are assembled with the VDLO operator of Zhang et al. [8]. Equilibrium, traction continuity, prescribed stress boundaries, and a conservative polygonal yield condition are imposed on element-local stresses. Complete-line inequalities are retained for stress projection, and a separate calculation is used to demonstrate why they cannot replace the pointwise yield constraints. Dual variables associated with those constraints drive conforming refinement and trace the two Mohr–Coulomb characteristic families. A conventional DLO upper-bound model is refined independently around its active points. The Prandtl punch and frictional-slope examples compare the resulting bounds and critical regions, while a plane-strain extrusion example additionally provides an exact slip-line reference. The verified lower-bound claim is limited to homogeneous, plane, associated perfect plasticity and T3 stress spaces for which equilibrium and yield admissibility hold throughout each element.

Section 2 gives the static formulation, line projection, and characteristic reconstruction. The adaptive discretizations are described in Section 3, followed by the numerical studies in Section 4. Section 5 and Section 6 close the paper.

2. DLO Lower-Bound Formulation

2.1. Problem Statement and Load Partition

Consider a plane-strain body $\Omega \subset {\mathbb{R}}^2$ with boundary $\partial \Omega ={\Gamma }_t\cup {\Gamma }_u$. Body forces and prescribed tractions are partitioned into dead and live contributions,

$$\mathbf{b}\left(\lambda \right)={\mathbf{b}}^d+\lambda {\mathbf{b}}^l,\overline{\mathbf{t}}\left(\lambda \right)={\overline{\mathbf{t}}}^d+\lambda {\overline{\mathbf{t}}}^l,$$

(1)

where $\lambda \ge 0$ is the load multiplier. Self-weight is placed in ${\mathbf{b}}^l$ when a gravity multiplier is sought; surcharge or other actions that remain fixed are assigned to the dead-load terms. The lower-bound problem maximizes $\lambda$ over all statically and plastically admissible stress fields. Tension is positive throughout.

Hereafter, e indexes elements and a an element-local stress degree of freedom. The element-local nodal stresses are written in tensor and Voigt form as

$${\sigma }_a^e=\left[\begin{array}{cc}{\sigma }_{xx,a}^e& {\tau }_{xy,a}^e\\ {\tau }_{xy,a}^e& {\sigma }_{yy,a}^e\end{array}\right],{\mathbf{s}}_a^e={\left[\begin{array}{ccc}{\sigma }_{xx,a}^e& {\sigma }_{yy,a}^e& {\tau }_{xy,a}^e\end{array}\right]}^{\mathsf{T}}.$$

(2)

Elements sharing a geometric node do not share the complete stress tensor. This permits statically admissible stress discontinuities, while the two traction components are matched explicitly on every internal edge.

2.2. Interpolated Stresses and Static Admissibility

Within element e, the tensor interpolation and its matrix counterpart are

$${\sigma }_h^e\left(\mathbf{x}\right)={\displaystyle \sum _{a=1}^{n_e}}{N}_a^e\left(\mathbf{x}\right){\sigma }_a^e,{\mathbf{s}}_h^e\left(\mathbf{x}\right)={\mathbf{N}}_{\sigma }^e\left(\mathbf{x}\right){\mathbf{s}}_e,$$

(3)

where T3 and Q4 denote three-node triangular and four-node quadrilateral stress elements, respectively, $n_e=3$ for T3 and $n_e=4$ for Q4, and

$${\mathbf{s}}_e={\left[\begin{array}{ccc}{\left({\mathbf{s}}_1^e\right)}^{\mathsf{T}}& \dots & {\left({\mathbf{s}}_{n_e}^e\right)}^{\mathsf{T}}\end{array}\right]}^{\mathsf{T}},{\mathbf{N}}_{\sigma }^e=\left[\begin{array}{ccc}{N}_1^e{\mathbf{I}}_3& \dots & N_{n_e}^e{\mathbf{I}}_3\end{array}\right].$$

(4)

where ${\mathbf{I}}_3$ is the $3\times 3$ identity matrix. The linear T3 stress space is denoted P1 below. The strong equilibrium equation is

$$\nabla ·{\sigma }_h^e+{\mathbf{b}}^d+\lambda {\mathbf{b}}^l=\mathbf{0}.$$

(5)

With the adopted Voigt ordering, Eq. (5) becomes

$$\mathbf{D}{\mathbf{N}}_{\sigma }^e{\mathbf{s}}_e+{\mathbf{b}}^d+\lambda {\mathbf{b}}^l=\mathbf{0},\mathbf{D}=\left[\begin{array}{ccc}{\partial }_x& 0& {\partial }_y\\ 0& {\partial }_y& {\partial }_x\end{array}\right].$$

(6)

For a T3 element, the divergence of Eq. (3) is constant, so Eq. (6) is imposed exactly with two rows. The Q4 implementation also supplies eight equilibrium rows evaluated at the four $2\times 2$ Gauss points. This collocation is exact for the affine Q4 stress tests used here, but it does not establish static admissibility throughout a generally distorted isoparametric element.

For unit normal $\mathbf{n}={[n_x,n_y]}^{\mathsf{T}}$, the traction vector has tensor and matrix forms

$${\mathbf{t}}_h\left(\mathbf{n}\right)={\sigma }_h\mathbf{n}=\mathbf{P}\left(\mathbf{n}\right){\mathbf{s}}_h,\mathbf{P}\left(\mathbf{n}\right)=\left[\begin{array}{ccc}{n}_x& 0& n_y\\ 0& n_y& n_x\end{array}\right].$$

(7)

Hence, on an interface ${\Gamma }_{e{e}^{\prime }}$, the endpoint equations

$$\left[\left[{\sigma }_h\mathbf{n}\right]\right]=\left[\left[\mathbf{P}\left(\mathbf{n}\right){\mathbf{s}}_h\right]\right]=\mathbf{0}$$

(8)

where $\left[[·]\right]={(·)}^e-{(·)}^{e^{\prime }}$ denotes the difference between the two element traces. These equations match both traction components. Because edge tractions are linearly interpolated, endpoint equality enforces traction continuity along a complete T3 or affine-Q4 edge. On ${\Gamma }_t$,

$$\mathbf{P}\left(\mathbf{n}\right){\mathbf{s}}_h={\overline{\mathbf{t}}}^d+\lambda {\overline{\mathbf{t}}}^l,$$

(9)

whereas stresses associated with reactions on ${\Gamma }_u$ remain free.

2.3. Pointwise Mohr–Coulomb Admissibility

For cohesion c and friction angle $\varphi$, the plane Mohr–Coulomb condition under the adopted sign convention is

$$\sqrt{{\left({\displaystyle \frac{{\sigma }_{xx}-{\sigma }_{yy}}{2}}\right)}^2+{\tau }_{xy}^2}+{\displaystyle \frac{{\sigma }_{xx}+{\sigma }_{yy}}{2}}sin\varphi \le ccos\varphi .$$

(10)

The circular section in Eq. (10) is replaced by an inward regular polygon with an even number m of sides. Defining ${\vartheta }_k=2\pi (k-1)/m$ and $r_m=cos(\pi /m)$ gives, at every element-local stress node, the matrix inequalities

$${\left({\mathbf{a}}_k^{\mathrm{mc}}\right)}^{\mathsf{T}}{\mathbf{s}}_a^e\le r_{m}ccos\varphi ,k=1,\dots ,m,$$

(11)

where

$${\mathbf{a}}_k^{\mathrm{mc}}={\displaystyle \frac{1}{2}}\left[\begin{array}{c}cos{\vartheta }_k+r_{m}sin\varphi \\ -cos{\vartheta }_k+r_{m}sin\varphi \\ 2sin{\vartheta }_k\end{array}\right].$$

(12)

The factor $r_m$ makes the polygon conservative with respect to the smooth surface. For a T3 element, the interpolated stress is a convex combination of its nodal stresses; nodal satisfaction of the convex polygon therefore implies satisfaction throughout the element.

2.4. Intersecting Candidate Lines and Stress Projection

Let $\mathcal{L}$ be a DLO-style set of candidate lines and ${\Gamma }_{\ell }=\ell \cap \Omega$ the in-domain part of line ℓ. A candidate may intersect other lines and traverse multiple stress elements without being constrained to element edges. Its unit tangent, unit normal, and in-domain length are ${\mathbf{t}}_{\ell }$, ${\mathbf{n}}_{\ell }$, and $L_{\ell }$. The normal and tangential stress components at $\mathbf{x}\in {\Gamma }_{\ell }$ are the tensor projections

$${\sigma }_{n\ell }=({\mathbf{n}}_{\ell }\otimes {\mathbf{n}}_{\ell }):{\sigma }_h,{\tau }_{nt\ell }=({\mathbf{n}}_{\ell }\otimes {\mathbf{t}}_{\ell }):{\sigma }_h.$$

(13)

Because ${\sigma }_h$ is symmetric, $({\mathbf{n}}_{\ell }\otimes {\mathbf{t}}_{\ell }):{\sigma }_h={({\mathbf{n}}_{\ell }\otimes {\mathbf{t}}_{\ell })}^{\mathrm{S}}:{\sigma }_h$. This is the local normal–tangential traction decomposition used in cracking-elements formulations [28,29,30]. In Voigt form,

$$\left[\begin{array}{c}{\sigma }_{n\ell }\\ {\tau }_{nt\ell }\end{array}\right]={\mathbf{Q}}_{\ell }{\mathbf{s}}_h,{\mathbf{Q}}_{\ell }=\left[\begin{array}{ccc}{n}_x^2& n_y^2& 2n_x{n}_y\\ n_x{t}_x& n_y{t}_y& n_x{t}_y+n_y{t}_x\end{array}\right],$$

(14)

where the subscript ℓ on the components of ${\mathbf{n}}_{\ell }$ and ${\mathbf{t}}_{\ell }$ is omitted inside ${\mathbf{Q}}_{\ell }$. Let ${\Omega }_e$ denote the domain of element e and let $\mathbf{s}$ collect all element vectors ${\mathbf{s}}_e$. Integrating Eq. (14) along the complete in-domain line gives

$${\mathbf{r}}_{\ell }=\left[\begin{array}{c}{N}_{\ell }\\ T_{\ell }\end{array}\right]=\sum _{e:{\Gamma }_{\ell }\cap {\Omega }_e\ne ⌀}\underset{{\displaystyle {\mathbf{H}}_{\ell e}}}{\underbrace{{\int }_{{\Gamma }_{\ell }\cap {\Omega }_e}{\mathbf{Q}}_{\ell }{\mathbf{N}}_{\sigma }^e\left(\mathbf{x}\right)\mathrm{d}s}}{\mathbf{s}}_e={\mathbf{H}}_{\ell }\mathbf{s}.$$

(15)

Thus, $N_{\ell }$ and $T_{\ell }$ are the line-integrated normal and tangential stress resultants, respectively. Equation (11) of VDLO extends the local traction relation in Eq. (13) to a line-integrated virtual-work row [8]. In the present formulation, the prescribed instantaneous field is replaced by the unknown stress vector $\mathbf{s}$. Clipping each line against the elements that it crosses yields the sparse matrix ${\mathbf{H}}_{\ell }$.

For homogeneous strength, two linear complete-line inequalities can be formed,

$$T_{\ell }+tan\varphi N_{\ell }\le c{L}_{\ell },-T_{\ell }+tan\varphi N_{\ell }\le c{L}_{\ell }.$$

(16)

These rows quantify line-level mobilization but do not replace Eq. (11). Integration can compensate a locally over-yielded segment with less mobilized portions of the same line. Conversely, once pointwise admissibility is imposed for a homogeneous convex strength domain, Eq. (16) is an integral consequence and is normally redundant. The candidate lines can thus measure and display mobilization, but the pointwise constraints establish the lower bound.

2.5. Sparse Linear Program

Collecting all element-local stress vectors in $\mathbf{s}$ and appending the multiplier gives $\mathbf{x}={[{\mathbf{s}}^{\mathsf{T}},\lambda ]}^{\mathsf{T}}$. The global matrix problem is

$$\begin{array}{cccc}\hfill & \underset{\mathbf{x}}{maximize}\hfill & \hfill & \lambda \hfill \\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill & {\mathbf{A}}_{\mathrm{eq}}\mathbf{x}={\mathbf{b}}_{\mathrm{eq}},\hfill \\ \hfill & \hfill & & {\mathbf{A}}_{\mathrm{mc}}\mathbf{x}\le {\mathbf{b}}_{\mathrm{mc}},\hfill \\ \hfill & \hfill & & {\mathbf{A}}_{\mathrm{line}}\mathbf{x}\le {\mathbf{b}}_{\mathrm{line}},\hfill \\ \hfill & \hfill & & \lambda \ge 0.\hfill \end{array}$$

(17)

Here ${\mathbf{A}}_{\mathrm{eq}}$ contains equilibrium, traction continuity, and prescribed stress boundaries; ${\mathbf{A}}_{\mathrm{mc}}$ contains the pointwise inward Mohr–Coulomb inequalities; and ${\mathbf{A}}_{\mathrm{line}}$ contains optional complete-line rows. The sparse problem is assembled once and passed to a linear-programming solver. No active-line preselection, nonlinear return mapping, or load stepping is used.

Figure 1 collects these operations on a schematic half-symmetric Prandtl domain. The Sloan-type mesh in panel (a) supplies the element-local stress variables and the equilibrium, traction-continuity, and pointwise-yield rows [3]. Independently, the DLO ground structure in panel (b) joins admissible node pairs by lines that may cross without creating intersection nodes [7]. Panel (c) combines the interpolated stress field and a candidate line through the VDLO projection [8], producing normal and tangential resultants and the line-level Mohr–Coulomb rows assembled in ${\mathbf{A}}_{\mathrm{line}}$. Panel (d) illustrates a critical line pattern formed from members already present in the candidate set. This panel is schematic, not a computed mechanism. Yield admissibility defines the bound; static characteristics are not velocity discontinuities.

2.6. Selection and Interpretation of a Non-Unique Optimal Stress Field

The primary linear program may have several admissible stress fields at the same optimal multiplier. Two secondary linear programs select one of them for reporting. RGB mapping and Mohr–Coulomb characteristic tracing are then applied to that selected field. The two selection problems retain the primary optimum; the subsequent plots involve no further optimization.

2.6.1. Lexicographic Selection of a Representative Stress Field

The primary linear program can admit multiple stress fields at the same optimal load. This non-uniqueness does not affect the lower-bound value, but an arbitrary optimizer vertex can contain pronounced element-scale variations in stress components that are not fixed by interface equilibrium. In particular, traction continuity gives $\left[\left[{\sigma }_h\right]\right]\mathbf{n}=\mathbf{0}$ on an internal edge, whereas the tangential–tangential component ${\mathbf{t}}^{\mathsf{T}}\left[\left[{\sigma }_h\right]\right]\mathbf{t}$ may remain discontinuous. A secondary optimization is therefore used to select, rather than smooth, one member of the primary optimal set.

Let $q^{*}$ denote the optimum of the primary objective and $\mathcal{F}\left(q^{*}\right)$ the original feasible set with that objective fixed at $q^{*}$. The discrete interelement-jump measure is

$$J_h\left(\mathbf{s}\right)=\sum _{\gamma \in {\mathcal{E}}_h^{\mathrm{int}}}{L}_{\gamma }{\displaystyle \sum _{r=1}^{m_{\gamma }}}{\omega }_{\gamma r}\left|{\mathbf{t}}_{\gamma }^{\mathsf{T}}\left[\left[{\sigma }_h\left({\mathbf{x}}_{\gamma r}\right)\right]\right]{\mathbf{t}}_{\gamma }\right|,$$

(18)

where ${\mathcal{E}}_h^{\mathrm{int}}$ is the set of internal edges; $L_{\gamma }$ and ${\mathbf{t}}_{\gamma }$ are the length and unit tangent of edge $\gamma$; and ${\{{\mathbf{x}}_{\gamma r},{\omega }_{\gamma r}\}}_{r=1}^{m_{\gamma }}$ are its sampling points and normalized quadrature weights. A P1 edge uses its two endpoints with ${\mathbf{\omega }}_{\gamma }=(1/2,1/2)$, whereas a P2 edge uses its two endpoints and geometric midpoint with ${\mathbf{\omega }}_{\gamma }=(1/6,4/6,1/6)$. At the midpoint, the quadratic Bernstein stress is evaluated from its edge controls with weights $(1/4,1/2,1/4)$; the middle control is not treated as a nodal stress. Introducing nonnegative epigraph variables for the absolute values gives the first selection problem

$$J_h^{*}=\underset{\mathbf{s}\in \mathcal{F}\left(q^{*}\right)}{min}{J}_h\left(\mathbf{s}\right).$$

(19)

All equilibrium, boundary-traction, and pointwise-yield rows of the primary problem are retained. Equation (19) consequently changes neither admissibility nor $q^{*}$.

Minimum jump alone can still leave an unnecessarily confined member of the optimal set. With tension positive, define the compressive mean-stress epigraph at local control a of element e by

$$p_a^e\ge -{\displaystyle \frac{{\sigma }_{xx,a}^e+{\sigma }_{yy,a}^e}{2}},p_a^e\ge 0.$$

(20)

The associated functional is $P_h\left(\mathbf{s}\right)={min}_{\mathbf{p}}{\sum }_e{\sum }_{a=1}^{n_e}{w}_{ea}{p}_a^e$ subject to Eq. (20), where $\mathbf{p}=\left\{p_a^e\right\}$, $w_{ea}=A_e/n_e$, $A_e$ is the element area, and $n_e$ is its number of stress controls (three for P1 T3 and six for P2 T3). A lexicographic second stage then solves

$$\underset{\mathbf{s}\in \mathcal{F}\left(q^{*}\right)}{min}{P}_h\left(\mathbf{s}\right)\mathrm{subject}\mathrm{to}{J}_h\left(\mathbf{s}\right)\le J_h^{*}+{\epsilon }_J.$$

(21)

The small tolerance ${\epsilon }_J$ prevents numerical loss of the first-stage optimum. Both stages are linear programs. Together they define a lexicographic gauge for the non-unique stress field; they are neither a constitutive update nor cross-element stress averaging.

The implementation uses three complete linear-program solves. The first solve stores the primary objective value $q^{*}$ and load multiplier ${\lambda }^{*}$. For the second solve, all primary variables, rows, and bounds are copied, and $\lambda$ is fixed to ${\lambda }^{*}$. When the primary objective is not $\lambda$ itself, the additional equality ${\mathbf{c}}^{\mathsf{T}}\mathbf{x}=q^{*}$ is imposed, where $\mathbf{c}$ is the coefficient vector of the primary linear objective. For every jump sample $j_{\gamma r}={\mathbf{t}}_{\gamma }^{\mathsf{T}}\left[\left[{\sigma }_h\left({\mathbf{x}}_{\gamma r}\right)\right]\right]{\mathbf{t}}_{\gamma }$, a nonnegative epigraph $z_{\gamma r}$ and the two rows $j_{\gamma r}\le z_{\gamma r}$ and $-j_{\gamma r}\le z_{\gamma r}$ are added before minimizing ${\sum }_{\gamma ,r}{L}_{\gamma }{\omega }_{\gamma r}{z}_{\gamma r}$. The third solve retains the same fixed primary optimum and appends the budget

$$\sum _{\gamma ,r}{L}_{\gamma }{\omega }_{\gamma r}{z}_{\gamma r}\le J_h^{*}+{10}^{-8}max(1,|J_h^{*}\left|\right)$$

(22)

before the pressure epigraphs in Eq. (20) are minimized. The reported calculations use the area weights $w_{ea}=A_e/n_e$ without a support-distance bias. A secondary solution is accepted only after the code has checked preservation of ${\lambda }^{*}$ to a relative scale of ${10}^{-8}$ and the jump budget to a relative scale of ${10}^{-6}$. These checks keep the reported stress field on the original optimal face.

2.6.2. RGB Representation of the Stress Tensor Field

For compact visualization, the three normalized components $\widehat{\mathbf{s}}={[{\sigma }_{xx},{\sigma }_{yy},{\tau }_{xy}]}^{\mathsf{T}}/c$ are assigned to the red, green, and blue channels,

$$C_j=clip\left({\displaystyle \frac{{\widehat{s}}_j-a_j}{b_j-a_j}},0,1\right),j=1,2,3.$$

(23)

Here $clip(v,0,1)=min\{1,max\{0,v\left\}\right\}$. The component limits $[a_j,b_j]$ are held fixed across every panel in a stated comparison. Colors are evaluated directly from the element-local stress interpolation; no global nodal averaging or image-space smoothing is applied.

Figure 2. Post-solution selection and RGB representation of one admissible stress field from the primary optimal set.

Figure 2. Post-solution selection and RGB representation of one admissible stress field from the primary optimal set.

Section 4 reports the corresponding stress maps, computational costs, and checks that the lower-bound value is unchanged.

2.6.3. Static Characteristic-Field Recovery

The optimization problem in Eq. (17) contains no displacement jump, and its active yield points should not be interpreted as a kinematic failure mechanism. For mechanical interpretation, a continuous static field is recovered from the two Mohr–Coulomb characteristic directions of the selected stress solution. Let the major-principal direction be

$$\psi ={\displaystyle \frac{1}{2}}atan2\left(2{\tau }_{xy},{\sigma }_{xx}-{\sigma }_{yy}\right).$$

(24)

The two unoriented characteristic tangents are

$${\theta }^{\alpha }=\psi -\left({\displaystyle \frac{\pi }{4}}+{\displaystyle \frac{\varphi }{2}}\right),{\theta }^{\beta }=\psi +\left({\displaystyle \frac{\pi }{4}}+{\displaystyle \frac{\varphi }{2}}\right).$$

(25)

The superscripts $\alpha$ and $\beta$ label the two conjugate, unoriented characteristic families; they denote neither bound types nor additional optimization variables. For $\varphi =0$, they reduce to the orthogonal maximum-shear directions $\psi \pm \pi /4$. For $\varphi >0$, their acute intersection angle is $\pi /2-\varphi$. Direct interpolation of $\theta$ is unreliable near the $0/\pi$ branch. For each family $q\in \{\alpha ,\beta \}$, the implementation interpolates the double-angle orientation vector

$${\mathbf{d}}_h^q\left(\mathbf{x}\right)=\sum _a{N}_a\left(\mathbf{x}\right){\{cos2{\theta }_a^q,sin2{\theta }_a^q\}}^{\mathsf{T}},{\theta }_h^q={\displaystyle \frac{1}{2}}atan2(d_{h,2}^q,d_{h,1}^q).$$

(26)

Here $N_a$ is the current element’s stress shape function and ${\theta }_a^q$ is the direction of family q evaluated from local stress control a. Characteristic curves then satisfy

$${\displaystyle \frac{\mathrm{d}\mathbf{x}}{\mathrm{d}s}}=\pm {\{cos{\theta }_h^q,sin{\theta }_h^q\}}^{\mathsf{T}}.$$

(27)

The sign is chosen after every element crossing to maximize alignment with the preceding direction. Seeds and stopping points are selected from the region identified by active dual variables of the pointwise yield constraints. The reconstruction uses only the selected stress field and the dual variables of its yield constraints. It adds no optimization variable or compatible velocity field, and it leaves both the stress solution and collapse multiplier unchanged.

3. Adaptive Lower-Bound Stress Discretization

The two bounds require different adaptive operations. On the lower-bound side, the admissible stress space is enriched by local mesh refinement (h-adaptivity) or by increasing the stress order (p-adaptivity). The conventional DLO upper bound instead inserts points near active velocity discontinuities and regenerates its line network [25]. It is refined independently for the comparison in Section 4. All three procedures retain a sparse linear-programming problem.

3.1. h-Adaptivity by Dual-Weighted Mesh Refinement

The accuracy of a lower-bound stress approximation depends on where stress discontinuities and additional stress degrees of freedom are introduced [6]. Here the dual multipliers of the pointwise Mohr–Coulomb inequalities are used to identify the elements that control the current optimum. For element e in the current mesh ${\mathcal{T}}_h$, the indicator is

$${\eta }_e={\displaystyle \sum _{a=1}^{n_e}}{\displaystyle \sum _{k=1}^m}\left|y_{eak}\right|,$$

(28)

where $y_{eak}$ is the optimizer dual associated with polygon side k at local stress node a. The resulting quantity measures sensitivity to the strength constraints; it is not a plastic strain or a slip-line density.

A Dörfler bulk criterion marks the smallest set $\mathcal{M}$ obtained by descending ${\eta }_e$ such that

$$\sum _{e\in \mathcal{M}}{\eta }_e\ge \vartheta \sum _{e\in {\mathcal{T}}_h}{\eta }_e,0<\vartheta \le 1.$$

(29)

For a full-domain problem with mirror-symmetric geometry, loading, and boundary conditions, the marking units are mirror orbits ${\mathcal{O}}_e=\{e,e^{★}\}$. Their indicators are ${\eta }_{{\mathcal{O}}_e}={\eta }_e+{\eta }_{e^{★}}$, and selecting an orbit marks both elements. Pairing prevents optimizer tie-breaking or a non-unique dual solution from biasing refinement toward one side of a symmetric model. Mirror partners are matched geometrically using the tolerance ${10}^{-10}max(1,L_x,L_y)$, where $L_x$ and $L_y$ are the domain extents. They are required to form a unique involution before marking proceeds.

Marked T3 elements are red-refined. Neighboring elements receive conforming green transitions; if a green template would create a child below the prescribed quality threshold, that parent is promoted to red refinement and edge closure is repeated. In the mirror-paired sequence, a closure element with two split edges is also promoted to red refinement. The promotion removes any choice of a three-child diagonal based on local node numbering and preserves mirror compatibility after every refinement. Exposed boundary tags and linearly varying boundary tractions are inherited by child edges. The child-quality measure used in this test is

$$q_K={\displaystyle \frac{4\sqrt{3}{A}_K}{{\ell }_{K1}^2+{\ell }_{K2}^2+{\ell }_{K3}^2}},0<q_K\le 1,$$

(30)

where $A_K$ is the area of triangle K and ${\ell }_{Ki}$ are its three side lengths. The reported calculations set the minimum predicted child quality to $q_{min}=0.30$.

In one adaptive cycle, the current sparse linear program is assembled and solved; the yield-row duals are accumulated into ${\eta }_e$; elements or mirror orbits are sorted in descending indicator order and marked by Eq. (29); all three edges of every marked parent are split; and conformity and quality closure are repeated until no new edge is introduced. One midpoint is then created per split edge, child connectivity and boundary records are generated, and a new sparse problem is assembled on the refined mesh. At least one marking unit is retained if the dual indicator is degenerate, so a requested refinement cycle cannot stall. The studies in Section 4 use $\vartheta =0.60$ and a prescribed number of levels rather than an estimator-tolerance stopping test; the last reported level is solved but not marked.

The previous feasible stress field is prolonged to the refined mesh. A parent’s linear stress is evaluated at every new child node, body loads are copied, and interface tractions are inherited. The prolongated field therefore satisfies child equilibrium, traction continuity, boundary conditions, and the convex pointwise yield condition whenever the parent field does. Hence, the feasible set at level h is embedded in the feasible set at level $h/2$, and

$${\lambda }_{h/2}^{\mathrm{LB}}\ge {\lambda }_h^{\mathrm{LB}}$$

(31)

up to optimization tolerances. This monotonicity is a property of the nested adaptive sequence; independently regenerated mesh families need not be nested and should be reported as a mesh-family envelope. This prolongation proves nestedness and supplies an explicit feasible child field. The timings reported below still include a full rebuild and solve at every level; no optimizer warm start is used.

3.2. p-Adaptivity by Bernstein Stress Enrichment

The second route raises the polynomial order of the stress approximation while retaining the affine T3 geometry. To distinguish the element coordinates from the load multiplier $\lambda$, let $\{{\zeta }_1,{\zeta }_2,{\zeta }_3\}$ denote the barycentric coordinates. The quadratic Bernstein basis is

$$\begin{array}{cc}\hfill {\mathbf{B}}^{\left(2\right)}=& \left[{\zeta }_1^2,{\zeta }_2^2,{\zeta }_3^2,2{\zeta }_1{\zeta }_2,2{\zeta }_2{\zeta }_3,2{\zeta }_3{\zeta }_1\right],\hfill \\ \hfill {\sigma }_h^e=& {\displaystyle \sum _{A=1}^6}{B}_A^{\left(2\right)}{\widehat{\sigma }}_A^e,{\mathbf{s}}_h^e={\mathbf{N}}_{\sigma }^{e,\left(2\right)}{\widehat{\mathbf{s}}}_e,\hfill \\ \hfill {\mathbf{N}}_{\sigma }^{e,\left(2\right)}=& \left[\begin{array}{cccc}{B}_1^{\left(2\right)}{\mathbf{I}}_3& B_2^{\left(2\right)}{\mathbf{I}}_3& \dots & B_6^{\left(2\right)}{\mathbf{I}}_3\end{array}\right].\hfill \end{array}$$

(32)

Here $B_A^{\left(2\right)}$ is the Ath scalar component of ${\mathbf{B}}^{\left(2\right)}$; its $3\times 3$ Voigt interpolation block is $B_A^{\left(2\right)}{\mathbf{I}}_3$. The stacked control vector is ${\widehat{\mathbf{s}}}_e={[{\left({\widehat{\mathbf{s}}}_1^e\right)}^{\mathsf{T}},\dots ,{\left({\widehat{\mathbf{s}}}_6^e\right)}^{\mathsf{T}}]}^{\mathsf{T}}\in {\mathbb{R}}^{18}$, with ${\widehat{\mathbf{s}}}_A^e$ the Voigt vector of ${\widehat{\sigma }}_A^e$. The six element-local stress controls are located at the three vertices and three edge midpoints. Because $B_A^{\left(2\right)}\ge 0$ and ${\sum }_A{B}_A^{\left(2\right)}=1$, imposing the convex inward yield polygon on every ${\widehat{\sigma }}_A^e$ guarantees pointwise yield admissibility throughout the element. The term edge-midpoint control describes a Bernstein control quantity, not the stress sampled at the geometric midpoint. For an edge with endpoint controls ${\widehat{\mathbf{s}}}_i,{\widehat{\mathbf{s}}}_j$ and edge control ${\widehat{\mathbf{s}}}_{ij}$, the midpoint value is ${\mathbf{s}}_h\left({\mathbf{x}}_{\mathrm{mid}}\right)=\frac{1}{4}{\widehat{\mathbf{s}}}_i+\frac{1}{2}{\widehat{\mathbf{s}}}_{ij}+\frac{1}{4}{\widehat{\mathbf{s}}}_j$.

The divergence of the quadratic stress field is linear. Exact element equilibrium is therefore obtained by enforcing

$$\mathbf{D}{\mathbf{N}}_{\sigma }^{e,\left(2\right)}\left({\mathbf{x}}_i\right){\widehat{\mathbf{s}}}_e+{\mathbf{b}}^d+\lambda {\mathbf{b}}^l=\mathbf{0},i=1,2,3,$$

(33)

at the three triangle vertices, giving six scalar rows per element. Quadratic edge traction has two endpoint controls and one midpoint control, so continuity requires

$$\mathbf{P}\left(\mathbf{n}\right)\left({\widehat{\mathbf{s}}}_{A_r}^e-{\widehat{\mathbf{s}}}_{A_r^{\prime }}^{e^{\prime }}\right)=\mathbf{0},r=1,2,3,$$

(34)

where $A_r$ and $A_r^{\prime }$ are the local control indices associated with the same physical endpoint or midpoint on the two sides of interface ${\Gamma }_{e{e}^{\prime }}$. The same three-control construction is used for prescribed boundary tractions. Candidate-line resultants are evaluated with three-point Gauss integration, which is exact for the quadratic variation along each straight line segment. Every coefficient remains linear in the stress controls and the load multiplier.

The row construction is explicit. Each P2 T3 contributes 18 stress unknowns, six scalar equilibrium equalities, and $6m$ pointwise polygon inequalities. Each internal edge contributes six traction-continuity equalities, namely two traction components at each of its three Bernstein edge controls. A prescribed boundary edge contributes three rows for every prescribed traction component; linearly varying prescribed tractions use their endpoint values and their average at the middle control.

Candidate lines are generated from unordered geometric-node pairs. A pair is discarded when another node lies on its open segment, using the tolerance ${10}^{-10}max(1,L_{\ell })$; this removes collinear redundancies without preventing different candidates from crossing. Every retained line is clipped against each convex element. The resulting in-element intervals are sorted, overlap is removed, and intervals shorter than ${10}^{-12}$ in normalized line coordinate are ignored. On every remaining interval, the three Gauss–Legendre abscissae $\{-\sqrt{3/5},0,\sqrt{3/5}\}$ with weights $\{5/9,8/9,5/9\}$ accumulate the coefficients of ${\mathbf{H}}_{\ell e}$ in Eq. (15). These operations are performed during sparse assembly; no stress field is sampled or transferred after the solve.

The P1 space is embedded in the P2 space by degree elevation. If ${\mathbf{s}}_1$, ${\mathbf{s}}_2$, and ${\mathbf{s}}_3$ are the P1 vertex stresses, with the element superscript omitted, the equivalent P2 controls are

$${\widehat{\mathbf{s}}}_a={\mathbf{s}}_a(a=1,2,3),{\widehat{\mathbf{s}}}_4={\displaystyle \frac{{\mathbf{s}}_1+{\mathbf{s}}_2}{2}},{\widehat{\mathbf{s}}}_5={\displaystyle \frac{{\mathbf{s}}_2+{\mathbf{s}}_3}{2}},{\widehat{\mathbf{s}}}_6={\displaystyle \frac{{\mathbf{s}}_3+{\mathbf{s}}_1}{2}}.$$

(35)

Degree elevation therefore gives, on the same mesh,

$${\lambda }_{h,2}^{\mathrm{LB}}\ge {\lambda }_{h,1}^{\mathrm{LB}}$$

(36)

up to optimization tolerances, where the second subscript denotes the stress order. The implementation uses global P1-to-P2 enrichment. Element-wise mixed-order marking is not used in the reported calculations; this restriction keeps P1–P2 interface treatment outside the current validation scope. For the combined $h+p$ sequences reported later, the order flag is set to P2 on every element of every adaptive mesh, and the P2 pointwise-yield duals are used in Eq. (28). The same marking, closure, boundary inheritance, and prescribed-level stopping rules as in Section 3.1 are then applied. Here $h+p$ denotes global quadratic stress enrichment followed by adaptive h-refinement, not local switching between P1 and P2 elements.

4. Numerical Verification

All sparse linear programs were solved with MOSEK. Unless stated otherwise, the lower-bound calculations used T3 stress elements and an inward 32-sided Mohr–Coulomb polygon. The companion upper bounds used conventional DLO jump variables and nonnegative plastic multipliers on all admissible point-to-point lines after removal of collinear redundancies. Lower-bound characteristics and upper-bound active lines are different quantities. The former are traced from a statically admissible stress field; the latter form a compatible kinematic mechanism. For all characteristic-field plots, seeds were drawn from elements whose maximum yield-row dual exceeded $1\%$ of the case maximum and were separated using the area-derived mesh scale $h=\left(2\right|\Omega |/N_{\mathrm{el}}{)}^{1/2}$, where $N_{\mathrm{el}}$ is the element count. The integration step and seed spacing were $0.18h$ and $1.6h$, respectively, and tracing stopped at the boundary or when the interpolated dual activity fell below the same display threshold. No geometric smoothing was applied. These settings affect only the plotted characteristic curves, not the optimization, adaptivity, or collapse multiplier. Optimizer times reported below exclude assembly and output, follow one unreported warm-up solve, and were measured on an Apple M5 processor with 10 cores and 24 GB of memory using MOSEK 11.2. They are single-run wall-clock measurements intended to compare formulations within this study, rather than hardware-independent performance benchmarks. The adaptive tables report the primary bound-optimization time only; auxiliary lexicographic stress-selection solves are excluded. When an exact value $q_{\mathrm{ex}}$ is available, the reported bracket is $100(q^{\mathrm{UB}}-q^{\mathrm{LB}})/q_{\mathrm{ex}}$; otherwise the gap is $100(q^{\mathrm{UB}}-q^{\mathrm{LB}})/q^{\mathrm{UB}}$.

4.1. Complete Symmetric Prandtl Punch

The first benchmark is a weightless Tresca material with $c=1$ beneath a centered rigid rough strip footing. The complete domain is $10b$ wide and $4b$ deep, and the footing occupies $4b\le x\le 6b$, giving a total footing width $2b$. The vertical sides and base are reaction boundaries in the lower-bound problem and fixed boundaries in the upper-bound problem. The footing segments share one vertical velocity and have zero horizontal velocity. The exact bearing-capacity factor is $N_c=2+\pi =5.14159265$. The lower-bound objective is the mean compressive reaction beneath the footing; the upper-bound problem minimizes dissipation under unit footing work.

Figure 3 shows the physical problem and the coarsest discretizations used in the convergence study. The lower-bound model contains 80 T3 stress elements on a checkerboard grid. The upper-bound model uses the same 55 geometric points as a point cloud; admissible point-to-point lines are generated from this cloud but are omitted from the model panel for clarity.

Table 1 reports six complete-domain discretizations, all using the 32-sided inward polygon specified above. The notation in the first column gives the matched half-domain resolution; the complete model contains twice as many cells in the horizontal direction. The lower bound increases monotonically, and the normalized bracket contracts at every level from $21.06\%$ to $4.14\%$. The upper-bound sequence approaches the exact solution from above, although the final value is slightly higher than the preceding one because independently generated DLO point sets are not nested. Complete- and half-domain factors agree to the reported precision for the even half grids. On odd half grids, checkerboard T3 diagonals are not mirrored; this topological difference does not invalidate the lower bound.

Figure 4 plots the same sequence against the number of complete-domain nodes. Every pair encloses the exact solution, and the distance between the curves decreases with refinement.

A separate line-only comparison shows why pointwise yield admissibility is necessary. On the coarse $5\times 4$ half grid, complete-line Tresca inequalities return $N_c=5.666667$, but the maximum local Tresca ratio is $9.648$. Adding the inward pointwise polygon gives $N_c=4.031043$ and a maximum ratio of one. Adding the line rows to this statically and plastically admissible model leaves the optimum unchanged. Line integration retains the DLO representation but cannot replace continuum admissibility.

Figure 5 compares the coarsest $10\times 4$ and a representative $40\times 16$ complete-domain solution. Because the symmetric linear program admits non-unique optimal stress fields, the reported static field is the convex average of the computed solution and its mirror image; this operation preserves equilibrium, convex yield admissibility, and the collapse multiplier. The lower-bound panels identify symmetric footing-edge critical regions and two traced characteristic families. The upper-bound panels independently recover a mirror-symmetric pair of collapse fans. Enforcing the rough-footing horizontal velocity condition is essential to this symmetry; locking only the vertical velocity permits an energetically equivalent one-sided optimum in the degenerate linear program.

Mesh-shape tests were also conducted without introducing another benchmark. Across four regular T3 diagonal patterns at a half-grid resolution of $40\times 32$, the lower bounds ranged from $4.629681$ to $4.997661$. Ten-percent interior-node perturbations reduced the sample standard deviation from $0.1701$ at $10\times 8$ to $0.0908$ at $20\times 16$. Refinement reduced, but did not eliminate, topology sensitivity. Affine Q4 tests were more restrictive, and the distorted-Q4 collocation did not preserve static admissibility throughout the element under perturbation. The Q4 results are accordingly treated as a supplementary assessment, not as lower-bound evidence.

The stress-field selection in Section 2.6.1 was first examined on the $10\times 8$ symmetric half-domain T3 discretization. For visual comparison in Figure 6, the half-domain solutions were reflected about the footing centerline; the reflected ${\tau }_{xy}$ component changes sign. The primary and lexicographically selected fields both gave $N_c=4.758506$ and a maximum Tresca ratio of one. The weighted tangential-jump measure decreased from $68.0098$ to $6.00367$, while the area-averaged compressive pressure decreased from $1.35192c$ to $1.22220c$. The maximum individual jump remained $3.14042c$, showing that the selection reduces the domain-wide jump measure without imposing artificial pointwise stress continuity.

The complete Prandtl model was then used to compare P1 h-adaptivity with a combined $h+p$ sequence. Both started from the same $10\times 4$ checker-A grid and used $\vartheta =0.6$, a minimum child quality of $0.30$, and the same 32-sided yield polygon. Mirror-pair indicators and symmetry-preserving closure were used in both sequences. In the first sequence, P1 yield duals drove mesh refinement. In the second, a global P2 Bernstein stress interpolation was used on every mesh and its own yield duals drove refinement. Thus, “$h+p$” denotes P2 enrichment combined with h-adaptivity, not element-wise mixing of P1 and P2 stresses. Every reported adaptive mesh remained exactly compatible with reflection about the footing centerline.

Table 2. Prandtl P1 h-adaptive and P2 $h+p$-adaptive sequences.

	P1 h-adaptive			P2 $\mathit{h}+\mathit{p}$-adaptive
Level	T3	$N_c^{\mathrm{LB}}$	Time (s)	T3	$N_c^{\mathrm{LB}}$	Time (s)
0	80	4.583742	0.018	80	4.726271	0.054
1	100	4.680908	0.021	100	4.842493	0.086
2	148	4.749794	0.047	174	4.951714	0.179
3	352	4.896535	0.137	406	5.026566	0.518
4	898	4.978773	0.517	1012	5.062206	30.204

Table 2. Prandtl P1 h-adaptive and P2 $h+p$-adaptive sequences.

	P1 h-adaptive			P2 $\mathit{h}+\mathit{p}$-adaptive
Level	T3	$N_c^{\mathrm{LB}}$	Time (s)	T3	$N_c^{\mathrm{LB}}$	Time (s)
0	80	4.583742	0.018	80	4.726271	0.054
1	100	4.680908	0.021	100	4.842493	0.086
2	148	4.749794	0.047	174	4.951714	0.179
3	352	4.896535	0.137	406	5.026566	0.518
4	898	4.978773	0.517	1012	5.062206	30.204

In Figure 7, both indicators concentrate refinement beneath the two footing edges and preserve the two reflected fans. The P2 sequence develops a broader refined region. Orange elements mark the next refinement; final panels are unmarked because each sequence stops there.

Both nested sequences were monotone. As shown in Figure 8, the P1 error decreased from $-10.85\%$ to $-3.17\%$, whereas the P2 error decreased from $-8.08\%$ to $-1.54\%$. At the final level, the combined P2 sequence raised the lower bound by $1.68\%$ relative to P1. It also increased the element count from 898 to 1012 and the measured optimizer time from $0.517$ to $30.204$ s. The enriched stress space reduced the remaining discretization error, but the largest P2 linear program incurred a substantial cost. The mirror constraint was introduced to preserve the benchmark symmetry, not as an efficiency or same-size accuracy improvement.

4.2. One Frictional Slope over Three Friction Angles

The second benchmark follows the slope studied by Wang and Zhang [31]. Its polygonal boundary has vertices $(0,0)$, $(20,0)$, $(20,13)$, $(12,13)$, $(2,3)$, and $(0,3)$ m. Cohesion and unit weight are fixed at $c=12.38$ kPa and $\gamma =20.0124$ kN/m${}^3$, while $\varphi ={10}^{\circ }$, ${20}^{\circ }$, or ${30}^{\circ }$. The base is fixed, the vertical sides are normal-fixed and tangentially smooth, and the ground surface is traction free. Gravity is the only live action, so ${\lambda }_g$ multiplies the gravitational acceleration.

All three material cases were solved on the same $h=0.5$ m T3 geometry with 817 points and 1520 stress elements. The upper-bound network contained 201,153 admissible lines. Both bounds increase monotonically with friction angle in Table 3 and Figure 9. The bracket remains below $8\%$ at all three values, although it widened modestly as $\varphi$ increased.

Figure 10 compares the recovered lower-bound characteristics with the upper-bound DLO mechanisms. Grey points in the top row mark pointwise constraints whose dual exceeds $1\%$ of the case maximum; blue and orange curves are the two continuously traced characteristic families. The dominant lower-bound family follows the same toe-to-crest region as the principal upper-bound chain, while the conjugate family crosses that region. As $\varphi$ increases, both constructions produce a narrower and steeper mobilized zone. The agreement is spatial rather than kinematic: no velocity or displacement is reconstructed from the lower-bound stress field.

At $\varphi ={20}^{\circ }$, the 817-point upper bound is $1.043328$, differing by $0.062\%$ from the published value $1.042682$ [31]. The same case was used to test paired adaptivity. The initial checker-A discretization had $h=1$ m. The lower-bound calculation used $m=32$, a bulk fraction $\vartheta =0.6$, and a minimum T3 child quality of $0.30$. The upper-bound cloud accepted at most 24 new points per level, using eight trial directions at $0.35$ times the local nearest-neighbor distance.

The adaptive sequence in Table 4 is monotone from both sides. The lower bound increased from $0.913233$ to $0.974991$, the upper bound decreased from $1.058189$ to $1.052160$, and the bracket contracted from $13.70\%$ to $7.33\%$. The independently refined 817-point upper bound is lower because its globally uniform point cloud is substantially larger than the final 315-point adaptive cloud; it is reported as a separate reference, not as a member of the nested adaptive sequence. Figure 12 displays levels 0, 2, and 4 of the paired sequence. The upper row overlays the blue and orange lower-bound characteristic families on the adapted T3 meshes. In the lower row, the upper-bound point cloud is overlaid by the active DLO slip lines, whose colour and width increase with the normalized plastic multiplier. Thus, the panels show how each independently adapted discretization and its associated critical pattern evolve together. In this and the subsequent adaptive table, $N_L$ is the number of lower-bound mesh nodes and $N_U$ is the number of upper-bound point-cloud points.

Figure 10. Lower-bound characteristic lines (top) and upper-bound active DLO lines (bottom) for three friction angles.

Figure 10. Lower-bound characteristic lines (top) and upper-bound active DLO lines (bottom) for three friction angles.

Figure 11. Opposite monotone adaptive bound sequences for the slope.

Figure 11. Opposite monotone adaptive bound sequences for the slope.

Figure 12. Evolution of the lower-bound T3 mesh with both characteristic families and the upper-bound DLO point cloud with active slip lines.

Figure 12. Evolution of the lower-bound T3 mesh with both characteristic families and the upper-bound DLO point cloud with active slip lines.

For the lower-bound comparison, the P1 sequence was retained through level 4, whereas the combined P2 sequence was evaluated through level 3. Each sequence used its own yield-dual indicator and therefore generated a different nested mesh family. Table 5 reports the measured optimizer cost together with the bound evolution.

At nearly equal element counts on levels 2 and 3, P2 raised the lower bound by $0.01534$ and $0.01409$, respectively, but required substantially longer optimizer times. The comparison changes when a target bound is considered: the 635-element P2 solution reached $0.976682$ in $40.474$ s, exceeding the 2001-element P1 value of $0.974991$ obtained in $266.639$ s. P2 costs more per element but is cheaper at this target accuracy.

Figure 14 pairs the P1 h-adaptive meshes with their lexicographically selected RGB stress fields at levels 0, 2, and 4. The primary weighted jump measure decreased by $97.4\%$, $95.2\%$, and $77.0\%$, respectively, while the load multiplier and maximum yield ratio were unchanged. The selected fields retain the curved toe-to-crest stress transfer as the lower bound increases. RGB values are evaluated directly from the element-local P1 interpolation under common limits, without nodal averaging or image-space smoothing. P1 is used for the selected-field panels; P2 is retained only for the bound-and-cost comparison in Table 5 and Figure 13.

4.3. Frictionless Plane-Strain Extrusion

The third benchmark is Alexander’s complete slip-line solution for frictionless plane-strain extrusion [32]. The lower half of the symmetric domain has upstream and downstream heights of 3 and 1, respectively, giving a 3:1 reduction. The die corner is at $x=6$, and both straight channel extensions have length 6. The material is weightless and rigid-perfectly plastic. Tresca yielding is recovered from Mohr–Coulomb plasticity by setting $\varphi =0$ and $c={\tau }_y$. The inlet ram, symmetry line, and die walls are tangentially smooth, the outlet is traction free, and the objective is the mean compressive inlet pressure p. Alexander’s exact value is

$${\displaystyle \frac{p}{c}}={\displaystyle \frac{4}{3}}\left(1+{\displaystyle \frac{\pi }{2}}\right)=3.427728.$$

(37)

Figure 15 summarizes the half-domain geometry and the boundary conditions shared by both bound calculations. The inlet segments carry the unknown normal reaction in the lower-bound problem and share one normal velocity in the upper-bound problem.

Both bounds started from the same 40-point geometry. The initial lower-bound mesh contained 48 checker-pattern T3 elements, whereas the initial upper-bound network contained 397 admissible lines. Twelve independent adaptive updates were completed. The upper-bound cloud accepted at most 20 new points per update. On the lower-bound side, repeated all-edge refinement was found to spend most new elements within the already localized die region. For this example, at most 80 elements were therefore marked per update and were refined by conforming longest-edge bisection. This calculation control changes neither the yield-dual indicator nor the static admissibility constraints, and the resulting meshes remain nested.

Table 6 lists alternate common states to keep the two sequences on the same update coordinate. The complete 13-state histories are plotted in Figure 16(a). The lower bound increased from 2.982165 to 3.385592, while the upper bound decreased from 3.833333 to 3.511024. At update 12, the two values lie $1.23\%$ below and $2.43\%$ above Eq. (37), respectively. The normalized bracket decreased from $22.20\%$ to $3.57\%$. Figure 16(b) reports cumulative core time after one unrecorded warm-up solve. LBD time contains the primary lower-bound optimization, while UBD time contains both candidate-line generation and upper-bound optimization. The optional fixed-bound stress-selection solves used only for visualization are excluded.

Figure 17 shows four of the common states. The top row superposes the blue $\alpha$ and orange $\beta$ characteristic families on the lower-bound meshes. Refinement remains concentrated near the die entry, while both channel far fields remain coarse. The middle row shows the independently adapted upper-bound points and active DLO lines in the same region. The bottom row gives the lexicographically selected RGB stress fields at the same lower-bound values. The selection reduced the weighted tangential-jump measure from 44.43 to 11.51, 39.46 to 10.17, 58.57 to 22.57, and 67.41 to 13.49 at updates 0, 4, 8, and 12, respectively, without changing $p/c$ or the maximum yield ratio. Thus, the bottom row is another optimum within the fixed-bound admissible set rather than a smoothed image. The lower-bound rows provide static information; the upper-bound row remains a kinematically admissible collapse mechanism.

5. Discussion

The lower-bound guarantee in DLO-LB follows from the stress discretization rather than the candidate-line network. A complete-line inequality controls only an integrated resultant and can admit local over-yielding, as the line-only Prandtl calculation shows. Reversing the DLO objective or replacing velocity jumps with line forces is therefore insufficient, consistent with earlier pseudo lower-bound formulations [26]. In DLO-LB, T3 equilibrium, prescribed and continuous tractions, and pointwise polygonal yield constraints are imposed on the stress field. The intersecting DLO lines serve stress projection and visualization.

DLO-LB and DLO-UB are independent workflows sharing geometry, material, loading, and boundary inputs. DLO-UB returns a kinematically admissible collapse mechanism, whereas DLO-LB returns a conservative capacity and an equilibrated, yield-admissible stress field. Yield-constraint duals trace the Mohr–Coulomb characteristic families and guide lower-bound refinement. The selection in Section 2.6.1 chooses a reproducible member of the possibly non-unique optimal stress set without changing the load multiplier. Similar critical regions support a common failure interpretation, but static characteristics are neither a velocity mechanism nor a constitutive loading path. VDLO instead projects an externally supplied instantaneous stress field onto candidate lines; the external analysis provides any load or time evolution.

Verification is currently limited to homogeneous, plane, associated perfect plasticity with T3 stress spaces. Q4 results remain supplementary because Gauss-point collocation does not ensure element-wide equilibrium, while large candidate-line sets may dominate the cost. Extensions should address equilibrated quadrilateral spaces, scalable line generation, heterogeneous or non-associated materials, three-dimensional equilibrium, finite deformation, and softening. Within this scope, the paired solutions provide conservative capacity, an explicit collapse mechanism, and a numerical bracket.

6. Conclusions

A lower-bound DLO formulation was derived for plane perfect plasticity. Its lower-bound status follows from element-local nodal stresses, exact T3 equilibrium, traction continuity, prescribed stress boundaries, and pointwise inward Mohr–Coulomb constraints. The intersecting DLO lines remain available for stress projection, but complete-line strength constraints alone do not guarantee a lower bound.

The Prandtl, slope, and plane-strain extrusion calculations produced increasing lower bounds and decreasing upper bounds under separate adaptive refinements. The extrusion sequence bracketed the exact slip-line pressure. The associated static characteristics and upper-bound mechanisms occupied similar critical regions, although they remain mechanically different fields. These results are limited to T3 stress spaces and homogeneous, plane, associated perfect plasticity. Equilibrated quadrilateral spaces, scalable line generation, and heterogeneous or non-associated materials remain open.