A Python script for discontinuity layout optimization

Discontinuity layout optimization (DLO) is a powerful numerical limit analysis technique that can be used to identify the collapse load and associated failure mechanism of a solid or structure. The method successfully automates the traditional ‘upper bound’ method of plasticity, with applications including metal extrusion problems, where die forces are sought, and geotechnical engineering problems, where the stability of foundations or retaining walls are to be established. Notably the basic DLO method uses the same underlying mathematical formulation as ‘ground structure’-based truss layout (or ‘topology’) optimization and is demonstrated in this contribution via a Python script capable of solving plane strain limit analysis problems. Extensions to the basic method are presented to allow treatment of larger-scale problems incorporating cohesive-frictional materials, and with self-weight treated in a new and conceptually elegant way. Finally, various examples are presented to illustrate the capabilities of DLO, with displacement vectors shown to aid interpretation.


Introduction
famously demonstrated that optimal truss structures contain an infinite number of bars with infinitesimal areas, so-called Michell continua. It was later observed that Michell continua are remarkably similar to the so-called slip-line fields associated with plane plasticity problems; both involve special geometric forms known as 'Hencky-Prandtl nets'. Consequently, theories developed for plane plasticity problems were transferred to optimal trusses (e.g., Hemp 1958;Prager 1959).
However, despite the known similarities between optimal truss forms and the forms of plane plasticity failure mechanisms, the efficient numerical "ground structure"based method of identifying optimal truss structures later developed by Dorn et al. (1964) was not applied to plastic analysis problems until comparatively recently, by Smith and Gilbert (2007), who developed the so-called discontinuity layout optimization (DLO) procedure. This is perhaps surprising, given that efficient numerical means of treating plastic limit analysis problems had been sought for many years. For example, the method of characteristics was proposed by Sokolovskii (1965), though this only provides incomplete lower-bound type solutions. Efforts have also been made to enhance limit equilibrium approaches, such as the method of slices for slope stability (e.g., Zhu et al. 2003), but these methods rely on a number of assumptions and the solutions obtained lack formal status. Over the past few decades, finite element limit analysis-based formulations have also been proposed by researchers (e.g., Lysmer 1970;Sloan 1988;Kobayashi 2005;Makrodimopoulos and Martin 2006). Finite element limit analysis formulations usually involve discretization of a body using deformable solid elements, though rigid elements can also be used in conjunction with interface elements placed between solid elements to permit jumps in the stress or strain rate fields (Alwis 2000). When using early finite element limit analysis methods, a priori knowledge of the failure mechanism was often required to obtain accurate solutions, with the user required to tailor the mesh discretization on a case-by-case basis. This was to allow singularities in the stress and/or displacement rate field to be treated in an accurate manner. While the use of adaptive mesh refinement overcomes this to an extent, this 1 3 152 Page 2 of 17 comes at the expense of some complexity and also questions over the remeshing criteria to use remain (e.g., Lyamin et al. 2005). Nevertheless, these methods have the significant advantage of being able to model the failure state without the need to resort to incremental solution schemes, required when using the discrete element method (Cundall and Strack 1979) or non-linear finite elements (De Borst et al. 2012).
On the other hand, by taking advantage of the analogy between slip-line fields and Michell continua, Smith and Gilbert (2007) demonstrated that the numerical layout optimization method developed for truss design problems could be modified so as to able to identify failure mechanisms in plane plasticity problems, with any singularities identified in an entirely natural manner. Also, as with truss layout optimization, use of an adaptive solution scheme (Gilbert and Tyas 2003) means that a very large set of potential discontinuities can be treated, such that highly accurate solutions (e.g., often with ≪ 1% error) can be obtained. To date, several DLO formulations have been proposed, including those capable of identifying in-plane rotational mechanisms (Gilbert et al. 2010;Smith and Gilbert 2013), out-of-plane rotational mechanisms (Gilbert et al. 2014), and three-dimensional translational mechanisms (Hawksbee et al. 2013;Zhang 2017). DLO has to date been applied to numerous applications, ranging from tunnels subject to fire loading (Sun et al. 2019), metal cutting processes (Pritchard et al. 2019), multi-scale masonry analysis (Valentino et al. 2023), and the bearing capacity of volcanic pyroclasts (Galindo et al. 2021), with commercial DLO software tools used not only by industry, but also by researchers seeking to better understand new and longstanding problems alike, e.g., see Leshchinsky and Ambauen (2015), Xie and Leshchinsky (2015), Wang et al. (2017), Liang and Knappett (2017), Zhou et al. (2018), and Zheng et al. (2020). Finally the first textbook on DLO was recently published by Zhang et al. (2022).
Although it has been shown that the DLO procedure can be used to obtain accurate solutions at modest computational cost for many problems, the power of the method still appears under-appreciated by the community. This appears in part to be due to a lack of accessible educational resources for use by researchers and practitioners alike. Although a MATLAB DLO script was previously presented at a specialist conference (Gilbert et al. 2010a), this was limited in that it could only treat small-scale problems involving rectangular problem geometries. Also, self-weight was treated in a somewhat complex manner, not taking advantage of recent research by Smith and Gilbert (2022) that enables complex domain geometries to be handled elegantly. (The alternative simplified approach to treating self-weight for complex domain geometries proposed by Salinas and Zegard (2022) involves a number of approximations, while the use of a supplementary FEM analysis to estimate the effect of body forces at each discontinuity line, recently proposed by Zhang et al. (2022), comes at the expense of significant additional complexity.) Finally, the open-source language Python has been quickly gaining popularity in recent years, particularly in industry circles. These issues are all addressed in the present contribution.
The paper is organized as follows: in Sect. 2, the analogy between truss layout optimization and DLO is discussed and the basic DLO formulation is presented; in Sect. 3 the most important sections of the script are explained; in Sect. 4, the formulation is expanded in order to deal with different yield surfaces, dead loads, and problems involving a large number of nodes, taking advantage of an adaptive solution procedure; in Sect. 5 the expanded formulation is used to solve benchmark problems, showing the potential of the method; finally, in Sect. 6 conclusions are drawn.
2 Analogy with optimal truss layout optimization

