# Integration of Polynomials Times Double Step Function in Quadrilateral Domains for XFEM Analysis
## 3 Materials and Methods
### Integration Algorithm for Single Discontinuity Problems
If such polynomial \(\widehat{H}(\mathbf{x})\) exists, the integration can be performed on the entire domain of the element using traditional Gauss quadrature [13]. Therefore, substituting the definition of \(H(x)\) from (4) into (9), a polynomial function \(p_{n}(\mathbf{x})\) can be integrated over the domain \(\Omega^{+}\) (Figure 1) through the integral:  
\[I=\int_{\Omega^{+}}p_{n}(\mathbf{x})\,\mathrm{d}\Omega=\int_{\Omega}H(\mathbf{x})p_{ n}(\mathbf{x})\,\mathrm{d}\Omega=\int_{\Omega}\widehat{H}(\mathbf{x})p_{n}(\mathbf{x})\, \mathrm{d}\Omega \tag{10}\]  
Thus, it is no longer required to partition the domain of integration, and the integral The element stiffness matrix (in the context of XFEM) has to be evaluated on each of the four subdomains \(\Omega_{A}\), \(\Omega_{B}\), \(\Omega_{C}\) or \(\Omega_{D}\). As seen in (8), such a matrix will contain discontinuous functions times polynomials that cannot be integrated using standard quadrature rules over the entire domain, \(\Omega\). It should be observed that the definition of the Heaviside step function (4) is such that the integrand function is zeroed on the subdomains with negative signed distance, so that a proper definition of the discontinuities normal vectors allows us to perform the quadrature on each subdomain. This also allows straightforward extension of the method to the generalised Heaviside function, having values \(+1\) and \(-1\) on the two sides of the discontinuity instead of \(+1\) and \(0\), as per definition (4). Assume an \(n\)th-degree polynomial \(p_{n}(\mathbf{x})\) to be integrated across the subdomains \(\Omega_{A}\), \(\Omega_{B}\), \(\Omega_{C}\) or \(\Omega_{D}\) obtained by partitioning a regular \(2\times 2\) square centred in \((0,0)\) on the \((x,y)\) reference system with two lines \(q\) and \(r\), so it is:  
\[I_{i}=\int_{\Omega_{i}}p_{n}(\mathbf{x})\,\mathrm{d}\Omega \tag{11}\]  
where: \(i=\{A,B,C,D\}\) and \(p_{n}(\mathbf{x})=k_{0}+k_{1}x+k_{2}x^{2}+k_{3}xy+k_{4}y+k_{5}y^{2}+...+k_{m-1}x^ {n}+k_{m}y^{n}\). To compute the integral in (11) for each value of \(i\) with standard quadrature approaches, the definition of each integration subdomain \(\Omega_{i}\) is required. Due to the possibility of the integration domains producing rather complex polygonal shapes depending on the slope of the lines, this is not always an easy task. The target is to use the equivalent polynomials for the two discontinuities \(p\) and \(r\) in order to allow, for each subdomain \(\Omega_{i}\), the integration over the entire element domain \(\Omega\) by means of the regular Gauss-Legendre quadrature rule, i.e., to compute (11) with integrations on the entire element domain \(\Omega\) instead of \(\Omega_{i}\). Let \(\hat{H}_{p^{+}}\) and \(\hat{H}_{r^{+}}\) be the equivalent polynomials related to the normals \(\mathbf{n}_{p^{+}}\) and \(\mathbf{n}_{r^{+}}\), and \(\hat{H}_{p^{-}}\) and \(\hat{H}_{r^{-}}\)  
Equations (12) to (15) give a system of four equations in the four unknowns \(I_{A},I_{B},I_{C}\) and \(I_{D}\) that, in general, can be proved to be indeterminate. In particular, it can be observed that if the intersection point between the two discontinuities \(q\) and \(r\) is external to the element domain or is on its boundary, then the solution will be unique. If the intersection point lies inside the element domain, the system will be indeterminate, with the system coefficient matrix having rank three. When the intersection point between the two discontinuities \(q\) and \(r\) is internal to the element domain, the above observation suggests the introduction of an auxiliary integration limit \(s\) along the abscissa axis of the element domain to eliminate indeterminacy. The line \(x=s\) contains the discontinuities intersection point \(P\) and is parallel to the vertical axis of the reference system. It should be noted that the introduction of the auxiliary integration limit \(s\) keeps the reduced integration domain rectangular, so that standard quadrature rules can be applied.