Quadrature module

The quadrature module defines how we perform numerical integration over finite elements.
In the FEM assembly, element integrals like $${\int }_{\hat{\Omega }}(\nabla N_a{)}^T(\nabla N_b)|J|d\hat{\Omega }\text{and}{\int }_{\hat{\Omega }}{N}_{a}f|J|d\hat{\Omega }$$ are approximated with quadrature rules (a set of points and weights) on the reference element. The mapping to physical space is handled elsewhere via the Jacobian.

Quadrature struct

A quadrature rule is defined by:

• points: the evaluation points in reference coordinates, and
• weights: the corresponding weights.

#![allow(unused)]
fn main() {
#[derive(Clone, Debug)]
pub struct QuadRule {
    pub points: Vec<Point2<f64>>,
    pub weights: Vec<f64>,
}
}

This lightweight container is used by the solver during element-level integration. For correctness, points.len() must equal weights.len().

Quadrature implementations

Two families of rules are provided:

• triangle(order): simple rules on the reference triangle.

  • order = 1: 1-point rule (centroid), total weight $\frac{1}{2}$ which matches the reference triangle area.
  • order = 2: 3-point rule, exact for linear fields.

• quadrilateral(n): tensor-product Gauss rules on the reference square $[-1,1]\times [-1,1]$.

  • n = 1: 1-point Gauss rule at the center with weight $4$.
  • n = 2: (2\times2) Gauss rule; points at $\pm \frac{1}{\sqrt{3}}$, each with weight $1$.

#![allow(unused)]
fn main() {
impl QuadRule {
    pub fn triangle(order: usize) -> Self {
        match order {
            1 => QuadRule {
                points: vec![Point2::new(1.0 / 3.0, 1.0 / 3.0)],
                weights: vec![0.5],
            },
            2 => QuadRule {
                points: vec![
                    Point2::new(1.0 / 6.0, 1.0 / 6.0),
                    Point2::new(2.0 / 3.0, 1.0 / 6.0),
                    Point2::new(1.0 / 6.0, 2.0 / 3.0),
                ],
                weights: vec![1.0 / 6.0; 3],
            },
            _ => panic!("triangle quadratule of order > 2 not implemented"),
        }
    }

    pub fn quadrilateral(n: usize) -> Self {
        match n {
            1 => QuadRule {
                points: vec![Point2::new(0.0, 0.0)],
                weights: vec![4.0],
            },
            2 => {
                let a = 1.0 / 3.0f64.sqrt();
                let pts = [-a, a];
                let mut points = Vec::with_capacity(4);
                let mut weights = Vec::with_capacity(4);
                for xi in pts {
                    for eta in pts {
                        points.push(Point2::new(xi, eta));
                        weights.push(1.0);
                    }
                }
                QuadRule { points, weights }
            }
            _ => panic!("quadrilateral quadrature with n > 2 points not implemented"),
        }
    }
}
}

Notes

• For triangles, the implementation panics for order > 2 and similarly, quads panic for n > 2. Extend these if you need higher-order elements or exactness.
• The quadrilateral rule builds the 2D grid from the 1D Gauss abscissae $\pm 1/\sqrt{3}$ when n = 2.

Simple tests

The tests check that the rule sizes match expectations for the provided configurations.

#![allow(unused)]
fn main() {
#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_triangle_quadrature() {
        let rule = QuadRule::triangle(2);
        assert_eq!(rule.points.len(), 3);
        assert_eq!(rule.weights.len(), 3);
    }

    #[test]
    fn test_quadrilateral_quadrature() {
        let rule = QuadRule::quadrilateral(2);
        assert_eq!(rule.points.len(), 4);
        assert_eq!(rule.weights.len(), 4);
    }
}
}