# Integration of Polynomials Times Double Step Function in Quadrilateral Domains for XFEM Analysis
## 3 Materials and Methods
### Integration Algorithm for Single Discontinuity Problems  
Consider the body \(\mathcal{B}\) described in Section 3.1. Suppose that it is meshed utilising bilinear quadrilateral finite elements, and analyse an element \(\Omega\) of the mesh. Assuming that \(\Omega\) is split in two parts by a discontinuity line \(d\) (defined by (5)), it is possible to define two subdomains on the two sides of the discontinuity: \(\Omega^{+}\) and \(\Omega^{-}\) (Figure 1). \[d:ax+by+c=0 \tag{5}\]  
The vector \(\mathbf{n}^{+}\) in Figure 1 defines the positive portion of the element domain \(\Omega\) and points in the direction of \(\Omega^{+}\). The coefficients \(a\) and \(b\) in (5) for the internal discontinuity line \(d\) are used to describe the components of the vector \(\mathbf{n}^{+}\), which is orthogonal to \(d\). The equivalent polynomial technique [13; 19] allows us to map the integrand terms containing the discontinuous Heaviside function onto an equivalent polynomial with the same integral. This allows us to map the integrand onto a polynomial function, so that standard quadrature rules (e.g., Gaussian) can be applied and partitioning into subdomains is no longer required. To describe the equivalent polynomial technique, the symbolic form of the finite element stiffness matrix is taken as a reference [13], but the general idea is that any integrand containing terms given by a polynomial times the Heaviside function can be exactly integrated by replacing the discontinuous Heaviside function with an equivalent polynomial counterpart. Figure 1: A 2D quadrangular domain \(\Omega\) crossed by a discontinuity line \(d\). With reference to a classical linear small displacement elasticity problem, the displacement and strain fields in a single enriched finite element can be written as:  
\[\mathbf{u}=\mathbf{N}\mathbf{u}_{\Omega}+\Psi_{\Omega}\mathbf{N}\mathbf{a}_{\Omega} \tag{6}\]  
\[\epsilon=\mathbf{B}\mathbf{u}_{\Omega}+\Psi_{\Omega}\mathbf{B}\mathbf{a}_{\Omega}+(\
abla_{ \epsilon}\Psi_{\Omega})\mathbf{N}\mathbf{a}_{\Omega} \tag{7}\]  
where \(\mathbf{u}_{\Omega}\) and \(\mathbf{a}_{\Omega}\) are the element standard and enriched nodal variables and \(\
abla_{\epsilon}\) is the symmetric gradient operator, so \(\mathbf{B}\mathbf{u}_{\Omega}=(\
abla_{\epsilon}\mathbf{N})\mathbf{u}_{\Omega}\). \(\mathbf{E}\) is the elastic operator, so that the stress \(\mathbf{\sigma}\) is given by \(\mathbf{\sigma}=\mathbf{E}\mathbf{\epsilon}\)[13]. Therefore, we consider that:  
* The strain matrix \(\mathbf{B}\) is defined by linear functions for the bilinear quadrilateral element [13];
* \(\
abla_{\epsilon}\Psi_{\Omega}\) is identically zero [13]. Then, the element stiffness matrix for the nodal variables \(\mathbf{u}_{\Omega}\) and \(\mathbf{a}_{\Omega}\) is given by:  
\[K_{\Omega}=\int_{\Omega} \begin{bmatrix}\mathbf{k}_{uu}&\mathbf{k}_{ua}\\ \mathbf{k}_{au}&\mathbf{k}_{aa}\end{bmatrix} \mathrm{d}\Omega=\int_{\Omega} \begin{bmatrix} \mathbf{B}^{T}\mathbf{E}\mathbf{B}&\mathbf{H}\mathbf{B}^{T}\mathbf{E}\mathbf{B}\\ \mathbf{H}\mathbf{B}^{T}\mathbf{E}\mathbf{B}&\mathbf{H}^{2}\mathbf{B}^{T}\mathbf{E}\mathbf{B}\end{bmatrix}  \mathrm{d}\Omega \tag{8}\]  
where the off-diagonal part \(\mathbf{H}\mathbf{B}^{T}\mathbf{E}\mathbf{B}\) is not integrable by Gauss quadrature in the entire domain \(\Omega\)[13]. The solution proposed by [PERSON] [13; 19] through equivalent polynomials is based on substituting the discontinuous Heaviside function \(H(\mathbf{x})\) with an equivalent polynomial function \(\widehat{H}(\mathbf{x})\), continuous over the entire domain, \(\Omega\), such that [20]:  
\[\int_{\Omega}\widehat{H}(\mathbf{x})p_{n}(\mathbf{x})\,\mathrm{d}\Omega=\int_{\Omega ^{\top}}H(\mathbf{x})p_{n}(\mathbf{x})\,\mathrm{d}\Omega+\int_{\Omega^{+}}H(\mathbf{x})p_ {n}(\mathbf{x})\,\mathrm{d}\Omega \tag{9}\]  
where \(p_{n}(\mathbf{x})\) is an \(n\)-degree polynomial function to be integrated (such as \(\mathbf{B}^{T}\mathbf{E}\mathbf{B}\) in (8)).