Linear-elastic cantilever

S Structural mechanics Generated 2026-11-06 17:21 UTC 15 plots

Place a fixed budget of material in a clamped beam to make it as stiff as possible, differentiating through a finite-element solve.

Designing a stiff structure. Given a fixed amount of material, where should it go to make a loaded beam as rigid as possible? This is topology optimization, solved by differentiating a finite-element stress analysis with respect to a per-element material density field \(\rho\).

We minimize the compliance (inverse stiffness) of a 3D linear-elastic cantilever beam under the SIMP density-penalization scheme (\(p=3\), \(E_\max = 70{,}000\) MPa). Each element’s stiffness follows the constitutive relation \(E_\text{eff}(\rho) = E_\min + (E_\max - E_\min)\,\rho^3\), and the global stiffness matrix \(K(\rho)\) couples every element to the displacement field. The objective \(C = \mathbf{F}^\top K(\rho)^{-1}\mathbf{F}\) (external work under load \(\mathbf{F}\)) is smooth but non-convex in \(\rho\), so gradient-based optimization drives the design toward sparse, near-binary \(0/1\) material layouts, and the gradient must stay reliable throughout.

How these results were produced

These are example results, produced automatically on GitHub-hosted runners and refreshed on every release. Each solver runs on its intended device: GPU-capable solvers on a Tesla T4 GPU node, CPU-only solvers (OpenFOAM, deal.II, FEniCS, Firedrake, TopOpt.jl) on a separate CPU node. Each benchmark job has its virtual machine to itself, but the hardware and the load on the host can change from run to run, so wall-clock times vary. Read them for relative scaling between solvers rather than as absolute timings. Accuracy and gradient metrics change little across hardware, although float32 arithmetic is not identical across devices and chaotic rollouts amplify small discrepancies. For numbers that reflect your setup, run the benchmark yourself on your target hardware.

Boundary conditions

3D cantilever beam on domain \([0,2]\times[0,1]\times[0,1]\) (HEX8 elements, 2:1:1 aspect). Dirichlet: all nodes at \(x=0\) have zero displacement (clamped wall). Neumann: a prescribed total force is applied to the right face (\(x=2\)), either a uniform downward traction or a concentrated upward corner load depending on the experiment (controlled by the corner_load flag).

Checks

The results overview explains the columns and states.

Solver	Forward	Gradient	Optimization	Forward cost	VJP cost	Time to plateau	Peak memory	Effort
GPU solversForward and VJP, multiples of the fastest GPU solver
JAX-FEM	Pass	Pass	Pass	—	—	—	—	L0
torch-fem	Pass	Pass	Pass	—	—	—	—	L0
CPU solversForward and VJP at N=128, multiples of the fastest CPU solver
FEniCS	Pass	Pass	Pass	3.05 s2.0×	6.02 sfastest	6.1 minfastest	1.07 GiB	L0
Firedrake	Pass	Pass	Pass	—	—	—	—	L0
TopOpt.jl modified setup	Pass	Pass	Pass	—	—	—	—	L0
deal.II	Pass	Not supported	Not supported	1.51 sfastest	—	—	763 MiB	L0

Modified setups (1)

These solvers run a substitute for part of the setup, so their checks hold for the modified setup.

TopOpt.jl: Applies the cantilever load as nodal forces. Effect: None: on the uniform mesh the nodal forces equal the integrated traction.

Flags and exclusions (16)

Firedrake on cost/spatial_cost: not run
Firedrake on cost/temporal_cost: not run
Firedrake on cost/vjp_cost/by_N: not run
Firedrake on cost/vjp_cost/by_steps: not run
JAX-FEM on cost/spatial_cost: not run
JAX-FEM on cost/temporal_cost: not run
JAX-FEM on cost/vjp_cost/by_N: not run
JAX-FEM on cost/vjp_cost/by_steps: not run
TopOpt.jl on cost/spatial_cost: not run
TopOpt.jl on cost/temporal_cost: not run
TopOpt.jl on cost/vjp_cost/by_N: not run
TopOpt.jl on cost/vjp_cost/by_steps: not run
torch-fem on cost/spatial_cost: not run
torch-fem on cost/temporal_cost: not run
torch-fem on cost/vjp_cost/by_N: not run
torch-fem on cost/vjp_cost/by_steps: not run

