Application of a vertex chain operation algorithm on topological analysis of three-dimensional fractured rock masses

Abstract

Identifying the morphology of rock blocks is vital to accurate modelling of rock mass structures. This paper applies the concepts of directed edges and vertex chain operations which are typical for block tracing approach to block assembling approach to construct the structure of three-dimensional fractured rock masses. Polygon subtraction and union algorithms that rely merely on vertex chain operation are proposed, which allow a fast and convenient construction of complex faces/loops. Apart from its robustness in dealing with finite discontinuities and complex geometries, the advantages of the current methodology in tackling some challenging issues associated with the morphological analysis of rock blocks are addressed. In particular, the identification of complex blocks with interior voids such as cavity, pit and torus can be readily achieved based on the number and the type of loops. The improved morphology visualization approach can benefit the pre-processing stage when analyzing the stability of rock masses subject to various engineering impacts using the block theory and the discrete element method.

Keywords

morphology / block assembling / vertex operation / discontinuities

Introduction

The natural rock masses normally contain various discontinuities that are generated subject to the long and complex tectonic history. These discontinuities separate the original intact rocks into many sub-blocks, creating a complex discontinuous system []. In particular, the failure mode of the rock masses is dominated in the forms of block falling, sliding, rotation as well as opening/closure of discontinuities. Therefore, the primary task to evaluate the stability of rock mass subject to engineering disturbances is to determine the topographical structure of the discontinuous system.

The practice of structure modelling of rock masses can date back to Warburton [] who proposed a rock block identification algorithm by assuming that all geological planes are infinite in size. This algorithm was later extended by Heliot [] taking into account the effects of tectonic history. Despite of its simplicity and ease for computational implementation, this approach could only handle convex blocks, and severely overestimate the connectivity and size of rock blocks. In the past decades, great efforts have been devoted to investigating the complex topology of real rock blocks with various approaches developed. These approaches, in general, can be categorized into the block tracing approach (BTA) and the block assembling approach (BAA):

• The block tracing approach (BTA) constructs a dimensional object from the bottom (i.e., points) up to the top (e.g., polyhedrons) based on the topology coherence theory. The implementation of BTA can be summarized as shown in Fig. 1(a): 1) generating the discontinuous network within the analysis domain; 2) tracing all the vertexes and edges from the intersections between geological planes; 3) determining the topological relationships between the vertexes and edges and identifying the possible loops (closed faces) based on the theory of topology; 4) constructing the possible polyhedrons from the faces and identifying their convexity. Steps 3 and 4 are achieved by combination of vertex chains. Successful applications of BTA to systematically identify the rock mass structure include Lin and Fairhurst [], Ikegawa and Hudson [], Jing [], Lu [] and Elmouttie et al. [–].

• The block assembling approach (BAA) identifies a dimensional object from the top (i.e., polyhedrons) down to the bottom (points) [–]. The general procedure of BAA normally involves four steps (see Fig. 1 (b)): 1) generating the fracture/discontinuous network within the analysis domain; 2) temporarily ignoring the finiteness of discontinuities by extending the discontinuities to intersect with the boundaries of the analysis domain; 3) identifying all the convex element blocks; 4) restoring the exact finiteness of the discontinuities and merging the element blocks to form the real complex blocks. The key parts of BAA are the detection and combination of element blocks.

Both methods have their intrinsic advantages and disadvantages. The concept of BTA is straightforward and adheres to the basic principles and procedures of topology theory, i.e., construct the space object, from vertexes to edges, from edges to faces, until to the final polyhedral model. However, complex topological analysis and mathematical calculation are required to identify closed loops. On the other hand, the BAA is based on the idea of spatial partition followed by block assemblage. The partition procedure of BAA is very simple but the assemblage and the topology analysis of element blocks are quite complex.

This paper extended the previous work of Zhang and Lei [] by incorporating the concepts of oriented edges [] in defining the loops of element and complex blocks, and boundary chain operations [] when constructing complex blocks from element blocks. Polygon subtraction and union algorithms which are capable of handling rock blocks with arbitrary shapes are detailed. The feasibility and correctness of the new approach are demonstrated through a generic example of rock systems.

Data structure of the method

A complete procedure of the new approach is described in Fig. 2. For the ease of operation, the data structures for the rock block system need to be built prior to the implementation procedure of BAA. Similar to the traditional BAA, the data structure of the current approach comprises the vertex class, the fracture class, the element face class, the element block class, the complex face class and the complex block class. The internal structure of these classes is similar to Zhang and Lei [], except that some topological indexes are incorporated for vertex chain operation and the vertex IDs are arranged in directed sequences. The purposes of incorporating these topological indexes will be introduced later in Section 3.

Vertex class

As shown in Fig. 3, the only attributes of vertex class are the three Cartesian coordinate variables: x, y and z.

Fracture class

For simplicity, the fractures are normally treated as circular finite planes in BAA. As Fig. 3 shows, in addition to a geometry component defining the geometry, a mechanical component defining the strength and a VertexLists component which contains the coordinates of boundary vertexes of the fracture plane, the internal data structure for the fracture class also comprises a boundary box data component (Xrange, Yrange, Zrange) and two spatial indexes (Upperindex and Lowerindex). The boundary box data component defines the spatial range of the fracture (Fig. 4), while the two spatial indexes reflect the spatial relationship between the element block and the fracture plane.