Truss layout optimization
To demonstrate the analogy between truss layout optimization and discontinuity layout optimization, it is useful to first revisit the basic truss layout optimization formulation. For a planar truss design problem involving n nodes and m potential truss bars connecting those nodes, the plastic minimum volume truss layout optimization (equilibrium) formulation for a single-load case problem can be written as follows: where V is the volume of the structure, q T = {q + 1 , is a vector containing tensile and compressive forces, each non-negative: c T = {l 1 ∕ 1 , l 1 ∕ 1 , l 2 ∕ 2 , l 2 ∕ 2 ...l m ∕ m } , where l i and i are respectively the length and yield stress of each bar i. B is a suitable (2n × 2m) equilibrium matrix and where f x j and f y j are the x and y components of the external load applied to node j (j = 1, 2...n) . The presence of supports at nodes can be accounted for by omitting the relevant terms from f , together with the corresponding rows from B. Figure 1a-d presents the steps involved in setting up and solving a simple truss layout optimization problem.

Discontinuity Layout Optimization (DLO)
The layout of discontinuities that form at failure in the case of a quasi-statically loaded perfectly cohesive body in plane strain has been demonstrated to be analogous to the layout of bars forming an optimal truss (Smith and Gilbert 2007). Thus, the 'kinematic' slip-line DLO formulation for a body discretized using m nodal connections (slip-line discontinuities) and n nodes can be written as follows: where E is the total internal energy dissipated due to shearing along the discontinuities, punch problem, which results in the same optimal layout as for the truss problem described in Fig. 1a-d; the analogy between truss equilibrium at a node and the compatibility of displacements of bodies sliding relative to each other is illustrated in Fig. 2.
The DLO formulation can, thus, be interpreted as finding the mechanism that dissipates the minimum internal energy for a given imposed displacement u on the system. Conservation of energy then means that the load(s) associated with this displacement can be determined by the equation T = T , albeit the interpretation of and u in this context is not particularly intuitive or convenient. It is, thus, helpful to reformulate formulation (2) to allow easier interpretation and more general usage, as will be described in the next section.
Note that, for convenience, the terms 'energy dissipation' and 'displacement' are herein used as shorthand for 'rate of energy dissipation' and 'displacement rate' (or 'velocity'), respectively.

General form of DLO
While formulation (2) demonstrates the analogy between the truss layout optimization and DLO problem formulations, it is useful to introduce a more general formulation for DLO, which for example also allows the potential for dilational displacements to occur along slip-line discontinuities (see Fig. 3): .N Lm } is a vector containing live loads acting on the discontinuities, d T = {s 1 , n 1 , s 2 , n 2 ...n m } contains relative shear and normal displacements along the discontinuities; is the load factor, such that f T L d in Eq. (3a) is the work done by live loads. Also, p is a vector of non-negative plastic multipliers describing the plastic flow at discontinuities, such that the right-hand side of Eq. (3a) is the internal energy dissipation. Therefore, the objective function Eq. (3a) identifies the minimum live load required to cause plastic flow of the structure (i.e., collapse).
While theoretically live loads can be applied to any discontinuity, in general, they will only be applied to discontinuities lying on free boundaries, such that S Li and N Li will be zero for any non-free boundary i.
It is also important to note that the displacements involved are all relative and that in this paper, the following sign convention is adopted: shear displacements s are taken as positive clockwise (as shown in Fig. 2b) and for normal displacements n, dilational displacements are taken as positive. Thus, 'inward' displacement into a body at a boundary corresponds to dilation at that boundary (where the dilation is acting relative to a fixed 'external' domain). Correspondingly, a normal load at a boundary is considered positive if it is applied inwards with respect to the domain boundary, such that it does positive work. So, for example, if the same positive load is applied to the upper boundary or to the lower boundary, it is oriented downwards or upwards, respectively. Similarly a boundary shear load is considered positive if it acts in an anticlockwise direction around the boundary.
In constraint (3b), the compatibility matrix B i of the ith discontinuity can be written as follows: where i and i are suitable direction cosines for this boundary. Constraint (3c) imposes a flow rule linking displacements s i and n i . For the example shown in Fig. 2b, only shear plastic flow is involved, such that the flow rule for the ith slip line would be written as follows: where the normal plastic flow is set to zero. Note that the flow rule constraint (5) is applied to all internal slip lines. For boundary slip lines, the flow rules need to be modified to satisfy the particular boundary conditions involved.
Note that the use of p 1,i + p 2,i in the work equation (3a), with p 1,i , p 2,i ≥ 0 , ensures that work done is always positive, regardless of the direction of displacement s i . Since Eq. (3a) is being minimized, the flow rule can be viewed as being equivalent p 1,i + p 2,i = |s i | ; this is illustrated in Table 1, which shows a range of possible p 1,i , p 2,i values for cases where s i = 10 or s i = −10 , indicating that the optimal (minimum) value will always occur when p 1,i + p 2,i = |s i | , where one of the values of p will always be zero.
Live load is applied directly on boundary discontinuities. For sake of simplicity, here only a unit (inward, compressive) normal load is used: when using the sign convention adopted in this paper, the indicated relative displacement jump occurs moving across the discontinuity from below to above Table 1 Examples of plastic multiplier values, p 1,i , p 2,i , showing that the optimal (minimum) work value coincides with where is a set containing the loaded boundary discontinuities, and l i is the length of the ith discontinuity. Note that this defines a 'flexible' load, i.e., the loaded boundary line is free to deform.

