Thermoelastic3D: port the Beams3D cleanups (non-cubic support, render, tidy-ups)

[arthurdrake1]

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.

[linear]

SOH-131

[g-braeunlich]

@arthurdrake1 is there a branch / PR for beams3D from which I can backport (assuming you want me to work on this as you assigned me)?

[arthurdrake1]

@arthurdrake1 is there a branch / PR for beams3D from which I can backport (assuming you want me to work on this as you assigned me)?

We now have #257 as you have seen, so you can refer to it for the implementation details. Updated the issue description to match this.

[g-braeunlich]

In #262, I addressed most of the above points.
I still left out volume_fraction -> volume_fraction_error.
Also I did not investigate the 2 "Worth checking" bullet points. I suggest to move those 2 points to a separate issue

[g-braeunlich]

Question (maybe also @markfuge ):
Do we really need a new problem version after renaming volume_fraction -> volume_fraction_error?
I was assuming that the only reason for changing a problem version would be another dataset.

[markfuge]

Do we really need a new problem version after renaming volume_fraction -> volume_fraction_error? I was assuming that the only reason for changing a problem version would be another dataset.

Fundamentally, if the change doesn't affect the solver or optimization routine, such that the dataset, if re-run using the updates, would not change, then you are correct and no new version is needed. However, if the changes would affect either simulate or optimize such that running the same data generation script would lead to different outputs, then we need to bump the problem version and also the dataset version. @arthurdrake1 can comment on the extent to which that is the case here

[arthurdrake1]

Do we really need a new problem version after renaming volume_fraction -> volume_fraction_error? I was assuming that the only reason for changing a problem version would be another dataset.

Fundamentally, if the change doesn't affect the solver or optimization routine, such that the dataset, if re-run using the updates, would not change, then you are correct and no new version is needed. However, if the changes would affect either simulate or optimize such that running the same data generation script would lead to different outputs, then we need to bump the problem version and also the dataset version. @arthurdrake1 can comment on the extent to which that is the case here

Apologies for getting back to this late. My concern was that the actual value of volume_fraction in the HF dataset is not the error, but rather the true converged volume fraction. As a result there's a mismatch with the code, which has volume_fraction as an objective to minimize. What it really should be is the volume_fraction_error, which is currently not reported in the dataset. If we wanted to add that column or change the existing one to reflect error, that would require a v1 if I'm not mistaken.

[g-braeunlich]

@arthurdrake1 So you confirm that running the same data generation script would lead to different outputs?