diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3b6605ea..daaff4f7 100644
--- a/CHANGELOG.md
+++ b/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
diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
index 3d8de213..19be30db 100644
--- a/docs/src/private_functions.md
+++ b/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
diff --git a/ext/VortexStepMethodMakieExt.jl b/ext/VortexStepMethodMakieExt.jl
index ebde9043..87115e04 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 91d0c6ac..b0ee6e2a 100644
--- a/src/airfoil_aero/airfoil_io.jl
+++ b/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
diff --git a/src/airfoil_aero/airfoil_solvers/common.jl b/src/airfoil_aero/airfoil_solvers/common.jl
index dc047053..3a99ed13 100644
--- a/src/airfoil_aero/airfoil_solvers/common.jl
+++ b/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}
 
diff --git a/src/airfoil_aero/geometry_gen.jl b/src/airfoil_aero/geometry_gen.jl
index c682de15..d9ccbfe6 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
@@ -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)
diff --git a/src/airfoil_aero/polar_export.jl b/src/airfoil_aero/polar_export.jl
index 9ee408e4..9e8c46d2 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..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
 
 """
@@ -218,6 +251,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 +383,24 @@ 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,
+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
 `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 +422,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 fd3f5b43..35abf5b6 100644
--- a/src/obj_adapter/obj_slice.jl
+++ b/src/obj_adapter/obj_slice.jl
@@ -344,26 +344,6 @@ function march_edges(vertices, faces; step)
             point=[r.point for r in rows], tangent=[r.tangent for r in rows])
 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
 
@@ -379,7 +359,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 8c8e17df..609e84b6 100644
--- a/src/obj_adapter/obj_to_yaml.jl
+++ b/src/obj_adapter/obj_to_yaml.jl
@@ -238,7 +238,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 520dfe19..da46cbe3 100644
--- a/src/surfplan_adapter/SurfplanAdapter.jl
+++ b/src/surfplan_adapter/SurfplanAdapter.jl
@@ -76,7 +76,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 e58eb583..a58befc0 100644
--- a/test/airfoil_aero/test_airfoil_aero.jl
+++ b/test/airfoil_aero/test_airfoil_aero.jl
@@ -1,10 +1,11 @@
 using Test
+using Logging
 using VortexStepMethod.AirfoilAero
 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
 using VortexStepMethod: SectionAero, section_surface, read_section_aero
 using VortexStepMethod.AirfoilAero: write_section_aero, write_aero_matrix
 using VortexStepMethod: load_polar_data
@@ -274,6 +275,59 @@ end
     end
 end
 
+@testset "re-wrapping a wrapped section keeps it a simple closed curve" 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 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
+end
+
+@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
+
+@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 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
     x, y = read_dat_coordinates(joinpath(@__DIR__, "data", "test_airfoil.dat"))
     csv = joinpath(mktempdir(), "polar.csv")
@@ -424,3 +478,24 @@ end
     @test maximum(abs, collect(extrema(written_y)) .-
                        collect(extrema(fitted_y))) < 1e-4
 end
+
+@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))
+    @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))
+
+    # 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