From acbbb8fd6db4976fc564c775aad6c6f33bfbad3e Mon Sep 17 00:00:00 2001
From: Bart <bart@vandelint.net>
Date: Thu, 7 May 2026 17:12:39 +0200
Subject: [PATCH 1/8] Lipo deform

---
 src/obj_geometry.jl                           |   6 +-
 src/solver.jl                                 |   2 +
 src/wing_geometry.jl                          | 285 ++++++++++++------
 .../test_yaml_wing_deformation.jl             | 100 ++----
 4 files changed, 236 insertions(+), 157 deletions(-)

diff --git a/src/obj_geometry.jl b/src/obj_geometry.jl
index fdaefd2c..4835216c 100644
--- a/src/obj_geometry.jl
+++ b/src/obj_geometry.jl
@@ -309,7 +309,10 @@ function refine_obj_wing!(wing::AbstractWing{T}; recompute_mapping=true) where {
     end
 
     # 5. Compute panel mapping and update non_deformed_sections
-    recompute_mapping && VortexStepMethod.compute_refined_panel_mapping!(wing)
+    if recompute_mapping
+        VortexStepMethod.compute_refined_panel_mapping!(wing)
+        VortexStepMethod.compute_refined_section_interpolation!(wing)
+    end
     VortexStepMethod.update_non_deformed_sections!(wing)
 
     return nothing
@@ -584,6 +587,7 @@ function ObjWing(
             MVector{3, Float64}(spanwise_direction),
             sections, Section{Float64}[], remove_nan, use_prior_polar, 0.0,  # billowing_percentage
             Int16[],  # refined_panel_mapping empty
+            Int16[], Float64[],  # refined_section interpolation cache empty
             Section{Float64}[], zeros(Float64, n_panels), zeros(Float64, n_panels),
             Float64(mass), Float64(gamma_tip), Matrix{Float64}(inertia_tensor),
             MVector{3, Float64}(T_cad_body), MMatrix{3, 3, Float64, 9}(R_cad_body),
diff --git a/src/solver.jl b/src/solver.jl
index eff11eaf..ede7c33b 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -1015,6 +1015,8 @@ function _wing_with_eltype(wing::Wing{P, Float64}, ::Type{TD}) where {P, TD}
         wing.use_prior_polar,
         wing.billowing_percentage,
         copy(wing.refined_panel_mapping),
+        copy(wing.refined_section_left_idx),
+        Vector{TD}(wing.refined_section_weight),
         Section{TD}[_section_with_eltype(s, TD) for s in wing.non_deformed_sections],
         Vector{TD}(wing.theta_dist),
         Vector{TD}(wing.delta_dist),
diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl
index d184cd56..ce4f38db 100644
--- a/src/wing_geometry.jl
+++ b/src/wing_geometry.jl
@@ -246,6 +246,12 @@ mutable struct Wing{P, T} <: AbstractWing{T}
     # Grouping
     refined_panel_mapping::Vector{Int16}  # Maps each refined panel index to unrefined section index (1 to n_unrefined_sections)
 
+    # Linear interpolation cache from unrefined sections to refined sections.
+    # For refined section i: value = weight * unrefined[left_idx] + (1 - weight) * unrefined[left_idx + 1].
+    # weight == 1 at unrefined endpoints, so refined section values match unrefined endpoints exactly.
+    refined_section_left_idx::Vector{Int16}  # Length: n_panels + 1
+    refined_section_weight::Vector{T}        # Length: n_panels + 1
+
     # Deformation fields
     non_deformed_sections::Vector{Section{T}}
     theta_dist::Vector{T}  # Length: n_panels (panel twist angles)
@@ -308,6 +314,8 @@ function Wing(n_panels::Int;
         Section{Float64}[], Section{Float64}[], remove_nan, use_prior_polar, Float64(billowing_percentage),
         # Grouping
         Int16[],
+        # Refined-section interpolation cache
+        Int16[], Float64[],
         # Deformation fields
         Section{Float64}[], zeros(max(0, n_panels)), zeros(max(0, n_panels)),
         # Physical properties (defaults for non-OBJ wings)
@@ -341,89 +349,77 @@ function reinit!(wing::AbstractWing)
 end
 
 """
-    unrefined_deform!(wing::Wing, theta_angles=nothing, delta_angles=nothing; smooth=false)
+    unrefined_deform!(wing::Wing, theta_angles=nothing, delta_angles=nothing)
 
-Apply deformation angles directly to unrefined wing sections.
+Apply deformation angles defined per unrefined section.
 
-For wings that support deformation (OBJ-based wings with non_deformed_sections), this
-applies theta_angles and delta_angles directly to unrefined sections and then applies deformation.
-For wings without deformation support (YAML-based), this is a no-op that only succeeds
-if both angle inputs are nothing.
+Refined-section twist and TE-deflection values are computed by linear interpolation
+from the unrefined-section inputs using the precomputed `refined_section_left_idx` /
+`refined_section_weight` cache built at refinement time. Endpoint refined sections
+take the unrefined endpoint values exactly. The panel-level `theta_dist` /
+`delta_dist` arrays are then filled by averaging adjacent refined-section values, so
+downstream consumers (solver, body aerodynamics) see a per-panel value.
 
 # Arguments
-- `wing::Wing`: The wing to deform
-- `theta_angles::AbstractVector`: Twist angles in radians for each unrefined section (or nothing).
-  Length must be `n_unrefined_sections`
-- `delta_angles::AbstractVector`: Trailing edge deflection angles in radians for each unrefined section (or nothing).
-  Length must be `n_unrefined_sections`
-- `smooth::Bool`: DEPRECATED - no longer used. Apply smoothing to input angles if needed.
-
-# Algorithm
-1. Copies theta_angles and delta_angles directly to wing.theta_dist and wing.delta_dist
-2. Calls deform! to update wing geometry and propagate to refined sections
-
-# Errors
-- Throws `ArgumentError` if wing doesn't support deformation but angles are provided
-- Throws `ArgumentError` if angle vectors don't match n_unrefined_sections
-
-# Returns
-- `nothing` (modifies wing in-place)
+- `wing::Wing`: Wing to deform (must have non_deformed_sections).
+- `theta_angles::AbstractVector`: Twist angles in radians, one per unrefined section.
+  Pass `nothing` to leave twist unchanged.
+- `delta_angles::AbstractVector`: TE deflection angles in radians, one per unrefined
+  section. Pass `nothing` to leave deflection unchanged.
 """
 function unrefined_deform!(wing::Wing, theta_angles=nothing, delta_angles=nothing;
                            smooth=false, smooth_window=nothing)
-    # Check if deformation is supported
     can_deform = !isempty(wing.non_deformed_sections)
-
-    # If no deformation requested, just return
     isnothing(theta_angles) && isnothing(delta_angles) && return nothing
 
-    # If deformation requested but not supported, throw error
     if !can_deform
         throw(ArgumentError("This Wing does not support deformation. Only OBJ-based wings created with ObjWing() can be deformed."))
     end
 
-    # Validate inputs
-    !isnothing(theta_angles) && length(theta_angles) != wing.n_unrefined_sections &&
-        throw(ArgumentError("theta_angles must have length n_unrefined_sections = $(wing.n_unrefined_sections), got $(length(theta_angles))"))
-    !isnothing(delta_angles) && length(delta_angles) != wing.n_unrefined_sections &&
-        throw(ArgumentError("delta_angles must have length n_unrefined_sections = $(wing.n_unrefined_sections), got $(length(delta_angles))"))
+    n_unref = wing.n_unrefined_sections
+    !isnothing(theta_angles) && length(theta_angles) != n_unref &&
+        throw(ArgumentError("theta_angles must have length n_unrefined_sections = $(n_unref), got $(length(theta_angles))"))
+    !isnothing(delta_angles) && length(delta_angles) != n_unref &&
+        throw(ArgumentError("delta_angles must have length n_unrefined_sections = $(n_unref), got $(length(delta_angles))"))
+
+    n_sec = wing.n_panels + 1
+    length(wing.refined_section_left_idx) == n_sec &&
+        length(wing.refined_section_weight) == n_sec ||
+        throw(ArgumentError(
+            "Refined-section interpolation cache is not initialized. Call refine!(wing) first."))
 
-    # Map unrefined sections → panels → sections (no smoothing yet)
     if !isnothing(theta_angles)
-        map_unrefined_to_sections!(wing.theta_dist, theta_angles, wing.refined_panel_mapping, wing.n_panels)
+        section_thetas = _interpolate_unrefined_to_refined(wing, theta_angles)
+        _apply_refined_section_thetas!(wing, section_thetas)
+        for i in 1:wing.n_panels
+            wing.theta_dist[i] = (section_thetas[i] + section_thetas[i + 1]) / 2
+        end
     end
     if !isnothing(delta_angles)
-        map_unrefined_to_sections!(wing.delta_dist, delta_angles, wing.refined_panel_mapping, wing.n_panels)
+        section_deltas = _interpolate_unrefined_to_refined(wing, delta_angles)
+        for i in 1:wing.n_panels
+            wing.delta_dist[i] = (section_deltas[i] + section_deltas[i + 1]) / 2
+        end
     end
-
-    # Apply deformation with optional smoothing
-    deform!(wing, smooth=smooth, smooth_window=smooth_window)
     return nothing
 end
 
 """
-    map_unrefined_to_sections!(panel_dist, unrefined_angles, panel_mapping, n_panels)
-
-Map angles from unrefined sections to panels.
-Steps: unrefined[1:n_unrefined] → panels[1:n_panels]
+    _interpolate_unrefined_to_refined(wing, unrefined_values) -> Vector
 
-# Arguments
-- `panel_dist::Vector{Float64}`: Output panel angles (length n_panels)
-- `unrefined_angles::Vector{Float64}`: Input unrefined section angles
-- `panel_mapping::Vector{Int16}`: Maps panel index to unrefined section index
-- `n_panels::Int`: Number of panels
-"""
-function map_unrefined_to_sections!(panel_dist,
-                                     unrefined_angles,
-                                     panel_mapping,
-                                     n_panels)
-    # Map unrefined sections to panels
-    for i in 1:n_panels
-        unrefined_idx = panel_mapping[i]
-        panel_dist[i] = unrefined_angles[unrefined_idx]
+Linearly interpolate `unrefined_values` (length `n_unrefined_sections`) to refined
+sections (length `n_panels + 1`) using the cached weights. Endpoints match exactly.
+"""
+function _interpolate_unrefined_to_refined(wing::Wing{P, T},
+                                            unrefined_values) where {P, T}
+    n_sec = wing.n_panels + 1
+    out = Vector{T}(undef, n_sec)
+    for i in 1:n_sec
+        l = wing.refined_section_left_idx[i]
+        wl = wing.refined_section_weight[i]
+        out[i] = wl * unrefined_values[l] + (one(T) - wl) * unrefined_values[l + 1]
     end
-
-    return nothing
+    return out
 end
 
 """
@@ -502,51 +498,70 @@ Updates wing.refined_sections based on wing.non_deformed_sections and stored dis
 function deform!(wing::Wing{P, T}; smooth=false, smooth_window=nothing) where {P, T}
     !isempty(wing.non_deformed_sections) || return nothing
 
-    # Apply smoothing if requested
     if smooth
-        # Default window size based on refinement ratio
         if isnothing(smooth_window)
             smooth_window = max(3, round(Int, wing.n_panels / wing.n_unrefined_sections))
         end
-
-        # Ensure window is odd for symmetric averaging
         if smooth_window % 2 == 0
             smooth_window += 1
         end
-
-        # Apply moving average to theta_dist and delta_dist (panel-level)
         smooth_distribution!(wing.theta_dist, smooth_window)
         smooth_distribution!(wing.delta_dist, smooth_window)
     end
 
+    n_sec = wing.n_panels + 1
+    section_thetas = Vector{T}(undef, n_sec)
+    _panel_thetas_to_section_thetas!(section_thetas, wing.theta_dist)
+    _apply_refined_section_thetas!(wing, section_thetas)
+    return nothing
+end
+
+"""
+    _panel_thetas_to_section_thetas!(section_thetas, panel_thetas)
+
+Convert panel-level twist angles (length `n`) to refined-section angles (length `n+1`).
+Interior sections take the average of adjacent panels. Boundary sections use linear
+extrapolation from the two nearest panels, so a linear input produces a linear output
+across the full refined-section range.
+"""
+function _panel_thetas_to_section_thetas!(section_thetas, panel_thetas)
+    n = length(panel_thetas)
+    n_sec = length(section_thetas)
+    n_sec == n + 1 || throw(ArgumentError(
+        "section_thetas length $(n_sec) must equal panel_thetas length + 1 = $(n + 1)"))
+    if n == 1
+        section_thetas[1] = panel_thetas[1]
+        section_thetas[2] = panel_thetas[1]
+        return nothing
+    end
+    section_thetas[1] = 1.5 * panel_thetas[1] - 0.5 * panel_thetas[2]
+    for i in 2:n
+        section_thetas[i] = (panel_thetas[i - 1] + panel_thetas[i]) / 2
+    end
+    section_thetas[n_sec] = 1.5 * panel_thetas[n] - 0.5 * panel_thetas[n - 1]
+    return nothing
+end
+
+"""
+    _apply_refined_section_thetas!(wing, section_thetas)
+
+Rotate each refined section's TE point around its LE point by the given per-section
+twist angle, using `wing.non_deformed_sections` as the reference geometry.
+"""
+function _apply_refined_section_thetas!(wing::Wing{P, T}, section_thetas) where {P, T}
     local_y = zeros(MVector{3, T})
     chord = zeros(MVector{3, T})
     normal = zeros(MVector{3, T})
 
-    # Process all refined sections (n_panels + 1)
-    # Convert panel angles to section angles by averaging
     for i in 1:(wing.n_panels + 1)
-        # Determine theta for this section by averaging adjacent panels
-        if i == 1
-            # First section: use panel 1 angle
-            theta = wing.theta_dist[1]
-        elseif i == wing.n_panels + 1
-            # Last section: use last panel angle
-            theta = wing.theta_dist[wing.n_panels]
-        else
-            # Middle sections: average of panels (i-1) and i
-            theta = (wing.theta_dist[i-1] + wing.theta_dist[i]) / 2.0
-        end
-
+        theta = section_thetas[i]
         section = wing.non_deformed_sections[i]
 
-        # Compute local coordinate system
         if i < wing.n_panels + 1
-            section2 = wing.non_deformed_sections[i+1]
+            section2 = wing.non_deformed_sections[i + 1]
             local_y .= normalize(section.LE_point - section2.LE_point)
         else
-            # For last section, use same local_y as previous
-            section_prev = wing.non_deformed_sections[i-1]
+            section_prev = wing.non_deformed_sections[i - 1]
             local_y .= normalize(section_prev.LE_point - section.LE_point)
         end
 
@@ -805,19 +820,25 @@ function refine!(wing::AbstractWing{T}; recompute_mapping=true, sort_sections=tr
         if isequal(wing.spanwise_distribution, UNCHANGED) ||
                length(wing.unrefined_sections) == n_sections
             copy_sections_to_refined!(wing; reuse_aero_data)
-            recompute_mapping && compute_refined_panel_mapping!(wing)
+            if recompute_mapping
+                compute_refined_panel_mapping!(wing)
+                compute_refined_section_interpolation!(wing)
+            end
             update_non_deformed_sections!(wing)
             return nothing
         else
             wing.refined_sections = Section{T}[Section{T}() for _ in 1:wing.n_panels+1]
         end
     end
-    
+
     # Handle special cases
     if isequal(wing.spanwise_distribution, UNCHANGED) ||
             length(wing.unrefined_sections) == n_sections
         copy_sections_to_refined!(wing; reuse_aero_data)
-        recompute_mapping && compute_refined_panel_mapping!(wing)
+        if recompute_mapping
+            compute_refined_panel_mapping!(wing)
+            compute_refined_section_interpolation!(wing)
+        end
         update_non_deformed_sections!(wing)
         return nothing
     end
@@ -834,7 +855,10 @@ function refine!(wing::AbstractWing{T}; recompute_mapping=true, sort_sections=tr
         reinit!(wing.refined_sections[2],
             s2.LE_point, s2.TE_point, s2.aero_model,
             reuse_aero_data ? nothing : s2.aero_data)
-        recompute_mapping && compute_refined_panel_mapping!(wing)
+        if recompute_mapping
+            compute_refined_panel_mapping!(wing)
+            compute_refined_section_interpolation!(wing)
+        end
         update_non_deformed_sections!(wing)
         return nothing
     end
@@ -855,7 +879,10 @@ function refine!(wing::AbstractWing{T}; recompute_mapping=true, sort_sections=tr
     end
 
     # Compute panel mapping by finding closest unrefined section for each refined panel
-    recompute_mapping && compute_refined_panel_mapping!(wing)
+    if recompute_mapping
+        compute_refined_panel_mapping!(wing)
+        compute_refined_section_interpolation!(wing)
+    end
 
     # Update n_unrefined_sections based on actual sections
     wing.n_unrefined_sections = Int16(length(wing.unrefined_sections))
@@ -936,6 +963,90 @@ function compute_refined_panel_mapping!(wing::AbstractWing)
     return nothing
 end
 
+"""
+    compute_refined_section_interpolation!(wing::AbstractWing)
+
+Compute per-refined-section linear-interpolation weights from the unrefined
+sections. For refined section i, the interpolated value is:
+
+    out[i] = weight[i] * unrefined[left_idx[i]] +
+             (1 - weight[i]) * unrefined[left_idx[i] + 1]
+
+Positions are quarter-chord arc-length along the unrefined and refined sections.
+Endpoint refined sections always get `weight == 1` so they take the unrefined
+endpoint value exactly.
+"""
+function compute_refined_section_interpolation!(wing::AbstractWing{T}) where {T}
+    n_unref = length(wing.unrefined_sections)
+    n_sections = wing.n_panels + 1
+
+    if length(wing.refined_section_left_idx) != n_sections
+        resize!(wing.refined_section_left_idx, n_sections)
+    end
+    if length(wing.refined_section_weight) != n_sections
+        resize!(wing.refined_section_weight, n_sections)
+    end
+
+    if n_unref < 2
+        for i in 1:n_sections
+            wing.refined_section_left_idx[i] = Int16(1)
+            wing.refined_section_weight[i] = one(T)
+        end
+        return nothing
+    end
+
+    @inline qc(s, j) = s.LE_point[j] + 0.25 * (s.TE_point[j] - s.LE_point[j])
+
+    s_unref = Vector{Float64}(undef, n_unref)
+    s_unref[1] = 0.0
+    for k in 2:n_unref
+        u_prev = wing.unrefined_sections[k - 1]
+        u_cur = wing.unrefined_sections[k]
+        d = 0.0
+        for j in 1:3
+            dq = qc(u_cur, j) - qc(u_prev, j)
+            d += dq * dq
+        end
+        s_unref[k] = s_unref[k - 1] + sqrt(d)
+    end
+
+    s_ref = Vector{Float64}(undef, n_sections)
+    s_ref[1] = 0.0
+    for i in 2:n_sections
+        r_prev = wing.refined_sections[i - 1]
+        r_cur = wing.refined_sections[i]
+        d = 0.0
+        for j in 1:3
+            dq = qc(r_cur, j) - qc(r_prev, j)
+            d += dq * dq
+        end
+        s_ref[i] = s_ref[i - 1] + sqrt(d)
+    end
+
+    for i in 1:n_sections
+        target = s_ref[i]
+        left = 1
+        for k in 1:(n_unref - 1)
+            if s_unref[k + 1] >= target || k == n_unref - 1
+                left = k
+                break
+            end
+        end
+        seg = s_unref[left + 1] - s_unref[left]
+        t = seg > 1e-30 ? (target - s_unref[left]) / seg : 0.0
+        t = clamp(t, 0.0, 1.0)
+        wing.refined_section_left_idx[i] = Int16(left)
+        wing.refined_section_weight[i] = T(1.0 - t)
+    end
+
+    wing.refined_section_left_idx[1] = Int16(1)
+    wing.refined_section_weight[1] = one(T)
+    wing.refined_section_left_idx[n_sections] = Int16(n_unref - 1)
+    wing.refined_section_weight[n_sections] = zero(T)
+
+    return nothing
+end
+
 
 """
     calculate_new_aero_data(sections, section_index, left_weight, right_weight)
diff --git a/test/yaml_geometry/test_yaml_wing_deformation.jl b/test/yaml_geometry/test_yaml_wing_deformation.jl
index 6194e676..c0d71ff7 100644
--- a/test/yaml_geometry/test_yaml_wing_deformation.jl
+++ b/test/yaml_geometry/test_yaml_wing_deformation.jl
@@ -215,98 +215,60 @@ using Test
         @test body_aero.panels[4].delta ≈ deg2rad(20.0) atol=1e-6
     end
 
-    @testset "unrefined_deform! Maps to Panels" begin
-        # Test that unrefined_deform! correctly maps unrefined sections to panels
-        # Use complex_wing which has 7 unrefined sections
+    @testset "unrefined_deform! Linear Interpolation" begin
+        # unrefined_deform! linearly interpolates angles from unrefined sections to
+        # refined sections; panel-level dist values are the average of adjacent
+        # refined-section values.
         complex_wing_file = test_data_path("yaml_geometry", "complex_wing.yaml")
         wing = Wing(complex_wing_file; n_panels=12)
         refine!(wing)
         body_aero = BodyAerodynamics([wing])
 
-        # Verify we have 7 unrefined sections
         @test wing.n_unrefined_sections == 7
 
-        # Create unrefined section angles (7 sections)
-        # These will be mapped to panels via refined_panel_mapping
         theta_unrefined = deg2rad.([10.0, 15.0, 20.0, 25.0, 20.0, 15.0, 10.0])
         delta_unrefined = deg2rad.([5.0, 7.5, 10.0, 12.5, 10.0, 7.5, 5.0])
 
-        # Apply using unrefined_deform!
         VortexStepMethod.unrefined_deform!(wing, theta_unrefined, delta_unrefined)
         VortexStepMethod.reinit!(body_aero)
 
-        # Each panel should have the delta from its mapped unrefined section
-        for i in 1:wing.n_panels
-            unrefined_idx = wing.refined_panel_mapping[i]
-            expected_delta = delta_unrefined[unrefined_idx]
-            @test body_aero.panels[i].delta ≈ expected_delta atol=1e-6
+        # Endpoint refined sections take exact unrefined endpoint values, so the
+        # first/last panel value is the average of an endpoint and its neighbor.
+        n_sec = wing.n_panels + 1
+        @test wing.refined_section_weight[1] ≈ 1.0
+        @test wing.refined_section_weight[n_sec] ≈ 0.0
+        @test wing.refined_section_left_idx[1] == 1
+        @test wing.refined_section_left_idx[n_sec] == wing.n_unrefined_sections - 1
+
+        # All panel deltas lie within the input range (no overshoot from interp).
+        delta_min = minimum(delta_unrefined)
+        delta_max = maximum(delta_unrefined)
+        for p in body_aero.panels
+            @test delta_min - 1e-9 ≤ p.delta ≤ delta_max + 1e-9
         end
     end
 
-    @testset "Smooth vs Non-Smooth Deformation" begin
-        # Create test wing with 2 unrefined sections, refined to 40 panels
+    @testset "unrefined_deform! Endpoint and Linearity" begin
         simple_wing_file = test_data_path("yaml_geometry", "simple_wing.yaml")
         wing = Wing(simple_wing_file; n_panels=40)
         refine!(wing)
         @test wing.n_unrefined_sections == 2
 
-        # Define varying input angles at unrefined section level
         delta_input = deg2rad.([0.0, 10.0])
-
-        # Test 1: Non-smooth deformation has step-wise discontinuities
-        VortexStepMethod.unrefined_deform!(wing, nothing, delta_input; smooth=false)
-        delta_nonsmooth = copy(wing.delta_dist)
-
-        # Verify step-wise pattern: panels in same unrefined section have identical angles
-        for i in 1:wing.n_panels
-            unrefined_idx = wing.refined_panel_mapping[i]
-            @test delta_nonsmooth[i] ≈ delta_input[unrefined_idx] atol=1e-10
-        end
-
-        # Verify discontinuities exist at unrefined section boundaries
-        # Find boundary indices (where panel mapping changes)
-        max_gradient_nonsmooth = 0.0
-        for i in 1:(wing.n_panels-1)
-            if wing.refined_panel_mapping[i] != wing.refined_panel_mapping[i+1]
-                gradient = abs(delta_nonsmooth[i+1] - delta_nonsmooth[i])
-                max_gradient_nonsmooth = max(max_gradient_nonsmooth, gradient)
-            end
+        VortexStepMethod.unrefined_deform!(wing, nothing, delta_input)
+        delta_panels = copy(wing.delta_dist)
+
+        # With two unrefined sections, the refined-section interpolation is exactly
+        # linear in arc-length, so panel-level values (averages of adjacent refined
+        # sections) are linear and monotonic — no discontinuities.
+        for i in 1:(wing.n_panels - 1)
+            @test delta_panels[i + 1] ≥ delta_panels[i] - 1e-12
         end
-        @test max_gradient_nonsmooth > deg2rad(5.0)  # Should have large jumps
-
-        # Test 2: Smooth deformation is continuous
-        VortexStepMethod.unrefined_deform!(wing, nothing, delta_input; smooth=true)
-        delta_smooth = copy(wing.delta_dist)
-
-        # Verify gradients between adjacent panels are small
-        max_gradient_smooth = maximum(abs.(diff(delta_smooth)))
-        @test max_gradient_smooth < deg2rad(3.0)  # Should be smooth
+        max_gradient = maximum(abs.(diff(delta_panels)))
+        @test max_gradient < (delta_input[2] - delta_input[1]) / wing.n_panels * 1.5
 
-        # Verify no sharp discontinuities
-        for i in 1:(wing.n_panels-1)
-            @test abs(delta_smooth[i+1] - delta_smooth[i]) < deg2rad(3.0)
-        end
-
-        # Test 3: Smoothing reduces maximum gradient
-        @test max_gradient_smooth < max_gradient_nonsmooth
-
-        # Test 4: Angles match input at unrefined section centers (both modes)
-        # For non-smooth: extract angle at center panel of each unrefined section
-        VortexStepMethod.unrefined_deform!(wing, nothing, delta_input; smooth=false)
-        for i in 1:wing.n_unrefined_sections
-            # Find panels belonging to this unrefined section
-            panel_indices = findall(==(i), wing.refined_panel_mapping)
-            center_panel_idx = panel_indices[div(length(panel_indices), 2) + 1]
-            @test wing.delta_dist[center_panel_idx] ≈ delta_input[i] atol=1e-10
-        end
-
-        # For smooth: angles at center should be close to input (tolerance larger due to smoothing)
-        VortexStepMethod.unrefined_deform!(wing, nothing, delta_input; smooth=true)
-        for i in 1:wing.n_unrefined_sections
-            panel_indices = findall(==(i), wing.refined_panel_mapping)
-            center_panel_idx = panel_indices[div(length(panel_indices), 2) + 1]
-            # Smoothing may shift values slightly, use absolute tolerance for small angles
-            @test wing.delta_dist[center_panel_idx] ≈ delta_input[i] atol=deg2rad(2.0)
-        end
+        # Refined-section endpoints match unrefined endpoints exactly.
+        @test wing.refined_section_weight[1] ≈ 1.0
+        @test wing.refined_section_weight[end] ≈ 0.0
     end
 end

From 48d83c868d7006ef49dd28510cd92041566b72b3 Mon Sep 17 00:00:00 2001
From: Bart <bart@vandelint.net>
Date: Fri, 8 May 2026 10:21:31 +0200
Subject: [PATCH 2/8] Exact middle twist

---
 src/wing_geometry.jl                    | 43 +++++++++++++++++++------
 test/ram_geometry/test_kite_geometry.jl | 15 +++++----
 2 files changed, 41 insertions(+), 17 deletions(-)

diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl
index ce4f38db..5c36851d 100644
--- a/src/wing_geometry.jl
+++ b/src/wing_geometry.jl
@@ -552,22 +552,45 @@ function _apply_refined_section_thetas!(wing::Wing{P, T}, section_thetas) where
     local_y = zeros(MVector{3, T})
     chord = zeros(MVector{3, T})
     normal = zeros(MVector{3, T})
+    n_sec = wing.n_panels + 1
 
-    for i in 1:(wing.n_panels + 1)
+    for i in 1:n_sec
         theta = section_thetas[i]
         section = wing.non_deformed_sections[i]
-
-        if i < wing.n_panels + 1
-            section2 = wing.non_deformed_sections[i + 1]
-            local_y .= normalize(section.LE_point - section2.LE_point)
-        else
-            section_prev = wing.non_deformed_sections[i - 1]
-            local_y .= normalize(section_prev.LE_point - section.LE_point)
+        le = section.LE_point
+
+        # Spanwise axis: average of unit vectors to the left and right neighbour LE
+        # points so curvature does not bias the twist axis. Boundaries fall back to
+        # whichever side exists.
+        local_y .= zero(T)
+        if i > 1
+            le_prev = wing.non_deformed_sections[i - 1].LE_point
+            dx = le_prev[1] - le[1]
+            dy = le_prev[2] - le[2]
+            dz = le_prev[3] - le[3]
+            inv_n = 1 / sqrt(dx * dx + dy * dy + dz * dz)
+            local_y[1] += dx * inv_n
+            local_y[2] += dy * inv_n
+            local_y[3] += dz * inv_n
+        end
+        if i < n_sec
+            le_next = wing.non_deformed_sections[i + 1].LE_point
+            dx = le[1] - le_next[1]
+            dy = le[2] - le_next[2]
+            dz = le[3] - le_next[3]
+            inv_n = 1 / sqrt(dx * dx + dy * dy + dz * dz)
+            local_y[1] += dx * inv_n
+            local_y[2] += dy * inv_n
+            local_y[3] += dz * inv_n
         end
+        inv_n = 1 / sqrt(local_y[1]^2 + local_y[2]^2 + local_y[3]^2)
+        local_y[1] *= inv_n
+        local_y[2] *= inv_n
+        local_y[3] *= inv_n
 
-        chord .= section.TE_point .- section.LE_point
+        chord .= section.TE_point .- le
         normal .= chord × local_y
-        @. wing.refined_sections[i].TE_point = section.LE_point +
+        @. wing.refined_sections[i].TE_point = le +
             cos(theta) * chord - sin(theta) * normal
     end
     return nothing
diff --git a/test/ram_geometry/test_kite_geometry.jl b/test/ram_geometry/test_kite_geometry.jl
index 7df146dd..31e921e1 100644
--- a/test/ram_geometry/test_kite_geometry.jl
+++ b/test/ram_geometry/test_kite_geometry.jl
@@ -190,22 +190,23 @@ using Serialization
         wing = ObjWing(test_obj_path, test_dat_path; remove_nan=true)
         body_aero = BodyAerodynamics([wing])
 
-        # Store original TE point for comparison
-        i = length(body_aero.panels) ÷ 2
+        # Sample the apex panel: with even n_panels, refined section (n_panels/2 + 1)
+        # sits exactly at γ=0, so panels[n_panels/2 + 1].TE_point_1 is the apex TE
+        # where the LE tangent is purely along y and twist preserves y exactly.
+        @test iseven(wing.n_panels)
+        i = wing.n_panels ÷ 2 + 1
         original_te_point = copy(body_aero.panels[i].TE_point_1)
 
-        # Apply deformation with non-zero angles
-        theta_dist = fill(deg2rad(30.0), wing.n_panels)  # 30 degrees twist for all panels
-        delta_dist = fill(deg2rad(5.0), wing.n_panels)   # 5 degrees TE deflection for all panels
+        theta_dist = fill(deg2rad(30.0), wing.n_panels)
+        delta_dist = fill(deg2rad(5.0), wing.n_panels)
 
         VortexStepMethod.deform!(wing, theta_dist, delta_dist)
         VortexStepMethod.reinit!(body_aero)
 
-        # Check if TE point changed after deformation
         deformed_te_point = copy(body_aero.panels[i].TE_point_1)
         @test !isapprox(original_te_point, deformed_te_point, atol=1e-2)
         @test deformed_te_point[3] < original_te_point[3] # right hand rule
-        @test deformed_te_point[2] ≈ original_te_point[2] atol=1e-2 # right hand rule
+        @test deformed_te_point[2] ≈ original_te_point[2] atol=1e-5 # right hand rule
         @test deformed_te_point[1] < original_te_point[1] # right hand rule
         @test body_aero.panels[i].delta ≈ deg2rad(5.0)
 

From b744bc877d58a900239d7eeb6b62d721d8f0d7e1 Mon Sep 17 00:00:00 2001
From: Bart <bart@vandelint.net>
Date: Tue, 12 May 2026 10:31:54 +0200
Subject: [PATCH 3/8] Update default manifests

---
 Manifest-v1.11.toml.default |  86 +++++++++++---------
 Manifest-v1.12.toml.default | 156 ++++++++++++++++++++----------------
 2 files changed, 135 insertions(+), 107 deletions(-)

diff --git a/Manifest-v1.11.toml.default b/Manifest-v1.11.toml.default
index 5deb9545..1fdccfe7 100644
--- a/Manifest-v1.11.toml.default
+++ b/Manifest-v1.11.toml.default
@@ -2,12 +2,12 @@
 
 julia_version = "1.11.9"
 manifest_format = "2.0"
-project_hash = "2208ea1100be74002e5a4dc77070f9be50832d18"
+project_hash = "1ad349ec1b7d8858e1bce3c1d99302d12b82e8e6"
 
 [[deps.ADTypes]]
-git-tree-sha1 = "f7304359109c768cf32dc5fa2d371565bb63b68a"
+git-tree-sha1 = "bbc22a9a08a0ef6460041086d8a7b27940ed4ffd"
 uuid = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
-version = "1.21.0"
+version = "1.22.0"
 weakdeps = ["ChainRulesCore", "ConstructionBase", "EnzymeCore"]
 
     [deps.ADTypes.extensions]
@@ -229,9 +229,9 @@ version = "1.15.1"
 
 [[deps.DifferentiationInterface]]
 deps = ["ADTypes", "LinearAlgebra"]
-git-tree-sha1 = "7ae99144ea44715402c6c882bfef2adbeadbc4ce"
+git-tree-sha1 = "2147a95a217cc8a78ec96ee03581adf129468e49"
 uuid = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
-version = "0.7.16"
+version = "0.7.18"
 
     [deps.DifferentiationInterface.extensions]
     DifferentiationInterfaceChainRulesCoreExt = "ChainRulesCore"
@@ -241,8 +241,9 @@ version = "0.7.16"
     DifferentiationInterfaceFiniteDiffExt = "FiniteDiff"
     DifferentiationInterfaceFiniteDifferencesExt = "FiniteDifferences"
     DifferentiationInterfaceForwardDiffExt = ["ForwardDiff", "DiffResults"]
-    DifferentiationInterfaceGPUArraysCoreExt = "GPUArraysCore"
+    DifferentiationInterfaceGPUArraysCoreExt = ["GPUArraysCore", "Adapt"]
     DifferentiationInterfaceGTPSAExt = "GTPSA"
+    DifferentiationInterfaceHyperHessiansExt = "HyperHessians"
     DifferentiationInterfaceMooncakeExt = "Mooncake"
     DifferentiationInterfacePolyesterForwardDiffExt = ["PolyesterForwardDiff", "ForwardDiff", "DiffResults"]
     DifferentiationInterfaceReverseDiffExt = ["ReverseDiff", "DiffResults"]
@@ -255,6 +256,7 @@ version = "0.7.16"
     DifferentiationInterfaceZygoteExt = ["Zygote", "ForwardDiff"]
 
     [deps.DifferentiationInterface.weakdeps]
+    Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
     ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
     DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
     Diffractor = "9f5e2b26-1114-432f-b630-d3fe2085c51c"
@@ -266,6 +268,7 @@ version = "0.7.16"
     ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
     GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
     GTPSA = "b27dd330-f138-47c5-815b-40db9dd9b6e8"
+    HyperHessians = "06b494a0-c8e0-40cc-ad32-d99506a00a6c"
     Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6"
     PolyesterForwardDiff = "98d1487c-24ca-40b6-b7ab-df2af84e126b"
     ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
@@ -298,9 +301,9 @@ uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 version = "1.0.7"
 
 [[deps.EnzymeCore]]
-git-tree-sha1 = "24bbb6fc8fb87eb71c1f8d00184a60fc22c63903"
+git-tree-sha1 = "c6ee69ee502060982d12dbaaf3d8fcb4e835a0d1"
 uuid = "f151be2c-9106-41f4-ab19-57ee4f262869"
-version = "0.8.19"
+version = "0.8.20"
 weakdeps = ["Adapt", "ChainRulesCore"]
 
     [deps.EnzymeCore.extensions]
@@ -323,9 +326,9 @@ version = "1.11.0"
 
 [[deps.FiniteDiff]]
 deps = ["ArrayInterface", "LinearAlgebra", "Setfield"]
-git-tree-sha1 = "73e879af0e767bd6dfade7c5b09d7b05657a8284"
+git-tree-sha1 = "f7017a4f337f8df189fcce98e32b67a1298a2115"
 uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
-version = "2.30.0"
+version = "2.31.0"
 
     [deps.FiniteDiff.extensions]
     FiniteDiffBandedMatricesExt = "BandedMatrices"
@@ -362,9 +365,9 @@ version = "1.1.3"
 
 [[deps.FunctionWrappersWrappers]]
 deps = ["FunctionWrappers", "PrecompileTools", "TruncatedStacktraces"]
-git-tree-sha1 = "ce6762f8f0e7542534f01523ae051e625cbf0468"
+git-tree-sha1 = "c1b0c3a166a2a393257aa888787ca817532e14ce"
 uuid = "77dc65aa-8811-40c2-897b-53d922fa7daf"
-version = "1.5.0"
+version = "1.8.0"
 
     [deps.FunctionWrappersWrappers.extensions]
     FunctionWrappersWrappersEnzymeExt = ["Enzyme", "EnzymeCore"]
@@ -436,9 +439,9 @@ version = "1.0.0"
 
 [[deps.JLLWrappers]]
 deps = ["Artifacts", "Preferences"]
-git-tree-sha1 = "0533e564aae234aff59ab625543145446d8b6ec2"
+git-tree-sha1 = "7204148362dafe5fe6a273f855b8ccbe4df8173e"
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.7.1"
+version = "1.8.0"
 
 [[deps.Krylov]]
 deps = ["LinearAlgebra", "Printf", "SparseArrays"]
@@ -493,9 +496,9 @@ version = "1.18.0+0"
 
 [[deps.LineSearch]]
 deps = ["ADTypes", "CommonSolve", "ConcreteStructs", "FastClosures", "LinearAlgebra", "MaybeInplace", "PrecompileTools", "SciMLBase", "SciMLJacobianOperators", "StaticArraysCore"]
-git-tree-sha1 = "69da095e4c24ed3c4a168bb76dc9c620a6d7239c"
+git-tree-sha1 = "fd58a77c92e7c8f1db25c9839127d52943a49349"
 uuid = "87fe0de2-c867-4266-b59a-2f0a94fc965b"
-version = "0.1.7"
+version = "0.1.9"
 
     [deps.LineSearch.extensions]
     LineSearchLineSearchesExt = "LineSearches"
@@ -510,9 +513,9 @@ version = "1.11.0"
 
 [[deps.LinearSolve]]
 deps = ["ArrayInterface", "ConcreteStructs", "DocStringExtensions", "EnumX", "GPUArraysCore", "InteractiveUtils", "Krylov", "Libdl", "LinearAlgebra", "MKL_jll", "Markdown", "OpenBLAS_jll", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLLogging", "SciMLOperators", "Setfield", "StaticArraysCore"]
-git-tree-sha1 = "42b5cb44317e89ef75dd841c9c8eba9045bf9ff0"
+git-tree-sha1 = "04940e970234c6ac650eec9adfc62677cca2328d"
 uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-version = "3.75.0"
+version = "3.76.0"
 
     [deps.LinearSolve.extensions]
     LinearSolveAMDGPUExt = "AMDGPU"
@@ -537,12 +540,12 @@ version = "3.75.0"
     LinearSolveKrylovKitExt = "KrylovKit"
     LinearSolveMetalExt = "Metal"
     LinearSolveMooncakeExt = "Mooncake"
-    LinearSolvePETScCSRExt = ["PETSc", "SparseArrays", "SparseMatricesCSR"]
-    LinearSolvePETScExt = ["PETSc", "SparseArrays"]
+    LinearSolvePETScExt = ["PETSc", "SparseArrays", "SparseMatricesCSR"]
     LinearSolvePETScMPIExt = ["PETSc", "PartitionedArrays", "SparseArrays", "SparseMatricesCSR"]
     LinearSolveParUExt = ["ParU_jll", "SparseArrays"]
     LinearSolvePardisoExt = ["Pardiso", "SparseArrays"]
     LinearSolveRecursiveFactorizationExt = "RecursiveFactorization"
+    LinearSolveSTRUMPACKExt = ["SparseArrays", "STRUMPACK_jll"]
     LinearSolveSparseArraysExt = "SparseArrays"
     LinearSolveSparspakExt = ["SparseArrays", "Sparspak"]
 
@@ -574,6 +577,7 @@ version = "3.75.0"
     Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
     PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
     RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4"
+    STRUMPACK_jll = "86fbd0b9-476f-557c-b766-62c724b42d8c"
     SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
     SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
     Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac"
@@ -661,9 +665,9 @@ version = "1.2.0"
 
 [[deps.NonlinearSolve]]
 deps = ["ADTypes", "ArrayInterface", "BracketingNonlinearSolve", "CommonSolve", "ConcreteStructs", "DifferentiationInterface", "FastClosures", "FiniteDiff", "ForwardDiff", "LineSearch", "LinearAlgebra", "LinearSolve", "NonlinearSolveBase", "NonlinearSolveFirstOrder", "NonlinearSolveQuasiNewton", "NonlinearSolveSpectralMethods", "PrecompileTools", "Preferences", "Reexport", "SciMLBase", "SciMLLogging", "Setfield", "SimpleNonlinearSolve", "StaticArraysCore", "SymbolicIndexingInterface"]
-git-tree-sha1 = "e88921859836899abe94d08ea0fd42137067280e"
+git-tree-sha1 = "ae0b0f875e4dfda538ab6977ec7939126ee9a1fe"
 uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-version = "4.17.1"
+version = "4.19.0"
 
     [deps.NonlinearSolve.extensions]
     NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt"
@@ -853,9 +857,9 @@ version = "1.3.4"
 
 [[deps.RecursiveArrayTools]]
 deps = ["Adapt", "ArrayInterface", "DocStringExtensions", "GPUArraysCore", "LinearAlgebra", "PrecompileTools", "RecipesBase", "StaticArraysCore", "SymbolicIndexingInterface"]
-git-tree-sha1 = "79a5a1a5b6294c54602d18a59026ec0b07471054"
+git-tree-sha1 = "57b6fb3932fc8d1fc911f840d2c9de5fe3ba5008"
 uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
-version = "4.2.0"
+version = "4.3.0"
 
     [deps.RecursiveArrayTools.extensions]
     RecursiveArrayToolsCUDAExt = "CUDA"
@@ -912,11 +916,16 @@ version = "0.5.18"
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 version = "0.7.0"
 
+[[deps.SafeTestsets]]
+git-tree-sha1 = "81ec49d645af090901120a1542e67ecbbe044db3"
+uuid = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+version = "0.1.0"
+
 [[deps.SciMLBase]]
 deps = ["ADTypes", "Accessors", "Adapt", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "PreallocationTools", "PrecompileTools", "Preferences", "Printf", "Random", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLLogging", "SciMLOperators", "SciMLPublic", "SciMLStructures", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface"]
-git-tree-sha1 = "4fdad3606c60fbbd52424737c31ec4141672c809"
+git-tree-sha1 = "fed3c3d9781685d982256f94ac85ad4ddf8ee78f"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-version = "3.3.0"
+version = "3.10.0"
 
     [deps.SciMLBase.extensions]
     SciMLBaseChainRulesCoreExt = "ChainRulesCore"
@@ -963,9 +972,9 @@ version = "0.1.13"
 
 [[deps.SciMLLogging]]
 deps = ["Logging", "LoggingExtras", "Preferences"]
-git-tree-sha1 = "0161be062570af4042cf6f69e3d5d0b0555b6927"
+git-tree-sha1 = "35bdd8e578af78030d98b262461e8eedc045d943"
 uuid = "a6db7da4-7206-11f0-1eab-35f2a5dbe1d1"
-version = "1.9.1"
+version = "1.10.0"
 
     [deps.SciMLLogging.extensions]
     SciMLLoggingTracyExt = "Tracy"
@@ -974,16 +983,21 @@ version = "1.9.1"
     Tracy = "e689c965-62c8-4b79-b2c5-8359227902fd"
 
 [[deps.SciMLOperators]]
-deps = ["Accessors", "ArrayInterface", "DocStringExtensions", "LinearAlgebra"]
-git-tree-sha1 = "234869cf9fee9258a95464b7a7065cc7be84db00"
+deps = ["Accessors", "Adapt", "ArrayInterface", "DocStringExtensions", "LinearAlgebra", "SafeTestsets"]
+git-tree-sha1 = "ef54b9bcfc6694b11aea8458619d1416885bec1e"
 uuid = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
-version = "1.16.0"
-weakdeps = ["SparseArrays", "StaticArraysCore"]
+version = "1.20.0"
 
     [deps.SciMLOperators.extensions]
+    SciMLOperatorsLoopVectorizationExt = "LoopVectorization"
     SciMLOperatorsSparseArraysExt = "SparseArrays"
     SciMLOperatorsStaticArraysCoreExt = "StaticArraysCore"
 
+    [deps.SciMLOperators.weakdeps]
+    LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
+    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+    StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+
 [[deps.SciMLPublic]]
 git-tree-sha1 = "0ba076dbdce87ba230fff48ca9bca62e1f345c9b"
 uuid = "431bcebd-1456-4ced-9d72-93c2757fff0b"
@@ -1089,9 +1103,9 @@ version = "7.7.0+0"
 
 [[deps.SymbolicIndexingInterface]]
 deps = ["Accessors", "ArrayInterface", "RuntimeGeneratedFunctions", "StaticArraysCore"]
-git-tree-sha1 = "94c58884e013efff548002e8dc2fdd1cb74dfce5"
+git-tree-sha1 = "173ecfe5f7c5a36043b5f2b8cecaa30a4fc958ef"
 uuid = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
-version = "0.3.46"
+version = "0.3.47"
 
     [deps.SymbolicIndexingInterface.extensions]
     SymbolicIndexingInterfacePrettyTablesExt = "PrettyTables"
@@ -1182,9 +1196,9 @@ version = "1.59.0+0"
 
 [[deps.oneTBB_jll]]
 deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl"]
-git-tree-sha1 = "1350188a69a6e46f799d3945beef36435ed7262f"
+git-tree-sha1 = "da8c1f6eee04831f14edcfa5dae611d309807e57"
 uuid = "1317d2d5-d96f-522e-a858-c73665f53c3e"
-version = "2022.0.0+1"
+version = "2022.3.0+0"
 
 [[deps.p7zip_jll]]
 deps = ["Artifacts", "Libdl"]
diff --git a/Manifest-v1.12.toml.default b/Manifest-v1.12.toml.default
index 8284840d..7f9885b1 100644
--- a/Manifest-v1.12.toml.default
+++ b/Manifest-v1.12.toml.default
@@ -1,13 +1,13 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.12.6"
+julia_version = "1.12.5"
 manifest_format = "2.0"
-project_hash = "ba02310a71881cfe5a42404e47e65d7141af4c4a"
+project_hash = "7523bf605a533ebe3221a67767656d2d45635bb2"
 
 [[deps.ADTypes]]
-git-tree-sha1 = "f7304359109c768cf32dc5fa2d371565bb63b68a"
+git-tree-sha1 = "bbc22a9a08a0ef6460041086d8a7b27940ed4ffd"
 uuid = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
-version = "1.21.0"
+version = "1.22.0"
 weakdeps = ["ChainRulesCore", "ConstructionBase", "EnzymeCore"]
 
     [deps.ADTypes.extensions]
@@ -234,9 +234,9 @@ version = "0.15.9"
 
 [[deps.Cairo_jll]]
 deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "Libdl", "Pixman_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
-git-tree-sha1 = "d0efe2c6fdcdaa1c161d206aa8b933788397ec71"
+git-tree-sha1 = "1fa950ebc3e37eccd51c6a8fe1f92f7d86263522"
 uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
-version = "1.18.6+0"
+version = "1.18.7+0"
 
 [[deps.ChainRulesCore]]
 deps = ["Compat", "LinearAlgebra"]
@@ -466,9 +466,9 @@ version = "1.15.1"
 
 [[deps.DifferentiationInterface]]
 deps = ["ADTypes", "LinearAlgebra"]
-git-tree-sha1 = "7ae99144ea44715402c6c882bfef2adbeadbc4ce"
+git-tree-sha1 = "2147a95a217cc8a78ec96ee03581adf129468e49"
 uuid = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
-version = "0.7.16"
+version = "0.7.18"
 
     [deps.DifferentiationInterface.extensions]
     DifferentiationInterfaceChainRulesCoreExt = "ChainRulesCore"
@@ -478,8 +478,9 @@ version = "0.7.16"
     DifferentiationInterfaceFiniteDiffExt = "FiniteDiff"
     DifferentiationInterfaceFiniteDifferencesExt = "FiniteDifferences"
     DifferentiationInterfaceForwardDiffExt = ["ForwardDiff", "DiffResults"]
-    DifferentiationInterfaceGPUArraysCoreExt = "GPUArraysCore"
+    DifferentiationInterfaceGPUArraysCoreExt = ["GPUArraysCore", "Adapt"]
     DifferentiationInterfaceGTPSAExt = "GTPSA"
+    DifferentiationInterfaceHyperHessiansExt = "HyperHessians"
     DifferentiationInterfaceMooncakeExt = "Mooncake"
     DifferentiationInterfacePolyesterForwardDiffExt = ["PolyesterForwardDiff", "ForwardDiff", "DiffResults"]
     DifferentiationInterfaceReverseDiffExt = ["ReverseDiff", "DiffResults"]
@@ -492,6 +493,7 @@ version = "0.7.16"
     DifferentiationInterfaceZygoteExt = ["Zygote", "ForwardDiff"]
 
     [deps.DifferentiationInterface.weakdeps]
+    Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
     ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
     DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
     Diffractor = "9f5e2b26-1114-432f-b630-d3fe2085c51c"
@@ -503,6 +505,7 @@ version = "0.7.16"
     ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
     GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
     GTPSA = "b27dd330-f138-47c5-815b-40db9dd9b6e8"
+    HyperHessians = "06b494a0-c8e0-40cc-ad32-d99506a00a6c"
     Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6"
     PolyesterForwardDiff = "98d1487c-24ca-40b6-b7ab-df2af84e126b"
     ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
@@ -521,9 +524,9 @@ version = "1.11.0"
 
 [[deps.Distributions]]
 deps = ["AliasTables", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns"]
-git-tree-sha1 = "12184a8cf11c7cbd90a4db8b2cb2f7b6f057cc46"
+git-tree-sha1 = "e421c1938fafab0165b04dc1a9dbe2a26272952c"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
-version = "0.25.124"
+version = "0.25.125"
 
     [deps.Distributions.extensions]
     DistributionsChainRulesCoreExt = "ChainRulesCore"
@@ -563,9 +566,9 @@ uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 version = "1.0.7"
 
 [[deps.EnzymeCore]]
-git-tree-sha1 = "24bbb6fc8fb87eb71c1f8d00184a60fc22c63903"
+git-tree-sha1 = "c6ee69ee502060982d12dbaaf3d8fcb4e835a0d1"
 uuid = "f151be2c-9106-41f4-ab19-57ee4f262869"
-version = "0.8.19"
+version = "0.8.20"
 weakdeps = ["Adapt", "ChainRulesCore"]
 
     [deps.EnzymeCore.extensions]
@@ -592,9 +595,9 @@ version = "0.1.11"
 
 [[deps.Expat_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "27af30de8b5445644e8ffe3bcb0d72049c089cf1"
+git-tree-sha1 = "8f05e9a2e7c2e3eb524102bb2926c5743c07fbe1"
 uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
-version = "2.7.3+0"
+version = "2.8.0+0"
 
 [[deps.ExprTools]]
 git-tree-sha1 = "27415f162e6028e81c72b82ef756bf321213b6ec"
@@ -608,9 +611,9 @@ version = "0.1.6"
 
 [[deps.FFMPEG_jll]]
 deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "PCRE2_jll", "Zlib_jll", "libaom_jll", "libass_jll", "libfdk_aac_jll", "libva_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
-git-tree-sha1 = "66381d7059b5f3f6162f28831854008040a4e905"
+git-tree-sha1 = "cac41ca6b2d399adfc95e51240566f8a60a80806"
 uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
-version = "8.0.1+1"
+version = "8.1.0+0"
 
 [[deps.FFTA]]
 deps = ["AbstractFFTs", "DocStringExtensions", "LinearAlgebra", "MuladdMacro", "Primes", "Random", "Reexport"]
@@ -625,9 +628,9 @@ version = "0.3.2"
 
 [[deps.FileIO]]
 deps = ["Pkg", "Requires", "UUIDs"]
-git-tree-sha1 = "6522cfb3b8fe97bec632252263057996cbd3de20"
+git-tree-sha1 = "8e9c059d6857607253e837730dbf780b6b151acd"
 uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
-version = "1.18.0"
+version = "1.19.0"
 weakdeps = ["HTTP"]
 
     [deps.FileIO.extensions]
@@ -679,9 +682,9 @@ weakdeps = ["PDMats", "SparseArrays", "StaticArrays", "Statistics"]
 
 [[deps.FiniteDiff]]
 deps = ["ArrayInterface", "LinearAlgebra", "Setfield"]
-git-tree-sha1 = "73e879af0e767bd6dfade7c5b09d7b05657a8284"
+git-tree-sha1 = "f7017a4f337f8df189fcce98e32b67a1298a2115"
 uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
-version = "2.30.0"
+version = "2.31.0"
 
     [deps.FiniteDiff.extensions]
     FiniteDiffBandedMatricesExt = "BandedMatrices"
@@ -753,9 +756,9 @@ version = "1.1.3"
 
 [[deps.FunctionWrappersWrappers]]
 deps = ["FunctionWrappers", "PrecompileTools", "TruncatedStacktraces"]
-git-tree-sha1 = "ce6762f8f0e7542534f01523ae051e625cbf0468"
+git-tree-sha1 = "c1b0c3a166a2a393257aa888787ca817532e14ce"
 uuid = "77dc65aa-8811-40c2-897b-53d922fa7daf"
-version = "1.5.0"
+version = "1.8.0"
 
     [deps.FunctionWrappersWrappers.extensions]
     FunctionWrappersWrappersEnzymeExt = ["Enzyme", "EnzymeCore"]
@@ -833,9 +836,9 @@ version = "3.7.0+0"
 
 [[deps.Git_jll]]
 deps = ["Artifacts", "Expat_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Libiconv_jll", "OpenSSL_jll", "PCRE2_jll", "Zlib_jll"]
-git-tree-sha1 = "dc34a3e3d96b4ed305b641e626dc14c12b7824b8"
+git-tree-sha1 = "0dd4cfb426924210c8f42742751cbde74b27bfa3"
 uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb"
-version = "2.53.0+0"
+version = "2.54.0+0"
 
 [[deps.Glib_jll]]
 deps = ["Artifacts", "GettextRuntime_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Zlib_jll"]
@@ -983,9 +986,9 @@ weakdeps = ["ForwardDiff", "Unitful"]
 
 [[deps.IntervalArithmetic]]
 deps = ["CRlibm", "CoreMath", "MacroTools", "OpenBLASConsistentFPCSR_jll", "Printf", "Random", "RoundingEmulator"]
-git-tree-sha1 = "f1c42fcaca2d8034fe392f3e86c2e0809f75b2a1"
+git-tree-sha1 = "3e6273749a2df3a5c9067657510ad01ba5039a92"
 uuid = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
-version = "1.0.6"
+version = "1.0.8"
 
     [deps.IntervalArithmetic.extensions]
     IntervalArithmeticArblibExt = "Arblib"
@@ -1066,15 +1069,15 @@ weakdeps = ["UnPack"]
 
 [[deps.JLLWrappers]]
 deps = ["Artifacts", "Preferences"]
-git-tree-sha1 = "0533e564aae234aff59ab625543145446d8b6ec2"
+git-tree-sha1 = "7204148362dafe5fe6a273f855b8ccbe4df8173e"
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.7.1"
+version = "1.8.0"
 
 [[deps.JSON]]
 deps = ["Dates", "Logging", "Parsers", "PrecompileTools", "StructUtils", "UUIDs", "Unicode"]
-git-tree-sha1 = "67c6f1f085cb2671c93fe34244c9cccde30f7a26"
+git-tree-sha1 = "fe23330af47b8ab4e135b2ff65f7398c3a2bfc65"
 uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
-version = "1.5.0"
+version = "1.5.2"
 
     [deps.JSON.extensions]
     JSONArrowExt = ["ArrowTypes"]
@@ -1216,9 +1219,9 @@ version = "2.42.0+0"
 
 [[deps.LineSearch]]
 deps = ["ADTypes", "CommonSolve", "ConcreteStructs", "FastClosures", "LinearAlgebra", "MaybeInplace", "PrecompileTools", "SciMLBase", "SciMLJacobianOperators", "StaticArraysCore"]
-git-tree-sha1 = "69da095e4c24ed3c4a168bb76dc9c620a6d7239c"
+git-tree-sha1 = "fd58a77c92e7c8f1db25c9839127d52943a49349"
 uuid = "87fe0de2-c867-4266-b59a-2f0a94fc965b"
-version = "0.1.7"
+version = "0.1.9"
 
     [deps.LineSearch.extensions]
     LineSearchLineSearchesExt = "LineSearches"
@@ -1233,9 +1236,9 @@ version = "1.12.0"
 
 [[deps.LinearSolve]]
 deps = ["ArrayInterface", "ConcreteStructs", "DocStringExtensions", "EnumX", "GPUArraysCore", "InteractiveUtils", "Krylov", "Libdl", "LinearAlgebra", "MKL_jll", "Markdown", "OpenBLAS_jll", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLLogging", "SciMLOperators", "Setfield", "StaticArraysCore"]
-git-tree-sha1 = "42b5cb44317e89ef75dd841c9c8eba9045bf9ff0"
+git-tree-sha1 = "04940e970234c6ac650eec9adfc62677cca2328d"
 uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-version = "3.75.0"
+version = "3.76.0"
 
     [deps.LinearSolve.extensions]
     LinearSolveAMDGPUExt = "AMDGPU"
@@ -1260,12 +1263,12 @@ version = "3.75.0"
     LinearSolveKrylovKitExt = "KrylovKit"
     LinearSolveMetalExt = "Metal"
     LinearSolveMooncakeExt = "Mooncake"
-    LinearSolvePETScCSRExt = ["PETSc", "SparseArrays", "SparseMatricesCSR"]
-    LinearSolvePETScExt = ["PETSc", "SparseArrays"]
+    LinearSolvePETScExt = ["PETSc", "SparseArrays", "SparseMatricesCSR"]
     LinearSolvePETScMPIExt = ["PETSc", "PartitionedArrays", "SparseArrays", "SparseMatricesCSR"]
     LinearSolveParUExt = ["ParU_jll", "SparseArrays"]
     LinearSolvePardisoExt = ["Pardiso", "SparseArrays"]
     LinearSolveRecursiveFactorizationExt = "RecursiveFactorization"
+    LinearSolveSTRUMPACKExt = ["SparseArrays", "STRUMPACK_jll"]
     LinearSolveSparseArraysExt = "SparseArrays"
     LinearSolveSparspakExt = ["SparseArrays", "Sparspak"]
 
@@ -1297,6 +1300,7 @@ version = "3.75.0"
     Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
     PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
     RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4"
+    STRUMPACK_jll = "86fbd0b9-476f-557c-b766-62c724b42d8c"
     SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
     SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
     Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac"
@@ -1466,9 +1470,9 @@ version = "1.3.0"
 
 [[deps.NonlinearSolve]]
 deps = ["ADTypes", "ArrayInterface", "BracketingNonlinearSolve", "CommonSolve", "ConcreteStructs", "DifferentiationInterface", "FastClosures", "FiniteDiff", "ForwardDiff", "LineSearch", "LinearAlgebra", "LinearSolve", "NonlinearSolveBase", "NonlinearSolveFirstOrder", "NonlinearSolveQuasiNewton", "NonlinearSolveSpectralMethods", "PrecompileTools", "Preferences", "Reexport", "SciMLBase", "SciMLLogging", "Setfield", "SimpleNonlinearSolve", "StaticArraysCore", "SymbolicIndexingInterface"]
-git-tree-sha1 = "e88921859836899abe94d08ea0fd42137067280e"
+git-tree-sha1 = "ae0b0f875e4dfda538ab6977ec7939126ee9a1fe"
 uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-version = "4.17.1"
+version = "4.19.0"
 
     [deps.NonlinearSolve.extensions]
     NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt"
@@ -1577,9 +1581,9 @@ version = "1.3.6+0"
 
 [[deps.OpenBLASConsistentFPCSR_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "f2b3b9e52a5eb6a3434c8cca67ad2dde011194f4"
+git-tree-sha1 = "3287ec88df50429a934ebc6cf14606215e27b987"
 uuid = "6cdc7f73-28fd-5e50-80fb-958a8875b1af"
-version = "0.3.30+0"
+version = "0.3.33+0"
 
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
@@ -1684,15 +1688,15 @@ version = "0.12.3"
 
 [[deps.Parsers]]
 deps = ["Dates", "PrecompileTools", "UUIDs"]
-git-tree-sha1 = "7d2f8f21da5db6a806faf7b9b292296da42b2810"
+git-tree-sha1 = "5d5e0a78e971354b1c7bff0655d11fdc1b0e12c8"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "2.8.3"
+version = "2.8.4"
 
 [[deps.Pixman_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl"]
-git-tree-sha1 = "db76b1ecd5e9715f3d043cec13b2ec93ce015d53"
+git-tree-sha1 = "e4a6721aa89e62e5d4217c0b21bd714263779dda"
 uuid = "30392449-352a-5448-841d-b1acce4e97dc"
-version = "0.44.2+0"
+version = "0.46.4+0"
 
 [[deps.Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "Random", "SHA", "TOML", "Tar", "UUIDs", "p7zip_jll"]
@@ -1856,9 +1860,9 @@ version = "1.3.4"
 
 [[deps.RecursiveArrayTools]]
 deps = ["Adapt", "ArrayInterface", "DocStringExtensions", "GPUArraysCore", "LinearAlgebra", "PrecompileTools", "RecipesBase", "StaticArraysCore", "SymbolicIndexingInterface"]
-git-tree-sha1 = "79a5a1a5b6294c54602d18a59026ec0b07471054"
+git-tree-sha1 = "57b6fb3932fc8d1fc911f840d2c9de5fe3ba5008"
 uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
-version = "4.2.0"
+version = "4.3.0"
 
     [deps.RecursiveArrayTools.extensions]
     RecursiveArrayToolsCUDAExt = "CUDA"
@@ -1950,11 +1954,16 @@ git-tree-sha1 = "e24dc23107d426a096d3eae6c165b921e74c18e4"
 uuid = "fdea26ae-647d-5447-a871-4b548cad5224"
 version = "3.7.2"
 
+[[deps.SafeTestsets]]
+git-tree-sha1 = "81ec49d645af090901120a1542e67ecbbe044db3"
+uuid = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+version = "0.1.0"
+
 [[deps.SciMLBase]]
 deps = ["ADTypes", "Accessors", "Adapt", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "PreallocationTools", "PrecompileTools", "Preferences", "Printf", "Random", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLLogging", "SciMLOperators", "SciMLPublic", "SciMLStructures", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface"]
-git-tree-sha1 = "4fdad3606c60fbbd52424737c31ec4141672c809"
+git-tree-sha1 = "fed3c3d9781685d982256f94ac85ad4ddf8ee78f"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-version = "3.3.0"
+version = "3.10.0"
 
     [deps.SciMLBase.extensions]
     SciMLBaseChainRulesCoreExt = "ChainRulesCore"
@@ -2001,9 +2010,9 @@ version = "0.1.13"
 
 [[deps.SciMLLogging]]
 deps = ["Logging", "LoggingExtras", "Preferences"]
-git-tree-sha1 = "0161be062570af4042cf6f69e3d5d0b0555b6927"
+git-tree-sha1 = "35bdd8e578af78030d98b262461e8eedc045d943"
 uuid = "a6db7da4-7206-11f0-1eab-35f2a5dbe1d1"
-version = "1.9.1"
+version = "1.10.0"
 
     [deps.SciMLLogging.extensions]
     SciMLLoggingTracyExt = "Tracy"
@@ -2012,16 +2021,21 @@ version = "1.9.1"
     Tracy = "e689c965-62c8-4b79-b2c5-8359227902fd"
 
 [[deps.SciMLOperators]]
-deps = ["Accessors", "ArrayInterface", "DocStringExtensions", "LinearAlgebra"]
-git-tree-sha1 = "234869cf9fee9258a95464b7a7065cc7be84db00"
+deps = ["Accessors", "Adapt", "ArrayInterface", "DocStringExtensions", "LinearAlgebra", "SafeTestsets"]
+git-tree-sha1 = "ef54b9bcfc6694b11aea8458619d1416885bec1e"
 uuid = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
-version = "1.16.0"
-weakdeps = ["SparseArrays", "StaticArraysCore"]
+version = "1.20.0"
 
     [deps.SciMLOperators.extensions]
+    SciMLOperatorsLoopVectorizationExt = "LoopVectorization"
     SciMLOperatorsSparseArraysExt = "SparseArrays"
     SciMLOperatorsStaticArraysCoreExt = "StaticArraysCore"
 
+    [deps.SciMLOperators.weakdeps]
+    LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
+    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+    StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+
 [[deps.SciMLPublic]]
 git-tree-sha1 = "0ba076dbdce87ba230fff48ca9bca62e1f345c9b"
 uuid = "431bcebd-1456-4ced-9d72-93c2757fff0b"
@@ -2035,9 +2049,9 @@ version = "1.10.0"
 
 [[deps.ScopedValues]]
 deps = ["HashArrayMappedTries", "Logging"]
-git-tree-sha1 = "ac4b837d89a58c848e85e698e2a2514e9d59d8f6"
+git-tree-sha1 = "67a144433c4ce877ee6d1ada69a124d6b1ecf7be"
 uuid = "7e506255-f358-4e82-b7e4-beb19740aa63"
-version = "1.6.0"
+version = "1.6.2"
 
 [[deps.Scratch]]
 deps = ["Dates"]
@@ -2243,9 +2257,9 @@ version = "0.2.3"
 
 [[deps.StructUtils]]
 deps = ["Dates", "UUIDs"]
-git-tree-sha1 = "aab80fbf866600f3299dd7f6656d80e7be177cfe"
+git-tree-sha1 = "dd974aefe288ef2898733aecf40858dc86742d74"
 uuid = "ec057cc2-7a8d-4b58-b3b3-92acb9f63b42"
-version = "2.7.2"
+version = "2.8.1"
 
     [deps.StructUtils.extensions]
     StructUtilsMeasurementsExt = ["Measurements"]
@@ -2272,9 +2286,9 @@ version = "7.8.3+2"
 
 [[deps.SymbolicIndexingInterface]]
 deps = ["Accessors", "ArrayInterface", "RuntimeGeneratedFunctions", "StaticArraysCore"]
-git-tree-sha1 = "94c58884e013efff548002e8dc2fdd1cb74dfce5"
+git-tree-sha1 = "173ecfe5f7c5a36043b5f2b8cecaa30a4fc958ef"
 uuid = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
-version = "0.3.46"
+version = "0.3.47"
 weakdeps = ["PrettyTables"]
 
     [deps.SymbolicIndexingInterface.extensions]
@@ -2407,10 +2421,10 @@ uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
 version = "1.3.0"
 
 [[deps.VortexStepMethod]]
-deps = ["Colors", "DefaultApplication", "DelimitedFiles", "DifferentiationInterface", "FiniteDiff", "Interpolations", "LaTeXStrings", "LinearAlgebra", "Logging", "Measures", "NonlinearSolve", "Parameters", "Pkg", "PreallocationTools", "PrecompileTools", "RecursiveArrayTools", "SciMLBase", "Serialization", "StaticArrays", "Statistics", "StructMapping", "Timers", "Xfoil", "YAML"]
+deps = ["Colors", "DefaultApplication", "DelimitedFiles", "DifferentiationInterface", "FiniteDiff", "ForwardDiff", "Interpolations", "LaTeXStrings", "LinearAlgebra", "Logging", "Measures", "NonlinearSolve", "Parameters", "Pkg", "PreallocationTools", "PrecompileTools", "RecursiveArrayTools", "SciMLBase", "Serialization", "StaticArrays", "Statistics", "StructMapping", "Timers", "Xfoil", "YAML"]
 path = "."
 uuid = "ed3cd733-9f0f-46a9-93e0-89b8d4998dd9"
-version = "3.0.1"
+version = "3.3.0"
 weakdeps = ["ControlPlots", "Makie"]
 
     [deps.VortexStepMethod.extensions]
@@ -2425,9 +2439,9 @@ version = "1.24.0+0"
 
 [[deps.WeakRefStrings]]
 deps = ["DataAPI", "InlineStrings", "Parsers"]
-git-tree-sha1 = "b1be2855ed9ed8eac54e5caff2afcdb442d52c23"
+git-tree-sha1 = "0716e01c3b40413de5dedbc9c5c69f27cddfddfc"
 uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
-version = "1.4.2"
+version = "1.4.3"
 
 [[deps.WebP]]
 deps = ["CEnum", "ColorTypes", "FileIO", "FixedPointNumbers", "ImageCore", "libwebp_jll"]
@@ -2520,9 +2534,9 @@ version = "0.9.12+0"
 
 [[deps.Xorg_libpciaccess_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Zlib_jll"]
-git-tree-sha1 = "4909eb8f1cbf6bd4b1c30dd18b2ead9019ef2fad"
+git-tree-sha1 = "58972370b81423fc546c56a60ed1a009450177c3"
 uuid = "a65dc6b1-eb27-53a1-bb3e-dea574b5389e"
-version = "0.18.1+0"
+version = "0.19.0+0"
 
 [[deps.Xorg_libxcb_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libXau_jll", "Xorg_libXdmcp_jll"]
@@ -2544,9 +2558,9 @@ version = "1.4.7+0"
 
 [[deps.Xorg_xkeyboard_config_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_xkbcomp_jll"]
-git-tree-sha1 = "00af7ebdc563c9217ecc67776d1bbf037dbcebf4"
+git-tree-sha1 = "429722587208f02b1cecbddcd20133df2f1ed796"
 uuid = "33bec58e-1273-512f-9401-5d533626f822"
-version = "2.44.0+0"
+version = "2.47.0+0"
 
 [[deps.Xorg_xtrans_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -2649,9 +2663,9 @@ version = "1.64.0+1"
 
 [[deps.oneTBB_jll]]
 deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl"]
-git-tree-sha1 = "1350188a69a6e46f799d3945beef36435ed7262f"
+git-tree-sha1 = "da8c1f6eee04831f14edcfa5dae611d309807e57"
 uuid = "1317d2d5-d96f-522e-a858-c73665f53c3e"
-version = "2022.0.0+1"
+version = "2022.3.0+0"
 
 [[deps.p7zip_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]

From 1019d826fdd52a0bf03a65679b4d75370c720e68 Mon Sep 17 00:00:00 2001
From: Bart <bart@vandelint.net>
Date: Tue, 12 May 2026 17:45:40 +0200
Subject: [PATCH 4/8] Reduce relaxation factor

---
 test/solver/test_forwarddiff.jl | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/test/solver/test_forwarddiff.jl b/test/solver/test_forwarddiff.jl
index b3a3a9b5..600513d0 100644
--- a/test/solver/test_forwarddiff.jl
+++ b/test/solver/test_forwarddiff.jl
@@ -21,7 +21,7 @@ using Test
 
     @testset "AutoForwardDiff matches AutoFiniteDiff (LOOP, INVISCID)" begin
         # FD warm-starts otherwise bias the reference Jacobian.
-        solver = Solver(body_aero; use_gamma_prev=false)
+        solver = Solver(body_aero; use_gamma_prev=false, relaxation_factor=0.01)
 
         jac_fwd, _, fwd_converged = VortexStepMethod.linearize(
             solver, body_aero, y0;
@@ -36,6 +36,7 @@ using Test
             backend=AutoFiniteDiff(absstep=1e-5, relstep=1e-5))
         @test fd_converged
 
+        @info "INVISCID linearize jacobian norms" norm_fwd=norm(jac_fwd) norm_fd=norm(jac_fd)
         rel_err = maximum(abs.(jac_fwd .- jac_fd)) / maximum(abs, jac_fwd)
         @test rel_err < 1e-4
     end
@@ -73,6 +74,7 @@ using Test
             rtol=1e-7,
             solver_type=LOOP,
             use_gamma_prev=false,
+            relaxation_factor=0.01,
         )
 
         v_a = 15.0
@@ -94,6 +96,7 @@ using Test
             backend=AutoFiniteDiff(absstep=1e-5, relstep=1e-5))
         @test conv_fd
 
+        @info "POLAR_MATRICES linearize jacobian norms" norm_fwd=norm(jac_fwd) norm_fd=norm(jac_fd)
         rel_err = maximum(abs.(jac_fwd .- jac_fd)) / maximum(abs, jac_fwd)
         @test rel_err < 1e-3
     end

From 43a663786d2d8a0781e2535637111b3c3c384564 Mon Sep 17 00:00:00 2001
From: Bart <bart@vandelint.net>
Date: Tue, 12 May 2026 17:47:08 +0200
Subject: [PATCH 5/8] Update default manifest

---
 Manifest-v1.12.toml.default | 34 +++++++++++++++++-----------------
 1 file changed, 17 insertions(+), 17 deletions(-)

diff --git a/Manifest-v1.12.toml.default b/Manifest-v1.12.toml.default
index 7f9885b1..91468f4c 100644
--- a/Manifest-v1.12.toml.default
+++ b/Manifest-v1.12.toml.default
@@ -166,10 +166,10 @@ uuid = "18cc8868-cbac-4acf-b575-c8ff214dc66f"
 version = "1.3.2"
 
 [[deps.BenchmarkTools]]
-deps = ["Compat", "JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
-git-tree-sha1 = "6876e30dc02dc69f0613cb6ece242144f2ca9e56"
+deps = ["Compat", "JSON", "Logging", "PrecompileTools", "Printf", "Profile", "Statistics", "UUIDs"]
+git-tree-sha1 = "9670d3febc2b6da60a0ae57846ba74670290653f"
 uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
-version = "1.7.0"
+version = "1.8.0"
 
 [[deps.BitFlags]]
 git-tree-sha1 = "0691e34b3bb8be9307330f88d1a3c3f25466c24d"
@@ -228,9 +228,9 @@ version = "1.1.1"
 
 [[deps.CairoMakie]]
 deps = ["CRC32c", "Cairo", "Cairo_jll", "Colors", "FileIO", "FreeType", "GeometryBasics", "LinearAlgebra", "Makie", "PrecompileTools"]
-git-tree-sha1 = "fa072933899aae6dc61dde934febed8254e66c6a"
+git-tree-sha1 = "bf2d9cd1ec0c4ce3e0b5aaad192074969413f626"
 uuid = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
-version = "0.15.9"
+version = "0.15.10"
 
 [[deps.Cairo_jll]]
 deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "Libdl", "Pixman_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
@@ -374,9 +374,9 @@ version = "0.6.3"
 
 [[deps.ControlPlots]]
 deps = ["JLD2", "LaTeXStrings", "Pkg", "Printf", "PyPlot", "StaticArraysCore"]
-git-tree-sha1 = "e2df85353a8457a0a18d9248a50a13bef3287389"
+git-tree-sha1 = "fb1cf8d1213c9f2908f8682f58f7ec252f6cc171"
 uuid = "23c2ee80-7a9e-4350-b264-8e670f12517c"
-version = "0.2.14"
+version = "0.2.15"
 
     [deps.ControlPlots.extensions]
     ControlPlotsExt = "ControlSystemsBase"
@@ -408,9 +408,9 @@ version = "1.16.0"
 
 [[deps.DataFrames]]
 deps = ["Compat", "DataAPI", "DataStructures", "Future", "InlineStrings", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrecompileTools", "PrettyTables", "Printf", "Random", "Reexport", "SentinelArrays", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
-git-tree-sha1 = "d8928e9169ff76c6281f39a659f9bca3a573f24c"
+git-tree-sha1 = "5fab31e2e01e70ad66e3e24c968c264d1cf166d6"
 uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-version = "1.8.1"
+version = "1.8.2"
 
 [[deps.DataStructures]]
 deps = ["OrderedCollections"]
@@ -788,9 +788,9 @@ version = "3.4.1+1"
 
 [[deps.GLMakie]]
 deps = ["ColorTypes", "Colors", "FileIO", "FixedPointNumbers", "FreeTypeAbstraction", "GLFW", "GeometryBasics", "LinearAlgebra", "Makie", "Markdown", "MeshIO", "ModernGL", "Observables", "PrecompileTools", "Printf", "ShaderAbstractions", "StaticArrays"]
-git-tree-sha1 = "1e0d427d2c73eb5a7564394df2c9fec8b85e7805"
+git-tree-sha1 = "da0780bbf5f0faa1cdd1567d2dbee0cf841557a7"
 uuid = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
-version = "0.13.9"
+version = "0.13.10"
 
 [[deps.GPUArraysCore]]
 deps = ["Adapt"]
@@ -1356,9 +1356,9 @@ version = "0.5.16"
 
 [[deps.Makie]]
 deps = ["Animations", "Base64", "CRC32c", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "ComputePipeline", "Contour", "Dates", "DelaunayTriangulation", "Distributions", "DocStringExtensions", "Downloads", "FFMPEG_jll", "FileIO", "FilePaths", "FixedPointNumbers", "Format", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "GridLayoutBase", "ImageBase", "ImageIO", "InteractiveUtils", "Interpolations", "IntervalSets", "InverseFunctions", "Isoband", "KernelDensity", "LaTeXStrings", "LinearAlgebra", "MacroTools", "Markdown", "MathTeXEngine", "Observables", "OffsetArrays", "PNGFiles", "Packing", "Pkg", "PlotUtils", "PolygonOps", "PrecompileTools", "Printf", "REPL", "Random", "RelocatableFolders", "Scratch", "ShaderAbstractions", "Showoff", "SignedDistanceFields", "SparseArrays", "Statistics", "StatsBase", "StatsFuns", "StructArrays", "TriplotBase", "UnicodeFun", "Unitful"]
-git-tree-sha1 = "68af66ec16af8b152309310251ecb4fbfe39869f"
+git-tree-sha1 = "0708c6a1f3cb18ba6482c4174058084c8d6deaf4"
 uuid = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
-version = "0.24.9"
+version = "0.24.10"
 
     [deps.Makie.extensions]
     MakieDynamicQuantitiesExt = "DynamicQuantities"
@@ -1406,9 +1406,9 @@ version = "1.1.10"
 
 [[deps.MbedTLS_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "ff69a2b1330bcb730b9ac1ab7dd680176f5896b8"
+git-tree-sha1 = "926c6af3a037c68d02596a44c22ec3595f5f760b"
 uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
-version = "2.28.1010+0"
+version = "2.28.6+0"
 
 [[deps.Measures]]
 git-tree-sha1 = "b513cedd20d9c914783d8ad83d08120702bf2c77"
@@ -1748,9 +1748,9 @@ version = "1.2.0"
 
 [[deps.PrecompileTools]]
 deps = ["Preferences"]
-git-tree-sha1 = "07a921781cab75691315adc645096ed5e370cb77"
+git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f"
 uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-version = "1.3.3"
+version = "1.2.1"
 
 [[deps.Preferences]]
 deps = ["TOML"]

From d5b169728297279ef70106b9cc6777939e042c5c Mon Sep 17 00:00:00 2001
From: Bart <bart@vandelint.net>
Date: Tue, 12 May 2026 17:52:52 +0200
Subject: [PATCH 6/8] Address copilot comments

---
 src/wing_geometry.jl | 21 ++++++++++++++++-----
 1 file changed, 16 insertions(+), 5 deletions(-)

diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl
index 5c36851d..c572f337 100644
--- a/src/wing_geometry.jl
+++ b/src/wing_geometry.jl
@@ -248,7 +248,8 @@ mutable struct Wing{P, T} <: AbstractWing{T}
 
     # Linear interpolation cache from unrefined sections to refined sections.
     # For refined section i: value = weight * unrefined[left_idx] + (1 - weight) * unrefined[left_idx + 1].
-    # weight == 1 at unrefined endpoints, so refined section values match unrefined endpoints exactly.
+    # First refined section is pinned to weight==1, left_idx==1 (takes unrefined[1]); last refined
+    # section is pinned to weight==0, left_idx==n_unref-1 (takes unrefined[end]).
     refined_section_left_idx::Vector{Int16}  # Length: n_panels + 1
     refined_section_weight::Vector{T}        # Length: n_panels + 1
 
@@ -361,11 +362,17 @@ take the unrefined endpoint values exactly. The panel-level `theta_dist` /
 downstream consumers (solver, body aerodynamics) see a per-panel value.
 
 # Arguments
-- `wing::Wing`: Wing to deform (must have non_deformed_sections).
+- `wing::Wing`: Wing to deform (must have non_deformed_sections, populated by
+  `refine!` for manual/YAML wings or by OBJ refinement for OBJ-based wings).
 - `theta_angles::AbstractVector`: Twist angles in radians, one per unrefined section.
   Pass `nothing` to leave twist unchanged.
 - `delta_angles::AbstractVector`: TE deflection angles in radians, one per unrefined
   section. Pass `nothing` to leave deflection unchanged.
+
+# Keyword arguments
+- `smooth`, `smooth_window`: accepted for backwards compatibility with callers of
+  `deform!`, but ignored here — the linear interpolation between unrefined sections
+  is already smooth, so no post-hoc smoothing is applied.
 """
 function unrefined_deform!(wing::Wing, theta_angles=nothing, delta_angles=nothing;
                            smooth=false, smooth_window=nothing)
@@ -373,7 +380,10 @@ function unrefined_deform!(wing::Wing, theta_angles=nothing, delta_angles=nothin
     isnothing(theta_angles) && isnothing(delta_angles) && return nothing
 
     if !can_deform
-        throw(ArgumentError("This Wing does not support deformation. Only OBJ-based wings created with ObjWing() can be deformed."))
+        throw(ArgumentError(
+            "This Wing has no non_deformed_sections to deform from. " *
+            "Call refine!(wing) (manual/YAML wings) or construct via ObjWing() " *
+            "before calling unrefined_deform!."))
     end
 
     n_unref = wing.n_unrefined_sections
@@ -996,8 +1006,9 @@ sections. For refined section i, the interpolated value is:
              (1 - weight[i]) * unrefined[left_idx[i] + 1]
 
 Positions are quarter-chord arc-length along the unrefined and refined sections.
-Endpoint refined sections always get `weight == 1` so they take the unrefined
-endpoint value exactly.
+The first refined section is pinned to `left_idx == 1`, `weight == 1` (returns
+`unrefined[1]` exactly) and the last refined section to `left_idx == n_unref - 1`,
+`weight == 0` (returns `unrefined[end]` exactly).
 """
 function compute_refined_section_interpolation!(wing::AbstractWing{T}) where {T}
     n_unref = length(wing.unrefined_sections)

From d431816209dcd1f0ba1927ac6c283bbabbfff7ab Mon Sep 17 00:00:00 2001
From: Bart <bart@vandelint.net>
Date: Tue, 12 May 2026 19:09:11 +0200
Subject: [PATCH 7/8] Use elliptic start

---
 test/solver/test_forwarddiff.jl | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/test/solver/test_forwarddiff.jl b/test/solver/test_forwarddiff.jl
index 600513d0..e17ba29f 100644
--- a/test/solver/test_forwarddiff.jl
+++ b/test/solver/test_forwarddiff.jl
@@ -20,8 +20,9 @@ using Test
     y0 = [va; omega]
 
     @testset "AutoForwardDiff matches AutoFiniteDiff (LOOP, INVISCID)" begin
-        # FD warm-starts otherwise bias the reference Jacobian.
-        solver = Solver(body_aero; use_gamma_prev=false, relaxation_factor=0.01)
+        solver = Solver(body_aero;
+            use_gamma_prev=false,
+            type_initial_gamma_distribution=ELLIPTIC)
 
         jac_fwd, _, fwd_converged = VortexStepMethod.linearize(
             solver, body_aero, y0;
@@ -50,8 +51,6 @@ using Test
     end
 
     @testset "AutoForwardDiff matches AutoFiniteDiff (LOOP, POLAR_MATRICES)" begin
-        # Exercises the mixed-eltype path: calculate_stall_angle_list! calls
-        # calculate_cl with Float64 alpha against Dual-typed panels.
         data_dir = joinpath(dirname(dirname(@__DIR__)), "data", "ram_air_kite")
         body_path = joinpath(tempdir(), "ram_air_kite_body.obj")
         foil_path = joinpath(tempdir(), "ram_air_kite_foil.dat")
@@ -74,7 +73,6 @@ using Test
             rtol=1e-7,
             solver_type=LOOP,
             use_gamma_prev=false,
-            relaxation_factor=0.01,
         )
 
         v_a = 15.0

From 8f256f832802f1cf697fd5a3d734106118251caa Mon Sep 17 00:00:00 2001
From: Bart <bart@vandelint.net>
Date: Tue, 12 May 2026 19:15:47 +0200
Subject: [PATCH 8/8] Use smooth sqrt

---
 src/solver.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/solver.jl b/src/solver.jl
index ede7c33b..4cd4bc79 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -656,6 +656,8 @@ function solve_base!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma
     nothing
 end
 
+@inline smooth_sqrt(x) = sqrt(x + 1e-30)
+
 @inline function update_gamma_candidate!(
     gamma_out,
     gamma_in,
@@ -718,11 +720,11 @@ end
     solver.lr.alpha_dist .= atan.(v_normal_array, v_tangential_array)
 
     @inbounds for i in 1:n_panels
-        solver.lr.v_a_dist[i] = sqrt(
+        solver.lr.v_a_dist[i] = smooth_sqrt(
             relative_velocity_crossz[i,1]^2 +
             relative_velocity_crossz[i,2]^2 +
             relative_velocity_crossz[i,3]^2)
-        va_magw_array[i] = sqrt(
+        va_magw_array[i] = smooth_sqrt(
             v_acrossz_array[i,1]^2 +
             v_acrossz_array[i,2]^2 +
             v_acrossz_array[i,3]^2)