Element module

The element module defines the finite element types and the associated reference elements used by the solver. It is responsible for everything specific to elements: types, connectivity, shape functions, gradients, and the jacobian needed to map between reference and physical coordinates.

This module is consumed by higher-level parts of the crate (mesh, quadrature, solver).

Element struct and types

The crate currently supports two classical 2D elements:

P1: 3-node linear triangle (Tri3),
Q1: 4-node bilinear quadrilateral (Quad4).

ElementType encodes which type is used, and Element stores the connectivity via global node indices.

#![allow(unused)]
fn main() {
#[derive(Clone, PartialEq, Eq, Debug)]
pub enum ElementType {
    /// 3-node triangle
    P1,
    /// 4-node quadrangle
    Q1,
}

/// An element stores a vector containing its global indices.
#[derive(Clone, Debug)]
pub struct Element {
    pub indices: Vec<usize>,
}
}

Notes

ElementType derives Clone, PartialEq, Eq, and Debug, which makes pattern matching and testing straightforward.
Element is a tiny container that holds the indices of the mesh vertices forming each element: geometry (actual coordinates) lives in the mesh.

Reference element enum

The reference element encodes the canonical (parameter-space) version of each element type:

Tri3: the unit reference triangle,
Quad4: the unit reference square $[- 1, 1] \times [- 1, 1]$ .

The method num_nodes() returns the number of nodes for each reference element.

#![allow(unused)]
fn main() {
#[derive(Debug, Clone)]
pub enum ReferenceElement {
    /// 3-node reference triangle
    Tri3,
    /// 4-node reference quadrangle
    Quad4,
}

impl ReferenceElement {
    pub fn num_nodes(&self) -> usize {
        match self {
            ReferenceElement::Tri3 => 3,
            ReferenceElement::Quad4 => 4,
        }
    }
}
}

This abstraction separates element formulas (defined on reference coordinates) from the actual geometry (physical coordinates from the mesh).

Reference element implementations

This block implements the core finite element kinematics on the reference element:

shape_functions(&Point2<f64>) -> Vec<f64>
Returns the values of the shape functions ${N_{a}}$ at given local coordinates $(ξ, η)$ .
These are used to interpolate fields (e.g., $u$ ) inside an element:
$u (ξ, η) = a \sum N_{a} (ξ, η) u_{a} .$
shape_gradients(&Point2<f64>) -> Vec<Vector2<f64>>
Returns the gradients on the reference element, i.e., $\partial N_{a} / \partial (ξ, η)$ .
They are combined with the inverse Jacobian to obtain physical gradients $\nabla N_{a}$ during assembly.
jacobian(vertices_coordinates, local_coordinates) -> Matrix2<f64>
Computes the Jacobian $J$ of the mapping from reference to physical coordinates.
For Tri3 this reduces to a constant matrix built from vertex differences; for Quad4 it is obtained by summing contributions of the shape function gradients weighted by vertex coordinates.

#![allow(unused)]
fn main() {
impl ReferenceElement {
    pub fn shape_functions(&self, local_coordinates: &Point2<f64>) -> Vec<f64> {
        match self {
            ReferenceElement::Tri3 => {
                let xi = local_coordinates.x;
                let eta = local_coordinates.y;
                vec![1.0 - xi - eta, xi, eta]
            }
            ReferenceElement::Quad4 => {
                let xi = local_coordinates.x;
                let eta = local_coordinates.y;
                let n1 = 0.25 * (1.0 - xi) * (1.0 - eta);
                let n2 = 0.25 * (1.0 + xi) * (1.0 - eta);
                let n3 = 0.25 * (1.0 + xi) * (1.0 + eta);
                let n4 = 0.25 * (1.0 - xi) * (1.0 + eta);
                vec![n1, n2, n3, n4]
            }
        }
    }

    pub fn shape_gradients(&self, local_coordinates: &Point2<f64>) -> Vec<Vector2<f64>> {
        match self {
            ReferenceElement::Tri3 => {
                vec![
                    Vector2::new(-1.0, 1.0),
                    Vector2::new(1.0, 0.0),
                    Vector2::new(0.0, 1.0),
                ]
            }
            ReferenceElement::Quad4 => {
                let xi = local_coordinates.x;
                let eta = local_coordinates.y;
                let dn1_dxi = -0.25 * (1.0 - eta);
                let dn1_deta = -0.25 * (1.0 - xi);
                let dn2_dxi = 0.25 * (1.0 - eta);
                let dn2_deta = -0.25 * (1.0 + xi);
                let dn3_dxi = 0.25 * (1.0 + eta);
                let dn3_deta = 0.25 * (1.0 + xi);
                let dn4_dxi = -0.25 * (1.0 + eta);
                let dn4_deta = 0.25 * (1.0 - xi);
                vec![
                    Vector2::new(dn1_dxi, dn1_deta),
                    Vector2::new(dn2_dxi, dn2_deta),
                    Vector2::new(dn3_dxi, dn3_deta),
                    Vector2::new(dn4_dxi, dn4_deta),
                ]
            }
        }
    }

    pub fn jacobian(
        &self,
        vertices_coordinates: &[Point2<f64>],
        local_coordinates: &Point2<f64>,
    ) -> Matrix2<f64> {
        match self {
            ReferenceElement::Tri3 => {
                let v0 = vertices_coordinates[0];
                let v1 = vertices_coordinates[1];
                let v2 = vertices_coordinates[2];
                let dx_dxi = v1.x - v0.x;
                let dy_dxi = v1.y - v0.y;
                let dx_deta = v2.x - v0.x;
                let dy_deta = v2.y - v0.y;
                Matrix2::new(dx_dxi, dx_deta, dy_dxi, dy_deta)
            }
            ReferenceElement::Quad4 => {
                let grads = self.shape_gradients(local_coordinates);
                let mut jac = Matrix2::zeros();
                for (grad, vertex) in grads.iter().zip(vertices_coordinates.iter()) {
                    jac[(0, 0)] += grad.x * vertex.x;
                    jac[(1, 0)] += grad.x * vertex.y;
                    jac[(0, 1)] += grad.y * vertex.x;
                    jac[(1, 1)] += grad.y * vertex.y;
                }
                jac
            }
        }
    }
}
}

How it fits into assembly:

Evaluate shape functions $N_{I}$ (and gradients) at quadrature points in reference space.
Build the Jacobian $J$ and its determinant $∣ J ∣$ .
Map reference gradients to physical gradients via $\nabla N_{I} = J^{- T} \nabla_{ξ η} N_{I}$ .
Accumulate local stiffness and load contributions.

Simple unit test

This smoke test checks that:

the number of nodes reported by each reference element is correct,
the shape function vectors have matching sizes.

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

    #[test]
    fn test_reference_element() {
        let tri3 = ReferenceElement::Tri3;
        let quad4 = ReferenceElement::Quad4;

        assert_eq!(tri3.num_nodes(), 3);
        assert_eq!(quad4.num_nodes(), 4);

        let local_coords = Point2::new(0.5, 0.5);
        let tri_shape_funcs = tri3.shape_functions(&local_coords);
        let quad_shape_funcs = quad4.shape_functions(&local_coords);

        assert_eq!(tri_shape_funcs.len(), 3);
        assert_eq!(quad_shape_funcs.len(), 4);
    }
}
}

Keyboard shortcuts

Learn Rust with Applied Sciences

Element module

Element struct and types

Reference element enum

Reference element implementations

Simple unit test