There are several points needed to be clarified when constructing the fracture data array:

1) The fracture plane is mathematically determined by:

(1)Ax+By+Cz=D,

where A=\sin \alpha \sin \beta, B=\sin \alpha \cos \beta, C=\cos \alpha, D={Ax}_{\mathrm{0}}+{By}_{\mathrm{0}}+{Cz}_{\mathrm{0}} with (x₀, y₀, z₀) being the coordinate of the center point of the fracture plane and the definitions of α and β refer to Fig. 3.

2) For computational convenience, the normal vector is always pointing to the upper half space of the fracture plane.

3) If the angle between the line connecting the vertex and the center of the fracture plane and the normal of the fracture is less than 90°, the vertex is in the upper space of the fracture; if this angle is equal to 90°, the vertex is on the fracture plane; otherwise, the vertex is beneath the fracture plane.

Elementface class

The faces created by assuming infinite size of fracture planes are defined as element face. An element face is always convex and lies on a single fracture plane. Hence, the geometrical and strength properties of an element face are inherited from its parent fracture plane through the FractureID. As shown in Fig. 5, the element face class also includes a VertexLists attribute which lists the coordinates of the boundary vertexes and a Topo attribute which stores the spatial relationship between the element face and the element block that contains the element face, an Area attribute (the area of the elementface), and a boundary box attribute (Xrange, Yrange, Zrange) defining the spatial range of the element face (see Fig. 4). Note that vertexes in the VertexLists are stored in an anticlockwise order with respect to the normal direction. The advantage of using such a convention will be explained later. Calculating the area of an element face can refer to Ref. [].

Elementblock class

The blocks identified with the assumption of infinite fracture planes are called element block which is convex in nature. Since each element block is made up of simple element faces, the internal structure of the element block data is based on the internal structure of element face data. As shown in Fig. 5, the internal structure of the Elementblock class contains: 1) ElementFaceLists directing to the data structure of the element faces comprising the element block; 2) Boundary Box data (Xrange, Yrange, Zrange) and 3) volume data (volume of the element block). Determination of the volume of an element block can refer to Ref. []. The Topo value of the element face mentioned in Section 2.3 is defined as: when the normal vector of the element face pointing to the inner space of the element block, the Topo value of the element face is+1; otherwise, it will be-1. For example, in the case shown in Fig. 5, the normal of F5 pointing to the inner space of the element block and the topo value of F5 is+1, while the Topo value of F3 is-1 as the normal vector of F3 pointing to the outer space of the element block.

Complexface class

When restoring the finiteness of the fracture planes, some of the element faces may be merged to form larger faces. These faces may have abnormal or irregular shapes and are called complex faces. As illustrated in Fig. 6, the geometry of a complex face could be much more complicated than an element face with some regular or irregular “holes” (empty space) inside due to the finiteness of the fracture plane. Therefore, following Lu [], apart from the basic features of the element face class, two additional attributes are introduced to describe a complex face: a) OuterLists which lists the vertexes of the outer boundary in an anticlockwise manner, and b) InnerLists which sorts the vertexes of the inner boundary clockwisely. It should be noted that there should only be one set of vertex list for the outer boundary, while the number of inner boundary vertex list could be any number in accordance with the number of internal “holes”.

Complexblock class

The last step of BAA is to assemble the real blocks from the “virtual” element blocks accounting for the exact finiteness of fracture planes. These blocks may own complex geometries and are thus categorized as complex block. As Fig. 6 shows, the internal structure of the Complexblock class is composed of: 1) ComplexFaceList referring to the data of all the complex faces forming the complex block; 2) OuterList pointing to the data of the complex faces forming the outer surface of the complex block; 3) InnerList referring to the data of the complex faces forming the inner surface of the complex block; 4) SubElementblockId stores the ID of the element blocks comprising the complex block and 5) Volume of the complex block which can be calculated from the volume summation of its element blocks.

Implementation of an improved morphological visualization method

Implementation of a BAA analysis usually involves: 1) identification of element blocks; 2) union operation of the element blocks; 3) construction of complex blocks and 4) modification of complex blocks. The last step is invoked only when excavation is considered.

Identification of element blocks

(a) Basic idea

The first step in a BAA analysis is to identify the element blocks by temporarily neglecting the finiteness of the fracture planes. The main procedure of the identification of element blocks is similar to traditional BAA analysis as illustrated in Fig. 7. Two strategies are employed to enhance the computational efficiency:

a) instantiate the computation domain into an assembly of convex element blocks. In the case when the computation domain is convex, it is taken as the first element block (EB1); otherwise, it could be decomposed into a number of convex subdomains and taken as an assembly of convex element blocks;

b) successively employ the fracture planes in a sequence of descending area (e.g., J2, J1, J3, and J4 in Fig. 7) to cut the existing element blocks and identify the intersection between the fracture plane and the existing element blocks. If the fracture plane is infinitely large, its area is taken as the part within the analysis domain. By adopting a descending area sequence, the number of generated element blocks and the times to perform intersection calculation can be reduced. This can be illustrated by a simple 2D example as shown in Fig. 8. A rectangular domain is cut by four fracture planes which are different in scale. If fracture cutting is performed in an arbitrary order shown in Fig. 8 (a), there will be nine element blocks generated at the end and the total times of intersection identification will be 1+2+4+7=14. Nevertheless, if fracture cutting is conducted in a sequence of descending area (Fig. 8 (b)), the total number of generated blocks will be reduced to seven and the total times of intersection calculation will be reduced to 1+2+3+5=11.

