Gauss quadrature over a bilinearly mapped rectangle: 22-Mec-B10 Finite Element Analysis, December 2016, Q1

Question 1: Gauss quadrature over a bilinearly mapped rectangle (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. A single bilinear (four-node) element covers the whole rectangle, so the geometry is described by the four printed shape functions acting on the corner coordinates:

Given data — Question 1
Quantity	Value
Integrand	$f(x,y)=x^{2}(y^{3}+1)$
Domain in $x$	$2 \le x \le 6$ (width 4)
Domain in $y$	$1 \le y \le 7$ (height 6)
Geometric interpolation	bilinear, $N_i=\tfrac14(1\pm\xi)(1\pm\eta)$
Parent domain	$-1 \le \xi,\eta \le 1$
Quoted exact value	$g_{\text{exact}}=42016$

Find. The value of $g$ obtained by Gauss–Legendre quadrature on the parent square, together with the smallest rule that is sufficient, and a reasoned comparison against the exact value.

Figure 1.1 — The rectangular domain 2 ≤ x ≤ 6, 1 ≤ y ≤ 7 and its bilinear map onto the parent square. The four teal markers are the 2 × 2 Gauss stations at ξ, η = ±1/√3, shown in both domains.

Approach. Map the rectangle onto the parent square with the given bilinear shape functions, show that the resulting Jacobian is constant, pick the Gauss order from the polynomial degree of the mapped integrand, and sum the weighted station values.

Map the physical rectangle onto the parent square. Interpolating the corner coordinates with the given shape functions, $x=\sum N_i x_i$ and $y=\sum N_i y_i$. For the corner set $(x_i)=(2,6,2,6)$ and $(y_i)=(1,1,7,7)$ this collapses to the affine relations $$x=\frac{2+6}{2}+\frac{6-2}{2}\,\xi = 4+2\xi,\qquad y=\frac{1+7}{2}+\frac{7-1}{2}\,\eta = 4+3\eta$$ because the element edges are parallel to the global axes and no cross term survives.
Evaluate the Jacobian of the mapping. Differentiating the mapping, $$[J]=\begin{bmatrix}\dfrac{\partial x}{\partial \xi} & \dfrac{\partial y}{\partial \xi}\\[6pt] \dfrac{\partial x}{\partial \eta} & \dfrac{\partial y}{\partial \eta}\end{bmatrix}=\begin{bmatrix}2 & 0\\ 0 & 3\end{bmatrix}$$ so that $$\boxed{\,|J| = (2)(3) - (0)(0) = 6\ \text{(constant over the element)}\,}$$ The constancy matters: it is what allows the quadrature order to be chosen from the degree of the integrand alone, with no rational function introduced by the mapping.
Rewrite the integral over the parent domain. With $dx\, dy = |J|\, d\xi\, d\eta$, $$g=\int_{-1}^{1}\!\!\int_{-1}^{1} (4+2\xi)^{2}\big[(4+3\eta)^{3}+1\big]\,|J|\;d\xi\, d\eta$$ Expanding the bracketed factors, the mapped integrand is a polynomial of degree 2 in $\xi$ and degree 3 in $\eta$.
Select the quadrature order. An $n$-point Gauss–Legendre rule integrates a polynomial of degree up to $2n-1$ exactly. The requirement in each direction is therefore $$2n-1 \ge 2 \;\Rightarrow\; n_\xi = 2,\qquad 2n-1 \ge 3 \;\Rightarrow\; n_\eta = 2$$ so a $2\times 2$ rule — four sampling points — is both necessary and sufficient. Its abscissae and weights are $$\xi_i,\eta_j = \pm\frac{1}{\sqrt{3}} = \pm 0.5773503,\qquad W_i = W_j = 1$$
Locate the stations in physical coordinates. Substituting the abscissae into the mapping, $$x = 4 \pm 2(0.5773503) \Rightarrow x = 2.845299,\; 5.154701;\qquad y = 4 \pm 3(0.5773503) \Rightarrow y = 2.267949,\; 5.732051$$

Accumulate the weighted station values. The quadrature sum is $$g \approx \sum_{i=1}^{2}\sum_{j=1}^{2} W_i W_j\,|J|\, x_{i}^{2}\big(y_{j}^{3}+1\big)$$ Each contribution is $6\, x_i^{2}(y_j^{3}+1)$ because both weights are unity:

Gauss station contributions ($|J|=6$, $W_iW_j=1$)
$(\xi_i,\eta_j)$	$x_i$	$y_j$	$x_i^{2}$	$y_j^{3}+1$	$6\, x_i^{2}(y_j^{3}+1)$
$(-,-)$	2.845299	2.267949	8.095729	12.665409	615.2143
$(-,+)$	2.845299	5.732051	8.095729	189.334591	9196.8093
$(+,-)$	5.154701	2.267949	26.570938	12.665409	2019.1907
$(+,+)$	5.154701	5.732051	26.570938	189.334591	30184.7857

Adding the four contributions, $615.2143+9196.8093+2019.1907+30184.7857$, gives $$\boxed{\, g_{2\times 2} = 42016.0\,}$$

Part (b) asks for the comparison, and the striking feature is that there is no discrepancy at all to explain away: the four-point rule returns the exact value to every digit carried. Two independent conditions make this happen simultaneously. First, the mapping is affine, so $|J|$ is a constant rather than a function of $\xi$ and $\eta$; a general distorted quadrilateral would put a rational function under the integral sign and no polynomial rule could then be exact. Second, the mapped integrand is of degree 2 in $\xi$ and 3 in $\eta$, and the two-point rule is exact through degree $2n-1=3$ in each direction. The quadrature error therefore vanishes identically, not merely to plotting accuracy.

The contrast worth drawing for the marker is what a cheaper rule would have produced. A single-point (centroid) rule samples $f$ at $x=4$, $y=4$ with weight $W=2$ in each direction: $$g_{1\times 1} = (2)(2)(6)\,(4)^{2}\big[(4)^{3}+1\big] = 24\,960$$ an error of $-40.6\%$. The one-point rule is exact only through degree 1, and the integrand is cubic in $y$, so its inability to see the curvature of $y^{3}$ is precisely the missing 40 per cent. Raising the rule to $3\times 3$ would return $42016$ again, at more than twice the cost per element for no gain — the practical lesson behind the question is that the correct quadrature order is the lowest one that is exact for the integrand at hand, since anything lower is inaccurate and anything higher is wasted work.

Question 1 — final results
Quantity	Result
Mapping	$x = 4+2\xi$, $y = 4+3\eta$
Jacobian determinant	$\|J\| = 6$ (constant)
Required rule	$2\times 2$ Gauss–Legendre ($\pm 1/\sqrt{3}$, $W=1$)
Quadrature result (a)	$g = 42016.0$
Exact value	$g_{\text{exact}} = 42016$
Difference (b)	Zero — the rule is exact for this integrand
One-point rule, for contrast	$24\,960$, i.e. $-40.6\%$ error

Topic. Gauss–Legendre quadrature in the isoparametric finite element formulation, and the choice of integration order.

Key relations. $\int_{-1}^{1}\!\!\int_{-1}^{1} F(\xi,\eta)\, d\xi\, d\eta \approx \sum_i\sum_j W_iW_j F(\xi_i,\eta_j)$; exactness of an $n$-point rule for polynomials of degree $\le 2n-1$; $dx\, dy = |J|\, d\xi\, d\eta$ with $|J|=\det[\partial(x,y)/\partial(\xi,\eta)]$; two-point abscissae $\pm 1/\sqrt3$ with unit weights.

Why this works. Gauss–Legendre quadrature is built by choosing both the abscissae and the weights to annihilate the error for the highest-degree polynomial the $2n$ free parameters can accommodate, which is degree $2n-1$. That is why a two-point rule beats a three-point Newton–Cotes (Simpson) rule despite using fewer evaluations. In finite element practice the integrand is a product of shape-function derivatives, material properties and the Jacobian, so the degree is known in advance and the rule can be pre-selected: $2\times 2$ for a four-node quadrilateral stiffness matrix, $3\times 3$ for an eight- or nine-node element. On an affinely mapped element the Jacobian leaves the parent-domain integrand polynomial, which is the situation exploited here; on a distorted element the Jacobian inverse makes the integrand rational and every rule becomes approximate, which is one reason badly shaped elements degrade accuracy.

Common pitfall. Choosing the order from the number of nodes rather than from the degree of the integrand — and, more seriously, taking the number of Gauss points needed for the geometry as the number needed for the integrand. A second frequent slip is forgetting the factor $|J|$ altogether, which here would report $g/6 = 7002.7$ and look like a gross modelling error rather than a missing scale factor.

Source. Logan, A First Course in the Finite Element Method, 6th ed., §10.4 (Gaussian quadrature, one- and two-dimensional rules); Reddy, 4th ed., §7.1.5; Bathe, Finite Element Procedures, 2nd ed., §5.5 (numerical integration and the required order).

Question 1 of 7: Gauss quadrature over a bilinearly mapped rectangle

Question 1: Gauss quadrature over a bilinearly mapped rectangle (20 marks)