Not supported (1 solver)

deal.II, 5 experiments. The deal.II structural tesseract is forward-only. It exposes no AD path (no VJP endpoint, no differentiable inputs), so cost/vjp and gradient experiments are categorically out of scope.; The deal.II structural tesseract has no AD path, and topology optimization (compliance gradient-based ρ updates) requires a differentiable compliance output.

Initial conditions

The initial conditions and fields that this setup’s experiments start from. They are plotted without running any solver.

Settings

two_density_bumps

{
  "nx": 16,
  "ny": 2,
  "nz": 8
}

uniform

{
  "rho_0": 0.5,
  "nx": 16
}

Forward

Is the prediction right? Forward experiments compare each solver’s output with a trusted reference, an analytic solution where one exists, and the other solvers.

Agreement between solvers

Structural compliance \(C = \mathbf{F}^\top \mathbf{u}\) vs density \(\rho_0\) at fixed mesh, sweeping uniform density to span the SIMP stiffness regime.

Settings

Sweeps rho_0 ∈ {0.2, 0.4, 0.5, 0.7, 0.9}

{
  "ic": {
    "name": "uniform",
    "seed": 0
  },
  "physics": {
    "nx": 8,
    "ny": 2,
    "nz": 4,
    "Lx": 2.0,
    "Ly": 1.0,
    "Lz": 1.0,
    "F_total": 1.0,
    "corner_load": false,
    "rho_0": 0.2
  },
  "sweep": {
    "key": "rho_0",
    "values": [
      0.2,
      0.4,
      0.5,
      0.7,
      0.9
    ]
  }
}

Agreement with the reference

Structural compliance \(C = \mathbf{F}^\top \mathbf{u}\) vs mesh resolution \(N\) for each solver, uniform density \(\rho_0=0.5\), full-face downward load.

Settings

Sweeps N ∈ {4, 6, 8, 12, 16}

{
  "ic": {
    "name": "uniform",
    "seed": 0
  },
  "physics": {
    "N": 4,
    "ny": 2,
    "nz": 4,
    "Lx": 2.0,
    "Ly": 1.0,
    "Lz": 1.0,
    "F_total": 1.0,
    "corner_load": false
  },
  "sweep": {
    "key": "N",
    "values": [
      4,
      6,
      8,
      12,
      16
    ]
  }
}

Conservation and decay laws

Diagnostic functionals (compliance, total displacement) vs total load \(F_\mathrm{total}\), validating linearity of the SIMP response.

Settings

Sweeps F_total ∈ {0.25, 0.5, 1.0, 2.0, 4.0}

{
  "ic": {
    "name": "uniform",
    "seed": 0
  },
  "physics": {
    "nx": 8,
    "ny": 2,
    "nz": 4,
    "Lx": 2.0,
    "Ly": 1.0,
    "Lz": 1.0,
    "corner_load": false,
    "rho_0": 0.5,
    "F_total": 0.25
  },
  "sweep": {
    "key": "F_total",
    "values": [
      0.25,
      0.5,
      1.0,
      2.0,
      4.0
    ]
  },
  "diagnostics": {
    "compliance": "<callable _get_compliance>"
  }
}

Results by solver

Solver	Mean rel. error
deal.II	0
FEniCS	0
Firedrake	0
TopOpt.jl	0
JAX-FEM	0
torch-fem	0

From Agreement between solvers.

Cost

What does it cost? Wall-clock time and peak memory of the forward and VJP passes as the problem size \(N\) and the number of steps grow. See the note on how these results were produced before reading absolute numbers.

Settings

spatial_cost

Sweeps nx ∈ {4, 6, 8, 12, 16, 24, 32, 48, 64, 96, 128}

{
  "physics": {
    "Lx": 2.0,
    "Ly": 1.0,
    "Lz": 1.0,
    "F_total": 1.0,
    "rho_0": 0.5,
    "corner_load": false,
    "steps": 1,
    "nx": 4
  },
  "cost": {
    "n_trials": 10
  },
  "sweep": {
    "key": "nx",
    "values": [
      4,
      6,
      8,
      12,
      16,
      24,
      32,
      48,
      64,
      96,
      128
    ]
  }
}

temporal_cost

Sweeps steps ∈ {1}