c) check whether the boundary box of the element block and that of the fracture plane intersect. If not (Fig. 9(a)), go to the next fracture plane.

A storage space will be allocated to the newly-created element blocks, while the information of the reduced existing element blocks will be updated in the Elementblock data array. Only if intersection calculation has been carried out with respect to all the existing element blocks, fracture cutting can be continued for the next fracture until all fracture planes have been considered.

If the fracture plane intersects the element block (Fig. 9(b)), further judgement is required to identify the spatial relationship between the fracture plane the element block based on the relationship matrix of the vertexes of the element block with respect to the fracture plane:

(2)\left[R\right]=\left[\begin{array}{c}\hfill A(x_{\mathrm{1}}-x_{\mathrm{0}})+B(y_{\mathrm{1}}-y_{\mathrm{0}})+C(z_{\mathrm{1}}-z_{\mathrm{0}})\hfill \\ \hfill A(x_{\mathrm{2}}-x_{\mathrm{0}})+B(y_{\mathrm{2}}-y_{\mathrm{0}})+C(z_{\mathrm{2}}-z_{\mathrm{0}})\hfill \\ \hfill M\hfill \\ \hfill A(x_n-x_{\mathrm{0}})+B(y_n-y_{\mathrm{0}})+C(z_n-z_{\mathrm{0}})\hfill \end{array}\right]=\left[\begin{array}{c}\hfill R_{\mathrm{1}}\hfill \\ \hfill R_{\mathrm{2}}\hfill \\ \hfill M\hfill \\ \hfill R_n\hfill \end{array}\right]

where the physical meaning and expression for A, B, C , x_{\mathrm{0}}, y_{\mathrm{0}} and z_{\mathrm{0}} can refer to Eq. 1, x_i, y_i and z_i are the coordinates of the ith vertex of the element block, R_i represents the spatial relationship between the ith vertex and the fracture plane. Note that the fracture plane is assumed to be infinite in Eq. 2. If all components of R are greater than or equal to 0 (EB1 in Fig. 9 (c)), the element block is in the upper space of the fracture plane; if all components ofR are less than or equal to 0 (EB2 in Fig. 9 (c)), the element block is in the lower space of the fracture plane; ifR contains both negative and positive values, the fracture plane cuts through the element block (EB3 in Fig. 9 (d)). In the case shown in Fig. 9 (d), the element block EB3 is further divided into two element blocks by the fracture plane and thus the attributes in the Element block class should be updated accordingly.

(b) Algorithm

Since the element block is made up of element faces, the cutting operations of fracture plane with respect to the element block can be transformed to the cutting operations of fracture plane with respect to the element faces. The algorithm of fracture cutting the element block can be summarized as follows:

• Determine the existing element blocks that will be cut through by the fracture plane by comparing the bounding box of the element blocks and that of the fracture plane.

• Secondly, set two temporary variables, UpperElementBlock and LowerElementBlock, for each fracture plane. Employ the fracture plane to cut the element face of each existing element block identified in the first step and analyze the topological relationship between the resultant faces and the fracture plane: if the vertex series of the resultant element face is above the fracture plane, the element face is stored in the UpperElementBlock; while if the vertex series of the resultant element face is below the fracture plane, the element face is stored in the LowerElementBlock. After all the element faces of the existing block have been cut by the fracture plane, the resultant element blocks can be determined from the UpperElementBlock and LowerElementBlock.

• Note that when the fracture plane cuts through each element face, new vertexes will be generated. These newly-generated vertexes are numbered continuously following the sequence of existing vertexes. Tolerance control is performed to merge the overlapping vertexes.

• The boundary vertex list of the element face will be updated obeying an anticlockwise order. To achieve this, as shown in Fig. 10 the global vertex coordinates are firstly transformed to the local coordinate system of the fracture plane. Then, the anticlockwise order can be obtained by sorting the vertex in an ascending order of the angle between the line connecting the vertex and the center of the fracture face with respect to thex axis of the local coordinate system.

The above procedure can be illustrated by a simple case shown in Fig. 11: an element block composed of five element faces (F1(V₁-V₂-V₃), F2(V₄-V₅-V₆), F3(V₃-V₅-V₂-V₆), F4(V₁-V₂-V₅-V₄) and F5(V₁-V₃-V₆-V₄)), is cut by a fracture plane which passes through the vertex V₂. Fracture plane cutting is performed sequentially from F1 to F5. Assuming that the fracture plane is infinitely large, F1 will be divided into two sub-facesV₈-V₂-V₃ and V₁-V₇-V₈ where V₇ and V₈ are the two newly-created vertexes. V₈-V₂-V₃ is stored into the UpperElementBlock variable while V₁-V₇-V₈ is stored into the LowerElementBlock variable. F2 is also cut through by the fracture plane generating two new vertexes V₉ and V₁₀ , and two sub-faces V₉-V₅-V₆-V₁₀ and V₄-V₉-V₁₀. V₉-V₅-V₆-V₁₀ is stored into the UpperElementBlock set while V₄-V₉-V₁₀ is stored into the LowerElementBlock set. F3 has no intersection with the extended fracture plane and thus its vertex series V₃-V₅-V₂-V₆ is stored into the UpperBlockElement set as it is above the fracture plane. For F4, two new vertexes V₁₁ and V₁₂ and one UpperElementBlock V₁₁-V₅-V₁₂ and one LowerElementBlock V₁-V₂-V₁₂-V₄ are created. Similarly, two new vertexes V₁₃ and V₁₄ and one UpperElementBlock V₁₃-V₃-V₆-V₁₄ and one LowerElementBlock V₁-V₁₃-V₁₄-V₄ are generated for F5. Distance calculation identifies four groups of overlapping vertexes: (V₂, V₇, V₁₁), (V₈, V₁₃), (V₉, V₁₂) and (V₁₀, V₁₄). These overlapping vertexes are merged into four vertexes, i.e., V₂, V₇, V₈ and V₉, respectively and the vertex numbers are updated accordingly in the UpperElementBlock and LowerElementBlock as shown in Fig. 11.

