# Integration of Polynomials Times Double Step Function in Quadrilateral Domains for XFEM Analysis
## 3 Materials and Methods
### Algorithm Description  
The purpose of the proposed algorithm is to deliver the expression for the equivalent polynomial function \(\hat{H}_{i}(\mathbf{x})\), in order to compute the integral in Equation (16) without splitting the integration domain. The usefulness of the proposed algorithm is presented by way of a generic example. Let us consider a polynomial \(p_{n}(\mathbf{x})\) to be integrated over a subdomain \(\hat{\Omega}_{A}\), generated by dividing a parallelogram \(\hat{\Omega}\) with two lines \(\vec{q}\) and \(\vec{r}\), as shown in Figure 4a. The problem is defined in the global coordinate system \(\mathbf{x}=(x,y)\). Applying Equation (16) after the equivalent polynomial function \(\hat{H}_{A}(\mathbf{x})\) determines:  
\[I=\int_{\hat{\Omega}_{A}}p_{n}(\mathbf{x})\,\mathrm{d}\hat{\Omega}=\int_{\hat{ \Omega}}\hat{H}_{A}(\mathbf{x})p_{n}(\mathbf{x})\,\mathrm{d}\hat{\Omega} \tag{22}\]  
Beforehand, the problem has to be mapped to a standard quadrature domain. Therefore, a change in variables from the \((x,y)\) coordinate system to the parent coordinate system \((\vec{x},\eta)\) is employed in order to compute the integral over a standard and regular domain, as illustrated in Figure 4b. By means of such a procedure, single parent geometry can be used to address various scenarios. Thus, parallelograms having any position and size in the global coordinate system can be mapped to the square parent geometry in the local coordinate system \((\vec{x},\eta)\in[-1,+1]\) illustrated in Figure 4b. The mathematical concept used in the proposed algorithm (as well as in the _DD_EQP_ library) for this purpose is _isoparametric mapping_, which is commonly employed in the FEM [9; 49]. Let \(P(\vec{x},\eta)\in\Omega\) be a generic point in the parent reference system, corresponding to the point \(\bar{P}(x,y)\in\hat{\Omega}\) in the global reference system. The mapping of \(P(\vec{x},\eta)\) onto \(\bar{P}(x,y)\) is described by:  
\[x=\sum_{i=1}^{v}N_{i}(\vec{x},\eta)x_{i} \tag{23a}\] \[y=\sum_{i=1}^{v}N_{i}(\vec{x},\eta)y_{i} \tag{23b}\]  
Figure 4: Isoparametric mapping of a quadrangular element. (**a**) Element configuration in the global coordinate system. (**b**) Element configuration in the parent coordinate system. where \(v\) is the number of nodes of the geometric element denoted by coordinates \((x_{i},y_{i})\) in the global reference system, and \(N_{i}(\xi,\eta)\) is the shape function in terms of local coordinates for the parent element \(i\)-th node [20]. Likewise, the discontinuities \(\bar{q}(\mathbf{x})\) and \(r(\mathbf{x})\) equations, defined in the global reference system, have to be mapped onto \(q(\xi)\) and \(r(\xi)\), defined in the parent coordinate system (Figure 4). This is achieved by:  
* Calculating the signed distances (\(D_{i}\)) in the global coordinate system between each discontinuity and each node of the integration domain;
* Writing the discontinuity coefficients (\(a\), \(b\) and \(c\)) in the parent coordinate system as a function of \(D_{i}\) by solving a linear equations system;
* Substituting the variables \(x\) and \(y\) in \(\bar{q}(\mathbf{x})\) and \(\bar{r}(\mathbf{x})\) by means of Equation (23), so that \(q(\xi)\) and \(r(\xi)\) are obtained in terms of the coefficients \(a^{\prime}\), \(b^{\prime}\) and \(c^{\prime}\) dependent on \(D_{i}\). For a 2D square parent element, the coefficients are:  
\[a^{\prime}=\frac{D_{2}-D_{1}}{2} \tag{24a}\] \[b^{\prime}=\frac{D_{3}-D_{1}}{2}\] (24b) \[c^{\prime}=\frac{D_{2}+D_{3}}{2} \tag{24c}\]  
After \(\bar{q}(\mathbf{x})\) and \(\bar{r}(\mathbf{x})\) are transformed into \(q(\xi)\) and \(r(\xi)\), the correct expression for the equivalent polynomial function \(\bar{H}_{i}(\mathbf{x})\) with respect to the parent domain coordinate system can be generated.