with the Voigt form of the strain tensor. The test functions and their corresponding derivatives can be obtained from the above expressions by replacing $\bm{u}_I$ and $\phi_I$ with $\delta\bm{u}_I$ and $\delta\phi_I$ respectively. Due to the arbitrariness of the latter two quantities, numerical approximation of the weak form in yields the following nonlinear system of equations at each node $I$:

$$\begin{aligned}
        \bm{r}^{\bm{u}}_I &= \int_\Omega g \left( \phi \right) \bm{B}_I^{\bm{u}T} \bm{\sigma} {\,\mathrm{d}}\Omega - \int_\Omega \bm{N}_I^{\bm{u}T} \bm{b} {\,\mathrm{d}}\Omega - \int_{\partial\Omega} \bm{N}_I^{\bm{u}T} \bm{t} {\,\mathrm{d}}S = \bm{0} \label{eq:fem_linearMomentum} \\
        r^\phi_I &= \int_\Omega \left[ \mathcal{G}_c\ell\, \bm{B}_I^{\phi T} \nabla\phi + \frac{\mathcal{G}_c}{\ell} N_I \phi \right] {\,\mathrm{d}}\Omega + \int_\Omega N_I g^\prime \left( \phi \right) \mathcal{H} {\,\mathrm{d}}\Omega = 0 \label{eq:fem_phaseField}
    \end{aligned}$$

where $\mathcal{H}$ is the threshold-based history variable defined in . Due to the non-convexity of , the above coupled system is solved using the alternate minimization algorithm outlined in [@Bourdin2000]. This involves cycling between and : the linear momentum equation is solved first using values of $\phi$ from the previous iteration. Next, the phase-field evolution equation is solved using the newly obtained values of $\bm{u}$. The iterations are carried out repeatedly until the prescribed criteria on the size of residuals and inter-iteration corrections on the unknowns are met.

For 