Boundary conditions
Boundary conditions affect the flow rule and work terms associated with the relevant discontinuities. For a fixed boundary, no additional conditions are required, since the relevant discontinuities have the same properties as internal ones. For discontinuities at free boundaries, the flow rules in Eq. (3c) do not need to be applied, since there is no requirement that the displacements (i.e., s i and n i ) are coupled. In addition, since no internal energy is dissipated on free boundaries, plastic multiplier terms are not included in the work calculation, Eq. (3a). Appendix 1 shows all boundary conditions considered in this work.
In the examples considered in this paper, each boundary is visually represented as follows: -Free boundary: line only -Fixed boundary: cross hatch -Symmetry boundary: dot-dash line -Loaded boundary: directional load arrows

Python implementation of basic DLO formulation
The formulation described in Sect. 2.3 has been programmed in the Python script dlo_basic.

Program code excerpts
Python is an open-source high-level interpreted programming language that is becoming an increasingly popular tool when solving scientific and engineering problems.
The key parts of the formulation described in Sect. 2.3 are now associated with the corresponding program code. Specifically, function DLO performs the high-level steps required to solve a DLO problem.
First, a polygonal problem domain is created: where here the geometrical library shapely is used to generate a polygon using its vertices vt. Nodes and discontinuities are then created: where Nd is a (n × 2) array of nodes, with rows defining the x and y coordinates of nodes. Cn is a (m × 3) array of discontinuities, with each row defining the indices of the connected nodes and the discontinuity length. Note that only nodes and discontinuities lying entirely inside the polygonal domain are created. Also, to remove redundant collinear discontinuities, overlapping connections are filtered out by imposing the following condition: where gcd is a function that finds the greatest common divisor of the x and y increments, dx and dy. Note that for sake of simplicity, only regular Cartesian nodal grids are considered here; alternative nodal grid and connection schemes can, however, be employed if required, e.g., see Zegard and Paulino (2014).
Boundary conditions are then defined: which generates an array bd defining boundary conditions for all discontinuities. A DLO problem is then set up and solved: and results are displayed graphically: Further details of the key steps involved are now presented.

