Merge pull request #376 from OpenSourceAWE/agent/361-march-edges-gives-each-station-the-previ

march_edges takes each station's tangent from its neighbours, so no cut near a closing tip grazes the surface

Author: 1-Bart-1

CHANGELOG.md

@@ -82,6 +82,10 @@
   from its own wing's `spanwise_direction`, not the first wing's. `solve` computes
   `wing_span` and `aspect_ratio_projected` from the extent of all wings along the
   first wing's span, through the new `calculate_span(wings, spanwise_direction)`.
+- `perpendicular_sections` and `obj_to_yaml` orient each cut by the leading edge through
+  the neighbouring stations, not by the step that reached the station, which near a
+  closing tip grazed the surface and returned an 80%-thick section. Every section moves
+  slightly; regenerate existing geometry to benefit.
 - `ELLIPTIC` initial circulation works on a body with more than one wing, where it threw
   an `ArgumentError`: each wing gets an ellipse over its own span, along its own
   `spanwise_direction` and centred on its own mid-span, also for a single wing off y = 0.

docs/src/private_functions.md

@@ -260,6 +260,9 @@ densify_contour
 create_interpolations
 find_circle_center_and_radius
 march_edges
+cut_station
+tip_station
+march_stations
 airfoils_from_yaml
 write_geometry_yaml
 resolve_aero_geometry

src/obj_adapter/obj_slice.jl

@@ -269,79 +269,100 @@ function slice_mesh_at_plane(vertices, faces, point, normal; tol=1e-6)
     return segments
 end
 
+"""
+    cut_station(vertices, faces, point, tangent) -> (; le, te) or nothing
+
+Cut the mesh with the vertical spanwise plane through `point` whose normal is the
+spanwise component of `tangent`, and return the crossings with the least and greatest
+`x` as the leading and trailing edge; `nothing` where the plane misses the mesh.
+"""
+function cut_station(vertices, faces, point, tangent)
+    segments = slice_mesh_at_plane(vertices, faces, point,
+                                   normalize([0.0, tangent[2], 0.0]))
+    isempty(segments) && return nothing
+    crossings = Vector{Float64}[]
+    for (start_point, end_point) in segments
+        push!(crossings, start_point)
+        push!(crossings, end_point)
+    end
+    return (; le=argmin(crossing -> crossing[1], crossings),
+            te=argmax(crossing -> crossing[1], crossings))
+end
+
+"""
+    tip_station(vertices, faces, prev_le, tangent, step) -> (; le, te, point) or nothing
+
+Bisect the step of length `step` from `prev_le` along `tangent` for the outermost cut
+that still yields an airfoil section whose leading edge stays within `step` of
+`prev_le` chordwise; `nothing` where no such cut exists.
+"""
+function tip_station(vertices, faces, prev_le, tangent, step)
+    lo, hi, tip = 0.0, step, nothing
+    for _ in 1:24
+        mid = (lo + hi) / 2
+        probe = prev_le .+ mid .* tangent
+        here = cut_station(vertices, faces, probe, tangent)
+        valid = here !== nothing &&
+                abs(here.le[1] - prev_le[1]) < step &&
+                build_section(vertices, faces, here.le, here.te, probe, tangent) !== nothing
+        valid ? (lo = mid; tip = (; here.le, here.te, point=probe)) : (hi = mid)
+    end
+    return tip
+end
+
+"""
+    march_stations(vertices, faces, start_le, step, direction) -> Vector{NamedTuple}
+
+March the leading edge from `start_le` in steps of arc length `step` towards `+y` for
+`direction = 1.0` or `-y` for `-1.0`, returning one `(; le, te, point)` per station,
+the tip station last. Stops when the leading edge stops advancing spanwise.
+"""
+function march_stations(vertices, faces, start_le, step, direction)
+    rows = NamedTuple[]
+    prev_le = start_le
+    tangent = [0.0, direction, 0.0]
+    for _ in 1:10_000
+        probe = prev_le .+ step .* tangent
+        found = cut_station(vertices, faces, probe, tangent)
+        if found === nothing
+            tip = tip_station(vertices, faces, prev_le, tangent, step)
+            tip === nothing || push!(rows, tip)
+            break
+        end
+        le_step = found.le .- prev_le
+        (norm(le_step) < 1e-9 || le_step[2] * direction <= 0.0) && break
+        push!(rows, (; found.le, found.te, point=probe))
+        tangent = normalize(le_step)
+        prev_le = found.le
+    end
+    return rows
+end
+
 """
     march_edges(vertices, faces; step) -> (; le, te, point, tangent)
 
-March the leading edge outward from mid-span in both directions in steps of arc
-length `step`. Each cut is a vertical spanwise plane (both the chordwise and vertical
-components of the running LE tangent dropped from the normal) so a tip that curls
-downward can't tilt the plane toward horizontal, where its min-chord "LE" pick would
-jump across the wing. Marching stops when the leading edge stops advancing spanwise.
-Cuts sample mesh *edges*, so the picks are robust to vertex density. Returns, ordered
-along the span, the LE/TE points and each cut's plane origin and tangent. Build the
+Cut the mesh at mid-span and march the leading edge outward from there in both
+directions with `march_stations`. Returns, ordered along the span, the LE/TE points,
+each cut's plane origin, and the LE tangent: the central difference over the
+neighbouring stations, where the end stations reuse their inner neighbour's. Build the
 airfoil for a chosen station with `build_section`.
 """
 function march_edges(vertices, faces; step)
     ys = [v[2] for v in vertices]
     y_min, y_max = extrema(ys)
     y_mid = (y_min + y_max) / 2
 
