Question 3 of 7: Gauss quadrature over a mapped rectangular domain

Question 3: Gauss quadrature over a mapped rectangular domain (20 marks)

Question text not reproduced: the examination questions are © Engineers and Geoscientists BC. Open the official past paper (linked at the top of this page) to read the question, then follow the worked solution below.

Given. The integrand and domain data are collected below; the four bilinear shape functions printed in the question fix the node ordering, since $N_1$ vanishes at $\xi=1$ and at $\eta=1$ and therefore belongs to the corner $(-1,-1)$, and so on around the element counter-clockwise.

Find. The value of $g$ by Gauss–Legendre quadrature, presented with the explicit Jacobian and in the vector–matrix form used in finite element codes; then an explanation of how that value compares with the exact 3768.

Approach. Interpolate the geometry with the given $N_i$ to obtain $x(\xi)$ and $y(\eta)$, form the Jacobian matrix and its determinant, count the polynomial degree of the mapped integrand to select the quadrature order, then evaluate the double sum in vector–matrix form and compare with the exact result.

1. Interpolate the geometry. With $x=\sum_{i=1}^{4}N_i x_i$ and $y=\sum_{i=1}^{4}N_i y_i$ and the nodal coordinates above, $$x(\xi,\eta)=\tfrac14\Big[(1-\xi)(1-\eta)(1)+(1+\xi)(1-\eta)(7)+(1+\xi)(1+\eta)(7)+(1-\xi)(1+\eta)(1)\Big],$$ which collapses because the two nodes on each vertical side share an $x$ value: $$x = 4+3\xi, \qquad y = 3+2\eta .$$ These are exactly the mid-point-plus-half-width forms $x=\tfrac{x_1+x_2}{2}+\tfrac{x_2-x_1}{2}\xi$ that the affine map of a rectangle must produce, and they confirm the node ordering read off the shape functions.
2. Form the Jacobian matrix. By definition $$[J]=\begin{bmatrix}\dfrac{\partial x}{\partial \xi} & \dfrac{\partial y}{\partial \xi}\\[8pt] \dfrac{\partial x}{\partial \eta} & \dfrac{\partial y}{\partial \eta}\end{bmatrix} =\begin{bmatrix} \sum \dfrac{\partial N_i}{\partial \xi}x_i & \sum \dfrac{\partial N_i}{\partial \xi}y_i\\[8pt] \sum \dfrac{\partial N_i}{\partial \eta}x_i & \sum \dfrac{\partial N_i}{\partial \eta}y_i \end{bmatrix}.$$ Differentiating $x=4+3\xi$ and $y=3+2\eta$, $$[J]=\begin{bmatrix}3 & 0\\ 0 & 2\end{bmatrix},\qquad \left|J\right| = (3)(2)-(0)(0)=6 .$$ The determinant is constant, which is the signature of an affine (rectangular or parallelogram) element and the reason the exercise closes so cleanly. Dimensionally it is the area ratio: the parent square has area 4 and the physical rectangle $6\times4=24$, and $24/4=6$ as required.
3. Transform the integral. The change of variables gives $$g=\int_{1}^{5}\!\!\int_{1}^{7} f(x,y)\,dx\,dy =\int_{-1}^{1}\!\!\int_{-1}^{1} f\big(x(\xi),y(\eta)\big)\,\left|J\right|\,d\xi\,d\eta,$$ with $$f\big(x(\xi),y(\eta)\big)=\left(4+3\xi\right)^{2}\Big[\left(4+3\xi\right)+\left(3+2\eta\right)\Big].$$
4. Select the quadrature order. Expanding, the mapped integrand is a polynomial of degree 3 in $\xi$ and degree 1 in $\eta$. An $n$-point Gauss–Legendre rule integrates a polynomial of degree $2n-1$ exactly, so $2n-1\ge 3$ requires $n\ge 2$ in $\xi$; $n=1$ would suffice in $\eta$, but a $2\times2$ rule is used in both directions as is standard practice for bilinear elements. $$\boxed{\text{a } 2\times2 \text{ rule: } \xi_i,\eta_j=\pm\tfrac{1}{\sqrt3}=\pm0.577350,\quad w_i=w_j=1}$$
5. Map the four Gauss stations. Substituting the station coordinates into $x=4+3\xi$ and $y=3+2\eta$,

   Gauss stations and integrand values

   $(\xi,\eta)$	$x$	$y$	$f=x^{2}(x+y)$
   $(-0.577350,\,-0.577350)$	2.267949	1.845299	21.15688
   $(-0.577350,\,+0.577350)$	2.267949	4.154701	33.03550
   $(+0.577350,\,-0.577350)$	5.732051	1.845299	248.96450
   $(+0.577350,\,+0.577350)$	5.732051	4.154701	324.84312
   Sum			628.00000

