diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0596fbfd..bd09addb 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,49 @@
 # Changelog
 
+## Unreleased
+
+### Changed
+
+- A live polar refresh writes into the panel's own storage instead of rebuilding
+  its interpolations. `Interpolations` hands back a gridded interpolation that
+  reads the very arrays it was given when it is built with `interpolate!`, so
+  once a panel's table reads its own knots and values, writing them is the whole
+  update: `set_polar!` on an unchanged table shape allocates nothing and leaves
+  the interpolation objects in place. A table that changes shape is rebuilt
+  once.
+- `refresh_live_polars!` now allocates nothing at all, where a refresh over 60
+  panels cost 15.3 MB before. On one core that buys little — the pass is
+  arithmetic-bound, not allocation-bound — but a sweep running one model per
+  core was spending half its wall time in the garbage collector, which stops
+  every worker. Eight workers refreshing concurrently went from 4.49 to 2.02 ms
+  a refresh, and the collector from 51% of the run to nothing. NeuralFoil's
+  forward pass was the bulk of the garbage: every layer built a fresh matrix for
+  its product, another for its activation, and the symmetry embedding did the
+  whole thing twice more for the flipped case. `NeuralFoilWorkspace` holds one
+  set of layer activations per symmetry and `fused_output!` runs the pass
+  through them, with the bias and activation folded into one in-place sweep, the
+  Mahalanobis penalty subtracted from the confidence logit case by case rather
+  than through a vector of distances, and the flip undone while averaging
+  instead of into a third matrix. `LivePolars` keeps that workspace, the network
+  itself, and a buffer for every per-panel quantity a refresh touches, so the
+  shapes deform in place (`deform_kulfan!`), the coefficients decode in place
+  (`decode_coefficients!`) and the panel takes them as views. `polar_drift`,
+  `live_surface_friction!` and `live_shape_offset!` allocate nothing either.
+- `refresh_live_pressure!` allocates 96% less (3.7 MB to 146 kB a refresh over
+  60 panels): it shares the polar refresh's workspace, reads the edge velocities
+  out of it into storage it holds, and reconstructs each panel's `Cp` through
+  one set of reused buffers (`contour_arc!`, `contour_pressure!`). What is left
+  is the monotone interpolation object itself, one per panel.
+- `set_polar!` no longer allocates when handed views rather than vectors.
+- `KulfanParameters` is now mutable, so a live polar source rewrites one shape
+  per panel every solve instead of building a new one, and the panel pointing at
+  it follows without being told. As an immutable it was boxed on its way into
+  the panel's `live_shape` field, which cost a refresh 48 bytes a panel even
+  when the shape had not changed; a panel's `live_shape` is now the very object
+  it was sampled from rather than a copy that compares equal to it. Two
+  separately built `KulfanParameters` therefore no longer compare `===` or `==`
+  on identical contents.
+
 ## VortexStepMethod v4.3.0 2026-08-31
 
 ### Fixed
diff --git a/docs/make.jl b/docs/make.jl
index 2688d0c7..fd446529 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -14,7 +14,7 @@ makedocs(;
              VortexStepMethod.ObjAdapter],
     authors="Uwe Fechner <uwe.fechner.msc@gmail.com>, Bart van de Lint <bart@vandelint.net> and contributors",
     sitename="VortexStepMethod.jl",