-    # Vertical spanwise plane (z dropped): a tip curl mustn't tilt the cut horizontal.
-    cut(point, tangent) = begin
-        segments = slice_mesh_at_plane(vertices, faces, point,
-                                       normalize([0.0, tangent[2], 0.0]))
-        isempty(segments) && return nothing
-        crossings = Vector{Float64}[]
-        for (start_point, end_point) in segments
-            push!(crossings, start_point)
-            push!(crossings, end_point)
-        end
-        return (; le=argmin(pt -> pt[1], crossings), te=argmax(pt -> pt[1], crossings))
-    end
-
-    mid_cut = cut([0.0, y_mid, 0.0], [0.0, 1.0, 0.0])
+    mid_cut = cut_station(vertices, faces, [0.0, y_mid, 0.0], [0.0, 1.0, 0.0])
     mid_cut === nothing && error("Could not slice mid-span; check mesh / rotation.")
 
-    march(direction) = begin
-        rows = NamedTuple[]
-        prev_le = mid_cut.le
-        tangent = [0.0, direction, 0.0]
-        for _ in 1:10_000
-            probe = prev_le .+ step .* tangent
-            found = cut(probe, tangent)
-            if found === nothing
-                # Bisect the last step to land the tip station on the outermost cut
-                # that still yields a valid (non-degenerate) airfoil section.
-                lo, hi, tip = 0.0, step, nothing
-                for _ in 1:24
-                    mid = (lo + hi) / 2
-                    probe_mid = prev_le .+ mid .* tangent
-                    here = cut(probe_mid, tangent)
-                    # Reject a near-degenerate tip slice whose min-chord "LE" has jumped
-                    # chordwise (an artifact that would misplace the tip station).
-                    valid = here !== nothing &&
-                            abs(here.le[1] - prev_le[1]) < step &&
-                            build_section(vertices, faces, here.le, here.te,
-                                          probe_mid, tangent) !== nothing
-                    valid ? (lo = mid; tip = (; here.le, here.te, point=probe_mid)) :
-                            (hi = mid)
-                end
-                tip === nothing || push!(rows, (; tip.le, tip.te, tip.point, tangent))
-                break
-            end
-            le_step = found.le .- prev_le
-            (norm(le_step) < 1e-9 || le_step[2] * direction <= 0.0) && break
-            push!(rows, (; found.le, found.te, point=probe, tangent))
-            tangent = normalize(le_step)
-            prev_le = found.le
-        end
-        return rows
-    end
-
-    center = (; mid_cut.le, mid_cut.te, point=[0.0, y_mid, 0.0], tangent=[0.0, 1.0, 0.0])
-    rows = vcat(reverse(march(-1.0)), [center], march(1.0))
-    return (; le=[r.le for r in rows], te=[r.te for r in rows],
-            point=[r.point for r in rows], tangent=[r.tangent for r in rows])
+    center = (; mid_cut.le, mid_cut.te, point=[0.0, y_mid, 0.0])
+    negative_y = march_stations(vertices, faces, mid_cut.le, step, -1.0)
+    positive_y = march_stations(vertices, faces, mid_cut.le, step, 1.0)
+    rows = vcat(reverse(negative_y), [center], positive_y)
+    le = [row.le for row in rows]
+    middle = clamp.(eachindex(le), 2, length(le) - 1)
+    tangent = [normalize(le[j + 1] .- le[j - 1]) for j in middle]
+    return (; le, te=[row.te for row in rows], point=[row.point for row in rows], tangent)
 end
 
 """

test/obj_adapter/test_obj_adapter.jl

@@ -29,6 +29,22 @@ obj_path = normpath(joinpath(@__DIR__, "..", "..",
         end
     end
 
+    @testset "every marched station cuts across the surface, closing tips included" begin
+        vertices, faces = read_faces(obj_path)
+        ys = [v[2] for v in vertices]
+        for n_bins in (60, 90)
+            march = ObjAdapter.march_edges(vertices, faces;
+                                           step=(maximum(ys) - minimum(ys)) / n_bins)
+            for i in eachindex(march.le)
+                section = ObjAdapter.build_section(vertices, faces, march.le[i],
+                    march.te[i], march.point[i], march.tangent[i])
+                depth = section === nothing ? Inf :
+                        maximum(section.y_airfoil) - minimum(section.y_airfoil)
+                @test depth < 0.3
+            end
+        end
+    end
+
     @testset "station_indices spreads sections evenly in span past a raked tip" begin
         # Straight LE over |y| ≤ 1, then tip caps whose LE runs 40 mm aft per 3 mm span.
         cap = [(1.0 + 0.003k, 0.04k) for k in 1:20]