As the complexity and challenges in the field of geomechanics rise, reducing computational costs has become a major theme to be investigated. This thesis evaluates the p-adaptive mesh optimization method’s performance for problems in the 3D finite element stress analysis of underground excavations with prismatic cross-sections. The p-adaptivity changes the element formulation within the finite element mesh, based on the concept of excavation disturbed zone. The changes in element formulation require the use of transition elements to connect linear hexahedral finite elements with quadratic ones. The forthcoming research was conducted using sim|FEM (Zsaki 2010), a research computer code intended for excavation design, which solves 2D plane strain and 3D problems. It was written in C++ and its key feature is its capability to test models with transition elements, which is not found in other analysis software. The research project started with an overview of 3D element formulations, both normal and transitional, which were implemented in the code then simple, yet practical models were tested and the results were analyzed. For some models, the results were compared to commercial software to prove that a correct behavior of the elements tested was obtained. Finally, the p-adaptive method was developed for this class of excavations and it was applied to a linear elastic medium with two circular excavations in a triaxial stress field representing a practical scenario of perhaps a transportation tunnel with a service tunnel running beside it excavated in intact rock. Mesh optimized and non-optimized models were compared and the results showed that optimization results in a reduction of the global stiffness matrix size on average by 82 percent and a reduction of solution times by about 82 percent for optimized models tested using 12-node and 16-node transitional elements, respectively, as compared to non-optimized models.