Union-Find of element blocks data

After obtaining the element block data, the next step is to determine the assembly relationship between the element blocks when the finiteness of the fracture plane is restored.

As shown in Fig. 12, if two element blocks can be merged, the following conditions must be satisfied:

1) the two element blocks should be located in different half spaces of the fracture plane;

2) there should be overlapping between the element faces of the element blocks on the extended fracture plane (i.e., beyond the actual area of the fracture plane).

Based on the aforementioned criteria, the union-find algorithm of the element blocks is performed by the following steps (Fig. 13) :

• Firstly, the spatial relationship between the element faces of each element block and the fracture plane is sequentially obtained from the UpperElementBlock and LowerElementBlock variables determined in Section 3.1. For the element face which is partially within the fracture plane, if the element block is in the upper half space of the fracture plane, store the element face ID and the relevant element block ID into the UpperIndex of the fracture plane’s data structure (see Section 2.2); otherwise, they should be stored into the LowerIndex of the fracture’s data structure. For each fracture plane, this should be performed with respect to all the element blocks.

• Secondly, for each fracture plane, the separating axes method [] is used to determine whether the element faces within the UpperIndex data set are overlapping with the element faces within the LowerIndex data set. If so, group the ID of the two element blocks. Such operation must be performed for all the fracture planes. As Fig. 13 shows, two ID groups {EB1, EB8} and {EB3, EB4} were obtained for J1, no ID groups were obtained for J2, two ID groups {EB1, EB2} and {EB9, EB10} were identified for J3, one ID group {EB5, EB6} was identified for J4 and one ID group {EB4, EB5} was determined for J6.

• Finally, the union-find algorithm is employed to merge the ID groups with common element face IDs. The merged ID groups ({EB1, EB2, EB8}, {EB3, EB4, EB5 EB6} and {EB9, EB10}) are used in the later step to construct the complex blocks.

Construction of complex blocks

After obtaining the grouped ID of EBs that can be merged, the next step is to construct the final complex blocks.

(a) Basic concept

The construction of complex blocks follows the steps below:

• Firstly, use each EB ID group and relevant element block data to instantiate a complex block data array. The element face (EF) data with the same parent fracture ID are stored into the same temporary data array. Considering the EB group of {EB1, EB2} in Fig. 14, the data of EF7 (the front face for EB1) and EF1 (the front face for EB2) should be stored in the same temporary complex face data array CF1 as the parent fracture plane of both EF7 and EF1 is the fracture with ID 3. According to the same rule, other complex face (CF) data array can be obtained, i.e., CF2 (EF2), CF3 (EF3 and EF9), CF4 (EF4 and EF10), CF5 (EF5), CF6 (EF6 and EF11), CF7 (EF8) and CF8 (EF12).

• Secondly, using vertex chain operation to analyze each temporary complex face data array and obtain the final complex face data array as well as the final number of complex faces on the same fracture plane. As shown in Fig. 15, there are two typical relationships for element faces on the same fracture plane. The element faces with the same topo value should be combined, while for the element faces with the opposite topo values their common part needs to be subtracted.

• Thirdly, according to the previously identified complex face data, construct the enclosed space and analyze the topology of the enclosed space.

(b) polygon–subtraction-and-union method

Several unfavorable topological identification problems may occur in the second step of complex face identification. The identification results may be different using different bool operation algorithms due to their intrinsic accuracy or lack of robustness [–]. An inaccurate identification of the geometry of the complex face will lead to errors in the succeeding topological analysis of complex blocks. For example, for the vertex profile presented in Fig. 16, the loop identification may have three different results: (a) an outer loop with three inner loops, (b) two outer loops one of which contains two inner loops), or (c) three outer loops with one outer loop having an inner loops) due to different self-intersection identification algorithms. The Greiner and Hormann’s [] and Vatti’s [] algorithms usually lead to the former two cases, while obviously, the third case should be the correct one. In particular, when the number of polygon edges and the number of degenerate cases are increased, the Greiner and Hormann’s and Vatti’s algorithms are more likely to fail to obtain the correct results of polygon bool operations.

To avoid the misinterpretation of the geometry of complex faces, a new complex face identification algorithm is proposed, named as the polygon–subtraction-and-union method. This method is based on the idea of decomposing the polygons into several sub simpler polygons and then combine them. The polygon subtraction and union should obey the following principles:

• Polygons subtraction: the subtracted polygons must belong to different spaces (either upper or lower) of the same fracture plane. For example, EF1, EF2 and EF3 in Fig. 15 Case B should be subtracted.

