Re-wrap a deformed section with the ball it was first wrapped with, and cut the loops out of the clearance offset

CHANGELOG.md

@@ -67,6 +67,9 @@
 - `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.
 - BREAKING: `obj_to_yaml` and `perpendicular_sections` spread the sections evenly over
   the span, measured along the quarter-chord line without its chordwise component,
   instead of over leading-edge arc length, and `wingtip_distance` is that spanwise
@@ -77,6 +80,14 @@
 
 ### 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` 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.
 - On a body whose wings span different directions, such as a wing and a vertical fin,
   `solve!`, `solve` and `linearize` take each panel's lift, drag and side directions
   from its own wing's `spanwise_direction`, not the first wing's. `solve` computes

docs/src/private_functions.md

@@ -177,12 +177,22 @@ pivot_step
 push_arc!
 edge_normal
 pivot_contour
+cut_loops
 largest_linking_gap
+densify_contour
 enforce_min_spacing!
 resample_arc
 smoothed_curvature
 ```
 
+### Section contour checks
+```@docs
+crossing_panels
+segments_cross
+extents_overlap
+side_of_line
+```
+
 ### NeuralFoil network
 ```@docs
 load_neuralfoil_model
@@ -256,7 +266,6 @@ build_section
 contour_to_airfoil
 plane_contour_to_airfoil
 reorder_airfoil_selig
-densify_contour
 create_interpolations
 find_circle_center_and_radius
 march_edges

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),

src/airfoil_aero/airfoil_io.jl

@@ -29,29 +29,31 @@ 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 table, CSV or
 Arrow as the suffix of `output_path` says. 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 [`XFoilSolver`](@ref) to pick the backend. With
-`delta_range === nothing` the sweep is over `alpha_range` only and written as a
-`POLAR_VECTORS` table (returns the `Vector{SectionSolution}`); pass a `delta_range` of
-trailing-edge deflections to sweep `(alpha, delta)` and write a long-format
-`POLAR_MATRICES` table (returns the `(cl, cd, cm)` matrices). Both angle ranges are in
-degrees. `crease_frac` is the chordwise hinge location (0–1) about which each
-`delta_range` deflection pivots. With `dat_prefix` set, each deflected shape is also
-written to `{dat_prefix}_{delta_suffix(δ)}.dat`.
+Wrap a raw or open single-membrane slice with [`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` table (returns the
+`Vector{SectionSolution}`); pass a `delta_range` of trailing-edge deflections to sweep
+`(alpha, delta)` and write a long-format `POLAR_MATRICES` table (returns the `(cl, cd,
+cm)` matrices). Both angle ranges are in degrees. `crease_frac` is the chordwise hinge
+location (0–1) about which each `delta_range` deflection pivots. With `dat_prefix` set,
+each deflected shape is also written to `{dat_prefix}_{delta_suffix(δ)}.dat`.
 """
 function generate_polar_from_coordinates(x::Vector, y::Vector, output_path::String;
                                          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(output_path, sols)
         return sols
@@ -61,7 +63,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(output_path, alphas, deltas, cl, cd, cm)
     return (cl, cd, cm)
 end

src/airfoil_aero/airfoil_solvers/common.jl

@@ -81,8 +81,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)
@@ -92,26 +91,73 @@ 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(; 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)
 end
 
+"""
+    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.
+"""
+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, 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
+
+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
+        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}
 

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
@@ -29,7 +30,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"))
@@ -46,7 +47,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)

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)

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)