5  Exact vs. Gauss Integration on Reference Elements

This notebook explores numerical quadrature on reference elements. In FEM, we often use Gauss quadrature to approximate integrals over reference elements. The choice of quadrature order affects both accuracy and computational cost.

import sympy as sp
from IPython.display import Math, display
from symbolic_fem_workbench.quadrature import (
    integrate_reference_triangle_exact,
    integrate_reference_quadrilateral_exact,
    triangle_one_point_rule,
    triangle_three_point_rule,
    quadrilateral_gauss_rule,
)

5.1 Part A: Quadrilateral integration

We test Gauss quadrature on the reference square \([-1, 1]^2\) against exact symbolic integration. Different Gauss orders provide exact integration for different polynomial degrees.

xi, eta = sp.symbols("xi eta")

# A polynomial integrand
f1 = 1 + xi + eta + xi*eta
print("Integrand:", f1)

exact = integrate_reference_quadrilateral_exact(f1, xi, eta)
gauss1 = quadrilateral_gauss_rule(f1, xi, eta, order=1)
gauss2 = quadrilateral_gauss_rule(f1, xi, eta, order=2)

print(f"Exact integral:     {exact}")
print(f"Gauss order 1:      {gauss1}")
print(f"Gauss order 2:      {gauss2}")
print(f"Match (order 2)?    {sp.simplify(exact - gauss2) == 0}")

Integrand: eta*xi + eta + xi + 1
Exact integral:     4
Gauss order 1:      4
Gauss order 2:      4
Match (order 2)?    True

5.2 Part B: Triangle integration

The reference triangle uses different quadrature rules. We compare the one-point rule (centroid) and three-point rule (edge midpoints) against exact integration.

# Linear polynomial on triangle
f2 = 1 + 2*xi + 3*eta
exact_tri = integrate_reference_triangle_exact(f2, xi, eta)
one_pt = triangle_one_point_rule(f2, xi, eta)
three_pt = triangle_three_point_rule(f2, xi, eta)

print(f"Integrand: {f2}")
print(f"Exact:      {exact_tri}")
print(f"1-point:    {one_pt}")
print(f"3-point:    {three_pt}")
print(f"1-point matches exact: {sp.simplify(exact_tri - one_pt) == 0}")
print(f"3-point matches exact: {sp.simplify(exact_tri - three_pt) == 0}")

Integrand: 3*eta + 2*xi + 1
Exact:      4/3
1-point:    4/3
3-point:    4/3
1-point matches exact: True
3-point matches exact: True

5.3 Part C: When quadrature fails

Lower-order quadrature rules cannot exactly integrate higher-degree polynomials. This demonstrates why choosing the correct quadrature order is critical for element accuracy.

# Higher-degree polynomial
f3 = xi**2 * eta**2
exact_high = integrate_reference_quadrilateral_exact(f3, xi, eta)
gauss1_high = quadrilateral_gauss_rule(f3, xi, eta, order=1)
gauss2_high = quadrilateral_gauss_rule(f3, xi, eta, order=2)
gauss3_high = quadrilateral_gauss_rule(f3, xi, eta, order=3)

print(f"Integrand: {f3}")
print(f"Exact:      {exact_high}")
print(f"Gauss 1:    {gauss1_high}  (error: {sp.simplify(exact_high - gauss1_high)})")
print(f"Gauss 2:    {gauss2_high}  (error: {sp.simplify(exact_high - gauss2_high)})")
print(f"Gauss 3:    {gauss3_high}  (error: {sp.simplify(exact_high - gauss3_high)})")

Integrand: eta**2*xi**2
Exact:      4/9
Gauss 1:    0  (error: 4/9)
Gauss 2:    4/9  (error: 0)
Gauss 3:    4/9  (error: 0)

5.4 Part D: Element stiffness integration (teaching insight)

For a P1 triangle element with constant gradients, the stiffness integrand is proportional to \(\text{det}(J)\), which is constant. This means 1-point quadrature is sufficient—no higher-order quadrature is needed!

print("For P1 triangle stiffness: integrand ~ constant * det(J)")
print("Constant polynomial exact integration:")
f_const = sp.Integer(1)
exact_const = integrate_reference_triangle_exact(f_const, xi, eta)
one_pt_const = triangle_one_point_rule(f_const, xi, eta)
print(f"  Exact: {exact_const}")
print(f"  1-point: {one_pt_const}")
print(f"  Match: {exact_const == one_pt_const}")

For P1 triangle stiffness: integrand ~ constant * det(J)
Constant polynomial exact integration:
  Exact: 1/2
  1-point: 1/2
  Match: True

5.5 Exercise

Question: For a P1 triangle element, the stiffness integrand involves products of constant gradients times \(\text{det}(J)\). What is the minimum quadrature order needed for exact integration? Test your answer using the workbench.

Then consider a P2 element: the shape functions are quadratic, so gradients are linear. What quadrature order is needed then?

# Exercise: For a P1 triangle element, the stiffness integrand involves
# products of constant gradients times det(J). What is the minimum
# quadrature order needed for exact integration?
# Test your answer using the workbench.