## 2 GPGPU-Parallelized 3D Y-HFDEM IDE  
In this study, an FDEM implemented in 3D Y-HFDEM IDE code was applied to simulate the fracturing of a direct shear test of a rough concrete-rock joint under shearing. Originally, the code was developed using the object-oriented programming tool (i.e., Visual C++) [44; 45; 46], and it was based on non-parallelized open-source Y-2D/3D libraries with OpenGL [37; 38; 36]. Using this code, the various rock engineering problems have been successfully solved by simulating rock fracturing [44; 45; 46]. However, due to the non-parallelized nature of the code, its previous applications were limited to small-scale 2D problems using relatively rough meshes. Recently, to overcome these limitations, a parallel programming scheme using GPGPU controlled by compute unified device architecture (CUDA) C/C++ was incorporated to develop the GPGPU-parallelized Y-HFDEM IDE code by the authors [47; 48; 49; 50; 51] for 2D and 3D FDEM modeling. Since the full details of the Y-HFDEM IDE code can be found in the literature [49; 50], we focus in this section on the essential parts of the 3D FDEM in order to model the fracture process in the direct shearing of the rough concrete-rock joint. In the remaining sections of this paper, the 3D FDEM implemented in the GPGPU-parallelized Y-HFDEM 3D IDE code is simply called the 3D FDEM. In the 3D FDEM, the continuous, i.e., elastic, deformations of rock and concrete are modeled using continuum solid elements (i.e., 4-noded tetrahedral elements, TET4s) [50]. The fracturing in the rock and concrete are modeled by an intrinsic cohesive zone model (ICZM) with the concept of smeared crack [55] representing opening and sliding cracks as well as mixed-mode I-II fracturing. In this model, 6-node initially zero-thickness cohesive elements (CE6) are inserted into all the boundaries of the TET4s at the beginning of the FDEM simulation to model the fracture of the rock and concrete (Figure 2). The opening and sliding of the faces in each CE6 are governed by tensile and shear softening curves as a function of the crack opening (\(w_{o}\)) and sliding displacements (\(w_{s}\)) (Figure 3). In the tensile and shear softening curves, the normal and shear cohesive tractions (i.e., \(\sigma_{\text{coh}}\) and \(\tau_{\text{coh}}\), respectively) acting on each side of CE6 are calculated using Equations (1) and (2), respectively:  
\[\sigma_{coh}=\left\{ \begin{array}{ll}\frac{2w_{o}}{w_{o}^{p}}T_{s}&if\ w_{o}<0 \\ \left[\frac{2w_{o}}{w_{o}^{p}}-\left(\frac{w_{o}}{w_{o}^{p}}\right)^{2}\right] f(D)T_{s}&if\ 0\leq w_{o}\leq w_{o}^{p}\\ f(D)T_{s}&if\ w_{o}^{p}<0\end{array} \right. \tag{1}\]  
\[\tau_{coh}=\left\{ \begin{array}{ll}\left[\frac{2|w_{o}|}{w_{s}^{p}}-\left( \frac{|w_{s}|}{w_{s}^{p}}\right)^{2}\right](-\sigma_{coh}\tan(\phi)+f(D)c)&if\ 0\leq |w_{s}|\leq w_{s}^{p}\\ -\sigma_{coh}\tan(\phi)+f(D)c&if\ w_{s}^{p}<|w_{s}|\end{array} \right. \tag{2}\]  
where \(w_{o}^{p}\) and \(w_{s}^{p}\) are the artificial elastic limits of \(w_{o}\) and \(w_{s}\), respectively. \(w_{o}^{overlap}\) is the representative overlap controlling the artificial elastic behavior of CE6s when \(o\) is negative, i.e., overlapping. \(T_{s}\), \(c\), and \(\phi\) are the microscopic tensile strength, cohesion, and internal friction angle, respectively, of the CE6. The \(w_{o}^{p}\), \(w_{s}^{p}\), and \(w_{o}^{overlap}\) are defined as \(w_{o}^{p}=2 ht_{s}/P_{open}\), \(w_{s}^{p}=2 hc/P_{tun}\), and \(w_{o}^{overlap}=2 ht_{s}/P_{overlap}\), respectively, where \(h\) is the element length, and \(P_{open}\), \(P_{tun}\) and \(P_{overlap}\) are the artificial penalty terms of the CE6 for opening, sliding, and overlapping, respectively. Equation (2) corresponds to the Mohr-Coulomb (MC) shear strength model with tension cut-off. The function \(f(D)\) is the characteristic function representing the tensile and shear softening curves (Figure 3), and it changes according to a damage value \(D\) of CE6.