4. Graph neural networks for unstructured meshes

Finite-element and finite-volume solvers rarely operate on regular grids. Their domains are discretized into unstructured meshes whose nodes sit at arbitrary positions and whose connectivity varies from one region to another — refined near stress concentrations, coarse in the far field. A graph neural network (GNN) is the natural learning architecture for such data: the mesh is a graph, with nodes as vertices and mesh edges as graph edges. This section explains how a GNN processes an unstructured mesh, how message passing and aggregation produce node-level field predictions, and how to design node and edge features so that the network can actually learn stress and displacement fields rather than memorizing geometry.

From mesh to graph: the data contract

Before any learning happens, the mesh must be expressed as a graph G = (V, E). Each node i ∈ V carries a feature vector x_i; each directed edge (i, j) ∈ E carries an edge feature e_ij. For a 2D triangular mesh with N nodes, a typical encoding is:

• Node features: nodal coordinates (x, y), boundary-condition flags (fixed / loaded / free), material identifiers, and any known input field such as applied traction or initial temperature.
• Edge features: the relative position vector r_ij = x_j − x_i, its length ‖r_ij‖, and optionally the edge's incident triangle areas or a boundary-edge flag.
• Targets: per-node displacement (u_x, u_y) and per-node stress components (σ_xx, σ_yy, σ_xy), or per-element stress if the mesh stores stress at element centroids.

A crucial convention must be fixed at the outset: all coordinates and displacements use the same global frame, and edge direction is defined as pointing from the source node i to the target node j. Every feature, message, and aggregation in this section respects that convention. If you later rotate the mesh for augmentation, rotate node coordinates and displacement targets together; never rotate one without the other.

Message passing: the core computation

A GNN layer updates each node by collecting information from its neighbors. The standard formulation has three stages. First, a message function computes what each neighbor sends:

m_ij = φ_m( x_i, x_j, e_ij )

Second, an aggregation function combines all incoming messages at node i:

a_i = □_{j ∈ N(i)} m_ij

where □ is a permutation-invariant operator — sum, mean, or max. Third, an update function produces the new node state:

x_i' = φ_u( x_i, a_i )

Stacking L such layers gives each node a receptive field of L hops. For a mesh, this matters physically: a node's stress depends on the displacement field in its neighborhood, and a deeper GNN can propagate boundary information (a fixed support, a point load) across the domain. In practice, 4–8 layers are enough for meshes with a few thousand nodes; beyond that, oversmoothing flattens node features and predictions lose local contrast.

The message function φ_m and update function φ_u are small MLPs. A practical choice for φ_m is a two-layer MLP with 64–128 hidden units and SiLU activation, taking the concatenation [x_i ‖ x_j ‖ e_ij] as input. The aggregation operator should be sum when the physical quantity is additive (e.g., contributions to a nodal force) and mean when you want scale invariance across nodes with different valence. Max aggregation is rarely appropriate for continuous field quantities because it discards magnitude information.

Aggregation and the readout

After L message-passing layers, each node holds a latent vector h_i. A small readout MLP maps h_i to the predicted field values at that node. Because the task is regression on continuous quantities, use a linear output layer and train with mean-squared error on normalized targets. Normalize each field channel independently: compute the mean and standard deviation of u_x, u_y, σ_xx, σ_yy, σ_xy over the training set, and standardize both inputs and targets. Without this, the loss is dominated by whichever channel has the largest raw magnitude — typically displacement in meters versus stress in pascals differ by many orders of magnitude.

For element-level stress, one common pattern is to predict nodal stress and then average to element centroids using the element's shape-function weights, or to attach a second readout that pools the three (or four) corner nodes of each element. The first approach keeps a single node-level output head and is simpler to implement; the second gives the network direct access to element topology.

Designing node and edge features that carry physics

The single largest determinant of GNN accuracy on mesh data is feature design, not architecture depth. Raw coordinates alone force the network to rediscover differential operators from scratch. Three design choices consistently pay off:

• Relative position on edges. Store r_ij = x_j − x_i rather than absolute coordinates on the edge. This makes messages translation-invariant and lets the MLP learn directional kernels analogous to finite-difference stencils.
• Boundary-condition encoding. A one-hot or integer flag per node (0 = free, 1 = fixed, 2 = prescribed displacement, 3 = loaded) tells the network where constraints live. Without it, a fixed support and an interior node with identical coordinates are indistinguishable.
• Material and geometry scalars. Young's modulus E, Poisson's ratio ν, and element area or edge length give the network the stiffness scale. For heterogeneous meshes, attach material ID per node or per edge.

A reproducible training loop

Consider a 2D cantilever plate discretized into roughly 3,000 triangular elements, with the left edge fixed and a downward traction on the right edge. Generate 200 load cases by varying traction magnitude and direction; solve each with a linear elastic FEM solver to obtain nodal displacement and stress. Split 160/20/20 into train/validation/test.

# Pseudocode — one training step
for batch in loader:                       # batch of graphs
    h = node_encoder(batch.x)              # [N, 64]
    for layer in gnn_layers:               # 6 message-passing layers
        m = mlp_msg(h[src], h[dst], batch.edge_attr)
        a = scatter_sum(m, dst, dim=0)     # sum aggregation
        h = mlp_upd(h, a)
    pred = readout(h)                      # [N, 5] → u_x, u_y, σ_xx, σ_yy, σ_xy
    loss = mse(pred, batch.y_norm)
    loss.backward(); optimizer.step()

Use Adam with learning rate 1e-3, batch size 8–16 graphs, and 200–500 epochs. Monitor validation MSE per channel separately; a single scalar loss hides the case where displacement converges but stress does not. On the test set, report relative L2 error ‖pred − true‖ / ‖true‖ per field. For this class of problem, a well-tuned GNN typically reaches 1–5% relative error on displacement and 5–15% on stress; stress is harder because it is a derivative quantity and amplifies high-frequency error.

What to verify before trusting a prediction

Beyond aggregate error metrics, inspect the predicted field visually on the mesh. A correct displacement field should be smooth and satisfy the prescribed boundary values to within numerical tolerance; a correct stress field should show the expected concentration at the fixed corner and decay toward the free end. If the network predicts a plausible-looking but spatially shifted stress peak, the likely culprit is an edge-direction convention mismatch between the data generator and the model — verify that r_ij is computed as x_j − x_i consistently in both. Finally, test on a mesh with different refinement than the training set: a GNN that generalizes across discretizations has learned the underlying operator, while one that only interpolates the training mesh has memorized it.