[Bug]: Heat of mixing fails to discretise with a quadratic or quartic particle profile #5876

PyBaMM Version

26.9.0.1.dev22+gf669fd8d9 (main)

Describe the bug

"heat of mixing": "true" with a polynomial particle profile ("quadratic profile" or "quartic profile") fails to discretise with a ShapeError. The same models without the heat of mixing build fine. This happens on main and is not changed by #5769, which only special-cases "uniform profile" (its heat of mixing is zero).

BaseThermal._heat_of_mixing integrates D |grad c|^2 dU/dc over the particle's r domain, using the particle concentration variable. For the polynomial profiles that concentration is reconstructed from the average and surface values instead of being a finite-volume state. Its grad and its r mesh do not match, so the matmul fails (21 vs 19 rows for the SPM, 420 vs 380 for the DFN).

Possible fixes: evaluate the heat of mixing analytically from the polynomial profile, whose dc/dr is known in closed form, or raise a clear OptionError for this combination until that is done.

Steps to Reproduce

import pybamm

for model in [
    pybamm.lithium_ion.SPM({"particle": "quadratic profile", "thermal": "lumped", "heat of mixing": "true"}),
    pybamm.lithium_ion.DFN({"particle": "quartic profile", "thermal": "lumped", "heat of mixing": "true"}),
]:
    pybamm.Simulation(model).build()

Relevant log output

pybamm.expression_tree.exceptions.ShapeError: Cannot find shape (original error: matmul: dimension mismatch with signature (n,k=21),(k=19,m)->(n,m))