-    format = Documenter.HTML(prettyurls = haskey(ENV, "CI")),
+    format = Documenter.HTML(prettyurls = haskey(ENV, "CI"), size_threshold = 400_000),
     pages=[
         "Home" => "index.md",
         "How it works" => "explanation.md",
diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
index a9ef5612..197866a0 100644
--- a/docs/src/private_functions.md
+++ b/docs/src/private_functions.md
@@ -99,8 +99,10 @@ update_non_deformed_sections!
 ```@docs
 calculate_new_aero_data
 set_polar!
-rebuild_polar
+refresh_polar!
+polar_column!
 polar_model
+reads_from
 same_knots
 polar_knots
 window_alpha
@@ -133,6 +135,8 @@ CurrentModule = VortexStepMethod.AirfoilAero
 
 ### Kulfan CST parametrization
 ```@docs
+deform_kulfan!
+chord_residual!
 bernstein_basis
 class_function
 leading_edge_basis
@@ -160,14 +164,27 @@ load_neuralfoil_model
 neuralfoil_section
 neuralfoil_fused_output
 fused_output
+fused_output!
 decode_surface_velocity
+decode_surface_velocity!
+surface_velocity_rows
 decode_coefficients
+decode_coefficients!
 nn_forward
+nn_forward!
+add_bias!
+layer_buffers
+case_capacity
 prepare_inputs
 fill_case_input!
 flip_inputs
+flip_inputs!
 flip_outputs
+flipped_row
+fuse_flipped!
 squared_mahalanobis_distance
+mahalanobis_case
+penalize_confidence!
 swish
 sigmoid
 ```
@@ -186,10 +203,16 @@ write_node_table
 flat_plate_cf
 neuralfoil_contour_solution
 contour_arc
+contour_arc!
 arc_at_chord
+arc_at_chord!
+sorted_surface!
+sort_pairs!
 trailing_edge_speed
 velocity_knots
+velocity_knots!
 contour_pressure
+contour_pressure!
 fill_node_nans!
 live_xfoil_solver
 xfoil_ramp
diff --git a/docs/src/private_types.md b/docs/src/private_types.md
index ff99a2c5..0941c427 100644
--- a/docs/src/private_types.md
+++ b/docs/src/private_types.md
@@ -31,6 +31,8 @@ LivePolarSettings
 LivePolars
 NeuralFoilModel
 NeuralFoilResult
+NeuralFoilWorkspace
+ContourPressureScratch
 ```
 ```@meta
 CurrentModule = VortexStepMethod
diff --git a/src/airfoil_aero/AirfoilAero.jl b/src/airfoil_aero/AirfoilAero.jl
index 140da02a..16f10813 100644
--- a/src/airfoil_aero/AirfoilAero.jl
+++ b/src/airfoil_aero/AirfoilAero.jl
@@ -16,11 +16,11 @@ include("kulfan.jl")
 include("deform.jl")
 include("shrink_wrap.jl")
 include("neuralfoil.jl")
-include("live_polar.jl")
 include("poly.jl")
 include("airfoil_solvers/common.jl")
 include("airfoil_solvers/xfoil_solver.jl")
 include("airfoil_solvers/neuralfoil_solver.jl")
+include("live_polar.jl")
 include("airfoil_io.jl")
 include("polar_gen.jl")
 include("polar_export.jl")
diff --git a/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl b/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl
index f49eed1e..54f3ba1e 100644
--- a/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl
+++ b/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl
@@ -49,8 +49,17 @@ along the upper surface toward the trailing edge, negative along the lower. The
 natural coordinate across a blunt nose, where a small chordwise step is a long one
 along the skin.
 """
-function contour_arc(x, y, le)
-    arc = zeros(Float64, length(x))
+contour_arc(x, y, le) = contour_arc!(zeros(Float64, length(x)), x, y, le)
+
+"""
+    contour_arc!(arc, x, y, le) -> arc
+
+[`contour_arc`](@ref) written into `arc`, which is resized to the contour. A live
+pressure refresh walks panel after panel through one buffer this way.
+"""
+function contour_arc!(arc::Vector{Float64}, x, y, le)
+    resize!(arc, length(x))
+    arc[le] = 0.0
     for k in (le - 1):-1:1
         arc[k] = arc[k + 1] + hypot(x[k] - x[k + 1], y[k] - y[k + 1])
     end
@@ -60,6 +69,31 @@ function contour_arc(x, y, le)
     return arc
 end
 
+"""
+    ContourPressureScratch()
+
+The buffers a surface-pressure reconstruction reuses, so [`contour_pressure!`](@ref)
+allocates nothing per panel or per frame. Each grows to whatever contour and station
+count it is first handed and stays that size.
+"""
+struct ContourPressureScratch
+    "One surface's node chord fractions, ascending."
+    surface_x::Vector{Float64}
+    "The same nodes' signed arc length."
+    surface_arc::Vector{Float64}
+    "Arc length of the network's stations on the upper surface."
+    upper_arc::Vector{Float64}
+    "Arc length of the network's stations on the lower surface."
+    lower_arc::Vector{Float64}
+    "Arc length of every knot of the joined edge-velocity curve."
+    knots::Vector{Float64}
+    "Signed edge velocity at those knots."
+    values::Vector{Float64}
+end
+
+ContourPressureScratch() = ContourPressureScratch(Float64[], Float64[], Float64[],
+                                                  Float64[], Float64[], Float64[])
+
 """
     arc_at_chord(x, arc, indices, fractions) -> Vector
 
@@ -67,14 +101,66 @@ Signed arc length at each chord fraction, read off the contour nodes `indices` 
 one surface the fractions belong to.
 """
 function arc_at_chord(x, arc, indices, fractions)
-    nodes = collect(indices)
-    order = sortperm(x[nodes])
-    xs, as = x[nodes][order], arc[nodes][order]
-    return [begin
+    scratch = ContourPressureScratch()
+    return arc_at_chord!(zeros(length(fractions)), scratch, x, arc, indices, fractions)
+end
+
+"""
+    arc_at_chord!(out, scratch, x, arc, indices, fractions) -> out
+
+[`arc_at_chord`](@ref) written into `out`, taking the sorted copy of the surface from
+`scratch` rather than allocating one per call.
+"""
+function arc_at_chord!(out::Vector{Float64}, scratch::ContourPressureScratch, x, arc,
+                       indices, fractions)
+    xs, as = sorted_surface!(scratch, x, arc, indices)
+    resize!(out, length(fractions))
+    @inbounds for (i, f) in enumerate(fractions)
         j = clamp(searchsortedfirst(xs, f), 2, length(xs))
         gap = max(xs[j] - xs[j - 1], eps())
-        as[j - 1] + (f - xs[j - 1]) / gap * (as[j] - as[j - 1])
-    end for f in fractions]
+        out[i] = as[j - 1] + (f - xs[j - 1]) / gap * (as[j] - as[j - 1])
+    end
+    return out
+end
+
+"""
+    sorted_surface!(scratch, x, arc, indices) -> (surface_x, surface_arc)
+
+One surface's nodes as chord fraction and signed arc length, ascending in chord, in the
+scratch's own storage. Sorted by insertion, which a surface listed from one end to the
+other is either already in or exactly reversed from, and which needs no scratch of its
+own.
+"""
+function sorted_surface!(scratch::ContourPressureScratch, x, arc, indices)
+    surface_x, surface_arc = scratch.surface_x, scratch.surface_arc
+    n = length(indices)
+    resize!(surface_x, n)
+    resize!(surface_arc, n)
+    @inbounds for (k, node) in enumerate(indices)
+        surface_x[k] = x[node]
+        surface_arc[k] = arc[node]
+    end
+    sort_pairs!(surface_x, surface_arc)
+    return (surface_x, surface_arc)
+end
+
+"""
+    sort_pairs!(sorted, carried)
+
+Insertion-sort `sorted` ascending, carrying `carried` along with it. Stable, in place,
+and linear on input that is already ordered, which both of its callers hand it.
+"""
+function sort_pairs!(sorted::Vector{Float64}, carried::Vector{Float64})
+    @inbounds for k in 2:length(sorted)
+        key, along = sorted[k], carried[k]
+        j = k - 1
+        while j >= 1 && sorted[j] > key
+            sorted[j + 1], carried[j + 1] = sorted[j], carried[j]
+            j -= 1
+        end
+        sorted[j + 1], carried[j + 1] = key, along
+    end
+    return nothing
 end
 
 """
@@ -111,20 +197,51 @@ than extrapolated off the end of one. The stagnation point enters as its own kno
 where `ue` interpolated linearly between those two stations reaches zero — the
 stagnation-point-flow result, `ue` being linear in arc length there.
 """
-function velocity_knots(station_x, ue_upper, ue_lower, x, arc, le)
-    upper_arc = arc_at_chord(x, arc, 1:le, station_x)
-    lower_arc = arc_at_chord(x, arc, le:length(x), station_x)
+velocity_knots(station_x, ue_upper, ue_lower, x, arc, le) =
+    velocity_knots!(ContourPressureScratch(), station_x, ue_upper, ue_lower, x, arc, le)
+
+"""
+    velocity_knots!(scratch, station_x, ue_upper, ue_lower, x, arc, le)
+        -> (knots, values)
+
+[`velocity_knots`](@ref) assembled in `scratch`. The two vectors returned are the
+scratch's own and are overwritten by the next call.
+"""
+function velocity_knots!(scratch::ContourPressureScratch, station_x, ue_upper, ue_lower,
+                         x, arc, le)
+    n_stations = length(station_x)
+    upper_arc = arc_at_chord!(scratch.upper_arc, scratch, x, arc, 1:le, station_x)
+    lower_arc = arc_at_chord!(scratch.lower_arc, scratch, x, arc, le:length(x),
+                              station_x)
     speed = trailing_edge_speed(ue_upper, ue_lower, upper_arc, lower_arc,
                                 arc[1], arc[end])
     lower_first, upper_first = abs(ue_lower[1]), abs(ue_upper[1])
     stagnation = lower_arc[1] + lower_first / max(lower_first + upper_first, eps()) *
                  (upper_arc[1] - lower_arc[1])
-    knots = [arc[end]; reverse(lower_arc); stagnation; upper_arc; arc[1]]
-    values = [-speed; reverse(-abs.(ue_lower)); 0.0; abs.(ue_upper); speed]
-    order = sortperm(knots)
-    knots, values = knots[order], values[order]
-    keep = [1; [k for k in 2:length(knots) if knots[k] > knots[k - 1] + 1e-9]]
-    return knots[keep], values[keep]
+    knots = resize!(scratch.knots, 2 * n_stations + 3)
+    values = resize!(scratch.values, 2 * n_stations + 3)
+    knots[1], values[1] = arc[end], -speed
+    for k in 1:n_stations
+        knots[1 + k], values[1 + k] = lower_arc[n_stations + 1 - k],
+                                      -abs(ue_lower[n_stations + 1 - k])
+    end
+    knots[n_stations + 2], values[n_stations + 2] = stagnation, 0.0
+    for k in 1:n_stations
+        knots[n_stations + 2 + k], values[n_stations + 2 + k] = upper_arc[k],
+                                                                abs(ue_upper[k])
+    end
+    knots[end], values[end] = arc[1], speed
+    sort_pairs!(knots, values)
+    kept, previous = 1, knots[1]
+    for k in 2:length(knots)
+        current = knots[k]
+        if current > previous + 1e-9
+            kept += 1
+            knots[kept], values[kept] = current, values[k]
+        end
+        previous = current
+    end
+    return (resize!(knots, kept), resize!(values, kept))
 end
 
 """
@@ -135,10 +252,24 @@ end
 one-curve edge velocity ([`velocity_knots`](@ref)), interpolates it with a
 shape-preserving monotone cubic in arc length, and squares it into `Cp = 1 - ue²`.
 """
-function contour_pressure(station_x, ue_upper, ue_lower, x, arc, le)
-    knots, values = velocity_knots(station_x, ue_upper, ue_lower, x, arc, le)
+contour_pressure(station_x, ue_upper, ue_lower, x, arc, le) =
+    contour_pressure!(zeros(length(x)), ContourPressureScratch(), station_x, ue_upper,
+                      ue_lower, x, arc, le)
+
+"""
+    contour_pressure!(cp, scratch, station_x, ue_upper, ue_lower, x, arc, le) -> cp
+
+[`contour_pressure`](@ref) written into `cp`, with `scratch` carrying the edge-velocity
+curve. Only the monotone interpolation itself still allocates, once per contour.
+"""
+function contour_pressure!(cp, scratch::ContourPressureScratch, station_x, ue_upper,
+                           ue_lower, x, arc, le)
+    knots, values = velocity_knots!(scratch, station_x, ue_upper, ue_lower, x, arc, le)
     speed = interpolate(knots, values, FritschButlandMonotonicInterpolation())
-    return [1 - speed(clamp(a, first(knots), last(knots)))^2 for a in arc]
+    @inbounds for k in eachindex(cp, arc)
+        cp[k] = 1 - speed(clamp(arc[k], first(knots), last(knots)))^2
+    end
+    return cp
 end
 
 """
diff --git a/src/airfoil_aero/deform.jl b/src/airfoil_aero/deform.jl
index a911a565..30814056 100644
--- a/src/airfoil_aero/deform.jl
+++ b/src/airfoil_aero/deform.jl
@@ -69,16 +69,39 @@ they are the two shape freedoms a chord-referenced deflection cannot resolve.
 function deform_kulfan(basis::KulfanBasis, base::KulfanParameters,
                        upper_deflection::AbstractVector,
                        lower_deflection::AbstractVector)
+    out = KulfanParameters(copy(base.upper_weights), copy(base.lower_weights),
+                           base.leading_edge_weight, base.TE_thickness)
+    return deform_kulfan!(out, basis, base, upper_deflection, lower_deflection,
+                          similar(basis.x))
+end
+
+"""
+    deform_kulfan!(out, basis, base, upper_deflection, lower_deflection, residual)
+        -> out
+
+[`deform_kulfan`](@ref) written into the shape `out` already is, with `residual` as the
+scratch [`chord_residual!`](@ref) fills — the form a live polar refreshes a panel's shape
+with, since it neither allocates nor replaces the object the panel points at.
+"""
+function deform_kulfan!(out::KulfanParameters, basis::KulfanBasis,
+                        base::KulfanParameters, upper_deflection::AbstractVector,
+                        lower_deflection::AbstractVector, residual::AbstractVector)
     length(base.upper_weights) == basis.n_weights || throw(ArgumentError(
         "KulfanBasis has $(basis.n_weights) weights, airfoil has " *
         "$(length(base.upper_weights))."))
     (length(upper_deflection) == length(basis.x) &&
      length(lower_deflection) == length(basis.x)) || throw(ArgumentError(
         "Deflections must be sampled on the basis' $(length(basis.x)) stations."))
-    return KulfanParameters(
-        base.upper_weights .+ basis.projection * chord_residual(basis, upper_deflection),
-        base.lower_weights .+ basis.projection * chord_residual(basis, lower_deflection),
-        base.leading_edge_weight, base.TE_thickness)
+    out.leading_edge_weight = base.leading_edge_weight
+    out.TE_thickness = base.TE_thickness
+    for (weights, base_weights, deflection) in
+            ((out.upper_weights, base.upper_weights, upper_deflection),
+             (out.lower_weights, base.lower_weights, lower_deflection))
+        chord_residual!(residual, basis, deflection)
+        mul!(weights, basis.projection, residual)
+        weights .+= base_weights
+    end
+    return out
 end
 
 """
@@ -94,9 +117,19 @@ magnitude past the ones it is correcting, and the airfoil that comes back is not
 Removes the straight line through the deflection's own endpoints, which
 [`chord_line`](@ref) returns, and gives back what is left.
 """
-function chord_residual(basis::KulfanBasis, deflection::AbstractVector)
+chord_residual(basis::KulfanBasis, deflection::AbstractVector) =
+    chord_residual!(similar(basis.x), basis, deflection)
+
+"""
+    chord_residual!(residual, basis, deflection) -> residual
+
+[`chord_residual`](@ref) written into storage the caller owns.
+"""
+function chord_residual!(residual::AbstractVector, basis::KulfanBasis,
+                         deflection::AbstractVector)
     offset, slope = chord_line(basis, deflection)
-    return deflection .- (offset .+ slope .* basis.x)
+    residual .= deflection .- (offset .+ slope .* basis.x)
+    return residual
 end
 
 """
@@ -119,6 +152,10 @@ end
 deform_kulfan(basis::KulfanBasis, base::KulfanParameters,
               camber::AbstractVector) = deform_kulfan(basis, base, camber, camber)
 
+deform_kulfan!(out::KulfanParameters, basis::KulfanBasis, base::KulfanParameters,
+               camber::AbstractVector, residual::AbstractVector) =
+    deform_kulfan!(out, basis, base, camber, camber, residual)
+
 """
     control_point_deflection(basis, fractions, deflections) -> Vector{Float64}
 
diff --git a/src/airfoil_aero/live_polar.jl b/src/airfoil_aero/live_polar.jl
index 215eccbc..ee98551b 100644
--- a/src/airfoil_aero/live_polar.jl
+++ b/src/airfoil_aero/live_polar.jl
@@ -32,9 +32,10 @@ end
     LivePolars(base; settings=LivePolarSettings(), n_stations=60)
 
 Live polar source for a wing whose panels each carry one undeformed airfoil in `base`.
-Holds the fixed CST basis, the per-panel reference angles and the network input scratch.
-A refresh is dominated by the forward pass itself; the deformation is a matvec against
-the constant basis, about half a microsecond a panel. Drive it with
+Holds the fixed CST basis, the per-panel reference angles, the network and every buffer a
+refresh needs, so a refresh writes through storage that is already there and allocates
+nothing at all. A refresh is dominated by the forward pass itself; the deformation is a
+matvec against the constant basis, about half a microsecond a panel. Drive it with
 [`refresh_live_polars!`](@ref).
 """
 mutable struct LivePolars
@@ -44,7 +45,7 @@ mutable struct LivePolars
     basis::KulfanBasis
     "Undeformed Kulfan parameters per panel."
     base::Vector{KulfanParameters}
-    "Deformed Kulfan parameters per panel, rewritten every refresh."
+    "Deformed Kulfan parameters per panel, rewritten in place every refresh."
     deformed::Vector{KulfanParameters}
     "Reference angle [rad] the last refresh sampled about, per panel."
     alpha_ref::Vector{Float64}
@@ -56,6 +57,34 @@ mutable struct LivePolars
     pressure_inputs::Matrix{Float32}
     "Lowest confidence NeuralFoil reported over each panel's samples, last refresh."
     confidence::Vector{Float64}
+    "The network the settings name, held so a refresh never goes to the model cache."
+    model::NeuralFoilModel
+    "Forward-pass scratch, sized for the sampled batch and so wide enough for either pass."
+    work::NeuralFoilWorkspace
+    "Per-panel Reynolds number of the current pass."
+    reynolds::Vector{Float64}
+    "Per-panel angle [rad] the surface-pressure pass was evaluated at."
+    alpha_at::Vector{Float64}
+    "Lift coefficient of every sample of the last refresh."
+    sample_cl::Vector{Float64}
+    "Drag coefficient of every sample of the last refresh."
+    sample_cd::Vector{Float64}
+    "Moment coefficient of every sample of the last refresh."
+    sample_cm::Vector{Float64}
+    "Analysis confidence of every sample of the last refresh."
+    sample_confidence::Vector{Float64}
+    "Upper-surface edge velocity per station and panel, last surface-pressure pass."
+    ue_upper::Matrix{Float64}
+    "Lower-surface edge velocity per station and panel, last surface-pressure pass."
+    ue_lower::Matrix{Float64}
+    "Signed arc length of the contour one panel's pressure is being spread over."
+    arc::Vector{Float64}
+    "Buffers the pressure reconstruction walks from panel to panel through."
+    pressure_scratch::ContourPressureScratch
+    "The representable part of one panel's deflection."
+    residual::Vector{Float64}
+    "One panel's weight change off its base shape."
+    weight_delta::Vector{Float64}
 end
 
 function LivePolars(base::AbstractVector{KulfanParameters};
@@ -69,13 +98,23 @@ function LivePolars(base::AbstractVector{KulfanParameters};
     first(offsets) <= 0 <= last(offsets) || throw(ArgumentError(
         "LivePolars sample offsets must straddle zero, so the reference angle lies " *
         "inside the sampled range; got $(rad2deg.(extrema(offsets))) deg."))
-    n_panels = length(base)
-    return LivePolars(settings, KulfanBasis(; n_stations,
-                                            n_weights=length(base[1].upper_weights)),
-                      collect(base), collect(base), zeros(n_panels),
-                      zeros(length(offsets)),
-                      zeros(Float32, 25, n_panels * length(offsets)),
-                      zeros(Float32, 25, n_panels), zeros(n_panels))
+    n_panels, n_samples = length(base), length(offsets)
+    n_weights = length(base[1].upper_weights)
+    n_cases = n_panels * n_samples
+    model = load_neuralfoil_model(settings.model_size; weights_dir=settings.weights_dir)
+    work = NeuralFoilWorkspace(model, n_cases)
+    n_network_stations = length(work.stations)
+    deformed = [KulfanParameters(copy(p.upper_weights), copy(p.lower_weights),
+                                 p.leading_edge_weight, p.TE_thickness) for p in base]
+    return LivePolars(settings, KulfanBasis(; n_stations, n_weights),
+                      collect(base), deformed, zeros(n_panels), zeros(n_samples),
+                      zeros(Float32, 25, n_cases), zeros(Float32, 25, n_panels),
+                      zeros(n_panels), model, work, zeros(n_panels), zeros(n_panels),
+                      zeros(n_cases), zeros(n_cases), zeros(n_cases), zeros(n_cases),
+                      zeros(n_network_stations, n_panels),
+                      zeros(n_network_stations, n_panels),
+                      Float64[], ContourPressureScratch(), zeros(n_stations),
+                      zeros(n_weights))
 end
 
 """
@@ -114,21 +153,32 @@ function polar_drift(live::LivePolars, alpha::AbstractVector)
         "polar_drift: $(length(alpha)) angles for $(length(live.alpha_ref)) panels."))
     offsets = live.settings.offsets
     reach = min(-first(offsets), last(offsets))
-    return maximum(abs.(alpha .- live.alpha_ref)) / reach
+    return maximum(i -> abs(alpha[i] - live.alpha_ref[i]),
+                   eachindex(alpha, live.alpha_ref)) / reach
 end
 
 """
     deform_live_shapes!(live::LivePolars, deflection) -> Vector{KulfanParameters}
 
 Deform every base airfoil by its camber increment and store the result in
-`live.deformed`, which is returned. `nothing` leaves the base shapes in place.
-The deformation is an analytic perturbation of one fixed weight vector, so it
-never refits and never inherits the non-uniqueness of a fit.
+`live.deformed`, which is returned. `nothing` writes the base shapes back. The
+deformation is an analytic perturbation of one fixed weight vector, so it never refits
+and never inherits the non-uniqueness of a fit.
+
+Each entry keeps the object it already was, its weights overwritten, so a panel holding
+one from an earlier refresh follows the current shape and no frame allocates a new one.
 """
 function deform_live_shapes!(live::LivePolars, deflection)
     for i in eachindex(live.base)
-        live.deformed[i] = isnothing(deflection) ? live.base[i] :
-            deform_kulfan(live.basis, live.base[i], deflection[i])
+        shape, base = live.deformed[i], live.base[i]
+        if isnothing(deflection)
+            shape.upper_weights .= base.upper_weights
+            shape.lower_weights .= base.lower_weights
+            shape.leading_edge_weight = base.leading_edge_weight
+            shape.TE_thickness = base.TE_thickness
+        else
+            deform_kulfan!(shape, live.basis, base, deflection[i], live.residual)
+        end
     end
     return live.deformed
 end
@@ -162,9 +212,10 @@ deform the base airfoil by `deflection` (a chord-normalized deflection on `live.
 or `nothing` to keep the base shape), evaluate NeuralFoil at `alpha_ref .+ offsets`, and
 hand those values straight to the panel.
 
