Geometric deep learning for computational mechanics Part I: Anisotropic Hyperelasticity

In plain words

Engineers predict how a material bends and stretches with a formula for the stored energy. For a material made of many small crystal grains, the grain shapes, sizes and orientations change the answer, but a few averaged numbers cannot capture how the grains touch each other.

This paper represents each grain as a dot, and each pair of touching grains as a line between dots. The dots carry the grain's volume and orientation. A neural network reads this graph and summarises it in a short list of numbers. A second small network combines that summary with the strain and outputs the energy.

The network is also trained on the stress, which is what the energy's slope gives you. Fitting the slope as well as the value keeps the learned energy smooth, so it does not produce false wiggles in the stress. The authors use the trained model in a crack simulation and report that the grain graph changes where cracks branch. The numbers come from the authors' own simulations, and the paper does not report any error values in the text.

TL;DR

• Input: a right Cauchy-Green strain tensor C, plus a polycrystal grain-contact graph whose nodes carry volume and three Euler angles §2.2 §3.2.
• Model: two graph-convolution layers, two dense layers, and a code of 9 values, concatenated with C into a two-layer ELU MLP that outputs the energy ψ §3.3 §7.1.
• Loss: the H1 Sobolev loss fits ψ and ∂ψ/∂C, which gives the stress, with no weights stated eq. 18 eq. 19.
• Training data come from FFT simulations of 150 random RVEs with 40 to 50 grains each §6.3 §7.
• All results are author-reported and qualitative in the text. The only quantitative results are eCDF curves fig. 12 fig. 13.

Problem

Classical constitutive models describe grain effects with a few averaged state variables, such as void fraction or orientation averages, and these cannot describe how grains connect §1. Neural networks that learn an energy from data can fit a single microstructure, but a model trained on a learned energy alone can have non-smooth derivatives, so the predicted stress can oscillate §4. The paper also assumes that the model should respect frame indifference, isotropy and convexity, and it checks these §5 §7.2. It applies to elastic (path-independent) behaviour only; the authors say history-dependent plasticity is future work §3.3.

Method

Problem setup. The paper models an elastic stored energy ψ(F, G), where F is the deformation gradient and G is a graph of the microstructure. The first Piola-Kirchhoff stress is P = ∂ψ/∂F eq. 1. Because the network takes C = FᵀF as input, frame indifference holds by construction §5.1 eq. 24 eq. 25.

Graph representation. Each grain is a node, and an edge joins each pair of grains that share a surface. The graph is binary and undirected, and only node weights are used §2.2. The node feature matrix X has one row per grain. In the experiments, D = 4 features are used: volume and the three Euler angles of each grain §3.2 §7.1.

Graph Laplacian. The paper builds the unnormalised Laplacian L = D − A, with D the degree matrix and A the adjacency matrix, and then the symmetric normalised form L^sym = D^(−1/2) L D^(−1/2) = I − D^(−1/2) A D^(−1/2) §2.1 eq. 4 eq. 5. Where a node has no neighbours, the entries are undefined; the paper does not say how to handle this eq. 6.

Graph convolution (encoder branch). A layer computes h_next = ReLU(L^sym · h · W + b) §3.2 eq. 16. The paper says it follows Kipf and Welling, with self-loops included, and says L^sym was chosen to avoid numerical instability §3.2. The first layer takes h⁰ = X. Two such layers are used, then two dense layers produce the code §3.3. The code uses a linear or softmax final activation; the authors report similar results for both §3.3.

Hybrid model (energy head). The code is concatenated with C in Voigt form. Two dense hidden layers use the ELU activation §3.1 eq. 13. The output is linear and is the energy ψ̂ §3.1.

Sobolev loss. The training loss has two terms, averaged over N samples. The first is the squared error of the energy eq. 17. The H1 version adds the squared error of ∂ψ/∂C, which is one half of the second Piola-Kirchhoff stress S = 2 ∂ψ/∂C eq. 18 eq. 19. The paper does not add a weight between the terms; both have weight 1 in eq. 18.

# Hybrid energy model  [§3.2] [§3.3] [fig. 6]
# Input:  lap (N, N) matrix from the Laplacian of eq. 5
#         feats (N, 4) volume and Euler angles per grain  [§3.2]
#         C (6,) right Cauchy-Green tensor in Voigt form  [§3.3]
# Output: psi (1,) predicted energy

