Multi‐slicing strategy for the three‐dimensional discontinuity layout optimization (3D DLO)

Summary Discontinuity layout optimization (DLO) is a recently presented topology optimization method for determining the critical layout of discontinuities and the associated upper bound limit load for plane two‐dimensional and three‐dimensional (3D) problems. The modelling process (pre‐processing) for DLO includes defining the discontinuities inside a specified domain and building the target function and the global constraint matrix for the optimization solver, which has great influence on the the efficiency of the computation processes and the reliability of the final results. This paper focuses on efficient and reliable pre‐processing of the discontinuities within the 3D DLO and presents a multi‐slicing strategy, which naturally avoids the overlapping and crossing of different discontinuities. Furthermore, the formulation of the 3D discontinuity considering a shape of an arbitrary convex polygon is introduced, permitting the efficient assembly of the global constraint matrix. The proposed method eliminates unnecessary discontinuities in 3D DLO, making it possible to apply 3D DLO for solving large‐scale engineering problems such as those involving landslides. Numerical examples including a footing test, a 3D landslide and a punch indentation are considered, illustrating the effectiveness of the presented method. © 2016 The Authors. International Journal for Numerical and Analytical Methods in Geomechanics published by John Wiley & Sons Ltd.


INTRODUCTION
Topology optimization is a type of mathematical approach for searching for the optimum material layout considering a given space. In engineering fields, topology optimization is commonly used for obtaining the best design considering one or several given target variables such as weight/cost and a series of constraints such as specific loads and supports [1]. With the progress of computing power in the last few decades, presently, it is possible to efficiently analyse large-scale problems with millions of variables on normal personal computers [2][3][4].
Belonging to the family of topology optimization, discontinuity layout optimization (DLO) is a recently presented, elegant and promising numerical method for determining (i) the critical layout of discontinuities in a body at failure and (ii) the associated upper bound limit load [5] for plasticity problems [6][7][8][9]. Introducing large numbers of potential discontinuities permitted to crossover one another in a given domain and specifying the loading and boundary conditions (live loads, dead loads and supports), DLO can automatically generate a set of activated discontinuities corresponding to the upper bound limit load. This method has been successfully used to determine the slip-line mechanisms of several geotechnical problems [7,10,11], the yield-line failure of slabs [12,13], failure mechanisms of masonry structures [14] and of retaining walls in seismic conditions [15].

489
In two-dimensional (2D) DLO, the discontinuities are geometrically represented by line segments, simplifying the pre-processing. Introducing large numbers of grid points into a specific domain, the potential discontinuities are built by connecting each pair of two grid points. In threedimensional (3D) DLO, on the other hand, the discontinuities are geometrically represented by polygons. Overlapping and crossing of different 3D discontinuities should be avoided to save computing resources [8,9], as an important pre-processing task. For the potential wide application of 3D DLO, a highly efficient pre-processing method is necessary.
In this paper, the author presents a multi-slicing strategy for the pre-processing of discontinuities for 3D DLO, naturally avoiding the aforementioned overlap and crossing of different 3D discontinuities. Furthermore, a 3D discontinuity formulation with the shape of an arbitrary convex polygon is introduced, assuring the efficient assembly of the global constraint matrix. Computing effort during pre-processing is saved while the number of unnecessary discontinuities is greatly reduced, paving the way for the analysis of large-scale engineering problems such as landslides with 3D DLO. This paper is organized as follows: in Section 2, DLO formulations for 2D and 3D problems, including the standard pre-processing steps, will be presented, highlighting the necessity for a new method for pre-processing of 3D DLO problems. The multi-slicing strategy and the formulation of 3D discontinuity considering a shape of an arbitrary convex polygon will be presented in Section 3. In Section 4, a footing test, a 3D landslide problem and a punch indentation problem will be investigated, proving the effectiveness of the presented method. This paper closes with concluding remarks given in Section 5.

