diff --git a/CHANGELOG.md b/CHANGELOG.md index 3b6605ea..14df8d25 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 3d8de213..97a2bee5 100644 --- a/docs/src/private_functions.md +++ b/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 diff --git a/src/obj_adapter/obj_slice.jl b/src/obj_adapter/obj_slice.jl index fd3f5b43..2472c3a3 100644 --- a/src/obj_adapter/obj_slice.jl +++ b/src/obj_adapter/obj_slice.jl @@ -269,16 +269,82 @@ 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) @@ -286,62 +352,17 @@ function march_edges(vertices, faces; step) 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 """ diff --git a/test/obj_adapter/test_obj_adapter.jl b/test/obj_adapter/test_obj_adapter.jl index 19a2b147..468931d1 100644 --- a/test/obj_adapter/test_obj_adapter.jl +++ b/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]