Live polars: regenerate a panel's polar from its deformed shape

.gitignore

@@ -59,3 +59,6 @@ lib/**/Manifest.toml
 
 # Generated test geometries (regenerated on demand by ram_air_matrix_wing)
 test/generated/
+
+# Generated aero tables and run logs
+*.arrow

docs/src/functions.md

@@ -41,6 +41,19 @@ deform_section
 analyze_section
 analyze_sweep
 neuralfoil_aero
+deform_kulfan
+chord_residual
+chord_line
+control_point_deflection
+panel_kulfan_parameters
+refresh_live_polars!
+deform_live_shapes!
+apply_live_shapes!
+polar_drift
+refresh_live_pressure!
+live_surface_friction!
+contour_shape_matrix
+live_shape_offset!
 generate_aero_matrices
 generate_polar_from_coordinates
 generate_polar_from_dat

docs/src/private_functions.md

@@ -98,6 +98,10 @@ update_non_deformed_sections!
 ### Aerodynamic data and Cp
 ```@docs
 calculate_new_aero_data
+set_sampled_polar!
+decode_surface_pressure
+contour_pressure
+sampled_value
 assemble_polar_matrix
 load_matrix_polar_data
 read_aero_matrix
@@ -153,8 +157,11 @@ smooth_turning!
 load_neuralfoil_model
 neuralfoil_section
 neuralfoil_fused_output
+fused_output
+decode_coefficients
 nn_forward
 prepare_inputs
+fill_case_input!
 flip_inputs
 flip_outputs
 squared_mahalanobis_distance
@@ -221,6 +228,7 @@ map_airfoil_3d
 fitted_airfoil_3d
 generated_slices
 airfoil_skin_geometry
+panel_contour
 panel_normal
 plate_hinge_local
 panel_plate_geometry

docs/src/private_types.md

@@ -7,6 +7,7 @@ CurrentModule = VortexStepMethod
 ### Wing Geometry, Panel and Aerodynamics
 ```@docs
 Panel
+KulfanParameters
 PanelProperties
 Filament
 BoundFilament
@@ -23,7 +24,9 @@ LEI_AIRFOIL_BREUKELS
 CurrentModule = VortexStepMethod.AirfoilAero
 ```
 ```@docs
-KulfanParameters
+KulfanBasis
+LivePolarSettings
+LivePolars
 NeuralFoilModel
 NeuralFoilResult
 ```

ext/VortexStepMethodMakieExt.jl

@@ -146,15 +146,34 @@ function Makie.plot!(ax, panel::VortexStepMethod.Panel; color=(:red, 0.2), R_b_w
     return plots
 end
 
+"""
+    panel_contour(panel, section) -> (x, y)
+
+The airfoil contour to draw a panel with: coordinates of its `live_shape` when a live
+polar source has put one there, else the section's tabulated contour at the panel's
+`delta`. The live branch reads the stored shape object itself, so what is drawn is what
+was flown.
+"""
+function panel_contour(panel, section)
+    shape = panel.live_shape
+    isnothing(shape) || return VortexStepMethod.AirfoilAero.kulfan_to_coordinates(shape)
+    return VortexStepMethod.section_surface(section.section_aero, 0.0, panel.delta)[1:2]
+end
+
 """
     airfoil_skin_geometry(body; R_b_w=nothing, T_b_w=nothing) -> (vertices, faces, ribs)
 
-Lofted airfoil skin of a `BodyAerodynamics`: each section's deflected contour
-(`section_surface` at the panel's `delta`) is fitted between the panel's `corner_points`
-by a 2D similarity, TE pinned so a deflection bulges the fore body up. The skin reflects
-`delta` only when the geometry carries per-`delta` slices (`obj_to_yaml` with a
-`delta_range`); with δ=0-only data it renders undeflected, a deliberate cue that the
-deflected slices are missing. Transformed to world by `R_b_w`/`T_b_w`.
+Lofted airfoil skin of a `BodyAerodynamics`: each section's contour is fitted between
+the panel's `corner_points` by a 2D similarity, TE pinned so a deflection bulges the
+fore body up.
+
+The contour is the panel's own `live_shape` when it has one — the very
+[`KulfanParameters`](@ref) object the live polar source deformed and handed to the
+airfoil solver, not a re-derivation of it, so a deformation bug shows up in the picture
+instead of being papered over. Otherwise it is `section_surface` at the panel's `delta`,
+which reflects `delta` only when the geometry carries per-`delta` slices (`obj_to_yaml`
+with a `delta_range`); with δ=0-only data it renders undeflected, a deliberate cue that
+the deflected slices are missing. Transformed to world by `R_b_w`/`T_b_w`.
 `vertices`/`faces` triangulate the skin between consecutive equal-node sections; `ribs` is
 one closed contour polyline per section. Sections without contour data are skipped.
 """
@@ -169,10 +188,10 @@ function airfoil_skin_geometry(body; R_b_w=nothing, T_b_w=nothing)
         n_panels = n - 1
         n_panels < 1 && continue
         for (i, section) in enumerate(sections)
-            isnothing(section.section_aero) && continue
             panel_idx = panel_offset + min(i, n_panels)
             panel_idx <= length(body.panels) || continue
             panel = body.panels[panel_idx]
+            (isnothing(section.section_aero) && isnothing(panel.live_shape)) && continue
             corners = panel.corner_points
             corner1 = Point3f(corners[:, 1]); corner3 = Point3f(corners[:, 3])
             corner2 = Point3f(corners[:, 2]); corner4 = Point3f(corners[:, 4])
@@ -183,8 +202,7 @@ function airfoil_skin_geometry(body; R_b_w=nothing, T_b_w=nothing)
             chord_len = norm(chord)
             chord_len < 1e-9 && continue
             up = panel_normal(panel)
-            xs, ys, _, _ = VortexStepMethod.section_surface(section.section_aero,
-                                                            0.0, panel.delta)
+            xs, ys = panel_contour(panel, section)
             le_i = argmin(xs)
             le_x = xs[le_i]; le_y = ys[le_i]
             te_x = 0.5 * (xs[1] + xs[end]); te_y = 0.5 * (ys[1] + ys[end])

src/VortexStepMethod.jl

@@ -36,7 +36,8 @@ export calculate_projected_area, calculate_span
 export MVec3
 
 export LLT, Model, VSM
-export AeroModel, INVISCID, POLY, LEI_AIRFOIL_BREUKELS, POLAR_MATRICES, POLAR_VECTORS
+export AeroModel, INVISCID, POLY, LEI_AIRFOIL_BREUKELS, POLAR_MATRICES, POLAR_VECTORS, SAMPLED
+export KulfanParameters
 export BILLOWING, COSINE, LINEAR, PanelDistribution, SPLIT_PROVIDED, UNCHANGED
 export ELLIPTIC, InitialGammaDistribution, ZEROS
 export FAILURE, FEASIBLE, INFEASIBLE, SolverStatus
@@ -237,7 +238,7 @@ Enumeration of the implemented wing types.
 @enum WingType  RECTANGULAR CURVED ELLIPTICAL
 
 """
-   AeroModel `POLY` `POLAR_VECTORS` `POLAR_MATRICES` `INVISCID`
+   AeroModel `POLY` `POLAR_VECTORS` `POLAR_MATRICES` `INVISCID` `SAMPLED`
 
 Enumeration of the implemented aerodynamic models. See also: [AeroData](@ref)
 
@@ -247,6 +248,11 @@ Enumeration of the implemented aerodynamic models. See also: [AeroData](@ref)
 - `POLAR_VECTORS`: Polar vectors as function of alpha (lookup tables with interpolation)
 - `POLAR_MATRICES`: Polar matrices as function of alpha and delta (lookup tables with interpolation)
 - INVISCID
+- `SAMPLED`: cl/cd/cm sampled at ascending angles of attack per panel, interpolated
+  between them and held flat past either end. This is what a live polar source writes
+  each solve, see [`refresh_live_polars!`](@ref
+  VortexStepMethod.AirfoilAero.refresh_live_polars!). Samples are the polar, so a stall
+  knee inside the sampled range is represented rather than smoothed. Ignores `delta`.
 
 `LEI_AIRFOIL_BREUKELS` is a deprecated alias of `POLY`.
 
@@ -257,6 +263,7 @@ where `alpha` is the angle of attack, `delta` is trailing edge angle.
    POLAR_VECTORS
    POLAR_MATRICES
    INVISCID
+   SAMPLED
 end
 
 """
@@ -339,6 +346,9 @@ Union of different definitions of the aerodynamic properties of a wing section.
   - (`alpha_range`, `cl_vector`, `cd_vector`, `cm_vector`) for `POLAR_VECTORS`
   - (`alpha_range`, `delta_range`, `cl_matrix`, `cd_matrix`, `cm_matrix`) for `POLAR_MATRICES`
 
+`SAMPLED` carries no section-level data: it is written onto a panel at run time by a
+live polar source, never read from a section.
+
 where `alpha` is the angle of attack [rad], `delta` is trailing edge angle [rad], `cl` the lift coefficient,
 `cd` the drag coefficient and `cm` the pitching moment coefficient. The camber of a kite refers to 
 the curvature of its airfoil shape. The camber is typically measured as the maximum distance 

src/airfoil_aero/AirfoilAero.jl

@@ -8,11 +8,14 @@ using NPZ
 using Xfoil
 using Printf: @sprintf
 using ..VortexStepMethod: SectionAero, interpolate_matrix_nans!, delta_suffix,
-                          write_node_rows
+                          write_node_rows, section_surface, set_sampled_polar!,
+                          KulfanParameters
 
 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")
@@ -29,6 +32,12 @@ export ShrinkWrap, shrink_wrap
 export fit_kulfan_parameters, kulfan_to_coordinates
 export NeuralFoilModel, NeuralFoilResult, load_neuralfoil_model
 export neuralfoil_aero, neuralfoil_section
+export KulfanBasis, deform_kulfan, control_point_deflection
+export chord_residual, chord_line
+export LivePolarSettings, LivePolars, panel_kulfan_parameters
+export refresh_live_polars!, refresh_live_pressure!, live_surface_friction!
+export contour_shape_matrix, live_shape_offset!
+export polar_drift
 export AbstractAirfoilSolver, XFoilSolver, NeuralFoilSolver
 export SectionSolution, DeformedSection, deform_section, analyze_section, analyze_sweep
 export create_2d_polars, generate_aero_matrices, generate_section_aero, lei_poly_coeffs

src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl

@@ -47,9 +47,8 @@ Assemble a full-contour [`SectionSolution`](@ref) for case `i` of a NeuralFoil s
 `(x, y)` split at the leading edge `le`, carrying the precomputed `cf`.
 """
 function neuralfoil_contour_solution(alpha, res, i, x, y, le, cf)
-    up = linear_interpolation(res.x, res.cp_upper[:, i]; extrapolation_bc=Line())
-    lo = linear_interpolation(res.x, res.cp_lower[:, i]; extrapolation_bc=Line())
-    cp = [(k <= le ? up : lo)(clamp(x[k], 0.0, 1.0)) for k in eachindex(x)]
+    cp = contour_pressure(res.x, view(res.cp_upper, :, i),
+                          view(res.cp_lower, :, i), x, le)
     return SectionSolution(alpha, res.cl[i], res.cd[i], res.cm[i],
                            res.confidence[i], x, y, cp, cf)
 end

src/airfoil_aero/deform.jl

@@ -0,0 +1,145 @@
+"""
+    KulfanBasis(; n_stations=60, n_weights=8, ridge=1e-3)
+
+The fixed CST basis a shape deformation is projected onto: chord stations `x` in
+`[0, 1]` and the ridge-regularised inverse `projection` of the class-times-Bernstein
+matrix `C(x)·B(x)` those stations span.
+
+`ridge` is the Tikhonov weight, relative to the basis' own largest singular value, that
+keeps the projection from answering a deflection it cannot represent with weights far
+larger than the airfoil they correct. A plain pseudoinverse has no such bound: a kinked
+deflection — a strut buckling is one — lands on the basis' weakest directions and comes
+back amplified a hundredfold and alternating in sign, which is not an airfoil. The ridge
+trades a small, measured under-response on deflections the basis *can* hold for a bounded
+answer on the ones it cannot.
+
+CST is linear in its weights, so a surface displacement is a matvec against this
+constant matrix — never a refit. Refitting inside a loop is not an option: the Kulfan
+fit is non-unique, so the same shape fitted twice returns weight vectors differing by
+more than the deformation signal, which reaches the polars as frame-to-frame jitter.
+
+`C(1) = 0`, so the parameterisation cannot put the trailing edge off the chord line.
+Deflections must therefore be measured **against the deformed chord**, with the chord
+rotation and stretch taken from the leading- and trailing-edge points; what is left is
+the representable residual.
+"""
+struct KulfanBasis
+    "Chord stations the deflection is sampled on, ascending in `[0, 1]`."
+    x::Vector{Float64}
+    "`n_weights × length(x)` regularised inverse of `C(x)·B(x)`, mapping deflection to weights."
+    projection::Matrix{Float64}
+    "Number of CST weights per surface, matching the airfoil being deformed."
+    n_weights::Int
+    "Tikhonov weight the projection was built with, relative to the basis' largest singular value."
+    ridge::Float64
+end
+
+function KulfanBasis(; n_stations::Int=60, n_weights::Int=8, ridge::Real=1e-3)
+    n_stations > n_weights || throw(ArgumentError(
+        "KulfanBasis needs more stations than weights; got $n_stations and $n_weights."))
+    ridge >= 0 || throw(ArgumentError("KulfanBasis ridge must not be negative."))
+    theta = range(0, pi, n_stations)
+    x = @. (1 - cos(theta)) / 2
+    shape = class_function(x) .* bernstein_basis(x, n_weights - 1)
+    scale = maximum(svdvals(shape))^2
+    projection = iszero(ridge) ? pinv(shape) :
+        (shape' * shape + ridge * scale * I) \ shape'
+    return KulfanBasis(collect(x), projection, n_weights, Float64(ridge))
+end
+
+"""
+    deform_kulfan(basis, base, upper_deflection, lower_deflection) -> KulfanParameters
+    deform_kulfan(basis, base, camber) -> KulfanParameters
+
+Add a surface deflection to a fixed set of Kulfan parameters, both deflections sampled
+on `basis.x` and normalized by chord. The three-argument form applies one camber
+deflection to both surfaces, leaving the thickness distribution untouched; the
+four-argument form deforms the surfaces independently, which is what a double-skin
+membrane with its own upper and lower control points needs.
+
+Each deflection is first reduced to the part the basis can carry, see
+[`chord_residual`](@ref): the straight line through its own endpoints is a chord
+rotation and translation, which the basis cannot express and `pinv` answers with
+runaway weights. Take that line from [`chord_line`](@ref) if the frame the deflection
+was measured in has not already absorbed it.
+
+The leading-edge weight and the trailing-edge thickness are carried over unchanged:
+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)
+    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)
+end
+
+"""
+    chord_residual(basis, deflection) -> Vector{Float64}
+    chord_line(basis, deflection) -> (offset, slope)
+
+The part of a deflection the CST basis can represent, and the part it cannot. `C(x)`
+vanishes at both chord ends, so the basis can put neither the leading nor the trailing
+edge off the chord line: a deflection that does not end at zero is asking for a chord
+**rotation and translation**, not a camber change. Projecting it anyway is not merely
+inexact, it is unstable — `pinv` answers an unrepresentable end displacement with
+weights an order of magnitude past the ones it is correcting, and the airfoil that
+comes back is not one.
+
+[`chord_residual`](@ref) removes the straight line through the deflection's own
+endpoints and returns what is left. [`chord_line`](@ref) returns that line as
+`(offset, slope)`, both over chord: `offset` moves the leading edge and `slope` is the
+chord rotation `atan(slope)` [rad] the caller owes its angle of attack, when the frame
+it measured the deflection in has not already absorbed it.
+"""
+function chord_residual(basis::KulfanBasis, deflection::AbstractVector)
+    offset, slope = chord_line(basis, deflection)
+    return deflection .- (offset .+ slope .* basis.x)
+end
+
+function chord_line(basis::KulfanBasis, deflection::AbstractVector)
+    length(deflection) == length(basis.x) || throw(ArgumentError(
+        "Deflection must be sampled on the basis' $(length(basis.x)) stations."))
+    span = basis.x[end] - basis.x[1]
+    abs(span) < eps() && return (deflection[1], 0.0)
+    slope = (deflection[end] - deflection[1]) / span
+    return (deflection[1] - slope * basis.x[1], slope)
+end
+
+deform_kulfan(basis::KulfanBasis, base::KulfanParameters,
+              camber::AbstractVector) = deform_kulfan(basis, base, camber, camber)
+
+"""
+    control_point_deflection(basis, fractions, deflections) -> Vector{Float64}
+
+Resample a deflection given at arbitrary chord `fractions` onto `basis.x`, by linear
+interpolation, held flat outside the sampled range. This is the generic bridge from
+control points — beam nodes now, membrane nodes later — to the CST basis: the caller
+only has to say where along the chord each point sits and how far off the chord line it
+moved, both normalized by the local chord.
+
+`fractions` need not be sorted, but must not repeat a station.
+"""
+function control_point_deflection(basis::KulfanBasis, fractions::AbstractVector,
+                                  deflections::AbstractVector)
+    length(fractions) == length(deflections) || throw(ArgumentError(
+        "control_point_deflection: $(length(fractions)) fractions for " *
+        "$(length(deflections)) deflections."))
+    order = sortperm(collect(float.(fractions)))
+    knots = collect(float.(fractions))[order]
+    values = collect(float.(deflections))[order]
+    if length(knots) < 2
+        return fill(isempty(values) ? 0.0 : values[1], length(basis.x))
+    end
+    allunique(knots) || throw(ArgumentError(
+        "control_point_deflection: repeated chord fraction in $knots."))
+    interp = linear_interpolation(knots, values; extrapolation_bc=Flat())
+    return [interp(xi) for xi in basis.x]
+end