## 2 Methods
### Determination of **n**, _Crack Propagation_, and _Initiation_  
A local criterion was first proposed in [54] for determining **n**, in which **n** is assumed to be the first eigenvector of the total strain \(\widehat{\boldsymbol{\varepsilon}^{(e)}}\), which is determined by  
\[\widehat{\boldsymbol{\varepsilon}^{(e)}}=\left[\begin{array}{c}\widehat{ \boldsymbol{\varepsilon}}_{x}^{(e)}\\ \widehat{\boldsymbol{\varepsilon}}_{y}^{(e)}\\ \widetilde{\gamma}_{xy}^{(e)}\end{array}\right]=\widehat{\
abla\mathbf{N}^{( e)}}\left[\begin{array}{c}\mathbf{u}_{1}\\ \vdots\\ \mathbf{u}_{n}\end{array}\right], \tag{14}\]  
which is independent of \(\zeta_{n}\) and \(\zeta_{t}\). After solving for the eigenvalue and eigenvector, we obtain  
\[\left[\begin{array}{c}n_{x}^{(e)}\\ n_{y}^{(e)}\end{array}\right]=\left[\begin{array}{c}a\;/\;\sqrt{a^{2}+b^{2}} \\ b\;/\;\sqrt{a^{2}+b^{2}}\end{array}\right]\]  
where  
\[a =\frac{\widehat{\gamma}_{xy}^{(e)}}{2} \tag{15}\] \[b =\frac{\widehat{\boldsymbol{\varepsilon}}_{y}^{(e)}-\widehat{ \boldsymbol{\varepsilon}}_{x}^{(e)}+\sqrt{\left(\widehat{\boldsymbol{ \varepsilon}}_{x}^{(e)}-\widehat{\boldsymbol{\varepsilon}}_{y}^{(e)}\right)^{2 }+\left(\widehat{\gamma}_{xy}^{(e)}\right)^{2}}}{2}.\]  
In the framework of the CEM, the elements will experience cracking one after another. Crack propagation is always checked first; then, crack initiation is considered. This procedure is based on the idea that, when the extensions of the existing cracks cannot fully release the extra stress, new cracks will appear. This procedure was used to capture the dynamic crack propagation in [55] where the branching and intersection of cracks are successfully captured.  
First, an index \(\phi_{RK}^{(e)}\) is introduced as  
\[\phi_{RK}^{(e)}=\left[\begin{array}{c}n_{x}^{(e)}\cdot n_{x}^{(e)}\\ n_{y}^{(e)}\cdot n_{y}^{(e)}\\ 2\ n_{x}^{(e)}\cdot n_{y}^{(e)}\end{array}\right]^{T}\mathbf{C}^{(e)}\left[ \begin{array}{c}\tilde{\varepsilon}_{x}^{(e)}\\ \tilde{\varepsilon}_{y}^{(e)}\\ \tilde{\gamma}_{xy}^{(e)}\end{array}\right]-f_{t}^{(e)} \tag{16}\]  
Then, the computing procedure below is conducted.  
1. The uncracked domain is divided into two subdomains, (i) a propagation domain and (ii) a crack root domain, based on a simple rule: when an element shares at least one edge with a cracking element, it belongs to the propagation domain; otherwise, it is in the crack root domain. Clearly, in the beginning, the whole domain is the crack root domain.
2. Find \(\left\{\max\left\{\phi_{RK}^{(e)}\right\}\right\}\) in the propagation domain. If \(\left\{\max\left\{\phi_{RK}^{(e)}\right\}\right\}>0\), the element becomes a cracking element; then, the two subdomains will be updated, the N-R iteration will be run, and this step will be conducted again. If \(\left\{\max\left\{\phi_{RK}^{(e)}\right\}\right\}\leq 0\), go to the next step.
3. Find \(\left\{\max\left\{\phi_{RK}^{(e)}\right\}\right\}\) in the crack root domain. If \(\left\{\max\left\{\phi_{RK}^{(e)}\right\}\right\}>0\), the element becomes a cracking element; then, the two subdomains will be updated, the N-R iteration will be run and step 2 will be conducted again. If \(\left\{\max\left\{\phi_{RK}^{(e)}\right\}\right\}\leq 0\), this loading step is considered to converge. A detailed flowchart can be found in [56].