BASICS OF DISCONTINUITY LAYOUT OPTIMIZATION
From this section on, the author denotes the local variable (vector, matrix) with a subscript, for example d i represents the displacement jump of the discontinuity "i". Meanwhile the symbol without subscript is used to represent a global variable (vector, matrix) such as d D OEd i ; ; d m T being the displacement jump of all the discontinuities in the domain. Furthermore, for better understanding, the dimensions of the equations are provided in the form of L for length and F for force (mass length time 2 ).

Kinematics
The detailed formulation of the DLO was first presented in [7]. For a domain with multiple potential discontinuities shown in Figure 1, several discontinuities become activated (i.e. displacement jump ¤ 0) during loading. When assuming rigid body displacements of the subdomains separated by the activated discontinuities, based on energy conservation, the following relationship is obtained: where W .l/ is the work performed by the live load, W .d / is the work performed by the dead load such as self weight and E is the energy dissipated at the discontinuities. When minimizing the live load W .l/ in the course of limit analysis, the target function of the optimization problem is obtained as In case several discontinuities are connected to one node, because of the rigid displacement of the subdomains, compatibility relationships are formulated at this node. For example, for the three 2D discontinuities i, j and k connected to the node P , with corresponding displacement jump ‡ in the normal and shear directions § d i D OEs i ; n i T , d j D s j ; n j T , and d k D OEs k ; n k T (Figure 2), following relationship is obtained: where X P is the relative displacement of node P , .B P;i d i / is the contribution of discontinuity i to the relative displacement of node P . ¶ B P;i is a coordinate transformation matrix (2 2 for 2D, see [7] for details and 3 3 for 3D, see Appendix A for details) transforming the local displacement jump d i (given in local coordinates) to the global coordinate system [7,8]. In 3D DLO, this compatibility relationship is formulated for every edge [8,9], such as edge PQ illustrated in Figure 3. The compatibility relationships are built at every node (for 2D DLO) or edge (for 3D DLO) in the whole domain. A global compatibility constraint is obtained as ‡ The author would like to emphasize that d i should not be considered as 'displacements of the discontinuity i'. d i is the 'displacement jump' between the two subdomains separated by the discontinuity i. § Rotational failure mechanisms are not considered in this paper, the formulation of which was given in [10]. ¶ For better understanding of Equation (3), the author suggests reverse process considering Figure 2(a)-(b). Because P 1 , P 2 and P 3 in Figure 2(b) come into the same position as P in Figure 2(a), the overall relative displacement of the node P contributed from the discontinuities i, j and k is equal to zero. ; m are the discontinuities of the domain. B is a sparse matrix with row number 2 np (for 2D DLO, np is the total number of nodes) or 3 ne (for 3D DLO, ne is the total number of edges) and column number 2 nd (for 2D DLO) or 3 nd (for 3D DLO), with nd being the total number of discontinuities.

Flow rule considering Mohr-Coulomb failure
Considering the geometrical positions of the discontinuities, there are two types of discontinuities in the domain as (i) boundary discontinuity and (ii) inner discontinuity || (Figure 4).
Considering a discontinuity j with displacement jump as d j D s j ; n j T for 2D or d j D s j ; t j ; n j T for 3D, if j is a boundary discontinuity, d j equals to the absolute displacement. The activations of the boundary discontinuities do not contribute to the dissipated energy E (Figure 1).
If j is an inner discontinuity, on the other hand, d j ¤ 0 indicates the sliding and opening of the subdomains separated by j . Considering an associated flow rule and Mohr-Coulomb failure, and the assumption of the rigid displacement of the subdomains, the following relationship is obtained as 2D: n j tan j js j j D 0 where j is the angle of friction (or dilation) of the inner discontinuity j . ** The authors of [7,8] introduced a plastic multiplier p j D p j;1 ; p j;2 T (for 2D) or p j D OEp j (for 3D) for transforming Equation (5) into the formulations respectively suitable for solution via linear programming (LP) and second-order cone programming (SOCP) problems as follows: 2D: with p j;1 0; p j;2 0: 3D: ; .L/ || It was pointed out in [7] that damage could happen at the supports as well. Therefore, [7] considered discontinuities at supports as inner discontinuities. In this paper, this type of damage is ignored for simplicity. Hence, the types of discontinuities depend only on their geometrical positions. ** Equation (5) indicates n j 0, assuring the two subdomains separated by the discontinuity will not intrude into each other. 492 Y. ZHANG where p j q s 2 j C t 2 j is the conic part in 3D DLO, and where the number of the cones equals the number of inner discontinuities.
Equation (6) indicates .n j D p j tan j / for 3D DLO, as also given in [8]. By substituting n j with . p j tan j / for every inner discontinuities j , great computing resources are saved (less unknowns regarding linear dependencies), bringing the formulation of 3D DLO later in Section 2.4.
For all inner discontinuities in the domain, the following constraint is obtained:

