improve handling of polynomial terms

[Davide-Miotti]

At the moment polynomials are handled well in the non-hermite case, but implementing the hermite version brings new challenges, the code does not look great now due to some temporary workarounds.

[Davide-Miotti]

An example of problematic code is this:

bc_i = data.boundary_conditions[i]
        if is_dirichlet(bc_i)
            # Dirichlet boundary: standard polynomial evaluation
            mon(a, xi)
        else
            # Neumann/Robin: α*P + β*∂ₙP
            ni = data.normals[i]
            nmon = length(a)

            # Evaluate polynomial and its normal derivative
            T = eltype(a)
            poly_vals = zeros(T, nmon)
            deriv_vals = zeros(T, nmon)

            mon(poly_vals, xi)
            ∂_normal(mon, ni)(deriv_vals, xi)

            # Apply boundary condition
            @inbounds for k in 1:nmon
                a[k] = α(bc_i) * poly_vals[k] + β(bc_i) * deriv_vals[k]
            end

Inside solve_hermite.jl

[kylebeggs]

was this solved with #53 ?

[Davide-Miotti]

no, still need to be addressed