Letting the structure design itself.
We reproduced the classic topology-optimization method — the one behind those organic, bone-like optimal structures — on the physicsbase engine. The FE solve is ours; the loop around it is the algorithm. Out came the iconic optimized cantilever.
Give a mechanical engineer a rectangle of steel, a place to bolt it down, and a load to carry, and eventually they'll ask the question that launched a whole field: where should the material actually go? Not "is this beam strong enough" — but "of all the ways to spend this much metal, which layout is stiffest?" That question is topology optimization, and its answers look uncannily like nature: branching struts, arches, and trusses that no one drew by hand.
The workhorse method is SIMP — Solid Isotropic Material with Penalization — and it is, at heart, a finite-element solve run over and over inside an optimization loop. That structure makes it a perfect stress test for physicsbase, whose entire premise is being the FE kernel an agent calls repeatedly. So we implemented SIMP on top of the engine and pointed it at two textbook problems. Here's how it works, what came out, and why the "solve-in-a-loop" shape matters.
The problem, precisely
Take a design domain and discretise it into finite elements. Give every element e a density variable \(\rho_e\) between 0 (void) and 1 (solid). The stiffness of that element scales with its density through a penalized power law:
The exponent is the trick. With p = 3, a half-dense element gets only ⅛ of the stiffness for half the material — intermediate "grey" densities are a bad deal, so the optimizer is pushed toward crisp black-and-white designs. We then minimise compliance (the structure's strain energy, \(c = \mathbf{F}^{\mathsf T}\mathbf{u}\) — low compliance means high stiffness) subject to a fixed material budget, say 40% of the domain.
At the boundary: what actually crosses the wire
Here is the whole relationship in one exchange. Each iteration, the agent hands physicsbase a description of the current design — geometry, which material each element is made of, where it's held, and what's pushing on it — as a single JSON object. Nothing else. This is a real request against the engine (trimmed to a two-element patch so it fits on the page):
"nodes": [[0,0],[1,0],[2,0],[0,1],[1,1],[2,1]], "materials": { "steel": { "E": 200e9, "nu": 0.3 } }, "elements": [ { "type":"quad4", "nodes":[0,1,4,3], "material":"steel", "section":{"t":0.01}, "id":"e0" }, { "type":"quad4", "nodes":[1,2,5,4], "material":"steel", "section":{"t":0.01}, "id":"e1" } ], "supports": [ {"node":0,"dofs":["ux","uy"]}, {"node":3,"dofs":["ux","uy"]} ], "loads": [ {"node":2,"fy":-5000}, {"node":5,"fy":-5000} ]
physicsbase solves it and returns the answer — the deformation at every node, the reactions at the supports, and the stress in every element. These are the numbers the optimizer reads; this is the real response (rounded):
"displacements": [ { "ux": 0.0, "uy": 0.0 }, // node 0 — held { "ux": -3.03e-5, "uy": -4.33e-5 }, { "ux": -4.04e-5, "uy": -1.271e-4 }, // node 2 — loaded corner /* … nodes 3–5 … */ ], "reactions": [ { "fx": 625.0, "fy": 5000.0 }, /* … */ ], // supports carry the 10 kN "elements": [ { "id": "e0", "stress": [0, 0, -1.00e6], "von_mises": 1.73e6 }, { "id": "e1", "stress": [0, 0, -1.00e6], "von_mises": 1.73e6 } ], "summary": { "n_dof": 12, "max_abs_displacement": 1.27e-4, "solve_ms": 0.9 }
Turning the response back into the next request
This is where the loop lives. From the response, the optimizer computes each element's strain energy — how hard that element is working — directly from the returned displacements. Elements that store a lot of energy are load paths worth reinforcing; elements sitting idle are material worth spending elsewhere. So the agent rewrites the design and sends it straight back:
# each element's material stiffness now reflects its density: # busy elements -> denser (stiffer material) in the next request # idle elements -> lighter (softer material) in the next request response = physicsbase.solve(design) # ← displacements, stresses energy = strain_energy_from(response) # read the numbers design = redistribute(design, energy) # → new material per element # ...repeat ~50 times. every step stands on one real solve.
Fifty requests in, the material has flowed to exactly where the responses said it was needed, and the shapes below appear.
Under the hood: physicsbase as the finite-element kernel
Every one of those requests is a linear FE solve of \(\mathbf{K}(\rho)\,\mathbf{u} = \mathbf{F}\), where the global stiffness is assembled from element stiffnesses scaled by each element's current density. That's exactly what physicsbase does internally — its plane-stress quad4 element and assembly are validated against closed-form solutions. For the tight 50-solve loop in the reference script we call the element formulation directly (skipping the JSON round-trip for speed), but the object crossing the boundary is the same one shown above:
# the FE kernel comes from physicsbase — a validated element from femengine.elements import Quad4 from femengine.materials import Material # 8x8 unit-square plane-stress stiffness, E = 1 KE = Quad4.stiffness([[0,0],[1,0],[1,1],[0,1]], Material(E=1.0, nu=0.3), {"t":1.0, "kind":"plane_stress"}) for it in range(max_iter): E = Emin + x**penal * (E0 - Emin) # SIMP interpolation K = assemble(KE, E, edof) # scale & scatter each element U = solve(K, F, free) # the physicsbase FE solve ce = einsum('ej,jk,ek->e', Ue, KE, Ue) # element strain energy dc = -penal * x**(penal-1) * (E0-Emin) * ce # compliance sensitivity dc = filter(dc) # density filter (kills checkerboards) x = oc_update(x, dc, volfrac) # optimality-criteria step
The three ingredients that make it converge, rather than degenerate into a pixelated checkerboard, are all standard: the compliance sensitivity \(\frac{\partial c}{\partial \rho_e} = -p\,\rho_e^{\,p-1}(E_0 - E_{\min})\,\mathbf{u}_e^{\mathsf T}\mathbf{k}_0\mathbf{u}_e\) (each element's contribution to stiffness change), a density filter that blurs the sensitivities over a small radius to enforce a minimum feature size, and the Optimality-Criteria update, a beautifully simple fixed-point rule that reshuffles material toward where it does the most good while a bisection holds the volume constraint.
The one call that is the physics
It's worth being precise about how little of this is physicsbase and how much of it is. The engine's entire contribution is a single function call:
>>> KE = Quad4.stiffness([[0,0],[1,0],[1,1],[0,1]], ... Material(E=1.0, nu=0.3), ... {"t": 1.0, "kind": "plane_stress"}) >>> KE.shape (8, 8) # 4 corners × 2 DOFs (ux, uy) = 8
That 8×8 matrix is the exact plane-stress stiffness of one element, and it's the same validated quad4 that passes the uniform-strain patch test in the engine's own suite. Two knobs on it matter for real work: section["kind"] chooses plane_stress (a thin sheet, out-of-plane stress zero) or plane_strain (a thick body, out-of-plane strain zero); section["t"] is the thickness. We evaluate it once, on a unit square with E = 1, because on a regular grid every element is geometrically identical — so scaling that single matrix by each element's SIMP modulus is the assembly. There is no other physics in the program.
Turning a response into a sensitivity
The optimizer asks exactly one thing of every solve: the strain energy stored in each element, \(\mathbf{u}_e^{\mathsf T}\mathbf{k}_0\mathbf{u}_e\). Everything needed to compute it is in the response. Gather each element's eight displacement components out of the returned displacements, and contract them with the same KE physicsbase handed you:
# response["displacements"] -> a flat DOF vector U U = np.array([[d["ux"], d["uy"]] for d in resp["displacements"]]).ravel() Ue = U[edof] # (n_elem, 8): each element's DOFs ce = np.einsum('ej,jk,ek->e', Ue, KE, Ue) # strain energy per element # physicsbase also hands back per-element stress & strain directly, # so the same quantity is available as an energy density with no re-derivation: # resp["elements"][e]["stress"], resp["elements"][e]["strain"]
This is the crux of the "usage" story: physicsbase returns fields the agent can act on. The displacements array is indexable by node; the elements array carries stress, strain, and von_mises per element keyed by the id you sent in. Nothing has to be re-derived — the sensitivity is a dot product over numbers the engine already computed.
In-process or over the wire — the same object
There are two ways to reach that solve, and they return identical numbers. For the tight fifty-iteration loop, the reference script imports Quad4 and assembles in-process — no serialization, fastest. An agent orchestrating from elsewhere POSTs the JSON problem shown earlier to /solve and reads the same displacements and elements back. The physics, the element, and the returned fields are the same; you choose the transport that fits where your loop runs.
What came out
We ran two of the canonical benchmarks. First, a cantilever: clamp the left edge, hang a load off the middle of the right edge, spend 40% of the box. The optimizer was handed a uniform grey slab and left to think for 50 iterations (about 15 seconds).
That shape is not designed; it is discovered. The optimizer has independently rebuilt the load path an engineer would draw — a top and bottom chord in tension and compression, diagonals carrying shear, closed triangular cells for stiffness — purely from "minimise compliance." Run the classic MBB beam instead (a simply-supported span pushed down at the top centre) and you get the other signature result:
And the optimization behaves. Compliance falls steeply for the first dozen iterations, then flattens as the topology sets — the design has found its structure and is only refining edges after that.
Why the loop is the point
Topology optimization is the purest example of the pattern physicsbase exists for. Nobody solves this once. It's solve → read the strain energy → redistribute material → solve again, fifty times over, each step depending entirely on the real numbers from the last. Hand an agent an FE library and it has to re-derive assembly, boundary conditions, and sensitivities every time — and it will get one of them subtly wrong. Hand it an FE solve as an API call and the loop becomes trivial to drive: the agent owns the design logic, the engine owns the physics, and every iteration stands on a validated solve.
That's the same argument as sizing a beam or checking a factor of safety, just with fifty times more solves and a far prettier answer.
The complete code
Here is the entire optimizer, exactly as it ran to produce the figures above
— one file, standing on physicsbase's Quad4 for the physics. No
hidden helpers.
"""SIMP compliance topology optimization, built on the physicsbase engine.
physicsbase supplies the finite-element kernel: the 8x8 plane-stress Quad4
element stiffness is taken straight from femengine.elements.Quad4 — the same
validated element used elsewhere in the engine. Everything else (the SIMP
interpolation, compliance sensitivities, density filter, and the
Optimality-Criteria update) is the optimization loop wrapped around it.
"""
import sys
from pathlib import Path
import numpy as np
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import spsolve
sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
from femengine.elements import Quad4 # physicsbase FE kernel
from femengine.materials import Material
def unit_element_stiffness(nu=0.3):
"""8x8 plane-stress stiffness of a unit square Quad4 with E = 1."""
coords = [[0, 0], [1, 0], [1, 1], [0, 1]] # CCW
return Quad4.stiffness(coords, Material(E=1.0, nu=nu),
{"t": 1.0, "kind": "plane_stress"})
def build_filter(nelx, nely, rmin):
"""Linear (cone) density filter weights as a sparse operator."""
nel = nelx * nely
rows, cols, vals = [], [], []
R = int(np.ceil(rmin)) - 1
for i in range(nelx):
for j in range(nely):
e = i * nely + j
for ii in range(max(i - R, 0), min(i + R + 1, nelx)):
for jj in range(max(j - R, 0), min(j + R + 1, nely)):
w = rmin - np.hypot(i - ii, j - jj)
if w > 0:
rows.append(e); cols.append(ii * nely + jj); vals.append(w)
H = coo_matrix((vals, (rows, cols)), shape=(nel, nel)).tocsr()
Hs = np.asarray(H.sum(axis=1)).ravel()
return H, Hs
def optimize(nelx=120, nely=60, volfrac=0.4, penal=3.0, rmin=2.4,
load="cantilever", max_iter=90, move=0.2, tol=0.01):
ndof = 2 * (nelx + 1) * (nely + 1)
KE = unit_element_stiffness()
# element -> global DOF map (node i,j -> id i*(nely+1)+j)
edof = np.zeros((nelx * nely, 8), dtype=int)
for i in range(nelx):
for j in range(nely):
n1 = i * (nely + 1) + j # bottom-left
n2 = (i + 1) * (nely + 1) + j # bottom-right
n3 = (i + 1) * (nely + 1) + j + 1 # top-right
n4 = i * (nely + 1) + j + 1 # top-left
edof[i * nely + j] = [2*n1, 2*n1+1, 2*n2, 2*n2+1,
2*n3, 2*n3+1, 2*n4, 2*n4+1]
iK = np.kron(edof, np.ones((8, 1), int)).flatten()
jK = np.kron(edof, np.ones((1, 8), int)).flatten()
# boundary conditions and unit load
F = np.zeros(ndof); fixed = []
if load == "cantilever": # clamp left edge, load right-mid
for j in range(nely + 1):
fixed += [2*j, 2*j+1]
ntip = nelx * (nely + 1) + (nely // 2)
F[2*ntip + 1] = -1.0
else: # MBB half-beam
for j in range(nely + 1):
fixed.append(2*j) # symmetry: left edge ux = 0
fixed.append(2*(nelx*(nely+1)) + 1) # bottom-right roller: uy = 0
F[2*nely + 1] = -1.0 # load at top-left
fixed = np.array(sorted(set(fixed)), dtype=int)
free = np.setdiff1d(np.arange(ndof), fixed)
H, Hs = build_filter(nelx, nely, rmin)
x = np.full(nelx * nely, volfrac)
for it in range(max_iter):
Emin, E0 = 1e-9, 1.0
E = Emin + x**penal * (E0 - Emin) # SIMP interpolation
sK = (KE.flatten()[np.newaxis]).T * E
K = coo_matrix((sK.flatten(order='F'), (iK, jK)),
shape=(ndof, ndof)).tocsc()
K = 0.5 * (K + K.T)
U = np.zeros(ndof)
U[free] = spsolve(K[free][:, free], F[free]) # the FE solve
Ue = U[edof]
ce = np.einsum('ej,jk,ek->e', Ue, KE, Ue) # element strain energy
comp = float(np.sum(E * ce))
dc = -penal * x**(penal - 1) * (E0 - Emin) * ce # compliance sensitivity
dv = np.ones(nelx * nely)
dc = np.asarray(H @ (dc / Hs)) # density filter
dv = np.asarray(H @ (dv / Hs))
l1, l2 = 0.0, 1e9 # OC update (bisection)
while (l2 - l1) / (l1 + l2 + 1e-30) > 1e-4:
lm = 0.5 * (l1 + l2)
xn = np.clip(x * np.sqrt(np.maximum(-dc, 0) / (dv*lm + 1e-30)),
np.maximum(0.0, x - move), np.minimum(1.0, x + move))
xn = np.clip(xn, 0.001, 1.0)
l1, l2 = (lm, l2) if xn.mean() > volfrac else (l1, lm)
change = float(np.max(np.abs(xn - x)))
x = xn
print(f"it {it:3d} compliance {comp:10.4f} vol {x.mean():.3f}")
if change < tol:
break
return x.reshape(nelx, nely).T
if __name__ == "__main__":
dens = optimize()
np.save("topopt_density.npy", dens)
print("final volume fraction:", dens.mean())That's the whole thing. The one import that matters — from femengine.elements import Quad4 — is where physicsbase does the heavy lifting; the rest is the classic algorithm. Run it, and the structures above appear.
Honest footnotes
What we implemented is the standard density-based method as described across the references below — the mathematics is common to all of them and to any graduate course on the subject. The code is our own, written around physicsbase's element; we did not reuse any of the published implementations. Our version is 2D minimum-compliance with a linear density filter and an OC update; it does not include the parallel solvers, 3D extension, or advanced projection filters that make the cited "large-scale" codes special. The point here was not to beat them — it was to show that a general-purpose FE engine, called in a loop, reproduces the canonical result cleanly. It does.
References
- M. P. Bendsøe & O. Sigmund, Topology Optimization: Theory, Methods, and Applications, Springer (2003) — the standard text for the SIMP method.
- O. Sigmund, "A 99 line topology optimization code written in Matlab," Structural and Multidisciplinary Optimization 21 (2001), 120–127.
- F. Ferrari & O. Sigmund, "A new generation 99 line Matlab code for compliance topology optimization and its extension to 3D," Struct. Multidisc. Optim. (2020). arXiv:2005.05436
- A. Gupta, R. Chowdhury, A. Chakrabarti & T. Rabczuk, "A 55-line code for large-scale parallel topology optimization in 2D and 3D" (2020). arXiv:2012.08208
examples/topopt.py in the engine repo, built entirely on physicsbase's quad4. Point it at your own domain, loads, and volume budget and watch a structure appear. Read the docs →