3. Message passing neural networks for mesh-based simulation

Finite element crash models are, at their core, attributed graphs. Every node carries a feature vector (position, velocity, acceleration, material state), every element edge carries geometric and constitutive information, and the physics of large deformation and contact is expressed as local interactions that propagate across the mesh. Message passing neural networks (MPNNs) formalize this intuition: a learned operator aggregates information from a node's neighborhood, updates the node state, and repeats. This chapter examines the four operator families most relevant to crash simulation — GraphSAGE, GCN, GAT, and custom operators designed for large deformation and contact — and shows how each maps onto the mechanics of an explicit FE solver.

3.1 The Message Passing Abstraction

A single message passing layer computes, for each node i with neighbors j ∈ N(i):

m_ij = φ_msg(h_i, h_j, e_ij)          # edge message
m_i  = ⊕_{j∈N(i)} m_ij                 # permutation-invariant aggregation
h_i' = φ_upd(h_i, m_i)                 # node update

Here h denotes node features, e_ij edge features (rest length, current length, element type), ⊕ is a permutation-invariant aggregator (sum, mean, max), and φ_msg, φ_upd are typically small MLPs. The permutation invariance is not a convenience — it is the discrete analogue of the fact that a finite element assembly is independent of element ordering. Any operator that violates this property will fail to generalize across mesh renumbering.

For crash simulation, the edge features matter as much as the node features. A shell element edge in a Belytschko–Tsay formulation carries membrane strain, bending curvature, and transverse shear; a solid element edge carries the full strain tensor. Encoding these into e_ij is what allows the network to distinguish, for example, a hinge forming in a thin-walled tube from a membrane stretch in the same geometry.

3.2 GraphSAGE: Inductive Aggregation for Varying Meshes

GraphSAGE (Hamilton et al., 2017) computes a fixed-size sample of neighbors and applies a learned aggregator. Its key property for crash models is inductiveness: the operator is defined by the aggregation function, not by a fixed adjacency matrix, so it can be applied to a mesh with different connectivity — a refined crush zone, a remeshed contact patch — without retraining.

The mean aggregator variant is the most natural fit for FE meshes:

h_i' = σ( W · [h_i ‖ mean_{j∈N(i)}(h_j) ] )

where ‖ denotes concatenation. For a crash model with ~10⁶ nodes and ~6×10⁶ edges, a two-layer GraphSAGE with hidden width 128 and mean aggregation is a reasonable starting point. The receptive field after L layers covers the L-hop neighborhood; with L=3, a node sees roughly 3 element rings, which is comparable to the stencil of a single explicit time step in a central-difference scheme.

3.3 GCN: Spectral Smoothing and Its Limits

The graph convolutional network (Kipf & Welling, 2017) applies a first-order Chebyshev approximation of a spectral filter:

H' = σ( D̃^{-1/2} Ã D̃^{-1/2} H W )

with à = A + I and D̃ the degree matrix. This is a low-pass filter on the graph: it smooths node features across neighbors. In a crash model, that smoothing is physically meaningful for displacement and velocity fields, which are continuous, but it is wrong for contact and fracture, where the solution is discontinuous across element boundaries. A GCN layer will smear a contact force across the interface, producing non-physical penetration.

GCN is therefore best used as a feature extractor in the early layers of a hybrid architecture, where it captures smooth deformation modes, and replaced by a discontinuity-aware operator near contact interfaces. Its other limitation is transductivity: the normalized adjacency is tied to a fixed graph, so a GCN trained on one mesh cannot be evaluated on another without recomputing the propagation rule.

3.4 GAT: Attention for Localized Nonlinearity

The graph attention network (Veličković et al., 2018) replaces fixed aggregation weights with learned attention coefficients:

α_ij = softmax_j( LeakyReLU( a^T [W h_i ‖ W h_j] ) )
h_i' = σ( Σ_j α_ij W h_j )

Attention is attractive for crash simulation because it lets the network modulate the influence of neighbors based on local state. A node near a plastic hinge can up-weight the neighbors that carry the bending moment; a node in a contact patch can down-weight neighbors that are separating. In practice, multi-head attention with 4–8 heads and hidden width 64 per head works well on shell meshes with ~10⁵ nodes.

The cost is quadratic in the number of neighbors per node, which is acceptable for FE meshes (typical valence 4–8 for shells, 8–27 for solids) but becomes a bottleneck if the graph is densified by contact search. A common compromise is to apply GAT only within a radius r_c of the contact interface and use a cheaper aggregator elsewhere.

3.5 Custom Operators for Large Deformation and Contact

Standard MPNN operators assume a fixed graph topology. Crash simulation violates this assumption in two ways: large deformation changes edge lengths and orientations continuously, and contact creates and destroys edges as surfaces approach and separate. A custom operator must therefore be geometry-aware and topology-dynamic.

A practical formulation for large deformation uses edge features that include the current length l_ij, the rest length l_ij⁰, and the direction cosine n_ij:

e_ij = [ l_ij, l_ij - l_ij⁰, n_ij, ε_ij ]   # ε_ij = (l_ij - l_ij⁰)/l_ij⁰
m_ij = MLP( h_i, h_j, e_ij )
m_i  = Σ_j m_ij
h_i' = h_i + MLP( h_i, m_i )

The residual update h_i' = h_i + Δ is important: it preserves the identity mapping when the message is zero, which stabilizes training on deep networks and mirrors the incremental nature of an explicit solver.

For contact, the operator must handle a dynamic neighbor set. A workable approach is to maintain a candidate edge list from a broad-phase search (e.g., a bounding-volume hierarchy or a uniform grid) and let the network learn a gating function:

g_ij = σ( MLP( h_i, h_j, d_ij ) )      # d_ij = gap distance
m_ij = g_ij · φ_msg( h_i, h_j, e_ij )

When d_ij exceeds a cutoff, g_ij → 0 and the edge contributes nothing, so the effective topology changes smoothly. This avoids the discontinuity of hard edge insertion/deletion and keeps gradients well-defined during training.

3.6 Choosing an Operator: A Practical Comparison

The table is a starting point, not a prescription. In practice, the strongest crash surrogates use a stacked architecture: a GCN or GraphSAGE layer to extract smooth fields, a GAT layer to localize nonlinearity, and a custom geometry-aware layer near contact. The operator choice is driven by the failure mode you are trying to capture, not by benchmark accuracy on a generic graph dataset.

The diagram above traces the full data path: node and edge features enter one of four operator families, are aggregated by a permutation-invariant reduction, and are combined with the previous node state through a residual update. The contact gating inset shows how a dynamic neighbor set is handled without hard edge insertion or deletion. The key engineering takeaway is that the operator is not a black box — each choice (mean vs. attention vs. gating) corresponds to a specific assumption about the physics, and that assumption should be validated against the failure mode the surrogate is meant to predict.