From 4a926cb7f54267567e9c02aac868cfe610dc5480 Mon Sep 17 00:00:00 2001 From: Bart Date: Fri, 14 Aug 2026 14:55:04 +0200 Subject: [PATCH 1/8] Keep refined surface tables through a coarsening remesh compute_refined_section_interpolation! reblended every refined section's SectionAero from the unrefined sections unconditionally, while aero_data was preserved under use_prior_polar. A wing rebuilt onto fewer structural stations therefore kept full-resolution polars and lost the contour/Cp/cf tables that pressure integration reads: on the TU Delft V3, 37 generated tables collapsed to the 10 sampled onto its struts, and both the pressure loads and the lofted airfoil plot were interpolations of those 10. The reblend is now skipped whenever the polars are being preserved and the refined sections already carry tables, so it still fills them the first time round. obj_to_yaml gains geometry_path, naming the geometry YAML instead of always writing output_dir/geometry.yaml. Emitted table references carry the path from the YAML's directory to output_dir, which is what the loader resolves them against, so a generated dataset can keep its bulk in a subdirectory while the geometry sits with the hand-written ones. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 13 ++++++ src/obj_adapter/obj_to_yaml.jl | 57 +++++++++++++++++++++--- src/wing_geometry.jl | 35 ++++++++++++--- test/obj_adapter/test_obj_adapter.jl | 23 ++++++++++ test/wing_geometry/test_wing_geometry.jl | 54 ++++++++++++++++++++++ 5 files changed, 170 insertions(+), 12 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2e73b25a..8ad3593c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,12 @@ and `surfplan_to_aero_yaml`: `:csv` (default, readable) or `:arrow` (binary, ~40× faster to load and 2.6× smaller). `read_section_aero` detects the format from the file suffix, so a geometry YAML can reference either. +- `geometry_path` keyword on `obj_to_yaml`, naming the geometry YAML itself + instead of always writing `output_dir/geometry.yaml`. Point it outside the + table directory and the emitted table references carry the path from the YAML's + directory to `output_dir`, which is what the geometry loader resolves them + against — so a generated dataset can keep its bulk in a subdirectory while the + geometry sits with the hand-written ones. - `convert_node_table` and `write_node_rows` rewrite a per-node table in the format the destination suffix names. `obj_to_yaml` migrates an existing dataset with them when `table_format` differs from what the directory holds, so a dataset @@ -18,6 +24,13 @@ instead of its own (`ZEROS`, `0.05`), and never said what `core_radius_fraction` measures. It now documents the `Solver` defaults and cites Damiani et al. (2019) for the 0.05 cut-off. +- A remesh under `use_prior_polar` no longer resamples the refined sections' + `SectionAero` surface tables down to whatever unrefined sections survive it. + `compute_refined_section_interpolation!` reblended contour, `cp` and `cf` from the + unrefined sections unconditionally while `aero_data` was preserved, so a wing rebuilt + onto fewer structural stations kept full-resolution polars but lost the surface tables + pressure integration reads. The reblend is now skipped when the polars are preserved + and the refined sections already carry tables. ### Changed - `SolverSettings` now defaults to the same values as `Solver`: `core_radius_fraction` diff --git a/src/obj_adapter/obj_to_yaml.jl b/src/obj_adapter/obj_to_yaml.jl index 78425319..27496fcd 100644 --- a/src/obj_adapter/obj_to_yaml.jl +++ b/src/obj_adapter/obj_to_yaml.jl @@ -54,8 +54,13 @@ to use XFoil instead. Each section's polar is written as `POLAR_VECTORS`. `crease_frac` is the chordwise hinge location (0–1) about which each `delta_range` trailing-edge deflection pivots. -With `force=false` (default) an existing `geometry.yaml` in `output_dir` is reused; -`force=true` regenerates it (e.g. after changing `delta_range` or the mesh). +With `force=false` (default) an existing geometry YAML is reused; `force=true` +regenerates it (e.g. after changing `delta_range` or the mesh). + +`geometry_path` names the YAML itself, `output_dir/geometry.yaml` by default. Point it +elsewhere to keep the YAML out of the table directory — the emitted table references +then carry the path from the YAML's directory to `output_dir`, which is what the +geometry loader resolves them against. `wrap_method` ([`ShrinkWrap`](@ref)) wraps each slice's point cloud into a clean closed airfoil, robust both to the noisy interior-structure points (ribs, spars) of a ram-air @@ -99,6 +104,9 @@ Rewrite a generated dataset's per-node `Cp`/`cf` tables in `table_format` and po an existing directory change format without re-running the airfoil solver that produced it — the polars are the slow part and they are untouched. The source tables are left in place. + +`output_dir` is what the YAML's relative table references resolve against, which is +its own directory — the rule the geometry loader follows. """ function migrate_node_tables(yaml_path::String, output_dir::String, table_format::Symbol; verbose::Bool=true) @@ -138,6 +146,41 @@ function migrate_node_tables(yaml_path::String, output_dir::String, return yaml_path end +""" + table_path_prefix(geometry_path, output_dir) -> String + +Path from the geometry YAML's directory to the table directory, empty when they are +the same. Table references resolve against the YAML's own directory, so a YAML +written outside `output_dir` has to carry this hop. +""" +function table_path_prefix(geometry_path::String, output_dir::String) + yaml_dir = dirname(abspath(geometry_path)) + tables = abspath(output_dir) + yaml_dir == tables && return "" + return relpath(tables, yaml_dir) +end + +""" + prefix_table_paths!(airfoil_rows, prefix) -> airfoil_rows + +Prepend `prefix` to every relative table reference in the `info_dict` of each row, so +a geometry YAML written outside the table directory still resolves them. A no-op on +an empty prefix. +""" +function prefix_table_paths!(airfoil_rows, prefix::String) + isempty(prefix) && return airfoil_rows + for row in airfoil_rows + info = row[end] + info isa AbstractDict || continue + for (key, value) in info + (endswith(String(key), "_file") || key == "csv_file_path") || continue + value isa AbstractString && !isabspath(value) && + (info[key] = joinpath(prefix, value)) + end + end + return airfoil_rows +end + function obj_to_yaml(obj_path::String, output_dir::String; n_sections::Int, Re::Real, alpha_range=-180:1:180, delta_range=nothing, @@ -146,16 +189,19 @@ function obj_to_yaml(obj_path::String, output_dir::String; reuse_valid_airfoils::Bool=true, max_thickness_ratio::Real=2.0, spanwise_direction=[0.0, 1.0, 0.0], rotation=I, wingtip_distance=0.05, crease_frac=0.75, force::Bool=false, - verbose::Bool=true, table_format::Symbol=:csv) + verbose::Bool=true, table_format::Symbol=:csv, + geometry_path::String=joinpath(output_dir, "geometry.yaml")) (!endswith(obj_path, ".obj")) && (obj_path *= ".obj") isfile(obj_path) || error("OBJ file not found: $obj_path") !isapprox(spanwise_direction, [0.0, 1.0, 0.0]) && throw(ArgumentError("Spanwise direction has to be [0.0, 1.0, 0.0]")) - yaml_path = joinpath(output_dir, "geometry.yaml") + yaml_path = geometry_path + mkpath(dirname(abspath(yaml_path))) if !force && isfile(yaml_path) verbose && @info "Reusing existing geometry (force=true to regenerate)" yaml_path - migrate_node_tables(yaml_path, output_dir, table_format; verbose) + migrate_node_tables(yaml_path, dirname(abspath(yaml_path)), table_format; + verbose) return yaml_path end @@ -188,6 +234,7 @@ function obj_to_yaml(obj_path::String, output_dir::String; delta_range, aero_solver, reuse_valid_airfoils, crease_frac, verbose, table_format) isempty(ok) && error("No section produced a valid polar in $obj_path") + prefix_table_paths!(airfoil_rows, table_path_prefix(yaml_path, output_dir)) # Each section uses its nearest airfoil that actually produced a polar — covering # both too-thick (degenerate) fits and sections the solver could not converge. diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index 24a3556f..25c0251b 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -812,6 +812,18 @@ end return all(_has_initialized_section_aero_data, wing.refined_sections) end +""" + _can_reuse_prior_refined_surface_tables(wing) -> Bool + +Whether every refined section already carries a [`SectionAero`](@ref), so a +remesh can keep them instead of reblending from the unrefined sections. False on +a wing that has none, where the blend is what fills them in the first place. +""" +@inline function _can_reuse_prior_refined_surface_tables(wing::AbstractWing) + isempty(wing.refined_sections) && return false + return all(s -> !isnothing(s.section_aero), wing.refined_sections) +end + """ copy_sections_to_refined!(wing; reuse_aero_data=false) @@ -924,7 +936,7 @@ function refine!(wing::AbstractWing{T}; recompute_mapping=true, sort_sections=tr copy_sections_to_refined!(wing; reuse_aero_data) if recompute_mapping compute_refined_panel_mapping!(wing) - compute_refined_section_interpolation!(wing) + compute_refined_section_interpolation!(wing; reuse_aero_data) end update_non_deformed_sections!(wing) return nothing @@ -939,7 +951,7 @@ function refine!(wing::AbstractWing{T}; recompute_mapping=true, sort_sections=tr copy_sections_to_refined!(wing; reuse_aero_data) if recompute_mapping compute_refined_panel_mapping!(wing) - compute_refined_section_interpolation!(wing) + compute_refined_section_interpolation!(wing; reuse_aero_data) end update_non_deformed_sections!(wing) return nothing @@ -959,7 +971,7 @@ function refine!(wing::AbstractWing{T}; recompute_mapping=true, sort_sections=tr reuse_aero_data ? nothing : s2.aero_data) if recompute_mapping compute_refined_panel_mapping!(wing) - compute_refined_section_interpolation!(wing) + compute_refined_section_interpolation!(wing; reuse_aero_data) end update_non_deformed_sections!(wing) return nothing @@ -983,7 +995,7 @@ function refine!(wing::AbstractWing{T}; recompute_mapping=true, sort_sections=tr # Compute panel mapping by finding closest unrefined section for each refined panel if recompute_mapping compute_refined_panel_mapping!(wing) - compute_refined_section_interpolation!(wing) + compute_refined_section_interpolation!(wing; reuse_aero_data) end # Update n_unrefined_sections based on actual sections @@ -1066,7 +1078,8 @@ function compute_refined_panel_mapping!(wing::AbstractWing) end """ - compute_refined_section_interpolation!(wing::AbstractWing) + compute_refined_section_interpolation!(wing::AbstractWing; + reuse_aero_data=false) Compute per-refined-section linear-interpolation weights from the unrefined sections. For refined section i, the interpolated value is: @@ -1078,8 +1091,15 @@ Positions are quarter-chord arc-length along the unrefined and refined sections. The first refined section is pinned to `left_idx == 1`, `weight == 1` (returns `unrefined[1]` exactly) and the last refined section to `left_idx == n_unref - 1`, `weight == 0` (returns `unrefined[end]` exactly). + +`reuse_aero_data` keeps the surface tables the refined sections already hold +instead of reblending them from the unrefined ones, the same preservation +[`refine!`](@ref) applies to `aero_data` under `use_prior_polar`. A remesh that +replaces the unrefined sections with a coarser set would otherwise resample the +surface tables down to that set even while the polars stay at full resolution. """ -function compute_refined_section_interpolation!(wing::AbstractWing{T}) where {T} +function compute_refined_section_interpolation!(wing::AbstractWing{T}; + reuse_aero_data::Bool=false) where {T} n_unref = length(wing.unrefined_sections) n_sections = wing.n_panels + 1 @@ -1147,7 +1167,8 @@ function compute_refined_section_interpolation!(wing::AbstractWing{T}) where {T} wing.refined_section_left_idx[n_sections] = Int16(n_unref - 1) wing.refined_section_weight[n_sections] = zero(T) - interpolate_section_aero_to_refined!(wing) + keep = reuse_aero_data && _can_reuse_prior_refined_surface_tables(wing) + keep || interpolate_section_aero_to_refined!(wing) return nothing end diff --git a/test/obj_adapter/test_obj_adapter.jl b/test/obj_adapter/test_obj_adapter.jl index 322113ed..e2f0ae91 100644 --- a/test/obj_adapter/test_obj_adapter.jl +++ b/test/obj_adapter/test_obj_adapter.jl @@ -90,6 +90,29 @@ obj_path = normpath(joinpath(@__DIR__, "..", "..", :csv) end + @testset "geometry_path writes the YAML outside the table directory" begin + root = mktempdir() + tables = joinpath(root, "tables") + yaml_path = joinpath(root, "nf_aero_geometry.yaml") + written = obj_to_yaml(obj_path, tables; n_sections=3, Re=5e5, + verbose=false, geometry_path=yaml_path) + @test written == yaml_path + @test isfile(yaml_path) + @test !isfile(joinpath(tables, "geometry.yaml")) + + info = Dict(YAML.load_file(yaml_path)["wing_airfoils"]["data"][1][3]) + for key in ("csv_file_path", "dat_file") + @test startswith(info[key], "tables/") + @test isfile(joinpath(dirname(yaml_path), info[key])) + end + # The loader resolves references against the YAML's own directory, so the + # prefixed paths have to be what makes this work. + @test Wing(yaml_path; n_panels=4) isa Wing + + @test ObjAdapter.table_path_prefix(joinpath(tables, "geometry.yaml"), + tables) == "" + end + @testset "center_to_com! rejects non-triangular faces" begin verts = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 1.0, 0.0]] @test_throws ArgumentError center_to_com!(verts, [[1, 2, 3, 4]]; prn=false) diff --git a/test/wing_geometry/test_wing_geometry.jl b/test/wing_geometry/test_wing_geometry.jl index 14f92408..c78a5dba 100644 --- a/test/wing_geometry/test_wing_geometry.jl +++ b/test/wing_geometry/test_wing_geometry.jl @@ -433,6 +433,60 @@ end @test wing_no_reuse.refined_sections[3].aero_data != no_reuse_baseline end + @testset "Surface tables across a coarsening remesh" begin + surface_table(scale) = VortexStepMethod.SectionAero( + deg2rad.([-5.0, 0.0, 5.0]), [0.0], + reshape(collect(range(0.0, 1.0; length=4)), 4, 1), + fill(0.05, 4, 1), fill(scale, 4, 3, 1), fill(0.1 * scale, 4, 3, 1)) + + alpha = deg2rad.([-5.0, 0.0, 5.0]) + polar = (collect(alpha), [1.0, 1.2, 1.4], [0.1, 0.11, 0.12], + [0.01, 0.02, 0.03]) + span_ys = [1.0, 0.5, 0.0, -0.5, -1.0] + scales = [1.0, 2.0, 10.0, 4.0, 5.0] + + function fine_wing(; use_prior_polar, tables=true) + wing = Wing(8; spanwise_distribution=LINEAR, use_prior_polar) + for (y, scale) in zip(span_ys, scales) + add_section!(wing, [0.0, y, 0.0], [1.0, y, 0.0], POLAR_VECTORS, + polar, tables ? surface_table(scale) : nothing) + end + refine!(wing) + return wing + end + + function keep_endpoints!(wing) + wing.unrefined_sections = [wing.unrefined_sections[1], + wing.unrefined_sections[end]] + wing.n_unrefined_sections = Int16(2) + refine!(wing; recompute_mapping=true, sort_sections=false) + return wing + end + + peak(wing) = maximum(maximum(s.section_aero.cp) + for s in wing.refined_sections) + + wing = fine_wing(use_prior_polar=true) + baseline = [copy(s.section_aero.cp) for s in wing.refined_sections] + @test peak(wing) > 5.0 + + keep_endpoints!(wing) + @test [s.section_aero.cp for s in wing.refined_sections] == baseline + # Only the dropped mid section carries a cp above 5; blending the two + # surviving endpoints cannot reach one. + @test peak(wing) > 5.0 + + @test peak(keep_endpoints!(fine_wing(use_prior_polar=false))) <= 5.0 + 1e-9 + + bare = fine_wing(use_prior_polar=true, tables=false) + @test all(isnothing(s.section_aero) for s in bare.refined_sections) + for (section, scale) in zip(bare.unrefined_sections, scales) + section.section_aero = surface_table(scale) + end + refine!(bare; recompute_mapping=true, sort_sections=false) + @test all(!isnothing(s.section_aero) for s in bare.refined_sections) + @test peak(bare) > 5.0 + end @testset "Refined panel mapping" begin # Test that refined panel mapping actually maps each panel to its closest unrefined panel From 912c86b18ffc2e9070257694f71356777b8459e9 Mon Sep 17 00:00:00 2001 From: Bart Date: Sun, 16 Aug 2026 12:57:52 +0200 Subject: [PATCH 2/8] Decide panel spanwise orientation once per wing MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit reinit! re-derived each panel's flip from dot(y_airf, spanwise_direction) on every call. A deforming wing can carry one panel's y_airf across that plane while its neighbours stay put, so the panel's z_airf inverted by 180° in a single step — and often back again the next one, chattering across the threshold. The traction pattern AeroPressure freezes is built on z_airf, so the aero force on that panel reversed with it. The orientation follows from the order a wing's sections run in, which is a property of its definition, so wing_span_flip computes it at BodyAerodynamics construction and every panel of that wing is reinitialized with it. Panels of one wing can no longer disagree and the value cannot change mid-run. reinit!(panel, …) drops spanwise_direction and takes flip::Bool instead: it no longer decides anything, it applies what the caller decided. Also carries the obj_to_yaml wingtip fix, the read_node_table speedup and the named-window plotting change that were already in the tree. Co-Authored-By: Claude Opus 5 (1M context) --- .gitignore | 1 + CHANGELOG.md | 23 ++++ docs/src/airfoil_pipeline.md | 9 +- docs/src/private_functions.md | 2 + examples/V3_kite.jl | 186 +++++++++++--------------------- ext/VortexStepMethodMakieExt.jl | 63 +++++++---- src/body_aerodynamics.jl | 30 ++++-- src/obj_adapter/obj_slice.jl | 33 +++--- src/obj_adapter/obj_to_yaml.jl | 7 +- src/panel.jl | 17 ++- test/panel/test_panel.jl | 3 +- 11 files changed, 193 insertions(+), 181 deletions(-) diff --git a/.gitignore b/.gitignore index a2266a9a..d7c105ab 100644 --- a/.gitignore +++ b/.gitignore @@ -46,6 +46,7 @@ LocalPreferences.toml # Pipeline-generated polars/dat/cp (regenerate via AirfoilAero/ObjAdapter) data/**/polars_neuralfoil/ data/**/polars_xfoil/ +data/**/generated_neuralfoil/ data/TUDELFT_V3_KITE/aero_geometry_neuralfoil.yaml # Local planning docs (not tracked) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8ad3593c..ad1c2d0d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,6 +32,22 @@ pressure integration reads. The reblend is now skipped when the polars are preserved and the refined sections already carry tables. +- `obj_to_yaml` no longer places sections on a wingtip that has closed to a point. + `station_indices` spreads its targets over the stations that still have a chord, + so the outermost section lands on the last sliceable one. A V3 mesh sliced with + the default `wingtip_distance` used to put a zero-chord section at each tip, + whose polar was `NaN` and took the whole solve with it; working around it meant + guessing a `wingtip_distance` large enough to skip past the tip. That workaround + is no longer the default: `wingtip_distance` is now `0.0`, an inset on top of the + trim for meshes whose slices just short of the tip are still too thin to analyse. + +- A panel's spanwise orientation is decided once per wing when `BodyAerodynamics` is + constructed (`wing_span_flip`) rather than per panel on every `reinit!`. Deciding it + from live geometry meant a deforming wing could carry one panel's `y_airf` across + `spanwise_direction` while its neighbours stayed put, inverting that panel's `z_airf` + by 180° in a single step — and inverting it back the next. Panels of one wing can no + longer disagree, and the orientation cannot change mid-run. + ### Changed - `SolverSettings` now defaults to the same values as `Solver`: `core_radius_fraction` `1e-20` → `0.05` and `type_initial_gamma_distribution` `ELLIPTIC` → `ZEROS`. The 0.05 @@ -50,9 +66,16 @@ already stores, take the core-radius cutoff from `|r1.r0|/|r0|` without forming the perpendicular vector, and defer the cross products that only one branch reads. Output is bit-identical; `solve!` is a further 1.47-1.55x faster. +- BREAKING: `reinit!(panel, …)` no longer takes `spanwise_direction`; it takes a + `flip::Bool` keyword instead, which the caller owns and must not derive from the + current geometry. - `read_node_table` parses into a preallocated matrix instead of `reduce(vcat, …)` over a generator, which was quadratic in the row count: ~21× faster on a 16 MB surface table (2.49 s → 0.12 s), benefiting every existing dataset. +- `is_show=true` draws into a window named after the plot title instead of into + whichever window the backend last used, so a script showing several plots gets + one window each and re-running it redraws them in place. `show_plot` takes the + window `name` as a keyword. ## VortexStepMethod v4.0.0 2026-08-03 diff --git a/docs/src/airfoil_pipeline.md b/docs/src/airfoil_pipeline.md index 0ad05192..63fe42d8 100644 --- a/docs/src/airfoil_pipeline.md +++ b/docs/src/airfoil_pipeline.md @@ -100,9 +100,12 @@ For each unique airfoil id `j`, `obj_to_yaml` writes into `output_dir`: - `geometry.yaml` — `wing_sections` (leading/trailing-edge points) plus `wing_airfoils` (each section's `type` and the `.dat`/`.csv` paths above) -A near-vanishing wingtip slice can shrink-wrap to an implausibly thick blob; such a -degenerate section reuses its nearest valid neighbour's airfoil and polar while keeping its -own edge positions, and a warning lists the reuse. All floats are rounded to millimetre +A tip that tapers to a point has no airfoil to slice, so the outermost stations stop at +the last slice that still has a chord; `wingtip_distance` moves them a further arc length +inboard when the slices just short of the tip are still too thin to analyse. A +near-vanishing slice that does get through can shrink-wrap to an implausibly thick blob; +such a degenerate section reuses its nearest valid neighbour's airfoil and polar while +keeping its own edge positions, and a warning lists the reuse. All floats are rounded to millimetre precision by the single [`write_yaml`](@ref VortexStepMethod.ObjAdapter.write_yaml) writer, so generated geometry files stay diff-friendly and consistent. diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 94c0f3e1..b0b0b82f 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -21,6 +21,7 @@ calculate_relative_alpha_and_velocity calculate_relative_alpha_and_relative_velocity update_effective_angle_of_attack! calculate_stall_angle_list +wing_span_flip calculate_circulation_distribution_elliptical_wing _compute_reference_velocity_from_distribution smooth_circulation! @@ -185,6 +186,7 @@ migrate_node_tables CurrentModule = Base.get_extension(VortexStepMethod, :VortexStepMethodMakieExt) ``` ```@docs +display_named create_geometry_plot_makie plot_line_segment_makie! set_axes_equal_makie! diff --git a/examples/V3_kite.jl b/examples/V3_kite.jl index 2c09c54b..4a7b831c 100644 --- a/examples/V3_kite.jl +++ b/examples/V3_kite.jl @@ -2,16 +2,21 @@ using Pkg if Base.active_project() != joinpath(@__DIR__, "Project.toml") Pkg.activate(@__DIR__) end -using LinearAlgebra using GLMakie using MakieControlPlots using VortexStepMethod +using VortexStepMethod.ObjAdapter +using VortexStepMethod.AirfoilAero: ShrinkWrap, NeuralFoilSolver PLOT = true SAVE_ALL = false USE_TEX = false DEFORM = false +NEURALFOIL = true +# Rolling-ball radius of the shrink wrap; fillets the concave tube-canopy junction. +MIN_CONCAVE_RADIUS = 0.4 OUTPUT_DIR = joinpath(dirname(@__DIR__), "output") +REFERENCE_POINT = [0.422646, 0.0, 9.3667] project_dir = dirname(@__DIR__) literature_paths = [ @@ -24,43 +29,21 @@ literature_paths = [ joinpath(project_dir, "data", "TUDELFT_V3_KITE", "literature_results", "windtunnel_alpha_sweep_beta_00_0_Poland_2025_Rey_5e5.csv"), ] -labels = [ - "VSM Julia Re=5e5", - "CFD Re=5e5", - "CFD Re=10e5", #with struts - "VSM Python Re=5e5", - "WindTunnel Re=5e5" #with struts -] beta_literature_paths = [ joinpath(project_dir, "data", "TUDELFT_V3_KITE", "literature_results", "windtunnel_beta_sweep_alpha_07_4_Poland_2025_Rey_5e5.csv"), ] -beta_labels = [ - labels[1], - "Wind Tunnel Re=5e5 beta sweep alpha=7.4", -] - -# Load YAML settings directly -settings_path = joinpath( - project_dir, "data", "TUDELFT_V3_KITE", "vsm_settings.yaml") -settings_data = VortexStepMethod.YAML.load_file(settings_path) -condition_cfg = settings_data["condition"] -wing_cfg = settings_data["wings"][1] -solver_cfg = settings_data["solver_settings"] -# Create wing, body_aero, and solver objects using settings -wing = Wing( - joinpath(project_dir, wing_cfg["geometry_file"]); - n_panels=wing_cfg["n_panels"], - spanwise_distribution=getproperty( - VortexStepMethod, - Symbol(wing_cfg["spanwise_panel_distribution"])), - spanwise_direction=Float64.(wing_cfg["spanwise_direction"]), - remove_nan=wing_cfg["remove_nan"], -) +settings = VSMSettings(joinpath(project_dir, "data", "TUDELFT_V3_KITE", + "vsm_settings.yaml"); data_prefix=false) +settings.wings[1].geometry_file = joinpath(project_dir, + settings.wings[1].geometry_file) +wing = Wing(settings) refine!(wing) body_aero = BodyAerodynamics([wing]) VortexStepMethod.reinit!(body_aero) +solver = Solver(body_aero, settings) +solver.reference_point .= REFERENCE_POINT if DEFORM VortexStepMethod.unrefined_deform!( @@ -72,52 +55,56 @@ if DEFORM VortexStepMethod.reinit!(body_aero; init_aero=false) end -# Construct Solver using keyword arguments from solver settings -solver = Solver(body_aero; - solver_type=(solver_cfg["solver_type"] == "NONLIN" ? NONLIN : LOOP), - aerodynamic_model_type=getproperty( - VortexStepMethod, - Symbol(solver_cfg["aerodynamic_model_type"])), - density=solver_cfg["density"], - max_iterations=solver_cfg["max_iterations"], - rtol=solver_cfg["rtol"], - tol_reference_error=solver_cfg["tol_reference_error"], - relaxation_factor=solver_cfg["relaxation_factor"], - is_with_artificial_damping=solver_cfg["artificial_damping"], - artificial_damping=(k2=solver_cfg["k2"], k4=solver_cfg["k4"]), - type_initial_gamma_distribution=getproperty( - VortexStepMethod, - Symbol(solver_cfg["type_initial_gamma_distribution"])), - use_gamma_prev=get(solver_cfg, "use_gamma_prev", - get(solver_cfg, "use_gamme_prev", true)), - core_radius_fraction=solver_cfg["core_radius_fraction"], - mu=solver_cfg["mu"], - is_only_f_and_gamma_output=get( - solver_cfg, "calc_only_f_and_gamma", false), - correct_aoa=get(solver_cfg, "correct_aoa", false), - reference_point=get(solver_cfg, "reference_point", - [0.422646, 0.0, 9.3667]), -) - -# Extract values for plotting -wind_speed = condition_cfg["wind_speed"] -angle_of_attack_deg = condition_cfg["alpha"] -sideslip_deg = condition_cfg["beta"] -yaw_rate = condition_cfg["yaw_rate"] - -# Set flight conditions from settings -α0 = deg2rad(angle_of_attack_deg) -β0 = deg2rad(sideslip_deg) -set_va!(body_aero, - wind_speed .* [cos(α0) * cos(β0), sin(β0), sin(α0) * cos(β0)]) +# Second sweep on generated polars: slice V3_25.obj, shrink-wrap every section into a +# closed airfoil and sweep it with NeuralFoil, so the same kite flies on polars derived +# from its own CAD surface instead of the checked-in CFD tables. +if NEURALFOIL + obj_file = joinpath(project_dir, "data", "TUDELFT_V3_KITE", "V3_25.obj") + generated_dir = joinpath(project_dir, "data", "TUDELFT_V3_KITE", + "generated_neuralfoil") + nf_yaml = obj_to_yaml(obj_file, generated_dir; + n_sections=settings.wings[1].n_panels, Re=1e6, force=false, + aero_solver=NeuralFoilSolver(model_size="large", n_crit=4.0, + xtr_upper=0.05, xtr_lower=0.05), + wrap_method=ShrinkWrap(clearance=0.0, + min_concave_radius=MIN_CONCAVE_RADIUS), + ) + settings_nf = deepcopy(settings) + settings_nf.wings[1].geometry_file = nf_yaml + wing_nf = Wing(settings_nf) + refine!(wing_nf) + body_nf = BodyAerodynamics([wing_nf]) + VortexStepMethod.reinit!(body_nf) + solver_nf = Solver(body_nf, settings_nf) + solver_nf.reference_point .= REFERENCE_POINT + + # Reading the generated directory instead of the OBJ shows the airfoils the polar + # pipeline actually analysed. Hover a slice to inspect its 2D fit. + PLOT && plot_slices_3d(generated_dir; obj_path=obj_file) +end -# Solve +solvers = NEURALFOIL ? [solver, solver_nf] : [solver] +bodies = NEURALFOIL ? [body_aero, body_nf] : [body_aero] +solver_labels = NEURALFOIL ? ["VSM Julia CFD", "VSM Julia NeuralFoil"] : + ["VSM Julia CFD"] +labels = [solver_labels; + ["CFD Re=5e5", + "CFD Re=10e5", #with struts + "VSM Python Re=5e5", + "WindTunnel Re=5e5"]] #with struts +beta_labels = [solver_labels; ["Wind Tunnel Re=5e5 beta sweep alpha=7.4"]] + +wind_speed = settings.condition.wind_speed +angle_of_attack_deg = settings.condition.alpha +sideslip_deg = settings.condition.beta +yaw_rate = settings.condition.yaw_rate + +set_va!(body_aero, settings) results = VortexStepMethod.solve(solver, body_aero; log=true) -# Plotting polars with moment coefficients PLOT && plot_polars( - [solver], - [body_aero], + solvers, + bodies, labels, literature_path_list=literature_paths, angle_range=range(-5, 25, length=31), @@ -130,7 +117,8 @@ PLOT && plot_polars( is_save=false || SAVE_ALL, is_show=true, use_tex=USE_TEX, - show_moments=true + show_moments=false, + cl_over_cd=true ) # Plotting geometry @@ -160,56 +148,10 @@ PLOT && plot_distribution( use_tex=USE_TEX ) -# --- Dual solver comparison: NONLIN vs LOOP --- -solver_cfg["solver_type"] = "LOOP" -solver_loop = Solver(body_aero; - solver_type=LOOP, - aerodynamic_model_type=getproperty( - VortexStepMethod, - Symbol(solver_cfg["aerodynamic_model_type"])), - density=solver_cfg["density"], - max_iterations=solver_cfg["max_iterations"], - rtol=solver_cfg["rtol"], - tol_reference_error=solver_cfg["tol_reference_error"], - relaxation_factor=solver_cfg["relaxation_factor"], - is_with_artificial_damping=solver_cfg["artificial_damping"], - artificial_damping=(k2=solver_cfg["k2"], k4=solver_cfg["k4"]), - type_initial_gamma_distribution=getproperty( - VortexStepMethod, - Symbol(solver_cfg["type_initial_gamma_distribution"])), - use_gamma_prev=get(solver_cfg, "use_gamma_prev", - get(solver_cfg, "use_gamme_prev", true)), - core_radius_fraction=solver_cfg["core_radius_fraction"], - mu=solver_cfg["mu"], - is_only_f_and_gamma_output=get( - solver_cfg, "calc_only_f_and_gamma", false), - correct_aoa=get(solver_cfg, "correct_aoa", false), - reference_point=get(solver_cfg, "reference_point", - [0.422646, 0.0, 9.3667]), -) - -PLOT && plot_polars( - [solver_loop], - [body_aero], - labels; - literature_path_list=literature_paths, - angle_range=range(-5, 20, step=1), - angle_type="angle_of_attack", - angle_of_attack=angle_of_attack_deg, - side_slip=sideslip_deg, - v_a=wind_speed, - title="LOOP solver", - show_moments=true, - save_path=OUTPUT_DIR, - is_save=false || SAVE_ALL, - is_show=true, - use_tex=USE_TEX -) - # --- Beta sweep --- PLOT && plot_polars( - [solver_loop], - [body_aero], + solvers, + bodies, beta_labels; literature_path_list=beta_literature_paths, angle_range=range(0, 12, step=1), @@ -225,4 +167,4 @@ PLOT && plot_polars( use_tex=USE_TEX ) -nothing \ No newline at end of file +nothing diff --git a/ext/VortexStepMethodMakieExt.jl b/ext/VortexStepMethodMakieExt.jl index dc6686f5..5e83bdb1 100644 --- a/ext/VortexStepMethodMakieExt.jl +++ b/ext/VortexStepMethodMakieExt.jl @@ -12,6 +12,32 @@ export plot_geometry, plot_distribution, plot_polars, save_plot, show_plot, const PANEL_MESH_OBSERVABLES = Ref{Union{Nothing,Dict}}(nothing) # Global storage for airfoil-skin observables, keyed by body objectid. const AIRFOIL_SKIN_OBSERVABLES = Ref{Union{Nothing,Dict}}(nothing) +const SCREENS = Dict{String,Any}() + +""" + display_named(fig, name) -> fig + +Show `fig` in the window registered under `name`, opening a new titled window for +a name that has none yet: plotting the same title again replaces its predecessor, +while a new title gets its own window. Backends without titled screens, such as +CairoMakie, fall back to a plain `display`. +""" +function display_named(fig::Makie.Figure, name::AbstractString) + screen = get(SCREENS, name, nothing) + if !isnothing(screen) && isopen(screen) + display(screen, fig) + return fig + end + screen = try + Makie.current_backend().Screen(; title=String(name)) + catch + display(fig) + return fig + end + SCREENS[name] = screen + display(screen, fig) + return fig +end """ PLATE_FACES @@ -465,18 +491,19 @@ function VortexStepMethod.save_plot(fig::Makie.Figure, save_path, title; data_ty end """ - show_plot(fig; dpi=130) + show_plot(fig; name="", dpi=130) -Display a Makie figure. +Display a Makie figure in an interactive session, in the window named `name`. # Arguments - `fig`: Makie Figure object # Keyword arguments +- `name`: Window to draw in; reusing a name reuses its window (default: "") - `dpi`: Dots per inch for the figure (default: 130) - currently unused in Makie """ -function VortexStepMethod.show_plot(fig::Makie.Figure; dpi=130) - isinteractive() && display(fig) +function VortexStepMethod.show_plot(fig::Makie.Figure; name="", dpi=130) + isinteractive() && display_named(fig, name) end """ @@ -671,9 +698,7 @@ function VortexStepMethod.plot_geometry(body_aero::BodyAerodynamics, title; fig = create_geometry_plot_makie(body_aero, title, view_elevation, view_azimuth) - if is_show && isinteractive() - display(fig) - end + is_show && show_plot(fig; name=title) return fig end @@ -801,9 +826,7 @@ function VortexStepMethod.plot_distribution(y_coordinates_list, results_list, la save_plot(fig, save_path, title, data_type=data_type) end - if is_show && isinteractive() - display(fig) - end + is_show && show_plot(fig; name=title) return fig end @@ -1023,9 +1046,7 @@ function VortexStepMethod.plot_polars( save_plot(fig, save_path, main_title; data_type) end - if is_show && isinteractive() - display(fig) - end + is_show && show_plot(fig; name=title) return fig end @@ -1082,9 +1103,7 @@ function VortexStepMethod.plot_polar_data(body_aero::BodyAerodynamics; color=:blue, linewidth=0.5, transparency=true) end - if is_show && isinteractive() - display(fig) - end + is_show && show_plot(fig; name="polar data") return fig else throw(ArgumentError( @@ -1454,9 +1473,7 @@ function VortexStepMethod.plot_combined_analysis( colsize!(fig.layout, 1, Relative(0.6)) colsize!(fig.layout, 2, Relative(0.4)) - if is_show && isinteractive() - display(fig) - end + is_show && show_plot(fig; name=title) return fig end @@ -1651,7 +1668,7 @@ function ObjAdapter.plot_slices_3d(path::String; n_slices::Int=10, rotation=I, m = ObjAdapter.march_edges(vertices, faces; step=span / n_bins) le = reduce(hcat, m.le) te = reduce(hcat, m.te) - idx = ObjAdapter.station_indices(m.arclen, n_slices; wingtip_distance) + idx = ObjAdapter.station_indices(m, n_slices; wingtip_distance) secs = filter(!isnothing, [ObjAdapter.build_section(vertices, faces, m.le[i], m.te[i], m.point[i], m.tangent[i]) for i in idx]) @@ -1726,7 +1743,7 @@ function ObjAdapter.plot_slices_3d(path::String; n_slices::Int=10, rotation=I, sel[] = best end - is_show && display(fig) + is_show && display_named(fig, "slices: $(basename(rstrip(path, '/')))") return fig end @@ -1780,7 +1797,7 @@ function ObjAdapter.plot_airfoil_fit(x::Vector, y::Vector; title::String="Airfoi axislegend(ax; position=:rt) - is_show && display(fig) + is_show && display_named(fig, title) return fig, params end @@ -1845,7 +1862,7 @@ function ObjAdapter.plot_airfoils(geometry_file::String; if is_save && !isnothing(save_path) VortexStepMethod.save_plot(fig, save_path, "airfoils"; data_type) end - is_show && display(fig) + is_show && display_named(fig, title) return fig end diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 87b514ec..629cd9ee 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -21,6 +21,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru - `projected_area::Float64` = 1.0: The area projected onto the xy-plane of the kite body reference frame [m²] - `c_ref::Float64` = 1.0: Reference chord length (max panel chord) [m] - `y::MVector{P, Float64}` = MVector{P,Float64}(zeros(P)) +- `span_flip::Vector{Int8}` = Int8[]: one orientation per wing, see [`wing_span_flip`](@ref) - `cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}}` = [LazyBufferCache() for _ in 1:15] """ @with_kw mutable struct BodyAerodynamics{P, W<:AbstractWing, T, PN<:Panel{T}} @@ -40,6 +41,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru projected_area::T = one(T) c_ref::T = one(T) y::MVector{P, T} = zeros(MVector{P, T}) + span_flip::Vector{Int8} = Int8[] cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:15] end @@ -113,11 +115,24 @@ function BodyAerodynamics( end end - body_aero = BodyAerodynamics{length(panels), W, T, eltype(panels)}(; panels, wings) + span_flip = Int8[wing_span_flip(wing) for wing in wings] + body_aero = BodyAerodynamics{length(panels), W, T, eltype(panels)}(; + panels, wings, span_flip) reinit!(body_aero; va, omega) return body_aero end +""" + wing_span_flip(wing) -> Int8 + +`-1` when `wing`'s sections run against its `spanwise_direction`, `+1` otherwise. The +`flip` every panel of the wing is reinitialized with ([`reinit!`](@ref)), decided once +here because deriving it per panel from live geometry inverts normals mid-run. +""" +wing_span_flip(wing) = + dot(first(wing.refined_sections).LE_point - last(wing.refined_sections).LE_point, + wing.spanwise_direction) < 0 ? Int8(-1) : Int8(1) + function Base.getproperty(obj::BodyAerodynamics, sym::Symbol) if sym === :va if getfield(obj, :has_distributed_va) @@ -248,12 +263,13 @@ function reinit!(body_aero::BodyAerodynamics{P, W, T}; ) where {P, W, T} idx = 1 vec = zeros(MVector{3, T}) - for wing in body_aero.wings + for (wing_idx, wing) in enumerate(body_aero.wings) reinit!(wing) validate_section_aero(wing.refined_sections) panel_props = wing.panel_props wing_init_aero = init_aero && !_can_skip_panel_aero_reinit(wing, body_aero.panels, idx) - + wing_flip = body_aero.span_flip[wing_idx] == -1 + # Create panels for i in 1:wing.n_panels if length(wing.delta_dist) > 0 @@ -263,7 +279,7 @@ function reinit!(body_aero::BodyAerodynamics{P, W, T}; delta = zero(T) end @views reinit!( - body_aero.panels[idx], + body_aero.panels[idx], wing.refined_sections[i], wing.refined_sections[i+1], panel_props.aero_centers[i, :], @@ -274,10 +290,10 @@ function reinit!(body_aero::BodyAerodynamics{P, W, T}; panel_props.y_airf[i, :], panel_props.z_airf[i, :], delta, - vec, - wing.spanwise_direction; + vec; remove_nan=wing.remove_nan, - init_aero=wing_init_aero + init_aero=wing_init_aero, + flip=wing_flip ) body_aero.panels[idx].crease_frac = wing.crease_frac idx += 1 diff --git a/src/obj_adapter/obj_slice.jl b/src/obj_adapter/obj_slice.jl index 592f958a..729895f7 100644 --- a/src/obj_adapter/obj_slice.jl +++ b/src/obj_adapter/obj_slice.jl @@ -433,6 +433,9 @@ built (`build_section`) at the marched station nearest each equal **leading-edge arc-length** target. Each section is `(; LE_point, TE_point, span_dir, contour3d, x_airfoil, y_airfoil)`. +Stations closed to a point at the tips are skipped ([`station_indices`](@ref)), and +`wingtip_distance` insets the outermost sections a further arc length. + The slicer assumes `x` = chordwise, `y` = spanwise, `z` = up. Pass a `3×3` rotation matrix to reorient a mesh stored in another convention before slicing. """ @@ -442,7 +445,7 @@ function perpendicular_sections(vertices, faces, n_sections; n_bins=60, rotation ys = [v[2] for v in vertices] m = march_edges(vertices, faces; step=(maximum(ys) - minimum(ys)) / n_bins) out = NamedTuple[] - for i in station_indices(m.arclen, n_sections; wingtip_distance) + for i in station_indices(m, n_sections; wingtip_distance) sec = build_section(vertices, faces, m.le[i], m.te[i], m.point[i], m.tangent[i]) sec === nothing || push!(out, sec) end @@ -450,17 +453,21 @@ function perpendicular_sections(vertices, faces, n_sections; n_bins=60, rotation end """ - station_indices(arclen, n; wingtip_distance=0.0) -> Vector{Int} - -Indices of the marched stations nearest `n` targets spread over the leading-edge arc -length. The first and last targets sit `wingtip_distance` (arc length) in from the -tips; with the default `0.0` they land exactly on the tips. Inset the tips a little -(e.g. `0.1` m) to avoid degenerate near-zero-chord tip sections that some solvers -(XFoil) cannot analyse. + station_indices(march, n; wingtip_distance=0.0, min_chord_frac=0.01) -> Vector{Int} + +Indices of the [`march_edges`](@ref) stations nearest `n` targets spread over the +leading-edge arc length. Stations whose chord has closed to less than +`min_chord_frac` of the longest one are left out of that range first, so a wing +tapering to a point puts its outermost sections on the last stations that still +have an airfoil to slice rather than on the point itself. The remaining first and +last targets sit a further `wingtip_distance` (arc length) inboard. """ -function station_indices(arclen, n; wingtip_distance=0.0) - total = arclen[end] - d = clamp(wingtip_distance, 0.0, total / 2) - n == 1 && return [argmin(abs.(arclen .- total / 2))] - return [argmin(abs.(arclen .- t)) for t in range(d, total - d, n)] +function station_indices(march, n; wingtip_distance=0.0, min_chord_frac=0.01) + chords = [norm(te .- le) for (le, te) in zip(march.le, march.te)] + usable = findall(≥(min_chord_frac * maximum(chords)), chords) + arclen = march.arclen + inner, outer = arclen[first(usable)], arclen[last(usable)] + d = clamp(wingtip_distance, 0.0, (outer - inner) / 2) + n == 1 && return [argmin(abs.(arclen .- (inner + outer) / 2))] + return [argmin(abs.(arclen .- t)) for t in range(inner + d, outer - d, n)] end diff --git a/src/obj_adapter/obj_to_yaml.jl b/src/obj_adapter/obj_to_yaml.jl index 27496fcd..5e8a58eb 100644 --- a/src/obj_adapter/obj_to_yaml.jl +++ b/src/obj_adapter/obj_to_yaml.jl @@ -45,7 +45,10 @@ Convert a 3D wing `.obj` mesh to the native YAML geometry route. Stations are placed at equal leading-edge arc-length intervals and sliced perpendicular to the local span (see [`perpendicular_sections`](@ref)), which keeps the airfoil undistorted near curved tips; each shape is then shrink-wrapped -into a clean airfoil and evaluated with `aero_solver`. +into a clean airfoil and evaluated with `aero_solver`. A tip that tapers to a +point carries no airfoil, so the outermost stations stop at the last slice that +still has a chord ([`station_indices`](@ref)); `wingtip_distance` moves them a +further arc length inboard. `aero_solver` selects the 2D-airfoil backend: [`NeuralFoilSolver`](@ref) (default, fast) or [`XFoilSolver`](@ref) (viscous panel code); pass `aero_solver=XFoilSolver()` @@ -188,7 +191,7 @@ function obj_to_yaml(obj_path::String, output_dir::String; wrap_method::ShrinkWrap=ShrinkWrap(), reuse_valid_airfoils::Bool=true, max_thickness_ratio::Real=2.0, spanwise_direction=[0.0, 1.0, 0.0], rotation=I, - wingtip_distance=0.05, crease_frac=0.75, force::Bool=false, + wingtip_distance=0.0, crease_frac=0.75, force::Bool=false, verbose::Bool=true, table_format::Symbol=:csv, geometry_path::String=joinpath(output_dir, "geometry.yaml")) (!endswith(obj_path, ".obj")) && (obj_path *= ".obj") diff --git a/src/panel.jl b/src/panel.jl index 7dd980f7..782ef330 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -227,14 +227,14 @@ end """ reinit!(panel, section_1, section_2, aero_center, control_point, bound_point_1, - bound_point_2, x_airf, y_airf, z_airf, delta, vec, spanwise_direction; kwargs...) + bound_point_2, x_airf, y_airf, z_airf, delta, vec; kwargs...) Reinitialize a panel's geometry, horseshoe filaments and aerodynamic interpolations. -The panel is oriented so its `y_airf` (and the bound vortex `bound_2 -> bound_1`) points -along `+spanwise_direction`, with `z_airf` pointing to the airfoil upper surface. This -makes the aero independent of section ordering: a reversed order would otherwise flip the -normal and make the panel look up its polar at a negated angle of attack. +`flip` reverses the section order so `y_airf` points along `+spanwise_direction` and +`z_airf` to the airfoil upper surface, making the aero independent of section ordering. +The caller owns it and must not derive it from the live geometry, which would invert +normals mid-run; [`reinit!(::BodyAerodynamics)`](@ref) decides it once per wing. """ function reinit!( panel::Panel, @@ -248,12 +248,11 @@ function reinit!( y_airf, z_airf, delta, - vec, - spanwise_direction; + vec; init_aero = true, - remove_nan = true + remove_nan = true, + flip::Bool = false ) - flip = dot(y_airf, spanwise_direction) < 0 if flip section_1, section_2 = section_2, section_1 bound_point_1, bound_point_2 = bound_point_2, bound_point_1 diff --git a/test/panel/test_panel.jl b/test/panel/test_panel.jl index 7dc014cb..eaeb6f8a 100644 --- a/test/panel/test_panel.jl +++ b/test/panel/test_panel.jl @@ -48,8 +48,7 @@ function create_panel(section1::Section, section2::Section) y_airf, z_airf, 0.0, - zeros(MVec3), - MVec3([0.0, 1.0, 0.0]) + zeros(MVec3) ) return panel end From 0d2120f5abf644a5f78dfc58660364548a9c09fa Mon Sep 17 00:00:00 2001 From: Bart Date: Sun, 16 Aug 2026 12:59:43 +0200 Subject: [PATCH 3/8] Drop the BREAKING note on the panel reinit! signature The panel method of reinit! takes internal panel_props slices and is not user-facing; users call reinit!(wing) and reinit!(body_aero). An internal signature change needs no changelog entry of its own. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 3 --- 1 file changed, 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ad1c2d0d..5964a86e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -66,9 +66,6 @@ already stores, take the core-radius cutoff from `|r1.r0|/|r0|` without forming the perpendicular vector, and defer the cross products that only one branch reads. Output is bit-identical; `solve!` is a further 1.47-1.55x faster. -- BREAKING: `reinit!(panel, …)` no longer takes `spanwise_direction`; it takes a - `flip::Bool` keyword instead, which the caller owns and must not derive from the - current geometry. - `read_node_table` parses into a preallocated matrix instead of `reduce(vcat, …)` over a generator, which was quadratic in the row count: ~21× faster on a 16 MB surface table (2.49 s → 0.12 s), benefiting every existing dataset. From 2ac530b508d1a9cb8390359ce3da997f244aee5a Mon Sep 17 00:00:00 2001 From: Bart Date: Sun, 16 Aug 2026 17:14:35 +0200 Subject: [PATCH 4/8] Normalize wing section order to +y -> -y on load MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Panel y_airf and z_airf are derived from the order sections are stored in, but that order was never an invariant: refine!'s sort was the only normalization, and callers that build particle wings pass sort_sections=false. A geometry file written -y -> +y therefore inverted every panel normal, and a wing whose sections were replaced after its panels were built inverted them mid-run, one panel at a time as each crossed spanwise_direction. normalize_span_order! settles it at the boundaries: on YAML load, before any structural pairing exists, and in refine! for wings assembled through add_section!. The refine! hook is gated on recompute_mapping so the per-step call that only updates positions still never reorders. sort_sections stays for a scrambled list. obj_to_yaml and surfplan_to_aero_yaml emit +y -> -y as well. Not breaking: files of either order load the same now. With the order canonical, nothing has to decide a per-panel or per-wing flip from live geometry, so BodyAerodynamics drops the span_flip field it cached — which serialized into model bins and went stale — and reinit! reads wing_span_flip directly as a guard for hand-built wings. Measured on a V3 beam wing whose aero geometry was written the wrong way round: 72 of 72 panel normals inverted after init! and 56 after bridle relaxation, both now 0, with the sections staying monotonic through relaxation where they did not before. Spanwise distribution plots put +y on the left to match, through one span_axis helper replacing nine repeated Axis constructions. Also carries the shrink_wrap rolling-ball rewrite, the V3_neuralfoil and ram_air_kite example updates already in the tree. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 26 +- docs/src/private_functions.md | 16 +- examples/V3_neuralfoil.jl | 4 +- examples/ram_air_kite.jl | 7 +- ext/VortexStepMethodMakieExt.jl | 39 +-- src/airfoil_aero/shrink_wrap.jl | 373 +++++++++--------------- src/body_aerodynamics.jl | 14 +- src/obj_adapter/obj_to_yaml.jl | 2 +- src/surfplan_adapter/SurfplanAdapter.jl | 2 +- src/wing_geometry.jl | 17 ++ src/yaml_geometry.jl | 5 +- 11 files changed, 222 insertions(+), 283 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5964a86e..55377a57 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -41,12 +41,15 @@ is no longer the default: `wingtip_distance` is now `0.0`, an inset on top of the trim for meshes whose slices just short of the tip are still too thin to analyse. -- A panel's spanwise orientation is decided once per wing when `BodyAerodynamics` is - constructed (`wing_span_flip`) rather than per panel on every `reinit!`. Deciding it - from live geometry meant a deforming wing could carry one panel's `y_airf` across - `spanwise_direction` while its neighbours stayed put, inverting that panel's `z_airf` - by 180° in a single step — and inverting it back the next. Panels of one wing can no - longer disagree, and the orientation cannot change mid-run. +- Wing sections are normalized to `+y` to `-y` order on load (`normalize_span_order!`), + and by `refine!` for wings built through `add_section!`. Panel `y_airf` and `z_airf` + follow the order sections are stored in, so a geometry file written the other way + round inverted every panel normal, and a wing whose sections were replaced after its + panels were built (a structural remesh) inverted them mid-run, one panel at a time as + each crossed `spanwise_direction`. `obj_to_yaml` and `surfplan_to_aero_yaml` emit that + order too; files of either order keep loading the same. +- Spanwise distribution plots put `+y` on the left, matching that order and the kite + seen from the front. ### Changed - `SolverSettings` now defaults to the same values as `Solver`: `core_radius_fraction` @@ -66,6 +69,17 @@ already stores, take the core-radius cutoff from `|r1.r0|/|r0|` without forming the perpendicular vector, and defer the cross products that only one branch reads. Output is bit-identical; `solve!` is a further 1.47-1.55x faster. +- `shrink_wrap` traces the rolling ball exactly — pivoting it around the cloud and + emitting the arcs its contact side sweeps (`pivot_contour`) — instead of thresholding + and marching-squares-tracing a distance field, so there is no grid resolution left to + set. `ShrinkWrap`'s `cell_size` is accordingly named `min_clearance`, still accepted + under the old name, and only floors `clearance` — which is what it already did for a + single-membrane slice. The wrap now + sits at exactly `clearance` from the cloud instead of a cell over it, so a V3 canopy's + aft strip comes out `2 * clearance` thick where the grid gave `3.5 * cell_size`. + `min_concave_radius` also stops costing anything, having padded the grid in both + directions before: one V3 slice at radius 0.4 drops from 173 ms to 10 ms, and at the + default radius from 15 ms to 3 ms. - `read_node_table` parses into a preallocated matrix instead of `reduce(vcat, …)` over a generator, which was quadratic in the row count: ~21× faster on a 16 MB surface table (2.49 s → 0.12 s), benefiting every existing dataset. diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index b0b0b82f..5ce84e26 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -45,7 +45,7 @@ update_panel_properties! build_interps panel_interp_types reinit!(wing::AbstractWing) -reinit!(panel::Panel, section_1::Section, section_2::Section, aero_center, control_point, bound_point_1, bound_point_2, x_airf, y_airf, z_airf, delta, vec, spanwise_direction) +reinit!(panel::Panel, section_1::Section, section_2::Section, aero_center, control_point, bound_point_1, bound_point_2, x_airf, y_airf, z_airf, delta, vec) rotated_te calculate_filaments_for_plotting ``` @@ -62,6 +62,7 @@ _apply_refined_section_thetas! _panel_thetas_to_section_thetas! _interpolate_unrefined_to_refined _section_sort_key +normalize_span_order! refine_mesh_for_linear_cosine_distribution! refine_mesh_by_splitting_provided_sections! refine_mesh_with_billowing! @@ -112,13 +113,13 @@ get_lower_upper turn_trailing_edge! ``` -### Shrink-wrap distance field +### Shrink-wrap rolling ball ```@docs -distance_parabolas! -squared_distance_transform! -grid_sampler -flood_outside -trace_level_set +point_buckets +turn_measure +pivot_step +push_arc! +pivot_contour largest_linking_gap resample_arc smooth_turning! @@ -187,6 +188,7 @@ CurrentModule = Base.get_extension(VortexStepMethod, :VortexStepMethodMakieExt) ``` ```@docs display_named +span_axis create_geometry_plot_makie plot_line_segment_makie! set_axes_equal_makie! diff --git a/examples/V3_neuralfoil.jl b/examples/V3_neuralfoil.jl index 2acedbca..88064a5b 100644 --- a/examples/V3_neuralfoil.jl +++ b/examples/V3_neuralfoil.jl @@ -46,8 +46,8 @@ ROTATION = I # Marched leading-edge stations across the span (finer => smoother edge trace). N_BINS = 100 -# Shrink-wrap each raw slice into a clean closed airfoil: the distance-field wrap -# hugs the cloud at `clearance` (floored at one grid `cell_size`); a single-skin +# Shrink-wrap each raw slice into a clean closed airfoil: the rolling-ball wrap +# hugs the cloud at `clearance` (floored at `min_clearance`); a single-skin # canopy becomes a thin capsule. `min_concave_radius` sets the fillet bridging the # concave tube-canopy junction: 0.2 smooths the neck over entirely; the default # (0.02) traces it, which NeuralFoil also handles fine. diff --git a/examples/ram_air_kite.jl b/examples/ram_air_kite.jl index 635b4fee..2eabe7fe 100644 --- a/examples/ram_air_kite.jl +++ b/examples/ram_air_kite.jl @@ -44,12 +44,11 @@ alpha_range = -8:2:26 delta_range = 0:1 N_SECTIONS = 10 -# Inset the tip stations 10 cm along the leading edge so they avoid the near-zero-chord -# wingtips, which slice to degenerate airfoils that XFoil cannot analyse. +# At 0.0 XFoil only converges above 16 deg on the near-tip slice this mesh still has. WINGTIP_DISTANCE = 0.1 -# Shrink-wrap each raw slice into a clean closed airfoil: the distance-field wrap -# hugs the cloud at `clearance` (floored at one grid `cell_size`) with a round nose. +# Shrink-wrap each raw slice into a clean closed airfoil: the rolling-ball wrap +# hugs the cloud at `clearance` (floored at `min_clearance`) with a round nose. WRAP = ShrinkWrap(clearance=0.0) # 3D slice diagnostic (live preview): mesh + LE/TE curves + section contours, their diff --git a/ext/VortexStepMethodMakieExt.jl b/ext/VortexStepMethodMakieExt.jl index 5e83bdb1..1356df73 100644 --- a/ext/VortexStepMethodMakieExt.jl +++ b/ext/VortexStepMethodMakieExt.jl @@ -703,6 +703,15 @@ function VortexStepMethod.plot_geometry(body_aero::BodyAerodynamics, title; return fig end +""" + span_axis(position, title, ylabel) -> Axis + +Axis for a spanwise distribution, `+y` on the left, matching the kite seen from the +front and the `+y` to `-y` order its sections are stored in. +""" +span_axis(position, title, ylabel) = + Axis(position; title, xlabel="Spanwise Position y/b", ylabel, xreversed=true) + """ plot_distribution(y_coordinates_list, results_list, label_list; title="spanwise_distribution", data_type=nothing, @@ -740,28 +749,22 @@ function VortexStepMethod.plot_distribution(y_coordinates_list, results_list, la Label(fig[0, :], title, fontsize=20) # Row 1: CL, CD, Gamma - ax_cl = Axis(fig[1, 1], title="CL Distribution", - xlabel="Spanwise Position y/b", ylabel="Lift Coefficient CL") - ax_cd = Axis(fig[1, 2], title="CD Distribution", - xlabel="Spanwise Position y/b", ylabel="Drag Coefficient CD") - ax_gamma = Axis(fig[1, 3], title="Γ Distribution", - xlabel="Spanwise Position y/b", ylabel="Circulation Γ") + ax_cl = span_axis(fig[1, 1], "CL Distribution", "Lift Coefficient CL") + ax_cd = span_axis(fig[1, 2], "CD Distribution", "Drag Coefficient CD") + ax_gamma = span_axis(fig[1, 3], "Γ Distribution", "Circulation Γ") # Row 2: Alpha geometric, alpha at ac, alpha uncorrected - ax_alpha_geo = Axis(fig[2, 1], title="α Geometric", - xlabel="Spanwise Position y/b", ylabel="Angle of Attack α (deg)") - ax_alpha_ac = Axis(fig[2, 2], title="α result (corrected to aerodynamic center)", - xlabel="Spanwise Position y/b", ylabel="Angle of Attack α (deg)") - ax_alpha_unc = Axis(fig[2, 3], title="α Uncorrected (if VSM, at control point)", - xlabel="Spanwise Position y/b", ylabel="Angle of Attack α (deg)") + alpha_label = "Angle of Attack α (deg)" + ax_alpha_geo = span_axis(fig[2, 1], "α Geometric", alpha_label) + ax_alpha_ac = span_axis(fig[2, 2], + "α result (corrected to aerodynamic center)", alpha_label) + ax_alpha_unc = span_axis(fig[2, 3], + "α Uncorrected (if VSM, at control point)", alpha_label) # Row 3: Force components - ax_fx = Axis(fig[3, 1], title="Force in x direction", - xlabel="Spanwise Position y/b", ylabel="Fx") - ax_fy = Axis(fig[3, 2], title="Force in y direction", - xlabel="Spanwise Position y/b", ylabel="Fy") - ax_fz = Axis(fig[3, 3], title="Force in z direction", - xlabel="Spanwise Position y/b", ylabel="Fz") + ax_fx = span_axis(fig[3, 1], "Force in x direction", "Fx") + ax_fy = span_axis(fig[3, 2], "Force in y direction", "Fy") + ax_fz = span_axis(fig[3, 3], "Force in z direction", "Fz") # Plot CL for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) diff --git a/src/airfoil_aero/shrink_wrap.jl b/src/airfoil_aero/shrink_wrap.jl index 944d9a9b..45b8bdac 100644 --- a/src/airfoil_aero/shrink_wrap.jl +++ b/src/airfoil_aero/shrink_wrap.jl @@ -1,212 +1,176 @@ """ - ShrinkWrap(; clearance=0.006, min_concave_radius=0.02, cell_size=0.001, + ShrinkWrap(; clearance=0.006, min_concave_radius=0.02, min_clearance=0.001, n_points=120, curvature_weight=0.05) -Distance-field shrink wrap: the rolling-ball offset of a raw slice point cloud, -extracted as one closed airfoil contour. A grid distance field is thresholded at -the rolling-ball radius (bridging cloud gaps and crevices narrower than about -twice `min_concave_radius`), flood-filled, eroded back to `clearance`, and the -resulting level set is traced with marching squares, faired (each pass clamped so -the contour keeps `clearance`) and resampled. Because the contour is -parameterized by arclength rather than `x`, the leading edge comes out genuinely -round — the offset of the cloud nose — the blunt trailing edge is capped by an arc -of radius `clearance`, and a single-membrane cloud becomes a thin capsule (at least -one `cell_size` half-thickness). +Rolling-ball shrink wrap: the offset of a raw slice point cloud, traced exactly as +one closed airfoil contour. A disk of radius `min_concave_radius` is pivoted around +the outside of the cloud ([`pivot_contour`](@ref)); the wrap is what its contact +side sweeps, so it consists of arcs of radius `clearance` about the points the disk +touches, joined by arcs of the disk radius across the gaps it cannot enter. Being +built from arcs rather than sampled on a grid, the wrap has no resolution floor: +the leading edge comes out genuinely round — the offset of the cloud nose — the +blunt trailing edge is capped by an arc of radius `clearance`, and a +single-membrane cloud becomes a capsule of exactly that half-thickness. # Fields -- `clearance`: offset the contour keeps outside every cloud point; floored at one - `cell_size` (the grid cannot represent a tighter wrap). Also the radius every - convex corner is rounded at. +- `clearance`: offset the contour keeps outside every cloud point, and the radius + every convex corner is rounded at; floored at `min_clearance`. - `min_concave_radius`: rolling-ball radius — concave features narrower than about twice this are bridged by a fillet of roughly this radius; convex geometry is unaffected. Auto-raised so the ball can neither fall through the cloud's largest - point gap nor pinch off between points during erosion, so sparse clouds get a - correspondingly looser, smoother wrap. -- `cell_size`: distance-field grid resolution as a chord fraction; sets the - geometric fidelity of the wrap. + point gap nor dip below `clearance` between neighbouring points, so sparse clouds + get a correspondingly looser, smoother wrap. +- `min_clearance`: floor on `clearance`, so `clearance=0` still leaves a single + membrane a capsule with two distinguishable sides instead of a bare curve. Also + accepted under its former name `cell_size`, which set the resolution of the + distance field this used to be traced on. - `n_points`: output stations per surface; the contour has `2*n_points - 1` points, cosine-clustered in arclength at the leading and trailing edges. - `curvature_weight`: extra sampling measure per radian of contour turning (chord fraction), concentrating output points into corners so XFoil's spline can follow them; `0` gives plain cosine-in-arclength sampling. """ -@with_kw struct ShrinkWrap - clearance::Float64 = 0.006 - min_concave_radius::Float64 = 0.02 - cell_size::Float64 = 0.001 - n_points::Int = 120 - curvature_weight::Float64 = 0.05 +struct ShrinkWrap + clearance::Float64 + min_concave_radius::Float64 + min_clearance::Float64 + n_points::Int + curvature_weight::Float64 +end + +function ShrinkWrap(; clearance=0.006, min_concave_radius=0.02, min_clearance=0.001, + cell_size=nothing, n_points=120, curvature_weight=0.05) + return ShrinkWrap(clearance, min_concave_radius, + isnothing(cell_size) ? min_clearance : cell_size, + n_points, curvature_weight) end """ - distance_parabolas!(d, f, n, v, z) + point_buckets(x, y, reach) -> Dict -One pass of the Felzenszwalb–Huttenlocher distance transform: writes into `d` the -lower envelope of the parabolas `(q - p)^2 + f[p]` over `p`, for `q in 1:n`. -`v` and `z` are scratch (length `n` and `n + 1`). +Sparse uniform grid of the point indices, one bucket per `reach`-sized cell, so every +point within `reach` of a query lies in one of the nine buckets around it. """ -function distance_parabolas!(d, f, n, v, z) - k = 1 - v[1] = 1 - z[1] = -Inf - z[2] = Inf - for q in 2:n - s = ((f[q] + q * q) - (f[v[k]] + v[k] * v[k])) / (2 * (q - v[k])) - while s <= z[k] - k -= 1 - s = ((f[q] + q * q) - (f[v[k]] + v[k] * v[k])) / (2 * (q - v[k])) - end - k += 1 - v[k] = q - z[k] = s - z[k+1] = Inf - end - k = 1 - for q in 1:n - while z[k+1] < q - k += 1 - end - d[q] = (q - v[k])^2 + f[v[k]] +function point_buckets(x, y, reach) + buckets = Dict{NTuple{2,Int},Vector{Int}}() + for k in eachindex(x) + key = (floor(Int, x[k] / reach), floor(Int, y[k] / reach)) + push!(get!(buckets, key, Int[]), k) end - return d + return buckets end """ - squared_distance_transform!(field) -> field + turn_measure(cross, dot) -> Float64 -In-place 2D squared Euclidean distance transform in grid-index units. On input -`field` holds `0` (or a local squared offset) at seeds and a large finite value -elsewhere; on output every node holds its squared distance to the nearest seed. +Pseudo-angle in `[0, 4)` of the direction with the given `cross` and `dot` against a +reference, rising monotonically with the counterclockwise angle it stands for. Ordering +directions by it is ordering them by angle, without the `atan` that dominates the +pivot's inner loop. """ -function squared_distance_transform!(field::Matrix{Float64}) - nx, ny = size(field) - m = max(nx, ny) - d, f = Vector{Float64}(undef, m), Vector{Float64}(undef, m) - v, z = Vector{Int}(undef, m), Vector{Float64}(undef, m + 1) - for i in 1:nx - f[1:ny] .= @view field[i, :] - distance_parabolas!(d, f, ny, v, z) - field[i, :] .= @view d[1:ny] - end - for j in 1:ny - f[1:nx] .= @view field[:, j] - distance_parabolas!(d, f, nx, v, z) - field[:, j] .= @view d[1:nx] +function turn_measure(cross, dot) + if cross >= 0 + return dot >= 0 ? cross / (cross + dot) : 1 + (-dot) / (cross - dot) end - return field + return dot <= 0 ? 2 + (-cross) / (-cross - dot) : 3 + dot / (dot - cross) end """ - flood_outside(blocked) -> BitMatrix + pivot_step(x, y, buckets, r, p, centre) -> (q, centre) or nothing -Mark the nodes reachable from the grid border without entering `blocked` -(4-connected flood fill); unreached nodes are the solid plus its enclosed holes. +Rotate the empty disk of radius `r` counterclockwise about point `p`, starting from +`centre`, until its rim reaches a second point. Returns that point and the disk +centre touching both, or `nothing` when nothing lies within `2r` of `p`. """ -function flood_outside(blocked::AbstractMatrix{Bool}) - nx, ny = size(blocked) - outside = falses(nx, ny) - stack = Int[] - visit(i, j) = if !blocked[i, j] && !outside[i, j] - outside[i, j] = true - push!(stack, (j - 1) * nx + i) - end - for i in 1:nx - visit(i, 1) - visit(i, ny) - end - for j in 1:ny - visit(1, j) - visit(nx, j) - end - while !isempty(stack) - idx = pop!(stack) - i, j = (idx - 1) % nx + 1, (idx - 1) ÷ nx + 1 - i > 1 && visit(i - 1, j) - i < nx && visit(i + 1, j) - j > 1 && visit(i, j - 1) - j < ny && visit(i, j + 1) +function pivot_step(x, y, buckets, r, p, centre) + px, py = x[p], y[p] + ux, uy = centre[1] - px, centre[2] - py + best_turn, best_q, best_centre = Inf, 0, centre + i0, j0 = floor(Int, px / 2r), floor(Int, py / 2r) + for i in i0-1:i0+1, j in j0-1:j0+1 + bucket = get(buckets, (i, j), nothing) + bucket === nothing && continue + for q in bucket + q == p && continue + dx, dy = x[q] - px, y[q] - py + span = dx * dx + dy * dy + (span < 1e-24 || span > 4r * r) && continue + out = sqrt(max(r * r / span - 0.25, 0.0)) + mx, my = dx / 2, dy / 2 + ox, oy = -dy * out, dx * out + for (vx, vy) in ((mx + ox, my + oy), (mx - ox, my - oy)) + (vx - ux)^2 + (vy - uy)^2 < 1e-24 && continue + turn = turn_measure(ux * vy - uy * vx, ux * vx + uy * vy) + (turn < 1e-9 || turn >= best_turn) && continue + best_turn, best_q, best_centre = turn, q, (px + vx, py + vy) + end + end end - return outside + best_q == 0 && return nothing + return best_q, best_centre end """ - trace_level_set(field, level) -> Vector{Vector{NTuple{2,Float64}}} + push_arc!(px, py, cx, cy, radius, from, sweep) -Marching-squares contours of `field .== level`, each returned as a closed loop of -grid-frame vertices (unit = one cell, node 1 at 1.0). Saddle cells are resolved by -the cell-center average. Loops touching the grid border are dropped. +Append the arc about `(cx, cy)` that starts at angle `from` and turns by the signed +angle `sweep`, stepped fine enough to stay within `1e-6` of the true arc. The closing +point is left off so consecutive arcs concatenate without duplicates. """ -function trace_level_set(field::Matrix{Float64}, level::Float64) - nx, ny = size(field) - inside(v) = v > level - frac(a, b) = (level - a) / (b - a) - points = Dict{NTuple{3,Int},NTuple{2,Float64}}() - links = Dict{NTuple{3,Int},Vector{NTuple{3,Int}}}() - connect(a, b) = begin - push!(get!(links, a, NTuple{3,Int}[]), b) - push!(get!(links, b, NTuple{3,Int}[]), a) +function push_arc!(px, py, cx, cy, radius, from, sweep) + steps = radius < 1e-12 ? 1 : + clamp(ceil(Int, abs(sweep) / sqrt(8e-6 / radius)), 1, 400) + for k in 0:steps-1 + angle = from + sweep * k / steps + push!(px, cx + radius * cos(angle)) + push!(py, cy + radius * sin(angle)) end - for j in 1:ny-1, i in 1:nx-1 - v00, v10 = field[i, j], field[i+1, j] - v01, v11 = field[i, j+1], field[i+1, j+1] - b00, b10, b01, b11 = inside(v00), inside(v10), inside(v01), inside(v11) - b00 == b10 == b01 == b11 && continue - bottom, top = (i, j, 0), (i, j + 1, 0) - left, right = (i, j, 1), (i + 1, j, 1) - b00 != b10 && (points[bottom] = (i + frac(v00, v10), Float64(j))) - b01 != b11 && (points[top] = (i + frac(v01, v11), Float64(j + 1))) - b00 != b01 && (points[left] = (Float64(i), j + frac(v00, v01))) - b10 != b11 && (points[right] = (Float64(i + 1), j + frac(v10, v11))) - crossed = [e for (e, c) in ((bottom, b00 != b10), (top, b01 != b11), - (left, b00 != b01), (right, b10 != b11)) if c] - if length(crossed) == 4 - if inside((v00 + v10 + v01 + v11) / 4) == b00 - connect(bottom, right) - connect(top, left) - else - connect(bottom, left) - connect(top, right) - end - else - connect(crossed[1], crossed[2]) - end - end - visited = Set{NTuple{3,Int}}() - loops = Vector{Vector{NTuple{2,Float64}}}() - for start in keys(links) - start in visited && continue - loop = NTuple{2,Float64}[] - prev, current = start, start - closed = false - while true - push!(visited, current) - push!(loop, points[current]) - nbrs = links[current] - length(nbrs) == 2 || break - prev, current = current, (nbrs[1] == prev ? nbrs[2] : nbrs[1]) - current == start && (closed = true; break) - end - closed && push!(loops, loop) - end - return loops + return nothing end """ - grid_sampler(field, x0, y0, cell) -> f(px, py) + pivot_contour(x, y, r, clearance) -> (px, py) -Bilinear interpolant of `field` at physical points; node `(i, j)` sits at -`(x0 + (i-1)*cell, y0 + (j-1)*cell)`. Queries are clamped to the grid. +Trace the rolling-ball closing of the cloud `(x, y)` offset outward by `clearance`, +as a dense closed polyline. A disk of radius `r` is pivoted around the outside of the +cloud from its leftmost point; each point the disk touches contributes an arc of +radius `clearance` about itself, and each pivot an arc of radius `r - clearance` +about the empty disk's centre, spanning the gap the disk could not enter. Those arcs +are the exact wrap boundary, so nothing here has a resolution to choose. """ -function grid_sampler(field::Matrix{Float64}, x0, y0, cell) - nx, ny = size(field) - return function (px, py) - gi = clamp((px - x0) / cell + 1, 1.0, nx - 1e-9) - gj = clamp((py - y0) / cell + 1, 1.0, ny - 1e-9) - i, j = floor(Int, gi), floor(Int, gj) - ti, tj = gi - i, gj - j - return (field[i, j] * (1 - ti) + field[i+1, j] * ti) * (1 - tj) + - (field[i, j+1] * (1 - ti) + field[i+1, j+1] * ti) * tj +function pivot_contour(x, y, r, clearance) + start = argmin(x) + buckets = point_buckets(x, y, 2r) + contacts, centres = Int[], NTuple{2,Float64}[] + p, centre = start, (x[start] - r, y[start]) + # A single-membrane stretch is touched once from each side, so the cycle closes on + # the (point, centre) pair rather than on the point alone. + for _ in 1:4*length(x)+100 + step = pivot_step(x, y, buckets, r, p, centre) + step === nothing && break + next, centre = step + !isempty(contacts) && p == contacts[1] && + hypot(centre[1] - centres[1][1], centre[2] - centres[1][2]) < 1e-12 && + break + push!(contacts, p) + push!(centres, centre) + p = next + end + px, py = Float64[], Float64[] + isempty(contacts) && return px, py + for (k, i) in enumerate(contacts) + incoming = centres[mod1(k - 1, length(centres))] + outgoing = centres[k] + from = atan(incoming[2] - y[i], incoming[1] - x[i]) + to = atan(outgoing[2] - y[i], outgoing[1] - x[i]) + push_arc!(px, py, x[i], y[i], clearance, from, mod(to - from, 2pi)) + j = contacts[mod1(k + 1, length(contacts))] + from = atan(y[i] - outgoing[2], x[i] - outgoing[1]) + to = atan(y[j] - outgoing[2], x[j] - outgoing[1]) + push_arc!(px, py, outgoing[1], outgoing[2], r - clearance, from, + rem(to - from, 2pi, RoundNearest)) end + return px, py end """ @@ -303,11 +267,11 @@ end shrink_wrap(x, y, method::ShrinkWrap) -> (x, y) Wrap the point cloud `(x, y)` into a clean closed airfoil in Selig order (TE upper → -LE → TE lower), following [`ShrinkWrap`](@ref): distance field on a `cell_size` -grid, closing with the rolling ball (`min_concave_radius`), offset outward by -`clearance`, traced as a single -closed contour and resampled to cosine panels in a curvature-weighted arclength -measure. The first and last point coincide at the +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). 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 @@ -315,73 +279,14 @@ normalized frame of the input cloud (chord slightly longer than 1, nose apex nea """ function shrink_wrap(x, y, method::ShrinkWrap) xn, yn, _ = normalize_airfoil(collect(float.(x)), collect(float.(y))) - cell = method.cell_size - gap = max(method.clearance, cell) + gap = max(method.clearance, method.min_clearance) linking = largest_linking_gap(xn, yn) - ball = max(method.min_concave_radius, gap + 3cell, + ball = max(method.min_concave_radius, gap, 1.01 * linking / 2, 1.5 * (linking^2 / 4 + gap^2) / (2gap)) - pad = ball + 3cell - x0, y0 = -pad, minimum(yn) - pad - nx = ceil(Int, (1.0 + pad - x0) / cell) + 1 - ny = ceil(Int, (maximum(yn) + pad - y0) / cell) + 1 - - far = Float64(nx^2 + ny^2) - field = fill(far, nx, ny) - for (px, py) in zip(xn, yn) - gi, gj = (px - x0) / cell + 1, (py - y0) / cell + 1 - for i in (floor(Int, gi), floor(Int, gi) + 1), - j in (floor(Int, gj), floor(Int, gj) + 1) - - (1 <= i <= nx && 1 <= j <= ny) || continue - field[i, j] = min(field[i, j], (gi - i)^2 + (gj - j)^2) - end - end - squared_distance_transform!(field) - cloud_dist = sqrt.(field) .* cell - - outside = flood_outside(cloud_dist .<= ball) - depth = [outside[i, j] ? 0.0 : far for i in 1:nx, j in 1:ny] - squared_distance_transform!(depth) - depth .= sqrt.(depth) .* cell - - loops = trace_level_set(depth, ball - gap) - isempty(loops) && error("ShrinkWrap traced no contour; the cloud may be too" * - " sparse or thin for cell_size=$(cell).") - shoelace(loop) = abs(sum(p[1] * loop[mod1(k + 1, length(loop))][2] - - loop[mod1(k + 1, length(loop))][1] * p[2] - for (k, p) in enumerate(loop))) / 2 - areas = shoelace.(loops) - main = loops[argmax(areas)] - if length(loops) > 1 && sort(areas)[end-1] > 0.01 * maximum(areas) - @warn "ShrinkWrap discarded a comparable secondary contour; increase" * - " min_concave_radius to bridge cloud gaps." - end - px = [x0 + (p[1] - 1) * cell for p in main] - py = [y0 + (p[2] - 1) * cell for p in main] - - dist_at = grid_sampler(cloud_dist, x0, y0, cell) - keep_clear(cx, cy) = begin - d = dist_at(cx, cy) - d >= gap && return (cx, cy) - h = cell / 2 - gx = (dist_at(cx + h, cy) - dist_at(cx - h, cy)) / cell - gy = (dist_at(cx, cy + h) - dist_at(cx, cy - h)) / cell - len = hypot(gx, gy) - len < 1e-9 && return (cx, cy) - return (cx + (gap - d) * gx / len, cy + (gap - d) * gy / len) - end + px, py = pivot_contour(xn, yn, ball, gap) + length(px) < 3 && error("ShrinkWrap traced no contour; the cloud may be too" * + " sparse or thin to pivot a ball of $(ball) around.") m = length(px) - for pass in 0:20 - ox, oy = copy(px), copy(py) - for k in 1:m - if pass > 0 - a, b = mod1(k - 1, m), mod1(k + 1, m) - px[k] = ox[k] + 0.4 * (ox[a] - 2ox[k] + ox[b]) - py[k] = oy[k] + 0.4 * (oy[a] - 2oy[k] + oy[b]) - end - px[k], py[k] = keep_clear(px[k], py[k]) - end - end anchor = argmax(px) px, py = circshift(px, 1 - anchor), circshift(py, 1 - anchor) diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 629cd9ee..a25543a7 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -21,7 +21,6 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru - `projected_area::Float64` = 1.0: The area projected onto the xy-plane of the kite body reference frame [m²] - `c_ref::Float64` = 1.0: Reference chord length (max panel chord) [m] - `y::MVector{P, Float64}` = MVector{P,Float64}(zeros(P)) -- `span_flip::Vector{Int8}` = Int8[]: one orientation per wing, see [`wing_span_flip`](@ref) - `cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}}` = [LazyBufferCache() for _ in 1:15] """ @with_kw mutable struct BodyAerodynamics{P, W<:AbstractWing, T, PN<:Panel{T}} @@ -41,7 +40,6 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru projected_area::T = one(T) c_ref::T = one(T) y::MVector{P, T} = zeros(MVector{P, T}) - span_flip::Vector{Int8} = Int8[] cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:15] end @@ -115,9 +113,7 @@ function BodyAerodynamics( end end - span_flip = Int8[wing_span_flip(wing) for wing in wings] - body_aero = BodyAerodynamics{length(panels), W, T, eltype(panels)}(; - panels, wings, span_flip) + body_aero = BodyAerodynamics{length(panels), W, T, eltype(panels)}(; panels, wings) reinit!(body_aero; va, omega) return body_aero end @@ -125,9 +121,9 @@ end """ wing_span_flip(wing) -> Int8 -`-1` when `wing`'s sections run against its `spanwise_direction`, `+1` otherwise. The -`flip` every panel of the wing is reinitialized with ([`reinit!`](@ref)), decided once -here because deriving it per panel from live geometry inverts normals mid-run. +`-1` when `wing`'s sections run against its `spanwise_direction`, `+1` otherwise: the +`flip` every panel of the wing is reinitialized with ([`reinit!`](@ref)). One answer per +wing, so neighbouring panels cannot disagree and invert a single normal by 180°. """ wing_span_flip(wing) = dot(first(wing.refined_sections).LE_point - last(wing.refined_sections).LE_point, @@ -268,7 +264,7 @@ function reinit!(body_aero::BodyAerodynamics{P, W, T}; validate_section_aero(wing.refined_sections) panel_props = wing.panel_props wing_init_aero = init_aero && !_can_skip_panel_aero_reinit(wing, body_aero.panels, idx) - wing_flip = body_aero.span_flip[wing_idx] == -1 + wing_flip = wing_span_flip(wing) == -1 # Create panels for i in 1:wing.n_panels diff --git a/src/obj_adapter/obj_to_yaml.jl b/src/obj_adapter/obj_to_yaml.jl index 5e8a58eb..bc05908f 100644 --- a/src/obj_adapter/obj_to_yaml.jl +++ b/src/obj_adapter/obj_to_yaml.jl @@ -252,7 +252,7 @@ function obj_to_yaml(obj_path::String, output_dir::String; s.TE_point[1], s.TE_point[2], s.TE_point[3]]) end - sort!(section_rows; by = row -> row[3]) # clean spanwise order (by LE_y) + sort!(section_rows; by = row -> row[3], rev = true) # +y to -y, by LE_y sort!(airfoil_rows; by = row -> row[1]) # airfoils by id write_geometry_yaml(yaml_path, section_rows, airfoil_rows) verbose && @info "Wrote geometry to $yaml_path ($(length(section_rows)) sections)" diff --git a/src/surfplan_adapter/SurfplanAdapter.jl b/src/surfplan_adapter/SurfplanAdapter.jl index e1cbe40c..a723daab 100644 --- a/src/surfplan_adapter/SurfplanAdapter.jl +++ b/src/surfplan_adapter/SurfplanAdapter.jl @@ -83,7 +83,7 @@ function surfplan_to_aero_yaml(adapter_dir::AbstractString, output_dir::Abstract push!(section_rows, Any[fid, s.LE_point[1], s.LE_point[2], s.LE_point[3], s.TE_point[1], s.TE_point[2], s.TE_point[3]]) end - sort!(section_rows; by = row -> row[3]) + sort!(section_rows; by = row -> row[3], rev = true) # +y to -y, by LE_y sort!(airfoil_rows; by = row -> row[1]) write_geometry_yaml(yaml_path, section_rows, airfoil_rows) verbose && @info "Wrote pressure geometry to $yaml_path " * diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index 25c0251b..8245671e 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -56,6 +56,21 @@ Function to update a [Section](@ref) in place. """ @inline _section_sort_key(s::Section) = s.LE_point[2] +"""Spanwise coordinate [`normalize_span_order!`](@ref) orders on.""" +_span_order_key(s::Section) = _section_sort_key(s) + +""" + normalize_span_order!(sections) -> sections + +Reverse `sections` if they do not already run `+y` to `-y`, the order panel normals +are built from. Use `refine!`'s `sort_sections` for a scrambled list. +""" +function normalize_span_order!(sections) + length(sections) > 1 && _span_order_key(last(sections)) > + _span_order_key(first(sections)) && reverse!(sections) + return sections +end + function reinit!(section::Section, LE_point, TE_point, aero_model=nothing, aero_data=nothing, section_aero=nothing) section.LE_point .= LE_point @@ -925,6 +940,8 @@ function refine!(wing::AbstractWing{T}; recompute_mapping=true, sort_sections=tr end sorted || sort!(wing.unrefined_sections; by=_section_sort_key, rev=true) + elseif recompute_mapping + normalize_span_order!(wing.unrefined_sections) end n_sections = wing.n_panels + 1 diff --git a/src/yaml_geometry.jl b/src/yaml_geometry.jl index 94e101df..ee57f206 100644 --- a/src/yaml_geometry.jl +++ b/src/yaml_geometry.jl @@ -19,6 +19,8 @@ end TE_z::Float64 end +_span_order_key(section::WingSectionData) = section.LE_y + @with_kw struct WingAirfoilData airfoil_id::Int64 type::String @@ -247,7 +249,8 @@ function Wing( TE_z = section_dict["TE_z"] )) end - + normalize_span_order!(sections) + # Convert wing airfoils wing_airfoils_data = data["wing_airfoils"] airfoils = WingAirfoilData[] From 5f27ade020a2f31b731c9900bd16c7a777c5798f Mon Sep 17 00:00:00 2001 From: Bart Date: Mon, 17 Aug 2026 09:29:50 +0200 Subject: [PATCH 5/8] Different field for corrected vs uncorrected --- src/body_aerodynamics.jl | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index a25543a7..b9a56ca5 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -16,8 +16,11 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru - `alpha_dist::MVector{P, Float64}` = zeros(Float64, P) - `v_a_dist::MVector{P, Float64}` = zeros(Float64, P) - `work_vectors`::NTuple{10, MVec3} = ntuple(_ -> zeros(MVec3), 10) -- `AIC::Array{Float64, 3}` = zeros(P, P, 3): influence coefficients, component last so - that each `AIC[:, :, k]` slice is a contiguous BLAS matrix +- `AIC::Array{Float64, 3}` = zeros(P, P, 3): control-point influence coefficients, the + matrix the circulation is solved against; component last so that + each `AIC[:, :, k]` slice is a contiguous BLAS matrix +- `AIC_aero_center::Array{Float64, 3}` = zeros(P, P, 3): aerodynamic-centre (LLT) + influence coefficients, used only for the corrected angle of attack - `projected_area::Float64` = 1.0: The area projected onto the xy-plane of the kite body reference frame [m²] - `c_ref::Float64` = 1.0: Reference chord length (max panel chord) [m] - `y::MVector{P, Float64}` = MVector{P,Float64}(zeros(P)) @@ -37,6 +40,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru v_a_dist::MVector{P, T} = zeros(MVector{P, T}) work_vectors::NTuple{10, MVector{3, T}} = ntuple(_ -> zeros(MVector{3, T}), 10) AIC::Array{T, 3} = zeros(T, P, P, 3) + AIC_aero_center::Array{T, 3} = zeros(T, P, P, 3) projected_area::T = one(T) c_ref::T = one(T) y::MVector{P, T} = zeros(MVector{P, T}) @@ -382,7 +386,8 @@ Returns: nothing @inline function calculate_AIC_matrices!(body_aero::BodyAerodynamics{P, W, T}, model::Model, core_radius_fraction, va_norm_array::AbstractVector{T}, - va_unit_array::AbstractMatrix{T}) where {P, W, T} + va_unit_array::AbstractMatrix{T}, + target::AbstractArray{T, 3}=body_aero.AIC) where {P, W, T} # Determine evaluation point based on model evaluation_point = model == VSM ? :control_point : :aero_center evaluation_point_on_bound = model == LLT @@ -436,7 +441,7 @@ Returns: nothing velocity_induced .-= U_2D end @inbounds for k in 1:3 - body_aero.AIC[icp, jring, k] = velocity_induced[k] + target[icp, jring, k] = velocity_induced[k] end end end @@ -494,11 +499,14 @@ function update_effective_angle_of_attack!(alpha_corrected, va_norm_array, va_unit_array) - calculate_AIC_matrices!(body_aero, LLT, core_radius_fraction, va_norm_array, va_unit_array) + # Its own buffer: `AIC` holds the control-point matrix the circulation was solved + # against, so overwriting it here would leave post-solve readers on the LLT one. + calculate_AIC_matrices!(body_aero, LLT, core_radius_fraction, va_norm_array, + va_unit_array, body_aero.AIC_aero_center) induced_velocity = body_aero.cache[1][va_array] for k in 1:3 - mul!(view(induced_velocity, :, k), view(body_aero.AIC, :, :, k), gamma) + mul!(view(induced_velocity, :, k), view(body_aero.AIC_aero_center, :, :, k), gamma) end # In-place relative velocity calculation From 89641ddf5a928550b711d4df205cf25f42b8ba5f Mon Sep 17 00:00:00 2001 From: Bart Date: Mon, 17 Aug 2026 11:45:21 +0200 Subject: [PATCH 6/8] Document the new span-order and table-path helpers Documenter's missing_docs check fails the docs build on any docstring in a checked module that no @docs block lists, and four new ones were not listed. span_order_key subsumes _section_sort_key, which named the same coordinate for the same purpose, and both it and can_reuse_prior_refined_surface_tables lose the underscore prefix now that they are documented names. Co-Authored-By: Claude Opus 5 (1M context) --- docs/src/private_functions.md | 6 +++++- src/wing_geometry.jl | 32 ++++++++++++++++++-------------- src/yaml_geometry.jl | 2 +- 3 files changed, 24 insertions(+), 16 deletions(-) diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 5ce84e26..9ea565b8 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -45,6 +45,7 @@ update_panel_properties! build_interps panel_interp_types reinit!(wing::AbstractWing) +reinit!(section::Section, LE_point, TE_point, aero_model, aero_data, section_aero) reinit!(panel::Panel, section_1::Section, section_2::Section, aero_center, control_point, bound_point_1, bound_point_2, x_airf, y_airf, z_airf, delta, vec) rotated_te calculate_filaments_for_plotting @@ -61,8 +62,9 @@ copy_sections_to_refined! _apply_refined_section_thetas! _panel_thetas_to_section_thetas! _interpolate_unrefined_to_refined -_section_sort_key +span_order_key normalize_span_order! +can_reuse_prior_refined_surface_tables refine_mesh_for_linear_cosine_distribution! refine_mesh_by_splitting_provided_sections! refine_mesh_with_billowing! @@ -178,6 +180,8 @@ center_to_com! airfoils_from_yaml write_geometry_yaml resolve_aero_geometry +table_path_prefix +prefix_table_paths! plot_airfoil_fit migrate_node_tables ``` diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index 8245671e..247e5db6 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -50,14 +50,12 @@ function Section(LE_point, TE_point, aero_model, aero_data, section_aero=nothing end """ - reinit!(section::Section, LE_point::PosVector, TE_point::PosVector, aero_model=nothing, aero_data=nothing) + span_order_key(section) -> Float64 -Function to update a [Section](@ref) in place. +Spanwise coordinate that [`normalize_span_order!`](@ref) and `refine!`'s +`sort_sections` order sections on. """ -@inline _section_sort_key(s::Section) = s.LE_point[2] - -"""Spanwise coordinate [`normalize_span_order!`](@ref) orders on.""" -_span_order_key(s::Section) = _section_sort_key(s) +@inline span_order_key(s::Section) = s.LE_point[2] """ normalize_span_order!(sections) -> sections @@ -66,11 +64,17 @@ Reverse `sections` if they do not already run `+y` to `-y`, the order panel norm are built from. Use `refine!`'s `sort_sections` for a scrambled list. """ function normalize_span_order!(sections) - length(sections) > 1 && _span_order_key(last(sections)) > - _span_order_key(first(sections)) && reverse!(sections) + length(sections) > 1 && span_order_key(last(sections)) > + span_order_key(first(sections)) && reverse!(sections) return sections end +""" + reinit!(section::Section, LE_point, TE_point, aero_model=nothing, aero_data=nothing, + section_aero=nothing) + +Update a [Section](@ref) in place. +""" function reinit!(section::Section, LE_point, TE_point, aero_model=nothing, aero_data=nothing, section_aero=nothing) section.LE_point .= LE_point @@ -828,13 +832,13 @@ end end """ - _can_reuse_prior_refined_surface_tables(wing) -> Bool + can_reuse_prior_refined_surface_tables(wing) -> Bool Whether every refined section already carries a [`SectionAero`](@ref), so a remesh can keep them instead of reblending from the unrefined sections. False on a wing that has none, where the blend is what fills them in the first place. """ -@inline function _can_reuse_prior_refined_surface_tables(wing::AbstractWing) +@inline function can_reuse_prior_refined_surface_tables(wing::AbstractWing) isempty(wing.refined_sections) && return false return all(s -> !isnothing(s.section_aero), wing.refined_sections) end @@ -932,14 +936,14 @@ function refine!(wing::AbstractWing{T}; recompute_mapping=true, sort_sections=tr if sort_sections sorted = true for i in 1:(length(wing.unrefined_sections) - 1) - if _section_sort_key(wing.unrefined_sections[i]) < - _section_sort_key(wing.unrefined_sections[i+1]) + if span_order_key(wing.unrefined_sections[i]) < + span_order_key(wing.unrefined_sections[i+1]) sorted = false break end end sorted || sort!(wing.unrefined_sections; - by=_section_sort_key, rev=true) + by=span_order_key, rev=true) elseif recompute_mapping normalize_span_order!(wing.unrefined_sections) end @@ -1184,7 +1188,7 @@ function compute_refined_section_interpolation!(wing::AbstractWing{T}; wing.refined_section_left_idx[n_sections] = Int16(n_unref - 1) wing.refined_section_weight[n_sections] = zero(T) - keep = reuse_aero_data && _can_reuse_prior_refined_surface_tables(wing) + keep = reuse_aero_data && can_reuse_prior_refined_surface_tables(wing) keep || interpolate_section_aero_to_refined!(wing) return nothing end diff --git a/src/yaml_geometry.jl b/src/yaml_geometry.jl index ee57f206..11da48ac 100644 --- a/src/yaml_geometry.jl +++ b/src/yaml_geometry.jl @@ -19,7 +19,7 @@ end TE_z::Float64 end -_span_order_key(section::WingSectionData) = section.LE_y +span_order_key(section::WingSectionData) = section.LE_y @with_kw struct WingAirfoilData airfoil_id::Int64 From 20309082422dd3cfde03c8173904ab2c2acf2b98 Mon Sep 17 00:00:00 2001 From: Bart Date: Mon, 17 Aug 2026 12:20:30 +0200 Subject: [PATCH 7/8] Leave the shrink-wrap rewrite out of this branch The rolling-ball rewrite rode along in the section-order commit and is unrelated to it. It also regresses the XFoil path: the wrapped NACA0012 of test/solver/test_backend_comparison.jl no longer converges at any angle, so generate_aero_matrices has nothing to interpolate the NaNs from and throws. The distance-field wrap comes back for now. The rewrite continues on wrap/rolling-ball, together with the two resampling fixes it needs and a record of what still fails. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 11 - docs/src/private_functions.md | 12 +- examples/V3_neuralfoil.jl | 4 +- examples/ram_air_kite.jl | 4 +- src/airfoil_aero/shrink_wrap.jl | 373 ++++++++++++++++++++------------ 5 files changed, 244 insertions(+), 160 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 55377a57..646c9351 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -69,17 +69,6 @@ already stores, take the core-radius cutoff from `|r1.r0|/|r0|` without forming the perpendicular vector, and defer the cross products that only one branch reads. Output is bit-identical; `solve!` is a further 1.47-1.55x faster. -- `shrink_wrap` traces the rolling ball exactly — pivoting it around the cloud and - emitting the arcs its contact side sweeps (`pivot_contour`) — instead of thresholding - and marching-squares-tracing a distance field, so there is no grid resolution left to - set. `ShrinkWrap`'s `cell_size` is accordingly named `min_clearance`, still accepted - under the old name, and only floors `clearance` — which is what it already did for a - single-membrane slice. The wrap now - sits at exactly `clearance` from the cloud instead of a cell over it, so a V3 canopy's - aft strip comes out `2 * clearance` thick where the grid gave `3.5 * cell_size`. - `min_concave_radius` also stops costing anything, having padded the grid in both - directions before: one V3 slice at radius 0.4 drops from 173 ms to 10 ms, and at the - default radius from 15 ms to 3 ms. - `read_node_table` parses into a preallocated matrix instead of `reduce(vcat, …)` over a generator, which was quadratic in the row count: ~21× faster on a 16 MB surface table (2.49 s → 0.12 s), benefiting every existing dataset. diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 9ea565b8..662d5fa0 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -115,13 +115,13 @@ get_lower_upper turn_trailing_edge! ``` -### Shrink-wrap rolling ball +### Shrink-wrap distance field ```@docs -point_buckets -turn_measure -pivot_step -push_arc! -pivot_contour +distance_parabolas! +squared_distance_transform! +grid_sampler +flood_outside +trace_level_set largest_linking_gap resample_arc smooth_turning! diff --git a/examples/V3_neuralfoil.jl b/examples/V3_neuralfoil.jl index 88064a5b..2acedbca 100644 --- a/examples/V3_neuralfoil.jl +++ b/examples/V3_neuralfoil.jl @@ -46,8 +46,8 @@ ROTATION = I # Marched leading-edge stations across the span (finer => smoother edge trace). N_BINS = 100 -# Shrink-wrap each raw slice into a clean closed airfoil: the rolling-ball wrap -# hugs the cloud at `clearance` (floored at `min_clearance`); a single-skin +# Shrink-wrap each raw slice into a clean closed airfoil: the distance-field wrap +# hugs the cloud at `clearance` (floored at one grid `cell_size`); a single-skin # canopy becomes a thin capsule. `min_concave_radius` sets the fillet bridging the # concave tube-canopy junction: 0.2 smooths the neck over entirely; the default # (0.02) traces it, which NeuralFoil also handles fine. diff --git a/examples/ram_air_kite.jl b/examples/ram_air_kite.jl index 2eabe7fe..fe2705c5 100644 --- a/examples/ram_air_kite.jl +++ b/examples/ram_air_kite.jl @@ -47,8 +47,8 @@ N_SECTIONS = 10 # At 0.0 XFoil only converges above 16 deg on the near-tip slice this mesh still has. WINGTIP_DISTANCE = 0.1 -# Shrink-wrap each raw slice into a clean closed airfoil: the rolling-ball wrap -# hugs the cloud at `clearance` (floored at `min_clearance`) with a round nose. +# Shrink-wrap each raw slice into a clean closed airfoil: the distance-field wrap +# hugs the cloud at `clearance` (floored at one grid `cell_size`) with a round nose. WRAP = ShrinkWrap(clearance=0.0) # 3D slice diagnostic (live preview): mesh + LE/TE curves + section contours, their diff --git a/src/airfoil_aero/shrink_wrap.jl b/src/airfoil_aero/shrink_wrap.jl index 45b8bdac..944d9a9b 100644 --- a/src/airfoil_aero/shrink_wrap.jl +++ b/src/airfoil_aero/shrink_wrap.jl @@ -1,176 +1,212 @@ """ - ShrinkWrap(; clearance=0.006, min_concave_radius=0.02, min_clearance=0.001, + ShrinkWrap(; clearance=0.006, min_concave_radius=0.02, cell_size=0.001, n_points=120, curvature_weight=0.05) -Rolling-ball shrink wrap: the offset of a raw slice point cloud, traced exactly as -one closed airfoil contour. A disk of radius `min_concave_radius` is pivoted around -the outside of the cloud ([`pivot_contour`](@ref)); the wrap is what its contact -side sweeps, so it consists of arcs of radius `clearance` about the points the disk -touches, joined by arcs of the disk radius across the gaps it cannot enter. Being -built from arcs rather than sampled on a grid, the wrap has no resolution floor: -the leading edge comes out genuinely round — the offset of the cloud nose — the -blunt trailing edge is capped by an arc of radius `clearance`, and a -single-membrane cloud becomes a capsule of exactly that half-thickness. +Distance-field shrink wrap: the rolling-ball offset of a raw slice point cloud, +extracted as one closed airfoil contour. A grid distance field is thresholded at +the rolling-ball radius (bridging cloud gaps and crevices narrower than about +twice `min_concave_radius`), flood-filled, eroded back to `clearance`, and the +resulting level set is traced with marching squares, faired (each pass clamped so +the contour keeps `clearance`) and resampled. Because the contour is +parameterized by arclength rather than `x`, the leading edge comes out genuinely +round — the offset of the cloud nose — the blunt trailing edge is capped by an arc +of radius `clearance`, and a single-membrane cloud becomes a thin capsule (at least +one `cell_size` half-thickness). # Fields -- `clearance`: offset the contour keeps outside every cloud point, and the radius - every convex corner is rounded at; floored at `min_clearance`. +- `clearance`: offset the contour keeps outside every cloud point; floored at one + `cell_size` (the grid cannot represent a tighter wrap). Also the radius every + convex corner is rounded at. - `min_concave_radius`: rolling-ball radius — concave features narrower than about twice this are bridged by a fillet of roughly this radius; convex geometry is unaffected. Auto-raised so the ball can neither fall through the cloud's largest - point gap nor dip below `clearance` between neighbouring points, so sparse clouds - get a correspondingly looser, smoother wrap. -- `min_clearance`: floor on `clearance`, so `clearance=0` still leaves a single - membrane a capsule with two distinguishable sides instead of a bare curve. Also - accepted under its former name `cell_size`, which set the resolution of the - distance field this used to be traced on. + point gap nor pinch off between points during erosion, so sparse clouds get a + correspondingly looser, smoother wrap. +- `cell_size`: distance-field grid resolution as a chord fraction; sets the + geometric fidelity of the wrap. - `n_points`: output stations per surface; the contour has `2*n_points - 1` points, cosine-clustered in arclength at the leading and trailing edges. - `curvature_weight`: extra sampling measure per radian of contour turning (chord fraction), concentrating output points into corners so XFoil's spline can follow them; `0` gives plain cosine-in-arclength sampling. """ -struct ShrinkWrap - clearance::Float64 - min_concave_radius::Float64 - min_clearance::Float64 - n_points::Int - curvature_weight::Float64 -end - -function ShrinkWrap(; clearance=0.006, min_concave_radius=0.02, min_clearance=0.001, - cell_size=nothing, n_points=120, curvature_weight=0.05) - return ShrinkWrap(clearance, min_concave_radius, - isnothing(cell_size) ? min_clearance : cell_size, - n_points, curvature_weight) +@with_kw struct ShrinkWrap + clearance::Float64 = 0.006 + min_concave_radius::Float64 = 0.02 + cell_size::Float64 = 0.001 + n_points::Int = 120 + curvature_weight::Float64 = 0.05 end """ - point_buckets(x, y, reach) -> Dict + distance_parabolas!(d, f, n, v, z) -Sparse uniform grid of the point indices, one bucket per `reach`-sized cell, so every -point within `reach` of a query lies in one of the nine buckets around it. +One pass of the Felzenszwalb–Huttenlocher distance transform: writes into `d` the +lower envelope of the parabolas `(q - p)^2 + f[p]` over `p`, for `q in 1:n`. +`v` and `z` are scratch (length `n` and `n + 1`). """ -function point_buckets(x, y, reach) - buckets = Dict{NTuple{2,Int},Vector{Int}}() - for k in eachindex(x) - key = (floor(Int, x[k] / reach), floor(Int, y[k] / reach)) - push!(get!(buckets, key, Int[]), k) +function distance_parabolas!(d, f, n, v, z) + k = 1 + v[1] = 1 + z[1] = -Inf + z[2] = Inf + for q in 2:n + s = ((f[q] + q * q) - (f[v[k]] + v[k] * v[k])) / (2 * (q - v[k])) + while s <= z[k] + k -= 1 + s = ((f[q] + q * q) - (f[v[k]] + v[k] * v[k])) / (2 * (q - v[k])) + end + k += 1 + v[k] = q + z[k] = s + z[k+1] = Inf + end + k = 1 + for q in 1:n + while z[k+1] < q + k += 1 + end + d[q] = (q - v[k])^2 + f[v[k]] end - return buckets + return d end """ - turn_measure(cross, dot) -> Float64 + squared_distance_transform!(field) -> field -Pseudo-angle in `[0, 4)` of the direction with the given `cross` and `dot` against a -reference, rising monotonically with the counterclockwise angle it stands for. Ordering -directions by it is ordering them by angle, without the `atan` that dominates the -pivot's inner loop. +In-place 2D squared Euclidean distance transform in grid-index units. On input +`field` holds `0` (or a local squared offset) at seeds and a large finite value +elsewhere; on output every node holds its squared distance to the nearest seed. """ -function turn_measure(cross, dot) - if cross >= 0 - return dot >= 0 ? cross / (cross + dot) : 1 + (-dot) / (cross - dot) +function squared_distance_transform!(field::Matrix{Float64}) + nx, ny = size(field) + m = max(nx, ny) + d, f = Vector{Float64}(undef, m), Vector{Float64}(undef, m) + v, z = Vector{Int}(undef, m), Vector{Float64}(undef, m + 1) + for i in 1:nx + f[1:ny] .= @view field[i, :] + distance_parabolas!(d, f, ny, v, z) + field[i, :] .= @view d[1:ny] + end + for j in 1:ny + f[1:nx] .= @view field[:, j] + distance_parabolas!(d, f, nx, v, z) + field[:, j] .= @view d[1:nx] end - return dot <= 0 ? 2 + (-cross) / (-cross - dot) : 3 + dot / (dot - cross) + return field end """ - pivot_step(x, y, buckets, r, p, centre) -> (q, centre) or nothing + flood_outside(blocked) -> BitMatrix -Rotate the empty disk of radius `r` counterclockwise about point `p`, starting from -`centre`, until its rim reaches a second point. Returns that point and the disk -centre touching both, or `nothing` when nothing lies within `2r` of `p`. +Mark the nodes reachable from the grid border without entering `blocked` +(4-connected flood fill); unreached nodes are the solid plus its enclosed holes. """ -function pivot_step(x, y, buckets, r, p, centre) - px, py = x[p], y[p] - ux, uy = centre[1] - px, centre[2] - py - best_turn, best_q, best_centre = Inf, 0, centre - i0, j0 = floor(Int, px / 2r), floor(Int, py / 2r) - for i in i0-1:i0+1, j in j0-1:j0+1 - bucket = get(buckets, (i, j), nothing) - bucket === nothing && continue - for q in bucket - q == p && continue - dx, dy = x[q] - px, y[q] - py - span = dx * dx + dy * dy - (span < 1e-24 || span > 4r * r) && continue - out = sqrt(max(r * r / span - 0.25, 0.0)) - mx, my = dx / 2, dy / 2 - ox, oy = -dy * out, dx * out - for (vx, vy) in ((mx + ox, my + oy), (mx - ox, my - oy)) - (vx - ux)^2 + (vy - uy)^2 < 1e-24 && continue - turn = turn_measure(ux * vy - uy * vx, ux * vx + uy * vy) - (turn < 1e-9 || turn >= best_turn) && continue - best_turn, best_q, best_centre = turn, q, (px + vx, py + vy) - end - end +function flood_outside(blocked::AbstractMatrix{Bool}) + nx, ny = size(blocked) + outside = falses(nx, ny) + stack = Int[] + visit(i, j) = if !blocked[i, j] && !outside[i, j] + outside[i, j] = true + push!(stack, (j - 1) * nx + i) end - best_q == 0 && return nothing - return best_q, best_centre + for i in 1:nx + visit(i, 1) + visit(i, ny) + end + for j in 1:ny + visit(1, j) + visit(nx, j) + end + while !isempty(stack) + idx = pop!(stack) + i, j = (idx - 1) % nx + 1, (idx - 1) ÷ nx + 1 + i > 1 && visit(i - 1, j) + i < nx && visit(i + 1, j) + j > 1 && visit(i, j - 1) + j < ny && visit(i, j + 1) + end + return outside end """ - push_arc!(px, py, cx, cy, radius, from, sweep) + trace_level_set(field, level) -> Vector{Vector{NTuple{2,Float64}}} -Append the arc about `(cx, cy)` that starts at angle `from` and turns by the signed -angle `sweep`, stepped fine enough to stay within `1e-6` of the true arc. The closing -point is left off so consecutive arcs concatenate without duplicates. +Marching-squares contours of `field .== level`, each returned as a closed loop of +grid-frame vertices (unit = one cell, node 1 at 1.0). Saddle cells are resolved by +the cell-center average. Loops touching the grid border are dropped. """ -function push_arc!(px, py, cx, cy, radius, from, sweep) - steps = radius < 1e-12 ? 1 : - clamp(ceil(Int, abs(sweep) / sqrt(8e-6 / radius)), 1, 400) - for k in 0:steps-1 - angle = from + sweep * k / steps - push!(px, cx + radius * cos(angle)) - push!(py, cy + radius * sin(angle)) +function trace_level_set(field::Matrix{Float64}, level::Float64) + nx, ny = size(field) + inside(v) = v > level + frac(a, b) = (level - a) / (b - a) + points = Dict{NTuple{3,Int},NTuple{2,Float64}}() + links = Dict{NTuple{3,Int},Vector{NTuple{3,Int}}}() + connect(a, b) = begin + push!(get!(links, a, NTuple{3,Int}[]), b) + push!(get!(links, b, NTuple{3,Int}[]), a) end - return nothing + for j in 1:ny-1, i in 1:nx-1 + v00, v10 = field[i, j], field[i+1, j] + v01, v11 = field[i, j+1], field[i+1, j+1] + b00, b10, b01, b11 = inside(v00), inside(v10), inside(v01), inside(v11) + b00 == b10 == b01 == b11 && continue + bottom, top = (i, j, 0), (i, j + 1, 0) + left, right = (i, j, 1), (i + 1, j, 1) + b00 != b10 && (points[bottom] = (i + frac(v00, v10), Float64(j))) + b01 != b11 && (points[top] = (i + frac(v01, v11), Float64(j + 1))) + b00 != b01 && (points[left] = (Float64(i), j + frac(v00, v01))) + b10 != b11 && (points[right] = (Float64(i + 1), j + frac(v10, v11))) + crossed = [e for (e, c) in ((bottom, b00 != b10), (top, b01 != b11), + (left, b00 != b01), (right, b10 != b11)) if c] + if length(crossed) == 4 + if inside((v00 + v10 + v01 + v11) / 4) == b00 + connect(bottom, right) + connect(top, left) + else + connect(bottom, left) + connect(top, right) + end + else + connect(crossed[1], crossed[2]) + end + end + visited = Set{NTuple{3,Int}}() + loops = Vector{Vector{NTuple{2,Float64}}}() + for start in keys(links) + start in visited && continue + loop = NTuple{2,Float64}[] + prev, current = start, start + closed = false + while true + push!(visited, current) + push!(loop, points[current]) + nbrs = links[current] + length(nbrs) == 2 || break + prev, current = current, (nbrs[1] == prev ? nbrs[2] : nbrs[1]) + current == start && (closed = true; break) + end + closed && push!(loops, loop) + end + return loops end """ - pivot_contour(x, y, r, clearance) -> (px, py) + grid_sampler(field, x0, y0, cell) -> f(px, py) -Trace the rolling-ball closing of the cloud `(x, y)` offset outward by `clearance`, -as a dense closed polyline. A disk of radius `r` is pivoted around the outside of the -cloud from its leftmost point; each point the disk touches contributes an arc of -radius `clearance` about itself, and each pivot an arc of radius `r - clearance` -about the empty disk's centre, spanning the gap the disk could not enter. Those arcs -are the exact wrap boundary, so nothing here has a resolution to choose. +Bilinear interpolant of `field` at physical points; node `(i, j)` sits at +`(x0 + (i-1)*cell, y0 + (j-1)*cell)`. Queries are clamped to the grid. """ -function pivot_contour(x, y, r, clearance) - start = argmin(x) - buckets = point_buckets(x, y, 2r) - contacts, centres = Int[], NTuple{2,Float64}[] - p, centre = start, (x[start] - r, y[start]) - # A single-membrane stretch is touched once from each side, so the cycle closes on - # the (point, centre) pair rather than on the point alone. - for _ in 1:4*length(x)+100 - step = pivot_step(x, y, buckets, r, p, centre) - step === nothing && break - next, centre = step - !isempty(contacts) && p == contacts[1] && - hypot(centre[1] - centres[1][1], centre[2] - centres[1][2]) < 1e-12 && - break - push!(contacts, p) - push!(centres, centre) - p = next - end - px, py = Float64[], Float64[] - isempty(contacts) && return px, py - for (k, i) in enumerate(contacts) - incoming = centres[mod1(k - 1, length(centres))] - outgoing = centres[k] - from = atan(incoming[2] - y[i], incoming[1] - x[i]) - to = atan(outgoing[2] - y[i], outgoing[1] - x[i]) - push_arc!(px, py, x[i], y[i], clearance, from, mod(to - from, 2pi)) - j = contacts[mod1(k + 1, length(contacts))] - from = atan(y[i] - outgoing[2], x[i] - outgoing[1]) - to = atan(y[j] - outgoing[2], x[j] - outgoing[1]) - push_arc!(px, py, outgoing[1], outgoing[2], r - clearance, from, - rem(to - from, 2pi, RoundNearest)) +function grid_sampler(field::Matrix{Float64}, x0, y0, cell) + nx, ny = size(field) + return function (px, py) + gi = clamp((px - x0) / cell + 1, 1.0, nx - 1e-9) + gj = clamp((py - y0) / cell + 1, 1.0, ny - 1e-9) + i, j = floor(Int, gi), floor(Int, gj) + ti, tj = gi - i, gj - j + return (field[i, j] * (1 - ti) + field[i+1, j] * ti) * (1 - tj) + + (field[i, j+1] * (1 - ti) + field[i+1, j+1] * ti) * tj end - return px, py end """ @@ -267,11 +303,11 @@ end shrink_wrap(x, y, method::ShrinkWrap) -> (x, y) Wrap the point cloud `(x, y)` into a clean closed airfoil in Selig order (TE upper → -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 +LE → TE lower), following [`ShrinkWrap`](@ref): distance field on a `cell_size` +grid, closing with the rolling ball (`min_concave_radius`), offset outward by +`clearance`, traced as a single +closed contour and 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). 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 @@ -279,14 +315,73 @@ normalized frame of the input cloud (chord slightly longer than 1, nose apex nea """ function shrink_wrap(x, y, method::ShrinkWrap) xn, yn, _ = normalize_airfoil(collect(float.(x)), collect(float.(y))) - gap = max(method.clearance, method.min_clearance) + cell = method.cell_size + gap = max(method.clearance, cell) linking = largest_linking_gap(xn, yn) - ball = max(method.min_concave_radius, gap, 1.01 * linking / 2, + ball = max(method.min_concave_radius, gap + 3cell, 1.5 * (linking^2 / 4 + gap^2) / (2gap)) - px, py = pivot_contour(xn, yn, ball, gap) - length(px) < 3 && error("ShrinkWrap traced no contour; the cloud may be too" * - " sparse or thin to pivot a ball of $(ball) around.") + pad = ball + 3cell + x0, y0 = -pad, minimum(yn) - pad + nx = ceil(Int, (1.0 + pad - x0) / cell) + 1 + ny = ceil(Int, (maximum(yn) + pad - y0) / cell) + 1 + + far = Float64(nx^2 + ny^2) + field = fill(far, nx, ny) + for (px, py) in zip(xn, yn) + gi, gj = (px - x0) / cell + 1, (py - y0) / cell + 1 + for i in (floor(Int, gi), floor(Int, gi) + 1), + j in (floor(Int, gj), floor(Int, gj) + 1) + + (1 <= i <= nx && 1 <= j <= ny) || continue + field[i, j] = min(field[i, j], (gi - i)^2 + (gj - j)^2) + end + end + squared_distance_transform!(field) + cloud_dist = sqrt.(field) .* cell + + outside = flood_outside(cloud_dist .<= ball) + depth = [outside[i, j] ? 0.0 : far for i in 1:nx, j in 1:ny] + squared_distance_transform!(depth) + depth .= sqrt.(depth) .* cell + + loops = trace_level_set(depth, ball - gap) + isempty(loops) && error("ShrinkWrap traced no contour; the cloud may be too" * + " sparse or thin for cell_size=$(cell).") + shoelace(loop) = abs(sum(p[1] * loop[mod1(k + 1, length(loop))][2] - + loop[mod1(k + 1, length(loop))][1] * p[2] + for (k, p) in enumerate(loop))) / 2 + areas = shoelace.(loops) + main = loops[argmax(areas)] + if length(loops) > 1 && sort(areas)[end-1] > 0.01 * maximum(areas) + @warn "ShrinkWrap discarded a comparable secondary contour; increase" * + " min_concave_radius to bridge cloud gaps." + end + px = [x0 + (p[1] - 1) * cell for p in main] + py = [y0 + (p[2] - 1) * cell for p in main] + + dist_at = grid_sampler(cloud_dist, x0, y0, cell) + keep_clear(cx, cy) = begin + d = dist_at(cx, cy) + d >= gap && return (cx, cy) + h = cell / 2 + gx = (dist_at(cx + h, cy) - dist_at(cx - h, cy)) / cell + gy = (dist_at(cx, cy + h) - dist_at(cx, cy - h)) / cell + len = hypot(gx, gy) + len < 1e-9 && return (cx, cy) + return (cx + (gap - d) * gx / len, cy + (gap - d) * gy / len) + end m = length(px) + for pass in 0:20 + ox, oy = copy(px), copy(py) + for k in 1:m + if pass > 0 + a, b = mod1(k - 1, m), mod1(k + 1, m) + px[k] = ox[k] + 0.4 * (ox[a] - 2ox[k] + ox[b]) + py[k] = oy[k] + 0.4 * (oy[a] - 2oy[k] + oy[b]) + end + px[k], py[k] = keep_clear(px[k], py[k]) + end + end anchor = argmax(px) px, py = circshift(px, 1 - anchor), circshift(py, 1 - anchor) From 3d627b1620e480a8e01a33e22a15b1dc2d848e7b Mon Sep 17 00:00:00 2001 From: Bart Date: Mon, 17 Aug 2026 12:31:56 +0200 Subject: [PATCH 8/8] Leave the Section method of reinit! undocumented Its positional defaults make the docstring's signature a Union of four arities, which no call-form @docs entry can match, so Documenter counted the docstring as missing however it was listed. The text was an orphan attached to the sort-key helper before this branch, and the sibling refined-section overload carries none either. Verified with a local docs build. Co-Authored-By: Claude Opus 5 (1M context) --- docs/src/private_functions.md | 1 - src/wing_geometry.jl | 6 ------ 2 files changed, 7 deletions(-) diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 662d5fa0..8841f5e7 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -45,7 +45,6 @@ update_panel_properties! build_interps panel_interp_types reinit!(wing::AbstractWing) -reinit!(section::Section, LE_point, TE_point, aero_model, aero_data, section_aero) reinit!(panel::Panel, section_1::Section, section_2::Section, aero_center, control_point, bound_point_1, bound_point_2, x_airf, y_airf, z_airf, delta, vec) rotated_te calculate_filaments_for_plotting diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index 247e5db6..be6d5ebd 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -69,12 +69,6 @@ function normalize_span_order!(sections) return sections end -""" - reinit!(section::Section, LE_point, TE_point, aero_model=nothing, aero_data=nothing, - section_aero=nothing) - -Update a [Section](@ref) in place. -""" function reinit!(section::Section, LE_point, TE_point, aero_model=nothing, aero_data=nothing, section_aero=nothing) section.LE_point .= LE_point