-The write is in place — same knot count, same vectors — so a refresh every solve costs
-the forward pass and the three interpolation rebuilds `Interpolations` needs to take the
-new values, which it copies rather than references.
+The write is in place — same knot count, same vectors — and so is everything around it:
+the shapes, the network inputs, both symmetries of the forward pass and the decoded
+coefficients all land in storage `live` already holds, so a refresh every solve costs the
+forward pass and nothing else.
 
 Each panel keeps the deformed shape it was evaluated at as its `live_shape`, so a plot
 draws the airfoil the network actually saw.
@@ -182,9 +233,8 @@ function refresh_live_polars!(live::LivePolars, panels, alpha_ref, reynolds;
     n_panels = length(live.base)
     length(panels) == n_panels || throw(ArgumentError(
         "refresh_live_polars!: $(length(panels)) panels for $n_panels base airfoils."))
-    per_panel(v) = v isa Number ? fill(float(v), n_panels) : collect(float.(v))
-    alpha_vec, re_vec = per_panel(alpha_ref), per_panel(reynolds)
-    live.alpha_ref .= alpha_vec
+    live.alpha_ref .= alpha_ref
+    live.reynolds .= reynolds
 
     offsets = live.settings.offsets
     n_samples = length(offsets)
@@ -192,23 +242,24 @@ function refresh_live_polars!(live::LivePolars, panels, alpha_ref, reynolds;
     for i in 1:n_panels
         for k in 1:n_samples
             fill_case_input!(live.inputs, (i - 1) * n_samples + k, live.deformed[i],
-                             rad2deg(alpha_vec[i] + offsets[k]), re_vec[i],
+                             rad2deg(live.alpha_ref[i] + offsets[k]), live.reynolds[i],
                              live.settings.n_crit, 1.0, 1.0)
         end
     end
 
-    model = load_neuralfoil_model(live.settings.model_size;
-                                  weights_dir=live.settings.weights_dir)
-    cl, cd, cm, confidence = decode_coefficients(fused_output(live.inputs, model))
+    decode_coefficients!(live.sample_cl, live.sample_cd, live.sample_cm,
+                         live.sample_confidence,
+                         fused_output!(live.work, live.inputs, live.model))
 
     for i in 1:n_panels
         samples = ((i - 1) * n_samples + 1):(i * n_samples)
-        live.knots .= alpha_vec[i] .+ offsets
-        live.confidence[i] = minimum(view(confidence, samples))
-        set_polar!(panels[i], live.knots, view(cl, samples), view(cd, samples),
-                   view(cm, samples); shape=live.deformed[i])
+        live.knots .= live.alpha_ref[i] .+ offsets
+        live.confidence[i] = minimum(view(live.sample_confidence, samples))
+        set_polar!(panels[i], live.knots, view(live.sample_cl, samples),
+                   view(live.sample_cd, samples), view(live.sample_cm, samples);
+                   shape=live.deformed[i])
     end
-    return minimum(confidence)
+    return minimum(live.sample_confidence)
 end
 
 """
@@ -293,7 +344,8 @@ panel forces follow the deformed shape; without this the pattern that spreads th
 forces over the structure would still come from the undeformed section, so a
 deformation would change how hard a panel pulls but not where it pulls.
 
-One batched forward pass over all panels, at the converged angle of attack rather
+One batched forward pass over all panels — through the same workspace the polar refresh
+uses, whose inputs it has already consumed — at the converged angle of attack rather
 than at the sampled ones — a panel's `Cp` is wanted at exactly one angle, and
 evaluating there is both cheaper than storing the samples and exact. Call it after
 the solve has converged, with the contours the traction pattern is indexed on —
@@ -308,22 +360,20 @@ function refresh_live_pressure!(cp, live::LivePolars, contour_x, contour_y,
         length(leading_edge) == n_panels ||
         throw(ArgumentError("refresh_live_pressure!: $(length(cp)) pressure and " *
             "$(length(contour_x)) contour entries for $n_panels panels."))
-    per_panel(v) = v isa Number ? fill(float(v), n_panels) : collect(float.(v))
-    alpha_vec, re_vec = per_panel(alpha), per_panel(reynolds)
+    live.alpha_at .= alpha
+    live.reynolds .= reynolds
     for i in 1:n_panels
         fill_case_input!(live.pressure_inputs, i, live.deformed[i],
-                         rad2deg(alpha_vec[i]), re_vec[i], live.settings.n_crit,
-                         1.0, 1.0)
+                         rad2deg(live.alpha_at[i]), live.reynolds[i],
+                         live.settings.n_crit, 1.0, 1.0)
     end
-    model = load_neuralfoil_model(live.settings.model_size;
-                                  weights_dir=live.settings.weights_dir)
-    station_x, ue_upper, ue_lower = decode_surface_velocity(
-        fused_output(live.pressure_inputs, model))
+    decode_surface_velocity!(live.ue_upper, live.ue_lower,
+                             fused_output!(live.work, live.pressure_inputs, live.model))
     for i in 1:n_panels
-        arc = contour_arc(contour_x[i], contour_y[i], leading_edge[i])
-        cp[i] .= contour_pressure(station_x, view(ue_upper, :, i),
-                                  view(ue_lower, :, i), contour_x[i], arc,
-                                  leading_edge[i])
+        arc = contour_arc!(live.arc, contour_x[i], contour_y[i], leading_edge[i])
+        contour_pressure!(cp[i], live.pressure_scratch, live.work.stations,
+                          view(live.ue_upper, :, i), view(live.ue_lower, :, i),
+                          contour_x[i], arc, leading_edge[i])
     end
     return cp
 end
@@ -342,11 +392,12 @@ function live_surface_friction!(cf, contour_x, reynolds)
     length(cf) == length(contour_x) || throw(ArgumentError(
         "live_surface_friction!: $(length(cf)) friction and " *
         "$(length(contour_x)) contour entries."))
-    re_vec = reynolds isa Number ? fill(float(reynolds), length(cf)) :
-        collect(float.(reynolds))
     for i in eachindex(cf)
-        cf[i] .= (flat_plate_cf(clamp(x, 0.0, 1.0), re_vec[i])
-                  for x in contour_x[i])
+        re = reynolds isa Number ? float(reynolds) : float(reynolds[i])
+        nodes = contour_x[i]
+        @inbounds for k in eachindex(cf[i], nodes)
+            cf[i][k] = flat_plate_cf(clamp(nodes[k], 0.0, 1.0), re)
+        end
     end
     return cf
 end
@@ -380,8 +431,8 @@ function live_shape_offset!(offset, live::LivePolars, shape)
         "live_shape_offset!: $(length(offset)) offset and $(length(shape)) shape " *
         "entries for $(length(live.base)) panels."))
     for i in eachindex(live.base)
-        mul!(offset[i], shape[i],
-             live.deformed[i].upper_weights .- live.base[i].upper_weights)
+        live.weight_delta .= live.deformed[i].upper_weights .- live.base[i].upper_weights
+        mul!(offset[i], shape[i], live.weight_delta)
     end
     return offset
 end
diff --git a/src/airfoil_aero/neuralfoil.jl b/src/airfoil_aero/neuralfoil.jl
index 080a3919..4a229bdb 100644
--- a/src/airfoil_aero/neuralfoil.jl
+++ b/src/airfoil_aero/neuralfoil.jl
@@ -121,6 +121,50 @@ function load_neuralfoil_model(model_size::String="xlarge"; weights_dir=nothing)
     return model
 end
 
+"""
+    layer_buffers(model, n_cases) -> Vector{Matrix{Float32}}
+
+One activation matrix per network layer, `n_cases` wide: the storage a forward pass
+writes through. [`NeuralFoilWorkspace`](@ref) holds two sets of them, one per symmetry.
+"""
+layer_buffers(model::NeuralFoilModel, n_cases::Int) =
+    [zeros(Float32, size(w, 1), n_cases) for w in model.weights]
+
+"""
+    NeuralFoilWorkspace(model, n_cases)
+
+Scratch a symmetry-fused forward pass runs inside, so a caller evaluating a batch of the
+same width every solve allocates nothing at all. Sized once against `model` for at most
+`n_cases` cases; [`fused_output!`](@ref) takes any batch up to that width, which is what
+lets one workspace serve both a sampled polar batch and the narrower pressure batch.
+"""
+struct NeuralFoilWorkspace
+    "Layer activations of the direct pass."
+    direct::Vector{Matrix{Float32}}
+    "Layer activations of the top/bottom-flipped pass."
+    flipped::Vector{Matrix{Float32}}
+    "The flipped input batch."
+    inputs_flip::Matrix{Float32}
+    "One input case less the training mean."
+    centered::Vector{Float64}
+    "NeuralFoil's boundary-layer station x/c, the axis a surface reconstruction reads."
+    stations::Vector{Float64}
+end
+
+function NeuralFoilWorkspace(model::NeuralFoilModel, n_cases::Int)
+    n_cases >= 1 || throw(ArgumentError(
+        "A NeuralFoilWorkspace needs at least one case; got $n_cases."))
+    n_outputs = size(last(model.weights), 1)
+    return NeuralFoilWorkspace(layer_buffers(model, n_cases),
+                               layer_buffers(model, n_cases),
+                               zeros(Float32, model.n_inputs, n_cases),
+                               zeros(model.n_inputs),
+                               compute_optimal_x_points((n_outputs - 6) ÷ 6))
+end
+
+"How many cases a workspace was sized for."
+case_capacity(work::NeuralFoilWorkspace) = size(work.inputs_flip, 2)
+
 """
     nn_forward(x::AbstractMatrix, model::NeuralFoilModel)
 
@@ -133,36 +177,92 @@ Forward pass through the neural network.
 # Returns
 - Output matrix of shape (n_outputs, n_cases)
 """
-function nn_forward(x::AbstractMatrix{T}, model::NeuralFoilModel) where T
-    n_layers = length(model.weights)
+nn_forward(x::AbstractMatrix, model::NeuralFoilModel) =
+    Matrix(nn_forward!(layer_buffers(model, size(x, 2)), x, model))
 
-    for i in 1:n_layers
-        x = model.weights[i] * x .+ model.biases[i]
-        if i < n_layers  # No activation on last layer
-            x = swish(x)
-        end
+"""
+    nn_forward!(layers, x, model, n_cases=size(x, 2)) -> AbstractMatrix
+
+Forward pass writing each layer's activations into `layers`, returning a view of the last
+one over the `n_cases` columns used. The matrix products go through BLAS in place, so a
+pass over a batch the buffers were sized for allocates nothing.
+"""
+function nn_forward!(layers::Vector{Matrix{Float32}}, x::AbstractMatrix,
+                     model::NeuralFoilModel, n_cases::Int=size(x, 2))
+    n_layers = length(model.weights)
+    out = view(layers[1], :, 1:n_cases)
+    mul!(out, model.weights[1], view(x, :, 1:n_cases))
+    add_bias!(out, model.biases[1], n_layers > 1)
+    for i in 2:n_layers
+        out = view(layers[i], :, 1:n_cases)
+        mul!(out, model.weights[i], view(layers[i - 1], :, 1:n_cases))
+        add_bias!(out, model.biases[i], i < n_layers)
     end
+    return view(layers[n_layers], :, 1:n_cases)
+end
 
-    return x
+"""
+    add_bias!(out, bias, activate) -> out
+
+Add a layer's bias to its activations in place, passing them through [`swish`](@ref)
+unless this is the output layer, which carries no activation.
+"""
+function add_bias!(out::AbstractMatrix, bias::AbstractVector, activate::Bool)
+    @inbounds for case in axes(out, 2), k in eachindex(bias)
+        value = out[k, case] + bias[k]
+        out[k, case] = activate ? swish(value) : value
+    end
+    return out
 end
 
 """
     squared_mahalanobis_distance(x::AbstractMatrix, model::NeuralFoilModel)
 
-Compute squared Mahalanobis distance from training distribution.
+Compute squared Mahalanobis distance from training distribution, one per case.
 
 This is used to penalize predictions far from the training data.
 """
 function squared_mahalanobis_distance(x::AbstractMatrix, model::NeuralFoilModel)
-    n_cases = size(x, 2)
-    result = zeros(n_cases)
+    centered = zeros(model.n_inputs)
+    return [mahalanobis_case(x, case, model, centered) for case in axes(x, 2)]
+end
 
-    for i in 1:n_cases
-        x_minus_mean = x[:, i] .- model.mean_inputs
-        result[i] = dot(x_minus_mean, model.inv_cov_inputs * x_minus_mean)
+"""
+    mahalanobis_case(x, case, model, centered) -> Float64
+
+One case's squared Mahalanobis distance, taking `centered` as the scratch the centred
+input column is written into rather than allocating one per case.
+"""
+function mahalanobis_case(x::AbstractMatrix, case::Int, model::NeuralFoilModel,
+                          centered::Vector{Float64})
+    inv_cov = model.inv_cov_inputs
+    @inbounds for i in eachindex(centered)
+        centered[i] = x[i, case] - model.mean_inputs[i]
     end
+    distance = 0.0
+    @inbounds for j in eachindex(centered)
+        weighted = 0.0
+        for i in eachindex(centered)
+            weighted += inv_cov[i, j] * centered[i]
+        end
+        distance += weighted * centered[j]
+    end
+    return distance
+end
+
+"""
+    penalize_confidence!(y, x, model, centered) -> y
 
-    return result
+Subtract every case's Mahalanobis penalty from the confidence logit, row 1 of the network
+output. This is what makes a shape far from the training distribution report a low
+confidence, and it is applied to each symmetry before the two are fused.
+"""
+function penalize_confidence!(y::AbstractMatrix, x::AbstractMatrix,
+                              model::NeuralFoilModel, centered::Vector{Float64})
+    @inbounds for case in axes(y, 2)
+        y[1, case] -= mahalanobis_case(x, case, model, centered) / (2 * model.n_inputs)
+    end
+    return y
 end
 
 """
@@ -242,23 +342,27 @@ end
 
 Flip inputs for symmetry embedding (swap upper/lower, negate alpha).
 """
-function flip_inputs(x::AbstractMatrix{T}) where T
-    x_flip = copy(x)
-
-    # Swap upper and lower weights with sign flip
-    x_flip[1:8, :] .= -x[9:16, :]   # upper <- -lower
-    x_flip[9:16, :] .= -x[1:8, :]   # lower <- -upper
-
-    # Flip LE weight
-    x_flip[17, :] .= -x[17, :]
-
-    # Flip sin(2*alpha) (index 19)
-    x_flip[19, :] .= -x[19, :]
+flip_inputs(x::AbstractMatrix) = flip_inputs!(similar(x), x)
 
-    # Swap transition locations
-    x_flip[24, :] .= x[25, :]  # xtr_upper <- xtr_lower
-    x_flip[25, :] .= x[24, :]  # xtr_lower <- xtr_upper
+"""
+    flip_inputs!(x_flip, x) -> x_flip
 
+[`flip_inputs`](@ref) written into storage the caller owns: the upper and lower weights
+swap with a sign change, the leading-edge weight and `sin(2·alpha)` negate, and the two
+transition locations swap. Everything else carries over.
+"""
+function flip_inputs!(x_flip::AbstractMatrix, x::AbstractMatrix)
+    x_flip .= x
+    @inbounds for case in axes(x, 2)
+        for i in 1:8
+            x_flip[i, case] = -x[8 + i, case]
+            x_flip[8 + i, case] = -x[i, case]
+        end
+        x_flip[17, case] = -x[17, case]
+        x_flip[19, case] = -x[19, case]
+        x_flip[24, case] = x[25, case]
+        x_flip[25, case] = x[24, case]
+    end
     return x_flip
 end
 
@@ -267,30 +371,58 @@ end
 
 Flip outputs back after evaluating with flipped inputs.
 """
-function flip_outputs(y::AbstractMatrix{T}) where T
-    y_flip = copy(y)
-    N = (size(y, 1) - 6) ÷ 6
-
-    # CL (row 2) and CM (row 4) flip sign
-    y_flip[2, :] .= -y[2, :]
-    y_flip[4, :] .= -y[4, :]
-
-    # Transition locations swap (rows 5, 6)
-    y_flip[5, :] .= y[6, :]
-    y_flip[6, :] .= y[5, :]
-
-    if N > 0
-        # Swap upper/lower theta+H blocks (rows 7:6+2N ↔ 7+3N:6+5N)
-        y_flip[7:(6+2N), :] .= y[(7+3N):(6+5N), :]
-        y_flip[(7+3N):(6+5N), :] .= y[7:(6+2N), :]
-        # Swap upper/lower ue/vinf blocks with sign flip (velocity mirrors)
-        y_flip[(7+2N):(6+3N), :] .= -1 .* y[(7+5N):(6+6N), :]
-        y_flip[(7+5N):(6+6N), :] .= -1 .* y[(7+2N):(6+3N), :]
+function flip_outputs(y::AbstractMatrix)
+    n_stations = (size(y, 1) - 6) ÷ 6
+    y_flip = similar(y)
+    for row in axes(y, 1)
+        source, sign = flipped_row(row, n_stations)
+        @views y_flip[row, :] .= sign .* y[source, :]
     end
-
     return y_flip
 end
 
+"""
+    flipped_row(row, n_stations) -> (source, sign)
+
+Which row of a network output a top/bottom-flipped case reads `row` from, and with which
+sign. The single place the mirror symmetry of the output layout is written down, shared
+by [`flip_outputs`](@ref) and [`fuse_flipped!`](@ref).
+
+`CL` and `CM` change sign, the two transition locations swap, and of the six
+`n_stations`-long boundary-layer blocks the upper and lower halves swap — the momentum
+thickness and shape factor as they are, the edge velocity mirrored, so it changes sign.
+"""
+function flipped_row(row::Int, n_stations::Int)
+    (row == 2 || row == 4) && return (row, -1)
+    row == 5 && return (6, 1)
+    row == 6 && return (5, 1)
+    (row <= 6 || n_stations == 0) && return (row, 1)
+    block = (row - 7) ÷ n_stations
+    block in (0, 1) && return (row + 3n_stations, 1)
+    block == 2 && return (row + 3n_stations, -1)
+    block in (3, 4) && return (row - 3n_stations, 1)
+    block == 5 && return (row - 3n_stations, -1)
+    return (row, 1)
+end
+
+"""
+    fuse_flipped!(y, y_flip) -> y
+
+Average a direct output with its flipped counterpart in place, undoing the flip on the
+way ([`flipped_row`](@ref)). Reading `y_flip` while writing `y` is what lets the fusion
+land in storage that is already there instead of a third matrix.
+"""
+function fuse_flipped!(y::AbstractMatrix, y_flip::AbstractMatrix)
+    n_stations = (size(y, 1) - 6) ÷ 6
+    @inbounds for row in axes(y, 1)
+        source, sign = flipped_row(row, n_stations)
+        for case in axes(y, 2)
+            y[row, case] = (y[row, case] + sign * y_flip[source, case]) / 2
+        end
+    end
+    return y
+end
+
 """
     compute_optimal_x_points(n) -> Vector{Float64}
 
@@ -352,10 +484,25 @@ Turn a fused network output matrix into the integrated coefficients per case, un
 NeuralFoil's output scaling. The single place that scaling is written down.
 """
 function decode_coefficients(y::AbstractMatrix)
-    return (Vector{Float64}(y[2, :] ./ 2),
-            Vector{Float64}(clamp.(exp.((y[3, :] .- 2) .* 2), 0.0, 1.0)),
-            Vector{Float64}(y[4, :] ./ 20),
-            Vector{Float64}(sigmoid.(y[1, :])))
+    n_cases = size(y, 2)
+    return decode_coefficients!(zeros(n_cases), zeros(n_cases), zeros(n_cases),
+                                zeros(n_cases), y)
+end
+
+"""
+    decode_coefficients!(cl, cd, cm, confidence, y) -> (cl, cd, cm, confidence)
+
+[`decode_coefficients`](@ref) into vectors the caller owns, so a live polar reads a
+refresh out of the network without allocating four vectors for it.
+"""
+function decode_coefficients!(cl, cd, cm, confidence, y::AbstractMatrix)
+    @inbounds for case in axes(y, 2)
+        cl[case] = y[2, case] / 2
+        cd[case] = clamp(exp((y[3, case] - 2) * 2), 0.0, 1.0)
+        cm[case] = y[4, case] / 20
+        confidence[case] = sigmoid(y[1, case])
+    end
+    return (cl, cd, cm, confidence)
 end
 
 """
@@ -381,15 +528,29 @@ Symmetry-fused forward pass over a prepared input matrix, see
 [`neuralfoil_fused_output`](@ref). Takes the inputs already built so a caller that
 assembles its own batch does not go back through [`prepare_inputs`](@ref).
 """
-function fused_output(x::AbstractMatrix, model::NeuralFoilModel)
-    y = nn_forward(x, model)
-    y[1, :] .-= squared_mahalanobis_distance(x, model) ./ (2 * model.n_inputs)
+fused_output(x::AbstractMatrix, model::NeuralFoilModel) =
+    Matrix(fused_output!(NeuralFoilWorkspace(model, size(x, 2)), x, model))
 
-    x_flip = flip_inputs(x)
-    y_flip = nn_forward(x_flip, model)
-    y_flip[1, :] .-= squared_mahalanobis_distance(x_flip, model) ./ (2 * model.n_inputs)
+"""
+    fused_output!(work, x, model) -> AbstractMatrix
 
-    return (y .+ flip_outputs(y_flip)) ./ 2
+[`fused_output`](@ref) run entirely inside `work`, allocating nothing for a batch the
+workspace has room for. The result is a view of the workspace's own storage and stays
+valid until the next pass through it.
+"""
+function fused_output!(work::NeuralFoilWorkspace, x::AbstractMatrix,
+                       model::NeuralFoilModel)
+    n_cases = size(x, 2)
+    n_cases <= case_capacity(work) || throw(ArgumentError(
+        "A NeuralFoilWorkspace sized for $(case_capacity(work)) cases was handed " *
+        "$n_cases."))
+    x_flip = view(work.inputs_flip, :, 1:n_cases)
+    flip_inputs!(x_flip, x)
+    y = nn_forward!(work.direct, x, model, n_cases)
+    penalize_confidence!(y, x, model, work.centered)
+    y_flip = nn_forward!(work.flipped, x_flip, model, n_cases)
+    penalize_confidence!(y_flip, x_flip, model, work.centered)
+    return fuse_flipped!(y, y_flip)
 end
 
 """
@@ -431,10 +592,44 @@ fixed stations. Each `ue` matrix is `N × n_cases`. Interpolate these rather tha
 `Cp` is quadratic.
 """
 function decode_surface_velocity(y::AbstractMatrix)
-    N = (size(y, 1) - 6) ÷ 6
-    return (compute_optimal_x_points(N),
-            Matrix{Float64}(y[(7 + 2N):(6 + 3N), :]),
-            Matrix{Float64}(y[(7 + 5N):(6 + 6N), :]))
+    n_stations, upper, lower = surface_velocity_rows(y)
+    return (compute_optimal_x_points(n_stations),
+            Matrix{Float64}(y[upper, :]), Matrix{Float64}(y[lower, :]))
+end
+
+"""
+    decode_surface_velocity!(ue_upper, ue_lower, y) -> (ue_upper, ue_lower)
+
+[`decode_surface_velocity`](@ref) into matrices the caller owns. The station axis is
+fixed by the network and is carried by [`NeuralFoilWorkspace`](@ref) instead, so a live
+pressure refresh reads a pass out without allocating anything for it.
+"""
+function decode_surface_velocity!(ue_upper::AbstractMatrix, ue_lower::AbstractMatrix,
+                                  y::AbstractMatrix)
+    _, upper, lower = surface_velocity_rows(y)
+    @inbounds for case in axes(y, 2)
+        for (k, row) in enumerate(upper)
+            ue_upper[k, case] = y[row, case]
+        end
+        for (k, row) in enumerate(lower)
+            ue_lower[k, case] = y[row, case]
+        end
+    end
+    return (ue_upper, ue_lower)
+end
+
+"""
+    surface_velocity_rows(y) -> (n_stations, upper, lower)
+
+The number of boundary-layer stations a fused output matrix carries and the row ranges
+its upper and lower edge-velocity ratios sit in. Viewing those rows is how a caller that
+already holds the output reads the velocities without copying them out
+([`decode_surface_velocity`](@ref) is the copying form).
+"""
+function surface_velocity_rows(y::AbstractMatrix)
+    n_stations = (size(y, 1) - 6) ÷ 6
+    return (n_stations, (7 + 2n_stations):(6 + 3n_stations),
+            (7 + 5n_stations):(6 + 6n_stations))
 end
 
 """
diff --git a/src/panel.jl b/src/panel.jl
index 3565ce07..8ddcd826 100644
--- a/src/panel.jl
+++ b/src/panel.jl
@@ -297,19 +297,21 @@ Rewrite a panel's polar table with values at `alphas` [rad], ascending, and rebu
 interpolations. The table keeps the shape the panel was built with: a `POLAR_VECTORS`
 panel takes the angles as they are, and a `POLAR_MATRICES` panel takes them at every
 `delta` it already spans, since a regenerated polar carries its deflection as shape
-rather than as a flap angle and so says the same thing at each. The panel's `alpha_ref` and `alpha_window` are taken from the
-angles, so a table covering only a window around one angle is held at its ends rather than
-extrapolated past them (see [`window_alpha`](@ref)).
+rather than as a flap angle and so says the same thing at each. The panel's `alpha_ref`
+and `alpha_window` are taken from the angles, so a table covering only a window around one
+angle is held at its ends rather than extrapolated past them (see [`window_alpha`](@ref)).
 
 Written in place: the knots and the three value vectors reuse the panel's own
 `alpha_knots` and `cl_coeffs`/`cd_coeffs`/`cm_coeffs` storage whenever the sample count is
-unchanged, which is what lets a live polar source refresh every solve without growing the
-panel. The interpolations themselves are rebuilt, since `Interpolations` copies the values
-it is handed; they keep the extrapolation the panel was built with.
+unchanged ([`polar_column!`](@ref)), and the interpolations read that storage rather than
+a copy of it, so writing it is the whole update and the interpolation objects stay put
+([`refresh_polar!`](@ref)). Only a table that changes shape is rebuilt. Values may be
+views into a caller's own buffers; nothing here holds on to them.
 
 `shape` is the [`KulfanParameters`](@ref) the values were generated from, stored on the
 panel as `live_shape` so a panel's polar and the shape behind it are set together and
-cannot drift apart.
+cannot drift apart. The panel holds the shape itself, not a copy, so a source that
+rewrites one in place has already updated every panel flying it.
 """
 function set_polar!(panel::Panel, alphas, cl, cd, cm; shape=nothing)
     length(alphas) == length(cl) == length(cd) == length(cm) ||
@@ -324,40 +326,68 @@ function set_polar!(panel::Panel, alphas, cl, cd, cm; shape=nothing)
     panel.alpha_ref = (alphas[1] + alphas[end]) / 2
     panel.alpha_window = (alphas[end] - alphas[1]) / 2
     panel.live_shape = shape
-    for (dst_sym, src) in ((:alpha_knots, alphas), (:cl_coeffs, cl),
-                           (:cd_coeffs, cd), (:cm_coeffs, cm))
-        dst = getfield(panel, dst_sym)
-        if length(dst) == length(src)
-            dst .= src
-        else
-            setfield!(panel, dst_sym, collect(Float64, src))
-        end
-    end
-    panel.cl_interp = rebuild_polar(panel.cl_interp, panel.alpha_knots, panel.cl_coeffs)
-    panel.cd_interp = rebuild_polar(panel.cd_interp, panel.alpha_knots, panel.cd_coeffs)
-    panel.cm_interp = rebuild_polar(panel.cm_interp, panel.alpha_knots, panel.cm_coeffs)
+    panel.alpha_knots = polar_column!(panel.alpha_knots, alphas)
+    panel.cl_coeffs = polar_column!(panel.cl_coeffs, cl)
+    panel.cd_coeffs = polar_column!(panel.cd_coeffs, cd)
+    panel.cm_coeffs = polar_column!(panel.cm_coeffs, cm)
+    cl_interp = refresh_polar!(panel.cl_interp, panel.alpha_knots, panel.cl_coeffs)
+    cd_interp = refresh_polar!(panel.cd_interp, panel.alpha_knots, panel.cd_coeffs)
+    cm_interp = refresh_polar!(panel.cm_interp, panel.alpha_knots, panel.cm_coeffs)
+    cl_interp === panel.cl_interp || (panel.cl_interp = cl_interp)
+    cd_interp === panel.cd_interp || (panel.cd_interp = cd_interp)
+    cm_interp === panel.cm_interp || (panel.cm_interp = cm_interp)
     return nothing
 end
 
 """
-    rebuild_polar(old, knots, values) -> Extrapolation
+    polar_column!(stored, values) -> Vector{Float64}
+
+One column of a rewritten polar table in the panel's own `stored` vector, which is
+returned, or in a fresh vector when the sample count has changed — the one case a
+rewritten table has to grow, and the one that costs the panel's interpolations a rebuild.
+"""
+function polar_column!(stored::Vector{Float64}, values)
+    length(stored) == length(values) || return collect(Float64, values)
+    stored .= values
+    return stored
+end
+
+"""
+    refresh_polar!(old, knots, values) -> Extrapolation
 
-A 1D linear interpolation over `knots`/`values` carrying `old`'s extrapolation and knot
-container, so the rebuilt object has the type the panel's field was parameterised with.
-Only the values need rebuilding — `Interpolations` holds the knots by reference — but
-they are what it copies, so the object is rebuilt until it can take them in place.
+An interpolation reading `knots` and `values` themselves, carrying `old`'s extrapolation
+and knot container so it keeps the type the panel's field was parameterised with.
+`interpolate!` takes both arrays by reference, so once an interpolation reads a panel's
+own storage, writing that storage is the whole update and there is nothing to rebuild.
+A table that changes shape is rebuilt once, and reads in place from then on. The two
+dimensional form spreads the angles over every `delta` the panel spans.
 """
-rebuild_polar(old::Interpolations.Extrapolation{<:Any, 1}, knots, values) =
-    linear_interpolation(same_knots(old.itp.knots[1], knots), values;
-                         extrapolation_bc=old.et)
+function refresh_polar!(old::E, knots, values
+                       ) where {E <: Interpolations.Extrapolation{<:Any, 1}}
+    reads_from(old.itp.knots[1], knots) && old.itp.coefs === values && return old
+    return extrapolate(interpolate!((same_knots(old.itp.knots[1], knots),), values,
+                                    Gridded(Linear())), old.et)::E
+end
 
-function rebuild_polar(old::Interpolations.Extrapolation{<:Any, 2}, knots, values)
+function refresh_polar!(old::E, knots, values
+                       ) where {E <: Interpolations.Extrapolation{<:Any, 2}}
     deltas = old.itp.knots[2]
-    return linear_interpolation((same_knots(old.itp.knots[1], knots), deltas),
-                                repeat(values, 1, length(deltas));
-                                extrapolation_bc=old.et)
+    coefs = old.itp.coefs
+    if reads_from(old.itp.knots[1], knots) && size(coefs) == (length(knots), length(deltas))
+        for column in axes(coefs, 2)
+            @views coefs[:, column] .= values
+        end
+        return old
+    end
+    return extrapolate(interpolate!((same_knots(old.itp.knots[1], knots), deltas),
+                                    repeat(values, 1, length(deltas)),
+                                    Gridded(Linear())), old.et)::E
 end
 
+"Whether an interpolation's knot container is a view of `knots` rather than a copy."
+reads_from(container::ScanKnots, knots) = container.data === knots
+reads_from(container::AbstractVector, knots) = container === knots
+
 "The angles in the container the panel's interpolations were parameterised with."
 same_knots(::ScanKnots, knots) = ScanKnots(knots)
 same_knots(::AbstractVector, knots) = knots
diff --git a/src/section_aero.jl b/src/section_aero.jl
index 30dee3fd..eb84ebb0 100644
--- a/src/section_aero.jl
+++ b/src/section_aero.jl
@@ -9,13 +9,18 @@ plot or a traction pattern reads without knowing how the shape was produced. Fit
 with `AirfoilAero.fit_kulfan_parameters`, deform one with `AirfoilAero.deform_kulfan`
 and turn one into coordinates with `AirfoilAero.kulfan_to_coordinates`.
 
+Mutable, so a live polar source can rewrite one shape every solve instead of building a
+new one, and the panel pointing at it follows without being told
+(`AirfoilAero.deform_kulfan!`). That also makes a panel's `live_shape` the very object it
+was sampled from rather than a copy that compares equal to it.
+
 # Fields
 - `upper_weights::Vector{Float64}`: weights for upper surface
 - `lower_weights::Vector{Float64}`: weights for lower surface
 - `leading_edge_weight::Float64`: Leading edge modification weight
 - `TE_thickness::Float64`: Trailing edge thickness
 """
-struct KulfanParameters
+mutable struct KulfanParameters
     upper_weights::Vector{Float64}
     lower_weights::Vector{Float64}
     leading_edge_weight::Float64
diff --git a/test/airfoil_aero/test_live_polar.jl b/test/airfoil_aero/test_live_polar.jl
index 58afab15..9a799ec4 100644
--- a/test/airfoil_aero/test_live_polar.jl
+++ b/test/airfoil_aero/test_live_polar.jl
@@ -19,6 +19,22 @@ function polar_panels(n_panels; n_samples=5)
     end
 end
 
+"""One `refresh_live_polars!` at the reference condition, its confidence dropped."""
+function refresh_once!(live, panels, deflection)
+    refresh_live_polars!(live, panels, deg2rad(6.0), 3e6; deflection)
+    return nothing
+end
+
+"""
+Allocations of one `refresh_once!`, warmed up first. The confidence the refresh returns
+has to be dropped rather than measured: on Julia 1.11 `@allocated` boxes a `Float64`
+result, which reads as 16 bytes the refresh itself never allocated.
+"""
+function refresh_allocs(live, panels; deflection=nothing)
+    refresh_once!(live, panels, deflection)
+    return @allocated refresh_once!(live, panels, deflection)
+end
+
 @testset "Kulfan deformation" begin
     basis = KulfanBasis()
     base = KulfanParameters(fill(0.15, 8), fill(-0.05, 8), 0.0, 0.0)
@@ -121,6 +137,18 @@ end
                        [0.03, 0.02, 0.02, 0.03, 0.06], zeros(5))
     @test panel.alpha_knots === knots && panel.cl_coeffs === coeffs
 
+    # A refresh of the same shape is a write, not a rebuild: the interpolations read
+    # the panel's own storage, so the objects survive and nothing is allocated.
+    set_polar!(panel, alphas, [0.1, 0.4, 0.7, 1.0, 0.9], fill(0.02, 5), fill(-0.1, 5))
+    cl_interp = panel.cl_interp
+    set_polar!(panel, alphas, [0.2, 0.5, 0.8, 1.1, 1.0], fill(0.03, 5), fill(-0.2, 5))
+    @test panel.cl_interp === cl_interp
+    @test calculate_cl(panel, alphas[2]) ≈ 0.5
+    cl_next, cd_next, cm_next = [0.2, 0.5, 0.8, 1.1, 1.0], fill(0.03, 5), fill(-0.2, 5)
+    refresh(p, a, l, d, m) = set_polar!(p, a, l, d, m)
+    refresh(panel, alphas, cl_next, cd_next, cm_next)
+    @test (@allocated refresh(panel, alphas, cl_next, cd_next, cm_next)) == 0
+
     @test_throws ArgumentError set_polar!(panel, alphas, [0.0], [0.0], [0.0])
     @test_throws ArgumentError set_polar!(panel, reverse(alphas), zeros(5),
                                           zeros(5), zeros(5))
@@ -214,6 +242,13 @@ end
     @test panels[1].live_shape === live.deformed[1]
     @test panels[1].live_shape.upper_weights ≈ live.base[1].upper_weights
 
+    # Which makes a refresh the forward pass and nothing else. The shapes, the network
+    # inputs, both symmetries of the pass and the decoded coefficients all land in
+    # storage `LivePolars` already holds, so a solve loop refreshing every step adds
+    # nothing to the heap — the property the whole struct is shaped around.
+    @test refresh_allocs(live, panels; deflection=fill(camber, n_panels)) == 0
+    @test refresh_allocs(live, panels) == 0
+
     @test_throws ArgumentError refresh_live_polars!(live, panels[1:2], 0.0, 3e6)
     @test_throws ArgumentError LivePolars(fill(base, 2);
         settings=LivePolarSettings(; offsets=deg2rad.([1.0, 2.0, 3.0])))
@@ -311,7 +346,11 @@ end
     contour_x = [1.0, 0.5, 0.0, 0.5, 1.0]
     cf = [zeros(5), zeros(5)]
     same = cf[1]
-    live_surface_friction!(cf, [contour_x, contour_x], [3e6, 1e6])
+    contours, reynolds = [contour_x, contour_x], [3e6, 1e6]
+    live_surface_friction!(cf, contours, reynolds)
+    friction(c, x, re) = live_surface_friction!(c, x, re)
+    friction(cf, contours, reynolds)
+    @test (@allocated friction(cf, contours, reynolds)) == 0
     @test cf[1] === same
     @test cf[1] ≈ [AirfoilAero.flat_plate_cf(x, 3e6) for x in contour_x]
     # Lower Reynolds is more friction, and it is highest at the nose.