def encode_graph(lap, feats):
    # Two graph-convolution layers with ReLU  [§3.2] [eq. 16]
    hidden = relu(lap @ feats @ W1 + b1)
    hidden = relu(lap @ hidden @ W2 + b2)
    # Flatten, then two dense layers to a code  [§3.3]
    # The code length of 9 was chosen by the authors  [§7.1]
    code = dense2(dense1(flatten(hidden)))
    return code                              # (9,)

def energy(lap, feats, C):
    # Concatenate the code with C, then two ELU hidden layers  [§3.1] [eq. 13]
    code = encode_graph(lap, feats)
    x = concatenate(code, C)
    h = elu(x @ V1 + c1)
    h = elu(h @ V2 + c2)
    # Linear output  [§3.1]
    return h @ V3 + c3

# Sobolev H1 training loss  [eq. 17] [eq. 18] [eq. 19]
# Input:  C (N, 6) strains, psi_true (N,), S_true (N, 6) stresses
# Output: scalar loss

def h1_loss(model, C, psi_true, S_true):
    # Energy prediction; C must be differentiable  [§4]
    psi_pred = model(C)
    # Stress from the derivative of the prediction  [eq. 18]
    dpsi_pred = grad(sum(psi_pred), C)
    # Target derivative is one half of S  [eq. 19]
    dpsi_true = 0.5 * S_true
    # Both squared errors, averaged over the N samples  [eq. 18]
    energy_error = square(psi_true - psi_pred)
    stress_error = sum(square(dpsi_true - dpsi_pred), axis=-1)
    return mean(energy_error + stress_error)

Mesoscale data. An FFT solver computes the homogenised response of random RVEs. Each RVE is a 49 × 49 × 49 grid with 40 to 50 equiaxed grains from NEPER, and an initial orientation distribution built from uniform and unimodal parts in MTEX §6.1 §6.3 eq. 33 eq. 34 eq. 38. The material is an anisotropic Fung model with constants given in eq. 37 §6.2. Each RVE is loaded with 200 random average strains, with each strain component of F − I in [0, 0.1] §6.3.

Degraded stress for fracture. In the phase-field model, the stored energy is split into tensile ψ⁺ and compressive ψ⁻ parts, and only the tensile part is degraded by g(d) = (1 − d)² eq. 44 eq. 45 §7.3. The split uses J = det F. For J < 1, ψ⁻ = ψ(F_vol) and ψ⁺ = ψ(F) − ψ(F_vol); otherwise ψ⁺ = ψ(F) and ψ⁻ = 0 eq. 49 eq. 50. The degraded stress is P = 2F [g(d) ∂ψ̂⁺/∂C + ∂ψ̂⁻/∂C] eq. 51.

# Degraded first Piola-Kirchhoff stress  [eq. 45] [eq. 46] [eq. 49-51]
# Input:  F (3, 3), d (scalar), G (graph), model returns psi(C, G)
# Output: P (3, 3) stress

def degraded_stress(F, d, G, model):
    # Volumetric part of F  [eq. 46] [eq. 47]
    J = det(F)
    C = transpose(F) @ F
    # C of F_vol is J^(2/3) times identity (derived from eq. 46)
    C_vol = J ** (2 / 3) * identity(3)
    psi_vol = model(C_vol, G)
    # Tensile and compressive energies  [eq. 49] [eq. 50]
    if J >= 1:
        psi_plus = model(C, G)
        psi_minus = 0.0
    else:
        psi_plus = model(C, G) - psi_vol
        psi_minus = psi_vol
    # Degrade only the tensile part  [eq. 45] [eq. 51]
    g = (1 - d) ** 2
    return 2 * F @ (g * grad(psi_plus, C) + grad(psi_minus, C))

Phase-field equations. The paper uses the viscous regularised brittle phase-field model of Miehe et al. §7.3 eq. 40 eq. 41 eq. 42. The driving force H is the maximum of ψ⁺ over the history of the loading eq. 43. The paper does not write the residual parameter r except as approximately zero §7.3.

Hyperparameters

Values are cited to the paper. Where the paper does not give a number, it is "not stated". The text says that training hyperparameters were kept the same across models wherever possible §7.1, but it does not list them.

