Finite element mesh improvement using an a priori local p-refinement for stress analysis of underground excavations

Abstract As our understanding and modeling capabilities evolve, an ever-increasing complexity of models representing the behaviour of geologic medium are analyzed. One way to evaluate these substantial problems is to optimize the underlying discretization of the governing differential equations by concentrating finite elements where solution accuracy counts the most. This paper develops and evaluates the performance of a priori local p-refinement method for finite element mesh improvement for stress analysis of underground excavations. This type of refinement entails a mesh with higher-order elements near the region of interest and lower-order elements elsewhere. The focus of the paper is the automated insertion of transitional elements at the interface of the two regions. The method relies on transitional finite elements in order to connect a mesh of quadratic interpolation order elements with a mesh of linear interpolation order elements. Four types of transitional elements were considered (4-node and 5-node triangles, 5-node and 7-node quadrilaterals). These were incorporated into a finite element code, and their performance was tested using representative problems such as a pressurized cavity or tunnelling through rock. For these problems the global stiffness matrix size was reduced on average by 85% and by 81% for the models using triangles and quadrilaterals, respectively, as a result, the calculation times were considerably shortened as well. While the average percentage of error with respect to the models without improvement, measured at critical points, was 0.04% and 0.02% in the case of triangular and quadrilateral elements, respectively.

Abstract: As our understanding and modeling capabilities evolve, an ever-increasing complexity of models representing the behaviour of geologic medium are analyzed. One way to evaluate these substantial problems is to optimize the underlying discretization of the governing differential equations by concentrating finite elements where solution accuracy counts the most. This paper develops and evaluates the performance of a priori local p-refinement method for finite element mesh improvement for stress analysis of underground excavations. This type of refinement entails a mesh with higher-order elements near the region of interest and lower-order elements elsewhere. The focus of the paper is the automated insertion of transitional elements at the interface of the two regions. The method relies on transitional finite elements in order to connect a mesh of quadratic interpolation order elements with a mesh of linear interpolation order elements. Four types of transitional elements were considered (4-node and 5-node triangles, 5-node and 7-node quadrilaterals). These were incorporated into a finite element code, and their performance was tested using representative problems such as a pressurized cavity or tunnelling through rock. For these problems the global stiffness matrix size was reduced on average by 85% and by 81% for the models using ABOUT THE AUTHORS Mr. D. Garcia Rosero completed his Master of Applied Science graduate degree under the supervision of Dr. Zsaki. His research interests are computer applications in geomechanics and finite element analysis. Dr. A. M. Zsaki is an Associate Professor in the Department of Building, Civil and Environmental Engineering at Concordia University. He obtained his B.Eng. degree from Ryerson University and his M.A.Sc. and Ph.D. degrees in civil engineering from the University of Toronto. Dr. Zsaki's research is focused on modeling and computational aspects of geosciences with particular interest in multiphysics modeling of continuum and discontinuum. His other areas of interest are scientific computing, parallel computing, computer graphics and mesh generation. In addition to academia, Dr. Zsaki has worked in the industry as software developer and consultant for geomechanics analysis software, and lately on high-performance scientific computing applications for modeling continuum behaviour.

PUBLIC INTEREST STATEMENT
As our understanding and modeling capabilities evolve, an ever-increasing complexity of models representing the behaviour of geologic medium are analyzed. This paper develops and evaluates the performance of a novel method for finite element mesh improvement for stress analysis of underground excavations. The focus of the paper is the automated insertion of elements that connect higher order (more precise) and lower-order (less precise) elements at the interface of these two regions. Their performance was tested using representative problems such as a pressurized cavity or tunnelling through rock. For these problems, the problem size was reduced on average by 85%, as a result, the calculation times were considerably shortened as well. While the average percentage of error with respect to the models without improvement, was on average 0.03%.

