- Email: [email protected]
Systematic identification of polyhedral rock blocks with arbitrary joints and faults
Computers and Geotechnics 29 (2002) 49–72 www.elsevier.com/locate/compgeo
Systematic identification of polyhedral rock blocks with arbitrary joints and faults Jun Lu Department of Mechanical Engineering, Division of Engineering and Applied Science, California Institute of Technology, Pasadena, CA 91125, USA Received 9 October 2000; received in revised form 9 May 2001; accepted 14 May 2001
Abstract The problem of automatic identification of three dimensional rock blocks formed by discontinuities such as joints and faults is studied in this paper. An algorithm and its simulation results for identification of rock blocks are presented. A dynamic link list associated with the algorithm is employed to realize the representation of a polyhedron as well as its forming process by topological identification. Attention is focused on the simplification and the realization of the blocky structure recognition system. Two typical examples are given, one for a blocky structure with all convex polyhedra, and the other with polyhedra of convex and concave geometries. The rock block generator is served as a pre-processor for discrete element method (DEM) using polyhedral elements or for other methods such as block theory (BT) and discontinuous deformation analysis (DDA). The computer simulation can also be extended to rock block identification involving curved joints and faults. # 2001 Published by Elsevier Science Ltd. All rights reserved. Keywords: Rock blocks; system identification; Rock joints and discontinuities
1. Introduction Rock mechanics is facing many challenges since the rapid development in mining and construction engineering, especially fields involving stability of high rock slopes, high arch dam abutments and underground rock excavations. Rock blocks from a geological structure viewpoint can be regarded as a system of blocks cut by planes in space. The geological planes, which serve as discontinuities, could be faults, joints, or cracks. Major structural planes such as faults and large joints are the most important discontinuous planes for rock systems. The geometrical distribution and physical properties of discontinuities considerably affect the mechanical behavior of 0266-352X/02/$ - see front matter # 2001 Published by Elsevier Science Ltd. All rights reserved. PII: S0266-352X(01)00018-0
50
J. Lu / Computers and Geotechnics 29 (2002) 49–72
a jointed rock system [1–4]. Such rock systems are the subjects of study for several computational rock mechanics methods, including discrete element method [5–11], block theory [12] and discontinuous deformation analysis [13]. In practice, a threedimensional rock system is an assembly of many rock blocks of complicated geometries, separated by structural planes representing physical interfaces of the rock system, whose geometrical characteristics cannot be oversimplified. A thorough discontinuous rock mechanical analysis of a rock system requires its complete structure, with its various rock blocks identified. Generally speaking, separating the rock blocks by hand for a three-dimensional system is obviously impractical. Therefore, a computer implementation for identifying rock blocks is a necessary preprocessing tool for three-dimensional discontinuous rock mechanics analysis of systems consisting of complicated topological and geological information. Warburton [14] developed a rock block identification algorithm using infinitely large discontinuity faces, without taking the influence of finite size of faces into consideration. The simple algorithm was not able to treat concave polygons and polyhedra, only convex polyhedra generated after detection. The BGL (Block Generation Language) by Heliot [15] was also limited to convex block detection. The above approach, conceptually simple and easy to program, has been adopted by 3DEC code [16] as a pre-processor for generating rock blocks. The effects of discontinuous faults and joints may be exaggerated due to the assumption that faults and joints are of infinite sizes in space. Discontinuous rock mechanics analysis could give completely different system response if joints and faults of finite sizes are used instead. Lin et al. [17] first described a framework to identify three-dimensional block assemblages for jointed rock masses based on topological concepts. The topological concepts were also used by Jing and Stephansson [18] and Jing [19] for both two- and three-dimensional rock mass identification. Fracture surfaces were assumed planar and finite in three-dimensional space. Attentions were mainly focused on the analytical aspect [17– 19]. The introduction of concepts such as simplex and complex was helpful but not necessary from a practical viewpoint for computer implementation. In addition, no verification by computer simulation of rock masses was reported by the authors. The object of this study is to develop a computer simulation tool capable of generating rock blocks for rock systems with complicated geological and topological structures. To simplify the problem, few topological concepts are used. The algorithm description is simpler compared with that in [17–19]. A universal data structure most suitable for representing polyhedral structures is utilized to facilitate the computer implementation. Block systems generated by computer simulation are presented at the end of the paper.
2. Rock block system and its data structure 2.1. Structural plane representation for rock mass systems Rock mass systems, subject to changes during their long geological histories, evolve to form very complicated geological bodies. They incorporate rock block
J. Lu / Computers and Geotechnics 29 (2002) 49–72
51
assemblies partially or fully separated by some faults, joints and cracks of different sizes and shapes (Fig. 1). Geophysical methods such as borehole, exploratory shaft and test trench are generally used to investigate information about structural planes near or under exposed surfaces of rock systems. Accurate determination of structural surfaces is very limited, and the assumption of flat planes will not adversely affect the accuracy of analysis. In engineering practice, a structural surface is usually treated as a flat plane with or without finite size. In the following sections, structural planes are all treated as three-dimensional mathematical polygons (which may have holes) without thickness. If thickness needs to be considered in mechanical analysis, twin planes replacing the original plane can be used to generate rock systems, with rock blocks in between the twin planes representing materials of faults or joints. In this paper, the terms polyhedron and block are exchangeable. 2.2. The fundamental elements of a polyhedron with orientation Three-dimensional polyhedra can be treated according to their geometrical and topological information, as the geometrical data refer to the spatial positions of the elements of a polyhedron while topological data describe the relationships between various elements of the polyhedron. In this study, the basic elements of a polyhedron are defined including vertices, edges, loops and faces [Fig. 2(a)–(e)]. In addition, the concepts of orientation, associated with every element except vertices, will be widely utilized in the description of a polyhedron. 2.2.1. Directed polyhedra A directed polyhedron (denoted by P) is a finite space enclosed by directed faces. It is a non-null, finite, continuous and closed subset in three-dimensional space. All its elements except its vertices are also directed elements defined as follows.
Fig. 1. The continuous–discontinuous characteristics of rock system at different scales.
52
J. Lu / Computers and Geotechnics 29 (2002) 49–72
Fig. 2. Elements of a directed polyhedron: (a) a directed edges; (b) a directed exterior loop; (c) a directed interior loops; (d) a directed faces; (e) a directed polyhedra; (f) relation among elements of a directed polyhedron.
J. Lu / Computers and Geotechnics 29 (2002) 49–72
53
2.2.2. Directed faces A directed face consists of only one exterior directed loop and may incorporate one or more interior directed loops. It is a non-null, finite, continuous and closed subset in three-dimensional plane. It has a unique orientation given by the normal vector of its exterior loop according to the right-hand rule. Let FðPÞ ¼ ff1 ; f2 ; ; fm g be the set of all directed faces of a polyhedron P. It follows that the boundary of m S P is @P ¼ fi . The normal of a face f points outside of its polyhedron. The boundi¼1
ary @f is a union of directed edges of face f. 2.2.3. Directed loops A loop is said to be directed if its closed boundary incorporates ordered, directed edges such that the end vertex of any of its edge elements is also a starting vertex of its next neighboring edge element. In addition, all its edge components do not intersect with each other except where they meet together at starting and end vertices for neighboring edges in sequence. As shown in Fig. 2(b)–(d), loops can be classified as interior and exterior loops, the former being clockwise and the latter being counterclockwise, assuming the normal of their directed face pointing out of paper. A noteworthy property associated with the description is that if one proceeds along an edge of a directed loop (whether an interior or exterior loop), one will find that his left side always points to the interior of the loop while the right side to the exterior. Let f be a directed face, which incorporates a set of directed loops fc0 ; c1 ; cS cn g (n50), with c0 the prescribed exterior loop. Then it follows that (1) 2 ; :::;S S @f ¼ c0 c1 cn ; (2) f consists of all the points which are interior to c0 and ci ðn5i51Þ including all the points belong to the loops; (3) No intersection occurs between any two loops except at finite number of points. 2.2.4. Directed edges The concept of orientation gives rise to the term directed edge (denoted by e), which is associated with an edge of finite length. For a polyhedron, its edges are the intersection set of its adjacent faces. There are two vertices associated with an edge, one defined as the starting vertex and the other as the end vertex according to the orientation of the edge. Consider a polyhedron P and its vertices set V(P), the edge set E(P) of P is a set such that (1) the vertices of the edges of E(P) belong to V(P); (2) none of the points inside an edge of E(P) is an element of V(Q); (3) every point of an edge of E(P), except the two vertices, belongs to two and only two faces of its polyhedron. 2.2.5. Vertices Obviously, there is no orientation associated with vertices. A vertex can beTrepresentedT as the T intersection of two edges or three faces of a polyhedron, v ¼ e1 e2 or v ¼ f1 f2 f3 . Vertices, edges, loops and faces are closely related as shown in Fig. 2(f). Two useful definitions are given here. The twin face of a directed face f can be formed if all its loops are switched to the reversed directions while the other elements
54
J. Lu / Computers and Geotechnics 29 (2002) 49–72
remain unchanged. The twin edge of an edge e is formed by exchanging the starting and end vertices of e. Owing to their simplicity, their ready visualization and the ease with which they lead themselves to computational simulations, polyhedra have naturally attracted the attention of rock block representation scheme. The above-mentioned elements are the necessary and sufficient terms for the representation of a simply connected or multiconnected polyhedron. The orientation is the most important concept involved in understanding the topological relation between various elements of a polyhedron. 2.3. Data structure representation of spatial rock polyhedra with finite sizes Recall the previous definitions, one can readily see that all the directed elements are their natural extensions used in ordinary geometry. To realize the representation of a polyhedron as well as its forming process by topological identification, a linked list will be used to store its geometrical and topological information. Matrix computation method was utilized by Jing [18] to record polyhedral data information and to keep track of rock block identification process. The problem with this method is that polyhedral data and processing data do not lend themselves easily to matrix storage, and inefficient memory usage is unavoidable. In contrast, a link list of variable length is more efficient and also more straightforward conceptually, as illustrated in Fig. 3. The first stage link list is the one for polyhedra with every node of the list pointing to a polyhedral structure. The second stage link list stores nodes pointing to vertex link lists and face link list. Global coordinates are stored for every indexed vertex of the corresponding polyhedron in the vertex list. The face list has nodes referring to one or more loop link lists depending on whether there is one or more interior loops for the face. A loop list contains its indices of vertices and orientation information. Instead of using concepts of simplex and complex for n-dimensional polyhedra [16, 18], only the concept of orientation is introduced to associate with edges, loops, faces, and polyhedra in three-dimensional space. This simple but important concept is the key for algorithm description, data structure and computer implementation. It will be stated in detail in the following sections how twin edges and twin faces will be used for loop detection and body identification. In the following discussion, without loss of generality, a master rock block (polyhedron), either simply connected or multiconnected, is assumed as the original rock block to be cut. Cutting faces are simply connected or multi-connected spatial polygons (Fig. 4). Rock blocks are generated from the master block by the cutting faces.
3. Identification 3.1. Basic idea of rock mass identification To describe a polyhedron, one would start from its face components, then to its loops followed by edges and vertices. To detect a polyhedron, on the other hand, a reverse process is employed; beginning with vertices and directed edges, then directed
J. Lu / Computers and Geotechnics 29 (2002) 49–72
55
Fig. 3. Dynamic link list used to store topological and geometrical information of a directed polyhedron.
Fig. 4. Rock body and discontinuities: (a) a schematic diagram of a master rock block; (b) a schematic diagram of discontinuous cutting faces.
56
J. Lu / Computers and Geotechnics 29 (2002) 49–72
loops and directed faces, ending with individual directed polyhedral blocks. The general steps of the detection algorithm employed in the following sections are briefly given as follows: 1. Compute all new intersecting vertices from spatial faces and edges and store them into an auxiliary link list along with their original sources. Re-arrange the sequences of all points according to their edge–face and face–face relations. The new structure will serve as a necessary structure for detection in the next step. 2. Introduce a new auxiliary edge link list for each directed face. For instance, two directed edges ab and ba related to an undirected edge ab will be stored into two independent nodes of the edge link list. However, in the case of new vertices being inserted into a boundary edge of a cutting surface, only one directed edge with the orientation the same as that of the boundary edge is added into the link list. 3. Search for directed loops using the data in the auxiliary link list. Find all the possible directed faces thereafter. 4. Introduce another auxiliary link list used to store face information temporarily. A new face associated with a cutting face will be stored as twin directed faces with opposite orientations. In the mean time, any new face associated with a face of the master rock block before being cut will give rise to only one directed face with the same orientation as that of the original face. 5. Deleting any face outside the region of the master rock block before body detection. Then perform body identification, generating link lists and deleting all auxiliary link lists. The flow chart of the algorithm is shown in Fig. 5.
3.2. Formation of new vertices and new directed edges The first step of body identification is to generate vertices owing to three-face and face-to-edge intersections. Assume the polyhedron to be cut has n faces, fi ði ¼ 1; 2; :::nÞ, and there are m cutting planes, fj ðj ¼ n þ 1; :::; n þ mÞ, of finite sizes. 3.2.1. Three-face intersection Possible new points come from intersection among three faces. Fig. 6(a) shows three non-parallel faces fi , fj and fk , at least one of which does not belongs to the master T Trock polyhedron being cut. Then the coordinates of the intersection point fi fj fk is i1 h v ¼ n i p i ; n j p j ; n k p k n Ti ; n Tj ; n Tk
ð1Þ
where v ¼ ðvx ; vy ; vz Þ are the position vector of the intersection point; n i ; n j ; n k are the unit normal of the directed faces fi ; fj and fk respectively; ½n Ti ; n Tj ; n Tk 1 is the inverse of a 3 by 3 matrix whose first, second and third columns consist of three components of vectors n i ; n j and n k respectively; p i ; p j ; p k are position vectors of three points on fi , fj , fk respectively; ðn i p i ; n j p j ; n k p k Þ is a vector formed by inner products of the normal vectors and the position vectors.
J. Lu / Computers and Geotechnics 29 (2002) 49–72
57
Fig. 5. Flow chart of block identification algorithm.
A further step needs to be taken to check if the new intersection point is inside the three faces of finite sizes. A signed arc-length method is used in this study to judge whether a point on the plane of a loop is also interior to the same loop. As shown in Fig. 7, an exterior loop on a spatial face consists of directed edges a1 a2 ; a2 a3 ; :::; ai aiþ1 ; :::; an1 an ; an a1 . Let a0i a0iþ1 be the projective arc of a directed edge ai aiþ1 . Assume the orientation of arc a0i a0iþ1 to be consistent with that of directed edge ai aiþ1 , in other words, if a point q moves along ai aj from ai to aj , its corresponding projection point q0 also proceeds along the arc from a0i to a0j , where a0i and a0j are the intersecting points of lines pai , paj and the unit circle respectively. Therefore, the sum P of the signed arc length is l ¼ ni¼1 a0i a0iþ1 where anþ1 ¼ a1 . It has the following characteristics: (i) l ¼ 2 if p is interior to the loop; (ii) l ¼ 0 if p is exterior to the loop; (iii) l ¼ if p is on the boundary of the loop, i.e. on one of its edges. On the other hand, if an exterior loop with clockwise orientation is considered, assuming same notation for points, a conclusion can be made as follows: (1) l ¼ 2 if p is exterior to the loop; (2) l ¼ 0 if p is interior to the loop; (3) l ¼ if p is on the loop. Finally, a point on the plane of a face is said to be inside the face if it is interior to all the loops of the faces.
58
J. Lu / Computers and Geotechnics 29 (2002) 49–72
3.2.2. Edge-face intersection Owing to the finite sizes of the cutting faces and the faces of the master rock block, new intersecting points may appear along the boundaries of the faces. As shown in Fig. 6(b), the position vector v of a point of intersection between a directed edge e of face fq and face fp is given by *
*
* *
*
*
*
v ¼ v s þ ½ð p n v s n Þ=ð n e^Þ e^
Fig. 6. New vertex generation: (a) three-face intersection; (b) edge-face intersection.
ð2Þ
J. Lu / Computers and Geotechnics 29 (2002) 49–72
59
Fig. 7. Arc-length method used to detect point-loop relation ( an exterior loop is shown here): (a) The point p is interior to the loop; (b) the point p is exterior to the loop; (c) the point p is on the loop.
*
where v s is the position vector of the starting vertex of the edge e, p is any position vector of a point on face fp , e^ is the unit vector of edge e, n is the unit normal of face fp . It is also necessary to check if the point is inside the edge and the face. The procedure is similar and will not be stated in detail. 3.2.3. New vertex recording To record new intersection points, an auxiliary link list PLF is setup to store the points according to different face pairs. This link list was used for both three-face and edge-face intersections. For instance, consider three faces f1 , f2 and f3 , among which, if there exists a common point, then the point will be stored at three different locations having face pairs f1 f2 , f1 f3 , and f2 f3 . As for edge–face intersection, we first find either one or two faces containing the edge being used, then store the intersection point into three different pairsTof the faces involved. For example, let edge e belong to faces f1 and f2 , i.e. e ¼ f1 f2 , then the intersection between e and
60
J. Lu / Computers and Geotechnics 29 (2002) 49–72
T T T T T the third face f3 is v ¼ e f3 ¼ ðf1 f2 Þ f3 ¼ f1 f2 f3 . The last step of storing the point in PLF is the same as that of three-face intersection. 3.2.4. Ordering With the assistance of the PLF list, an unsigned (non-oriented) edge list can thus be generated for each face including the cutting faces and the original faces of the master rock block. So far the new intersecting points along a line contained in the PLF link list have not been properly ordered. Therefore, a simple ordering algorithm can be applied to all the points to obtain ordered lists of points, details are omitted here. 3.2.5. Edge link list setup The last step before performing loop detection is to build an auxiliary edge link list, which will be explained by an example. As shown in Fig. 8(a), face f, a polygon
Fig. 8. Generation of directed edges before loop detection: (a) an example of the directed edges of a face; (b) the auxiliary link list used to store the directed edges of directed faces.
J. Lu / Computers and Geotechnics 29 (2002) 49–72
61
originally having 4 vertices acdf, was cut by another 6 faces, resulting in 10 line segment groups and 13 new points of intersection. The total number of points of interest is 17. Since e and b are on the boundary of face f, ab, bc, cd, de, ef and fa are the directed edges generated by the boundary of face f. Their twin edges are not considered. On the other hand, directed edges formed by at least one interior point among g, h, i, j, k, l, m, n, o, p, q have their twin edges. For instance, directed edge hj has its own twin edge jh. The total number of directed edges is 30. The structure of the link list used to store the edges for all faces is illustrated in Fig. 8(b). The first level pointers are indices for faces. The second level pointers are for storage of the directed edges of each face. Infinite domain is not treated in this study. However, it is a trivial case for block identification and it can be added without difficulty if necessary. 3.3. Loop detection Loop detection uses the edge link list generated in last section. The principle of loop detection is demonstrated graphically by an example in Fig. 9(a)–(g). A directed face shown in Fig. 9(a) has an exterior directed loop abcehja and an interior directed loop qprsp. All corresponding directed edges after cutting are shown in the Fig. 9(b) with an arrow next to each directed edge. Walking along either an interior or exterior directed loop, one will find that the solid domain is kept at the left side of the advancing direction. Before the detection of loops, the flag attached to each directed edge is set to 0 (unmarked), which means that an edge has not been identified so far. If, during the detection process, a directed edge is selected, then the flag associated with the edge is assigned to 1 (marked), prohibiting the edge as a candidate for other loops thereafter. Without loss of generality, we can arbitrarily select a directed edge from the edge link list. For instance, assume that the first edge of the first loop is jk. The tracing proceeds from vertex j to vertex k. The second edge will be started at vertex k according to a maximum-right-handed-angle rule, which will be stated more accurately in the second part of this section. Starting from vertex k, two directed edges ka; kl are the candidate edges. Only one edge according to the maximum-right-handed-angle rule will be selected automatically. In this case, kl forms the largest right-handed angle so it is the second directed edge for the first loop. After setting the flag of the second edge, detection proceeds until the end vertex returns back to vertex j such that a closed directed loop c0 ¼ jk ln mnonij is formed as shown by the bold lines with orientations in Fig. 9(c). Having identified the first loop, a randomly selected edge, from the rest of the directed edges in the edge link list, bc for example, serves as the initial edge of the second loop, cd is the only next connecting edge so it is the second edge of the loop. According to the maximum-right-handed rule, one of the two directed edges started with vertex d is found to be the third edge of the loop. Keep doing this until the end vertex of the loop coincides with the starting vertex of the loop, such that a new loop is detected, c1 ¼ bcdyrsqxghino nmnlkab [Fig. 9(d)].
62
J. Lu / Computers and Geotechnics 29 (2002) 49–72
Fig. 9. Examples of loop detection: (a) an original directed face; (b) the original directed face after cutting; (c) loop c0 ¼ jklnmnonij; (d) loop c1 ¼ bcdyrsqxghinonmnlkab; (e) loop c2 ¼defvwrg0 f0 ug0 e0 g0 vtvfgxpyd; (f) loop c3 ¼ a0 b0 c0 d0 a0 (g) loop c4 ¼ a0 d0 c0 b0 a0 .
J. Lu / Computers and Geotechnics 29 (2002) 49–72
63
A similar searching technique can be applied to the remainder of the unused edges to find three other loops, c2 ¼defvwv0 f0 g0 ug0 e0 g0 vtvfgxpyd, c3 ¼ a0 b0 c0 d0 a0 and c4 ¼ a0 d0 c0 b0 a0 as depicted in Fig. 9(e)–(g). Note the different orientations associated with loop c3 and c4 . The maximum-right-handed-angle rule, a key criterion for loop detection, is given as follows. As shown in Fig. 10, the next possible edges for a loop after an edge e0 are edges ei (i=1, 2,. . .,n) and e00 , the twin edge of e0 with opposite orientation. (1) If n=0, i.e. the only edge e00 will be the next edge of the loop. (2) If n51, let i be the right-handed angle of ei with respect to e0 , then i ¼ cos1 ðe^0 e^i Þ i ¼ þ cos1 ðe^0 e^i Þ
*
if ðe^i e^0 Þn 50 * if ðe^i e^0 Þn < 0
ð3Þ
where n is the unit normal of the face f where the loop locates and is determined by the exterior loop of face f, pointing out of paper, e^0 , e^i are the unit vectors of edges e0 and ei respectively. Then the edge ek associated with max ¼ maxð1 ; 2 ; :::n Þ ¼ k is the next edge of the loop. A noteworthy property associated with the criterion is the case of a degenerate loop consisting of only twin edges. Due to the existence of such degenerate case, the selection of next edge from candidate edges may not be unique. This is because that a starting edge is randomly chosen. For example, for the loop in Fig. 9(e), if the initial edge were fv, then fvwvg0 f0 g0 ug0 e0 g0 vtvf would become one loop defvwvg0 f0 g0 e0 g0 vtvfgxpyd and defgxpyd would be another loop. If vg0 , followed by fv and de, were chosen as the initial edges, then the components and number of loops would change. Three loops will be formed—vg0 f0 g0 ug0 e0 g0 v, fvwvtvf and defgxpyd. However, all
Fig. 10. A schematic diagram of loop search criterion.
64
J. Lu / Computers and Geotechnics 29 (2002) 49–72
adjacent degenerate loops can be grouped together to form a loop of maximum nodes and then attached to a non-degenerate neighboring loop at a common vertex, such that the randomness of initial edges has no influence on any loop structure. For degenerate loops vg0 f0 g0 ug0 e0 g0 v and fvwvtvf shown in Fig. 9(e), loop defgxpyd and fvwvtvf have a common vertex f, so fvwvtvf can join the loop defgxpyd to form a larger loop, which in turn has a common vertex v with another loop vg0 f0 g0 ug0 e0 g0 v such that another larger loop can be formed, which incorporates the components of all degenerated loops. Thus, unless there is no connected neighboring loop for a degenerate loop, it will be joined with its non-degenerated neighbor to form a combined loop. The loop detection algorithm is summarized below: Algorithm (loop detection). Input: a link list of directed edges of a directed face; output: a link list of directed loops of a directed face. Begin: Unmark all edges. Let n be the normal of the directed face. Step 1: If all edges are marked, go to step 4; otherwise, randomly select an unmarked directed edge e0 whose unit directional vector is e^0 . Let e0 be the initial and current edge of the new loop to be detected. Then marked the edge e0 . Step 2: Search among all unmarked directed edges from the edge link list whose starting vertices coincide with the end vertex of current directed edge e0 . (1) The twin edge of edge e0 is denoted e00 if it exists. (2) If there exist n(n51) other non-twin edges, denote the edges ei ði ¼ 1; 2; :::; nÞ and calculate i ði ¼ 1; 2; :::; nÞ, the righthanded angle from the current edge e0 to the edge ei , using Eq. (3). Go to step 3. Step 3: If n51, then calculate max ¼ maxð1 ; :::; n Þ ¼ k . Marked the edge ek corresponds to max and let it be the current edge of the loop. If only the twin edge is the candidate, then marked the twin edge and let it be the current edge (denoted by e0 ) of the loop. If the end vertex of the current edge coincides with the starting vertex of the initial edge during the same loop detection, a new loop has been found, go to loop 1; otherwise, go to step 2; Step 4: Classify interior and exterior loops. For degenerate loops, try to combine them with their nearest neighbors if possible to give an unique solution. End 3.4. Face identification Generally speaking, for both simply connected and multi-connected domains, a face consists of one and only one exterior loop and possibly one or more interior loops. To detect a directed face, the relations between exterior and interior loops are required. The maximum-right-handed-angle rule ensures that there is no intersection among different loops and no self-intersection inside the same loop. It is clear that interior and exterior loops can be classified according to either the normal directions of loops by the right-hand rule. For instance, if the unit normal of a face before the cutting process is n f and the unit normal of a loop c on that face is n c , then we have the following: (1) c is an exterior loop if n f n c ¼ 1; (2) c is an interior loop if n f n c ¼ 1. Assume a face f has total s loops, which can be divided into group of exterior loops {c1 ; c1 ; ::::; cr }ðr4sÞ and group of interior loops {crþ1 ; crþ2 ; ::::; cs }, a
J. Lu / Computers and Geotechnics 29 (2002) 49–72
65
simple procedure is presented for face detection: (1) every ci ði ¼ 1; 2; :::; rÞ forms an exterior loop for a new directed face fi ði ¼ 1; 2; :::; rÞ. So the total number of new directed faces will be r; (2) if any interior loop cj ð j ¼ r þ 1; r þ 2; :::; sÞ is enclosed by any exterior loop ci ði ¼ 1; 2; :::; rÞ, then the loop becomes an interior loop of the face fi ði ¼ 1; 2; :::; rÞ. No numerical calculation is needed at this step except that for determination of enclosure relation among loops. The example illustrated in Fig. 9 shows that 5 loops have been detected. It can be seen that c0 , c4 , c4 and c0 form a new face. c1 , c2 and c3 form 3 new faces independently. 3.5. Body identification After spatial two-dimensional face detection has been completed, a link list of directed faces can then be created to assist in body identification. The body identification is the last major step in block cutting; the first two are directed edge formation and face generation. The complexity of the problem is reduced by a topological superposition process. So it is expected that the last step for body identification should be of the same difficulty in principle as that of the previous face detection process. In addition, the maximum-right-handed-angle criterion for the face detection offers a clue for the body-tracing step, i.e. the criterion can be generalized from the two-dimensional plane to three-dimensional space. The key element of face detection was the unit vectors of directed edges; for the threedimensional case, it is the unit normal of directed faces. Thus the three dimensional problem becomes similar to that of two-dimensional case. A criterion for body tracing is stated as follows: Let f0 be a current directed face during a body detection process. f 00 is the twin face of f0 with opposite normal. Let e be a directed edge of face f0 . Assume fi ði ¼ 1; 2; :::; nÞ has a twin edge of edge e. To find the next adjacent face among the possible fi s and f 00 , a generalized right-handed-angle criterion is given as follows (see Fig. 11): i ¼ cos1 ðn i n 0 Þ i ¼ þ cos1 ðn i n 0 Þ
if ðn 0 n i Þe < 0 if ðn 0 n i Þe50
ð4Þ
where n 0 is the unit normal of f0 and n i the unit normal of fi . 1. If n=0 which means f 00 is the only candidate, then f 00 belongs to the new polyhedron. 2. If n51; calculate max ¼ maxð1 ; 2 :::n Þ ¼ k . The corresponding face fk ð14k4nÞ belongs to the new polyhedron. The only difference between the face detection and the body detection is that for a face detection, an edge only has two connecting points, i.e. the starting vertex and the end vertex while for a body detection, a face has three or more edges which serve as the
66
J. Lu / Computers and Geotechnics 29 (2002) 49–72
Fig. 11. A schematic diagram of block search criterion.
common connecting units between different faces, resulting in a face search in multiple directions. The algorithm is universal for body tracing. For case of multi-connected blocks with cavities, enclosure relations among different closed surfaces need to be checked, which is a relatively simple step. Before proceed, a face link list is needed as an input data structure for the algorithm. The structure of the face link list is similar to that of edge link list. The algorithm is stated more rigorously in the following: Algorithm (polyhedral block detection). Input: link list of directed faces; output: directed polyhedra. Begin Step 1: All faces outside the master rock block are removed. Unmark all the directed faces and all the directed edges. Step 2: If all the faces are marked twice (which implies that all the directed edges are marked), go to step 8; otherwise, randomly select from the link list an unmarked directed face. Let the face be the current directed face f0 of the new polyhedron to be formed. Step 3: Randomly select an unmarked directed edge e of face f0 . Let the edge serve as the current edge. Mark the edge e. Search for directed faces having twin edges of the current edge e. (1) The twin face is denoted f00 if it exists. (2) If there exist n other non-twin faces, denote the faces fi ði ¼ 1; 2; :::; nÞ and calculate i ði ¼ 1; 2; :::; nÞ, the right-handed angle from the current face f0 to the face fi with respect to edge e using Eq. (4). Go to step 4. Step 4: If n51, then calculate max ¼ maxð1 ; :::; n Þ ¼ k . Mark once the face fk with respect to the edge e; otherwise, mark once the face f00 with respect to the edge e. If the newly marked face is new to the new polyhedron to be formed, then add the face to the polyhedral structure. Go to step 5.
J. Lu / Computers and Geotechnics 29 (2002) 49–72
67
Step 5: If not all the edges of face f0 are marked, go to step 3. Otherwise, if all the faces of the new polyhedron are marked twice (i.e. all the directed edges of those faces are marked), go to step 7; otherwise go to step 6. Step 6: Select a directed face that is not marked twice from the face family of the new polyhedron to be formed. Denote the face f0 . Go to step 3. Step 7: A new polyhedron has been identified. Go to step 2. Step 8: Check enclosure relation and multi-connected blocks. End An example is given in Fig. 12 to demonstrate how to remove extraneous fracture faces and crack lines, to which special care must be taken in the removal process to preserve directed faces. For a polyhedron cut by a family of finite-sized faces, ‘‘fracture planes’’ may appear inside blocks after body identification [Fig. 12(a)], a situation similar to that of a polygon cut by a group of line segments, which may lead to ‘‘line marks’’ inside the polygon after face tracing. Each ‘‘fracture plane’’ consists of twin faces of same size, shape and position but of opposite normal directions. At the mean time, every ‘‘line marks’’ incorporates twin edges [Fig. 12(b)]. These kinds of ‘‘fracture planes’’ and ‘‘line marks’’ may or may not be needed in some mechanical analysis. To perform the removal operation, faces inside blocks are removed first [Fig. 12(b)], leaving line marks on other faces which will be removed afterwards
Fig. 12. Removal of ‘‘fracture planes’’ and ‘‘line marks’’: (a) a fractured rock block; (b) the fractured rock block after fracture faces removed; (c) the fractured rock block after fracture faces and line marks removed; (d) process of removing line marks from the top face of the block.
68
J. Lu / Computers and Geotechnics 29 (2002) 49–72
[Fig. 12(c)]. For line mark such as ab in Fig. 12(d), a simple cleanup operation to remove the line and followed by union operation of the two connecting edges if they are linear. However, for a line mark such as cd in Fig. 12(d), after the removal process, the two co-linear faces must be first joined before the union operation of the edges.
4. Examples of convex and concave polyhedra A program has been developed to solve practical problems. The ability of the computer implementation to treat polyhedral rock structures with arbitrary joints and faults is demonstrated by two computer simulations.
Fig. 13. Example one of block identification: (a) a cubic master rock block to be cut by 12 faces of large sizes; (b) the master block and the faces after cutting; (c) a typical face of the master block; (d) a typical cutting face; (e) convex blocks after identification (operation of translational movement is performed for each block) (continued on next page).
J. Lu / Computers and Geotechnics 29 (2002) 49–72
69
Fig. 13. (continued)
4.1. Block identification—cutting a polyhedron by a family of faces of relative large sizes This first example is a special case of convex polyhedra. The original system consists of a cubic master rock block being cut by a group of 12 faces of large sizes depicted before cutting in Fig. 13(a) and after cutting in Fig. 13(b). It results in new block system of all convex polyhedra from the master block after block identification. The result would be the same if faces of infinite size were used because the sizes used here are large enough so that no original edge of the cutting faces is in the master block. A typical face of the master block cut by the group of cutting faces is shown in Fig. 13(c) and a typical cutting face cut by other faces is shown in Fig. 13(d). The convex blocks after identification are displayed in Fig. 13(e). To illustrate the blocks more clearly, the operation of translational movement has been applied to each block with respect to its original position. 4.2. Block identification—cutting a polyhedron by a family of faces of finite sizes The second example is a more general case of both convex and concave polyhedra. The original system consists of a cubic master block and a group of 12 faces of
70
J. Lu / Computers and Geotechnics 29 (2002) 49–72
relatively small and finite sizes as illustrated in Fig. 14(a). Fig. 14(b) shows the actual system after cutting. A typical face of the master block cut by the group of cutting faces is shown in Fig. 14(c) and a typical cutting face cut by other faces is shown in Fig. 14(d). As shown in Fig. 14(e)–(h), rock blocks consist of both convex and concave polyhedra. ‘‘Fracture planes’’ have been removed from the blocks in Fig. 14(g). Both ‘‘fracture planes’’ and ‘‘line marks’’ have been removed from the blocks in Fig. 14(h). The difference between the current cutting process and the previous cutting process in Section 4.1. is that not all the edges of the original cutting faces are outside the master block, resulting more complicated shapes in general. For the case of multi-connected polyhedra, no special treatment is needed for non-cavity multi-connected polyhedra, and only enclosure check is required for processing polyhedra with cavities. The latter is trivial for computer implementation.
Fig. 14. Example two of block identification: (a) a cubic master rock block to be cut by 12 faces of small sizes; (b) the master block and the faces after cutting; (c) a typical face of the master block; (d) a typical cutting face; (e) convex and concave blocks after identification; (f) operation of translational movement is performed for each block; (g) rock blocks whose ‘‘fracture faces’’ have been removed; (h) rock blocks whose ‘‘fracture faces’’ and ‘‘line marks’’ have been removed (continued on next page).
J. Lu / Computers and Geotechnics 29 (2002) 49–72
71
Fig. 14. (continued).
5. Conclusions This paper has presented an algorithm for systematic identification of rock blocks and its computer implementation. A new data structure is introduced with few topological concepts. The identification process can be decomposed into three major steps, namely, edge detection, loop (face) tracing and body identification. Two critical criteria, described in vector forms, are proposed for loop and body identifications; the latter is a generalization of the loop-tracing criterion. Degenerated loops and post-processing of ‘‘fracture faces’’ or ‘‘crack lines’’ are also considered. Not only simply connected blocks can be treated successfully, multi-connected blocks, which turns out to be a trivial case for computer implementation, can also be handled very easily without introducing more complicated concepts and techniques. Simulation examples verify the feasibility of the algorithm employed. The computer implementation, a pre-processor for discontinuous rock mechanics analysis, is therefore capable of dealing with any jointed rock system found in rock engineering, provided that the corresponding geological data are available. Discontinuous mechanical analysis of rock blocks of high arch dam abutments under seismic loading will be presented in another paper.
72
J. Lu / Computers and Geotechnics 29 (2002) 49–72
Acknowledgements The research was conducted with funding provided by the national eighth fiveyear plan of China under grant number 208-01-01-16. The author gratefully acknowledges many suggestions from Dr. Y. Hsu at California Institute of Technology.
References [1] Chappell BA. Deformational response in discontinua. J Rock Mech Min Sci and Geomech Abstr 1979;16(6):377–90. [2] Hoek E, Bray JW. Rock slope engineering. 2nd ed. London: The Institution of Mining and Metallurgy, 1977. [3] Dershowitz WS, Einstein HH. Characterizing rock joint geometry with joint system models. Rock Mech Rock Eng 1988;21(1):21–51. [4] Kulatilake PHSW, Wang S, Stephansson O. Effect of finite-size joints on the deformability of jointed rock in 3-dimensions. Int J Rock Mech Min Sci and Geomech Abstr 1993;30(5):479–501. [5] Cundall PA. Computer model for simulating progressive large scale movements in blocky systems. In: Proceedings of the Symposium of the International Society of Rock Mechanics, Nacy, France, Vol. 1, Paper No. II-8, 1971. [6] Cundall PA. UDEC-A generalized distinct element program for modelling jointed rock. European Research Office and US Army: Peter Associates, 1980. [7] Cundall PA, Hart RD. Development of generalized 2-D and 3-D distinct element programs for modelling jointed rock. ITASCA Consulting Group & US Army Corps of Engineers. Misc Paper SL-85-1, 1985. [8] Dowding CH, Belytschko TB, Yen HJ. A coupled finite element-rigid block method for transient analysis of rock caverns. Int J Num Analy Methods Geomech 1983;7(1):117–27. [9] Lorig LJ, Brady BHG, Cundall PA. Hybrid distinct element-boundary element analysis of jointed rock. Int J Rock Mech Min Sci and Geomech Abstr 1986;23(4):303–12. [10] Cundall PA. Formulation of three-dimensional distinct element model Part I—A scheme to detect and represent contacts in system composed of many polyhedral blocks. Int J Rock Mech Min Sci & Geomech Abstr 1988;25(3):107–16. [11] Hart R, Cundall PA, Lemos J. Formulation of three-dimensional distinct element model Part II— mechanical calculations of a system composed of many polyhedral blocks. Int J Rock Mech Min Sci and Geomech Abstr 1988;25(3):117–25. [12] Goodman RE, Shi GH. Block theory and its application to rock engineering. Prentice-Hall, 1985. [13] Shi GH, Goodman RE. Two dimensional discontinuous deformation analysis. Int J Num Analy Mech Geomech 1989;9(6):541–56. [14] Warburton PM. Application of a new computer model for reconstructing blocky block geometry analysis single block stability and identifying keystones. In: Proc 5th Int Congress on Rock Mechanics, Melbourne, 1983. p. F225–F230. [15] Heliot D. Generating a blocky rock mass. Int J Rock Mech Min Sci and Geomech Abstr 1988;25(3): 127–38. [16] The 3DEC code Manual. ITASCA Consulting Group Ltd, 1994. [17] Lin D, Fairhurst C, Starfield AM. Geometrical identification of three dimensional rock block system using topological techniques. Int J Rock & Geomech Abstr 1987;24(6):331–8. [18] Jing L, Stephansson O. Topological identification of blocky assemblages for jointed rock masses. Technical note. Int J Rock Mech & Geomech Abstr 1994;31(2):163–72. [19] Jing L. Block system construction for three-dimensional discrete element models of fractured rocks. Int J Rock Mech & Geomech Abstr 2000;37(4):645–59.
Systematic identification of polyhedral rock blocks with arbitrary joints and faults Jun Lu Department of Mechanical Engineering, Division of Engineering and Applied Science, California Institute of Technology, Pasadena, CA 91125, USA Received 9 October 2000; received in revised form 9 May 2001; accepted 14 May 2001
Abstract The problem of automatic identification of three dimensional rock blocks formed by discontinuities such as joints and faults is studied in this paper. An algorithm and its simulation results for identification of rock blocks are presented. A dynamic link list associated with the algorithm is employed to realize the representation of a polyhedron as well as its forming process by topological identification. Attention is focused on the simplification and the realization of the blocky structure recognition system. Two typical examples are given, one for a blocky structure with all convex polyhedra, and the other with polyhedra of convex and concave geometries. The rock block generator is served as a pre-processor for discrete element method (DEM) using polyhedral elements or for other methods such as block theory (BT) and discontinuous deformation analysis (DDA). The computer simulation can also be extended to rock block identification involving curved joints and faults. # 2001 Published by Elsevier Science Ltd. All rights reserved. Keywords: Rock blocks; system identification; Rock joints and discontinuities
1. Introduction Rock mechanics is facing many challenges since the rapid development in mining and construction engineering, especially fields involving stability of high rock slopes, high arch dam abutments and underground rock excavations. Rock blocks from a geological structure viewpoint can be regarded as a system of blocks cut by planes in space. The geological planes, which serve as discontinuities, could be faults, joints, or cracks. Major structural planes such as faults and large joints are the most important discontinuous planes for rock systems. The geometrical distribution and physical properties of discontinuities considerably affect the mechanical behavior of 0266-352X/02/$ - see front matter # 2001 Published by Elsevier Science Ltd. All rights reserved. PII: S0266-352X(01)00018-0
50
J. Lu / Computers and Geotechnics 29 (2002) 49–72
a jointed rock system [1–4]. Such rock systems are the subjects of study for several computational rock mechanics methods, including discrete element method [5–11], block theory [12] and discontinuous deformation analysis [13]. In practice, a threedimensional rock system is an assembly of many rock blocks of complicated geometries, separated by structural planes representing physical interfaces of the rock system, whose geometrical characteristics cannot be oversimplified. A thorough discontinuous rock mechanical analysis of a rock system requires its complete structure, with its various rock blocks identified. Generally speaking, separating the rock blocks by hand for a three-dimensional system is obviously impractical. Therefore, a computer implementation for identifying rock blocks is a necessary preprocessing tool for three-dimensional discontinuous rock mechanics analysis of systems consisting of complicated topological and geological information. Warburton [14] developed a rock block identification algorithm using infinitely large discontinuity faces, without taking the influence of finite size of faces into consideration. The simple algorithm was not able to treat concave polygons and polyhedra, only convex polyhedra generated after detection. The BGL (Block Generation Language) by Heliot [15] was also limited to convex block detection. The above approach, conceptually simple and easy to program, has been adopted by 3DEC code [16] as a pre-processor for generating rock blocks. The effects of discontinuous faults and joints may be exaggerated due to the assumption that faults and joints are of infinite sizes in space. Discontinuous rock mechanics analysis could give completely different system response if joints and faults of finite sizes are used instead. Lin et al. [17] first described a framework to identify three-dimensional block assemblages for jointed rock masses based on topological concepts. The topological concepts were also used by Jing and Stephansson [18] and Jing [19] for both two- and three-dimensional rock mass identification. Fracture surfaces were assumed planar and finite in three-dimensional space. Attentions were mainly focused on the analytical aspect [17– 19]. The introduction of concepts such as simplex and complex was helpful but not necessary from a practical viewpoint for computer implementation. In addition, no verification by computer simulation of rock masses was reported by the authors. The object of this study is to develop a computer simulation tool capable of generating rock blocks for rock systems with complicated geological and topological structures. To simplify the problem, few topological concepts are used. The algorithm description is simpler compared with that in [17–19]. A universal data structure most suitable for representing polyhedral structures is utilized to facilitate the computer implementation. Block systems generated by computer simulation are presented at the end of the paper.
2. Rock block system and its data structure 2.1. Structural plane representation for rock mass systems Rock mass systems, subject to changes during their long geological histories, evolve to form very complicated geological bodies. They incorporate rock block
J. Lu / Computers and Geotechnics 29 (2002) 49–72
51
assemblies partially or fully separated by some faults, joints and cracks of different sizes and shapes (Fig. 1). Geophysical methods such as borehole, exploratory shaft and test trench are generally used to investigate information about structural planes near or under exposed surfaces of rock systems. Accurate determination of structural surfaces is very limited, and the assumption of flat planes will not adversely affect the accuracy of analysis. In engineering practice, a structural surface is usually treated as a flat plane with or without finite size. In the following sections, structural planes are all treated as three-dimensional mathematical polygons (which may have holes) without thickness. If thickness needs to be considered in mechanical analysis, twin planes replacing the original plane can be used to generate rock systems, with rock blocks in between the twin planes representing materials of faults or joints. In this paper, the terms polyhedron and block are exchangeable. 2.2. The fundamental elements of a polyhedron with orientation Three-dimensional polyhedra can be treated according to their geometrical and topological information, as the geometrical data refer to the spatial positions of the elements of a polyhedron while topological data describe the relationships between various elements of the polyhedron. In this study, the basic elements of a polyhedron are defined including vertices, edges, loops and faces [Fig. 2(a)–(e)]. In addition, the concepts of orientation, associated with every element except vertices, will be widely utilized in the description of a polyhedron. 2.2.1. Directed polyhedra A directed polyhedron (denoted by P) is a finite space enclosed by directed faces. It is a non-null, finite, continuous and closed subset in three-dimensional space. All its elements except its vertices are also directed elements defined as follows.
Fig. 1. The continuous–discontinuous characteristics of rock system at different scales.
52
J. Lu / Computers and Geotechnics 29 (2002) 49–72
Fig. 2. Elements of a directed polyhedron: (a) a directed edges; (b) a directed exterior loop; (c) a directed interior loops; (d) a directed faces; (e) a directed polyhedra; (f) relation among elements of a directed polyhedron.
J. Lu / Computers and Geotechnics 29 (2002) 49–72
53
2.2.2. Directed faces A directed face consists of only one exterior directed loop and may incorporate one or more interior directed loops. It is a non-null, finite, continuous and closed subset in three-dimensional plane. It has a unique orientation given by the normal vector of its exterior loop according to the right-hand rule. Let FðPÞ ¼ ff1 ; f2 ; ; fm g be the set of all directed faces of a polyhedron P. It follows that the boundary of m S P is @P ¼ fi . The normal of a face f points outside of its polyhedron. The boundi¼1
ary @f is a union of directed edges of face f. 2.2.3. Directed loops A loop is said to be directed if its closed boundary incorporates ordered, directed edges such that the end vertex of any of its edge elements is also a starting vertex of its next neighboring edge element. In addition, all its edge components do not intersect with each other except where they meet together at starting and end vertices for neighboring edges in sequence. As shown in Fig. 2(b)–(d), loops can be classified as interior and exterior loops, the former being clockwise and the latter being counterclockwise, assuming the normal of their directed face pointing out of paper. A noteworthy property associated with the description is that if one proceeds along an edge of a directed loop (whether an interior or exterior loop), one will find that his left side always points to the interior of the loop while the right side to the exterior. Let f be a directed face, which incorporates a set of directed loops fc0 ; c1 ; cS cn g (n50), with c0 the prescribed exterior loop. Then it follows that (1) 2 ; :::;S S @f ¼ c0 c1 cn ; (2) f consists of all the points which are interior to c0 and ci ðn5i51Þ including all the points belong to the loops; (3) No intersection occurs between any two loops except at finite number of points. 2.2.4. Directed edges The concept of orientation gives rise to the term directed edge (denoted by e), which is associated with an edge of finite length. For a polyhedron, its edges are the intersection set of its adjacent faces. There are two vertices associated with an edge, one defined as the starting vertex and the other as the end vertex according to the orientation of the edge. Consider a polyhedron P and its vertices set V(P), the edge set E(P) of P is a set such that (1) the vertices of the edges of E(P) belong to V(P); (2) none of the points inside an edge of E(P) is an element of V(Q); (3) every point of an edge of E(P), except the two vertices, belongs to two and only two faces of its polyhedron. 2.2.5. Vertices Obviously, there is no orientation associated with vertices. A vertex can beTrepresentedT as the T intersection of two edges or three faces of a polyhedron, v ¼ e1 e2 or v ¼ f1 f2 f3 . Vertices, edges, loops and faces are closely related as shown in Fig. 2(f). Two useful definitions are given here. The twin face of a directed face f can be formed if all its loops are switched to the reversed directions while the other elements
54
J. Lu / Computers and Geotechnics 29 (2002) 49–72
remain unchanged. The twin edge of an edge e is formed by exchanging the starting and end vertices of e. Owing to their simplicity, their ready visualization and the ease with which they lead themselves to computational simulations, polyhedra have naturally attracted the attention of rock block representation scheme. The above-mentioned elements are the necessary and sufficient terms for the representation of a simply connected or multiconnected polyhedron. The orientation is the most important concept involved in understanding the topological relation between various elements of a polyhedron. 2.3. Data structure representation of spatial rock polyhedra with finite sizes Recall the previous definitions, one can readily see that all the directed elements are their natural extensions used in ordinary geometry. To realize the representation of a polyhedron as well as its forming process by topological identification, a linked list will be used to store its geometrical and topological information. Matrix computation method was utilized by Jing [18] to record polyhedral data information and to keep track of rock block identification process. The problem with this method is that polyhedral data and processing data do not lend themselves easily to matrix storage, and inefficient memory usage is unavoidable. In contrast, a link list of variable length is more efficient and also more straightforward conceptually, as illustrated in Fig. 3. The first stage link list is the one for polyhedra with every node of the list pointing to a polyhedral structure. The second stage link list stores nodes pointing to vertex link lists and face link list. Global coordinates are stored for every indexed vertex of the corresponding polyhedron in the vertex list. The face list has nodes referring to one or more loop link lists depending on whether there is one or more interior loops for the face. A loop list contains its indices of vertices and orientation information. Instead of using concepts of simplex and complex for n-dimensional polyhedra [16, 18], only the concept of orientation is introduced to associate with edges, loops, faces, and polyhedra in three-dimensional space. This simple but important concept is the key for algorithm description, data structure and computer implementation. It will be stated in detail in the following sections how twin edges and twin faces will be used for loop detection and body identification. In the following discussion, without loss of generality, a master rock block (polyhedron), either simply connected or multiconnected, is assumed as the original rock block to be cut. Cutting faces are simply connected or multi-connected spatial polygons (Fig. 4). Rock blocks are generated from the master block by the cutting faces.
3. Identification 3.1. Basic idea of rock mass identification To describe a polyhedron, one would start from its face components, then to its loops followed by edges and vertices. To detect a polyhedron, on the other hand, a reverse process is employed; beginning with vertices and directed edges, then directed
J. Lu / Computers and Geotechnics 29 (2002) 49–72
55
Fig. 3. Dynamic link list used to store topological and geometrical information of a directed polyhedron.
Fig. 4. Rock body and discontinuities: (a) a schematic diagram of a master rock block; (b) a schematic diagram of discontinuous cutting faces.
56
J. Lu / Computers and Geotechnics 29 (2002) 49–72
loops and directed faces, ending with individual directed polyhedral blocks. The general steps of the detection algorithm employed in the following sections are briefly given as follows: 1. Compute all new intersecting vertices from spatial faces and edges and store them into an auxiliary link list along with their original sources. Re-arrange the sequences of all points according to their edge–face and face–face relations. The new structure will serve as a necessary structure for detection in the next step. 2. Introduce a new auxiliary edge link list for each directed face. For instance, two directed edges ab and ba related to an undirected edge ab will be stored into two independent nodes of the edge link list. However, in the case of new vertices being inserted into a boundary edge of a cutting surface, only one directed edge with the orientation the same as that of the boundary edge is added into the link list. 3. Search for directed loops using the data in the auxiliary link list. Find all the possible directed faces thereafter. 4. Introduce another auxiliary link list used to store face information temporarily. A new face associated with a cutting face will be stored as twin directed faces with opposite orientations. In the mean time, any new face associated with a face of the master rock block before being cut will give rise to only one directed face with the same orientation as that of the original face. 5. Deleting any face outside the region of the master rock block before body detection. Then perform body identification, generating link lists and deleting all auxiliary link lists. The flow chart of the algorithm is shown in Fig. 5.
3.2. Formation of new vertices and new directed edges The first step of body identification is to generate vertices owing to three-face and face-to-edge intersections. Assume the polyhedron to be cut has n faces, fi ði ¼ 1; 2; :::nÞ, and there are m cutting planes, fj ðj ¼ n þ 1; :::; n þ mÞ, of finite sizes. 3.2.1. Three-face intersection Possible new points come from intersection among three faces. Fig. 6(a) shows three non-parallel faces fi , fj and fk , at least one of which does not belongs to the master T Trock polyhedron being cut. Then the coordinates of the intersection point fi fj fk is i1 h v ¼ n i p i ; n j p j ; n k p k n Ti ; n Tj ; n Tk
ð1Þ
where v ¼ ðvx ; vy ; vz Þ are the position vector of the intersection point; n i ; n j ; n k are the unit normal of the directed faces fi ; fj and fk respectively; ½n Ti ; n Tj ; n Tk 1 is the inverse of a 3 by 3 matrix whose first, second and third columns consist of three components of vectors n i ; n j and n k respectively; p i ; p j ; p k are position vectors of three points on fi , fj , fk respectively; ðn i p i ; n j p j ; n k p k Þ is a vector formed by inner products of the normal vectors and the position vectors.
J. Lu / Computers and Geotechnics 29 (2002) 49–72
57
Fig. 5. Flow chart of block identification algorithm.
A further step needs to be taken to check if the new intersection point is inside the three faces of finite sizes. A signed arc-length method is used in this study to judge whether a point on the plane of a loop is also interior to the same loop. As shown in Fig. 7, an exterior loop on a spatial face consists of directed edges a1 a2 ; a2 a3 ; :::; ai aiþ1 ; :::; an1 an ; an a1 . Let a0i a0iþ1 be the projective arc of a directed edge ai aiþ1 . Assume the orientation of arc a0i a0iþ1 to be consistent with that of directed edge ai aiþ1 , in other words, if a point q moves along ai aj from ai to aj , its corresponding projection point q0 also proceeds along the arc from a0i to a0j , where a0i and a0j are the intersecting points of lines pai , paj and the unit circle respectively. Therefore, the sum P of the signed arc length is l ¼ ni¼1 a0i a0iþ1 where anþ1 ¼ a1 . It has the following characteristics: (i) l ¼ 2 if p is interior to the loop; (ii) l ¼ 0 if p is exterior to the loop; (iii) l ¼ if p is on the boundary of the loop, i.e. on one of its edges. On the other hand, if an exterior loop with clockwise orientation is considered, assuming same notation for points, a conclusion can be made as follows: (1) l ¼ 2 if p is exterior to the loop; (2) l ¼ 0 if p is interior to the loop; (3) l ¼ if p is on the loop. Finally, a point on the plane of a face is said to be inside the face if it is interior to all the loops of the faces.
58
J. Lu / Computers and Geotechnics 29 (2002) 49–72
3.2.2. Edge-face intersection Owing to the finite sizes of the cutting faces and the faces of the master rock block, new intersecting points may appear along the boundaries of the faces. As shown in Fig. 6(b), the position vector v of a point of intersection between a directed edge e of face fq and face fp is given by *
*
* *
*
*
*
v ¼ v s þ ½ð p n v s n Þ=ð n e^Þ e^
Fig. 6. New vertex generation: (a) three-face intersection; (b) edge-face intersection.
ð2Þ
J. Lu / Computers and Geotechnics 29 (2002) 49–72
59
Fig. 7. Arc-length method used to detect point-loop relation ( an exterior loop is shown here): (a) The point p is interior to the loop; (b) the point p is exterior to the loop; (c) the point p is on the loop.
*
where v s is the position vector of the starting vertex of the edge e, p is any position vector of a point on face fp , e^ is the unit vector of edge e, n is the unit normal of face fp . It is also necessary to check if the point is inside the edge and the face. The procedure is similar and will not be stated in detail. 3.2.3. New vertex recording To record new intersection points, an auxiliary link list PLF is setup to store the points according to different face pairs. This link list was used for both three-face and edge-face intersections. For instance, consider three faces f1 , f2 and f3 , among which, if there exists a common point, then the point will be stored at three different locations having face pairs f1 f2 , f1 f3 , and f2 f3 . As for edge–face intersection, we first find either one or two faces containing the edge being used, then store the intersection point into three different pairsTof the faces involved. For example, let edge e belong to faces f1 and f2 , i.e. e ¼ f1 f2 , then the intersection between e and
60
J. Lu / Computers and Geotechnics 29 (2002) 49–72
T T T T T the third face f3 is v ¼ e f3 ¼ ðf1 f2 Þ f3 ¼ f1 f2 f3 . The last step of storing the point in PLF is the same as that of three-face intersection. 3.2.4. Ordering With the assistance of the PLF list, an unsigned (non-oriented) edge list can thus be generated for each face including the cutting faces and the original faces of the master rock block. So far the new intersecting points along a line contained in the PLF link list have not been properly ordered. Therefore, a simple ordering algorithm can be applied to all the points to obtain ordered lists of points, details are omitted here. 3.2.5. Edge link list setup The last step before performing loop detection is to build an auxiliary edge link list, which will be explained by an example. As shown in Fig. 8(a), face f, a polygon
Fig. 8. Generation of directed edges before loop detection: (a) an example of the directed edges of a face; (b) the auxiliary link list used to store the directed edges of directed faces.
J. Lu / Computers and Geotechnics 29 (2002) 49–72
61
originally having 4 vertices acdf, was cut by another 6 faces, resulting in 10 line segment groups and 13 new points of intersection. The total number of points of interest is 17. Since e and b are on the boundary of face f, ab, bc, cd, de, ef and fa are the directed edges generated by the boundary of face f. Their twin edges are not considered. On the other hand, directed edges formed by at least one interior point among g, h, i, j, k, l, m, n, o, p, q have their twin edges. For instance, directed edge hj has its own twin edge jh. The total number of directed edges is 30. The structure of the link list used to store the edges for all faces is illustrated in Fig. 8(b). The first level pointers are indices for faces. The second level pointers are for storage of the directed edges of each face. Infinite domain is not treated in this study. However, it is a trivial case for block identification and it can be added without difficulty if necessary. 3.3. Loop detection Loop detection uses the edge link list generated in last section. The principle of loop detection is demonstrated graphically by an example in Fig. 9(a)–(g). A directed face shown in Fig. 9(a) has an exterior directed loop abcehja and an interior directed loop qprsp. All corresponding directed edges after cutting are shown in the Fig. 9(b) with an arrow next to each directed edge. Walking along either an interior or exterior directed loop, one will find that the solid domain is kept at the left side of the advancing direction. Before the detection of loops, the flag attached to each directed edge is set to 0 (unmarked), which means that an edge has not been identified so far. If, during the detection process, a directed edge is selected, then the flag associated with the edge is assigned to 1 (marked), prohibiting the edge as a candidate for other loops thereafter. Without loss of generality, we can arbitrarily select a directed edge from the edge link list. For instance, assume that the first edge of the first loop is jk. The tracing proceeds from vertex j to vertex k. The second edge will be started at vertex k according to a maximum-right-handed-angle rule, which will be stated more accurately in the second part of this section. Starting from vertex k, two directed edges ka; kl are the candidate edges. Only one edge according to the maximum-right-handed-angle rule will be selected automatically. In this case, kl forms the largest right-handed angle so it is the second directed edge for the first loop. After setting the flag of the second edge, detection proceeds until the end vertex returns back to vertex j such that a closed directed loop c0 ¼ jk ln mnonij is formed as shown by the bold lines with orientations in Fig. 9(c). Having identified the first loop, a randomly selected edge, from the rest of the directed edges in the edge link list, bc for example, serves as the initial edge of the second loop, cd is the only next connecting edge so it is the second edge of the loop. According to the maximum-right-handed rule, one of the two directed edges started with vertex d is found to be the third edge of the loop. Keep doing this until the end vertex of the loop coincides with the starting vertex of the loop, such that a new loop is detected, c1 ¼ bcdyrsqxghino nmnlkab [Fig. 9(d)].
62
J. Lu / Computers and Geotechnics 29 (2002) 49–72
Fig. 9. Examples of loop detection: (a) an original directed face; (b) the original directed face after cutting; (c) loop c0 ¼ jklnmnonij; (d) loop c1 ¼ bcdyrsqxghinonmnlkab; (e) loop c2 ¼defvwrg0 f0 ug0 e0 g0 vtvfgxpyd; (f) loop c3 ¼ a0 b0 c0 d0 a0 (g) loop c4 ¼ a0 d0 c0 b0 a0 .
J. Lu / Computers and Geotechnics 29 (2002) 49–72
63
A similar searching technique can be applied to the remainder of the unused edges to find three other loops, c2 ¼defvwv0 f0 g0 ug0 e0 g0 vtvfgxpyd, c3 ¼ a0 b0 c0 d0 a0 and c4 ¼ a0 d0 c0 b0 a0 as depicted in Fig. 9(e)–(g). Note the different orientations associated with loop c3 and c4 . The maximum-right-handed-angle rule, a key criterion for loop detection, is given as follows. As shown in Fig. 10, the next possible edges for a loop after an edge e0 are edges ei (i=1, 2,. . .,n) and e00 , the twin edge of e0 with opposite orientation. (1) If n=0, i.e. the only edge e00 will be the next edge of the loop. (2) If n51, let i be the right-handed angle of ei with respect to e0 , then i ¼ cos1 ðe^0 e^i Þ i ¼ þ cos1 ðe^0 e^i Þ
*
if ðe^i e^0 Þn 50 * if ðe^i e^0 Þn < 0
ð3Þ
where n is the unit normal of the face f where the loop locates and is determined by the exterior loop of face f, pointing out of paper, e^0 , e^i are the unit vectors of edges e0 and ei respectively. Then the edge ek associated with max ¼ maxð1 ; 2 ; :::n Þ ¼ k is the next edge of the loop. A noteworthy property associated with the criterion is the case of a degenerate loop consisting of only twin edges. Due to the existence of such degenerate case, the selection of next edge from candidate edges may not be unique. This is because that a starting edge is randomly chosen. For example, for the loop in Fig. 9(e), if the initial edge were fv, then fvwvg0 f0 g0 ug0 e0 g0 vtvf would become one loop defvwvg0 f0 g0 e0 g0 vtvfgxpyd and defgxpyd would be another loop. If vg0 , followed by fv and de, were chosen as the initial edges, then the components and number of loops would change. Three loops will be formed—vg0 f0 g0 ug0 e0 g0 v, fvwvtvf and defgxpyd. However, all
Fig. 10. A schematic diagram of loop search criterion.
64
J. Lu / Computers and Geotechnics 29 (2002) 49–72
adjacent degenerate loops can be grouped together to form a loop of maximum nodes and then attached to a non-degenerate neighboring loop at a common vertex, such that the randomness of initial edges has no influence on any loop structure. For degenerate loops vg0 f0 g0 ug0 e0 g0 v and fvwvtvf shown in Fig. 9(e), loop defgxpyd and fvwvtvf have a common vertex f, so fvwvtvf can join the loop defgxpyd to form a larger loop, which in turn has a common vertex v with another loop vg0 f0 g0 ug0 e0 g0 v such that another larger loop can be formed, which incorporates the components of all degenerated loops. Thus, unless there is no connected neighboring loop for a degenerate loop, it will be joined with its non-degenerated neighbor to form a combined loop. The loop detection algorithm is summarized below: Algorithm (loop detection). Input: a link list of directed edges of a directed face; output: a link list of directed loops of a directed face. Begin: Unmark all edges. Let n be the normal of the directed face. Step 1: If all edges are marked, go to step 4; otherwise, randomly select an unmarked directed edge e0 whose unit directional vector is e^0 . Let e0 be the initial and current edge of the new loop to be detected. Then marked the edge e0 . Step 2: Search among all unmarked directed edges from the edge link list whose starting vertices coincide with the end vertex of current directed edge e0 . (1) The twin edge of edge e0 is denoted e00 if it exists. (2) If there exist n(n51) other non-twin edges, denote the edges ei ði ¼ 1; 2; :::; nÞ and calculate i ði ¼ 1; 2; :::; nÞ, the righthanded angle from the current edge e0 to the edge ei , using Eq. (3). Go to step 3. Step 3: If n51, then calculate max ¼ maxð1 ; :::; n Þ ¼ k . Marked the edge ek corresponds to max and let it be the current edge of the loop. If only the twin edge is the candidate, then marked the twin edge and let it be the current edge (denoted by e0 ) of the loop. If the end vertex of the current edge coincides with the starting vertex of the initial edge during the same loop detection, a new loop has been found, go to loop 1; otherwise, go to step 2; Step 4: Classify interior and exterior loops. For degenerate loops, try to combine them with their nearest neighbors if possible to give an unique solution. End 3.4. Face identification Generally speaking, for both simply connected and multi-connected domains, a face consists of one and only one exterior loop and possibly one or more interior loops. To detect a directed face, the relations between exterior and interior loops are required. The maximum-right-handed-angle rule ensures that there is no intersection among different loops and no self-intersection inside the same loop. It is clear that interior and exterior loops can be classified according to either the normal directions of loops by the right-hand rule. For instance, if the unit normal of a face before the cutting process is n f and the unit normal of a loop c on that face is n c , then we have the following: (1) c is an exterior loop if n f n c ¼ 1; (2) c is an interior loop if n f n c ¼ 1. Assume a face f has total s loops, which can be divided into group of exterior loops {c1 ; c1 ; ::::; cr }ðr4sÞ and group of interior loops {crþ1 ; crþ2 ; ::::; cs }, a
J. Lu / Computers and Geotechnics 29 (2002) 49–72
65
simple procedure is presented for face detection: (1) every ci ði ¼ 1; 2; :::; rÞ forms an exterior loop for a new directed face fi ði ¼ 1; 2; :::; rÞ. So the total number of new directed faces will be r; (2) if any interior loop cj ð j ¼ r þ 1; r þ 2; :::; sÞ is enclosed by any exterior loop ci ði ¼ 1; 2; :::; rÞ, then the loop becomes an interior loop of the face fi ði ¼ 1; 2; :::; rÞ. No numerical calculation is needed at this step except that for determination of enclosure relation among loops. The example illustrated in Fig. 9 shows that 5 loops have been detected. It can be seen that c0 , c4 , c4 and c0 form a new face. c1 , c2 and c3 form 3 new faces independently. 3.5. Body identification After spatial two-dimensional face detection has been completed, a link list of directed faces can then be created to assist in body identification. The body identification is the last major step in block cutting; the first two are directed edge formation and face generation. The complexity of the problem is reduced by a topological superposition process. So it is expected that the last step for body identification should be of the same difficulty in principle as that of the previous face detection process. In addition, the maximum-right-handed-angle criterion for the face detection offers a clue for the body-tracing step, i.e. the criterion can be generalized from the two-dimensional plane to three-dimensional space. The key element of face detection was the unit vectors of directed edges; for the threedimensional case, it is the unit normal of directed faces. Thus the three dimensional problem becomes similar to that of two-dimensional case. A criterion for body tracing is stated as follows: Let f0 be a current directed face during a body detection process. f 00 is the twin face of f0 with opposite normal. Let e be a directed edge of face f0 . Assume fi ði ¼ 1; 2; :::; nÞ has a twin edge of edge e. To find the next adjacent face among the possible fi s and f 00 , a generalized right-handed-angle criterion is given as follows (see Fig. 11): i ¼ cos1 ðn i n 0 Þ i ¼ þ cos1 ðn i n 0 Þ
if ðn 0 n i Þe < 0 if ðn 0 n i Þe50
ð4Þ
where n 0 is the unit normal of f0 and n i the unit normal of fi . 1. If n=0 which means f 00 is the only candidate, then f 00 belongs to the new polyhedron. 2. If n51; calculate max ¼ maxð1 ; 2 :::n Þ ¼ k . The corresponding face fk ð14k4nÞ belongs to the new polyhedron. The only difference between the face detection and the body detection is that for a face detection, an edge only has two connecting points, i.e. the starting vertex and the end vertex while for a body detection, a face has three or more edges which serve as the
66
J. Lu / Computers and Geotechnics 29 (2002) 49–72
Fig. 11. A schematic diagram of block search criterion.
common connecting units between different faces, resulting in a face search in multiple directions. The algorithm is universal for body tracing. For case of multi-connected blocks with cavities, enclosure relations among different closed surfaces need to be checked, which is a relatively simple step. Before proceed, a face link list is needed as an input data structure for the algorithm. The structure of the face link list is similar to that of edge link list. The algorithm is stated more rigorously in the following: Algorithm (polyhedral block detection). Input: link list of directed faces; output: directed polyhedra. Begin Step 1: All faces outside the master rock block are removed. Unmark all the directed faces and all the directed edges. Step 2: If all the faces are marked twice (which implies that all the directed edges are marked), go to step 8; otherwise, randomly select from the link list an unmarked directed face. Let the face be the current directed face f0 of the new polyhedron to be formed. Step 3: Randomly select an unmarked directed edge e of face f0 . Let the edge serve as the current edge. Mark the edge e. Search for directed faces having twin edges of the current edge e. (1) The twin face is denoted f00 if it exists. (2) If there exist n other non-twin faces, denote the faces fi ði ¼ 1; 2; :::; nÞ and calculate i ði ¼ 1; 2; :::; nÞ, the right-handed angle from the current face f0 to the face fi with respect to edge e using Eq. (4). Go to step 4. Step 4: If n51, then calculate max ¼ maxð1 ; :::; n Þ ¼ k . Mark once the face fk with respect to the edge e; otherwise, mark once the face f00 with respect to the edge e. If the newly marked face is new to the new polyhedron to be formed, then add the face to the polyhedral structure. Go to step 5.
J. Lu / Computers and Geotechnics 29 (2002) 49–72
67
Step 5: If not all the edges of face f0 are marked, go to step 3. Otherwise, if all the faces of the new polyhedron are marked twice (i.e. all the directed edges of those faces are marked), go to step 7; otherwise go to step 6. Step 6: Select a directed face that is not marked twice from the face family of the new polyhedron to be formed. Denote the face f0 . Go to step 3. Step 7: A new polyhedron has been identified. Go to step 2. Step 8: Check enclosure relation and multi-connected blocks. End An example is given in Fig. 12 to demonstrate how to remove extraneous fracture faces and crack lines, to which special care must be taken in the removal process to preserve directed faces. For a polyhedron cut by a family of finite-sized faces, ‘‘fracture planes’’ may appear inside blocks after body identification [Fig. 12(a)], a situation similar to that of a polygon cut by a group of line segments, which may lead to ‘‘line marks’’ inside the polygon after face tracing. Each ‘‘fracture plane’’ consists of twin faces of same size, shape and position but of opposite normal directions. At the mean time, every ‘‘line marks’’ incorporates twin edges [Fig. 12(b)]. These kinds of ‘‘fracture planes’’ and ‘‘line marks’’ may or may not be needed in some mechanical analysis. To perform the removal operation, faces inside blocks are removed first [Fig. 12(b)], leaving line marks on other faces which will be removed afterwards
Fig. 12. Removal of ‘‘fracture planes’’ and ‘‘line marks’’: (a) a fractured rock block; (b) the fractured rock block after fracture faces removed; (c) the fractured rock block after fracture faces and line marks removed; (d) process of removing line marks from the top face of the block.
68
J. Lu / Computers and Geotechnics 29 (2002) 49–72
[Fig. 12(c)]. For line mark such as ab in Fig. 12(d), a simple cleanup operation to remove the line and followed by union operation of the two connecting edges if they are linear. However, for a line mark such as cd in Fig. 12(d), after the removal process, the two co-linear faces must be first joined before the union operation of the edges.
4. Examples of convex and concave polyhedra A program has been developed to solve practical problems. The ability of the computer implementation to treat polyhedral rock structures with arbitrary joints and faults is demonstrated by two computer simulations.
Fig. 13. Example one of block identification: (a) a cubic master rock block to be cut by 12 faces of large sizes; (b) the master block and the faces after cutting; (c) a typical face of the master block; (d) a typical cutting face; (e) convex blocks after identification (operation of translational movement is performed for each block) (continued on next page).
J. Lu / Computers and Geotechnics 29 (2002) 49–72
69
Fig. 13. (continued)
4.1. Block identification—cutting a polyhedron by a family of faces of relative large sizes This first example is a special case of convex polyhedra. The original system consists of a cubic master rock block being cut by a group of 12 faces of large sizes depicted before cutting in Fig. 13(a) and after cutting in Fig. 13(b). It results in new block system of all convex polyhedra from the master block after block identification. The result would be the same if faces of infinite size were used because the sizes used here are large enough so that no original edge of the cutting faces is in the master block. A typical face of the master block cut by the group of cutting faces is shown in Fig. 13(c) and a typical cutting face cut by other faces is shown in Fig. 13(d). The convex blocks after identification are displayed in Fig. 13(e). To illustrate the blocks more clearly, the operation of translational movement has been applied to each block with respect to its original position. 4.2. Block identification—cutting a polyhedron by a family of faces of finite sizes The second example is a more general case of both convex and concave polyhedra. The original system consists of a cubic master block and a group of 12 faces of
70
J. Lu / Computers and Geotechnics 29 (2002) 49–72
relatively small and finite sizes as illustrated in Fig. 14(a). Fig. 14(b) shows the actual system after cutting. A typical face of the master block cut by the group of cutting faces is shown in Fig. 14(c) and a typical cutting face cut by other faces is shown in Fig. 14(d). As shown in Fig. 14(e)–(h), rock blocks consist of both convex and concave polyhedra. ‘‘Fracture planes’’ have been removed from the blocks in Fig. 14(g). Both ‘‘fracture planes’’ and ‘‘line marks’’ have been removed from the blocks in Fig. 14(h). The difference between the current cutting process and the previous cutting process in Section 4.1. is that not all the edges of the original cutting faces are outside the master block, resulting more complicated shapes in general. For the case of multi-connected polyhedra, no special treatment is needed for non-cavity multi-connected polyhedra, and only enclosure check is required for processing polyhedra with cavities. The latter is trivial for computer implementation.
Fig. 14. Example two of block identification: (a) a cubic master rock block to be cut by 12 faces of small sizes; (b) the master block and the faces after cutting; (c) a typical face of the master block; (d) a typical cutting face; (e) convex and concave blocks after identification; (f) operation of translational movement is performed for each block; (g) rock blocks whose ‘‘fracture faces’’ have been removed; (h) rock blocks whose ‘‘fracture faces’’ and ‘‘line marks’’ have been removed (continued on next page).
J. Lu / Computers and Geotechnics 29 (2002) 49–72
71
Fig. 14. (continued).
5. Conclusions This paper has presented an algorithm for systematic identification of rock blocks and its computer implementation. A new data structure is introduced with few topological concepts. The identification process can be decomposed into three major steps, namely, edge detection, loop (face) tracing and body identification. Two critical criteria, described in vector forms, are proposed for loop and body identifications; the latter is a generalization of the loop-tracing criterion. Degenerated loops and post-processing of ‘‘fracture faces’’ or ‘‘crack lines’’ are also considered. Not only simply connected blocks can be treated successfully, multi-connected blocks, which turns out to be a trivial case for computer implementation, can also be handled very easily without introducing more complicated concepts and techniques. Simulation examples verify the feasibility of the algorithm employed. The computer implementation, a pre-processor for discontinuous rock mechanics analysis, is therefore capable of dealing with any jointed rock system found in rock engineering, provided that the corresponding geological data are available. Discontinuous mechanical analysis of rock blocks of high arch dam abutments under seismic loading will be presented in another paper.
72
J. Lu / Computers and Geotechnics 29 (2002) 49–72
Acknowledgements The research was conducted with funding provided by the national eighth fiveyear plan of China under grant number 208-01-01-16. The author gratefully acknowledges many suggestions from Dr. Y. Hsu at California Institute of Technology.
References [1] Chappell BA. Deformational response in discontinua. J Rock Mech Min Sci and Geomech Abstr 1979;16(6):377–90. [2] Hoek E, Bray JW. Rock slope engineering. 2nd ed. London: The Institution of Mining and Metallurgy, 1977. [3] Dershowitz WS, Einstein HH. Characterizing rock joint geometry with joint system models. Rock Mech Rock Eng 1988;21(1):21–51. [4] Kulatilake PHSW, Wang S, Stephansson O. Effect of finite-size joints on the deformability of jointed rock in 3-dimensions. Int J Rock Mech Min Sci and Geomech Abstr 1993;30(5):479–501. [5] Cundall PA. Computer model for simulating progressive large scale movements in blocky systems. In: Proceedings of the Symposium of the International Society of Rock Mechanics, Nacy, France, Vol. 1, Paper No. II-8, 1971. [6] Cundall PA. UDEC-A generalized distinct element program for modelling jointed rock. European Research Office and US Army: Peter Associates, 1980. [7] Cundall PA, Hart RD. Development of generalized 2-D and 3-D distinct element programs for modelling jointed rock. ITASCA Consulting Group & US Army Corps of Engineers. Misc Paper SL-85-1, 1985. [8] Dowding CH, Belytschko TB, Yen HJ. A coupled finite element-rigid block method for transient analysis of rock caverns. Int J Num Analy Methods Geomech 1983;7(1):117–27. [9] Lorig LJ, Brady BHG, Cundall PA. Hybrid distinct element-boundary element analysis of jointed rock. Int J Rock Mech Min Sci and Geomech Abstr 1986;23(4):303–12. [10] Cundall PA. Formulation of three-dimensional distinct element model Part I—A scheme to detect and represent contacts in system composed of many polyhedral blocks. Int J Rock Mech Min Sci & Geomech Abstr 1988;25(3):107–16. [11] Hart R, Cundall PA, Lemos J. Formulation of three-dimensional distinct element model Part II— mechanical calculations of a system composed of many polyhedral blocks. Int J Rock Mech Min Sci and Geomech Abstr 1988;25(3):117–25. [12] Goodman RE, Shi GH. Block theory and its application to rock engineering. Prentice-Hall, 1985. [13] Shi GH, Goodman RE. Two dimensional discontinuous deformation analysis. Int J Num Analy Mech Geomech 1989;9(6):541–56. [14] Warburton PM. Application of a new computer model for reconstructing blocky block geometry analysis single block stability and identifying keystones. In: Proc 5th Int Congress on Rock Mechanics, Melbourne, 1983. p. F225–F230. [15] Heliot D. Generating a blocky rock mass. Int J Rock Mech Min Sci and Geomech Abstr 1988;25(3): 127–38. [16] The 3DEC code Manual. ITASCA Consulting Group Ltd, 1994. [17] Lin D, Fairhurst C, Starfield AM. Geometrical identification of three dimensional rock block system using topological techniques. Int J Rock & Geomech Abstr 1987;24(6):331–8. [18] Jing L, Stephansson O. Topological identification of blocky assemblages for jointed rock masses. Technical note. Int J Rock Mech & Geomech Abstr 1994;31(2):163–72. [19] Jing L. Block system construction for three-dimensional discrete element models of fractured rocks. Int J Rock Mech & Geomech Abstr 2000;37(4):645–59.