{
  "physics": {
    "Lx": 2.0,
    "Ly": 1.0,
    "Lz": 1.0,
    "F_total": 1.0,
    "rho_0": 0.5,
    "corner_load": false,
    "nx": 8,
    "steps": 1
  },
  "cost": {
    "n_trials": 10
  },
  "sweep": {
    "key": "steps",
    "values": [
      1
    ]
  }
}

Results by solver

Solver	Device	Forward time (N=128)	Forward first call	VJP time (N=128)	VJP / forward	Forward peak memory	VJP peak memory
deal.II	CPU	1.51 s	1.66 s	—	—	763 MiB	—
FEniCS	CPU	3.05 s	6.47 s	6.02 s	1.98×	1.02e+03 MiB	1.07 GiB

From Forward cost across resolution.

Gradient

Is the gradient right? Gradient experiments compare each solver’s AD or adjoint gradient with finite differences, in magnitude (relative \(L^2\) error) and direction (cosine similarity).

Gradient against finite differences

U-curves of finite-difference gradient error vs perturbation size \(\varepsilon\) with subspace cosine, validating VJP correctness on a random density.

Settings

{
  "ic": {
    "name": "random",
    "seed": 0
  },
  "physics": {
    "nx": 8,
    "ny": 2,
    "nz": 4,
    "Lx": 2.0,
    "Ly": 1.0,
    "Lz": 1.0,
    "F_total": 1.0,
    "corner_load": true
  },
  "fd": {
    "eps_values": [
      2.0,
      0.5,
      0.1,
      0.03,
      0.01,
      0.003,
      0.001,
      0.0003,
      0.0001
    ],
    "n_dirs": 6
  }
}

Gradient accuracy across the control parameter

Gradient norm, best-\(\varepsilon\) FD error, direction cosine, and U-curves vs uniform density \(\rho_0\).

Settings

Sweeps rho_0 ∈ {0.2, 0.4, 0.6, 0.8}

{
  "ic": {
    "name": "uniform",
    "seed": 0
  },
  "physics": {
    "nx": 8,
    "ny": 2,
    "nz": 4,
    "Lx": 2.0,
    "Ly": 1.0,
    "Lz": 1.0,
    "F_total": 1.0,
    "corner_load": true,
    "rho_0": 0.2
  },
  "fd": {
    "eps_values": [
      0.5,
      0.1,
      0.03,
      0.01,
      0.003,
      0.001,
      0.0003
    ],
    "n_dirs": 6
  },
  "sweep": {
    "key": "rho_0",
    "values": [
      0.2,
      0.4,
      0.6,
      0.8
    ]
  }
}

Results by solver

Solver	Best-ε FD error	1 − cosine
FEniCS	2.11e-05	1.97e-10
Firedrake	2.11e-05	1.97e-10
TopOpt.jl	2.11e-05	1.97e-10
JAX-FEM	2.11e-05	1.97e-10
torch-fem	2.11e-05	1.97e-10

From Gradient against finite differences.

Optimization

Can you optimize through it? Optimization experiments run a gradient-based optimizer on each solver’s own gradients. A gradient can pass the finite-difference check and still fail to drive an optimization loop.

Topology optimization

SIMP topology optimization on a \(16\times2\times8\) cantilever beam with Adam (lr=0.05): compliance \(C = \mathbf{F}^\top \mathbf{u}\) and density field evolution under a 50% volume-fraction constraint.

Settings

{
  "ic": {
    "name": "uniform",
    "seed": 0
  },
  "physics": {
    "nx": 16,
    "ny": 2,
    "nz": 8,
    "Lx": 2.0,
    "Ly": 1.0,
    "Lz": 1.0,
    "F_total": 1.0,
    "corner_load": true,
    "v_frac": 0.5,
    "compliance_key": "compliance",
    "penalty_weight": 50.0,
    "x_min": 0.001,
    "snap_interval": 5
  },
  "optim": {
    "lr": 0.05,
    "max_iters": 2500,
    "patience": 100
  }
}

Results by solver

Solver	Final compliance	Iterations to plateau	Time to plateau
JAX-FEM	1.00e-03	—	—
FEniCS	1.00e-03	1460	6.1 min
Firedrake	1.00e-03	—	—
TopOpt.jl	1.00e-03	—	—
torch-fem	1.00e-03	—	—

From Topology optimization.