Computing Scalar Integrals Over Triangles- Methods

What Is a Scalar Integral Over a Triangle?

A scalar integral over a triangle means computing the integral of some function f(x,y) across the area of a 2D triangle. This shows up constantly in finite element methods, computer graphics, and computational physics. If you're doing anything with mesh-based simulations, you'll run into this.

The integral looks like this:

∫∫_T f(x,y) dA

Where T is your triangle. The hard part is that there's no closed-form solution for arbitrary f(x,y). You need numerical methods.

Why This Matters

Triangle integrals are everywhere in practice:

• Computing mass matrices in FEM
• Calculating centroids and moments of inertia for mesh elements
• Rendering with area-weighted sampling
• Solving PDEs with Galerkin methods

If your mesh has millions of triangles, you need fast, accurate integration methods. The wrong approach will kill your performance.

The Three Main Approaches

1. Affine Coordinate Integration

Transform the triangle to a reference triangle using barycentric coordinates. Map point (ξ,η) in the reference triangle [0,1]×[0,1] to your physical triangle T. The Jacobian of this transformation gives you a constant factor.

This works well for polynomials. For a linear function f(x,y) = a + bx + cy, you can solve this analytically. The integral over any triangle becomes:

∫∫_T (a + bx + cy) dA = Area(T) × (a + b·x̄ + c·ȳ)

Where (x̄,ȳ) is the triangle centroid. That's it for linear functions.

2. Gaussian Quadrature

For anything more complex than a polynomial, use Gaussian quadrature. This approximates the integral as a weighted sum of function evaluations at specific points.

For a triangle, you need triangular Gaussian quadrature rules, not the standard 1D rules. These are derived from orthogonal polynomials on triangles.

Common rules:

• 3-point rule: 1st degree exact for polynomials up to degree 1
• 7-point rule: 2nd degree exact
• 12-point rule: 3rd degree exact
• 16-point rule: 4th degree exact

3. Recursive Subdivision

Split the triangle into 4 smaller triangles by connecting edge midpoints. Apply your integration rule to each subtriangle and sum the results.

This approach handles singularities better than fixed quadrature rules. It's adaptive—you can refine more where the function varies rapidly.

Comparison of Methods

How to Implement Gaussian Quadrature on a Triangle

Here's the practical setup. First, define your triangle vertices:

V0 = (x0, y0), V1 = (x1, y1), V2 = (x2, y2)

Map from reference coordinates (α,β) with α,β ≥ 0 and α+β ≤ 1:

P(α,β) = V0 + α·(V1-V0) + β·(V2-V0)

The area element transforms as:

dA = 2·Area(T)·dα·dβ

For the 7-point rule (2nd degree exact):

• Point 1: α=1/3, β=1/3, weight=9/20
• Points 2-7: α=0.6, β=0.2 and permutations, weight=1/20 each

Evaluate f at each mapped point, multiply by weights, sum everything up, multiply by 2×Area(T).

Watch Out For

Jacobian singularities: If your triangle is degenerate (zero area), the Jacobian becomes zero. Your integral will be zero, but division by zero in the mapping will break things. Check for degenerate triangles first.

Conditioning: Thin, stretched triangles cause numerical issues. The Jacobian determinant can become very small, amplifying floating-point errors. Consider transforming to a well-shaped reference triangle.

Function complexity: If f(x,y) has sharp gradients or kinks, low-order quadrature will fail. Adaptive subdivision or high-order rules are your options. High-order rules are usually faster if you know the function is smooth.

Which Method Should You Use?

Use analytic integration for linear polynomials. It's exact and trivial to implement.

Use 7-point Gaussian quadrature for general FEM assembly. It balances accuracy and cost for most engineering applications.

Use recursive subdivision when your integrand has localized features or when you need guaranteed accuracy within a tolerance.

Avoid high-order rules (12+ points) unless you have a specific reason. The performance cost is usually not worth it. Most engineering tolerances are satisfied at 2nd-degree exactness.