Setting up the DLO optimization problem
As successfully utilized by He et al. (2019), the convex optimization package cvxpy (Diamond and Boyd 2016) is here used to solve the minimization problem Eq. (3), processed via function solveLP. Firstly, all coefficient vectors and matrices in problem (3) are obtained: Boundary conditions are then considered: which generate two vectors activeN and activeG, containing binary data used later to impose boundary conditions when working with the flow rule matrix N and the energy dissipation vector g.
Optimization variables are then created: To improve readability, mathematical expressions are written in as natural a format as possible. Also note that matrix multiplications in cvxpy (from version 1.0) are defined using the symbol "@". The objective function Eq. (3a) is created using: where the array activeG is used as a mask to set certain coefficients in g to 0, to fit the boundary conditions involved.
All constraints in problem (3) are contained in a list cons: The optimization problem can now be created and solved: where here the LP problem is initially solved via the free ECOS solver (Domahidi et al. 2013), which is installed with the cvxpy package. The optimization variables can then be obtained:

Compatibility constraints
To improve computer memory efficiency, the compatibility matrix in Eq. (4) is stored in a sparse matrix. Therefore, it is necessary only to store the values and locations of nonzeros (i.e., row and column identifiers) in matrix B . Since numpy can handle element-wise calculations in arrays, it is convenient to define local compatibility matrices for all discontinuities: where here s and n are respectively index vectors of shear and normal displacements in d ; see Eq. (4). Also n1 and n2 are indices of the first and second nodes connected by discontinuities.
The above script collects all non-zeros in the B matrix and their corresponding locations. To create the sparse matrix, the following function is called: which creates a 2n × 2m sparse matrix using the non-zeros provided.

Flow rule
The flow rule matrix N in Eq. (5) is calculated as follows:

Applied loading
In accordance with Eq. (6), loads can be applied to discontinuities lying along boundaries using the following code fragment: Note that this applies a 'flexible' unit load (since no kinematic constraints have been created to link together the displacements of adjacent loaded segments).

Illustrative examples
The formulation described is now applied to simple literature problems. For all problems, a laptop PC equipped with an Intel I7-7700HQ CPU and running 64-bit Windows 10 was used. These problems are characterized by a pure cohesive model (i.e., one that can be described by the Tresca failure envelope, in which only shear plastic strains occur). Appendix 2 indicates the function name to be called to recreate the following and all other examples presented in the paper.

Metal extrusion
The first example is a classical metal extrusion problem considered by Hill (1950), in which metal is pushed through a rectangular die by a ram, leading to 'steady motion' metal extrusion (i.e., a uniform displacement (rate) field exists at the bottom boundary) (Fig. 4). Figure 5 presents results for three different domain heights. It is evident that a slip-line field similar to that obtained by Hill is only obtained when the loaded boundary is a sufficient distance from the opening. This is because the loading presented in Sect. 3.1.4 is 'flexible,' and does not ensure a uniform displacement field is present at the loaded boundary; this can be addressed by instead using a rigid load, as will be described in Sect. 4.2. Figure 6a shows a variant of the well-known Prandtl punch problem (Hill 1950). Taking symmetry conditions into account, a rectangular domain with 10 × 5 nodal divisions is used here; see Fig. 6b. The computed load factor = 5.222 , which is just 1.56% above the analytical solution of 2 + . An important benefit of the DLO method compared with comparable finite element analysis methods is that the singularity in the displacement field that occurs at the edge of the punch is identified automatically, without the need for e.g., tailored meshes or adaptive mesh refinement.

Extensions
A range of extensions to the basic DLO method are now considered, to increase its range of applicability.

Cohesive-frictional materials
By changing the flow rule matrix, it is possible to treat different convex yield surfaces. For example, it is straightforward to implement the Mohr-Coulomb model for handling cohesive-frictional materials. To achieve this, no changes to the compatibility conditions imposed in Eq. (3b) are required, since normal displacements are already represented in the basic DLO formulation (e.g., to permit the presence of normal displacements at boundaries). However, to implement the Mohr-Coulomb model, the flow rule constraint for the ith slip line needs to be modified to read as follows: where N i is the local plastic flow matrix, p i is a vector containing non-negative plastic multipliers, and is the angle of friction of the material.