• Polygons union: the added polygons must be located in the same half space (both upper or both lower) of the same fracture plane. For example, EF1, EF2 and EF3 in Fig. 15 Case A should be combined.

Polygons subtraction

As shown in Fig. 17, operations of polygon subtraction involve the following steps:

• Firstly, define the polygon having more vertexes as the subject polygon (the rectangle) and the polygon with fewer vertexes (the triangle) as the clipper polygon;

• Secondly, employ each edge of the two polygons to cut the two polygons into small sub-element faces ({EF1 EF3 EF4} for the subject polygon and {EF2 EF5 EF6} for the clipper polygon);

• Finally, remove the overlapping part EF3 and EF5.

Polygons union

After polygon subtraction, all polygons have been simplified as sub element faces and the overlapping element faces have been removed. Then polygon union can be carried out based on the vertex lists.

In the polygon union operations, if the vertex lists of the polygons satisfy any of the following three situations, the vertex lists of the two polygons could be merged.

Case 1: There are two vertexes in common for the two polygons and these two vertexes are successive in the vertex list (V₃-V₂ in Fig. 18 (a)).

For this case, taking Fig. 18 (a) as an example, the polygon union operations can be summarized as follows:

• Firstly, according to the sequential vertex order of the subject polygon (V₁-V₂-V₃-V₄-V₅-V₆), rearrange the vertex list order of the clipper polygon (V₂-V₃-V₇) with the vertex list starting from the first common vertex in the vertex list of the subject polygon (i.e., V₂) and ending at the second common vertex in the vertex list of the subject polygon (i.e.,V₃). Therefore, the vertex list of the clipper polygon after such a rearrangement becomes V₂-V₇-V₃.

• Then, break the vertex list of the subject polygon at the common vertexes, creating two vertex lists (i.e.,V₁-V₂ and V₃-V₄-V₅-V₆). Use the rearranged vertex list of the clipper polygon to chain these two vertex lists which leads toV₁-V₂- V₂- V₇-V₃- V₃-V₄-V₅-V₆. Remove the repeated common vertexes V₂ and V₃, and get the final vertex list for the complex face V₁-V₂-V₇-V₃-V₄-V₅-V₆. Note that after such an operation, the vertex list is still in an anticlockwise order, indicating that the loop is an outer boundary.

Case 2: There are more than two vertexes in common for the two polygons and all the vertexes are sequentially numbered in the vertexes list.

Polygon union operation for this case can be illustrated by Fig. 18 (b) in which the subject polygon (V₁-V₂-V₃-V₄-V₅-V₆-V₇-V₈) has three vertexes (V₂, V₃ and V₄) in common with the clipper polygon (V₂-V₃-V₄-V₉):

• Firstly, according to the sequential order of the subject polygon, remove the middle common vertexes (V₃ in Fig. 18 (b)) in the vertex list of the clipper polygon but retain the first (V₂) and the last common vertex (V₄). After such operation, the vertex list of the clipper polygon becomes V₂-V₄-V₉ and the problem is reduced to Case 1.

• Secondly, following the same procedure introduced for Case 1, combine the vertex lists of the two polygons and get the vertex list of the combined face, which isV₁-V₂-V₉-V₄-V₅-V₆-V₇-V₈ for the case shown in Fig. 18 (b).

Case 3: There are over two vertexes in common that are not in the sequential order for the two polygons (V₃-V₄, V₆-V₇, V₉-V₁₀). If at least two vertexes of these common vertexes are successive, the two polygons can still be merged.

Apparently, for Case 1 and Case 2 there will only be one chained vertex list for the resultant complex face. This chained vertex list should always represent an outer boundary and a solid face. However, vertex operations for Case 3 is more complicated as more than one chained vertex lists may be obtained and empty space with different shapes may be identified within the face. For the case illustrated in Fig. 18 (c), the subject polygon (V₁-V₂-V₃-V₄-V₅-V₆-V₇-V₈-V₉-V₁₀-V₁₁-V₁₂-V₁₃-V₁₄) has six vertexes (V₃, V₄, V₆, V₇, V₉ and V₁₀) in common with the clipper polygon (V₃-V₁₅-V₁₆-V₁₀-V₉-V₇-V₆-V₄). V₃ and V₄, V₆, and V₇ and V₉ and V₁₀ are the three pairs of successive vertexes. To deal with this situation,

• firstly, based on the sequential order of the common vertexes, separate the vertex lists of both the subject polygon and the clipper polygon into several sub vertex lists. If there are more than two vertexes in sequence, remove the middle vertexes and only retain the first and last vertexes in the common vertex list. After such an operation, we haveV₁-V₂-V₃, V₄-V₅-V₆, V₇-V₈-V₉ and V₁₀-V₁₁-V₁₂-V₁₃-V₁₄ for the subject polygon and V₃-V₁₅-V₁₆-V₁₀, V₉-V₇, V₆-V₄ for the clipper polygon.