Introduction
Using the current state of modeling capabilities in geomechanics, one can create complex models incorporating millions of degrees of freedom. Non-linear and time-dependent material constitutive behavior and elaborate loading paths can be easily incorporated into the models as well. However, these models still require a considerable amount of resources both in computers to analyze and human expertise to generate the models, while perhaps unintentionally introducing user-bias during the modeling process. If a continuum is the most appropriate idealization of a problem, the finite element method (FEM) is one of the most used discretizations of the governing differential equations. One approach to reduce the consumption of resources is to adapt a model such that the density of nodes (mesh density) is increased locally. Another approach is to place more accurate (higher-order) elements where they count the most; in zones of rapid changes in stresses, dictated by estimation of solution error such as proposed by Zienkiewicz and Zhu (1987) and others (Kardani et al., 2013;Yu et al., 2018;Zienkiewicz et al., 2005), or in areas where the analyst requires high solution quality.
Traditionally, mesh adaptivity or refinement methods have increased the number of degrees of freedom in each consecutive solution of the problem to arrive at a better solution (Ikhenazen et al., 2019). This required repeated analysis of a problem at the expense of a longer solution time. There are several adaptivity-based methods to optimize meshes in the context of placing elements or nodes where they count the most (Sapre et al., 2019). The most notable methods are the radaptivity (relocating nodes), h-adaptivity (locally refining or coarsening mesh density), p-adaptivity (changing order of basis functions) and hp-adaptivity (a combination of the previous two). Their essentials are summarized in Figure 1.
Somewhat different approaches to mesh improvement is to use a mesh that is created, prior to an actual solution of the problem, such that it has more elements or higher-order elements in areas where such elements are needed the most. In this paper, the concept of improvement refers to generating a mesh that is appropriate to the problem being solved, e.g. by considering the underlying physics by using an approximate solution (Kirsch solution), to aid the placement of elements and nodes at the mesh generation phase. In contrast to traditional geometry-based approaches, where more elements are placed close to geometric features, such as excavations, our approach of mesh enrichment or improvement considers the physical definition of the problem to add more or better elements. This optimized mesh generation is done before any analysis is made, in contrast to the above methods. Once the mesh is optimized, the solution of a problem can be carried out, and there is no need to re-analyze a problem like in h-adaptivity methods. This problem-driven selective mesh refinement/coarsening method was developed by the second author (Zsaki & Curran, 2005a, 2005b, 2005c) and subsequently expanded on by Zsaki (2010) and Hazegh and Zsaki (2013). Instead of relying upon successive iterations of mesh refinement, the mesh improvement method incorporates an a priori knowledge contained in the definition and analysis intent of a problem being developed. This is accomplished by either knowing where highstress gradients are expected (around excavations, for example) or generating a high-quality (quality refers to mesh quality, such as aspect ratios or other metrics) mesh where the analyst needs it (opening new regions of a mine for ore extraction or proposed critical infrastructure). However, the previous mesh improvement method (Zsaki & Curran, 2005a, 2005b, 2005c and its variants (Hazegh & Zsaki, 2013;Zsaki, 2010) were based essentially on r-refinement and h-refinement; locally coarsening or refining mesh density along with relocating nodes. This paper presents an alternative approach to mesh improvement; a local p-refinement driven by the excavation disturbed zone (EDZ) concept done before any solution takes place, without a need for re-analysis. Starting with any finite element mesh populated with higher-order elements, irrespective if it is a converged mesh or not, the application of the proposed method is capable of reducing the solution time while preserving the accuracy of a model. This is achieved by selectively, using the EDZ concept, changing the element interpolation order. The development of this method will be explained in the following.
Excavations or openings in soil or a rock mass disturb the in-situ stress field around them. However, their effect is localized to about a few excavation diameters, beyond which the in-situ stress field is not affected (Brady & Brown, 2006;Martin & Kozak, 1992). The proposed method approximately maps out this zone of disturbance using an analytical solution given by Kirsch (in Ramamurthy, 2007) for infinite plates with circular holes. The analytical solution gives a continuum mechanics-based means to automatically guide the mesh improvement, which, in the authors' opinion, is considerably better than letting the person carrying out a mesh improvement to define less significant zones manually perhaps introducing user-bias. Although the analytical solution is valid only for an isotropic, homogeneous elastic medium with the excavations approximated by circles, its use to find the extent of disturbance caused by an excavation is deemed appropriate considering the inherent uncertainty and difficulty in obtaining representative soil or rock constitutive properties (Zsaki & Curran, 2005a). For anisotropic rock masses, solutions similar to the Kirsch solution can be adapted instead (Brady & Brown, 2006;Brown, 1987). In either case, this step is solely used to determine the zone of transition in the element order. Even though the EDZ, as determined by the Kirsch solution, is approximate, it is capable of estimating the extent of disturbance since the effect of an excavation diminishes at a rate of about one over the square of distance (cf. Kirsch solution) (Ramamurthy, 2007). The mesh improvement presented in this paper was developed for two-dimensional analysis using a plane strain assumption. But with appropriate analytical solutions for spherical or ellipsoidal cavities, it can be readily extended into three dimensions.
The extent of the EDZ is estimated to be where the presence of an excavation, or its approximation, disturbs the values of principal in-situ stresses by at least 5%, as measured by the percent difference between the principal stresses obtained with the excavation present and the preexcavation in-situ principal stresses (Zsaki & Curran, 2005a). Elements beyond this zone have their basis function's order reduced in comparison to the ones within the EDZ. For example, if second-order elements (isoparametric 6-noded triangular elements and 8-noded quadrilateral elements) are used in the EDZ, then first-order elements (isoparametric 3-noded triangular elements and 4-noded quadrilateral elements) are employed outside the EDZ. This necessitates the use of transitional elements to connect the second and first-order elements in a band within a model, as shown in Figure 2. Literature reports a few element formulations for transitional elements (Gupta, 1978) such as 4 and 5-noded isoparametric triangles and 5 and 7-noded isoparametric quadrilaterals. More recently, Huang and Xie (2010) have used a 5-noded transitional element for the solution of Poisson's equation. However, this paper uses the transitional element formulation developed by Gupta (1978). The practical validation and performance testing of these elements can be found in the publications of Gupta (1978) and Garcia Rosero (2011). For the sake of completeness, the shape or basis functions of 4-noded and 5-noded transitional triangular elements are summarized in Table 1 while Figure 3 shows the location of nodes and degrees of freedom for these elements. Figure 4 contains the same information for 5noded and 7-noded transitional quadrilateral elements and Table 2 lists their basis functions. These transitional elements were implemented in sim|FEM, which is a geomechanics-specific research code for finite element analysis for excavation design (Zsaki, 2018).

Procedure for mesh improvement using local p-refinement
The improvement of a mesh using the proposed procedure starts with an existing geometric model/mesh generated by any of the appropriate modeling tools available. For the purpose of the research summarized in this paper, Phase 2 (recently renamed to RS2) was used (rocscience inc, 2014) to generate the models, though Phase 2 was not used for the analysis. However, any current modeling tool could be used that generates finite element input files. The mesh should be entirely   comprised of the highest order elements that the modeler desires to incorporate; 6-noded triangles for the case of triangular elements and 8-noded quadrilaterals for quadrilateral elements. The input file should contain a complete definition of a model including the location of nodes, the definition of elements, boundary conditions and materials properties. Most of these will remain intact with the exception of nodal locations and the element definitions with respect to indices of constituent nodes.
In reference to the flowchart shown in Figure 5 and the rest of this section, the mesh improvement method was implemented into a standalone computer program named sim|OPTIM. This program takes an input file describing a model generated by Phase 2 and performs the mesh improvement. The output of sim|OPTIM is an input file that sim|FEM, our finite element package, can use. sim|FEM was used to carry out all finite element analysis in this paper. The first step in the mesh improvement process is to approximate each excavation by a circle. A circle is fitted such that its center is found by averaging the excavation geometry's coordinates in x and y directions and its radius is the distance computed from the center to the farthest vertex in the excavation geometry. The resulting circle slightly over-estimates the size of an excavation for non-circular shapes, nevertheless, resulting in a perhaps more accurate final numerical solution since the zone of higher-order elements extends farther. For excavation geometries with high aspect ratios (greater than 1.5 (Zsaki, 2003)), the circle-fitting could greatly over-estimate the original excavation size. However, since most tunneling works in civil and mining engineering projects have aspect ratios close to unity (e.g. a circle), a circular approximation and the subsequent analytical solution for stresses using the Kirsch equations are deemed appropriate (Zsaki & Curran, 2005a). For high aspect ratio excavations, the fundamental analytical solution for a circle should be substituted by one for an ellipse (Zsaki, 2003) without affecting the rest of the mesh improvement process. For the ensuing discussion, it is assumed that an excavation can be approximated by a circle. Once each excavation is approximated, the zone of influence (the extent of the EDZ) of each individual excavation is determined using the Kirsch equations. The outer boundary of the EDZ is taken where the excavation disturbs the in-situ stress field by no less than 5%, as measured by comparing the values of pre-excavation and post-excavation principal stresses. The Kirsch solution for computing the principal stresses is evaluated on a regular grid in the domain. Once the 5% level is found, a close-fitting circle is fitted to these grid points using standard computational geometry procedures. Although the 5% is quite arbitrary, it is a good compromise since it enables the definition of a zone far enough from an excavation yet puts a finite bound on its effect (Zsaki & Curran, 2005a).
Once the boundaries of the EDZ are mapped out, all elements intersecting this boundary are flagged to be potentially transformed into transitional elements. Since both 4-and 5-noded triangles and 5-and 7-noded quadrilaterals can be generated, the algorithm selects one of the generating schemes based on mesh type (triangles or quadrilaterals in the original mesh) and user preference (number of nodes representing a transitional element). All elements within the EDZ keep their original interpolation order and all elements outsize an EDZ are reduced in order (e.g. 3noded triangles or 4-noded quadrilaterals). At this point, the remaining band of flagged elements needs to be transformed into transitional elements. All elements that share an edge with a higherorder element (6-noded triangles or 8-noded quadrilaterals) retain their mid-side nodes on that edge and the mid-side nodes on their edge common with lower-order elements are removed, depending on the transitional element type. Thus, by iterating over all flagged elements, the mesh is transformed into separate zones of second, transitional and first interpolation order elements. These steps are summarized in the following pseudo-code as a two-pass algorithm over the set of flagged transitional elements: for all elements flagged as transitional if the element shares an edge with an element outsize the EDZ -reduce the element interpolation order to either 4 or 5-noded triangle or 5 or 7-noded quadrilateral, depending on mesh and user preference -change element type -remove mid-side nodes -flag element as processed for all elements flagged as transitional and not flagged as processed for all edges that connect to a flagged transitional element -remove, if needed a mid-side node on the connecting edge -based on mid-side node removal, change element type to an appropriate transitional element -flag element as processed However, not all flagged elements can be or need to be transformed into transitional elements. Consider the situation shown in Figure 6, where on the left a structured mesh of triangles is shown, while on the right an unstructured triangular mesh is depicted. The EDZ, represented by the dashed line, intersects the elements shaded in gray, denoting the potential transitional elements. As all potential transitional elements with mid-side nodes shared by the second-order element are retained, thus only those shown in darker gray on Figure 6 are actually transitional elements while the rest of the lighter gray elements are converted into first-order elements by removing all of their mid-side nodes. The procedure for meshes with quadrilateral elements is essentially the same; however, in our experience, it is very difficult to have all transitional elements within a mesh to be exclusively 5-noded or exclusively 7-noded transitional quadrilaterals. Usually, the mesh will contain a mixture of both types of transitional elements. Thus, in the first pass, 7-noded quadrilaterals are inserted beside 8-noded quadrilaterals. Any element adjacent to an edge with a midside node of a newly inserted 7-noded quadrilateral is transformed into a 5-noded quadrilateral thus completing the transition of a mesh.
The elimination of mid-side nodes necessitates an update of the data structure holding the model's nodes using a re-numbering phase to reduce data storage requirements. Similarly, the nodal indices for each element have to be updated to reflect a removal of mid-side nodes and the subsequent renumbering of nodes. Finally, all affected boundary conditions of a domain need to be updated; either removed if referring to an eliminated mid-side node or linked to the now re-numbered node. At this point, an updated input file is generated reflecting the presence of transitional elements and the associated changes in nodes and boundary conditions. Subsequently, this file can be used in sim|FEM analyses. For models comprised of multiple excavations, the effect of each excavation, as manifested by its EDZ, is assessed individually using the analytical solution. All elements in the model are flagged as either; a) inside an EDZ of an excavation, b) outside all EDZs or c) potential transitional elements. Once all EDZs are mapped out, a logical "union" operation is performed merging all intersecting EDZs and forming a zone encompassing all of them. Then, in turn, all flagged elements are either transformed into transitional or left as first-order elements (3-noded triangles or 4noded quadrilaterals) according to the rules described above.

Verification of practical performance of transitional elements-pressurized cavity
Once a mesh is modified to reflect the presence of transitional elements by using an updated input file generated by sim|OPTIM, the numerical solution of a problem can be performed using sim|FEM. As a simple assessment of the performance of this type of mesh, a model comprised of a pressurized cavity in an infinite elastic medium was created. This relatively simple model was used to test the capabilities of transitional elements. This model is representative of a typical water conveyance tunnel in a hydroelectric power generation system (Hoek, 2007). As a basis of comparison, the original models with 3-and 6-noded triangles and 4-and 8-noded quadrilaterals were solved in sim|FEM and results of the original models and improved model were compared. In addition, results from an h-refinement-based adaptive solution were compared to the ones obtained using the transitional elements.
The test model comprised of a domain containing a pressurized cavity. Only one-quarter of the domain was modeled by exploiting symmetry along two axes to reduce the model size. The loading and geometry of the problem are shown in Figure 7. The material representing the domain was assumed to be linear elastic and isotropic with a modulus of elasticity of 1.0 × 10 8 kPa, and Poisson's ratio of 0.3, characterizing a high compressive strength rock such as a competent igneous rock typical of the Canadian Shield (Jackson et al., 1995).
Three distinct zones within the model were considered for the sake of testing transitional elements in the model. Figure 7 shows the zones at 7.5, 10 and 12.5 m from the center of a 5m radius opening. Out of these three zones, the one at 10 m corresponds to the calculated EDZ, based the 5% influence, as set forth in the preceding section. The reason for considering the 7.5 Figure 7. Load, boundary conditions and geometry for the pressurized cavity modeltransitional zones also shown. and 12.5-m radius zones was to investigate the effect of moving the transitional zone closer and farther from the excavation and its impact on solution quality. For the comparison of performance, the domain was discretized using 3-noded and 6-noded triangles and 4-and 8-noded quadrilaterals. The 6-noded triangular and 8-noded quadrilateral meshes were then optimized using the presented algorithm for the cases of an inner zone at 7.5 m, EDZ at 10 m and outer zone at 12.5 m. For the case of triangular elements, transitional zones were represented by 5-noded triangles, while the transitional zones for quadrilateral elements were filled by 5-noded elements. These six models were then analyzed in sim|FEM and the results visualized in sim|FEM using its advanced PostScript output capabilities. Figure 8 shows the meshes with 3-noded triangles and 8-noded quadrilaterals.

Discussion of results and observations
For the case of meshes composed of triangular elements, vertical displacements along the left boundary were plotted, as shown in Figure 9. Two trends were observed; 3-noded triangles tend to under-estimate displacements in comparison to 6-noded triangles, while for the three cases containing transitional elements, the obtained values closely follow the trend established by the  6-noded triangles. For these three cases, there is no noticeable deviation in the values of displacement in the zones where 3-noded elements are used outside the EDZ. Although not shown, the same observations were made regarding displacements along the bottom boundary, where nearly identical values were obtained in the analysis. As measured by an L ∞ -norm, the maximum percent difference in displacements between the all 6-noded triangular mesh and meshes with transitional elements over the whole domain was found to be 2.34% (inner zone at 7.5 m), 2.26% (EDZ at 10 m) and 1.76% (outer zone at 12.5 m), respectively. In all three cases, the locations of nodes with the maximum value of the L ∞ -norm were around the excavation in 6-noded triangular elements. Similarly, for the models with transitional elements, the stress response of the whole domain was in close agreement to the values obtained for 6-noded triangles, as shown in Figure  10 (left image), where the contours of normal stresses in the x-direction (σ xx ) are shown for the model with the transitional zone located at 7.5 m (Figure 10, right image).
Assessment of results for models comprised of quadrilateral elements led to similar conclusions; 4-noded quadrilaterals under-estimate the displacements along the vertical (left) boundary in comparison to the second-order 8-noded quadrilaterals. Similar to the triangular transitional elements, meshes with 5 and 7-noded transitional quadrilaterals gave results close to those obtained using 8-noded quadrilaterals with the exception of one model, as shown in Figure 11. The model with the transitional zone at 12.5 m slightly (by 1.7% in an L ∞ -norm) under-estimates the displacements at the boundary of the circular excavation and over-estimates the displacements in the neighborhood of the transition zone (1.1%). These were the largest differences among any of the transitional meshes. In summary, similar to the triangular meshes, the L ∞ -norm for percent difference in displacements of the all 8-noded quadrilateral mesh and the meshes with quadrilateral transitional elements yielded the following; 2.06% (inner zone at 7.5 m), 1.89% (EDZ at 10 m) and 1.7% (outer zone at 12.5 m), respectively. In all cases, the nodes contributing to the highest L ∞ -norm were located near the circular excavation. Generally, the response of the domain was comparable to the base case (second-order 8-noded quadrilaterals, Figure 12, left image), as judged by what is shown in Figure 12 (right image) for the contour plot of normal stresses in the vertical direction (σ yy ) for the model with transitional elements at 12.5 m.
Another approach of assessing the performance of the proposed mesh improvement method is to compare it with other mesh improvement techniques. The most prevalent method of mesh Figure 10. Comparison of contours of horizontal normal stress σ xx (in kPa) for the models with 6-noded triangles (left) and the inner zone at 7.5 m (right).
improvement, generally a refinement, is based on h-refinement, as briefly discussed in the introduction. Although the proposed mesh improvement method does not need repeated re-meshing and re-computation of the solution like h-refinement methods, it is worthwhile to compare the results of using transitional elements in sim|FEM to that of results obtained from successive hrefinement. Figure 13    incorporated into sim|FEM. Element quality was improved by local Laplacian smoothing and edge swapping, based on the element aspect ratio, in the vicinity of newly-split elements.
The results of models that were refined using the h-refinement method and those with the transitional elements were plotted in a similar manner as before (cf. Figures 9 and 11) showing displacements along the left boundary of the model. Results for meshes with triangular elements are shown in Figure 14 and for meshes with quadrilateral elements in Figure 15, respectively. It is apparent from Figure 14 that the mesh with the first-level h-refinement (TRI-HREF1) underestimates the displacements, while the mesh with second level h-refinement (TRI-HREF2) Figure 13. Models used in the comparison of transitional (EDZ at 10 m) and h-refined meshes. Figure 14. Vertical displacements at left boundary, comparison of transitional and hrefined triangular meshes. and the transitional elements (TRI-TRANS) gave results that were nearly identical (0.4% difference at the boundary of the excavation). Thus, two key observations can be made at this point. First, even though, as it will be shown later in the discussion about performance, the TRI-HREF1 results in fewer degrees of freedom than TRI-TRANS, the h-refined mesh does not contain enough elements to give an accurate result. Second, since TRI-HREF2 and TRI-TRANS resulted in values close to each other, it can be concluded that TRI-TRANS is a converged mesh, given that TRI-HREF2 contains a very dense mesh in the vicinity of the circular excavation (about 10 times denser than the transitional mesh, as measured by element surface area). This can be further supported by looking at Figure 15, where the results for quadrilateral elements for the first (QUAD-HREF1) and second level of refinement (QUAD-HREF2) have converged to the same displacement value as TRI-TRANS. The results for the models with quadrilateral elements, as shown in Figure 15, gave values for all meshes, including the one with transitional quadrilateral elements (QUAD-TRANS), that were nearly the same (0.23% difference in displacements at the excavation boundary between QUAD-TRANS and QUAD-HREF1 and 0.46% between QUAD-TRANS and QUAD-HREF2).
With the confirmation of an acceptable performance of models with meshes incorporating transitional elements with respect to the accuracy of results, the savings in computational resources can be assessed. Two measures were used in determining this performance; reduction in the size of the global stiffness matrix (summarized in Table 3) and the associated reduction in the time to solution (including the mesh improvement steps outlined in Figure 5), the results of which were summarized in Table 4. The reduction of global stiffness matrix size is a direct result of the elimination of nodes outside the EDZ by the improvement algorithm, which effectively reduced the number of active degrees of freedom. The time to solution, although dependent on the processor/memory combination, is a good measure of practical performance. It considers factors such as paging (accessing data stored on hard drives that could not fit into memory) while accessing memory and the general performance of a computing platform. Table 3 lists the number of nodes, elements, active degrees of freedom and the size of the global stiffness matrix for all cases used in modeling the pressurized cavity problem.
For the case of triangles, memory consumption by the global stiffness matrix ranged from 7.5 MB for the 6-noded triangles to 0.5 MB for the 3-noded triangles. Similarly, the number of degrees of freedom was 972 on the high end and 252 on the low end. As evident from the previous discussion, the introduction of transitional elements along with the 6-noded and 3-noded triangles in the same mesh exhibits a behavior very close to that of a mesh comprised entirely of 6-noded triangles (cf. Figure 9). Reduction in the number of degrees of freedom achieved was about threefold for the first transition zone at 7.5 m, 2.4-fold for the second zone and 2.07-fold for the third transition zone at 12.5 m. Similarly, savings in memory consumption with respect to the 6noded triangle mesh were 8.83 times for the transition zone at 7.5 m, 5.76 times for the transition zone at 10 m and 4.31 times for the transition zone at 12.5 m. For the quadrilateral elements, reduction in degrees of freedom was 2.36 times, 1.97 times and 1.66 times for the three transition  zones. Although somewhat less than that for the triangular meshes, the achieved reduction was comparable. It has to be noted that this paper uses a simple full-matrix storage for the stiffness matrix, which is far from optimum when it comes to other memory storage schemes such as the "skyline" method. Nevertheless, no matter what the storage scheme is, the number of degrees of freedom was substantially reduced while giving accurate results.
The second measure used in the assessment of the performance of the improvement algorithm was the reduction in computation time. Table 4 summarizes findings for the optimized and unoptimized meshes along with the number of iteration steps the solution took for convergence using a conjugate gradient-type solver. For meshes comprised of triangular elements, the achieved reduction in computation time ranged from 13.23-fold for the transition zone at 7.5 m to about 4.86-fold for the furthest location of the transition zone. These reductions directly translate into a considerable speedup, as compared to the un-optimized mesh with 6-noded triangles. Similarly, for the quadrilateral elements, reduction in time ranged from 8.41 to 2.83 for the closest-toexcavation to the farthest transition zone as compared to the 8-noded quadrilateral mesh. Although not as high as for the triangular meshes, the speedup in solution time was substantial considering that the accuracy of the solution was maintained.
Meshes with transitional elements were compared for both accuracy and efficiency with meshes comprised of only linear and quadratic elements. In addition, their accuracy was compared to meshes obtained using h-refinement. But h-refinement comes at a cost; the model needs to be remeshed and re-analyzed for every refinement level. Table 5 summarizes the key parameters for meshes with transitional elements and those obtained using h-refinement. Even though the first level of h-refinement of the base mesh with linear triangular elements (TRI-HREF1) resulted in fewer degrees of freedom and faster solution time than the transitional mesh (TRI-TRANS), its solution accuracy was lower than the one obtained from TRI-TRANS, as discussed earlier. The model in the second level of refinement (TRI-HREF2) gave results nearly identical to TRI-TRANS, however, this was achieved at a premium; almost tripling the computation time and an almost 12fold increase in matrix storage. For meshes with quadrilateral elements, the solution of QUAD-TRANS takes about half the time to solve as compared with the model with the first level of refinement (QUAD-HREF1). Since QUAD-HREF1 matches the solution of QUAD-TRANS and that of the mesh with second-order 8-noded quadrilaterals, signaling that the first refinement level is sufficiently accurate, the comparison of QUAD-TRANS and QUAD-HREF2 is not necessary. Judging from the statistics in Table 5, the proposed mesh improvement method was able to solve the problem with fewer degrees of freedom than the h-refined meshes and it required less computation time.

Application of an a priori local p-refinement mesh improvement to underground excavations
As a practical example of the mesh improvement algorithm, consider a case of two parallel tunnels, perhaps representing a main tunnel and a service tunnel used as part of a transportation infrastructure. The tunnels are to be excavated out of a homogeneous, isotropic rock mass often found in the Canadian Shield (Jackson et al., 1995), with the properties and the prevailing state of in-situ stress summarized in Table 6. The rock constitutive model was chosen to be represented by a Mohr-Coulomb material, allowing both elastic deformation and yield, while being general enough. It was assumed that the tunnels are deep enough that the free surface effects can be neglected. The in-situ vertical stress of 7.5 MPa represents approximately 275 m of overburden rock, while the horizontal to vertical stress ratio of 1.33 is characteristic of many rock masses around the world (Hoek, 2007). Figure 16 shows the geometry of tunnels along with the extent of  the domain modeled. The problem was discretized using both 6-noded triangles and 8-noded quadrilaterals to serve as un-optimized base cases. Meshes, which are converged meshes with respect to their individual solution accuracy, for these two models are shown in Figures 17 and 18, respectively. The mesh with triangular elements was comprised of 4609 elements and 18582 degrees of freedom, while the mesh with quadrilateral elements was discretized into 4502 elements with 27210 degrees of freedom. The resulting contours for the magnitude of displacements (u T ) around the tunnels are shown in Figure 19 for both meshes.
Since there were two excavations, the improvement process for this model identified two zones of transition (the extent of EDZ, one for each excavation) using the analytical solution; one with a radius of 10.28 m centered slightly above the crown of the smaller excavation and for the larger excavation, the zone was found to extend 21.6 m and centered near the lower-left corner of that excavation. However, since there was an overlap between the two zones, a logical "union"  operation was performed on the two zones of transition, leading to a combined zone surrounding the two excavations. Figure 20 shows the zones for both mesh types along with the circles denoting the EDZ of each excavation. The elements intersected by the EDZ were flagged according to the rules set out in Section 2 and these elements were transformed into transitional elements. The resulting model input files were analyzed using sim|FEM and a summary of the number of nodes, elements, degrees of freedom, size of global stiffness matrix and solution times (which include the mesh improvement steps outlined in Figure 5, as appropriate) are listed in Table 7.

Discussion of results and observations
The reduction in the number of degrees of freedom for the optimized models was 2.65-fold for the triangular meshes and 2.24-fold for the quadrilateral meshes, respectively. This translated into a direct reduction in solution time by a factor of 12.5 and 6.18 for these two mesh types. The disparity in ratios between the number of degrees of freedom and solution time was attributed to the memory requirements; about 2.7GB was needed for the 6-noded triangles and 5.9GB for the 8noded quadrilaterals. Thus, for the model with quadrilaterals, disk swapping (due to lack of available memory) occurred when the memory was exhausted resulting in a considerable time penalty. For other, more memory-conscious matrix storage schemes, the achievable reduction in solution time should be less, in line with what was found for the mesh with triangles. Nevertheless, no matter what solution scheme is used, the number of degrees of freedom was reduced by a factor of 2.65 and 2.24, respectively.
With the advantages in shorter solution times established, the accuracy of the solution was investigated. Figure 21 shows the contours of the magnitude of displacement (u T ) for both optimized models. Note that the extent of the transition zone, where the element types change, is visible on the same figure. In comparison to Figure 19, the solutions show considerable Figure 19. Results of the unoptimized models (contours of the magnitude of displacements (in μmeters)).
agreement. Note that in the zones of transition the contour lines are smooth, thus the presence of transitional elements is not detrimental to the analysis results. For a numerical assessment of results, displacements at four points (top, bottom, right, left) on the boundary of each excavation were compared. Furthermore, the displacements along a line connecting the two excavations were compared considering that most analyses will be concerned with the behavior of a pillar between the excavations. Differences in displacements between the un-optimized (UM) and optimized (OM) models were calculated using the following equation: % difference at node i ¼ displacement in UM at node i À displacement in OM at node i displacement in UM at node i Tables 8 and 9 summarize the differences in the solution for the triangular and quadrilateral meshes, respectively. It is noted that for the triangular mesh containing transitional elements, the maximum difference (as compared to an all 6-noded triangular mesh) in any displacement (either in x and y direction) on the excavation boundaries was about 0.2%, while the maximum difference in the L ∞ -norm was about 0.28% over the whole mesh. For the mesh composed of quadrilateral elements, the maximum difference in displacements (between the all 8-noded quadrilateral mesh and one with transitional elements) around the excavation boundaries was about 0.15% and the maximum difference in L ∞ -norm of the displacement solution (in x and y direction) was about 0.19% for the whole mesh, which is still quite acceptable knowing the inherent uncertainty in parameters defining rock mass properties or in-situ stresses. Lastly, displacements between the two excavations were assessed, resulting in a plot shown in Figure 22 for the triangular meshes and Figure 23 for the one with quadrilaterals. The line spans the distance from the point denoted "right" on the left excavation to the point "left" on the right excavation, as already presented in Figure 16.
The maximum difference in displacements along the diagonal line was 0.23% for the case of models with triangular elements, while on average the displacements differed about 0.04%. Similarly, for the quadrilateral elements, the maximum difference was 0.04% and the average difference in results was 0.02%.

Conclusions
With a more thorough understanding of the mechanical behavior of soils and rock masses arises the challenge of analyzing models of ever-increasing complexity. However, the uncertainty in characterizing rock mass properties or in-situ stresses presents a limit in the accuracy achievable for any modeling in geomechanics. For large models analyzed using the FEM, an additional limit is posed by the availability of computational resources that could hinder modeling all the complexity inherent in a geologic medium. By reducing the size of a model, as measured by the number of degrees of freedom or the size of a global stiffness matrix, while maintaining the accuracy of the solution where it matters the most, analysts can incorporate considerably more detail in a model for a solution using the same computational resource. One way to reduce the number of degrees of freedom is to selectively eliminate nodes in a model based on the assumption that higher interpolation order elements, with more nodes, are required close to excavations or where they matter the most. The algorithm presented in this paper removes higher interpolation order elements that are outside the EDZ of an excavation while keeping them inside an EDZ. However, the transition from higher to lower interpolation order elements necessitates the use of a special element type to connect the two, ensuring the validity of finite element meshes. The use of transitional elements presented in this paper (4-and 5-noded triangles and 5-and 7-noded quadrilaterals) enabled the formation of this transition zone. Since the hybrid mesh, consisting of three types of elements, could pose problems in a FEM solution, the performance of these meshes was tested using a typical problem of a pressurized cavity expansion. Meshes comprised of  Figure 21. Results of the optimized models (contours of the magnitude of displacements (in μmeters)) with the transitional and EDZ zone overlaid.    triangular and quadrilateral elements were tested and the solution of the problem using hybrid meshes closely matched the solution obtained from meshes of 6-noded triangles and 8-noded quadrilaterals. Subsequently, the mesh improvement algorithm was used to reduce the model size for a typical problem in geomechanics, where the response of a geologic medium is sought after between two parallel excavations. The results obtained showed a considerable agreement between the un-optimized and optimized (or improved) models. Thus, using an a priori mesh modification, for the examples presented, it was possible to considerably reduce model sizes. It is anticipated that the method has the potential to reduce problem sizes in general, while preserving the accuracy of the solution where it counts the most.