## 2 Development of the GPGPU-Accelerated 3D FDEM Code for the RC Structure  
### GPGPU-Accelerated 3D FDEM for Modeling Concrete Fracturing  
A 3D FDEM code has been developed to model the fracture of rock-like materials such as rocks and concrete. This code was originally developed in a C++ programming language and is based on the open-source Y2D/Y3D libraries. This code has been successfully used to deal with a variety of rock engineering issues by modeling the fracture of rock and rock-like material [26; 27; 28]. However, since the applied algorithms in the original Y-code were sequential, its prior applications were restricted to small-scale 2D issues with relatively coarse meshes. To further improve computation processing such as speed and capabilities, the FDEM code has incorporated a parallel computing using a GPGPU device controlled by computing unified device architecture (CUDA) C++ [29; 30]. The FDEM is primarily based on continuum mechanics, the cohesive zone model (CZM), and contact mechanics within an explicit FEM framework. In this section, the only necessary component applied in this paper is addressed, while other comprehensive concepts of the FDEM code can be found in literature [29; 30]. In 3D FDEM, the continuous deformation including the stress wave propagation in the intact regime of materials such as concrete is modeled using a continuum 4-node tetrahedral finite element model (TET4s) (Figure 1) [29]. The CZM is applied with the concept of smeared crack to describe the transition from continuum to discontinuum [31]. The tensile and shear softening curves are applied using the assembly of 6-node cohesive elements(CE6) (Figure 1), which are initially zero-thickness, and inserted into all boundaries of the TET4s at the start of the simulation to model the behavior of the fracture process zone (FPZ) in front of the crack tip. Tensile and shear softening curves as a function of crack opening (o) and sliding displacements (s) control the opening and sliding of the faces in each CE6 (Figure 2). The normal and shear cohesive tractions (\(\sigma_{coh}\) and \(\tau_{coh}\)) acting on each side of CE6 are computed using Equations (1) and (2), respectively. \[\sigma=\left\{ \begin{array}{ll}\frac{2o}{\sigma_{overlap}}T_{s},&\text{if $o<0$}\\ \left[\frac{2o}{o_{p}}-\left(\frac{o}{o_{p}}\right)^{2}\right]f(D)T_{s},&\text{ if $0\leq o\leq o_{p}$}\\ f(D)T_{s},&\text{if $o_{p}<0$}\end{array} \right. \tag{1}\]  
\[\tau_{coh}\left\{ \begin{array}{ll}\left[\frac{2|s|}{s_{p}}-\left(\frac{|s|}{ s_{p}}\right)^{2}\right](-\sigma_{coh}\tan\phi+f(D)c),&\text{if $o\leq|s|\leq s_{p}$}\\ -\sigma_{coh}\tan\phi+f(D)c,&\text{if $s_{p}<|s|$}\end{array} \right. \tag{2}\]  
where \(o_{p}\) and \(s_{p}\) denote the corresponding artificial elastic limits of \(o\) and \(s\). \(o_{overlap}\) is the representative overlap, which controls the artificial elastic behavior of CE6s when \(o\) is negative, i.e., overlapping. \(T_{s}\), \(c\) and \(\phi\) are the microscopic tensile strength, the microscopic cohesion and the microscopic internal friction angle of the CE6, respectively. The \(o_{p}\), \(s_{p}\) and \(o_{overlap}\) are specified as \(o_{p}=2 hT_{s}/P_{overlap}\), \(s_{p}=2 hc/P_{tan}\) and \(o_{overlap}=2 hT_{s}/P_{overlap}\), respectively. \(h\) denotes the size of the CE6, and \(P_{open}\), \(P_{tan}\) and \(P_{overlap}\) indicate the CE6's artificial penalty factors for opening, sliding and overlapping. Equation (2) is equivalent to the Mohr-Coulomb (MC) shear strength model with a tension cut-off. The function \(f(D)\) is the characteristic function for the tensile and shear softening curves (Figure 2), and it varies with respect to the damage value D of CE6. D not only reflects the opening and sliding of the faces of CE6, but also mixed-mode crack behavior according to the following equation [31]:  
\[f(D)= \begin{array}{ll}\left[1-\frac{A+B-1}{A+B}\exp\left(D\frac{A+CB}{(A+B)( 1-A-B)}\right)\right]\\ \left[A(1-D)+B(1-D)^{C}\right]\end{array}