4. Contact graph construction

In a physics-aware GNN for crash simulation, the contact graph is the structural backbone that tells the network which parts of the geometry can physically interact and how they are oriented relative to one another. Unlike a proximity graph built purely from Euclidean distance, a contact graph must encode the mechanics of contact: the signed gap along a surface normal, the normal direction itself, and the distinction between two separate bodies approaching versus a single body folding onto itself. This section covers the three construction stages — candidate pair detection, edge feature assembly, and self-contact handling — with concrete parameters and implementation logic you can reproduce.

Stage 1: Detecting Candidate Contact Pairs

Brute-force all-pairs testing is infeasible for meshes with 105–106 nodes. Candidate detection therefore uses a broad-phase spatial query followed by a narrow-phase geometric filter. A practical and widely used broad phase is a uniform spatial hash grid sized to the maximum expected contact search radius.

Define the search radius from the physics, not arbitrarily. If the smallest shell thickness in the model is t_min = 1.2 mm and you want to capture contacts that open up to roughly one element edge before impact, a search radius of r_search = 3.0 mm is reasonable. Cell size should be ≥ r_search so that a query only needs to inspect the 27 neighboring cells (3×3×3) around each node.

# Broad phase: spatial hash candidate generation (illustrative)
cell = r_search                      # 3.0 mm
grid = defaultdict(list)
for nid, x in nodes.items():
    grid[tuple((x // cell).astype(int))].append(nid)

candidates = set()
for key, bucket in grid.items():
    neigh = [grid.get((key[0]+i, key[1]+j, key[2]+k))
             for i in (-1,0,1) for j in (-1,0,1) for k in (-1,0,1)]
    for a in bucket:
        for cell_nodes in neigh:
            if not cell_nodes: continue
            for b in cell_nodes:
                if a < b:
                    candidates.add((a, b))

The narrow phase then rejects pairs that cannot be in contact. Two filters matter most. First, a normal-compatibility test: compute the angle between the surface normals n_a and n_b; if n_a · n_b > cos(60°) ≈ 0.5, the surfaces face the same direction and cannot form a contact pair (they are back-to-back or coplanar-facing). Second, a signed gap test: project the separation vector onto the average normal and discard pairs whose gap exceeds r_search. After filtering, a typical automotive B-pillar model with ~80k nodes yields on the order of 2–5×105 candidate pairs per frame — manageable for edge feature computation.

Tune r_search against the coarsest expected deformation increment. If a node can move more than r_search between two saved frames, contacts will be missed; either increase the radius or sub-sample the simulation more densely.

Stage 2: Building Edges with Gap and Normal Features

Each surviving candidate becomes a directed edge (or an undirected edge with symmetric features). The edge feature vector should carry enough geometric information for the message-passing layer to infer contact state without recomputing geometry. A compact, effective feature set is:

FeatureSymbolDefinitionTypical range
Signed gapg(x_b − x_a) · n̄, negative = penetration−2 to +3 mm
Normal alignmentc_nn_a · n_b−1 to 1
Relative normalΔnn_b − n_a (3 components)−2 to 2
Tangential offsetd_tseparation projected onto tangent plane0 to r_search
Edge lengthL‖x_b − x_a‖0 to r_search

The sign convention must be fixed once and documented: here n̄ = normalize(n_a + n_b) points from body A toward body B, and g < 0 means the surfaces have interpenetrated. This convention propagates into the loss function and the contact-force readout, so inconsistency between training and inference is a common source of silent error.

# Narrow phase + edge features (illustrative)
for (a, b) in candidates:
    d = x[b] - x[a]
    nbar = normalize(n[a] + n[b])
    g = dot(d, nbar)
    if g > r_search: continue
    if dot(n[a], n[b]) > 0.5: continue      # same-facing, reject
    dt = norm(d - g * nbar)
    edge_feat[(a,b)] = [g, dot(n[a], n[b]), *(n[b]-n[a]), dt, norm(d)]

Note that g is computed from the current deformed coordinates, so the graph must be rebuilt (or at least re-featured) at every saved frame. A useful engineering shortcut is to keep the topology from a reference frame and only refresh features, rebuilding topology every N frames or when the maximum nodal displacement exceeds 0.5 · r_search.

Stage 3: Self-Contact in Folded and Crushed Geometries

Self-contact is the hard case. When a thin-walled structure crushes, a single connected component folds back onto itself: the two contacting surfaces share the same mesh, often the same material region, and may be only a few element layers apart in the original topology. Naive candidate generation then produces two failure modes — false positives from topologically adjacent nodes that are close in the reference configuration but never actually touch, and missed contacts when the fold brings distant-in-topology but near-in-space surfaces together.

The standard remedy is a topological exclusion mask combined with a geodesic distance gate. Exclude pairs whose graph-hop distance is below a threshold (e.g., 2 hops, to remove immediate neighbors and their neighbors), and additionally exclude pairs whose geodesic distance along the surface is below d_geo_min. For a shell with element size 2 mm, d_geo_min = 6 mm removes the self-pairs created by a single fold while preserving genuine contact between opposing walls of a crushed tube.

Do not use Euclidean distance alone to suppress self-pairs. In a tight fold, two genuinely contacting surfaces can be closer than two adjacent elements on the same wall. The geodesic gate is what separates "near because folded" from "near because adjacent."

For folded geometries, also track the fold line: nodes on opposite sides of a fold have normals that rotate rapidly along the surface. A useful edge feature is the normal curvature across the pair, approximated as ‖Δn‖ / L. Values above roughly 0.3 mm⁻¹ flag a fold-adjacent pair, which the GNN can learn to weight differently from a flat-on-flat impact.

Contact Graph Construction Pipeline 1 · Candidate Detection Spatial hash, cell = r_search 27-cell neighborhood query Normal test: n_a·n_b ≤ 0.5 Gap test: g ≤ r_search ~2–5×10⁵ pairs / frame 2 · Edge Features g = (x_b − x_a)·n̄ (signed) c_n = n_a·n_b, Δn = n_b − n_a d_t = tangential offset L = ‖x_b − x_a‖ rebuild per frame or every N 3 · Self-Contact Hop-distance mask (≥2 hops) Geodesic gate d_geo ≥ 6 mm Fold curvature ‖Δn‖/L > 0.3 Suppress adjacent false pairs preserve true fold contacts Self-contact in a crushed thin-wall section (side view) true contact (g≈0) adjacent pair (masked) n_a n_b outer wall folded inner wall excluded self-pair retained contact edge

Putting It Together: A Reproducible Check

To verify the construction on a real model, run the pipeline on a reference frame and inspect three diagnostics. First, the candidate count per node should be roughly uniform (mean ≈ 6–10 for a shell mesh at r_search = 3 mm); a heavy tail indicates a mis-sized grid. Second, the gap histogram of retained edges should show a clear mode near g ≈ 0 for contacting regions and a decaying tail toward r_search; a spike at g < −t_min signals penetration that the narrow phase should have flagged. Third, for a known fold, confirm that the geodesic gate removed the adjacent false pairs while the true fold contact at g ≈ 0 survived — the SVG above shows exactly this discrimination. If the true contact is missing, lower d_geo_min or check that the fold's two walls are not connected by a short topological path through a shared edge.

Once these diagnostics pass, the contact graph is ready to feed the message-passing layers: each edge carries a signed gap, a normal alignment, and a fold-curvature signal, and the self-contact mask has already removed the topologically adjacent impostors that would otherwise dominate the neighborhood aggregation.