Numerical Integration

For higher-order elements like Q4, the B-matrix varies with position, making the stiffness integral impossible to evaluate analytically:

$$[K_e] = \int_{-1}^{1} \int_{-1}^{1} [B]^T [D] [B] \, |J| \, d\xi \, d\eta$$

We need numerical integration — specifically Gauss quadrature — to evaluate these integrals efficiently and accurately.

Why Numerical Integration?

The Problem

For a Q4 element, the integrand contains:

• $[B]$ — varies linearly with $\xi$, $\eta$
• $|J|$ — varies with position for non-rectangular elements
• Their product is a polynomial that's tedious or impossible to integrate by hand

The Solution

Approximate the integral as a weighted sum:

$$\int_{-1}^{1} f(\xi) \, d\xi \approx \sum_{i=1}^{n} w_i \, f(\xi_i)$$

Where:

• $\xi_i$ = integration points (sampling locations)
• $w_i$ = weights (how much each point contributes)
• $n$ = number of integration points

The genius of Gauss quadrature: By choosing $\xi_i$ and $w_i$ optimally, we can integrate polynomials exactly up to degree $2n-1$.

Gauss-Legendre Quadrature

1-Point Rule

$$\int_{-1}^{1} f(\xi) \, d\xi \approx 2 \cdot f(0)$$

• Point: $\xi_1 = 0$
• Weight: $w_1 = 2$
• Exact for: polynomials up to degree 1 (linear)

2-Point Rule

$$\int_{-1}^{1} f(\xi) \, d\xi \approx 1 \cdot f(-\frac{1}{\sqrt{3}}) + 1 \cdot f(\frac{1}{\sqrt{3}})$$

• Points: $\xi_{1,2} = \pm \frac{1}{\sqrt{3}} \approx \pm 0.5774$
• Weights: $w_1 = w_2 = 1$
• Exact for: polynomials up to degree 3 (cubic)

3-Point Rule

• Points: $\xi_1 = 0$, $\xi_{2,3} = \pm\sqrt{0.6} \approx \pm 0.7746$
• Weights: $w_1 = \frac{8}{9}$, $w_{2,3} = \frac{5}{9}$
• Exact for: polynomials up to degree 5 (quintic)

General Pattern

2D Integration

For 2D integrals, apply 1D quadrature in each direction:

$$\int_{-1}^{1} \int_{-1}^{1} f(\xi, \eta) \, d\xi \, d\eta \approx \sum_{i=1}^{n} \sum_{j=1}^{m} w_i w_j \, f(\xi_i, \eta_j)$$

2×2 Gauss (Q4 Standard)

4 integration points at:

$$(\xi, \eta) = \left(\pm\frac{1}{\sqrt{3}}, \pm\frac{1}{\sqrt{3}}\right)$$

All weights = 1.

This integrates the Q4 stiffness matrix exactly for rectangular elements.

3×3 Gauss (Q8 Standard)

9 integration points — needed for quadratic elements.

Triangle Integration

Triangles use different rules based on area coordinates:

Full vs Reduced Integration

Full Integration

Use enough Gauss points to integrate the stiffness matrix exactly:

• Q4: 2×2 = 4 points
• Q8: 3×3 = 9 points
• CST: 1 point (constant integrand)
• LST: 3 points

Reduced Integration

Use fewer points than "full":

• Q4 reduced: 1 point (at center)
• Q8 reduced: 2×2 = 4 points

• Computational savings: Fewer evaluations of $[B]^T[D][B]$
• Avoids shear locking: Full integration can make elements too stiff in bending
• Better accuracy in some cases (counterintuitive!)

The Locking Problem

For a beam in bending:

• True behavior: $\gamma_{xy} = 0$ (no shear strain)
• Q4 with full integration: Forces $\gamma_{xy} \neq 0$ (artificial stiffness)
• Q4 with reduced integration: Can represent $\gamma_{xy} = 0$ correctly

Hourglassing

The downside of reduced integration: hourglass modes.

With 1-point integration, the element can deform in a zero-energy "hourglass" pattern:

This deformation produces zero strain at the single integration point!

• Hourglass control (artificial stiffness)
• Selective reduced integration
• Use more elements
• Use higher-order elements

Selective Reduced Integration

A clever compromise: Use different integration for different terms.

For Q4 in bending:

• Full integration (2×2) for volumetric/normal terms
• Reduced integration (1 point) for shear terms

This avoids locking while maintaining stability.

Integration Requirements by Element

Computing the Stiffness Matrix

Here's the algorithm for Q4 with 2×2 Gauss:

Initialize K = zeros(8, 8)

# Gauss points and weights
points = [(-1/√3, -1/√3), (1/√3, -1/√3),
          (1/√3, 1/√3), (-1/√3, 1/√3)]
weights = [1, 1, 1, 1]

for each Gauss point (ξᵢ, ηᵢ) with weight wᵢ:
    # Compute shape function derivatives
    dN/dξ, dN/dη = shape_derivatives(ξᵢ, ηᵢ)

    # Compute Jacobian
    J = compute_jacobian(dN/dξ, dN/dη, node_coords)
    detJ = det(J)

    # Transform derivatives to physical coordinates
    dN/dx, dN/dy = J⁻¹ @ [dN/dξ, dN/dη]

    # Build B-matrix
    B = build_B_matrix(dN/dx, dN/dy)

    # Accumulate stiffness
    K += wᵢ * B.T @ D @ B * detJ * thickness

Stress Recovery

After solving for displacements, we compute stresses at integration points:

$$\{\sigma\}_i = [D][B(\xi_i, \eta_i)]\{u_e\}$$

• Stresses are most accurate there (superconvergent)
• Extrapolating to nodes can introduce errors

Use shape function interpolation in reverse to get nodal stress values from Gauss point values.

Common Mistakes

• Wrong number of points: Using 1-point for Q4 without hourglass control
• Negative Jacobian: Check element orientation and shape quality
• Mismatched integration: Different rules for stiffness vs mass matrix
• Over-integration: Too many points wastes computation without benefit
• Stress at nodes: Extrapolated values less accurate than Gauss point values

Key Takeaways

• Gauss quadrature approximates integrals as weighted sums at optimal points
• n points integrate polynomials up to degree 2n-1 exactly
• 2D integration: Apply 1D rule in each direction (2×2, 3×3, etc.)
• Full integration: Enough points for exact stiffness (2×2 for Q4)
• Reduced integration: Fewer points, avoids locking but risks hourglassing
• Selective integration: Best of both — full for volumetric, reduced for shear
• Stresses most accurate at integration points (superconvergent)

What's Next

We've used natural coordinates $(\xi, \eta)$ and the Jacobian throughout. The next lesson dives deep into isoparametric formulation — the unified framework that makes all this work elegantly for arbitrary element shapes.