Results

Author-reported. The paper gives no numeric error values in the text. Quantitative results appear only as eCDF curves in figures, and their values are not stated. Setup: FFT data from anisotropic Fung RVEs, K-fold splits, and a phase-field fracture study §6 §7.

Sources of the baseline labels: M^L2_mlp and M^H1_mlp are feed-forward MLPs with L2 or H1 loss; M^H1_hybrid is the graph model with H1; M^H1_reg adds regularisation tab. 1.

Paper vs code

No official code. The paper gives no code link, and its text contains no project page. I searched GitHub four ways. The search for "hyperelastic graph sobolev" and "geometric learning hyperelasticity polycrystal" returned no repositories. The search for "vlassis" returned unrelated repositories and some third-party projects, none by the authors. The search for "WaiChing Sun" returned only a third-party Claude skill. Nothing was found, so code is omitted.

Pitfalls

• Laplacian choice. Eq. 16 multiplies by L^sym, which is I − D^(−1/2) A D^(−1/2) eq. 5. The text says the GCN uses self-loops, which means  = A + I §3.2 §2.1. The LaTeX source also has a commented-out line for the form D^(−1/2)  D^(−1/2) next to eq. 16 §3.2. The paper does not say which of these the filter is. Implement the self-loop version as the text describes, and check it numerically against the output scale.
• Nodes with no neighbours. A grain with no contacts has degree zero, so D^(−1/2) is undefined eq. 5. The paper does not say how to handle it.
• Missing training details. The optimiser, learning rate, batch size, epochs, dropout rate and L2 coefficient are not stated. The paper says they were kept fixed across models §7.1. A reimplementation cannot reproduce the reported curves from the text.
• Loss scaling. The energy and stress terms have unit weight in eq. 18. The scaled MSE in eq. 20 is used for evaluation. The paper does not say whether training data were scaled, so the two terms may have very different magnitudes.
• Stress convention. Eq. 19 uses S = 2 ∂ψ/∂C, and eq. 18 compares ∂ψ/∂C against S/2. The paper does not say how the 6 Voigt components are normalised in the norm.
• Split derivative. In the compressive branch, ψ(F_vol) depends on C through J. Eq. 51 does not say whether the derivative of ψ(F_vol) includes this dependence. The full chain rule is the consistent choice. This is derived, not stated.
• Extra network calls. Each evaluation of the degraded stress calls the network at F_vol, in addition to F, when J < 1 eq. 49 eq. 50. This cost is derived from the formulas. The paper does not report run time.
• Dataset sizes. The text gives 150 RVEs in §7 and 100 RVEs for the K-fold study §7.1. The fig. 13 caption gives 150 with a 100/50 split. The relation between these sets is not stated. The number of samples is 200 strains per RVE §6.3, which is 30,000 samples for 150 RVEs (derived).
• Convexity is only a necessary check. The paper says the sampled pairwise test is necessary but not sufficient, and it is checked on grid points §7.2.2 eq. 31. Use it only as a sanity check.
• Isotropy is qualitative. The paper shows that rotated RVEs give different responses. It gives no metric §7.2.1 fig. 14.
• Crack results are qualitative. The paper gives crack patterns but no crack-length or branching numbers §7.3 fig. 17 fig. 18. The mesh is coarse and the authors say so §7.3.
• Missing phase-field constants. g_c and η are not stated eq. 41. The residual r is only approximately zero §7.3.
• Equation numbers. The numbers in this card were counted by hand from the LaTeX source, in order of appearance. Check them against the PDF before citing.
• Scope. The paper covers elastic energies only. The history-dependent case is the subject of Part II §3.3 §7.3.

Lineage

• Builds on: Sobolev training (Czarnecki et al., 2017) for the H1 loss §4. Graph convolution in the style of Kipf and Welling (2017) §3.2. Viscous phase-field fracture of Miehe et al. (2010) §7.3. Generalised Fung elasticity for the mesoscale data §6.2. FFT homogenisation via the Lippmann-Schwinger equation §6.1.
• Compared against: M^L2_mlp (L2 loss MLP) and M^H1_mlp (H1 loss MLP) §7.1 tab. 1.
• Superseded by: not stated.