## 3 Integration of Equations
### Finite Difference Method  
The time-stepping algorithm has two stages. The first finds the values of the damage at the new time step, and the second finds the values of the displacement and velocity fields at the new time step. Specifically, denoting by \(n\) the time step, we have  
1. Solve explicitly for \(\Phi^{j,n+1}=\mathcal{F}_{1}(\Phi^{j,n};U^{j,n},V^{j,n})\) using the first order accurate forward Euler method, where \(\mathcal{F}_{1}\) is \(\mathcal{F}\) when \(\dot{u}\) and \(\dot{v}\) are set to zero.  
Specifically, the system to be discretized is:  
\[M^{0}\dot{\Phi}^{j}=F_{\phi}. \tag{31}\]  
Applying the forward Euler method, one arrives at the discrete formulation:  
\[M^{0}_{ij}\Phi^{j,n+1} =M^{0}_{ij}\Phi^{j,n}+\Delta tF^{n}_{\phi} \tag{32}\] \[=M^{0}_{ij}\Phi^{j,n}+\Delta t\left(\eta^{n}_{c}\bigg{/}\frac{ \beta}{2}(1-\phi^{n}_{h})^{\beta-1}\mathsf{C}[\
abla\mathbf{u}^{n}_{h}]\cdot \
abla\mathbf{u}^{n}_{h}-G^{n}_{cr}\bigg{)},w^{i}\right). \tag{33}\]  
Everything on the right-hand side is known, and thus we can solve for the new values of the damage.  
1. Solve implicitly for \([U^{j,n+1},V^{j,n+1}]=\mathcal{F}_{2}\left(U^{j,n+1},V^{j,n+1};U^{j,n},V^{j,n},\Phi^{j,n+1},\Phi^{j,n}\right)\), where \(\mathcal{F}_{2}\) is \(\mathcal{F}\) when \(\dot{\phi}\) is set to zero.  
Note that since we already have the values of \(\Phi^{j,n+1}\), the resulting system is linear. However, in such a system the matrices and right-hand-side are not constant, so that some additional care must be exercised when implementing this process into a computer code. The system we wish to solve is:  
\[\begin{bmatrix}M&0\\ 0&M^{\rho}\end{bmatrix}\begin{bmatrix}\dot{U}^{j}\\ \dot{V}^{j}\end{bmatrix}+\begin{bmatrix}0&-M\\ A^{E}&0\end{bmatrix}\begin{bmatrix}U^{j}\\ V^{j}\end{bmatrix}=\begin{bmatrix}0\\ F_{v}\end{bmatrix}, \tag{34}\]  
which can be more compactly written as:  
\[\mathcal{M}\dot{\xi}+\mathcal{K}\xi=\mathcal{F}, \tag{35}\]  
where \(\mathcal{M}\) is the mass matrix on the left-hand side of (34), \(\mathcal{K}\) is the stiffness matrix on the left-hand side of (34), and where \(\mathcal{F}=[0,F_{v}]^{T}\). Applying the Crank-Nicholson scheme (see, e.g., [13]) to the above yields the discrete form of (34),  
\[\left(\mathcal{M}+\frac{\Delta t}{2}\mathcal{K}^{n+1}\right)\xi^{n+1}=\left( \mathcal{M}-\frac{\Delta t}{2}\mathcal{K}^{n}\right)\xi^{n}+\frac{\Delta t}{ 2}\left(\mathcal{F}^{n+1}+\mathcal{F}^{n}\right). \tag{36}\]  
Note that the stiffness matrix (\(\mathcal{K}\)) is dependent on the time step only through the damage variable, which, due to the staggered solution algorithm, is now simply a parameter. The time-step is chosen based on the concepts of the Courant-Friedrichs-Lewy (CFL) limit (again, see [13]); we employ the wavespeed to determine the time-step via  
\[\Delta t=h_{\text{min}}\sqrt{\frac{\rho}{E}}, \tag{37}\]  
where \(h_{\text{min}}\) denotes the diameter of the smallest element in the grid.  
When the damage is constant, the choice of \(\mathcal{Q}^{2}\) finite elements and a Crank-Nicholson finite difference scheme results in a formulation which is second-order accurate in both time and space for the elastic unknowns; however, when damage is evolving we make no claim of the order of the accuracy of the overall algorithm. We have implemented the numerical methods presented in this section in a custom C++ program, which relies heavily on the deal.ii finite element library (see, for reference, [14]).