## 3 Integration of Equations
### Finite Element Method  
Let \(\xi(t)=[\mathbf{u}(t),\mathbf{v}(t),\phi(t)]^{T}\in\mathcal{V}\times[0,T]\) be a solution of our system, where \(\mathcal{V}=H^{1}\times H^{1}\times L^{2}\) is infinite dimensional, \(T\) is some end time, \(H^{1}\) is selected as the appropriate function space for the variables \(\mathbf{u}\) and \(\mathbf{v}\), and where \(L^{2}\) is the chosen function space for the damage field.1 If we view the system in (18)-(20) in the abstract form  
Footnote 1: The notation \(H^{1}\) is standard in the field of partial differential equations to denote the Sobolev space of Lebesgue integrable functions whose weak derivatives of order one are also integrable. The notation \(L^{2}\) denotes the space of square (Lebesgue) integrable functions. Although we have not explicitly indicated it, it is understood that the \(H^{1}\) spaces in question conform to the prescribed Dirichlet boundary conditions for the fields of interest. \[\mathcal{F}(\xi(t))=0, \tag{21}\]  
we can then concisely write the problem's weak form as: Find \(\xi(t)\in\mathcal{V}\times[0,T]\) such that:  
\[(\mathcal{F}(\xi(t)),\zeta):=\int_{\Omega}\mathcal{F}(\xi(t))\zeta\,d\Omega=0 \quad\forall\zeta\in\mathcal{V},\quad\forall t\in[0,T], \tag{22}\]  
where \(\Omega\) is the spatial domain over which the governing equations of the problem are defined.2 As is customary with the Ritz-Galerkin Method, we now introduce a finite dimensional space \(\mathcal{V}_{h}\subset\mathcal{V}\), and let \(\xi_{h}\in\mathcal{V}_{h}\) denote an element from this space. We let \(\Xi=\{[\mathbf{w}^{i},\mathbf{w}^{i},w^{i}]^{T}\}\) denote a basis for \(\mathcal{V}_{h}\) (i.e., \(\mathrm{span}(\Xi)=\mathcal{V}_{h}\)) such that  
Footnote 2: \(\Omega\) is the domain of the problem. \[\mathbf{w}^{i}\in\big{(}\mathcal{Q}^{2}(K)\big{)}^{3}\quad\text{and}\quad w^{ i}\in\mathcal{P}^{0}(K), \tag{23}\]  
where \(K\) is any element in the triangulation \(\mathcal{T}_{h}\) of the domain \(\Omega\), \(\mathcal{Q}^{2}(K)\) is the space of quadratic Lagrange polynomials with support \(K\) and that are globally continuous over the triangulation, and \(\mathcal{P}^{0}(K)\) is the space of piecewise constant functions over the triangulation. We can now reformulate the weak statement of the problem in the usual way: Find \(\xi_{h}(t)\in\mathcal{V}_{h}\times[0,T]\) such that:  
\[\big{(}\mathcal{F}(\xi_{h}(t)),\Xi^{i}\big{)}=0\quad\forall\Xi^{i}\in\Xi, \quad\forall t\in[0,T]. \tag{24}\]  
Substituting (18)-(20) into the above, selecting the following representation of our trial solutions,  
\[u_{h}=U^{j}(t)\mathbf{w}^{j}(x),\;\;\;v_{h}=V^{j}(t)\mathbf{w}^{j}(x),\;\;\; \phi_{h}=\Phi^{j}(t)w^{j}(x), \tag{25}\]  
and applying integration by parts, yields a matrix form of (24),  
\[ \begin{bmatrix}M&0&0\\ 0&M^{\rho}&0\\ 0&0&M^{0}\end{bmatrix}\begin{bmatrix}\dot{U}^{j}\\ \dot{V}^{j}\\ \dot{\Phi}^{j}\end{bmatrix} + \begin{bmatrix}0&-M&0\\ A^{E}&0&0\\ 0&0&0\end{bmatrix}\begin{bmatrix}U^{j}\\ V^{j}\\ \Phi^{j}\end{bmatrix} = \begin{bmatrix}0\\ F_{v}\\ F_{\phi}\end{bmatrix} , \tag{26}\]  
where the elements of the block-matrices on the left-hand side of (26) are defined as  
\[M_{ij}=(\mathbf{w}^{j},\mathbf{w}^{i}),\qquad M_{ij}^{\rho}=( \rho\mathbf{w}^{j},\mathbf{w}^{i}), \tag{27}\] \[M_{ij}^{0}=(w^{j},w^{i}),\qquad\;\;A_{ij}^{E}=\big{(}(1-\phi)^{2 }\mathsf{C}_{0}[\
abla\mathbf{w}^{j}],\
abla\mathbf{w}^{i}\big{)}\,, \tag{28}\]and where the elements of the right-hand side of (26) are  
\[F_{v} =\left(\mathbf{b}_{\kappa},\mathbf{w}^{i}\right)+\left(\mathbf{s}_{ \kappa},\mathbf{w}_{i}\right)_{\partial\mathcal{B}}, \tag{29}\] \[F_{\phi} =\left(\eta_{c}\left\langle\frac{\beta}{2}(1-\phi_{h})^{\beta-1} \mathsf{C}_{0}[\
abla\mathbf{u}_{h}]\cdot\
abla\mathbf{u}_{h}-G_{cr}\right\rangle,w^{i}\right). \tag{30}\]  
Here \(\mathbf{s}_{\kappa}\) is the traction force applied to the boundary of the domain in the reference configuration. The finite element scheme just defined approximates the damage variable field as a piecewise constant function.