Rigid loads
In contrast to the 'flexible' loads defined using Eq. (6), it is possible to specify loads that are rigid, i.e., such that the shape of a given loaded boundary line remains fixed, with discontinuities belonging to the boundary all displacing the same amount. This can be implemented by introducing additional equality constraints that link the displacement variables involved: where linkN and linkS are arrays containing the indices of discontinuities to be linked, considering normal and shear displacements, respectively; also sL and nL are arrays of the corresponding displacement variables. Appendix 1 presents further details of how various types of load can be represented.

Theory
In previous work (e.g Smith and Gilbert 2007), the work done by body forces was implemented by considering the work done moving a column of material that lies, e.g., vertically above a given slip-line discontinuity. For simple examples involving domains with a flat uppermost boundary (e.g., see Fig. 7a), it is relatively easy to calculate the gravity load imposed by materials lying above any discontinuity. However, for general cases, the calculations can become rather complex. For example, in Fig. 7b, since the uppermost edge is non-smooth, any algorithm developed to calculate the gravity load would need to first identify intersection points and vertices along this edge in order to obtain the polygonal areas above a given underlying discontinuity line. Due to the requirement to calculate intersection points, this process can also become computationally expensive when a large number of discontinuities are present in a given DLO problem. For this reason, handling distributed body forces was  Fig. 7 Calculating the body force G associated with the shaded area for a discontinuity AB using the approach described in Smith and Gilbert (2007): a if the top edge is flat then it is straightforward to calculate the area 'above' the discontinuity, e.g., see Gilbert et al. (2010); b however, if a non-smooth top edge is present, then a more complex algorithm is required, as the profile of the top edge must be considered identified as being somewhat challenging when using DLO by He and Gilbert (2016). However, following recent work by Smith and Gilbert (2022), a much simpler and more elegant approach becomes possible. Here, this will be described from a conceptual standpoint.
Consider first a body containing only a non-dilational material within which a translational mechanism is formed (Fig. 8a). Due to conservation of volume, all normal displacements at the domain boundary must sum to zero. To an external observer, material that is displaced at one boundary discontinuity by a normal displacement must 'reappear' at one or more other boundary discontinuities as illustrated for example in Fig. 8b. Here the volume of material pushed downwards at the top of the slope must equal the volume of material pushed outwards on the slope face, i.e., n 1 l 1 + n 5 l 5 = 0 . Given that adopted sign convention for n is dilation positive, the interpretation at boundaries is thus that n 1 is positive and n 5 is negative.
To compute the work done by body forces, it is therefore not necessary to track the movement of material throughout a body (since this is done implicitly by enforcing compatibility elsewhere in the DLO formulation), but to simply note the potential of material that vanishes (positive normal displacement), or appears (negative normal displacement) at a boundary, and to sum these to form the body force work term. Since shear displacements do not affect volume, these need not be considered. Additionally, since the calculation is in terms of a body force energy potential, the computation is independent of the direction of the normal displacement. Hence in Fig. 8b, under normal 1 g gravity, the potential energy change of this global movement of material is equal to − (l 1 n 1 h 1 + l 5 n 5 h 5 ) , where h is the height from an arbitrary datum to the centroid of the discontinuity that is undergoing normal displacement. Since displacement n is assumed to be small, the centroid can be assumed to remain at the mid-point of the discontinuity.
For a material that undergoes volume change on deformation, e.g., dilation, the argument can be extended to include volume generation (or loss) internal to the body, as illustrated in Fig. 8c, with the same principles applying as for the non-dilational material. Here the volume of material pushed downwards on the soil surface plus the volume of dilation on interfaces 2, 3 and 4, must equal the volume of material pushed upwards at the soil surface, i.e., n 1 l 1 + n 2 l 2 + n 3 l 3 + n 4 l 4 + n 5 l 5 = 0 . Under normal 1 g gravity, the potential energy change of this global movement of material is thus equal to -(l 1 n 1 h 1 + l 2 n 2 h 2 + l 3 n 3 h 3 + l 4 n 4 h 4 + l 5 n 5 h 5 ).