• secondly, chain the sub vertex lists to get the vertex lists of the complex face. Chaining of the vertex lists starts from the sub vertex lists of the subject polygon in an order of ascending vertex number. The sub vertexes from the clipper polygon and the subject polygon at the common vertexes are alternatively joined until there is no more remaining sub vertex list initiating from the last vertex of the joined vertex list. According to this rule, three vertex lists can be generated. The first joined vertex list starts fromV₁-V₂-V₃ of the subject polygon, becomes V₁-V₂-V₃-V₃-V₁₅-V₁₆-V₁₀ after joined with V₃-V₁₅-V₁₆-V₁₀ of the clipper polygon and finally becomes V₁-V₂-V₃-V₃-V₁₅-V₁₆-V₁₀-V₁₀-V₁₁-V₁₂-V₁₃-V₁₄ after joined with V₁₀-V₁₁-V₁₂-V₁₃-V₁₄ of the subject polygon. This vertex list becomes V₁- V₂-V₃-V₁₅-V₁₆-V₁₀-V₁₁-V₁₂-V₁₃-V₁₄ after removing the repeated common vertexes V₃ and V₁₀. The second chained vertex list starts from V₄-V₅-V₆ and becomes V₄-V₅-V₆-V₆-V₄ after joined with the sub vertex list V₆-V₄ of the clipper polygon. After removing the repeated common vertexes (V₄ and V₆), the second vertex list becomes V₄-V₅-V₆. Similarly, we can obtain the third vertex list, which is V₇-V₈-V₉. Note that the first joined vertex list is in an anticlockwise order, while the second and the third vertex lists are in a clockwise order. According to the rule defined for the Complexface Class in Section 2.5, the first joined vertex list will be stored in the OuterLists indicating an outer periphery, while the second and the third vertex lists should be stored in the InnerLists and represent two inner boundaries. The merged complex face is illustrated in Fig. 18 (c) which shows thatV₁-V₂-V₃-V₁₅-V₁₆-V₁₀-V₁₁-V₁₂-V₁₃-V₁₄ is indeed an outer circumference and V₄-V₅-V₆ and V₇-V₈-V₉ are two inner boundaries, indicating that the current algorithm is correct and robust.

Application of the polygon union algorithms is not restricted to the aforementioned scenarios in which both faces are solid. They are also applicable to circumstances when the faces that need to be merged contain empty inner space such as the case shown in Fig. 19 (a). Of particular interest is that when the common vertexes coincide with the vertexes of the inner loop of the subject polygon, the complex face becomes a solid face (Fig. 19 (b)). There may be more than one complex face on the same fracture and it is necessary to calculate the number of complex faces on the same parent fracture.

After obtaining the topology of the enclosed space, the topological characteristics of the complex block can be readily identified: if the complex block only has one enclosed space, the complex block should be either a simple convex block (Fig. 20 (a)) or a simple concave block (Fig. 20 (b)). If there is any vertex not in the enclosed space made up of the outer loops, it must contain a cavity (Fig. 20 (c)). For the enclosed space containing inner loops, if it has only one loop on the external profile face, it must be a pit (Fig. 20 (d)), otherwise it must be a torus (Fig. 20 (e)).

Excavation simulation

In the current approach, excavation simulation is performed after the complex block construction has been completed for the original model. Excavation operation includes the following steps:

• Generation of the excavation region: based on the geometry of the excavation space, the excavation region is instantiated to be a series of element blocks. Similar to the generation of computation domain, if the geometry of the excavation space is too complex, it could be divided into several convex regions. The excavation faces are instantiated to be in the fracture plane class.

• Identification of the influenced blocks: due to the complex topological feature of complex blocks, following Zhang and Lei (2014) and Wu et al. (2015) the spatial relationship between the complex blocks and the excavation planes is evaluated considering the spatial relationship between the element blocks (ElementBlock 1 and 2 in Fig. 21) and the excavation plane. Only those complex blocks that intersect with the excavation element blocks will be considered for the later analysis.

• Generation of the new element blocks: considering the influenced complex blocks and following the methodology described in Section 3.1, determine the new element block data array (ElementBlock 2, 3 and 4 in Fig. 21).

• Removal of excavated element blocks: since the newly-generated element blocks could only be above or below the excavation planes, if the centroid of a newly-generated element block is within the excavation domain, it could be removed from the overall list of element blocks (ElementBlock 4 in Fig. 21).

• Reconstruction of the complex blocks: based on the principles given in Sections 3.2 and 3.3, reconstruct the final complex block data.

Validity checking

After the morphological analysis, it is necessary to check whether the topological analysis of the complex blocks is correct or not. There are generally two ways for correctness checking:

• Topological information check

Following Shi [], the Euler-Poincaré formula can be described as:

(3)N_v-N_e+N_f-(N_l-N_f)-\mathrm{2}( N_s-N_g)\mathrm{ =0},

where N_v, N_f, N_e, N_g, N_l, N_s are respectively the number of vertexes, faces, edges, holes, loops and shells.

• Block volume summation check

The total block volumes should be equal to the volume of the computation domain or the remaining volume of the computation domain after subtracting the excavation volume.

(4){\displaystyle \sum _{i=\mathrm{1}}^{n}V\left(B_i\right)}=V\left(T\right)\hspace{0.17em}\mathrm{or}\hspace{0.17em}{\displaystyle \sum _{i=\mathrm{1}}^{n}V\left(B_i\right)}=V\left(T\right)-V\left(E\right),

where n is the total number of complex blocks, V(B_i) is the volume of the ith complex block, V(T) is the total volume of the analysis domain and V(E) is the volume of the excavation area.

Application examples

To evaluate the effectiveness of the developed approach, two generic examples with different shapes of excavation are created.

Excavation of a cubic tunnel within jointed rock masses