6. Cast the evaluation in vector–matrix form. Collecting the four integrand values in the matrix $[F]$ with $F_{ij}=f(\xi_i,\eta_j)$ and the weights in the vector $\{w\}$, the double sum is a quadratic form: $$g=\left|J\right|\sum_{i=1}^{2}\sum_{j=1}^{2}w_i\,w_j\,F_{ij} =\left|J\right|\,\{w\}^{\mathsf T}\,[F]\,\{w\},\qquad \{w\}=\begin{Bmatrix}1\\ 1\end{Bmatrix},\quad [F]=\begin{bmatrix}21.15688 & 33.03550\\ 248.96450 & 324.84312\end{bmatrix}.$$ This is exactly the structure a finite element code uses when it accumulates an element matrix over the integration loop, with $|J|$ playing the role of the differential volume.
7. Evaluate. Since all weights are unity the quadratic form is simply the sum of the entries of $[F]$: $$\{w\}^{\mathsf T}[F]\{w\}=21.15688+33.03550+248.96450+324.84312=628.00000,$$ $$g = \left|J\right|\times 628.00000 = 6\times 628.00000$$ $$\boxed{\,g = 3768.00\,}$$
8. Part (b) — compare with the exact value. The exact integral is obtained directly: $$\int_{1}^{7}x^{3}dx=\frac{7^{4}-1}{4}=600,\qquad \int_{1}^{7}x^{2}dx=\frac{7^{3}-1}{3}=114,$$ $$g_{\text{exact}}=600\int_{1}^{5}dy+114\int_{1}^{5}y\,dy=600(4)+114(12)=2400+1368=3768,$$ confirming the quoted value. The numerical and exact results are therefore identical, not merely close — the difference is zero to machine precision, not a small truncation error.
9. Explain why they coincide. Two facts combine. First, the physical element is a rectangle with sides parallel to the axes, so the isoparametric map is affine: $x$ depends on $\xi$ alone, $y$ on $\eta$ alone, and $|J|=6$ is constant and can be taken outside the sum. Second, with a constant Jacobian the mapped integrand remains a polynomial, of degree 3 in $\xi$ and 1 in $\eta$, and a two-point Gauss rule is exact for polynomials up to degree $2n-1=3$ in each variable. Gauss quadrature is not an approximation here; it reproduces the integral identically. Had the element been a general (distorted) quadrilateral, $|J|$ would be a function of $\xi$ and $\eta$ and the integrand would become a rational function, for which no finite Gauss rule is exact.
10. Show the contrast with a lower-order rule. To demonstrate that the agreement is earned rather than automatic, apply the one-point rule ($\xi=\eta=0$, $w=2$ in each direction), which is exact only to degree 1: $$g_{1\text{-pt}} = |J|\,(2)(2)\,f(4,3) = 6\times 4\times \left(4^{2}\right)(4+3)=6\times4\times112 = 2688,$$ an error of $\dfrac{2688-3768}{3768}=-28.66\%$. The one-point rule under-integrates badly because it cannot see the cubic growth of $x^{3}$ across the element, and the same deficiency in a stiffness integral is what produces rank deficiency and hourglass modes.