General equations
Taking the origin as datum, for an individual discontinuity i, the loss of body force potential P i due to a normal displacement n i is given by: where (x m , y m ) is the coordinate of the mid-point of the discontinuity, is the material unit weight; l i is the discontinuity length; k v and k h are respectively the vertical and horizontal body force accelerations acting in the positive x and y directions (e.g., for normal gravity set k v = −1 ). Smith and Gilbert (2022) provide an analytical proof of this formulation Calculating body force potential energy changes using the current approach based on normal displacements: a example two wedge slope problem, showing dimensions; b displaced wedges for problem with non-dilational material with normal displacements at the boundaries only; c displaced wedges for problem with dilational material with internal and boundary normal displacements. Centroid of dilating volumes are marked with a dot. Normal displacements are exaggerated for illustrative purposes but are considered infinitesimal in the analysis (hence, for example, the centroid positions for interfaces 1 and 5 coincide with the soil surface) in terms of the integral of a body force potential function within the dual equilibrium form of DLO, discussed later in this paper.
When combined with discrete loads applied to any discontinuity, this gives the general equation for loading (live or dead) on a discontinuity i as follows: where S i and N i are respectively shear and normal loads applied to the discontinuity. Body forces may be applied as either live or dead loads. In this paper, only gravity loads are considered as dead loads, and therefore Eq. (9) can be simplified to which involves significantly simpler computations than those associated with the strip model shown in Fig. 7.
To accommodate such body forces it is necessary to extend Eq. (3a) to now include dead loads (D) as follows:

Alternative LP solvers
By default, the LP problem (3) is solved using the opensource solver ECOS (Domahidi et al. 2013) that is distributed with cvxpy. However, cvxpy also supports many other, potentially more efficient, solvers-albeit these need to be installed separately by users. For example, to use the MOSEK solver (MOSEK 2019), the solve command in function solveLP is replaced with the following: where the MOSEK parameter "MSK_IPAR_INTPNT_ BASIS" disables the unnecessary basis identification step to improve speed.

Adaptive solution procedure
Similar to the adaptive 'member adding' procedure proposed for numerical truss layout optimization problems (Gilbert and Tyas 2003;He et al. 2019), an adaptive solution scheme can also be employed when solving DLO problems, significantly improving computational efficiency. Figure 9 shows how a solution is obtained for a Prandtl punch problem when using the adaptive solution process. For sake of simplicity, in this work, all potential discontinuities are created (Fig. 9c), though only small subsets of these are actually used to solve problem (3), e.g., Fig. 9d-i. (Note however that, to improve the memory efficiency further, the step in Fig. 9c can alternatively be omitted, with only required subsets stored.) In the adaptive solution procedure, the dual problem of Eq. (3) is examined. Using duality theory (e.g., see He et al. 2019), the dual problem of Eq. (3), extended in Eq. (11), can be derived as follows: . 9 Solving a DLO problem via the adaptive solution scheme: a problem specification; b node discretization ( 10 × 5 divisions); c create set of potential discontinuities by interconnecting all node pairs (total: 1361); d activate subset of potential discontinuities (215); e iteration 1, solve problem (3), = 6.000 ; f check dual violation using (13), and activate most violating discontinuities via (14); g iteration 2, = 5.333 , with dual violation; h iteration 4, = 5.259 , with dual violation; i final iteration, = 5.222 , no dual violation (total activated discontinuities: 248); j discontinuities with non-zero energy dissipation highlighted and showing displacement vectors obtained using algorithm presented in Sect. 4.6 (taking cohesion c = 1) where t is a vector containing nodal forces, q = [S 1 , N 1 , S 2 , N 2 , ..., S m , N m ] T is a vector of discontinuity forces, and where S i and N i are shear and normal forces acting on discontinuity i of m discontinuities. Note that the inequality constraint Eq. (12c) now defines the Mohr-Coulomb yield-criterion, e.g., for discontinuity i: When primal problem (3) is solved via the primal-dual interior point method, the nodal forces t at every node are obtained from the compatibility constraints. Therefore, the discontinuity force vector q for all potential discontinuities (i.e., discontinuities that are not included in the primal problem) can be obtained from equality constraint Eq. (12b): Since the potential discontinuities are not included in the primal problem, yield condition Eq. (12c) is unlikely to be satisfied for all these discontinuities. However, the degree to which the yield condition is violated can be calculated from where, for each discontinuity, two violation numbers, corresponding to the two inequality constraints in (13), are obtained. These are extracted using The violation for all potential discontinuities can now be sorted, in descending order: The solution obtained in the next iteration is clearly likely to be improved if discontinuities where violation of the yield condition has been identified are added to the primal problem. However, in the interests of computational efficiency, only discontinuities where violation is greatest should be added initially. Taking m v to be the number of potential (i.e., currently non-activated) discontinuities where violation has been found, and m p as the total number of potential discontinuities, the following selection criteria are used to obtain the set of potential discontinuities to be added to the problem at the next iteration: where Δm is the number of discontinuities to be added; 1 and 2 are coefficients determining the percentage of new discontinuities to be added; in this work, both coefficients are taken as 0.05. In Eq. (14), if the number of violated discontinuities is relatively large (i.e., m v ≥ 2 m p ), only a small proportion of the violated discontinuities are added, Eq. (14a). This prevents the problem from growing very rapidly during early iterations of the adaptive solution process, when a large number of dual violations are expected.
On the other hand, if the number of violated discontinuities is relatively small (i.e., m v ≤ 1 2 m p ), all violated discontinuities are added, Eq. (14b), since adding these will only increase the size of the problem slightly. Alternatively, if the number of violated discontinuities is neither large nor small, a fixed proportion of currently non-activated discontinuities are added, Eq. (14c). Finally, if no violation is detected, considering all potential discontinuities (i.e., Δm = m v = 0 ), then no discontinuities need to be added to the problem and the adaptive procedure terminates. Note that although the selection criteria given in Eq. (14) is based on heuristics, and different heuristic strategies are possible, the adaptive procedure will always obtain the same load factor as that obtained by solving the full problem, regardless of the specific heuristics used.