As shown in Fig. 22. A cubic analysis domain with a dimension of 100 m×100 m×100 m is created and the geological information of the six boundary planes are listed in Table 1. Five fracture planes exist in the analysis domain and their geological information are given in Table 2. A tunnel with a rectangular cross-section of 30 m×30 m is directed in they direction through the analysis domain with the tunnel axis originated at the point (50, 0, 50).

Identifying the structure of rock mass after excavation contains two parts:

Stage I: Construction of block system prior to excavation

According to the procedure described in Section 3.4, the first step is to construct the block system prior to excavation. As Fig. 22 shows, 216 element blocks were created assuming that all fracture planes are infinite in size. Then union-operation was performed to cluster these element blocks by restoring the exact size of each fracture plane. Finally, block subtraction and block union were carried out to merge the element blocks in each cluster to form the final complex blocks. In total, twelve complex blocks were identified with the summed block volume being 999999.9985 m³ which is almost identical to the volume of the computational domain (1×10⁶m³).

Stage II: Construction of block system after excavation

After obtaining the block information before excavation, the four excavation planes defining the excavation domain were considered (Fig. 23). The spatial relationship between the complex blocks identified in Stage I and the excavation planes was evaluated by considering the spatial relationship between the element blocks comprising each complex block and the excavation planes. The complex blocks with ID 4–12 had intersection with the excavation planes, and the element blocks of these blocks were considered in the later analysis, while the element blocks of the complex blocks with ID 1–3 which were not affected by the excavation planes were not considered in the later union-find procedure. In total 351 element blocks were identified after introducing the four excavation planes, of which 57 were located within the excavation domain and were removed prior to the construction of complex blocks. The properties of each element block were updated in the Elementblock Class array, and they were assembled by block subtraction and block union to form the new complex blocks. The information of the previous defined complex block in the Complexblock Class array were updated accordingly.

The geometries of the final twelve complex blocks after excavation are illustrated in Fig. 23. Blocks 4, 5, 6, 8, 11 were concave, while the remaining blocks were convex. The total volume of these complex blocks were 909999.9928m³, which approximates the difference between the total volume of the analysis domain (1×10⁶m³) and the excavated volume (9×10⁴ m³). The geometry parameters of the final twelve complex blocks are given in Table 4.

As suggested by Jing [] and Shi [], the correctness of morphological analysis can be evaluated using the Euler-Poincaré formula. The numbers of vertexes, edges, faces, loops, shells and holes for each complex block are listed in Table 4. As can be calculated from Table 4, Eq. 4 is satisfied for each complex block. Considering the whole block system, there are in total 146 vertexes, 219 edges, 95 faces, 95 loops, 12 shell and 1 hole. Therefore, the Euler-Poincaré formula for the whole system becomes:

(5)\mathrm{146}-\mathrm{219}+\mathrm{95}-(\mathrm{95}-\mathrm{95})-\mathrm{2}(\mathrm{12}-\mathrm{1})=\mathrm{0}.

Excavation of a circular tunnel in jointed rock masses

In engineering practice, the geometry of an excavation can be quite complex. Topologically any shape can be treated as a composition of straight lines in 2D or planar faces in 3D. This rule is employed in the current study. Fig. 24 illustrates an example of identifying the topological feature of complex blocks after excavation of a circular tunnel in jointed rock mass. A cubic analysis domain with a dimension of 50 m×100 m×50 m is created and the geological information of the six boundary planes are listed in Table 5. Three fracture planes exist in the analysis domain and their geological information are given in Table 6. A circular tunnel with a radius of 7 m is directed in they direction through the analysis domain with the tunnel axis originated at the point (25, 0, 25). The procedure of constructing the complex block system is similar to Section 4.1 and is thereby omitted herein. To identify the complex blocks after excavation, the circular tunnel is replaced by a polyhedron composed of 15 planar planes with a stepwise change of 24° in dip angle (see Table 7). The closeness between a polyhedron and a circular tunnel can be improved by increasing the number of planar planes composing the polyhedral. However, by so doing, the computational cost will be increased accordingly as the number of times to evaluate the spatial relationship between element blocks and the excavation planes is proportional to the number of excavation planes. The total volume of the polyhedron is around 14947.5716 m³, which is very close to the total volume of the circular tunnel, i.e., 15393.804 m³. The spatial relationship between the complex blocks identified prior to excavation and the excavation planes was evaluated by considering the spatial relationship between the element blocks comprising each complex block and the excavation planes. In total 15 element blocks were identified after introducing the 15 excavation planes, from which eight complex blocks were identified after removing the element blocks within the excavation area. The geometrical parameters of each complex block after excavation are given in Table 8. The total volume of complex blocks after excavation is 235052.4283 m³, which is almost the same as the difference between the volume of the analysis domain (2.5×10⁵ m³) and the polyhedra. The Euler parameters for each complex block are also listed in Table 8. One can easily check that the Euler-Poincaré formula is satisfied for all the complex blocks.

Conclusions

An improved morphological visualization method is developed based on the block assembling approach to determine the morphological information of rock blocks. The new approach incorporates the concept of directed edge into the storage of vertex chain data. Different orientations are used to distinguish the internal loops from the external loops. Vertex chain operations are performed when assembling element blocks into complex blocks. The data structure of the new approach was presented and the algorithms for different stages of analysis were introduced. The proposed approach has following main merits in constructing complex rock blocks:

a) Application of the vertex chain operation in block assemblage analysis has its advantage as the vertex list of each element face/loop will be known after fracture cutting;

b) Block assembling is achieved by employing the polygon subtraction and union algorithms that are merely based on union operation of vertex chains with no need to calculate the intersection between complex faces;

c) The topology and convexity of the final complex blocks can be readily identified from the number and the type of their pertaining loops.

The proposed new algorithm has potential to be an efficient and convenient pre-processing tool prior to the kinematic analysis using block theory, and model initiation and contact detection in discrete element simulations.

References

[1] Lei Q, Wang  X. Tectonic interpretation of the scaling properties of a multiscale fracture system in limestone: structure and connectivity. Geophysical Research Letters, 2016, 43: 1551–1558

[2] Warburton P. A computer program for reconstructing blocky rock geometry and analyzing single block stability. Computers & Geosciences, 1985, 11(6): 707–712

[3]   Heliot D. Generating a blocky rock mass. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts: Elsevier; 1988, 25(3): 127–138

[4]   Lin D, Fairhurst  C, Starfield A . Geometrical identification of three-dimensional rock block systems using topological techniques. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts: Elsevier; 1987, 24(6): 331–338

[5] Ikegawa Y, Hudson  J A. A novel automatic identification system for three-dimensional multi-block systems. Engineering Computations, 1992, 9(2): 169–179

[6] Jing L. Block system construction for three-dimensional discrete element models of fractured rocks. International Journal of Rock Mechanics and Mining Sciences, 2000, 37(4): 645–659

[7] Lu J. Systematic identification of polyhedral rock blocks with arbitrary joints and faults. Computers and Geotechnics, 2002, 29(1): 49–72

[8] Elmouttie M, Krähenbühl  G, Poropat G . Robust algorithms for polyhedral modelling of fractured rock mass structure. Computers and Geotechnics, 2013, 53: 83–94

[9] Elmouttie M, Poropat  G, Krähenbühl G. Polyhedral modelling of rock mass structure. International Journal of Rock Mechanics and Mining Sciences, 2010, 47(4): 544–552

[10] Elmouttie M, Poropat  G, Krähenbühl G. Polyhedral modelling of underground excavations. Computers and Geotechnics, 2010, 37(4): 529–535

[11]   Cundall PA . Formulation of a three-dimensional distinct element model—Part I. A scheme to detect and represent contacts in a system composed of many polyhedral blocks. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts: Elsevier; 1988, 25(3): 107–116

[12]   Hart R, Cundall  P, Lemos J . Formulation of a three-dimensional distinct element model—Part II. Mechanical calculations for motion and interaction of a system composed of many polyhedral blocks. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts: Elsevier; 1988, 25(3): 117–125

[13] Yu Q C, Ohnishi  Y, Xue G F ,  Chen D. A generalized procedure to identify three‐dimensional rock blocks around complex excavations. International Journal for Numerical and Snalytical Methods in Geomechanics., 2009, 33(3): 355–375

[14] Zhang Y, Xiao  M, Chen J . A new methodology for block identification and its application in a large scale underground cavern complex. Tunnelling and Underground Space Technology, 2010, 25(2): 168–180

[15] Zhang Z, Lei  Q. Object-oriented modeling for three-dimensional multi-block systems. Computers and Geotechnics, 2013, 48: 208–227

[16] Zhang Z, Lei  Q. A morphological visualization method for removability analysis of blocks in discontinuous rock masses. Rock Mechanics and Rock Engineering, 2014, 47(4): 1237–1254

[17] Wu J, Zhang  Z, Kwok C . Stability analysis of rock blocks around a cross-harbor tunnel using the improved morphological visualization method. Engineering Geology, 2015, 187: 10–31

[18] Hao J, Shi  K B, Chen  G M, Bai  X J. Block theory of limited trace lengths and its application to probability analysis of block sliding of surrounding rock. Chinese Journal of Rock Mechanics and Engineering., 2014, 33(7): 1471–1477

[19] Liu X G, Zhu  H H, Liu  X Z, Wu  W. Improvement of contact detection algorithm of three-dimensional blocks. Chinese Journal of Rock Mechanics and Engineering., 2015, 34(3): 489–497

[20] Bourke P. Calculating the area and centroid of a polygon. 1988. Aavilable online:

[21] Büeler B, Enge  A, Fukuda K . Exact Volume Computation for Polytopes: A Practical Study. In: Kalai G, Ziegler G, editors. Polytopes — Combinatorics and Computation: Birkhäuser Basel; 2000, 29(6): 131–154

[22]   Greiner G ,  Hormann K . Efficient clipping of arbitrary polygons[J]. ACM Transactions on Graphics (TOG), 1998, 17(2): 71–83

[23] Vatti B R. A generic solution to polygon clipping. Communications of the ACM, 1992, 35(7): 56–63

[24] Rivero M, Feito  F R. Boolean operations on general planar polygons. Computers & Graphics, 2000, 24(6): 881–896

[25] Martínez F, Rueda  A J, Feito  F R. A new algorithm for computing Boolean operations on polygons. Computers & Geosciences, 2009, 35(6): 1177–1185

[26]   Shi G.Producing joint polygons, cutting joint blocks and finding key blocks for general free surfaces. Chinese Journal of Rock Mechanics and Engineering. 2006, 25(11): 2161–2170

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg