We're adding a beams3d problem (PR #257), and a number of the things we fixed there turned out to be present in thermoelastic3d too — the two share most of their code. This is a checklist for bringing the same fixes over; once the beams3d PR lands it's the reference for each. The biggest one, the filter bug, has its own issue (#255).

A note on non-cubic grids: simulate() and optimize() already run on non-cubic grids if you hand them properly-shaped inputs (I checked on 4×6×3), and simulate() gives the right answer there. Apart from the filter bug, the two items marked [non-cubic] below are what's left before arbitrary grid sizes are usable through the normal API.

Functional

• [non-cubic] There's no __init__, so you can't configure a problem. ThermoElastic3D doesn't define __init__, so it uses the base Problem.__init__(self, seed=0). That means ThermoElastic3D(config={"nelx": 8, ...}) raises TypeError and self.config is never set. The only way to run a non-cubic grid today is to pass all five boundary masks (shaped (nelx+1, nely+1, nelz+1)) on every call and let the solver infer the grid from them. It also means optimize() falls back to self.Config.max_iter (the class default) instead of an instance value — you can still override per call with optimize(config={"max_iter": N}), but you can't set it at construction. The fix is the __init__(self, seed=0, config=None) that the beams3d PR adds, which sets up self.config, resizes the default masks, and rebuilds design_space/dataset_id for the chosen grid.
• [non-cubic] design_space has the axes in the wrong order (and is stuck at 16³). It's declared spaces.Box(shape=(NELX, NELY, NELZ)), i.e. [x, y, z], but the solver stores the design as (nely, nelx, nelz) and indexes it x[ely, elx, elz]. And since there's no __init__, it stays (16, 16, 16) no matter what grid you run. So a design built from design_space.shape has the wrong shape/axes on any non-cubic grid, and check_constraints, design_space.contains, and random_design are all wrong there. The beams3d PR uses (NELY, NELX, NELZ) and rebuilds it per instance.
• render() needs a Qt backend, and it screenshots the window after it's closed. The thermoelastic3d extra in pyproject.toml lists plain napari with no Qt binding, so a fresh pip install engibench[thermoelastic3d] followed by render() dies with qtpy.QtBindingsNotFoundError. On top of that, render() calls viewer.export_figure() after napari.run() returns, so when you close the window (open_window=True) the canvas is already gone and you get a "C++ object already deleted" error. The beams3d PR uses a plain matplotlib voxel plot with no Qt; the minimal fix here is to add a Qt binding (e.g. napari[pyqt6]) and grab the figure before calling napari.run().
• The volume_fraction objective actually reports the error, not the fraction. It returns abs(mean(x) - volfrac), which heads toward 0, not ~0.3. thermoelastic2d already renamed this to volume_fraction_error (fix(thermoelastic2d): align column names with v1 dataset #247). Changing the objective name needs a new problem version.
• Worth checking: the thermal matrix kth ends up non-symmetric but is solved as if it weren't. Heat-sink boundary conditions are applied by overwriting just the row of each fixed DOF (kth.rows[d]=[d]) and leaving the column alone, which breaks symmetry. It's then solved with solve_spd_with_amg (CG plus smoothed-aggregation AMG), both of which assume a symmetric positive-definite matrix — for the temperature solve and for the thermal adjoint. Someone should confirm this converges to the right field; if not, zero the column too (symmetric elimination) or solve the free-DOF block the way the mechanical side does. This applies to every grid size, not just non-cubic.
• Worth checking: the two compliances come from one combined-load solve. When 0 < weight < 1, the mechanical right-hand side is fp + feps, so the single displacement um (and therefore the reported structural_compliance) includes the thermal load's contribution. The two reported numbers aren't the separate quantities the names suggest. At a minimum this should be documented; better would be to evaluate each compliance under its own load case.

Cosmetic

• The elasticity matrix is written wrong and fixed on the next line. In fem_matrix_builder.py, row 2 of the 6×6 D is written as [lam, lam+2*mu, lam+2*mu, ...] and then corrected with d[2, 1] = lam. The result is correct, but it's easy to break if someone edits the literal and misses the correction. The beams3d PR writes the row correctly.
• Some docstrings still describe the 3D arrays as 2D. A few were copied from the 2D problem — e.g. fe_mthm_bc_3d documents x as a "2D array … (nely, nelx)", and run calls the masks an "NxN binary array". They're really 3D (nely, nelx, nelz) fields and (nelx+1, nely+1, nelz+1) node masks.