Graphical display of mechanism kinematics
To provide a visual indication of the movements associated with any given failure mechanism, the latter can be overlain with a grid of displacement vectors. This is here achieved by passing horizontal rays across the domain at regular y-intervals. For each ray, the absolute displacement is first set to zero, and then, as the ray passes across the domain from left to right, the absolute displacement is updated as each discontinuity is crossed (by adding the relative displacement occurring at that discontinuity), as shown in Fig. 10. It is convenient to then display vectors of absolute displacement at regular x-and y-intervals; this grid need not coincide with the grid of DLO nodes.

Numerical examples
In this section, a range of examples that illustrate the extensions described in Sect. 4 are presented. All applied loads are 'rigid' unless stated otherwise.

Metal extrusion example revisited
Taking advantage of features introduced in Sect. 4, the metal extrusion example described in Sect. 3.2.1 is now revisited.
First, the ram can now be modeled as a rigid load, ensuring 'steady motion' without the need to set the loaded boundary far from the opening. Also, to obtain more accurate solutions than shown in Fig. 5c, the number of nodal divisions considered can be increased by taking advantage of the adaptive solution scheme described in Sect. 4.5. Figure 11 shows results obtained using 30 × 30 and 60 × 60 nodal divisions while Table 2 provides a summary of the solutions obtained with varying numbers of nodal divisions.
Note that if the adaptive solution procedure is not employed, the open-source ECOS solver fails to obtain solutions when the number of nodal divisions is increased from 15 × 15 to 30 × 30 , though solutions to this problem can be obtained using the more robust MOSEK solver. On the other hand, the ECOS solver is able to solve the problem with 60 × 60 nodal divisions when the adaptive solution scheme is used. Also, increasingly significant reductions in CPU time are observed as problem size increases, when using the MOSEK solver. It is important to note that, although for each problem the CPU cost may vary markedly depending on the solution strategy (i.e., solving  the full problem directly vs using the adaptive solution scheme) and solver (i.e., using ECOS vs using MOSEK), the reported load factors remain the same for any given nodal discretization.