Loading and boundary conditions
As mentioned in Section 2.1, the displacement field of the domain is described by d, with W .d / obtained as where f i is the load vector (dead load) applied on discontinuity i [7]. For work performed by the live load W .l/, the unit displacement constraint [7] is introduced as where F is the unknown limit load .F L 2 /, † † A i is a vector taking into account the area of the discontinuity i (when i is loaded by the live load) and the projection of unit vector parallel to the live load on the discontinuity i. Dissipated energy E is obtained by the plastic multiplier vector as E D ep D e j e k 2 6 4 p j : : : where e j is the dissipated energy vector of the inner discontinuity j . Considering c j as the cohesion of discontinuity j , e j D c j l j ; c j l j (for 2D) with l j being the length of discontinuity j or e j D c j a j (for 3D) with a j being the area of the discontinuity j . As mentioned in Section 2.2, for the boundary discontinuities, the displacement jump equals to the absolute displacement of the discontinuities, which is used for prescribing fixed supports, such as If the discontinuity i is fixed (supported) along the shear direction, set s i D 0 (for 2D) or s i D 0; t i D 0 (for 3D); If the discontinuity i is fixed (supported) along the normal direction, set n i D 0.
Another type of constraint commonly encountered is locked discontinuities. For example, if discontinuities i and j are locked together (d i D d j ), the following constraint is added into the optimization formulation: where I is the unit matrix.

Optimization formulation
Considering Equations (4), (7), (8), (9), (10) and (11), the optimization formulation (written in matrix form) for the 2D DLO problem is obtained as follows: As mentioned before, to save computing resources in 3D DLO, n j D p j tan j is considered for every inner discontinuity j . After redefining d j D s j ; t j ; p j T for every inner discontinuity j , the following optimization formulation (written in matrix form) for 3D DLO is obtained and p j q s 2 j C t 2 j I ! If discontinuity i is fixed along the shear direction: set s i D 0 and t i D 0I ! If discontinuity i is fixed along the normal direction: set n i D 0: (13) where f new and A new are transformed from OE f e and A by considering the redefinition of d j D s j ; t j ; p j T and n j D p j tan j for every inner discontinuity j . Equation (13) is equivalent to the original formulation given in [8].

Y. ZHANG
Equations (12) and (13) can be solved respectively using LP and SOCP solvers. As recommended by [7,8], the author also uses the software 'Mosek' for solving the optimization problems, see [16] for details.

Standard pre-processing in discontinuity layout optimization
A fully connected grid method is used for pre-processing of 2D DLO problems [7]. First, grid points are given in the 2D domain as nodes, with every two nodes connected as a discontinuity. The loading and boundary conditions will be set at the corresponding discontinuities. The activated discontinuities giving minimum limit load F will be automatically obtained after optimization, as illustrated in Figure 5.
This method is very direct and elegant in 2D DLO, with the maximum number of discontinuities equal to C .m; 2/ (m is the number of nodes) considering the mathematical combinations available. On the other hand, when using a similar method for 3D DLO, that is, considering nodes in a 3D domain and connecting every three nodes as a 3D discontinuity, the maximum number of discontinuities becomes C .m; 3/, which dramatically increases with increasing m (Figure 6).

