From b7a47bf899f2b4a8f40bf939454ac5ae22847f90 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Mon, 14 Sep 2026 18:52:44 +0200 Subject: [PATCH 1/8] Refuse a contour XFoil's panel code cannot take A self-crossing contour can reach a bare Fortran STOP inside XFoil, which ends the Julia process with exit code 0 and no exception. `analyze_sweep` now checks the contour before `set_coordinates` and throws `ArgumentError` instead, so a section XFoil has no solution for is one failed sweep rather than the end of the session. The same check covers a contour with more nodes than XFoil's panel arrays hold: ABCOPY refused it and left the previously loaded airfoil in place, so every angle came back solved on the wrong shape. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_018rHHzpuPmi8pBjVk3r5QZD --- CHANGELOG.md | 11 ++++++ docs/src/private_functions.md | 8 +++++ src/airfoil_aero/airfoil_solvers/common.jl | 35 +++++++++++++++++++ .../airfoil_solvers/xfoil_solver.jl | 21 +++++++++++ test/airfoil_aero/test_airfoil_aero.jl | 22 +++++++++++- 5 files changed, 96 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 532b0816..dd2f327d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,16 @@ # Changelog +## Unreleased + +### Fixed + +- `analyze_sweep(::XFoilSolver, ...)` throws `ArgumentError` for a contour XFoil's + panel code cannot take, rather than handing it to the Fortran. A self-crossing + contour could reach a bare `STOP` there, ending the Julia process with exit code 0 + and no exception; a contour with more nodes than XFoil's panel arrays hold was + refused by XFoil, which then solved every angle on whichever airfoil was loaded + before and returned it as this one. + ## VortexStepMethod v5.1.1 2026-09-12 ### Fixed diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index d7807c42..35cc72a8 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -167,6 +167,14 @@ resample_arc smoothed_curvature ``` +### Section solver preconditions +```@docs +validate_xfoil_contour +crossing_panels +segments_cross +side_of_line +``` + ### NeuralFoil network ```@docs load_neuralfoil_model diff --git a/src/airfoil_aero/airfoil_solvers/common.jl b/src/airfoil_aero/airfoil_solvers/common.jl index 60c79409..b4e7725e 100644 --- a/src/airfoil_aero/airfoil_solvers/common.jl +++ b/src/airfoil_aero/airfoil_solvers/common.jl @@ -116,6 +116,41 @@ function deform_section(x, y, delta; crease_frac=0.9, thickness_frac=1.0, return DeformedSection(kulfan, xd, yd) end +""" + side_of_line(a, b, p) -> Float64 + +Twice the signed area of the triangle `a`, `b`, `p`: positive with `p` left of the +line from `a` to `b`, negative right of it, zero on it. +""" +side_of_line(a, b, p) = (b[1] - a[1]) * (p[2] - a[2]) - (b[2] - a[2]) * (p[1] - a[1]) + +""" + segments_cross(p, q, r, s) -> Bool + +Whether the segments `p`-`q` and `r`-`s` cross properly, each strictly separating +the other's endpoints. Touching at an endpoint or lying along each other does not +count. +""" +segments_cross(p, q, r, s) = + side_of_line(p, q, r) * side_of_line(p, q, s) < 0 && + side_of_line(r, s, p) * side_of_line(r, s, q) < 0 + +""" + crossing_panels(x, y) -> Tuple{Int,Int} or nothing + +The first pair of non-neighbouring panels of the closed contour `(x, y)` that +cross, or `nothing` when the contour is a simple closed curve. +""" +function crossing_panels(x, y) + nodes = collect(zip(x, y)) + last_panel = length(nodes) - 1 + for i in 1:last_panel, j in (i + 2):last_panel + i == 1 && j == last_panel && continue + segments_cross(nodes[i], nodes[i+1], nodes[j], nodes[j+1]) && return (i, j) + end + return nothing +end + """ analyze_sweep(solver, def, alpha_range, Re) -> Vector{SectionSolution} diff --git a/src/airfoil_aero/airfoil_solvers/xfoil_solver.jl b/src/airfoil_aero/airfoil_solvers/xfoil_solver.jl index 28ff516e..a5b9204c 100644 --- a/src/airfoil_aero/airfoil_solvers/xfoil_solver.jl +++ b/src/airfoil_aero/airfoil_solvers/xfoil_solver.jl @@ -28,6 +28,23 @@ runs (`ncrit=9`, `max_iter=100`, incompressible) so the two backends are compara repanel::Bool = false end +""" + validate_xfoil_contour(def::DeformedSection) + +Enforce what XFoil's panel code needs of `def`'s coordinates: no more nodes than its +panel arrays hold, and no two panels crossing. Throws `ArgumentError` on a violation. +""" +function validate_xfoil_contour(def::DeformedSection) + max_nodes = Xfoil.IQX - 5 + length(def.x) <= max_nodes || throw(ArgumentError( + "XFoil holds $max_nodes panel nodes; this contour has $(length(def.x)).")) + crossing = crossing_panels(def.x, def.y) + isnothing(crossing) || throw(ArgumentError( + "XFoil needs a contour that does not cross itself; panels $(crossing[1]) " * + "and $(crossing[2]) of this one do.")) + return nothing +end + """ analyze_sweep(solver::XFoilSolver, def, alpha_range, Re) -> Vector{SectionSolution} @@ -37,8 +54,12 @@ reinit at each side for convergence. Each converged angle reads the surface pres (`Xfoil.cpdump`) and the boundary layer (`Xfoil.bldump`, giving `cf` and the node coordinates) at the same panel nodes. Non-converged angles yield empty node arrays and `NaN` confidence. + +Throws `ArgumentError` for a contour XFoil's panel code cannot take, see +[`validate_xfoil_contour`](@ref). """ function analyze_sweep(solver::XFoilSolver, def::DeformedSection, alpha_range, Re) + validate_xfoil_contour(def) Xfoil.set_coordinates(def.x, def.y) solver.repanel && Xfoil.pane(npan=solver.npan) sols = Vector{SectionSolution}(undef, length(alpha_range)) diff --git a/test/airfoil_aero/test_airfoil_aero.jl b/test/airfoil_aero/test_airfoil_aero.jl index ac4b9363..90700a41 100644 --- a/test/airfoil_aero/test_airfoil_aero.jl +++ b/test/airfoil_aero/test_airfoil_aero.jl @@ -4,7 +4,8 @@ import VortexStepMethod using VortexStepMethod.AirfoilAero: KulfanParameters, LeastSquaresFit, ShrinkWrap, shrink_wrap, fit_kulfan_parameters, kulfan_to_coordinates, neuralfoil_aero, class_function, bernstein_basis, - leading_edge_basis, normalize_airfoil + leading_edge_basis, normalize_airfoil, crossing_panels, + DeformedSection, XFoilSolver, analyze_sweep, Xfoil using VortexStepMethod: SectionAero, section_surface, read_section_aero using VortexStepMethod.AirfoilAero: write_section_aero @@ -299,3 +300,22 @@ end @test maximum(abs, collect(extrema(written_y)) .- collect(extrema(fitted_y))) < 1e-4 end + +@testset "XFoil refuses a contour it has no solution for" begin + clean = KulfanParameters(fill(0.15, 8), fill(-0.15, 8), 0.0, 0.0) + x, y = collect.(kulfan_to_coordinates(clean; n_points=60)) + alphas = deg2rad.([0.0]) + @test isnothing(crossing_panels(x, y)) + + # the upper surface driven through the lower one over a stretch of the chord + folded = copy(y) + folded[20:40] .= -3 .* folded[20:40] + @test !isnothing(crossing_panels(x, folded)) + @test_throws ArgumentError analyze_sweep(XFoilSolver(), + DeformedSection(clean, x, folded), alphas, 1e6) + + crowded_x, crowded_y = collect.(kulfan_to_coordinates(clean; n_points=200)) + @test length(crowded_x) > Xfoil.IQX - 5 + @test_throws ArgumentError analyze_sweep(XFoilSolver(), + DeformedSection(clean, crowded_x, crowded_y), alphas, 1e6) +end From d1ffa0c03f094a49905755b380a3f4b92fe63ada Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Sun, 20 Sep 2026 22:08:40 +0200 Subject: [PATCH 2/8] Test the panel that closes the contour back to node 1 Co-Authored-By: Claude Opus 5 --- src/airfoil_aero/airfoil_solvers/common.jl | 9 +++++---- test/airfoil_aero/test_airfoil_aero.jl | 15 ++++++++++++--- 2 files changed, 17 insertions(+), 7 deletions(-) diff --git a/src/airfoil_aero/airfoil_solvers/common.jl b/src/airfoil_aero/airfoil_solvers/common.jl index b4e7725e..8d713e94 100644 --- a/src/airfoil_aero/airfoil_solvers/common.jl +++ b/src/airfoil_aero/airfoil_solvers/common.jl @@ -117,7 +117,7 @@ function deform_section(x, y, delta; crease_frac=0.9, thickness_frac=1.0, end """ - side_of_line(a, b, p) -> Float64 + side_of_line(a, b, p) Twice the signed area of the triangle `a`, `b`, `p`: positive with `p` left of the line from `a` to `b`, negative right of it, zero on it. @@ -138,14 +138,15 @@ segments_cross(p, q, r, s) = """ crossing_panels(x, y) -> Tuple{Int,Int} or nothing -The first pair of non-neighbouring panels of the closed contour `(x, y)` that -cross, or `nothing` when the contour is a simple closed curve. +The first pair of non-neighbouring panels of the contour `(x, y)` that cross, or +`nothing` when the contour is a simple closed curve. A contour that does not end on its +first node is closed by a panel back to it, which carries the highest panel index. """ function crossing_panels(x, y) nodes = collect(zip(x, y)) + last(nodes) == first(nodes) || push!(nodes, first(nodes)) last_panel = length(nodes) - 1 for i in 1:last_panel, j in (i + 2):last_panel - i == 1 && j == last_panel && continue segments_cross(nodes[i], nodes[i+1], nodes[j], nodes[j+1]) && return (i, j) end return nothing diff --git a/test/airfoil_aero/test_airfoil_aero.jl b/test/airfoil_aero/test_airfoil_aero.jl index f6600ca3..37cedff9 100644 --- a/test/airfoil_aero/test_airfoil_aero.jl +++ b/test/airfoil_aero/test_airfoil_aero.jl @@ -5,7 +5,8 @@ using VortexStepMethod.AirfoilAero: KulfanParameters, LeastSquaresFit, ShrinkWra shrink_wrap, fit_kulfan_parameters, kulfan_to_coordinates, neuralfoil_aero, class_function, bernstein_basis, leading_edge_basis, normalize_airfoil, crossing_panels, - DeformedSection, XFoilSolver, analyze_sweep, Xfoil + validate_xfoil_contour, DeformedSection, XFoilSolver, + analyze_sweep, Xfoil using VortexStepMethod: SectionAero, section_surface, read_section_aero using VortexStepMethod.AirfoilAero: write_section_aero @@ -326,11 +327,11 @@ end collect(extrema(fitted_y))) < 1e-4 end -@testset "XFoil refuses a contour it has no solution for" begin +@testset "a contour XFoil has no solution for is refused before XFoil sees it" begin clean = KulfanParameters(fill(0.15, 8), fill(-0.15, 8), 0.0, 0.0) x, y = collect.(kulfan_to_coordinates(clean; n_points=60)) alphas = deg2rad.([0.0]) - @test isnothing(crossing_panels(x, y)) + @test isnothing(validate_xfoil_contour(DeformedSection(clean, x, y))) # the upper surface driven through the lower one over a stretch of the chord folded = copy(y) @@ -339,6 +340,14 @@ end @test_throws ArgumentError analyze_sweep(XFoilSolver(), DeformedSection(clean, x, folded), alphas, 1e6) + gapped = KulfanParameters(fill(0.15, 8), fill(-0.15, 8), 0.0, 0.02) + gapped_x, gapped_y = collect.(kulfan_to_coordinates(gapped; n_points=60)) + @test (gapped_x[1], gapped_y[1]) != (gapped_x[end], gapped_y[end]) + @test isnothing(validate_xfoil_contour(DeformedSection(gapped, gapped_x, gapped_y))) + + # a contour whose only crossing is the panel closing it back to node 1 + @test crossing_panels([0.0, 1.0, 1.0, 2.0], [0.0, 2.0, -2.0, 1.0]) == (2, 4) + crowded_x, crowded_y = collect.(kulfan_to_coordinates(clean; n_points=200)) @test length(crowded_x) > Xfoil.IQX - 5 @test_throws ArgumentError analyze_sweep(XFoilSolver(), From 85652b397dbd454cb0fa361550a92ec8818aadd0 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Mon, 21 Sep 2026 13:14:55 +0200 Subject: [PATCH 3/8] Re-wrap a deformed section with the ball it was first wrapped with, on its polygon deform_section took a fixed ShrinkWrap(clearance=0.0) for the re-wrap, so a section wrapped at min_concave_radius 0.2 was re-wrapped at 0.02. shrink_wrap rolled that ball on the input's nodes only, and across the long panels of a thin wrapped canopy it touched the other skin and came back crossing itself. deform_section now takes the wrap the section came from and re-wraps with its ball at zero clearance; the wrap is threaded through generate_airfoils, generate_airfoil_aero, generate_aero_matrices, generate_polar_from_coordinates and the Makie previews. shrink_wrap densifies a closed input to edges of at most min(0.01, radius/2) so the ball rolls on the polygon, and warns when the contour it returns crosses itself. densify_contour moves from ObjAdapter to the shrink wrap it now serves. Co-Authored-By: Claude Opus 5 --- docs/src/private_functions.md | 2 +- ext/VortexStepMethodMakieExt.jl | 6 ++- src/airfoil_aero/airfoil_io.jl | 12 +++--- src/airfoil_aero/airfoil_solvers/common.jl | 18 ++++----- src/airfoil_aero/geometry_gen.jl | 7 ++-- src/airfoil_aero/polar_export.jl | 10 +++-- src/airfoil_aero/section_aero_gen.jl | 29 +++++++-------- src/airfoil_aero/shrink_wrap.jl | 43 ++++++++++++++++++---- src/obj_adapter/obj_slice.jl | 22 +---------- src/obj_adapter/obj_to_yaml.jl | 2 +- src/surfplan_adapter/SurfplanAdapter.jl | 2 +- test/airfoil_aero/test_airfoil_aero.jl | 30 +++++++++++++++ 12 files changed, 113 insertions(+), 70 deletions(-) diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index a29c9f24..314a8789 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -169,6 +169,7 @@ push_arc! edge_normal pivot_contour largest_linking_gap +densify_contour enforce_min_spacing! resample_arc smoothed_curvature @@ -255,7 +256,6 @@ build_section contour_to_airfoil plane_contour_to_airfoil reorder_airfoil_selig -densify_contour create_interpolations find_circle_center_and_radius march_edges diff --git a/ext/VortexStepMethodMakieExt.jl b/ext/VortexStepMethodMakieExt.jl index 3cd2569a..c5c092cf 100644 --- a/ext/VortexStepMethodMakieExt.jl +++ b/ext/VortexStepMethodMakieExt.jl @@ -1542,7 +1542,8 @@ function fitted_airfoil_3d(s, wrap_method; delta=0.0, crease_frac=0.75) collect(pz ./ chord), wrap_method) maximum(abs, yf) > 1.0 && return nothing if !iszero(delta) - def = AirfoilAero.deform_section(xf, yf, deg2rad(delta); crease_frac) + def = AirfoilAero.deform_section(xf, yf, deg2rad(delta); crease_frac, + wrap_method) xf, yf = def.x, def.y end return reduce(hcat, [s.LE_point .+ x_af .* (x0 + xf[i] * chord) .+ z_af .* (yf[i] * chord) @@ -1672,7 +1673,8 @@ function ObjAdapter.plot_slices_3d(path::String; n_slices::Int=10, rotation=I, xf, yf = AirfoilAero.shrink_wrap(collect(Float64, s.x_airfoil), collect(Float64, s.y_airfoil), wrap_method) def = iszero(delta) ? nothing : - AirfoilAero.deform_section(xf, yf, deg2rad(delta); crease_frac) + AirfoilAero.deform_section(xf, yf, deg2rad(delta); crease_frac, + wrap_method) d2 = (; raw=Point2f.(s.x_airfoil, s.y_airfoil), fit=Point2f.(xf, yf), fit_kulfan=fit_pts(xf, yf), def=def === nothing ? Point2f[] : Point2f.(def.x, def.y), diff --git a/src/airfoil_aero/airfoil_io.jl b/src/airfoil_aero/airfoil_io.jl index 10cefff5..ef04bd3d 100644 --- a/src/airfoil_aero/airfoil_io.jl +++ b/src/airfoil_aero/airfoil_io.jl @@ -33,12 +33,13 @@ end """ generate_polar_from_coordinates(x, y, output_path; Re, alpha_range=-180:1:180, solver=NeuralFoilSolver(), delta_range=nothing, - crease_frac=0.75, dat_prefix=nothing) + crease_frac=0.75, dat_prefix=nothing, + wrap_method=ShrinkWrap()) Sweep `solver` over the airfoil coordinates `(x, y)` and write the polar CSV. XFoil uses the coordinates directly; NeuralFoil fits [`LeastSquaresFit`](@ref) Kulfan parameters ([`deform_section`](@ref)). Wrap a raw or open single-membrane slice with -[`shrink_wrap`](@ref) before calling this. Pass a [`NeuralFoilSolver`](@ref) or +[`shrink_wrap`](@ref) before calling this, and pass the same `wrap_method`. Pass a [`NeuralFoilSolver`](@ref) or [`XFoilSolver`](@ref) to pick the backend. With `delta_range === nothing` the sweep is over `alpha_range` only and written as a `POLAR_VECTORS` CSV (returns the `Vector{SectionSolution}`); pass a `delta_range` of trailing-edge deflections to sweep @@ -51,10 +52,11 @@ function generate_polar_from_coordinates(x::Vector, y::Vector, output_path::Stri Re::Real, alpha_range=-180:1:180, solver::AbstractAirfoilSolver=NeuralFoilSolver(), delta_range=nothing, crease_frac=0.75, - dat_prefix=nothing) + dat_prefix=nothing, + wrap_method::ShrinkWrap=ShrinkWrap()) alphas = deg2rad.(collect(Float64, alpha_range)) if delta_range === nothing - def = deform_section(x, y, 0.0) + def = deform_section(x, y, 0.0; wrap_method) sols = analyze_sweep(solver, def, alphas, Re) write_polar_csv(output_path, sols) return sols @@ -64,7 +66,7 @@ function generate_polar_from_coordinates(x::Vector, y::Vector, output_path::Stri (d, xd, yd) -> write_dat("$(dat_prefix)_$(delta_suffix(d)).dat", "deflection", xd, yd) cl, cd, cm = generate_aero_matrices(solver, x, y; - alpha_range=alphas, delta_range=deltas, Re, crease_frac, on_deform) + alpha_range=alphas, delta_range=deltas, Re, crease_frac, on_deform, wrap_method) write_polar_matrix_csv(output_path, alphas, deltas, cl, cd, cm) return (cl, cd, cm) end diff --git a/src/airfoil_aero/airfoil_solvers/common.jl b/src/airfoil_aero/airfoil_solvers/common.jl index 8d713e94..5eb0fc8e 100644 --- a/src/airfoil_aero/airfoil_solvers/common.jl +++ b/src/airfoil_aero/airfoil_solvers/common.jl @@ -85,8 +85,7 @@ end """ deform_section(x, y, delta; crease_frac=0.9, thickness_frac=1.0, - flip_thickness_neg=true, - wrap_method=ShrinkWrap(clearance=0.0)) + flip_thickness_neg=true, wrap_method=ShrinkWrap()) -> DeformedSection Deform the airfoil coordinates `(x, y)` by trailing-edge deflection `delta` (radians) @@ -96,22 +95,21 @@ into clean cosine panels and fit [`LeastSquaresFit`](@ref) Kulfan parameters to XFoil consumes the coordinates directly, NeuralFoil the Kulfan parameters. `flip_thickness_neg` folds a soft membrane about its lower surface for negative `delta`. -The re-wrap uses zero clearance (it hugs the deflected shape at `min_clearance`); -the rolling-ball wrap bridges the crease with a `min_concave_radius` fillet instead -of the overlapping panels that XFoil's own repaneling can hit there. The wrap runs -for every `delta` including `0`, so all deflections share the same node count -(`2·n_points - 1`). +`wrap_method` is the wrap `(x, y)` came from: the re-wrap rolls the same ball at zero +clearance, bridging the crease with a `min_concave_radius` fillet. It runs for every +`delta` including `0`, so all deflections share the same node count (`2·n_points - 1`). """ function deform_section(x, y, delta; crease_frac=0.9, thickness_frac=1.0, - flip_thickness_neg=true, - wrap_method::ShrinkWrap=ShrinkWrap(clearance=0.0)) + flip_thickness_neg=true, wrap_method::ShrinkWrap=ShrinkWrap()) xd, yd = collect(float.(x)), collect(float.(y)) if !iszero(delta) pivot = flip_thickness_neg && delta < 0 ? 1 - thickness_frac : thickness_frac lower, upper = get_lower_upper(xd, yd, crease_frac) turn_trailing_edge!(delta, xd, yd, lower, upper, crease_frac; thickness_frac=pivot) end - xd, yd = shrink_wrap(xd, yd, wrap_method) + rewrap = ShrinkWrap(0.0, wrap_method.min_concave_radius, wrap_method.min_clearance, + wrap_method.n_points, wrap_method.curvature_weight) + xd, yd = shrink_wrap(xd, yd, rewrap) kulfan = fit_kulfan_parameters(xd, yd, LeastSquaresFit()) return DeformedSection(kulfan, xd, yd) end diff --git a/src/airfoil_aero/geometry_gen.jl b/src/airfoil_aero/geometry_gen.jl index 6244664d..6c00038d 100644 --- a/src/airfoil_aero/geometry_gen.jl +++ b/src/airfoil_aero/geometry_gen.jl @@ -2,7 +2,7 @@ generate_airfoils(airfoils, output_dir; Re, alpha_range=-180:1:180, delta_range=nothing, aero_solver=NeuralFoilSolver(), reuse_valid_airfoils=true, crease_frac=0.75, verbose=true, - table_format=:csv) -> (airfoil_rows, ok) + table_format=:csv, wrap_method=ShrinkWrap()) -> (airfoil_rows, ok) Run the 2D solver over a set of already-shrink-wrapped airfoils and write the per section files each geometry route references. Shared by the `.obj` and Surfplan @@ -12,6 +12,7 @@ placement; this writes the surface pressure/friction tables, polars and airfoil `airfoils` is a vector of `(; id, x_fit, y_fit, x_raw, y_raw)`: `x_fit`/`y_fit` is the wrapped airfoil the solver analyses; `x_raw`/`y_raw` the raw points it enclosed. +`wrap_method` is the [`ShrinkWrap`](@ref) that produced `x_fit`/`y_fit`. Writes into `output_dir` (indexed by `id`), one directory per file kind: `airfoils/{id}.dat` (wrapped shape), `airfoils/{id}_{delta_suffix(δ)}.dat` (per @@ -28,7 +29,7 @@ function generate_airfoils(airfoils, output_dir::String; Re::Real, alpha_range=-180:1:180, delta_range=nothing, aero_solver::AbstractAirfoilSolver=NeuralFoilSolver(), reuse_valid_airfoils::Bool=true, crease_frac=0.75, verbose::Bool=true, - table_format::Symbol=:csv) + table_format::Symbol=:csv, wrap_method::ShrinkWrap=ShrinkWrap()) mkpath(joinpath(output_dir, "airfoils")) mkpath(joinpath(output_dir, "polars")) mkpath(joinpath(output_dir, "pressure")) @@ -45,7 +46,7 @@ function generate_airfoils(airfoils, output_dir::String; aero, sols = generate_airfoil_aero(aero_solver, fit_kulfan_parameters(af.x_fit, af.y_fit, LeastSquaresFit()); alpha_range=alphas, delta_range=deltas, - reynolds_number=Float64(Re), crease_frac) + reynolds_number=Float64(Re), crease_frac, wrap_method) clvals = collect(sol.cl for sol in sols[1]) all(isnan, clvals) && error("solver produced no converged points") if isnothing(delta_range) diff --git a/src/airfoil_aero/polar_export.jl b/src/airfoil_aero/polar_export.jl index 40b15b02..428c6f64 100644 --- a/src/airfoil_aero/polar_export.jl +++ b/src/airfoil_aero/polar_export.jl @@ -21,7 +21,8 @@ end """ generate_aero_matrices(solver, x, y; alpha_range, delta_range, Re, - crease_frac=0.75, remove_nan=true, on_deform=nothing) + crease_frac=0.75, remove_nan=true, on_deform=nothing, + wrap_method=ShrinkWrap()) -> (cl, cd, cm) Build `(alpha × delta)` coefficient matrices for a base airfoil given as coordinates @@ -30,18 +31,19 @@ deflected shape is then swept over `alpha_range` (radians) with `solver` — any [`AbstractAirfoilSolver`](@ref), so this works identically for XFoil and NeuralFoil. `Re` is the Reynolds number. With `remove_nan` the (non-converged) `NaN` entries are interpolated away. `on_deform(delta, x, y)`, if given, is called with each deflected -shape's coordinates (e.g. to write a per-deflection `.dat`). +shape's coordinates (e.g. to write a per-deflection `.dat`). `wrap_method` is the +[`ShrinkWrap`](@ref) `(x, y)` was wrapped with. """ function generate_aero_matrices(solver::AbstractAirfoilSolver, x, y; alpha_range, delta_range, Re, crease_frac=0.75, remove_nan=true, - on_deform=nothing) + on_deform=nothing, wrap_method::ShrinkWrap=ShrinkWrap()) na, nd = length(alpha_range), length(delta_range) cl = fill(NaN, na, nd) cd = fill(NaN, na, nd) cm = fill(NaN, na, nd) alphas = collect(Float64, alpha_range) for (j, delta) in enumerate(delta_range) - def = deform_section(x, y, delta; crease_frac) + def = deform_section(x, y, delta; crease_frac, wrap_method) on_deform === nothing || on_deform(delta, def.x, def.y) sols = analyze_sweep(solver, def, alphas, Re) for (i, s) in enumerate(sols) diff --git a/src/airfoil_aero/section_aero_gen.jl b/src/airfoil_aero/section_aero_gen.jl index 9c515f88..2e0aff5e 100644 --- a/src/airfoil_aero/section_aero_gen.jl +++ b/src/airfoil_aero/section_aero_gen.jl @@ -1,21 +1,24 @@ """ generate_airfoil_aero(solver, base; alpha_range, delta_range, reynolds_number, - crease_frac=0.9, remove_nan=true) -> (SectionAero, sols) + crease_frac=0.9, remove_nan=true, wrap_method=ShrinkWrap()) + -> (SectionAero, sols) Run **one** solver sweep of a base airfoil (Kulfan) over the `(alpha, delta)` grid (radians) and return both the [`SectionAero`](@ref) (contour + `Cp` + `cf` per node) and the raw `sols::Vector{Vector{SectionSolution}}` (one inner vector per delta). The `sols` also carry `cl/cd/cm`, so a caller can write the polar from the same sweep — this is how `obj_to_yaml` avoids a second sweep. Non-converged points stay `NaN` and, when -`remove_nan`, are filled per node with `interpolate_matrix_nans!`. +`remove_nan`, are filled per node with `interpolate_matrix_nans!`. `wrap_method` is the +[`ShrinkWrap`](@ref) the base airfoil was wrapped with ([`deform_section`](@ref)). """ function generate_airfoil_aero(solver::AbstractAirfoilSolver, base::KulfanParameters; - alpha_range, delta_range, reynolds_number, crease_frac=0.9, remove_nan=true) + alpha_range, delta_range, reynolds_number, crease_frac=0.9, remove_nan=true, + wrap_method::ShrinkWrap=ShrinkWrap()) x0, y0 = kulfan_to_coordinates(base) n_alpha, n_delta = length(alpha_range), length(delta_range) sols = Vector{Vector{SectionSolution}}(undef, n_delta) for (jd, delta) in enumerate(delta_range) - def = deform_section(x0, y0, delta; crease_frac) + def = deform_section(x0, y0, delta; crease_frac, wrap_method) sols[jd] = analyze_sweep(solver, def, alpha_range, reynolds_number) end @@ -52,14 +55,14 @@ end """ generate_airfoil_aero(solver, x::Vector, y::Vector; kwargs...) -> (SectionAero, sols) -Convenience: [`shrink_wrap`](@ref) the coordinates and fit base Kulfan parameters -([`LeastSquaresFit`](@ref)) first. +Convenience: [`shrink_wrap`](@ref) the coordinates with `wrap_method` (keyword, +default `ShrinkWrap()`) and fit base Kulfan parameters ([`LeastSquaresFit`](@ref)) first. """ function generate_airfoil_aero(solver::AbstractAirfoilSolver, x::Vector, y::Vector; - kwargs...) - xw, yw = shrink_wrap(x, y, ShrinkWrap()) + wrap_method::ShrinkWrap=ShrinkWrap(), kwargs...) + xw, yw = shrink_wrap(x, y, wrap_method) return generate_airfoil_aero(solver, fit_kulfan_parameters(xw, yw, LeastSquaresFit()); - kwargs...) + wrap_method, kwargs...) end """ @@ -70,12 +73,8 @@ Just the [`SectionAero`](@ref) from [`generate_airfoil_aero`](@ref) (drops the r generate_section_aero(solver::AbstractAirfoilSolver, base::KulfanParameters; kwargs...) = generate_airfoil_aero(solver, base; kwargs...)[1] -function generate_section_aero(solver::AbstractAirfoilSolver, x::Vector, y::Vector; - kwargs...) - xw, yw = shrink_wrap(x, y, ShrinkWrap()) - return generate_section_aero(solver, fit_kulfan_parameters(xw, yw, LeastSquaresFit()); - kwargs...) -end +generate_section_aero(solver::AbstractAirfoilSolver, x::Vector, y::Vector; kwargs...) = + generate_airfoil_aero(solver, x, y; kwargs...)[1] """ fill_node_nans!(grid, i) diff --git a/src/airfoil_aero/shrink_wrap.jl b/src/airfoil_aero/shrink_wrap.jl index 14d2bd95..8e0ca4a7 100644 --- a/src/airfoil_aero/shrink_wrap.jl +++ b/src/airfoil_aero/shrink_wrap.jl @@ -218,6 +218,25 @@ function largest_linking_gap(x, y) return sqrt(gap2) end +""" + densify_contour(contour, max_edge) -> Vector + +Insert evenly spaced points along the edges of the point sequence `contour` longer +than `max_edge`, at most 100 pieces per edge. +""" +function densify_contour(contour, max_edge) + out = eltype(contour)[] + for i in 1:length(contour)-1 + p1, p2 = contour[i], contour[i+1] + n = clamp(ceil(Int, norm(p2 .- p1) / max_edge), 1, 100) + for k in 0:n-1 + push!(out, p1 .+ (p2 .- p1) .* (k / n)) + end + end + push!(out, contour[end]) + return out +end + """ smoothed_curvature(node_arclength, turn, band) -> Vector{Float64} @@ -331,18 +350,23 @@ LE → TE lower), following [`ShrinkWrap`](@ref): the rolling ball (`min_concave_radius`) is pivoted around the cloud, its contact side offset outward by `clearance` ([`pivot_contour`](@ref)), and the resulting arcs are resampled to cosine panels in a curvature-weighted arclength measure. The first and last point -coincide at the -trailing edge (the TE cap is part of the contour). A closed loop (first and last -cloud points coincident) keeps its true `clearance`, so `clearance=0` hugs the input -and leaves a sharp trailing edge sharp; an open single-membrane cloud is floored at -`min_clearance`. The output stays in the -normalized frame of the input cloud (chord slightly longer than 1, nose apex near +coincide at the trailing edge (the TE cap is part of the contour). A closed loop +(first and last cloud points coincident) is wrapped as the polygon through its points +([`densify_contour`](@ref)) and keeps its true `clearance`, so `clearance=0` hugs the +input and leaves a sharp trailing edge sharp; an open single-membrane cloud is floored +at `min_clearance`. Warns when the wrapped contour crosses itself. The output stays in +the normalized frame of the input cloud (chord slightly longer than 1, nose apex near `x = -clearance`) and is ready to write as a `.dat` or fit with [`LeastSquaresFit`](@ref). """ function shrink_wrap(x, y, method::ShrinkWrap) xn, yn, _ = normalize_airfoil(collect(float.(x)), collect(float.(y))) closed = hypot(xn[end] - xn[1], yn[end] - yn[1]) < 0.02 + if closed + max_edge = min(0.01, method.min_concave_radius / 2) + nodes = densify_contour(collect(zip(xn, yn)), max_edge) + xn, yn = first.(nodes), last.(nodes) + end gap = closed ? method.clearance : max(method.clearance, method.min_clearance) linking = largest_linking_gap(xn, yn) ball = max(method.min_concave_radius, 1.01 * linking / 2) @@ -364,5 +388,10 @@ function shrink_wrap(x, y, method::ShrinkWrap) method.curvature_weight) xl, yl = resample_arc(vcat(px[le:end], px[1]), vcat(py[le:end], py[1]), method.n_points, method.curvature_weight) - return vcat(xu, xl[2:end]), vcat(yu, yl[2:end]) + xo, yo = vcat(xu, xl[2:end]), vcat(yu, yl[2:end]) + crossing = crossing_panels(xo, yo) + isnothing(crossing) || @warn "shrink_wrap: panels $(crossing) of the wrapped " * + "contour cross; the rolling ball ($(round(ball; sigdigits=3))) or the " * + "clearance ($(gap)) does not fit this cloud" + return xo, yo end diff --git a/src/obj_adapter/obj_slice.jl b/src/obj_adapter/obj_slice.jl index 729895f7..f1ac8952 100644 --- a/src/obj_adapter/obj_slice.jl +++ b/src/obj_adapter/obj_slice.jl @@ -348,26 +348,6 @@ function march_edges(vertices, faces; step) point=[r.point for r in rows], tangent=[r.tangent for r in rows], arclen) end -""" - densify_contour(contour, max_edge) -> Vector{Vector{Float64}} - -Insert evenly spaced points along contour edges longer than `max_edge`. Coarse mesh -triangles otherwise leave large hops in the slice cloud, which force the shrink -wrap's auto-raised rolling ball far up and over-smooth the wrapped airfoil. -""" -function densify_contour(contour, max_edge) - out = Vector{Float64}[] - for i in 1:length(contour)-1 - p1, p2 = contour[i], contour[i+1] - n = clamp(ceil(Int, norm(p2 .- p1) / max_edge), 1, 100) - for k in 0:n-1 - push!(out, p1 .+ (p2 .- p1) .* (k / n)) - end - end - push!(out, contour[end]) - return out -end - """ build_section(vertices, faces, le, te, point, tangent) -> section or nothing @@ -383,7 +363,7 @@ function build_section(vertices, faces, le, te, point, tangent) isempty(segments) && return nothing contour = order_segments_to_contour(segments) length(contour) < 5 && return nothing - contour = densify_contour(contour, 0.005 * norm(te .- le)) + contour = AirfoilAero.densify_contour(contour, 0.005 * norm(te .- le)) af = plane_contour_to_airfoil(contour, le, x_af, z_af) af === nothing && return nothing return (; LE_point=le, TE_point=te, span_dir=y_af, contour3d=contour, diff --git a/src/obj_adapter/obj_to_yaml.jl b/src/obj_adapter/obj_to_yaml.jl index 10dc038f..ed61e13a 100644 --- a/src/obj_adapter/obj_to_yaml.jl +++ b/src/obj_adapter/obj_to_yaml.jl @@ -236,7 +236,7 @@ function obj_to_yaml(obj_path::String, output_dir::String; x_raw = stations[j].xa, y_raw = stations[j].ya) for j in unique(ids)] airfoil_rows, ok = generate_airfoils(airfoils, output_dir; Re, alpha_range, delta_range, aero_solver, reuse_valid_airfoils, crease_frac, verbose, - table_format) + table_format, wrap_method) isempty(ok) && error("No section produced a valid polar in $obj_path") prefix_table_paths!(airfoil_rows, table_path_prefix(yaml_path, output_dir)) diff --git a/src/surfplan_adapter/SurfplanAdapter.jl b/src/surfplan_adapter/SurfplanAdapter.jl index a723daab..935c41cb 100644 --- a/src/surfplan_adapter/SurfplanAdapter.jl +++ b/src/surfplan_adapter/SurfplanAdapter.jl @@ -73,7 +73,7 @@ function surfplan_to_aero_yaml(adapter_dir::AbstractString, output_dir::Abstract end airfoil_rows, ok = generate_airfoils(airfoils, output_dir; Re, alpha_range, - delta_range, aero_solver, crease_frac, verbose, table_format) + delta_range, aero_solver, crease_frac, verbose, table_format, wrap_method) isempty(ok) && error("No airfoil produced a valid polar from $adapter_dir") remap(id) = id in ok ? id : ok[argmin(abs.(ok .- id))] diff --git a/test/airfoil_aero/test_airfoil_aero.jl b/test/airfoil_aero/test_airfoil_aero.jl index 37cedff9..413e85e8 100644 --- a/test/airfoil_aero/test_airfoil_aero.jl +++ b/test/airfoil_aero/test_airfoil_aero.jl @@ -275,6 +275,36 @@ end end end +@testset "re-wrapping a wrapped section keeps it a simple closed curve" begin + shoelace(x, y) = abs(sum(x[i] * y[mod1(i + 1, end)] - x[mod1(i + 1, end)] * y[i] + for i in eachindex(x))) / 2 + obj = joinpath(pkgdir(VortexStepMethod), "data", "TUDELFT_V3_KITE", "V3_25.obj") + vertices, faces = VortexStepMethod.ObjAdapter.read_faces(obj) + canopy = VortexStepMethod.ObjAdapter.perpendicular_sections(vertices, faces, 18; + n_bins=100) + for radius in (0.02, 0.2) + wrap = ShrinkWrap(clearance=0.0, min_concave_radius=radius) + for section in canopy, delta in deg2rad.((0.0, 5.0)) + xw, yw = shrink_wrap(section.x_airfoil, section.y_airfoil, wrap) + def = deform_section(xw, yw, delta; wrap_method=wrap) + @test isnothing(crossing_panels(def.x, def.y)) + end + end + + # a clearance-padded wrap is re-wrapped at zero clearance, neither padded again + # nor collapsed by the ball falling between its long panels + xw, yw = shrink_wrap(read_dat_coordinates(joinpath(@__DIR__, "data", + "test_airfoil.dat"))..., + ShrinkWrap()) + def = deform_section(xw, yw, 0.0; wrap_method=ShrinkWrap()) + @test isnothing(crossing_panels(def.x, def.y)) + @test shoelace(def.x, def.y) ≈ shoelace(xw, yw) rtol = 0.05 + + @test_logs (:warn, r"cross") match_mode=:any shrink_wrap(canopy[2].x_airfoil, + canopy[2].y_airfoil, + ShrinkWrap()) +end + @testset "generate_polar_from_coordinates POLAR_VECTORS sweep" begin x, y = read_dat_coordinates(joinpath(@__DIR__, "data", "test_airfoil.dat")) csv = joinpath(mktempdir(), "polar.csv") From f016003a73f2fc6f092d04a90d0bc29be0063965 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Mon, 21 Sep 2026 13:18:02 +0200 Subject: [PATCH 4/8] Log the re-wrap fix and name the wrap generate_polar_from_coordinates takes Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 6 ++++++ src/airfoil_aero/airfoil_io.jl | 5 +++-- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 53b47333..f9bf11d4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,6 +32,12 @@ ### Fixed +- `deform_section` re-wraps a section with the rolling ball it was first wrapped with, + passed as `wrap_method` (also taken by `generate_airfoils`, `generate_airfoil_aero`, + `generate_aero_matrices` and `generate_polar_from_coordinates`), at zero clearance. + `shrink_wrap` rolls the ball on the polygon of a closed input rather than its nodes + and warns when the contour it returns crosses itself. A thin wrapped section came back + with its surfaces crossing. - `solve!` and `solve` throw a `DimensionMismatch` naming both sizes for a `body_aero` whose panel or unrefined-section count differs from the solver's, where they failed on a broadcast partway through or silently left section results at zero. diff --git a/src/airfoil_aero/airfoil_io.jl b/src/airfoil_aero/airfoil_io.jl index ef04bd3d..6792500e 100644 --- a/src/airfoil_aero/airfoil_io.jl +++ b/src/airfoil_aero/airfoil_io.jl @@ -39,8 +39,9 @@ end Sweep `solver` over the airfoil coordinates `(x, y)` and write the polar CSV. XFoil uses the coordinates directly; NeuralFoil fits [`LeastSquaresFit`](@ref) Kulfan parameters ([`deform_section`](@ref)). Wrap a raw or open single-membrane slice with -[`shrink_wrap`](@ref) before calling this, and pass the same `wrap_method`. Pass a [`NeuralFoilSolver`](@ref) or -[`XFoilSolver`](@ref) to pick the backend. With `delta_range === nothing` the sweep is +[`shrink_wrap`](@ref) before calling this and pass that [`ShrinkWrap`](@ref) as +`wrap_method`. Pass a [`NeuralFoilSolver`](@ref) or [`XFoilSolver`](@ref) to pick the +backend. With `delta_range === nothing` the sweep is over `alpha_range` only and written as a `POLAR_VECTORS` CSV (returns the `Vector{SectionSolution}`); pass a `delta_range` of trailing-edge deflections to sweep `(alpha, delta)` and write a long-format `POLAR_MATRICES` CSV (returns the `(cl, cd, From 7c5347447743560ddb20263f77086ce9ec440a77 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Mon, 21 Sep 2026 13:25:11 +0200 Subject: [PATCH 5/8] Name the re-wrap's fields, state the densify cap, test the warning on a built input The crossing warning is now tested on a wavy membrane wrapped at a clearance its waves cannot take, with a clean wrap of the same curve as the control, instead of on a V3 slice whose default-settings crossing is #363. Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 8 +++++--- src/airfoil_aero/airfoil_solvers/common.jl | 7 +++++-- src/airfoil_aero/shrink_wrap.jl | 5 +++-- test/airfoil_aero/test_airfoil_aero.jl | 19 ++++++++++++------- 4 files changed, 25 insertions(+), 14 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f9bf11d4..fd4963ed 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -29,15 +29,17 @@ - `Solver(body_aero; kwargs...)` and `Solver(body_aero, settings)` are deprecated and warn on use; build the solver with `Solver(settings)` or `Solver(n_panels, n_unrefined_sections)` instead. +- `shrink_wrap` splits the edges of a closed input to at most + `min(0.01, min_concave_radius/2)` chord before wrapping it, so the deflected sections + and polars `obj_to_yaml` generates from Kulfan contours move slightly. ### Fixed - `deform_section` re-wraps a section with the rolling ball it was first wrapped with, passed as `wrap_method` (also taken by `generate_airfoils`, `generate_airfoil_aero`, `generate_aero_matrices` and `generate_polar_from_coordinates`), at zero clearance. - `shrink_wrap` rolls the ball on the polygon of a closed input rather than its nodes - and warns when the contour it returns crosses itself. A thin wrapped section came back - with its surfaces crossing. + `shrink_wrap` warns when the contour it returns crosses itself. A thin wrapped section + came back with its surfaces crossing. - `solve!` and `solve` throw a `DimensionMismatch` naming both sizes for a `body_aero` whose panel or unrefined-section count differs from the solver's, where they failed on a broadcast partway through or silently left section results at zero. diff --git a/src/airfoil_aero/airfoil_solvers/common.jl b/src/airfoil_aero/airfoil_solvers/common.jl index 5eb0fc8e..b53654c3 100644 --- a/src/airfoil_aero/airfoil_solvers/common.jl +++ b/src/airfoil_aero/airfoil_solvers/common.jl @@ -107,8 +107,11 @@ function deform_section(x, y, delta; crease_frac=0.9, thickness_frac=1.0, lower, upper = get_lower_upper(xd, yd, crease_frac) turn_trailing_edge!(delta, xd, yd, lower, upper, crease_frac; thickness_frac=pivot) end - rewrap = ShrinkWrap(0.0, wrap_method.min_concave_radius, wrap_method.min_clearance, - wrap_method.n_points, wrap_method.curvature_weight) + rewrap = ShrinkWrap(; clearance=0.0, + min_concave_radius=wrap_method.min_concave_radius, + min_clearance=wrap_method.min_clearance, + n_points=wrap_method.n_points, + curvature_weight=wrap_method.curvature_weight) xd, yd = shrink_wrap(xd, yd, rewrap) kulfan = fit_kulfan_parameters(xd, yd, LeastSquaresFit()) return DeformedSection(kulfan, xd, yd) diff --git a/src/airfoil_aero/shrink_wrap.jl b/src/airfoil_aero/shrink_wrap.jl index 8e0ca4a7..31493137 100644 --- a/src/airfoil_aero/shrink_wrap.jl +++ b/src/airfoil_aero/shrink_wrap.jl @@ -351,8 +351,9 @@ LE → TE lower), following [`ShrinkWrap`](@ref): the rolling ball by `clearance` ([`pivot_contour`](@ref)), and the resulting arcs are resampled to cosine panels in a curvature-weighted arclength measure. The first and last point coincide at the trailing edge (the TE cap is part of the contour). A closed loop -(first and last cloud points coincident) is wrapped as the polygon through its points -([`densify_contour`](@ref)) and keeps its true `clearance`, so `clearance=0` hugs the +(first and last cloud points coincident) is wrapped as the polygon through its points, +its edges split to at most `min(0.01, min_concave_radius/2)` chord +([`densify_contour`](@ref)), and keeps its true `clearance`, so `clearance=0` hugs the input and leaves a sharp trailing edge sharp; an open single-membrane cloud is floored at `min_clearance`. Warns when the wrapped contour crosses itself. The output stays in the normalized frame of the input cloud (chord slightly longer than 1, nose apex near diff --git a/test/airfoil_aero/test_airfoil_aero.jl b/test/airfoil_aero/test_airfoil_aero.jl index 413e85e8..b2f0286d 100644 --- a/test/airfoil_aero/test_airfoil_aero.jl +++ b/test/airfoil_aero/test_airfoil_aero.jl @@ -1,4 +1,5 @@ using Test +using Logging using VortexStepMethod.AirfoilAero import VortexStepMethod using VortexStepMethod.AirfoilAero: KulfanParameters, LeastSquaresFit, ShrinkWrap, @@ -276,8 +277,6 @@ end end @testset "re-wrapping a wrapped section keeps it a simple closed curve" begin - shoelace(x, y) = abs(sum(x[i] * y[mod1(i + 1, end)] - x[mod1(i + 1, end)] * y[i] - for i in eachindex(x))) / 2 obj = joinpath(pkgdir(VortexStepMethod), "data", "TUDELFT_V3_KITE", "V3_25.obj") vertices, faces = VortexStepMethod.ObjAdapter.read_faces(obj) canopy = VortexStepMethod.ObjAdapter.perpendicular_sections(vertices, faces, 18; @@ -290,19 +289,25 @@ end @test isnothing(crossing_panels(def.x, def.y)) end end +end - # a clearance-padded wrap is re-wrapped at zero clearance, neither padded again - # nor collapsed by the ball falling between its long panels +@testset "re-wrapping a clearance-padded wrap keeps its area" begin + shoelace(x, y) = abs(sum(x[i] * y[mod1(i + 1, end)] - x[mod1(i + 1, end)] * y[i] + for i in eachindex(x))) / 2 xw, yw = shrink_wrap(read_dat_coordinates(joinpath(@__DIR__, "data", "test_airfoil.dat"))..., ShrinkWrap()) def = deform_section(xw, yw, 0.0; wrap_method=ShrinkWrap()) @test isnothing(crossing_panels(def.x, def.y)) @test shoelace(def.x, def.y) ≈ shoelace(xw, yw) rtol = 0.05 +end - @test_logs (:warn, r"cross") match_mode=:any shrink_wrap(canopy[2].x_airfoil, - canopy[2].y_airfoil, - ShrinkWrap()) +@testset "shrink_wrap warns when its contour crosses itself" begin + x = collect(range(0.0, 1.0, 400)) + y = 0.02 .* sin.(20pi .* x) .+ 0.05 .* sin.(pi .* x) + @test_logs (:warn, r"cross") match_mode=:any shrink_wrap(x, y, + ShrinkWrap(clearance=0.05)) + @test_logs min_level=Logging.Warn shrink_wrap(x, y, ShrinkWrap(clearance=0.0)) end @testset "generate_polar_from_coordinates POLAR_VECTORS sweep" begin From 0d77ab2dab50df4e43e689b61e2af9d319ccb45a Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Mon, 21 Sep 2026 14:21:16 +0200 Subject: [PATCH 6/8] Take #320's XFoil contour guard back out, keeping crossing_panels #320 closed unmerged, its fix having gone to xfoil_light instead. This branch carried it only as a base; the shrink wrap's warning still uses crossing_panels. Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 6 ------ docs/src/private_functions.md | 3 +-- .../airfoil_solvers/xfoil_solver.jl | 21 ------------------- test/airfoil_aero/test_airfoil_aero.jl | 21 +++---------------- 4 files changed, 4 insertions(+), 47 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index fd4963ed..45562af9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -56,12 +56,6 @@ panel's chordwise trailing segment were affected. - With `artificial_damping` on, an iteration whose circulation is already smooth no longer re-applies the previous iteration's damping correction. -- `analyze_sweep(::XFoilSolver, ...)` throws `ArgumentError` for a contour XFoil's - panel code cannot take, rather than handing it to the Fortran. A self-crossing - contour could reach a bare `STOP` there, ending the Julia process with exit code 0 - and no exception; a contour with more nodes than XFoil's panel arrays hold was - refused by XFoil, which then solved every angle on whichever airfoil was loaded - before and returned it as this one. ## VortexStepMethod v5.1.1 2026-09-12 diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 314a8789..8d8cbbaa 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -175,9 +175,8 @@ resample_arc smoothed_curvature ``` -### Section solver preconditions +### Section contour checks ```@docs -validate_xfoil_contour crossing_panels segments_cross side_of_line diff --git a/src/airfoil_aero/airfoil_solvers/xfoil_solver.jl b/src/airfoil_aero/airfoil_solvers/xfoil_solver.jl index a5b9204c..28ff516e 100644 --- a/src/airfoil_aero/airfoil_solvers/xfoil_solver.jl +++ b/src/airfoil_aero/airfoil_solvers/xfoil_solver.jl @@ -28,23 +28,6 @@ runs (`ncrit=9`, `max_iter=100`, incompressible) so the two backends are compara repanel::Bool = false end -""" - validate_xfoil_contour(def::DeformedSection) - -Enforce what XFoil's panel code needs of `def`'s coordinates: no more nodes than its -panel arrays hold, and no two panels crossing. Throws `ArgumentError` on a violation. -""" -function validate_xfoil_contour(def::DeformedSection) - max_nodes = Xfoil.IQX - 5 - length(def.x) <= max_nodes || throw(ArgumentError( - "XFoil holds $max_nodes panel nodes; this contour has $(length(def.x)).")) - crossing = crossing_panels(def.x, def.y) - isnothing(crossing) || throw(ArgumentError( - "XFoil needs a contour that does not cross itself; panels $(crossing[1]) " * - "and $(crossing[2]) of this one do.")) - return nothing -end - """ analyze_sweep(solver::XFoilSolver, def, alpha_range, Re) -> Vector{SectionSolution} @@ -54,12 +37,8 @@ reinit at each side for convergence. Each converged angle reads the surface pres (`Xfoil.cpdump`) and the boundary layer (`Xfoil.bldump`, giving `cf` and the node coordinates) at the same panel nodes. Non-converged angles yield empty node arrays and `NaN` confidence. - -Throws `ArgumentError` for a contour XFoil's panel code cannot take, see -[`validate_xfoil_contour`](@ref). """ function analyze_sweep(solver::XFoilSolver, def::DeformedSection, alpha_range, Re) - validate_xfoil_contour(def) Xfoil.set_coordinates(def.x, def.y) solver.repanel && Xfoil.pane(npan=solver.npan) sols = Vector{SectionSolution}(undef, length(alpha_range)) diff --git a/test/airfoil_aero/test_airfoil_aero.jl b/test/airfoil_aero/test_airfoil_aero.jl index b2f0286d..15aa3491 100644 --- a/test/airfoil_aero/test_airfoil_aero.jl +++ b/test/airfoil_aero/test_airfoil_aero.jl @@ -5,9 +5,7 @@ import VortexStepMethod using VortexStepMethod.AirfoilAero: KulfanParameters, LeastSquaresFit, ShrinkWrap, shrink_wrap, fit_kulfan_parameters, kulfan_to_coordinates, neuralfoil_aero, class_function, bernstein_basis, - leading_edge_basis, normalize_airfoil, crossing_panels, - validate_xfoil_contour, DeformedSection, XFoilSolver, - analyze_sweep, Xfoil + leading_edge_basis, normalize_airfoil, crossing_panels using VortexStepMethod: SectionAero, section_surface, read_section_aero using VortexStepMethod.AirfoilAero: write_section_aero @@ -362,29 +360,16 @@ end collect(extrema(fitted_y))) < 1e-4 end -@testset "a contour XFoil has no solution for is refused before XFoil sees it" begin +@testset "crossing_panels finds the first pair of crossing panels" begin clean = KulfanParameters(fill(0.15, 8), fill(-0.15, 8), 0.0, 0.0) x, y = collect.(kulfan_to_coordinates(clean; n_points=60)) - alphas = deg2rad.([0.0]) - @test isnothing(validate_xfoil_contour(DeformedSection(clean, x, y))) + @test isnothing(crossing_panels(x, y)) # the upper surface driven through the lower one over a stretch of the chord folded = copy(y) folded[20:40] .= -3 .* folded[20:40] @test !isnothing(crossing_panels(x, folded)) - @test_throws ArgumentError analyze_sweep(XFoilSolver(), - DeformedSection(clean, x, folded), alphas, 1e6) - - gapped = KulfanParameters(fill(0.15, 8), fill(-0.15, 8), 0.0, 0.02) - gapped_x, gapped_y = collect.(kulfan_to_coordinates(gapped; n_points=60)) - @test (gapped_x[1], gapped_y[1]) != (gapped_x[end], gapped_y[end]) - @test isnothing(validate_xfoil_contour(DeformedSection(gapped, gapped_x, gapped_y))) # a contour whose only crossing is the panel closing it back to node 1 @test crossing_panels([0.0, 1.0, 1.0, 2.0], [0.0, 2.0, -2.0, 1.0]) == (2, 4) - - crowded_x, crowded_y = collect.(kulfan_to_coordinates(clean; n_points=200)) - @test length(crowded_x) > Xfoil.IQX - 5 - @test_throws ArgumentError analyze_sweep(XFoilSolver(), - DeformedSection(clean, crowded_x, crowded_y), alphas, 1e6) end From 54ed7b518df68e7999341f381032165aafc432db Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Mon, 21 Sep 2026 14:29:41 +0200 Subject: [PATCH 7/8] Count collinear panels that do not overlap as not crossing Along a straight edge the side-of-line products are rounding noise of either sign, so two disjoint panels on it could be reported as crossing. That made the re-wrap test fail on CI's dependency versions near a V3 section's straight nose. segments_cross now also requires the two segments' extents to overlap. Co-Authored-By: Claude Opus 5 --- docs/src/private_functions.md | 1 + src/airfoil_aero/airfoil_solvers/common.jl | 11 ++++++++++- test/airfoil_aero/test_airfoil_aero.jl | 7 +++++++ 3 files changed, 18 insertions(+), 1 deletion(-) diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 6dc763da..bf78b7ec 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -179,6 +179,7 @@ smoothed_curvature ```@docs crossing_panels segments_cross +extents_overlap side_of_line ``` diff --git a/src/airfoil_aero/airfoil_solvers/common.jl b/src/airfoil_aero/airfoil_solvers/common.jl index b53654c3..b887a5a7 100644 --- a/src/airfoil_aero/airfoil_solvers/common.jl +++ b/src/airfoil_aero/airfoil_solvers/common.jl @@ -130,12 +130,21 @@ side_of_line(a, b, p) = (b[1] - a[1]) * (p[2] - a[2]) - (b[2] - a[2]) * (p[1] - Whether the segments `p`-`q` and `r`-`s` cross properly, each strictly separating the other's endpoints. Touching at an endpoint or lying along each other does not -count. +count, nor do segments whose coordinate extents do not overlap. """ segments_cross(p, q, r, s) = + extents_overlap(p[1], q[1], r[1], s[1]) && extents_overlap(p[2], q[2], r[2], s[2]) && side_of_line(p, q, r) * side_of_line(p, q, s) < 0 && side_of_line(r, s, p) * side_of_line(r, s, q) < 0 +""" + extents_overlap(a1, a2, b1, b2) -> Bool + +Whether the intervals spanned by `a1`, `a2` and by `b1`, `b2` overlap. +""" +extents_overlap(a1, a2, b1, b2) = + max(min(a1, a2), min(b1, b2)) <= min(max(a1, a2), max(b1, b2)) + """ crossing_panels(x, y) -> Tuple{Int,Int} or nothing diff --git a/test/airfoil_aero/test_airfoil_aero.jl b/test/airfoil_aero/test_airfoil_aero.jl index 15aa3491..7a1e7df2 100644 --- a/test/airfoil_aero/test_airfoil_aero.jl +++ b/test/airfoil_aero/test_airfoil_aero.jl @@ -372,4 +372,11 @@ end # a contour whose only crossing is the panel closing it back to node 1 @test crossing_panels([0.0, 1.0, 1.0, 2.0], [0.0, 2.0, -2.0, 1.0]) == (2, 4) + + # disjoint panels along one straight edge, whose side tests round to either sign + edge_x = [0.0008647734705084547, 0.000648580102881341, 0.00043238673525422734, + 0.00021619336762711367, 0.0] + edge_y = [0.0044110473582233455, 0.003308285518667509, 0.0022055236791116727, + 0.0011027618395558364, 0.0] + @test isnothing(crossing_panels([edge_x; 0.001], [edge_y; -0.002])) end From ee7a6b3c4ceb4454a60b33f5e2c505b6c15bfb04 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Mon, 21 Sep 2026 22:35:16 +0200 Subject: [PATCH 8/8] Cut the loops the clearance offset makes, so the default wrap is simple On a thin canopy the clearance offset of the contact polygon turns back on itself where contacts lie closer than the clearance, and shrink_wrap then returned a crossing contour and warned: 6 of 18 V3 sections at the default MeshSettings. pivot_contour now walks the offset and replaces each loop by its crossing point (cut_loops). Closes #363 Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 3 +++ docs/src/private_functions.md | 1 + src/airfoil_aero/shrink_wrap.jl | 37 ++++++++++++++++++++++++-- test/airfoil_aero/test_airfoil_aero.jl | 27 ++++++++++++++++--- 4 files changed, 62 insertions(+), 6 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5baef239..4b4b5039 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -70,6 +70,9 @@ `generate_aero_matrices` and `generate_polar_from_coordinates`), at zero clearance. `shrink_wrap` warns when the contour it returns crosses itself. A thin wrapped section came back with its surfaces crossing. +- `shrink_wrap` cuts out the loops its `clearance` offset makes on a thin canopy, so a + V3 section wrapped at the default `MeshSettings` is a simple closed curve and no longer + warns; 6 of 18 crossed themselves. - `get_lower_upper`, and with it the flap hinge in `deform_section`, takes the lower and upper surface heights where the contour crosses `x = crease_frac`. It took the nearest points below and above `y = 0`, which on a cambered section put the hinge near the diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 59c21f94..c2abbaf4 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -171,6 +171,7 @@ pivot_step push_arc! edge_normal pivot_contour +cut_loops largest_linking_gap densify_contour enforce_min_spacing! diff --git a/src/airfoil_aero/shrink_wrap.jl b/src/airfoil_aero/shrink_wrap.jl index 31493137..26282e7d 100644 --- a/src/airfoil_aero/shrink_wrap.jl +++ b/src/airfoil_aero/shrink_wrap.jl @@ -145,7 +145,8 @@ of the cloud from its leftmost point; the points it touches are the polygon's vertices in order, the ones inside a concavity narrower than the disk having been skipped and bridged straight. The offset rounds each convex vertex with an arc of radius `clearance` and chamfers each reflex one across the bisector, so it holds that -distance from every contact. +distance from every contact, and the loops the offset makes where contacts lie closer +than `clearance` are cut out ([`cut_loops`](@ref)). """ function pivot_contour(x, y, r, clearance) start = argmin(x) @@ -182,7 +183,39 @@ function pivot_contour(x, y, r, clearance) end from = to end - return px, py + return cut_loops(px, py) +end + +""" + cut_loops(px, py) -> (x, y) + +The closed polyline `(px, py)` walked from its first node with every loop it makes +by crossing itself cut out, the crossing point taking the loop's place, so what comes +back is a simple closed curve. +""" +function cut_loops(px, py) + n = length(px) + x, y = [px[1]], [py[1]] + for k in 2:n+1 + closing = k > n + p = (px[mod1(k, n)], py[mod1(k, n)]) + a = (x[end], y[end]) + for j in (closing ? 2 : 1):length(x)-2 + r, s = (x[j], y[j]), (x[j+1], y[j+1]) + segments_cross(r, s, a, p) || continue + t = side_of_line(r, s, a) / (side_of_line(r, s, a) - side_of_line(r, s, p)) + resize!(x, j) + resize!(y, j) + push!(x, a[1] + t * (p[1] - a[1])) + push!(y, a[2] + t * (p[2] - a[2])) + break + end + if !closing + push!(x, p[1]) + push!(y, p[2]) + end + end + return x, y end """ diff --git a/test/airfoil_aero/test_airfoil_aero.jl b/test/airfoil_aero/test_airfoil_aero.jl index e55a65cf..fd0665b0 100644 --- a/test/airfoil_aero/test_airfoil_aero.jl +++ b/test/airfoil_aero/test_airfoil_aero.jl @@ -301,12 +301,31 @@ end @test shoelace(def.x, def.y) ≈ shoelace(xw, yw) rtol = 0.05 end -@testset "shrink_wrap warns when its contour crosses itself" begin +@testset "shrink_wrap cuts the loops its clearance offset makes" begin + obj = joinpath(pkgdir(VortexStepMethod), "data", "TUDELFT_V3_KITE", "V3_25.obj") + vertices, faces = VortexStepMethod.ObjAdapter.read_faces(obj) + canopy = VortexStepMethod.ObjAdapter.perpendicular_sections(vertices, faces, 18; + n_bins=100) + for section in canopy + xw, yw = @test_logs min_level=Logging.Warn shrink_wrap(section.x_airfoil, + section.y_airfoil, + ShrinkWrap()) + @test isnothing(crossing_panels(xw, yw)) + end x = collect(range(0.0, 1.0, 400)) y = 0.02 .* sin.(20pi .* x) .+ 0.05 .* sin.(pi .* x) - @test_logs (:warn, r"cross") match_mode=:any shrink_wrap(x, y, - ShrinkWrap(clearance=0.05)) - @test_logs min_level=Logging.Warn shrink_wrap(x, y, ShrinkWrap(clearance=0.0)) + @test_logs min_level=Logging.Warn shrink_wrap(x, y, ShrinkWrap(clearance=0.05)) +end + +@testset "shrink_wrap warns when its contour crosses itself" begin + x = collect(range(1.0, 0.0, 800)) + camber = 0.05 .* sin.(pi .* x) .+ 0.02 .* sin.(20pi .* x) + loop(half) = (vcat(x, reverse(x)[2:end]), + vcat(camber .+ half, reverse(camber .- half)[2:end])) + @test_logs (:warn, r"cross") match_mode=:any shrink_wrap(loop(1e-4)..., + ShrinkWrap(clearance=0.0)) + @test_logs min_level=Logging.Warn shrink_wrap(loop(1e-3)..., + ShrinkWrap(clearance=0.0)) end @testset "generate_polar_from_coordinates POLAR_VECTORS sweep" begin