## 1 Introduction  
Quasi-brittle materials such as concrete and fiber-reinforced concrete are most commonly used in engineering structures. Predicting the limit load capacity and failure behavior is of vital importance in preventing the potential risks of these structures due to the extremely stochastic and nonlinear properties of quasi-brittle materials [1, 2]. The accurate analysis of failure in mechanical response, crack initiation and propagation via computational modeling is not only an indispensable tool for which systematic experiments are expensive, time-cost, difficult or even impractical, but it also provides insights into understanding the failure processes and mechanical behaviors of these materials, such as in the third numerical example provided in this work [3]. For quasi-brittle materials with strain softening and intrinsically heterogeneous properties, their fracture process and mechanical behavior cannot be analyzed simply by elastic or plastic mechanics methods. In the past decades, various approaches have been proposed to understand and model the mechanical response, crack initiation and propagation of quasi-brittle materials under different conditions [4, 5, 6]. These approaches can be divided into two groups: damage mechanics methods and fracture mechanics methods. However, conventional damage mechanics methods often suffer from mesh bias problems due to the strain localization phenomenon [4, 5] unless some additional methods are proposed, as in the nonlocal and gradient-enhanced models [6, 7, 8]. Similar problems exist in the smeared/diffusive and discrete methods of fracture mechanics. In fracture mechanics, there are many methods currently proposed for crack propagation modeling. These methods may be classified into discrete and smeared/diffusive methods. Discrete methods model cracks as real discontinuities. Representative methods are cohesive zone models(CZMs), XFEMs, etc. However, the CZM cannot easily obtain mesh-independent results [9; 10; 11]. The XFEM requires an explicit representation of the crack and additional criteria to achieve complex crack growth, and its convergence is not robust, especially in a 3D setting [12; 13]. On the other hand, the diffusive approach can be seen as a development of continuum damage mechanics, which treat the cracks as the result of damage accumulation. As the most popular diffuse approach, the phase field (PF) method has recently received more and more attention from researchers for its advantage of not requiring additional criteria to predict crack nucleation and propagation [1; 2; 14]. The method started with the linear-elastic fracture variational principle proposed by [PERSON] and [PERSON] [15]  
\[\Psi(\mathbf{u},\Gamma)=\Psi^{b}+\Psi^{\Gamma}=\int_{\Omega}\psi(\epsilon,\phi) \mathrm{d}\Omega+G_{c}\int_{\Gamma}\mathrm{d}\Gamma \tag{1}\]  
where the total potential energy functional \(\Psi(\mathbf{u},\Gamma)\) is assumed to include the surface energy \(\Psi^{\Gamma}\) of the crack and the bulk energy \(\Psi^{b}\); more details will be discussed in Section 2. Subsequently, [PERSON] [16], [PERSON] [17], [PERSON] [18] and some other scholars made further extensions to the PF model, making the phase field modeling for failure generalized. PF methods have been demonstrated to be a powerful modeling tool due to their elegant format for the prediction of failure in many engineering materials such as rock [19], asphalt concrete [20], fiber-reinforced composites [21; 22], piezoelectric composites [23] and shield tunnel lining [24], to cite a few. In 2017, [PERSON] discovered the commonality between the classical PF model and the nonlocal damage model and proposed a unified phase-field theory that is suitable for both brittle and cohesive failure [1; 25]. The model (e.g., PF-CZM) derived from the above theory has been validated to be able to overcome most of the drawbacks mentioned in the above models, such as mesh bias and length scale sensitivity, and even the size effect can be captured accurately [26; 27]. See [28; 29] for recent progress.