The overlap and crossing of discontinuities in discontinuity layout optimization
In the DLO, the overlap and crossing of the different discontinuities are commonly encountered, as illustrated in Figures 7 and 8.
Overlapping of discontinuities may give incorrect displacement jumps, ‡ ‡ because displacement jumps are 'shared' in overlapped discontinuities. Moreover, overlapping of discontinuities will introduce large numbers of unnecessary discontinuities into the domain. In 2D DLO, distance checking is used to avoid the overlap of discontinuities. Two nodes p 1 and p 2 will only be connected by a discontinuity if jl l 1 l 2 j > ( is a prescribed small value, and p o is all the other nodes except p 1 and p 2 ) for assuring no other nodes are located between p 1 and p 2 (Figure 9).
The presence of discontinuities, which cross over, does not affect the numerical stability of DLO. The fact that discontinuities can cross over one another is a great advantage of this method over other limit analysis methods (e.g. the finite element method [17], when using rigid elements) because of the resulting very wide search space. What should be borne in mind is that crossed discontinuities are not directly connected. Because there is no common node located at the intersection point, the geometrical compatibility relationship (Equation (3)) will not be built for the crossed discontinuities. Therefore, only limited failure patterns are available for crossed discontinuities ( Figure 10). Hence, in 2D DLO, to take into account large numbers of failure patterns, it is necessary to use a fine grid and ensure there are enough nodes and corresponding fully connected discontinuities [7]. Considering the 3D situation, obviously distance checking cannot be used to avoid the overlap of two 3D discontinuities. In addition, introducing more nodes into the 3D domain and building a set of fully connected 3D discontinuities will cause a huge increase in the number of discontinuities ( Figure 6). The only possible routine is to use limited nodes and properly build 3D discontinuities, which are not overlapped or crossed, thereby avoiding unnecessary discontinuities as well as taking into account a large number of possible failure patterns.
Presently, the standard method for dealing with the overlap and crossing of 3D discontinuities is to check every two 3D discontinuities in the domain then add more nodes and split the relevant discontinuities into more discontinuities [8,9] (Figure 11). With such type of pairwise checking, the number of times of this checking (one round) is about C .C .m; 3/ ; 2/, with m being the number of nodes. This process would occupy great computing resources to conduct on a personal computer when the number of the nodes exceeds 100 ( Figure 12). Furthermore, because new discontinuities are introduced during every round of the checking, several rounds of this checking process are necessary within every round, and more discontinuities need to be checked than in the last round. The standard pre-processing of 3D DLO greatly restrains the potential wide application of this promising method. § § A more efficient strategy for modelling 3D discontinuities is essential for 3D DLO.

Geometrical relationship between slicings and three-dimensional discontinuities
A slicing is a plane, cutting the whole domain into two subdomains. Assuming the outside surfaces of a domain are also slicings and considering a domain cut by multiple slicings into several subdomains, the section of the domain of a specific slicing is composed of several convex polygons ( Figure 13). Inspecting the section of a domain considering a slicing i as shown in Figure 14, every line is an intersection line of the slicing i and another slicing and every point is an intersection point of the slicing i and two other slicings. If the convex polygons formed by these lines and points are recognized as discontinuities, overlapping and crossing of the discontinuities are automatically avoided. Moreover, because every edge of the discontinuities is located on the intersection lines of the slicings, the displacement jump of the discontinuities on the corresponding slicings can be transferred to each other considering the geometrical compatibility relationship.
Considering recognizing the convex polygons on one plane from several intersected lines, some highly efficient algorithms are available such as in [18]. In this paper, on the other hand, the author uses a simple algorithm for this problem as illustrated in Figure 15. First, at every node, the intersected half lines are ordered with the polar angle (clockwise suggested by the author, bringing convenience for further building 3D discontinuities). Then, a start node and a search direction are specified, with the search direction turning to the neighbouring polar angle when coming to the next node. Once the search node comes to the start node again, a polygon is recognized. This process is repeated until every polygon has been recognized. § § Wisely, in the numerical examples shown in [8], the domain was first discretized with 3D tetrahedron elements (similar to a discretization used in finite element method), then the surfaces of the 3D elements were taken as discontinuities (the overlap and crossing of discontinuities are avoided) for 3D DLO. Such strategy is suitable for some special cases, and the reliability of the results depend on the discretization of the domain, see the example given in

Three-dimensional discontinuity considering shape of an arbitrary convex polygon
When a convex polygon is recognized, a corresponding 3D discontinuity is defined by a node loop such as .˛;ˇ; ; "/ (where,˛;ˇ; ; " are the id of the nodes), with the direction of the normal displacement jump n determined by the right hand rule (Figure 16). For 3D discontinuities on the same slicing, the author suggests a unified local coordinate system of n, s and t because of the associated easy coding during the programming processes (Appendix A). In addition, for discontinuities on the surfaces (boundary discontinuities), it is very important to assure the direction of n is defined from the outside of the domain point to the inside of the domain.
Considering geometrical compatibility in 3D DLO, there are two directions for each edge (such as in Figure 8, directions P ! Q and Q ! P ) [8]. Hence, the sign should be determined during the assembly of the compatibility matrix B (Section 2.1). In this paper, the sign is determined by

14)
The multi-slicing procedure for 3D DLO is summarized as follows: (1) Multiple slicings are prescribed; (2) Based on the multiple slicings, intersection points of every three non-parallel slicings are calculated and defined as nodes; (3) All nodes on one specific slicing are collected; then all convex polygons on this slicing formed by intersection lines are recognized as 3D discontinuities. This process is repeated until all slicings are checked; (4) With the recognized 3D discontinuities and the constraint given in Equation (14) assembled, the optimization formulation (SOCP problem, see Equation (13)) is built and sent to the optimization solver for solving.

NUMERICAL EXAMPLES
A personal computer with Intel Core i7-3930K 3.20G Hz processor and 16 GB memory was used to obtain the numerical results. The SOCP problem (Equation (13)) is built via the multi-slicing strategy presented then solved by Mosek 7.1 [16].

Footing problem
The weightless footing problem [19] is shown in Figure 17, with the live load 'F ' applied perpendicular on the marked region. When assuming the angle of friction D 0 ı and unified cohesion c for all discontinuities, the analytical solution for the limit loading 'F ' is equal to .2 C /c. Three cases with different slicings are considered ( Figure 18). Key points are the points prescribed for defining the slicings, with every three non-collinear points providing a slicing. ¶ ¶ The details ¶ ¶ The footing problem shown in this paper is actually a 2D problem and only the discontinuities perpendicular to y-ṕ lane will finally become activated. However, during the calculation, the author considers a general condition with all discontinuities given by the slicings taken into account.  of the calculations are given in Table I. |||| Activated discontinuities are shown in Figure 19. The] results indicate the applicability of the presented method. With the increasing number of slicings, the results are improved comparing to the analytical solution, with acceptable time cost for modelling. Theoretically, when introducing very large number of slicings, very accurate results would be obtained. Nevertheless, unlike the LP problem used in 2D DLO, the SOCP problem is much more difficult to solve, taking a long time (Table I) and occupying great computing resources (e.g. around 10 GB memory are occupied by Mosek when solving case III), which increases dramatically with increasing number of slicings. Besides, introducing more slicings will reduce the distance between nodes. When some nodes are 'too close' to each other considering the computational precision of the variables stored on the computer system, linear dependencies will be introduced in the global constraint matrix, greatly affecting the numerical stability of the optimization solver. *** Considering the computational cost and assuring the numerical stability of the solver, presently, a number of slicings less than 100 is suitable, which is unfortunately not enough to obtain a mighty accurate result. For better results, some techniques considering engineering properties of the case should be used, for example, taking advantage of symmetry and introducing more slicings in the region close to the live load.