Metal block compressed between rough platens
The second example is a rectangular metal block compressed vertically between two rough platens. The analytical solution to this problem has been presented by Chakrabarty (1991Chakrabarty ( , 2006. As observed by Smith and Gilbert (2007), when considering half the width of the metal block, the optimal layout of discontinuities is identical to the optimal layout of truss bars in a minimum volume cantilever truss, e.g., see Johnson (1961) and Lewiński et al. (1994).
Here, three different aspect ratios are considered, with sample graphical output shown in Fig. 12, and with higher resolution numerical solutions provided in Table 3 to permit comparison with analytical solutions. This indicates that highly accurate solutions can be obtained using the DLO procedure described herein.

Geotechnical examples
Availability of the Mohr-Coulomb failure criterion and the capability to apply body forces makes DLO formulation Eq.
(3) suitable for solving a range of geotechnical engineering problems, examples of which are now presented.

Bearing capacity problem with cohesive-frictional soil
In this example, the purely cohesive material present in the Prandtl punch problem considered in Sect. 3.2.2 is replaced with a cohesive-frictional material to represent a classical geotechnical bearing capacity problem. Here, the Mohr-Coulomb failure criterion is used, with the analytical solution for a weightless soil material obtainable from where c is the cohesion, is the angle of friction, and N q = exp( tan ) tan 2 (45 + ∕2) (see e.g., Smith 2005).
Taking advantage of symmetry, a half domain is here discretized using 48 × 16 nodal divisions, with a unit rigid load applied across 5 nodal divisions. Taking c = 1 and = 25 • , a value of 21.0124 can be computed using the DLO script presented herein, which differs from the exact analytical value of 20.721 by just 1.4%. The corresponding failure mechanism is presented in Fig. 13, which clearly shows how the introduction of friction affects the geometry of the failure mechanism.

Retaining wall
In the case of geotechnical engineering problems, selfweight is often important. Thus, the next example involves a smooth-faced rectangular retaining wall with self-weight taken into account. In this case, the 'passive' condition is considered, where the lateral resistance to movement of the wall into the retained soil body is to be determined. The soil cohesion and unit weight are taken as unity and the angle of friction of the soil = 20 • . The domain is discretized using 40× 20 nodal divisions and a unit load is applied on the left boundary.
In this case the computed load factor = 23.254 is just 0.007% higher than the theoretical value of 23.252 (calculated from 0 = 0.5K p H + 2c . 12 Metal block compressed between rough platens-DLO solutions for: a 10× 10 nodal divisions, aspect ratio = 1, = 2.000 ; b 36× 10 nodal divisions, aspect ratio = 3.6, = 3.325 ; c 67× 10 nodal divisions, aspect ratio = 6.7, = 4.978 (taking cohesion c = 1)  (1857) and Bell (1915)). The failure mechanism shown in Fig. 14 also matches perfectly the expected single 'wedge' failure mechanism for a smooth wall. It is straightforward to analyse the equivalent active condition, by reversing the direction of the live load, as defined in Eq. (6).

Terraced slope with crest surcharge
A key benefit of DLO and other generally applicable numerical methods over hand calculation methods is that they can be applied to problems with arbitrary geometry, where the mode of response is not known in advance. Thus, Fig. 15 depicts a more geometrically complex problem, a terraced slope with crest surcharge. The Mohr-Coulomb failure criterion is used, taking c = 1 and = 25 • , and the unit weight is taken as = 1 . The non-convex domain is discretized with 20 × 12 nodal divisions; nodes lying outside the domain are removed in function createNodes, and discontinuities intersecting concave edges on the top surface are removed in function createDiscontinuities. Since the slope has a non-smooth top surface, calculating the effects of the self-weight using previously described approaches (e.g., Smith and Gilbert 2007;Gilbert et al. 2010) becomes quite involved; see also Fig. 7. However, using the new body force formulation encapsulated in Eq. (8), the load effects can be calculated very efficiently. The computed DLO failure mechanism shown in Fig. 15 includes a piecewise linear-curved slip line, which converts the vertical displacement of the rigid load into horizontal movement of the underlying soil, as observed from the displacement vectors. The example, therefore, demonstrates the flexibility of the DLO method, and its ability to deal with complex geometries and boundary conditions that pose challenges when traditional hand analysis methods are used, given that it will often not be clear a priori whether failure will involve e.g., bearing capacity failure, slope failure, or a combined failure mechanism.

Concluding remarks
Discontinuity layout optimization (DLO) is a powerful numerical limit analysis technique that can be used to automatically identify the collapse load and associated failure mechanism of a solid or structure. Output is in the form of slip-line failure mechanisms that are familiar to engineers. DLO is a member of a wider family of layout optimization procedures that includes the well-known truss layout optimization method originated by Dorn et al. (1964). These methods rely on node and interconnecting lines, forming what is referred to as a 'ground structure' in the case of the truss layout optimization problem.
In this paper, a practical Python implementation of the DLO method is described. This allows the critical failure mechanism and associated load factor to be obtained using linear optimization for a wide range of translational plane To allow readers to run the example problems, an additional script, example.py, has been prepared. This includes functions to run each of the examples considered in this work (Table 5) and to solve all examples considered in the paper in sequence, readers can run the following command in a terminal: In addition, to obviate the need to install Python and associated packages, the scripts have also been made available via Google Colab (Bisong 2019), at: https:// colab. resea rch. google. com/ drive/ 1qgXh lD2JSC_ pKU_ PF0p3 kECFi M9gXT Yd. Entering this url in a web browser automatically runs all examples -scroll to the bottom of the page to view the text and graphical output generated for all problems, which will take a few seconds to appear. This link will work for as long as Google Colab services remain publicly available.  Extrusion12() Fig. 5b Extrusion15() Fig. 5c Extrusion18() Prandtl punch Fig. 6 Prandtl() Extrusion (rigid load) Fig. 11a Extrusion30() Fig. 11b Extrusion60() Metal block (compressed) Fig. 12a Metal10() Fig. 12b Metal36() Fig. 12c Metal67() Bearing capacity Fig. 13 Bearing() Retaining wall (passive failure) Fig. 14 Retaining() Slope (terraced with surcharge) Fig. 15 Slope()