Three-dimensional landslide
A weightless 3D landslide example is shown in Figure 20, with uniform cohesion c assumed. Owing to symmetry, only half of the problem is considered. The prescribed key points and the correspond-|||| The outside surfaces of the domain are also considered as slicings (six slicings in total) in Table I ing discontinuities are shown in Figure 21. Considering angle of frictions D 0 ı , D 15 ı and D 30 ı , the corresponding activated discontinuities and limit loads of the 3D landslide problem are shown in Figure 22, indicating a large failure zone and the possibility that the limit loads are overestimated.
Considering the failure pattern of the 3D landslide that many activated discontinuities appear at the loaded region, the author implemented a simple grid searching procedure into the multi-slicing strategy to obtain better results. Several standing points are prescribed, then a search point is added, which moves on the grid (Figure 23). The multiple slicings include the surfaces and the slicings defined by the search point and another two standing points (three points are non-collinear). In one search round, after all points on the grid are searched and calculated, the search point with minimum energy/limit load will be added as a new standing point; meanwhile, the corresponding slicings with activated discontinuities are also kept in the next search round. With two 50 50 grids ( Figure 23) and two search rounds on each grid, † † † considering D 30 ı , the obtained activated discontinuities, the standing points and search point in the end of every search round are shown in Figure 24, indicating the results are greatly improved with this search procedure. The final obtained minimum limit loads and failure patterns for the 3D landslide considering D 0 ı , D 15 ı and D 30 ı are shown in Figure 25, which are more reasonable than the results given in Figure 22.

Punch indentation
As a classic example, the punch indentation analysed in [8,20,21] is reanalysed here. A block with its bottom fixed and lateral surfaces restrained along the normal direction (n D 0) is loaded in the centre. Only 1=8 of the volume is considered, owing to symmetry (Figure 26), with uniform cohesion c considered.
As mentioned before, overlapping and crossing problems can also be avoided if discontinuities are transformed from the surfaces of 3D elements of a discretized domain. For testing this strategy, the model is discretized with tetrahedra (similar discretization is used in finite element method analysis in [22]) and the surfaces of the tetrahedra are transformed into discontinuities (Figure 27). The results are listed in Table II, ‡ ‡ ‡ regarding the edge length of the cube equals 1=4 (denoted as cube -1=4), 1=16 (denoted as cube -1=16) and 1=24 (denoted as cube -1=24). From the results, it can be found that one drawback of this strategy is that a finer mesh does not ensure a corresponding much better result, but does take much more computing time.
Then the multi-slicing method with grid searching procedure is used. The initial standing points for the search procedure are shown in Figure 26. In this example, grids similar to cobweb are † † † Two grids are searched, first grid 1 then grid 2. This process was run twice (two search rounds for each grid). ‡ ‡ ‡ Because the CPU time listed in Table II referring the work of [8] was obtained with a different computer, direct comparison could be unfair.  Table II, with the obtained activated discontinuities and deformation pattern regarding D 0 ı and 25 grid shown in Figure 29. It can be seen that using the multi-slicing method combing with the search procedure provides satisfying results even with 5 grid, which takes much less computing time than the case  [20] Multi-block analysis F D 6:561c -F D 104:019c - [21] FEM limit analysis F D 6:051c -F D 66:936c - [8] 3D cube -1=24. Furthermore, the results are continuously improved with finer search grid, indicating the effectiveness and robustness of this method.

CONCLUSION
In this paper, in order to avoid overlapping and crossing of 3D discontinuities, a multi-slicing strategy was presented, with 3D discontinuities being arbitrary convex polygons. The presented method has been proven to be effective, and numerical results are improved with increasing numbers of slicings. Although the number of slicings is still limited because of the complexity at solving large scale SOCP problems, considering the engineering properties of the cases, the author took the advantage of the symmetry and implemented a grid searching procedure, obtaining much better results. The work proposed in this paper opens the door to the application of the 3D DLO to some large scale engineering problems. As mentioned in Section 3, the author suggests a unified local coordinate system for discontinuities on the same slicing, which makes the compatibility matrix for all the discontinuities lying on this slicing the same, bringing great convenience for coding. Figure A.1 illustrates a 3D discontinuity o lying on slicing i with the corresponding local coordinate system s; t; n and global coordinate system x; y;´. When n being perpendicular to the slicing i, s is determined by two farthest nodes n 1 and n 2 lying on the slicing i. Assuming the unit vectors parallel to s; t; n are defined as vs; vt; vn, the coordinates of vs; vt; vn in x; y;´system are .vs x ; vs y ; vs´/, .vt x ; vt y ; vt´/ and .vn x ; vn y ; vn´/. Considering a change of basis, the compatibility matrix of o as B o (Equation (14)) is determined by