From f078c77e32118f36476cd565b5513ac031f13f15 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sat, 25 Oct 2025 22:36:47 +0200 Subject: [PATCH 01/22] Make sure 0 groups works --- src/ram_geometry.jl | 2 +- src/solver.jl | 53 +++++-- src/wing_geometry.jl | 2 +- src/yaml_geometry.jl | 2 +- test/data/solver/wings/solver_test_wing.yaml | 12 ++ test/solver/test_solver.jl | 148 +++++++++++++++++++ test/yaml_geometry/test_wing_constructor.jl | 7 +- 7 files changed, 206 insertions(+), 20 deletions(-) diff --git a/src/ram_geometry.jl b/src/ram_geometry.jl index e352c4e2..c267d727 100644 --- a/src/ram_geometry.jl +++ b/src/ram_geometry.jl @@ -415,7 +415,7 @@ function RamAirWing( interp_steps=n_sections # TODO: check if interpolations are still needed ) - !(n_panels % n_groups == 0) && throw(ArgumentError("Number of panels should be divisible by number of groups")) + !(n_groups == 0 || n_panels % n_groups == 0) && throw(ArgumentError("Number of panels should be divisible by number of groups")) !isapprox(spanwise_direction, [0.0, 1.0, 0.0]) && throw(ArgumentError("Spanwise direction has to be [0.0, 1.0, 0.0], not $spanwise_direction")) # Load or create polars diff --git a/src/solver.jl b/src/solver.jl index ecf34e3b..8b2d75fd 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -292,21 +292,29 @@ function solve!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution= moment_coeff_dist[i] = moment_dist[i] / (q_inf * projected_area) end - group_moment_dist = solver.sol.group_moment_dist - group_moment_coeff_dist = solver.sol.group_moment_coeff_dist - group_moment_dist .= 0.0 - group_moment_coeff_dist .= 0.0 - panel_idx = 1 - group_idx = 1 - for wing in body_aero.wings - panels_per_group = wing.n_panels ÷ wing.n_groups - for _ in 1:wing.n_groups - for _ in 1:panels_per_group - group_moment_dist[group_idx] += moment_dist[panel_idx] - group_moment_coeff_dist[group_idx] += moment_coeff_dist[panel_idx] - panel_idx += 1 + # Only compute group moments if there are groups + if length(solver.sol.group_moment_dist) > 0 + group_moment_dist = solver.sol.group_moment_dist + group_moment_coeff_dist = solver.sol.group_moment_coeff_dist + group_moment_dist .= 0.0 + group_moment_coeff_dist .= 0.0 + panel_idx = 1 + group_idx = 1 + for wing in body_aero.wings + if wing.n_groups > 0 + panels_per_group = wing.n_panels ÷ wing.n_groups + for _ in 1:wing.n_groups + for _ in 1:panels_per_group + group_moment_dist[group_idx] += moment_dist[panel_idx] + group_moment_coeff_dist[group_idx] += moment_coeff_dist[panel_idx] + panel_idx += 1 + end + group_idx += 1 + end + else + # Skip panels for wings with n_groups=0 + panel_idx += wing.n_panels end - group_idx += 1 end end @@ -689,8 +697,8 @@ jac, results = linearize( ) ``` """ -function linearize(solver::Solver, body_aero::BodyAerodynamics, y::Vector{T}; - theta_idxs=1:4, +function linearize(solver::Solver, body_aero::BodyAerodynamics, y::Vector{T}; + theta_idxs=1:4, delta_idxs=nothing, va_idxs=nothing, omega_idxs=nothing, @@ -700,6 +708,19 @@ function linearize(solver::Solver, body_aero::BodyAerodynamics, y::Vector{T}; !(length(body_aero.wings) == 1) && throw(ArgumentError("Linearization only works for a body_aero with one wing")) wing = body_aero.wings[1] + # Validate that theta_idxs and delta_idxs match the number of groups + if !isnothing(theta_idxs) && wing.n_groups > 0 + length(theta_idxs) != wing.n_groups && throw(ArgumentError( + "Length of theta_idxs ($(length(theta_idxs))) must match number of groups ($(wing.n_groups))")) + end + if !isnothing(delta_idxs) && wing.n_groups > 0 + length(delta_idxs) != wing.n_groups && throw(ArgumentError( + "Length of delta_idxs ($(length(delta_idxs))) must match number of groups ($(wing.n_groups))")) + end + if wing.n_groups == 0 && (!isnothing(theta_idxs) || !isnothing(delta_idxs)) + throw(ArgumentError("Cannot use theta_idxs or delta_idxs when wing has n_groups=0 (no group functionality)")) + end + init_va = body_aero.cache[1][body_aero.va] init_va .= body_aero.va if !isnothing(theta_idxs) diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index 506b91b4..bb1f5e06 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -235,7 +235,7 @@ function Wing(n_panels::Int; spanwise_distribution::PanelDistribution=LINEAR, spanwise_direction::PosVector=MVec3([0.0, 1.0, 0.0]), remove_nan=true) - !(n_panels % n_groups == 0) && throw(ArgumentError("Number of panels should be divisible by number of groups")) + !(n_groups == 0 || n_panels % n_groups == 0) && throw(ArgumentError("Number of panels should be divisible by number of groups")) panel_props = PanelProperties{n_panels}() Wing(n_panels, n_groups, spanwise_distribution, panel_props, spanwise_direction, Section[], Section[], remove_nan) end diff --git a/src/yaml_geometry.jl b/src/yaml_geometry.jl index 1abf3118..3ba618f5 100644 --- a/src/yaml_geometry.jl +++ b/src/yaml_geometry.jl @@ -190,7 +190,7 @@ function Wing( remove_nan=true, prn=false ) - !(n_panels % n_groups == 0) && throw(ArgumentError("Number of panels should be divisible by number of groups")) + !(n_groups == 0 || n_panels % n_groups == 0) && throw(ArgumentError("Number of panels should be divisible by number of groups")) !isapprox(spanwise_direction, [0.0, 1.0, 0.0]) && throw(ArgumentError("Spanwise direction has to be [0.0, 1.0, 0.0], not $spanwise_direction")) prn && @info "Reading YAML wing configuration from $geometry_file" diff --git a/test/data/solver/wings/solver_test_wing.yaml b/test/data/solver/wings/solver_test_wing.yaml index e69de29b..14969bce 100644 --- a/test/data/solver/wings/solver_test_wing.yaml +++ b/test/data/solver/wings/solver_test_wing.yaml @@ -0,0 +1,12 @@ +wing_sections: + headers: [airfoil_id, LE_x, LE_y, LE_z, TE_x, TE_y, TE_z] + data: + - [1, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0] + - [1, 0.0, -1.0, 0.0, 1.0, -1.0, 0.0] + +wing_airfoils: + alpha_range: [-10, 10, 1] + reynolds: 1000000 + headers: [airfoil_id, type, info_dict] + data: + - [1, polars, {csv_file_path: ""}] diff --git a/test/solver/test_solver.jl b/test/solver/test_solver.jl index 1d60061a..0dbb9971 100644 --- a/test/solver/test_solver.jl +++ b/test/solver/test_solver.jl @@ -29,4 +29,152 @@ using Test rm(settings_file; force=true) end end + + @testset "Solver with n_groups=0" begin + # Test that solver works correctly when n_groups=0 (no group functionality) + settings_file = create_temp_wing_settings("solver", "solver_test_wing.yaml"; + alpha=5.0, beta=0.0, wind_speed=10.0, n_groups=0) + + try + settings = VSMSettings(settings_file) + wing = Wing(settings) + @test wing.n_groups == 0 + + body_aero = BodyAerodynamics([wing]) + solver = Solver(body_aero, settings) + + # Verify solver has zero groups + @test length(solver.sol.group_moment_dist) == 0 + @test length(solver.sol.group_moment_coeff_dist) == 0 + + # Test that solve! works without errors + va = [10.0, 0.0, 0.0] + set_va!(body_aero, va) + sol = solve!(solver, body_aero) + + @test sol isa VSMSolution + @test sol.solver_status == SUCCESS + + # Verify that group moments are empty + @test length(sol.group_moment_dist) == 0 + @test length(sol.group_moment_coeff_dist) == 0 + + # But force and moment should still be computed + @test !all(sol.force .== 0.0) + @test norm(sol.force) > 0 + + finally + rm(settings_file; force=true) + end + end + + @testset "Linearize with n_groups=0" begin + # Test that linearize works correctly when n_groups=0 + settings_file = create_temp_wing_settings("solver", "solver_test_wing.yaml"; + alpha=5.0, beta=0.0, wind_speed=10.0, n_groups=0, n_panels=4) + + try + settings = VSMSettings(settings_file) + wing = Wing(settings) + @test wing.n_groups == 0 + + body_aero = BodyAerodynamics([wing]) + solver = Solver(body_aero, settings) + + # Set velocity + va = [10.0, 0.0, 0.0] + set_va!(body_aero, va) + + # Test linearize with only velocity (no theta or delta since n_groups=0) + y = va # Only velocity in input vector + jac, results = VortexStepMethod.linearize( + solver, + body_aero, + y; + theta_idxs=nothing, + delta_idxs=nothing, + va_idxs=1:3, + omega_idxs=nothing + ) + + # Results should only have 6 elements (force + moment, no group moments) + @test length(results) == 6 + @test size(jac) == (6, 3) # 6 outputs, 3 inputs (vx, vy, vz) + + # Verify forces are non-zero + @test norm(results[1:3]) > 0 + + # Test that using theta_idxs with n_groups=0 throws an error + @test_throws ArgumentError VortexStepMethod.linearize( + solver, + body_aero, + [0.0, 10.0, 0.0, 0.0]; # Invalid: trying to use theta + theta_idxs=1:1, + va_idxs=2:4 + ) + + # Test that using delta_idxs with n_groups=0 throws an error + @test_throws ArgumentError VortexStepMethod.linearize( + solver, + body_aero, + [0.0, 10.0, 0.0, 0.0]; # Invalid: trying to use delta + delta_idxs=1:1, + va_idxs=2:4 + ) + + finally + rm(settings_file; force=true) + end + end + + @testset "Linearize theta_idxs validation" begin + # Test that theta_idxs length must match n_groups + settings_file = create_temp_wing_settings("solver", "solver_test_wing.yaml"; + alpha=5.0, beta=0.0, wind_speed=10.0, n_groups=2, n_panels=4) + + try + settings = VSMSettings(settings_file) + wing = Wing(settings) + @test wing.n_groups == 2 + + body_aero = BodyAerodynamics([wing]) + solver = Solver(body_aero, settings) + + va = [10.0, 0.0, 0.0] + set_va!(body_aero, va) + + # Test with correct number of theta angles (2) + y_correct = [0.0, 0.0, 10.0, 0.0, 0.0] # 2 theta + 3 va + jac, results = VortexStepMethod.linearize( + solver, + body_aero, + y_correct; + theta_idxs=1:2, + va_idxs=3:5 + ) + @test size(jac, 1) == 8 # 6 + 2 group moments + + # Test with wrong number of theta angles (should throw error) + y_wrong = [0.0, 0.0, 0.0, 0.0, 10.0, 0.0, 0.0] # 4 theta + 3 va + @test_throws ArgumentError VortexStepMethod.linearize( + solver, + body_aero, + y_wrong; + theta_idxs=1:4, # Wrong: 4 angles but only 2 groups + va_idxs=5:7 + ) + + # Test with wrong number of delta angles + @test_throws ArgumentError VortexStepMethod.linearize( + solver, + body_aero, + y_wrong; + delta_idxs=1:4, # Wrong: 4 angles but only 2 groups + va_idxs=5:7 + ) + + finally + rm(settings_file; force=true) + end + end end diff --git a/test/yaml_geometry/test_wing_constructor.jl b/test/yaml_geometry/test_wing_constructor.jl index b7b30a45..19ec9869 100644 --- a/test/yaml_geometry/test_wing_constructor.jl +++ b/test/yaml_geometry/test_wing_constructor.jl @@ -151,7 +151,12 @@ wing_airfoils: # Test invalid n_panels/n_groups combination @test_throws ArgumentError Wing(test_yaml_path; n_panels=5, n_groups=2) - + + # Test n_groups=0 (no grouping functionality) + wing_no_groups = Wing(test_yaml_path; n_panels=4, n_groups=0) + @test wing_no_groups.n_groups == 0 + @test wing_no_groups.n_panels == 4 + # Test invalid spanwise direction @test_throws ArgumentError Wing(test_yaml_path; spanwise_direction=[1.0, 0.0, 0.0]) end From dba02d946901d56eb4e3141ab75554099d15e9e8 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sun, 26 Oct 2025 10:45:22 +0100 Subject: [PATCH 02/22] Unify Wing and RamAirWing --- docs/src/types.md | 10 +- examples/bench.jl | 6 +- examples/ram_air_kite.jl | 8 +- src/VortexStepMethod.jl | 4 +- src/body_aerodynamics.jl | 6 +- src/{ram_geometry.jl => obj_geometry.jl} | 239 ++++------------------- src/solver.jl | 8 +- src/wing_geometry.jl | 197 ++++++++++++++++++- test/body_aerodynamics/test_results.jl | 2 +- test/plotting/test_plotting.jl | 2 +- test/ram_geometry/test_kite_geometry.jl | 12 +- test/solver/test_solver.jl | 45 ++++- test/wake/test_wake.jl | 2 +- 13 files changed, 309 insertions(+), 232 deletions(-) rename src/{ram_geometry.jl => obj_geometry.jl} (66%) diff --git a/docs/src/types.md b/docs/src/types.md index 14e5cbdd..de083b35 100644 --- a/docs/src/types.md +++ b/docs/src/types.md @@ -24,17 +24,17 @@ AeroData ``` ## Wing Geometry, Panel and Aerodynamics -A body is constructed of one or more abstract wings. An abstract wing can be a Wing or a RamAirWing. -A Wing/ RamAirWing has one or more sections. +A body is constructed of one or more abstract wings. All wings are of type Wing. +A Wing has one or more sections and can be created from YAML files or OBJ geometry. ```@docs Section Section(LE_point::PosVector, TE_point::PosVector, aero_model) Wing Wing(n_panels::Int; spanwise_distribution::PanelDistribution=LINEAR, spanwise_direction::PosVector=MVec3([0.0, 1.0, 0.0])) -RamAirWing -RamAirWing(obj_path, dat_path; alpha=0.0, crease_frac=0.75, wind_vel=10., mass=1.0, - n_panels=54, n_sections=n_panels+1, spanwise_distribution=UNCHANGED, +ObjWing +ObjWing(obj_path, dat_path; alpha=0.0, crease_frac=0.75, wind_vel=10., mass=1.0, + n_panels=54, n_sections=n_panels+1, spanwise_distribution=UNCHANGED, spanwise_direction=[0.0, 1.0, 0.0]) BodyAerodynamics ``` diff --git a/examples/bench.jl b/examples/bench.jl index 89544278..90c41ff1 100644 --- a/examples/bench.jl +++ b/examples/bench.jl @@ -56,9 +56,9 @@ println("Rectangular wing, solve:") @time solve(vsm_solver, body_aero, nothing) # Create wing geometry -wing = RamAirWing( - joinpath("data", "ram_air_kite", "ram_air_kite_body.obj"), - joinpath("data", "ram_air_kite", "ram_air_kite_foil.dat"); +wing = ObjWing( + joinpath("data", "ram_air_kite", "ram_air_kite_body.obj"), + joinpath("data", "ram_air_kite", "ram_air_kite_foil.dat"); prn=false ) body_aero = BodyAerodynamics([wing]) diff --git a/examples/ram_air_kite.jl b/examples/ram_air_kite.jl index f10b8d5b..13e60c62 100644 --- a/examples/ram_air_kite.jl +++ b/examples/ram_air_kite.jl @@ -9,9 +9,9 @@ DEFORM = false LINEARIZE = false # Create wing geometry -wing = RamAirWing( - joinpath("data", "ram_air_kite", "ram_air_kite_body.obj"), - joinpath("data", "ram_air_kite", "ram_air_kite_foil.dat"); +wing = ObjWing( + joinpath("data", "ram_air_kite", "ram_air_kite_body.obj"), + joinpath("data", "ram_air_kite", "ram_air_kite_foil.dat"); prn=PRN ) body_aero = BodyAerodynamics([wing];) @@ -21,7 +21,7 @@ println("First init") if DEFORM # Linear interpolation of alpha from 10° at one tip to 0° at the other println("Deform") - @time VortexStepMethod.smooth_deform!(wing, deg2rad.([10,20,10,0]), deg2rad.([-10,0,-10,0])) + @time group_deform!(wing, deg2rad.([10,20,10,0]), deg2rad.([-10,0,-10,0]); smooth=true) println("Deform init") @time VortexStepMethod.reinit!(body_aero; init_aero=false) end diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index c852c48e..64a4a657 100644 --- a/src/VortexStepMethod.jl +++ b/src/VortexStepMethod.jl @@ -27,7 +27,7 @@ using Xfoil # Export public interface export VSMSettings, WingSettings, SolverSettings -export Wing, Section, RamAirWing, reinit! +export Wing, Section, ObjWing, reinit! export BodyAerodynamics export Solver, solve, solve_base!, solve!, VSMSolution, linearize export calculate_results @@ -272,7 +272,7 @@ end include("settings.jl") include("wing_geometry.jl") include("polars.jl") -include("ram_geometry.jl") +include("obj_geometry.jl") include("yaml_geometry.jl") include("filament.jl") include("panel.jl") diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 87e75fec..6444a598 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -5,7 +5,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru # Fields - panels::Vector{Panel}: Vector of [Panel](@ref) structs -- wings::Union{Vector{Wing}, Vector{RamAirWing}}: A vector of wings; a body can have multiple wings +- wings::Vector{Wing}: A vector of wings; a body can have multiple wings - `va::MVec3` = zeros(MVec3): A vector of the apparent wind speed, see: [MVec3](@ref) - `omega`::MVec3 = zeros(MVec3): A vector of the turn rates around the kite body axes - `gamma_distribution`=zeros(Float64, P): A vector of the circulation @@ -23,7 +23,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru """ @with_kw mutable struct BodyAerodynamics{P} panels::Vector{Panel} - wings::Union{Vector{Wing}, Vector{RamAirWing}} + wings::Vector{Wing} _va::MVec3 = zeros(MVec3) omega::MVec3 = zeros(MVec3) gamma_distribution::MVector{P, Float64} = zeros(P) @@ -142,7 +142,7 @@ function reinit!(body_aero::BodyAerodynamics; # Create panels for i in 1:wing.n_panels - if wing isa RamAirWing + if !isnothing(wing.delta_dist) delta = wing.delta_dist[i] else delta = 0.0 diff --git a/src/ram_geometry.jl b/src/obj_geometry.jl similarity index 66% rename from src/ram_geometry.jl rename to src/obj_geometry.jl index c267d727..dec38cb7 100644 --- a/src/ram_geometry.jl +++ b/src/obj_geometry.jl @@ -14,7 +14,7 @@ Read vertices and faces from an OBJ file. function read_faces(filename) vertices = [] faces = [] - + open(filename) do file for line in eachline(file) if startswith(line, "v ") && !startswith(line, "vt") && !startswith(line, "vn") @@ -115,13 +115,13 @@ function create_interpolations(vertices, circle_center_z, radius, gamma_tip, R=I gamma_range = range(-gamma_tip+1e-6, gamma_tip-1e-6, interp_steps) stepsize = gamma_range.step.hi vz_centered = [v[3] - circle_center_z for v in vertices] - + te_gammas = zeros(length(gamma_range)) le_gammas = zeros(length(gamma_range)) trailing_edges = zeros(3, length(gamma_range)) leading_edges = zeros(3, length(gamma_range)) areas = zeros(length(gamma_range)) - + for (j, gamma) in enumerate(gamma_range) trailing_edges[1, j] = -Inf leading_edges[1, j] = Inf @@ -164,7 +164,7 @@ function create_interpolations(vertices, circle_center_z, radius, gamma_tip, R=I te_interp = ntuple(i -> linear_interpolation(le_gammas, trailing_edges[i, :], extrapolation_bc=Line()), 3) area_interp = linear_interpolation(gamma_range, areas, extrapolation_bc=Line()) - + return (le_interp, te_interp, area_interp) end @@ -187,26 +187,26 @@ Calculate center of mass of a mesh and translate vertices so that COM is at orig function center_to_com!(vertices, faces; prn=true) area_total = 0.0 com = zeros(3) - + for face in faces if length(face) == 3 # Triangle case v1 = vertices[face[1]] v2 = vertices[face[2]] v3 = vertices[face[3]] - + # Calculate triangle area and centroid normal = cross(v2 - v1, v3 - v1) area = norm(normal) / 2 centroid = (v1 + v2 + v3) / 3 - + area_total += area com -= area * centroid else throw(ArgumentError("Triangulate faces in a CAD program first")) end end - + com = com / area_total !(abs(com[2]) < 0.01) && throw(ArgumentError("Center of mass $com of .obj file has to lie on the xz-plane.")) prn && @info "Centering vertices of .obj file to the center of mass: $com" @@ -220,7 +220,7 @@ end """ calculate_inertia_tensor(vertices, faces, mass, com) -Calculate the inertia tensor for a triangulated surface mesh, assuming a thin shell with uniform +Calculate the inertia tensor for a triangulated surface mesh, assuming a thin shell with uniform surface density. # Arguments @@ -244,23 +244,23 @@ function calculate_inertia_tensor(vertices, faces, mass, com) # Initialize inertia tensor I = zeros(3, 3) total_area = 0.0 - + for face in faces v1 = vertices[face[1]] .- com v2 = vertices[face[2]] .- com v3 = vertices[face[3]] .- com - + # Calculate triangle area normal = cross(v2 - v1, v3 - v1) area = norm(normal) / 2 total_area += area - + # Calculate contribution to inertia tensor for i in 1:3 for j in 1:3 # Vertices relative to center of mass points = [v1, v2, v3] - + # Calculate contribution to inertia tensor for p in points if i == j @@ -274,7 +274,7 @@ function calculate_inertia_tensor(vertices, faces, mass, com) end end end - + # Scale by mass/total_area to get actual inertia tensor return (mass / total_area) * I / 3 end @@ -291,14 +291,14 @@ function calc_inertia_y_rotation(I_b_tensor) # Transform inertia tensor I_rotated = R_y * I_b_tensor * R_y' # We want the off-diagonal xz elements to be zero - F[1] = I_rotated[1,3] + F[1] = I_rotated[1,3] end - + theta0 = [0.0] prob = NonlinearProblem(eq!, theta0, nothing) sol = NonlinearSolve.solve(prob, NewtonRaphson()) theta_opt = sol.u[1] - + R_b_p = [ cos(theta_opt) 0 sin(theta_opt); 0 1 0; @@ -312,67 +312,18 @@ end """ - RamAirWing <: AbstractWing - -A ram-air wing model that represents a curved parafoil with deformable aerodynamic surfaces. - -## Core Features -- Curved wing geometry derived from 3D mesh (.obj file) -- Aerodynamic properties based on 2D airfoil data (.dat file) -- Support for control inputs (twist angles and trailing edge deflections) -- Inertial and geometric properties calculation - -## Notable Fields -- `n_panels::Int16`: Number of panels in aerodynamic mesh -- `n_groups::Int16`: Number of control groups for distributed deformation -- `mass::Float64`: Total wing mass in kg -- `gamma_tip::Float64`: Angular extent from center to wing tip -- `inertia_tensor::Matrix{Float64}`: Full 3x3 inertia tensor in the kite body frame -- `T_cad_body::MVec3`: Translation vector from CAD frame to body frame -- `R_cad_body::MMat3`: Rotation matrix from CAD frame to body frame -- `radius::Float64`: Wing curvature radius -- `theta_dist::Vector{Float64}`: Panel twist angle distribution -- `delta_dist::Vector{Float64}`: Trailing edge deflection distribution - -See constructor `RamAirWing(obj_path, dat_path; kwargs...)` for usage details. -""" -mutable struct RamAirWing <: AbstractWing - n_panels::Int16 - n_groups::Int16 - spanwise_distribution::PanelDistribution - panel_props::PanelProperties - spanwise_direction::MVec3 - sections::Vector{Section} - refined_sections::Vector{Section} - remove_nan::Bool - - # Additional fields for RamAirWing - non_deformed_sections::Vector{Section} - mass::Float64 - gamma_tip::Float64 - inertia_tensor::Matrix{Float64} - T_cad_body::MVec3 - R_cad_body::MMat3 - radius::Float64 - le_interp::NTuple{3, Extrapolation} - te_interp::NTuple{3, Extrapolation} - area_interp::Extrapolation - theta_dist::Vector{Float64} - delta_dist::Vector{Float64} - cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} -end + ObjWing(obj_path, dat_path; kwargs...) -""" - RamAirWing(obj_path, dat_path; kwargs...) - -Create a ram-air wing model from 3D geometry and airfoil data files. +Create a deformable wing model from 3D geometry (.obj) and airfoil data (.dat) files. This constructor builds a complete aerodynamic model by: -1. Loading or generating wing geometry from the .obj file -2. Creating aerodynamic polars from the airfoil .dat file +1. Loading wing geometry from the .obj file +2. Creating aerodynamic polars from the airfoil .dat file (or loading existing) 3. Computing inertial properties and coordinate transformations 4. Setting up control surfaces and panel distribution +The resulting Wing supports deformation through group_deform! and deform! functions. + # Arguments - `obj_path`: Path to .obj file containing 3D wing geometry - `dat_path`: Path to .dat file containing 2D airfoil profile @@ -389,32 +340,34 @@ This constructor builds a complete aerodynamic model by: - `remove_nan=true`: Interpolate NaN values in aerodynamic data - `alpha_range=deg2rad.(-5:1:20)`: Angle of attack range for polars (rad) - `delta_range=deg2rad.(-5:1:20)`: Trailing edge deflection range for polars (rad) -- prn=true: if info messages shall be printed +- `prn=true`: Print informational messages # Returns -A fully initialized `RamAirWing` instance ready for aerodynamic simulation. +A fully initialized `Wing` instance ready for aerodynamic simulation with deformation support. # Example ```julia -# Create a ram-air wing from geometry files -wing = RamAirWing( +# Create a deformable wing from geometry files +wing = ObjWing( "path/to/wing.obj", "path/to/airfoil.dat"; mass=1.5, n_panels=40, n_groups=4 ) + +# Apply deformation +group_deform!(wing, deg2rad.([5, 10, 5, 0]), deg2rad.([-5, 0, -5, 0])) ``` """ -function RamAirWing( - obj_path, dat_path; - crease_frac=0.9, wind_vel=10., mass=1.0, - n_panels=56, n_sections=n_panels+1, n_groups=4, spanwise_distribution=UNCHANGED, +function ObjWing( + obj_path, dat_path; + crease_frac=0.9, wind_vel=10., mass=1.0, + n_panels=56, n_sections=n_panels+1, n_groups=4, spanwise_distribution=UNCHANGED, spanwise_direction=[0.0, 1.0, 0.0], remove_nan=true, align_to_principal=false, alpha_range=deg2rad.(-5:1:20), delta_range=deg2rad.(-5:1:20), prn=true, - interp_steps=n_sections # TODO: check if interpolations are still needed + interp_steps=n_sections ) - !(n_groups == 0 || n_panels % n_groups == 0) && throw(ArgumentError("Number of panels should be divisible by number of groups")) !isapprox(spanwise_direction, [0.0, 1.0, 0.0]) && throw(ArgumentError("Spanwise direction has to be [0.0, 1.0, 0.0], not $spanwise_direction")) @@ -436,7 +389,7 @@ function RamAirWing( if align_to_principal inertia_tensor, R_cad_body = calc_inertia_y_rotation(inertia_tensor) else - R_cad_body = I(3) + R_cad_body = Matrix{Float64}(I, 3, 3) end circle_center_z, radius, gamma_tip = find_circle_center_and_radius(vertices) le_interp, te_interp, area_interp = create_interpolations(vertices, circle_center_z, radius, gamma_tip, R_cad_body; interp_steps) @@ -446,7 +399,7 @@ function RamAirWing( if !ispath(cl_polar_path) || !ispath(cd_polar_path) || !ispath(cm_polar_path) width = 2gamma_tip * radius area = area_interp(gamma_tip) - create_polars(; dat_path, cl_polar_path, cd_polar_path, cm_polar_path, wind_vel, + create_polars(; dat_path, cl_polar_path, cd_polar_path, cm_polar_path, wind_vel, area, width, crease_frac, alpha_range, delta_range, remove_nan) end @@ -459,7 +412,7 @@ function RamAirWing( any(isnan.(cd_matrix)) && interpolate_matrix_nans!(cd_matrix; prn) any(isnan.(cm_matrix)) && interpolate_matrix_nans!(cm_matrix; prn) end - + # Create sections sections = Section[] refined_sections = Section[] @@ -474,12 +427,13 @@ function RamAirWing( end panel_props = PanelProperties{n_panels}() - cache = [LazyBufferCache()] + cache = [PreallocationTools.LazyBufferCache()] - RamAirWing(n_panels, n_groups, spanwise_distribution, panel_props, spanwise_direction, sections, - refined_sections, remove_nan, non_deformed_sections, + Wing(n_panels, n_groups, spanwise_distribution, panel_props, MVec3(spanwise_direction), + sections, refined_sections, remove_nan, + non_deformed_sections, zeros(n_panels), zeros(n_panels), mass, gamma_tip, inertia_tensor, T_cad_body, R_cad_body, radius, - le_interp, te_interp, area_interp, zeros(n_panels), zeros(n_panels), cache) + le_interp, te_interp, area_interp, cache) catch e if e isa BoundsError @@ -488,110 +442,3 @@ function RamAirWing( rethrow(e) end end - -""" - group_deform!(wing::RamAirWing, theta_angles::AbstractVector, delta_angles::AbstractVector) - -Distribute control angles across wing panels and apply smoothing using a moving average filter. - -# Arguments -- `wing::RamAirWing`: The wing to deform -- `theta_angles::AbstractVector`: Twist angles in radians for each control section -- `delta_angles::AbstractVector`: Trailing edge deflection angles in radians for each control section -- `smooth::Bool`: Wether to apply smoothing or not - -# Algorithm -1. Distributes each control input to its corresponding group of panels -2. Applies moving average smoothing with window size based on control group size - -# Errors -- Throws `ArgumentError` if panel count is not divisible by the number of control inputs - -# Returns -- `nothing` (modifies wing in-place) -""" -function group_deform!(wing::RamAirWing, theta_angles=nothing, delta_angles=nothing; smooth=false) - !isnothing(theta_angles) && !(wing.n_panels % length(theta_angles) == 0) && - throw(ArgumentError("Number of angles has to be a multiple of number of panels")) - !isnothing(delta_angles) && !(wing.n_panels % length(delta_angles) == 0) && - throw(ArgumentError("Number of angles has to be a multiple of number of panels")) - isnothing(theta_angles) && isnothing(delta_angles) && return nothing - - n_panels = wing.n_panels - theta_dist = wing.theta_dist - delta_dist = wing.delta_dist - n_angles = isnothing(theta_angles) ? length(delta_angles) : length(theta_angles) - - dist_idx = 0 - for angle_idx in 1:n_angles - for _ in 1:(wing.n_panels ÷ n_angles) - dist_idx += 1 - !isnothing(theta_angles) && (theta_dist[dist_idx] = theta_angles[angle_idx]) - !isnothing(delta_angles) && (delta_dist[dist_idx] = delta_angles[angle_idx]) - end - end - @assert (dist_idx == wing.n_panels) - - if smooth - window_size = wing.n_panels ÷ n_angles - if n_panels > window_size - smoothed = wing.cache[1][theta_dist] - - if !isnothing(theta_angles) - smoothed .= theta_dist - for i in (window_size÷2 + 1):(n_panels - window_size÷2) - @views smoothed[i] = mean(theta_dist[(i - window_size÷2):(i + window_size÷2)]) - end - theta_dist .= smoothed - end - - if !isnothing(delta_angles) - smoothed .= delta_dist - for i in (window_size÷2 + 1):(n_panels - window_size÷2) - @views smoothed[i] = mean(delta_dist[(i - window_size÷2):(i + window_size÷2)]) - end - delta_dist .= smoothed - end - end - end - deform!(wing) - return nothing -end - -""" - deform!(wing::RamAirWing, theta_dist::AbstractVector, delta_dist::AbstractVector; width) - -Deform wing by applying theta and delta distributions. - -# Arguments -- `wing`: RamAirWing to deform -- `theta_dist`: the angle distribution between of the kite and the body x-axis in radians of each panel -- `delta_dist`: the deformation of the trailing edges of each panel - -# Effects -Updates wing.sections with deformed geometry -""" -function deform!(wing::RamAirWing, theta_dist::AbstractVector, delta_dist::AbstractVector) - !(length(theta_dist) == wing.n_panels) && throw(ArgumentError("theta_dist and panels are of different lengths")) - !(length(delta_dist) == wing.n_panels) && throw(ArgumentError("delta_dist and panels are of different lengths")) - wing.theta_dist .= theta_dist - wing.delta_dist .= delta_dist - - deform!(wing) -end - -function deform!(wing::RamAirWing) - local_y = zeros(MVec3) - chord = zeros(MVec3) - normal = zeros(MVec3) - - for i in 1:wing.n_panels - section1 = wing.non_deformed_sections[i] - section2 = wing.non_deformed_sections[i+1] - local_y .= normalize(section1.LE_point - section2.LE_point) - chord .= section1.TE_point .- section1.LE_point - normal .= chord × local_y - @. wing.sections[i].TE_point = section1.LE_point + cos(wing.theta_dist[i]) * chord - sin(wing.theta_dist[i]) * normal - end - return nothing -end diff --git a/src/solver.jl b/src/solver.jl index 8b2d75fd..9bce66ae 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -654,15 +654,15 @@ function smooth_circulation!( end """ - linearize(solver::Solver, body_aero::BodyAerodynamics, wing::RamAirWing, y::Vector{T}; + linearize(solver::Solver, body_aero::BodyAerodynamics, wing::Wing, y::Vector{T}; theta_idxs=1:4, delta_idxs=nothing, va_idxs=nothing, omega_idxs=nothing, kwargs...) where T -Compute the Jacobian matrix for a ram air wing around an operating point using finite differences. +Compute the Jacobian matrix for a deformable wing around an operating point using finite differences. # Arguments - `solver`: VSM solver instance (must be initialized) - `body_aero`: Aerodynamic body representation -- `wing`: RamAirWing model to linearize +- `wing`: Wing model to linearize (must support deformation, i.e., created with ObjWing()) - `y`: Input vector at operating point, containing a combination of control angles and velocities # Keyword Arguments @@ -679,7 +679,7 @@ Compute the Jacobian matrix for a ram air wing around an operating point using f # Example ```julia # Initialize wing and solver -wing = RamAirWing("path/to/body.obj", "path/to/foil.dat") +wing = ObjWing("path/to/body.obj", "path/to/foil.dat") body_aero = BodyAerodynamics([wing]) solver = Solver(body_aero) diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index bb1f5e06..02f7a308 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -192,7 +192,7 @@ end Represents a wing composed of multiple sections with aerodynamic properties. -# Fields +# Core Fields (all wings) - `n_panels::Int16`: Number of panels in aerodynamic mesh - `n_groups::Int16`: Number of panel groups - `spanwise_distribution`::PanelDistribution: [PanelDistribution](@ref) @@ -201,6 +201,23 @@ Represents a wing composed of multiple sections with aerodynamic properties. - `refined_sections::Vector{Section}`: Vector of refined wing sections, see: [Section](@ref) - `remove_nan::Bool`: Wether to remove the NaNs from interpolations or not +# Deformation Fields (optional, for deformable wings) +- `non_deformed_sections::Vector{Section}`: Original undeformed sections +- `theta_dist::Vector{Float64}`: Panel twist angle distribution +- `delta_dist::Vector{Float64}`: Trailing edge deflection distribution + +# Physical Properties (optional, for OBJ-based wings) +- `mass::Float64`: Total wing mass in kg (0.0 if not applicable) +- `gamma_tip::Float64`: Angular extent from center to wing tip (0.0 if not applicable) +- `inertia_tensor::Matrix{Float64}`: 3x3 inertia tensor (empty if not applicable) +- `T_cad_body::MVec3`: Translation from CAD to body frame (zeros if not applicable) +- `R_cad_body::MMat3`: Rotation from CAD to body frame (identity if not applicable) +- `radius::Float64`: Wing curvature radius (0.0 if not applicable) +- `le_interp::Union{Nothing, NTuple{3, Extrapolation}}`: Leading edge interpolation +- `te_interp::Union{Nothing, NTuple{3, Extrapolation}}`: Trailing edge interpolation +- `area_interp::Union{Nothing, Extrapolation}`: Area interpolation +- `cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}}`: Preallocated buffers + """ mutable struct Wing <: AbstractWing n_panels::Int16 @@ -211,6 +228,23 @@ mutable struct Wing <: AbstractWing sections::Vector{Section} refined_sections::Vector{Section} remove_nan::Bool + + # Deformation fields + non_deformed_sections::Vector{Section} + theta_dist::Vector{Float64} + delta_dist::Vector{Float64} + + # Physical properties (OBJ-based wings) + mass::Float64 + gamma_tip::Float64 + inertia_tensor::Matrix{Float64} + T_cad_body::MVec3 + R_cad_body::MMat3 + radius::Float64 + le_interp::Union{Nothing, NTuple{3, Extrapolation}} + te_interp::Union{Nothing, NTuple{3, Extrapolation}} + area_interp::Union{Nothing, Extrapolation} + cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} end """ @@ -220,8 +254,8 @@ end spanwise_direction::PosVector=MVec3([0.0, 1.0, 0.0]), remove_nan::Bool=true) -Constructor for a [Wing](@ref) struct with default values that initializes the sections -and refined sections as empty arrays. +Constructor for a [Wing](@ref) struct with default values that initializes the sections +and refined sections as empty arrays. Creates a basic wing suitable for YAML-based construction. # Parameters - `n_panels::Int`: Number of panels in aerodynamic mesh @@ -237,12 +271,23 @@ function Wing(n_panels::Int; remove_nan=true) !(n_groups == 0 || n_panels % n_groups == 0) && throw(ArgumentError("Number of panels should be divisible by number of groups")) panel_props = PanelProperties{n_panels}() - Wing(n_panels, n_groups, spanwise_distribution, panel_props, spanwise_direction, Section[], Section[], remove_nan) + + # Initialize with default/empty values for optional fields + Wing( + n_panels, n_groups, spanwise_distribution, panel_props, spanwise_direction, + Section[], Section[], remove_nan, + # Deformation fields + Section[], zeros(n_panels), zeros(n_panels), + # Physical properties (defaults for non-OBJ wings) + 0.0, 0.0, zeros(0, 0), zeros(MVec3), Matrix{Float64}(I, 3, 3), + 0.0, nothing, nothing, nothing, + PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}[] + ) end function reinit!(wing::AbstractWing) refine_aerodynamic_mesh!(wing) - + # Calculate panel properties update_panel_properties!( wing.panel_props, @@ -252,6 +297,148 @@ function reinit!(wing::AbstractWing) return nothing end +""" + group_deform!(wing::Wing, theta_angles=nothing, delta_angles=nothing; smooth=false) + +Distribute control angles across wing panels and optionally apply smoothing. + +For wings that support deformation (OBJ-based wings with non_deformed_sections), this +distributes theta_angles and delta_angles to panel groups and applies deformation. +For wings without deformation support (YAML-based), this is a no-op that only succeeds +if both angle inputs are nothing. + +# Arguments +- `wing::Wing`: The wing to deform +- `theta_angles::AbstractVector`: Twist angles in radians for each control group (or nothing) +- `delta_angles::AbstractVector`: Trailing edge deflection angles in radians (or nothing) +- `smooth::Bool`: Whether to apply moving average smoothing + +# Algorithm +1. Distributes each control input to its corresponding group of panels +2. Optionally applies moving average smoothing with window based on group size +3. Calls deform! to update wing geometry + +# Errors +- Throws `ArgumentError` if wing doesn't support deformation but angles are provided +- Throws `ArgumentError` if panel count is not divisible by number of control inputs + +# Returns +- `nothing` (modifies wing in-place) +""" +function group_deform!(wing::Wing, theta_angles=nothing, delta_angles=nothing; smooth=false) + # Check if deformation is supported + can_deform = !isempty(wing.non_deformed_sections) + + # If no deformation requested, just return + isnothing(theta_angles) && isnothing(delta_angles) && return nothing + + # If deformation requested but not supported, throw error + if !can_deform + throw(ArgumentError("This Wing does not support deformation. Only OBJ-based wings created with ObjWing() can be deformed.")) + end + + # Validate inputs + !isnothing(theta_angles) && !(wing.n_panels % length(theta_angles) == 0) && + throw(ArgumentError("Number of theta_angles has to be a divisor of number of panels")) + !isnothing(delta_angles) && !(wing.n_panels % length(delta_angles) == 0) && + throw(ArgumentError("Number of delta_angles has to be a divisor of number of panels")) + + n_panels = wing.n_panels + theta_dist = wing.theta_dist + delta_dist = wing.delta_dist + n_angles = isnothing(theta_angles) ? length(delta_angles) : length(theta_angles) + + # Distribute angles to panels + dist_idx = 0 + for angle_idx in 1:n_angles + for _ in 1:(wing.n_panels ÷ n_angles) + dist_idx += 1 + !isnothing(theta_angles) && (theta_dist[dist_idx] = theta_angles[angle_idx]) + !isnothing(delta_angles) && (delta_dist[dist_idx] = delta_angles[angle_idx]) + end + end + @assert (dist_idx == wing.n_panels) + + # Apply smoothing if requested + if smooth + window_size = wing.n_panels ÷ n_angles + if n_panels > window_size + smoothed = wing.cache[1][theta_dist] + + if !isnothing(theta_angles) + smoothed .= theta_dist + for i in (window_size÷2 + 1):(n_panels - window_size÷2) + @views smoothed[i] = mean(theta_dist[(i - window_size÷2):(i + window_size÷2)]) + end + theta_dist .= smoothed + end + + if !isnothing(delta_angles) + smoothed .= delta_dist + for i in (window_size÷2 + 1):(n_panels - window_size÷2) + @views smoothed[i] = mean(delta_dist[(i - window_size÷2):(i + window_size÷2)]) + end + delta_dist .= smoothed + end + end + end + + deform!(wing) + return nothing +end + +""" + deform!(wing::Wing, theta_dist::AbstractVector, delta_dist::AbstractVector) + +Deform wing by applying theta and delta distributions directly. + +# Arguments +- `wing::Wing`: Wing to deform (must support deformation) +- `theta_dist::AbstractVector`: Twist angle in radians for each panel +- `delta_dist::AbstractVector`: Trailing edge deflection for each panel + +# Effects +Updates wing.sections with deformed geometry based on wing.non_deformed_sections +""" +function deform!(wing::Wing, theta_dist::AbstractVector, delta_dist::AbstractVector) + !isempty(wing.non_deformed_sections) || throw(ArgumentError("Wing does not support deformation")) + !(length(theta_dist) == wing.n_panels) && throw(ArgumentError("theta_dist and panels are of different lengths")) + !(length(delta_dist) == wing.n_panels) && throw(ArgumentError("delta_dist and panels are of different lengths")) + wing.theta_dist .= theta_dist + wing.delta_dist .= delta_dist + + deform!(wing) +end + +""" + deform!(wing::Wing) + +Apply stored theta_dist and delta_dist to deform the wing geometry. + +# Arguments +- `wing::Wing`: Wing to deform (must have non_deformed_sections) + +# Effects +Updates wing.sections based on wing.non_deformed_sections and stored distributions +""" +function deform!(wing::Wing) + !isempty(wing.non_deformed_sections) || return nothing + + local_y = zeros(MVec3) + chord = zeros(MVec3) + normal = zeros(MVec3) + + for i in 1:wing.n_panels + section1 = wing.non_deformed_sections[i] + section2 = wing.non_deformed_sections[i+1] + local_y .= normalize(section1.LE_point - section2.LE_point) + chord .= section1.TE_point .- section1.LE_point + normal .= chord × local_y + @. wing.sections[i].TE_point = section1.LE_point + cos(wing.theta_dist[i]) * chord - sin(wing.theta_dist[i]) * normal + end + return nothing +end + """ remove_vector_nans(aero_data) diff --git a/test/body_aerodynamics/test_results.jl b/test/body_aerodynamics/test_results.jl index 2b47dc5b..60d2d4cf 100644 --- a/test/body_aerodynamics/test_results.jl +++ b/test/body_aerodynamics/test_results.jl @@ -31,7 +31,7 @@ if !@isdefined ram_wing_results error("Required data files not found: $body_src or $foil_src") end - ram_wing = RamAirWing(body_path, foil_path; alpha_range=deg2rad.(-1:1), delta_range=deg2rad.(-1:1)) + ram_wing = ObjWing(body_path, foil_path; alpha_range=deg2rad.(-1:1), delta_range=deg2rad.(-1:1)) end @testset "Nonlinear vs Linear - Comprehensive Input Testing" begin diff --git a/test/plotting/test_plotting.jl b/test/plotting/test_plotting.jl index 385abcef..e1f82959 100644 --- a/test/plotting/test_plotting.jl +++ b/test/plotting/test_plotting.jl @@ -29,7 +29,7 @@ if !@isdefined ram_wing foil_src = joinpath(_ram_data_dir, "ram_air_kite_foil.dat") cp(body_src, body_path; force=true) cp(foil_src, foil_path; force=true) - ram_wing = RamAirWing(body_path, foil_path; alpha_range=deg2rad.(-1:1), delta_range=deg2rad.(-1:1)) + ram_wing = ObjWing(body_path, foil_path; alpha_range=deg2rad.(-1:1), delta_range=deg2rad.(-1:1)) end function create_body_aero() diff --git a/test/ram_geometry/test_kite_geometry.jl b/test/ram_geometry/test_kite_geometry.jl index 74b6fd9c..47dee1d4 100644 --- a/test/ram_geometry/test_kite_geometry.jl +++ b/test/ram_geometry/test_kite_geometry.jl @@ -165,9 +165,9 @@ using Serialization @test R_b_p2 ≈ I(3) end - @testset "RamAirWing Construction" begin - wing = RamAirWing(test_obj_path, test_dat_path; remove_nan=true) - + @testset "ObjWing Construction" begin + wing = ObjWing(test_obj_path, test_dat_path; remove_nan=true) + @test wing.n_panels == 56 # Default value @test wing.spanwise_distribution == UNCHANGED @test wing.spanwise_direction ≈ [0.0, 1.0, 0.0] @@ -179,15 +179,15 @@ using Serialization @test !isnan(wing.sections[1].aero_data[4][end]) @test !isnan(wing.sections[1].aero_data[5][end]) - wing = RamAirWing(test_obj_path, test_dat_path; remove_nan=false) + wing = ObjWing(test_obj_path, test_dat_path; remove_nan=false) @test isnan(wing.sections[1].aero_data[3][end]) @test isnan(wing.sections[1].aero_data[4][end]) @test isnan(wing.sections[1].aero_data[5][end]) end @testset "Wing Deformation" begin - # Create a RamAirWing for testing - wing = RamAirWing(test_obj_path, test_dat_path; remove_nan=true) + # Create an ObjWing for testing + wing = ObjWing(test_obj_path, test_dat_path; remove_nan=true) body_aero = BodyAerodynamics([wing]) # Store original TE point for comparison diff --git a/test/solver/test_solver.jl b/test/solver/test_solver.jl index 0dbb9971..084be994 100644 --- a/test/solver/test_solver.jl +++ b/test/solver/test_solver.jl @@ -53,7 +53,7 @@ using Test sol = solve!(solver, body_aero) @test sol isa VSMSolution - @test sol.solver_status == SUCCESS + @test sol.solver_status == FEASIBLE # Verify that group moments are empty @test length(sol.group_moment_dist) == 0 @@ -177,4 +177,47 @@ using Test rm(settings_file; force=true) end end + + @testset "Wing type cannot deform" begin + # Test that regular Wing type (YAML-based) cannot use group_deform! + settings_file = create_temp_wing_settings("solver", "solver_test_wing.yaml"; + alpha=5.0, beta=0.0, wind_speed=10.0, n_groups=2, n_panels=4) + + try + settings = VSMSettings(settings_file) + wing = Wing(settings) + @test wing isa Wing + @test wing.n_groups == 2 + + body_aero = BodyAerodynamics([wing]) + solver = Solver(body_aero, settings) + + va = [10.0, 0.0, 0.0] + set_va!(body_aero, va) + + # Test that trying to use theta angles with Wing throws an error + y = [0.0, 0.0, 10.0, 0.0, 0.0] # 2 theta + 3 va + @test_throws ArgumentError VortexStepMethod.linearize( + solver, + body_aero, + y; + theta_idxs=1:2, + va_idxs=3:5 + ) + + # But linearize should work fine with only velocity (no deformation) + y_velocity_only = [10.0, 0.0, 0.0] + jac, results = VortexStepMethod.linearize( + solver, + body_aero, + y_velocity_only; + theta_idxs=nothing, + va_idxs=1:3 + ) + @test size(jac, 1) == 8 # 6 + 2 group moments + + finally + rm(settings_file; force=true) + end + end end diff --git a/test/wake/test_wake.jl b/test/wake/test_wake.jl index 156e4ec9..b1c45f87 100644 --- a/test/wake/test_wake.jl +++ b/test/wake/test_wake.jl @@ -20,7 +20,7 @@ using VortexStepMethod try # Create wing and body aerodynamics with known good geometry - wing = RamAirWing(body_path, foil_path; n_panels=56) # Use default panels + wing = ObjWing(body_path, foil_path; n_panels=56) # Use default panels body_aero = BodyAerodynamics([wing]) # Test that frozen_wake! doesn't throw errors From 15811702ff937442d915c86bbe5b24e27a34b9c4 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sun, 26 Oct 2025 13:52:56 +0100 Subject: [PATCH 03/22] Restore test --- test/solver/test_solver.jl | 191 ------------------------------------- 1 file changed, 191 deletions(-) diff --git a/test/solver/test_solver.jl b/test/solver/test_solver.jl index 084be994..1d60061a 100644 --- a/test/solver/test_solver.jl +++ b/test/solver/test_solver.jl @@ -29,195 +29,4 @@ using Test rm(settings_file; force=true) end end - - @testset "Solver with n_groups=0" begin - # Test that solver works correctly when n_groups=0 (no group functionality) - settings_file = create_temp_wing_settings("solver", "solver_test_wing.yaml"; - alpha=5.0, beta=0.0, wind_speed=10.0, n_groups=0) - - try - settings = VSMSettings(settings_file) - wing = Wing(settings) - @test wing.n_groups == 0 - - body_aero = BodyAerodynamics([wing]) - solver = Solver(body_aero, settings) - - # Verify solver has zero groups - @test length(solver.sol.group_moment_dist) == 0 - @test length(solver.sol.group_moment_coeff_dist) == 0 - - # Test that solve! works without errors - va = [10.0, 0.0, 0.0] - set_va!(body_aero, va) - sol = solve!(solver, body_aero) - - @test sol isa VSMSolution - @test sol.solver_status == FEASIBLE - - # Verify that group moments are empty - @test length(sol.group_moment_dist) == 0 - @test length(sol.group_moment_coeff_dist) == 0 - - # But force and moment should still be computed - @test !all(sol.force .== 0.0) - @test norm(sol.force) > 0 - - finally - rm(settings_file; force=true) - end - end - - @testset "Linearize with n_groups=0" begin - # Test that linearize works correctly when n_groups=0 - settings_file = create_temp_wing_settings("solver", "solver_test_wing.yaml"; - alpha=5.0, beta=0.0, wind_speed=10.0, n_groups=0, n_panels=4) - - try - settings = VSMSettings(settings_file) - wing = Wing(settings) - @test wing.n_groups == 0 - - body_aero = BodyAerodynamics([wing]) - solver = Solver(body_aero, settings) - - # Set velocity - va = [10.0, 0.0, 0.0] - set_va!(body_aero, va) - - # Test linearize with only velocity (no theta or delta since n_groups=0) - y = va # Only velocity in input vector - jac, results = VortexStepMethod.linearize( - solver, - body_aero, - y; - theta_idxs=nothing, - delta_idxs=nothing, - va_idxs=1:3, - omega_idxs=nothing - ) - - # Results should only have 6 elements (force + moment, no group moments) - @test length(results) == 6 - @test size(jac) == (6, 3) # 6 outputs, 3 inputs (vx, vy, vz) - - # Verify forces are non-zero - @test norm(results[1:3]) > 0 - - # Test that using theta_idxs with n_groups=0 throws an error - @test_throws ArgumentError VortexStepMethod.linearize( - solver, - body_aero, - [0.0, 10.0, 0.0, 0.0]; # Invalid: trying to use theta - theta_idxs=1:1, - va_idxs=2:4 - ) - - # Test that using delta_idxs with n_groups=0 throws an error - @test_throws ArgumentError VortexStepMethod.linearize( - solver, - body_aero, - [0.0, 10.0, 0.0, 0.0]; # Invalid: trying to use delta - delta_idxs=1:1, - va_idxs=2:4 - ) - - finally - rm(settings_file; force=true) - end - end - - @testset "Linearize theta_idxs validation" begin - # Test that theta_idxs length must match n_groups - settings_file = create_temp_wing_settings("solver", "solver_test_wing.yaml"; - alpha=5.0, beta=0.0, wind_speed=10.0, n_groups=2, n_panels=4) - - try - settings = VSMSettings(settings_file) - wing = Wing(settings) - @test wing.n_groups == 2 - - body_aero = BodyAerodynamics([wing]) - solver = Solver(body_aero, settings) - - va = [10.0, 0.0, 0.0] - set_va!(body_aero, va) - - # Test with correct number of theta angles (2) - y_correct = [0.0, 0.0, 10.0, 0.0, 0.0] # 2 theta + 3 va - jac, results = VortexStepMethod.linearize( - solver, - body_aero, - y_correct; - theta_idxs=1:2, - va_idxs=3:5 - ) - @test size(jac, 1) == 8 # 6 + 2 group moments - - # Test with wrong number of theta angles (should throw error) - y_wrong = [0.0, 0.0, 0.0, 0.0, 10.0, 0.0, 0.0] # 4 theta + 3 va - @test_throws ArgumentError VortexStepMethod.linearize( - solver, - body_aero, - y_wrong; - theta_idxs=1:4, # Wrong: 4 angles but only 2 groups - va_idxs=5:7 - ) - - # Test with wrong number of delta angles - @test_throws ArgumentError VortexStepMethod.linearize( - solver, - body_aero, - y_wrong; - delta_idxs=1:4, # Wrong: 4 angles but only 2 groups - va_idxs=5:7 - ) - - finally - rm(settings_file; force=true) - end - end - - @testset "Wing type cannot deform" begin - # Test that regular Wing type (YAML-based) cannot use group_deform! - settings_file = create_temp_wing_settings("solver", "solver_test_wing.yaml"; - alpha=5.0, beta=0.0, wind_speed=10.0, n_groups=2, n_panels=4) - - try - settings = VSMSettings(settings_file) - wing = Wing(settings) - @test wing isa Wing - @test wing.n_groups == 2 - - body_aero = BodyAerodynamics([wing]) - solver = Solver(body_aero, settings) - - va = [10.0, 0.0, 0.0] - set_va!(body_aero, va) - - # Test that trying to use theta angles with Wing throws an error - y = [0.0, 0.0, 10.0, 0.0, 0.0] # 2 theta + 3 va - @test_throws ArgumentError VortexStepMethod.linearize( - solver, - body_aero, - y; - theta_idxs=1:2, - va_idxs=3:5 - ) - - # But linearize should work fine with only velocity (no deformation) - y_velocity_only = [10.0, 0.0, 0.0] - jac, results = VortexStepMethod.linearize( - solver, - body_aero, - y_velocity_only; - theta_idxs=nothing, - va_idxs=1:3 - ) - @test size(jac, 1) == 8 # 6 + 2 group moments - - finally - rm(settings_file; force=true) - end - end end From df97a2187674f727b2cbb82c1b851baa97a98f23 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Wed, 29 Oct 2025 14:52:52 +0100 Subject: [PATCH 04/22] Not passing bench --- src/body_aerodynamics.jl | 72 ++++++++++++++++++++++++++-------------- 1 file changed, 47 insertions(+), 25 deletions(-) diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 6444a598..9e09f54c 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -129,13 +129,13 @@ Initialize a BodyAerodynamics struct in-place by setting up panels and coefficie # Returns nothing """ -function reinit!(body_aero::BodyAerodynamics; +function reinit!(body_aero::BodyAerodynamics; init_aero=true, va=[15.0, 0.0, 0.0], omega=zeros(MVec3) ) idx = 1 - vec = zeros(MVec3) + vec = @MVector zeros(3) for wing in body_aero.wings reinit!(wing) panel_props = wing.panel_props @@ -168,8 +168,8 @@ function reinit!(body_aero::BodyAerodynamics; end # Initialize rest of the struct - body_aero.projected_area = sum(wing -> calculate_projected_area(wing), body_aero.wings) - body_aero.stall_angle_list .= calculate_stall_angle_list(body_aero.panels) + body_aero.projected_area = sum(calculate_projected_area, body_aero.wings) + calculate_stall_angle_list!(body_aero.stall_angle_list, body_aero.panels) body_aero.alpha_array .= 0.0 body_aero.v_a_array .= 0.0 body_aero.AIC .= 0.0 @@ -191,21 +191,28 @@ Returns: nothing """ @inline function calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::Model, core_radius_fraction, - va_norm_array, + va_norm_array, va_unit_array) # Determine evaluation point based on model evaluation_point = model == VSM ? :control_point : :aero_center evaluation_point_on_bound = model == LLT - - # Initialize AIC matrices - velocity_induced, tempvel, va_unit, U_2D = zeros(MVec3), zeros(MVec3), zeros(MVec3), zeros(MVec3) + + # Allocate work vectors for this function (separate from those used by child functions) + velocity_induced = @MVector zeros(3) + tempvel = @MVector zeros(3) + va_unit = @MVector zeros(3) + U_2D = @MVector zeros(3) # Calculate influence coefficients for icp in eachindex(body_aero.panels) - ep = getproperty(body_aero.panels[icp], evaluation_point) + panel_icp = body_aero.panels[icp] + ep = evaluation_point == :control_point ? panel_icp.control_point : panel_icp.aero_center for jring in eachindex(body_aero.panels) - va_unit .= @views va_unit_array[jring, :] - filaments = body_aero.panels[jring].filaments + panel_jring = body_aero.panels[jring] + @inbounds for k in 1:3 + va_unit[k] = va_unit_array[jring, k] + end + filaments = panel_jring.filaments va_norm = va_norm_array[jring] calculate_velocity_induced_single_ring_semiinfinite!( velocity_induced, @@ -222,8 +229,8 @@ Returns: nothing # Subtract 2D induced velocity for VSM if icp == jring && model == VSM - calculate_velocity_induced_bound_2D!(U_2D, body_aero.panels[jring], ep, body_aero.work_vectors) - velocity_induced .-= U_2D + calculate_velocity_induced_bound_2D!(U_2D, panel_jring, ep, body_aero.work_vectors) + velocity_induced .-= U_2D end body_aero.AIC[:, icp, jring] .= velocity_induced end @@ -276,19 +283,34 @@ function calculate_stall_angle_list(panels::Vector{Panel}; step_aoa=1.0, stall_angle_if_none_detected=50.0, cl_initial=-10.0) - - aoa_range = deg2rad.(range(begin_aoa, end_aoa, step=step_aoa)) - stall_angles = Float64[] - - for panel in panels + stall_angles = Vector{Float64}(undef, length(panels)) + calculate_stall_angle_list!(stall_angles, panels; + begin_aoa, end_aoa, step_aoa, + stall_angle_if_none_detected, cl_initial) + return stall_angles +end + +function calculate_stall_angle_list!(stall_angles::AbstractVector{Float64}, + panels::Vector{Panel}; + begin_aoa=9.0, + end_aoa=22.0, + step_aoa=1.0, + stall_angle_if_none_detected=50.0, + cl_initial=-10.0) + + # Pre-compute range values to avoid allocation + n_steps = Int(floor((end_aoa - begin_aoa) / step_aoa)) + 1 + + for (idx, panel) in enumerate(panels) # Default stall angle if none found panel_stall = stall_angle_if_none_detected - + # Start with minimum cl cl_old = cl_initial - + # Find stall angle - for aoa in aoa_range + for i in 0:(n_steps-1) + aoa = deg2rad(begin_aoa + i * step_aoa) cl = calculate_cl(panel, aoa) if cl < cl_old panel_stall = aoa @@ -296,11 +318,11 @@ function calculate_stall_angle_list(panels::Vector{Panel}; end cl_old = cl end - - push!(stall_angles, panel_stall) + + stall_angles[idx] = panel_stall end - - return stall_angles + + return nothing end """ From e28c0ced61539c023e5c2e74c22e9efbf407d42b Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Wed, 29 Oct 2025 15:10:21 +0100 Subject: [PATCH 05/22] Working grouping --- README.md | 4 +- bin/install | 2 +- docs/make.jl | 6 -- docs/src/glossary.md | 1 + docs/src/index.md | 4 +- docs/src/tips_and_tricks.md | 32 +++++++++ docs/src/types.md | 1 + examples/Project.toml | 7 ++ examples/bench.jl | 6 -- examples/menu.jl | 5 -- examples/stall_model.jl | 8 +-- scripts/Project.toml | 4 ++ scripts/build_docu.jl | 4 -- src/VortexStepMethod.jl | 12 ++++ src/obj_geometry.jl | 3 +- src/solver.jl | 24 +++++-- src/wing_geometry.jl | 114 +++++++++++++++++++++++++++++---- src/yaml_geometry.jl | 12 ++-- test/Aqua.jl | 5 -- test/Project.toml | 10 +++ test/bench.jl | 5 -- test/bench_solve.jl | 6 -- test/settings/test_settings.jl | 5 -- 23 files changed, 204 insertions(+), 76 deletions(-) create mode 100644 examples/Project.toml create mode 100644 scripts/Project.toml create mode 100644 test/Project.toml diff --git a/README.md b/README.md index de01a449..d93c3a0a 100644 --- a/README.md +++ b/README.md @@ -19,9 +19,9 @@ if you haven't already. On Linux, make sure that Python3 and Matplotlib are inst ``` sudo apt install python3-matplotlib ``` -Furthermore, the packages `TestEnv` and `ControlPlots` must be installed globally: +Furthermore, the package `ControlPlots` must be installed globally: ``` -julia -e 'using Pkg; Pkg.add("TestEnv"); Pkg.add("ControlPlots")' +julia -e 'using Pkg; Pkg.add("ControlPlots")' ``` Before installing this software it is suggested to create a new project, for example like this: diff --git a/bin/install b/bin/install index e4609513..49b359ad 100755 --- a/bin/install +++ b/bin/install @@ -20,7 +20,7 @@ fi export JULIA_PKG_SERVER_REGISTRY_PREFERENCE=eager -julia -e 'using Pkg; Pkg.add("TestEnv"); Pkg.add("ControlPlots")' +julia -e 'using Pkg; Pkg.add("ControlPlots")' julia --project -e 'include("bin/install.jl")' diff --git a/docs/make.jl b/docs/make.jl index 17d42209..677baf72 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,9 +1,3 @@ -using Pkg -if ("TestEnv" ∈ keys(Pkg.project().dependencies)) - if ! ("Documents" ∈ keys(Pkg.project().dependencies)) - using TestEnv; TestEnv.activate() - end -end using ControlPlots using VortexStepMethod using Documenter diff --git a/docs/src/glossary.md b/docs/src/glossary.md index a67c577d..940c6b5d 100644 --- a/docs/src/glossary.md +++ b/docs/src/glossary.md @@ -5,6 +5,7 @@ | AIC | Aerodynamic Influence Coefficient (AIC). The AIC matrix represents the relationship between the induced velocities or pressures on aerodynamic surfaces and the circulation strength or modal deformations of the lifting surfaces.| | inviscid | A fluid flow in which viscosity is considered negligible or zero. This means that there is no internal friction between the fluid layers, and the effects of viscosity on the flow are assumed to be insignificant. | | Panel | Flat surface element in 3D that approximate the contour of the aerodynamic body being studied.| +| Panel Group | A collection of panels whose aerodynamic forces and moments are summed together. Groups can be defined using EQUAL_SIZE (sequential grouping) or REFINE (based on original unrefined structure) methods.| | Section |A wing section, also known as an airfoil or aerofoil, is the cross-sectional shape of an aircraft wing.| | Span | Distance from one wing tip to the other wing tip. | | Polar | The polar typically plots the coefficient of lift (CL) against the coefficient of drag (CD), with the angle of attack as a parameter along the curve. | diff --git a/docs/src/index.md b/docs/src/index.md index fa6d54c1..3310d72c 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -18,9 +18,9 @@ if you haven't already. On Linux, make sure that Python3 and Matplotlib are inst ``` sudo apt install python3-matplotlib ``` -Furthermore, the packages `TestEnv` and `ControlPlots` must be installed globally: +Furthermore, the package `ControlPlots` must be installed globally: ``` -julia -e 'using Pkg; Pkg.add("TestEnv"); Pkg.add("ControlPlots")' +julia -e 'using Pkg; Pkg.add("ControlPlots")' ``` Before installing this software it is suggested to create a new project, for example like this: diff --git a/docs/src/tips_and_tricks.md b/docs/src/tips_and_tricks.md index 64995f34..63cf87d3 100644 --- a/docs/src/tips_and_tricks.md +++ b/docs/src/tips_and_tricks.md @@ -10,6 +10,38 @@ The following bodies can be simulated: To build the geometry of a RAM-air kite, a 3D .obj file can be used as input. In addition a `.dat` file is needed. It should have two columns, one for the `x` and one for the `y` coordinate of the 2D polar that is used. +## Panel Grouping Methods +When creating a wing, you can specify how panels should be grouped for moment and force calculations using the `grouping_method` parameter. Two methods are available: + +### EQUAL_SIZE (Default) +Divides refined panels into equally-sized sequential groups. This is the original behavior. + +```julia +wing = Wing(40; n_groups=4, grouping_method=EQUAL_SIZE) +``` + +In this example, with 40 panels and 4 groups, each group will contain 10 consecutive panels (panels 1-10, 11-20, 21-30, 31-40). + +### REFINE +Groups refined panels back to their original unrefined section. This is useful when you want group moments and forces to represent the original wing structure, regardless of panel refinement. + +```julia +# Create wing with 4 unrefined sections (3 panels) +wing = Wing(40; n_groups=3, grouping_method=REFINE) +add_section!(wing, [0, 5, 0], [1, 5, 0], INVISCID) # Section 1 +add_section!(wing, [0, 2.5, 0], [1, 2.5, 0], INVISCID) # Section 2 +add_section!(wing, [0, 0, 0], [1, 0, 0], INVISCID) # Section 3 +add_section!(wing, [0, -5, 0], [1, -5, 0], INVISCID) # Section 4 +``` + +**Important:** When using `REFINE`, `n_groups` must equal the number of unrefined panels (number of sections - 1). The solver will automatically map each refined panel to its closest original unrefined panel and sum their moments and forces accordingly. + +This is particularly useful for: +- LEI kites where you want loads per rib +- Wings with discrete control surfaces +- Cases where physical structure doesn't align with uniform panel distribution +- Dynamic simulations where you have fewer structural segments than panels needed for accurate VSM aerodynamics. For example, a 6-segment structural model can be combined with 40-panel aerodynamics by using `n_groups=6` and `grouping_method=REFINE` to map aerodynamic loads back to the structural segments. + ## RAM-air kite model If running the example `ram_air_kite.jl` fails, try to run the `cleanup.jl` script and then try again. Background: this example caches the calculated polars. Reading cached polars can fail after an update. diff --git a/docs/src/types.md b/docs/src/types.md index de083b35..9b76a7d0 100644 --- a/docs/src/types.md +++ b/docs/src/types.md @@ -7,6 +7,7 @@ Model WingType AeroModel PanelDistribution +PanelGroupingMethod InitialGammaDistribution SolverStatus ``` diff --git a/examples/Project.toml b/examples/Project.toml new file mode 100644 index 00000000..b94307b9 --- /dev/null +++ b/examples/Project.toml @@ -0,0 +1,7 @@ +[deps] +CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +ControlPlots = "23c2ee80-7a9e-4350-b264-8e670f12517c" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" +VortexStepMethod = "ed3cd733-9f0f-46a9-93e0-89b8d4998dd9" diff --git a/examples/bench.jl b/examples/bench.jl index 90c41ff1..6c2f8726 100644 --- a/examples/bench.jl +++ b/examples/bench.jl @@ -2,12 +2,6 @@ using LinearAlgebra using ControlPlots using VortexStepMethod -using Pkg - -if !("CSV" ∈ keys(Pkg.project().dependencies)) - using TestEnv - TestEnv.activate() -end # Step 1: Define wing parameters n_panels = 20 # Number of panels diff --git a/examples/menu.jl b/examples/menu.jl index 4454ccee..147990c3 100644 --- a/examples/menu.jl +++ b/examples/menu.jl @@ -1,8 +1,3 @@ -using Pkg -if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies)) - using TestEnv; TestEnv.activate() -end - using ControlPlots using VortexStepMethod using REPL.TerminalMenus diff --git a/examples/stall_model.jl b/examples/stall_model.jl index 6253faca..89f3c0f4 100644 --- a/examples/stall_model.jl +++ b/examples/stall_model.jl @@ -2,10 +2,6 @@ using ControlPlots using LinearAlgebra using VortexStepMethod -using Pkg -if ! ("CSV" ∈ keys(Pkg.project().dependencies)) - using TestEnv; TestEnv.activate() -end using CSV using DataFrames @@ -38,7 +34,9 @@ for row in eachrow(df) end # Create wing geometry -CAD_wing = Wing(n_panels; spanwise_distribution) +# Using REFINE grouping method: n_groups should equal number of unrefined panels (18 sections = 18 panels) +n_groups = length(rib_list) - 1 +CAD_wing = Wing(n_panels; spanwise_distribution, n_groups, grouping_method=REFINE) for rib in rib_list add_section!(CAD_wing, rib[1], rib[2], rib[3], rib[4]) end diff --git a/scripts/Project.toml b/scripts/Project.toml new file mode 100644 index 00000000..f020312d --- /dev/null +++ b/scripts/Project.toml @@ -0,0 +1,4 @@ +[deps] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589" +VortexStepMethod = "ed3cd733-9f0f-46a9-93e0-89b8d4998dd9" diff --git a/scripts/build_docu.jl b/scripts/build_docu.jl index f88f437c..f9398642 100644 --- a/scripts/build_docu.jl +++ b/scripts/build_docu.jl @@ -17,8 +17,4 @@ if !("LiveServer" in globaldependencies()) run(`julia -e 'using Pkg; Pkg.add("LiveServer")'`) end -if !("Documenter" ∈ keys(Pkg.project().dependencies)) - using TestEnv - TestEnv.activate() -end using LiveServer; servedocs(launch_browser=true) diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index 64a4a657..f3cb531a 100644 --- a/src/VortexStepMethod.jl +++ b/src/VortexStepMethod.jl @@ -37,6 +37,7 @@ export MVec3 export Model, VSM, LLT export AeroModel, LEI_AIRFOIL_BREUKELS, POLAR_VECTORS, POLAR_MATRICES, INVISCID export PanelDistribution, LINEAR, COSINE, COSINE_VAN_GARREL, SPLIT_PROVIDED, UNCHANGED +export PanelGroupingMethod, EQUAL_SIZE, REFINE export InitialGammaDistribution, ELLIPTIC, ZEROS export SolverStatus, FEASIBLE, INFEASIBLE, FAILURE export SolverType, LOOP, NONLIN @@ -138,6 +139,17 @@ Enumeration of the implemented panel distributions. UNCHANGED # Keep original sections end +""" + PanelGroupingMethod EQUAL_SIZE REFINE + +Enumeration of methods for grouping panels. + +# Elements +- EQUAL_SIZE: Divide panels into equally-sized sequential groups +- REFINE: Group refined panels back to their original unrefined section +""" +@enum PanelGroupingMethod EQUAL_SIZE REFINE + """ InitialGammaDistribution ELLIPTIC ZEROS diff --git a/src/obj_geometry.jl b/src/obj_geometry.jl index dec38cb7..1dae4cbc 100644 --- a/src/obj_geometry.jl +++ b/src/obj_geometry.jl @@ -366,7 +366,7 @@ function ObjWing( n_panels=56, n_sections=n_panels+1, n_groups=4, spanwise_distribution=UNCHANGED, spanwise_direction=[0.0, 1.0, 0.0], remove_nan=true, align_to_principal=false, alpha_range=deg2rad.(-5:1:20), delta_range=deg2rad.(-5:1:20), prn=true, - interp_steps=n_sections + interp_steps=n_sections, grouping_method::PanelGroupingMethod=EQUAL_SIZE ) !(n_groups == 0 || n_panels % n_groups == 0) && throw(ArgumentError("Number of panels should be divisible by number of groups")) !isapprox(spanwise_direction, [0.0, 1.0, 0.0]) && throw(ArgumentError("Spanwise direction has to be [0.0, 1.0, 0.0], not $spanwise_direction")) @@ -431,6 +431,7 @@ function ObjWing( Wing(n_panels, n_groups, spanwise_distribution, panel_props, MVec3(spanwise_direction), sections, refined_sections, remove_nan, + grouping_method, Int16[], non_deformed_sections, zeros(n_panels), zeros(n_panels), mass, gamma_tip, inertia_tensor, T_cad_body, R_cad_body, radius, le_interp, te_interp, area_interp, cache) diff --git a/src/solver.jl b/src/solver.jl index 9bce66ae..7f4daec3 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -302,14 +302,26 @@ function solve!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution= group_idx = 1 for wing in body_aero.wings if wing.n_groups > 0 - panels_per_group = wing.n_panels ÷ wing.n_groups - for _ in 1:wing.n_groups - for _ in 1:panels_per_group - group_moment_dist[group_idx] += moment_dist[panel_idx] - group_moment_coeff_dist[group_idx] += moment_coeff_dist[panel_idx] + if wing.grouping_method == EQUAL_SIZE + # Original method: divide panels into equally-sized sequential groups + panels_per_group = wing.n_panels ÷ wing.n_groups + for _ in 1:wing.n_groups + for _ in 1:panels_per_group + group_moment_dist[group_idx] += moment_dist[panel_idx] + group_moment_coeff_dist[group_idx] += moment_coeff_dist[panel_idx] + panel_idx += 1 + end + group_idx += 1 + end + elseif wing.grouping_method == REFINE + # REFINE method: group refined panels by their original unrefined section + for local_panel_idx in 1:wing.n_panels + original_section_idx = wing.refined_panel_mapping[local_panel_idx] + group_moment_dist[group_idx + original_section_idx - 1] += moment_dist[panel_idx] + group_moment_coeff_dist[group_idx + original_section_idx - 1] += moment_coeff_dist[panel_idx] panel_idx += 1 end - group_idx += 1 + group_idx += wing.n_groups end else # Skip panels for wings with n_groups=0 diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index 02f7a308..9536107f 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -229,6 +229,10 @@ mutable struct Wing <: AbstractWing refined_sections::Vector{Section} remove_nan::Bool + # Grouping + grouping_method::PanelGroupingMethod + refined_panel_mapping::Vector{Int16} # Maps each refined panel to its original unrefined section index + # Deformation fields non_deformed_sections::Vector{Section} theta_dist::Vector{Float64} @@ -252,7 +256,8 @@ end n_groups=n_panels, spanwise_distribution::PanelDistribution=LINEAR, spanwise_direction::PosVector=MVec3([0.0, 1.0, 0.0]), - remove_nan::Bool=true) + remove_nan::Bool=true, + grouping_method::PanelGroupingMethod=EQUAL_SIZE) Constructor for a [Wing](@ref) struct with default values that initializes the sections and refined sections as empty arrays. Creates a basic wing suitable for YAML-based construction. @@ -263,19 +268,28 @@ and refined sections as empty arrays. Creates a basic wing suitable for YAML-bas - `spanwise_distribution`::PanelDistribution = LINEAR: [PanelDistribution](@ref) - `spanwise_direction::MVec3` = MVec3([0.0, 1.0, 0.0]): Wing span direction vector - `remove_nan::Bool`: Wether to remove the NaNs from interpolations or not +- `grouping_method::PanelGroupingMethod` = EQUAL_SIZE: Method for grouping panels (EQUAL_SIZE or REFINE) """ function Wing(n_panels::Int; n_groups = n_panels, spanwise_distribution::PanelDistribution=LINEAR, spanwise_direction::PosVector=MVec3([0.0, 1.0, 0.0]), - remove_nan=true) - !(n_groups == 0 || n_panels % n_groups == 0) && throw(ArgumentError("Number of panels should be divisible by number of groups")) + remove_nan=true, + grouping_method::PanelGroupingMethod=EQUAL_SIZE) + # Validate grouping parameters + if grouping_method == EQUAL_SIZE + !(n_groups == 0 || n_panels % n_groups == 0) && throw(ArgumentError("With EQUAL_SIZE grouping, number of panels should be divisible by number of groups")) + end + # Note: For REFINE grouping, validation happens after refinement when we know the number of unrefined sections + panel_props = PanelProperties{n_panels}() # Initialize with default/empty values for optional fields Wing( n_panels, n_groups, spanwise_distribution, panel_props, spanwise_direction, Section[], Section[], remove_nan, + # Grouping + grouping_method, Int16[], # Deformation fields Section[], zeros(n_panels), zeros(n_panels), # Physical properties (defaults for non-OBJ wings) @@ -556,23 +570,25 @@ function refine_aerodynamic_mesh!(wing::AbstractWing) for i in eachindex(wing.sections) reinit!(wing.refined_sections[i], wing.sections[i]) end + compute_refined_panel_mapping!(wing) return nothing end @debug "Refining aerodynamic mesh from $(length(wing.sections)) sections to $n_sections sections." - + # Handle two-section case if n_sections == 2 reinit!(wing.refined_sections[1], LE[1,:], TE[1,:], aero_model[1], aero_data[1]) reinit!(wing.refined_sections[2], LE[end,:], TE[end,:], aero_model[end], aero_data[end]) + compute_refined_panel_mapping!(wing) return nothing end - + # Handle different distribution types if wing.spanwise_distribution == SPLIT_PROVIDED - return refine_mesh_by_splitting_provided_sections!(wing) + refine_mesh_by_splitting_provided_sections!(wing) elseif wing.spanwise_distribution in (LINEAR, COSINE, COSINE_VAN_GARREL) - return refine_mesh_for_linear_cosine_distribution!( + refine_mesh_for_linear_cosine_distribution!( wing, 1, wing.spanwise_distribution, @@ -585,6 +601,80 @@ function refine_aerodynamic_mesh!(wing::AbstractWing) else throw(ArgumentError("Unsupported spanwise panel distribution: $(wing.spanwise_distribution)")) end + + # Compute panel mapping by finding closest unrefined panel for each refined panel + compute_refined_panel_mapping!(wing) + + # Validate REFINE grouping method + if wing.grouping_method == REFINE && wing.n_groups > 0 + n_unrefined_panels = length(wing.sections) - 1 + if wing.n_groups != n_unrefined_panels + throw(ArgumentError( + "With REFINE grouping method, n_groups ($(wing.n_groups)) must equal " * + "the number of unrefined panels ($n_unrefined_panels). " * + "The wing has $(length(wing.sections)) unrefined sections, forming $n_unrefined_panels panels." + )) + end + end + + return nothing +end + + +""" + compute_refined_panel_mapping!(wing::AbstractWing) + +Compute the mapping from refined panels to unrefined panels by finding +the closest unrefined panel for each refined panel (based on panel center distance). +This is non-allocating and works after refinement is complete. +""" +function compute_refined_panel_mapping!(wing::AbstractWing) + n_unrefined_sections = length(wing.sections) + n_refined_panels = wing.n_panels + + # Handle case where no refinement occurred + if n_unrefined_sections == n_refined_panels + 1 + wing.refined_panel_mapping = Int16[i for i in 1:n_refined_panels] + return nothing + end + + # Ensure mapping array is allocated + if length(wing.refined_panel_mapping) != n_refined_panels + wing.refined_panel_mapping = zeros(Int16, n_refined_panels) + end + + # Compute centers of unrefined panels + n_unrefined_panels = n_unrefined_sections - 1 + unrefined_centers = Vector{MVec3}(undef, n_unrefined_panels) + for i in 1:n_unrefined_panels + le_mid = (wing.sections[i].LE_point + wing.sections[i+1].LE_point) / 2 + te_mid = (wing.sections[i].TE_point + wing.sections[i+1].TE_point) / 2 + unrefined_centers[i] = MVec3((le_mid + te_mid) / 2) + end + + # For each refined panel, find closest unrefined panel + for refined_panel_idx in 1:n_refined_panels + le_mid = (wing.refined_sections[refined_panel_idx].LE_point + + wing.refined_sections[refined_panel_idx+1].LE_point) / 2 + te_mid = (wing.refined_sections[refined_panel_idx].TE_point + + wing.refined_sections[refined_panel_idx+1].TE_point) / 2 + refined_center = MVec3((le_mid + te_mid) / 2) + + # Find closest unrefined panel + min_dist = Inf + closest_idx = 1 + for unrefined_panel_idx in 1:n_unrefined_panels + dist = norm(refined_center - unrefined_centers[unrefined_panel_idx]) + if dist < min_dist + min_dist = dist + closest_idx = unrefined_panel_idx + end + end + + wing.refined_panel_mapping[refined_panel_idx] = Int16(closest_idx) + end + + return nothing end @@ -736,7 +826,7 @@ function refine_mesh_for_linear_cosine_distribution!( target_length = target_lengths[i] # Find segment index - section_index = searchsortedlast(qc_cum_length, target_length) + section_index = searchsortedlast(qc_cum_length, target_length) section_index = clamp(section_index, 1, length(qc_cum_length)-1) # 4. Calculate weights @@ -747,7 +837,7 @@ function refine_mesh_for_linear_cosine_distribution!( right_weight = t # 5. Calculate quarter chord point - new_quarter_chord[i,:] = quarter_chord[section_index,:] + + new_quarter_chord[i,:] = quarter_chord[section_index,:] + t .* (quarter_chord[section_index+1,:] - quarter_chord[section_index,:]) # 6. Calculate chord vectors @@ -899,10 +989,10 @@ function refine_mesh_by_splitting_provided_sections!(wing::AbstractWing) # Add left section of pair reinit!(wing.refined_sections[idx], wing.sections[left_section_index]) idx += 1 - + # Calculate new sections for this pair num_new_sections = new_sections_per_pair + (left_section_index <= remaining ? 1 : 0) - + if num_new_sections > 0 # Prepare pair data LE_pair = hcat(LE[left_section_index], LE[left_section_index + 1])' @@ -915,7 +1005,7 @@ function refine_mesh_by_splitting_provided_sections!(wing::AbstractWing) aero_data[left_section_index], aero_data[left_section_index + 1] ] - + # Generate sections for this pair idx = refine_mesh_for_linear_cosine_distribution!( wing, diff --git a/src/yaml_geometry.jl b/src/yaml_geometry.jl index 3ba618f5..42a82199 100644 --- a/src/yaml_geometry.jl +++ b/src/yaml_geometry.jl @@ -188,7 +188,8 @@ function Wing( spanwise_distribution=LINEAR, spanwise_direction=[0.0, 1.0, 0.0], remove_nan=true, - prn=false + prn=false, + grouping_method::PanelGroupingMethod=EQUAL_SIZE ) !(n_groups == 0 || n_panels % n_groups == 0) && throw(ArgumentError("Number of panels should be divisible by number of groups")) !isapprox(spanwise_direction, [0.0, 1.0, 0.0]) && throw(ArgumentError("Spanwise direction has to be [0.0, 1.0, 0.0], not $spanwise_direction")) @@ -236,11 +237,12 @@ function Wing( end # Create Wing using the standard constructor - wing = Wing(n_panels; - n_groups=n_groups, + wing = Wing(n_panels; + n_groups=n_groups, spanwise_distribution=spanwise_distribution, - spanwise_direction=MVec3(spanwise_direction), - remove_nan=remove_nan + spanwise_direction=MVec3(spanwise_direction), + remove_nan=remove_nan, + grouping_method=grouping_method ) # Parse sections and populate wing diff --git a/test/Aqua.jl b/test/Aqua.jl index e16c0cc9..7a885062 100644 --- a/test/Aqua.jl +++ b/test/Aqua.jl @@ -1,8 +1,3 @@ -using Pkg -if ! ("Aqua" ∈ keys(Pkg.project().dependencies)) - using TestEnv; TestEnv.activate() -end - using Aqua, VortexStepMethod, Test @testset "Aqua.jl" begin Aqua.test_all( diff --git a/test/Project.toml b/test/Project.toml new file mode 100644 index 00000000..5354caa3 --- /dev/null +++ b/test/Project.toml @@ -0,0 +1,10 @@ +[deps] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +ControlPlots = "23c2ee80-7a9e-4350-b264-8e670f12517c" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +VortexStepMethod = "ed3cd733-9f0f-46a9-93e0-89b8d4998dd9" diff --git a/test/bench.jl b/test/bench.jl index 5acfb95f..88172d69 100644 --- a/test/bench.jl +++ b/test/bench.jl @@ -1,8 +1,3 @@ -using Pkg -if !("BenchmarkTools" ∈ keys(Pkg.project().dependencies)) - using TestEnv - TestEnv.activate() -end using BenchmarkTools using StaticArrays using VortexStepMethod diff --git a/test/bench_solve.jl b/test/bench_solve.jl index 1abe47d7..dc8c4590 100644 --- a/test/bench_solve.jl +++ b/test/bench_solve.jl @@ -6,12 +6,6 @@ using VortexStepMethod using BenchmarkTools using Test -using Pkg - -if !("CSV" ∈ keys(Pkg.project().dependencies)) - using TestEnv - TestEnv.activate() -end # Step 1: Define wing parameters n_panels = 20 # Number of panels diff --git a/test/settings/test_settings.jl b/test/settings/test_settings.jl index c7412894..33ccd7a7 100644 --- a/test/settings/test_settings.jl +++ b/test/settings/test_settings.jl @@ -1,8 +1,3 @@ -using Pkg -if ! ("Test" ∈ keys(Pkg.project().dependencies)) - using TestEnv; TestEnv.activate() -end - using VortexStepMethod using Test From c27d1572e3da954dfa43459bceec52c69a4e6070 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Wed, 29 Oct 2025 15:46:51 +0100 Subject: [PATCH 06/22] Working tests except aqua --- Project.toml | 22 ---- examples/stall_model.jl | 16 +-- test/Project.toml | 23 +++- test/wing_geometry/test_wing_geometry.jl | 153 +++++++++++++++++++++++ 4 files changed, 184 insertions(+), 30 deletions(-) diff --git a/Project.toml b/Project.toml index 6841b218..29000957 100644 --- a/Project.toml +++ b/Project.toml @@ -38,50 +38,28 @@ VortexStepMethodControlPlotsExt = "ControlPlots" VortexStepMethodMakieExt = "Makie" [compat] -Aqua = "0.8" -BenchmarkTools = "1" -CSV = "0.10" Colors = "0.13" -ControlPlots = "0.2.5" -DataFrames = "1.7" DefaultApplication = "1" DelimitedFiles = "1" DifferentiationInterface = "0.7.4" -Documenter = "1.8" FiniteDiff = "2.27.0" Interpolations = "0.15, 0.16" LaTeXStrings = "1" LinearAlgebra = "1" Logging = "1" -Makie = "0.24.6" Measures = "0.3" NonlinearSolve = "4.8.0" Parameters = "0.12" Pkg = "1" PreallocationTools = "0.4.31" PrecompileTools = "1.2.1" -Random = "1.10.0" RecursiveArrayTools = "3 - 3.36.0" SciMLBase = "2.77.0" Serialization = "1" StaticArrays = "1" Statistics = "1" StructMapping = "0.2.3" -Test = "1" Timers = "0.1" Xfoil = "1.1.0" YAML = "0.4.13" julia = "1.10, 1.11" - -[extras] -Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" -BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" -CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" -ControlPlots = "23c2ee80-7a9e-4350-b264-8e670f12517c" -DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" -Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" - -[targets] -test = ["Test", "DataFrames", "CSV", "Documenter", "BenchmarkTools", "ControlPlots", "Aqua", "Random"] diff --git a/examples/stall_model.jl b/examples/stall_model.jl index 89f3c0f4..78fb28fe 100644 --- a/examples/stall_model.jl +++ b/examples/stall_model.jl @@ -107,21 +107,23 @@ path_cfd_lebesque = joinpath( "V3_CL_CD_RANS_Lebesque_2024_Rey_300e4.csv" ) +# Only include literature data if file exists +literature_paths = isfile(path_cfd_lebesque) ? [path_cfd_lebesque] : String[] +labels = isfile(path_cfd_lebesque) ? + ["VSM CAD 19ribs", "VSM CAD 19ribs , with stall correction", "CFD_Lebesque Rey 30e5"] : + ["VSM CAD 19ribs", "VSM CAD 19ribs , with stall correction"] + PLOT && plot_polars( [vsm_solver, VSM_with_stall_correction], [body_aero, body_aero], - [ - "VSM CAD 19ribs", - "VSM CAD 19ribs , with stall correction", - "CFD_Lebesque Rey 30e5" - ]; - literature_path_list=[path_cfd_lebesque], + labels; + literature_path_list=literature_paths, angle_range=range(0, 25, length=25), angle_type="angle_of_attack", angle_of_attack=0, side_slip=0, v_a=10, - title="tutorial_testing_stall_model_n_panels_$(n_panels)_distribution_$(spanwise_distribution)", + title="tutorial_testing_stall_model_n_panels_$(n_panels)_distribution_$(spanwise_distribution)_grouping_$(CAD_wing.grouping_method)", data_type=".pdf", save_path=joinpath(save_folder, "polars"), is_save=true, diff --git a/test/Project.toml b/test/Project.toml index 5354caa3..32eda93f 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -5,6 +5,27 @@ CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" ControlPlots = "23c2ee80-7a9e-4350-b264-8e670f12517c" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -VortexStepMethod = "ed3cd733-9f0f-46a9-93e0-89b8d4998dd9" +YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6" + +[compat] +Aqua = "0.8" +BenchmarkTools = "1" +CSV = "0.10" +ControlPlots = "0.2.5" +DataFrames = "1.7" +Documenter = "1.8" +LinearAlgebra = "1" +Logging = "1" +Random = "1.10.0" +StaticArrays = "1" +Statistics = "1" +Test = "1" +YAML = "0.4.13" diff --git a/test/wing_geometry/test_wing_geometry.jl b/test/wing_geometry/test_wing_geometry.jl index d2d01f0f..ddb8ae25 100644 --- a/test/wing_geometry/test_wing_geometry.jl +++ b/test/wing_geometry/test_wing_geometry.jl @@ -332,4 +332,157 @@ end @test section.aero_model == INVISCID end end + + @testset "REFINE grouping panel mapping" begin + # Test that refined panel mapping actually maps each panel to its closest unrefined panel + + @testset "LINEAR distribution" begin + n_panels = 20 + span = 10.0 + + wing = Wing(n_panels; spanwise_distribution=LINEAR, n_groups=2, grouping_method=REFINE) + # 3 sections = 2 unrefined panels + add_section!(wing, [0.0, span/2, 0.0], [1.0, span/2, 0.0], INVISCID) + add_section!(wing, [0.0, 0.0, 0.0], [1.0, 0.0, 0.0], INVISCID) + add_section!(wing, [0.0, -span/2, 0.0], [1.0, -span/2, 0.0], INVISCID) + + refine_aerodynamic_mesh!(wing) + + @test length(wing.refined_panel_mapping) == n_panels + + # Manually verify each refined panel is mapped to its closest unrefined panel + n_unrefined_panels = length(wing.sections) - 1 + for refined_panel_idx in 1:n_panels + # Calculate refined panel center + le_mid = (wing.refined_sections[refined_panel_idx].LE_point + + wing.refined_sections[refined_panel_idx+1].LE_point) / 2 + te_mid = (wing.refined_sections[refined_panel_idx].TE_point + + wing.refined_sections[refined_panel_idx+1].TE_point) / 2 + refined_center = (le_mid + te_mid) / 2 + + # Find closest unrefined panel manually + min_dist = Inf + closest_idx = 1 + for unrefined_panel_idx in 1:n_unrefined_panels + le_mid_unref = (wing.sections[unrefined_panel_idx].LE_point + + wing.sections[unrefined_panel_idx+1].LE_point) / 2 + te_mid_unref = (wing.sections[unrefined_panel_idx].TE_point + + wing.sections[unrefined_panel_idx+1].TE_point) / 2 + unrefined_center = (le_mid_unref + te_mid_unref) / 2 + + dist = norm(refined_center - unrefined_center) + if dist < min_dist + min_dist = dist + closest_idx = unrefined_panel_idx + end + end + + # Verify mapping is correct + @test wing.refined_panel_mapping[refined_panel_idx] == closest_idx + end + end + + @testset "COSINE distribution" begin + n_panels = 30 + span = 20.0 + + wing = Wing(n_panels; spanwise_distribution=COSINE, n_groups=3, grouping_method=REFINE) + # 4 sections = 3 unrefined panels + add_section!(wing, [0.0, span/2, 0.0], [1.0, span/2, 0.0], INVISCID) + add_section!(wing, [0.0, span/6, 0.0], [1.0, span/6, 0.0], INVISCID) + add_section!(wing, [0.0, -span/6, 0.0], [1.0, -span/6, 0.0], INVISCID) + add_section!(wing, [0.0, -span/2, 0.0], [1.0, -span/2, 0.0], INVISCID) + + refine_aerodynamic_mesh!(wing) + + @test length(wing.refined_panel_mapping) == n_panels + + # Verify each panel is mapped to its closest unrefined panel + n_unrefined_panels = length(wing.sections) - 1 + for refined_panel_idx in 1:n_panels + # Calculate refined panel center + le_mid = (wing.refined_sections[refined_panel_idx].LE_point + + wing.refined_sections[refined_panel_idx+1].LE_point) / 2 + te_mid = (wing.refined_sections[refined_panel_idx].TE_point + + wing.refined_sections[refined_panel_idx+1].TE_point) / 2 + refined_center = (le_mid + te_mid) / 2 + + # Find closest unrefined panel manually + min_dist = Inf + closest_idx = 1 + for unrefined_panel_idx in 1:n_unrefined_panels + le_mid_unref = (wing.sections[unrefined_panel_idx].LE_point + + wing.sections[unrefined_panel_idx+1].LE_point) / 2 + te_mid_unref = (wing.sections[unrefined_panel_idx].TE_point + + wing.sections[unrefined_panel_idx+1].TE_point) / 2 + unrefined_center = (le_mid_unref + te_mid_unref) / 2 + + dist = norm(refined_center - unrefined_center) + if dist < min_dist + min_dist = dist + closest_idx = unrefined_panel_idx + end + end + + # Verify mapping is correct + @test wing.refined_panel_mapping[refined_panel_idx] == closest_idx + end + end + + @testset "SPLIT_PROVIDED distribution" begin + n_panels = 12 + + wing = Wing(n_panels; spanwise_distribution=SPLIT_PROVIDED, n_groups=3, grouping_method=REFINE) + # 4 sections = 3 unrefined panels + add_section!(wing, [0.0, 6.0, 0.0], [1.0, 6.0, 0.0], INVISCID) + add_section!(wing, [0.0, 2.0, 0.0], [1.0, 2.0, 0.0], INVISCID) + add_section!(wing, [0.0, -2.0, 0.0], [1.0, -2.0, 0.0], INVISCID) + add_section!(wing, [0.0, -6.0, 0.0], [1.0, -6.0, 0.0], INVISCID) + + refine_aerodynamic_mesh!(wing) + + @test length(wing.refined_panel_mapping) == n_panels + + # Verify each panel is mapped to its closest unrefined panel + n_unrefined_panels = length(wing.sections) - 1 + for refined_panel_idx in 1:n_panels + # Calculate refined panel center + le_mid = (wing.refined_sections[refined_panel_idx].LE_point + + wing.refined_sections[refined_panel_idx+1].LE_point) / 2 + te_mid = (wing.refined_sections[refined_panel_idx].TE_point + + wing.refined_sections[refined_panel_idx+1].TE_point) / 2 + refined_center = (le_mid + te_mid) / 2 + + # Find closest unrefined panel manually + min_dist = Inf + closest_idx = 1 + for unrefined_panel_idx in 1:n_unrefined_panels + le_mid_unref = (wing.sections[unrefined_panel_idx].LE_point + + wing.sections[unrefined_panel_idx+1].LE_point) / 2 + te_mid_unref = (wing.sections[unrefined_panel_idx].TE_point + + wing.sections[unrefined_panel_idx+1].TE_point) / 2 + unrefined_center = (le_mid_unref + te_mid_unref) / 2 + + dist = norm(refined_center - unrefined_center) + if dist < min_dist + min_dist = dist + closest_idx = unrefined_panel_idx + end + end + + # Verify mapping is correct + @test wing.refined_panel_mapping[refined_panel_idx] == closest_idx + end + end + + @testset "Validation: n_groups must equal unrefined panels" begin + wing = Wing(20; spanwise_distribution=LINEAR, n_groups=5, grouping_method=REFINE) + add_section!(wing, [0.0, 5.0, 0.0], [1.0, 5.0, 0.0], INVISCID) + add_section!(wing, [0.0, 0.0, 0.0], [1.0, 0.0, 0.0], INVISCID) + add_section!(wing, [0.0, -5.0, 0.0], [1.0, -5.0, 0.0], INVISCID) + + # Should throw error: 5 groups but only 2 unrefined panels + @test_throws ArgumentError refine_aerodynamic_mesh!(wing) + end + end end \ No newline at end of file From 35356430c0ef65814501f2fba2be71ea1a43a595 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Wed, 29 Oct 2025 15:57:35 +0100 Subject: [PATCH 07/22] Working tests with aqua --- Project.toml | 2 ++ test/Project.toml | 2 ++ 2 files changed, 4 insertions(+) diff --git a/Project.toml b/Project.toml index 29000957..7609394f 100644 --- a/Project.toml +++ b/Project.toml @@ -39,6 +39,7 @@ VortexStepMethodMakieExt = "Makie" [compat] Colors = "0.13" +ControlPlots = "0.2.5" DefaultApplication = "1" DelimitedFiles = "1" DifferentiationInterface = "0.7.4" @@ -47,6 +48,7 @@ Interpolations = "0.15, 0.16" LaTeXStrings = "1" LinearAlgebra = "1" Logging = "1" +Makie = "0.24.6" Measures = "0.3" NonlinearSolve = "4.8.0" Parameters = "0.12" diff --git a/test/Project.toml b/test/Project.toml index 32eda93f..9ede8c7a 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -22,9 +22,11 @@ CSV = "0.10" ControlPlots = "0.2.5" DataFrames = "1.7" Documenter = "1.8" +Interpolations = "0.15, 0.16" LinearAlgebra = "1" Logging = "1" Random = "1.10.0" +Serialization = "1" StaticArrays = "1" Statistics = "1" Test = "1" From 2b0301504a3ac6a2c14054df13f28bc6ee5713e8 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sat, 1 Nov 2025 23:31:49 +0100 Subject: [PATCH 08/22] Don't deform tuple --- examples/Project.toml | 3 + src/wing_geometry.jl | 8 +- test/yaml_geometry/test_yaml_geometry.jl | 1 + .../test_yaml_wing_deformation.jl | 175 ++++++++++++++++++ 4 files changed, 184 insertions(+), 3 deletions(-) create mode 100644 examples/Project.toml create mode 100644 test/yaml_geometry/test_yaml_wing_deformation.jl diff --git a/examples/Project.toml b/examples/Project.toml new file mode 100644 index 00000000..24966cd6 --- /dev/null +++ b/examples/Project.toml @@ -0,0 +1,3 @@ +[deps] +ControlPlots = "23c2ee80-7a9e-4350-b264-8e670f12517c" +VortexStepMethod = "ed3cd733-9f0f-46a9-93e0-89b8d4998dd9" diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index 02f7a308..4f00d2d4 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -44,10 +44,12 @@ function reinit!(section::Section, LE_point, TE_point, aero_model=nothing, aero_ section.TE_point .= TE_point (!isnothing(aero_model)) && (section.aero_model = aero_model) if !isnothing(aero_data) - if !isnothing(section.aero_data) - section.aero_data .= aero_data - else + # NTuple is immutable, so we must assign directly + # For mutable types (Vector, Matrix), we can broadcast for efficiency + if aero_data isa NTuple || isnothing(section.aero_data) section.aero_data = aero_data + else + section.aero_data .= aero_data end end nothing diff --git a/test/yaml_geometry/test_yaml_geometry.jl b/test/yaml_geometry/test_yaml_geometry.jl index 12daa54f..58924a91 100644 --- a/test/yaml_geometry/test_yaml_geometry.jl +++ b/test/yaml_geometry/test_yaml_geometry.jl @@ -4,4 +4,5 @@ using Test # Include specific test files for better organization include("test_load_polar_data.jl") include("test_wing_constructor.jl") + include("test_yaml_wing_deformation.jl") end diff --git a/test/yaml_geometry/test_yaml_wing_deformation.jl b/test/yaml_geometry/test_yaml_wing_deformation.jl new file mode 100644 index 00000000..da931177 --- /dev/null +++ b/test/yaml_geometry/test_yaml_wing_deformation.jl @@ -0,0 +1,175 @@ +using VortexStepMethod +using LinearAlgebra +using Test + +@testset "YAML Wing Deformation Tests" begin + @testset "Simple Wing Deformation" begin + # Load existing simple_wing.yaml + simple_wing_file = test_data_path("yaml_geometry", "simple_wing.yaml") + wing = Wing(simple_wing_file; n_panels=4, n_groups=2) + body_aero = BodyAerodynamics([wing]) + + # Store original TE point for comparison + i = length(body_aero.panels) ÷ 2 + original_te_point = copy(body_aero.panels[i].TE_point_1) + original_le_point = copy(body_aero.panels[i].LE_point_1) + + # Apply deformation with non-zero angles + theta_dist = fill(deg2rad(30.0), wing.n_panels) # 30 degrees twist + delta_dist = fill(deg2rad(5.0), wing.n_panels) # 5 degrees trailing edge deflection + + VortexStepMethod.deform!(wing, theta_dist, delta_dist) + VortexStepMethod.reinit!(body_aero) + + # Check if TE point changed after deformation + deformed_te_point = copy(body_aero.panels[i].TE_point_1) + deformed_le_point = copy(body_aero.panels[i].LE_point_1) + + # TE point should change significantly due to twist and deflection + @test !isapprox(original_te_point, deformed_te_point, atol=1e-2) + @test deformed_te_point[3] < original_te_point[3] # TE should move down with positive twist + + # LE point should also change due to twist + @test !isapprox(original_le_point, deformed_le_point, atol=1e-2) + + # Check delta is set correctly + @test body_aero.panels[i].delta ≈ deg2rad(5.0) + + # Reset deformation with zero angles + zero_theta_dist = zeros(wing.n_panels) + zero_delta_dist = zeros(wing.n_panels) + + VortexStepMethod.deform!(wing, zero_theta_dist, zero_delta_dist) + VortexStepMethod.reinit!(body_aero) + + # Check if TE point returned to original position + reset_te_point = copy(body_aero.panels[i].TE_point_1) + reset_le_point = copy(body_aero.panels[i].LE_point_1) + @test original_te_point ≈ reset_te_point atol=1e-4 + @test original_le_point ≈ reset_le_point atol=1e-4 + @test body_aero.panels[i].delta ≈ 0.0 atol=1e-4 + end + + @testset "Complex Wing Deformation" begin + # Load existing complex_wing.yaml with multiple sections + complex_wing_file = test_data_path("yaml_geometry", "complex_wing.yaml") + wing = Wing(complex_wing_file; n_panels=12, n_groups=3) + body_aero = BodyAerodynamics([wing]) + + # Store original points for multiple panels + original_points = [] + test_indices = [1, length(body_aero.panels) ÷ 2, length(body_aero.panels)] + for i in test_indices + push!(original_points, ( + LE=copy(body_aero.panels[i].LE_point_1), + TE=copy(body_aero.panels[i].TE_point_1) + )) + end + + # Apply spanwise-varying deformation + theta_dist = [deg2rad(10.0 * i / wing.n_panels) for i in 1:wing.n_panels] # Linear twist distribution + delta_dist = [deg2rad(-5.0 + 10.0 * i / wing.n_panels) for i in 1:wing.n_panels] # Varying deflection + + VortexStepMethod.deform!(wing, theta_dist, delta_dist) + VortexStepMethod.reinit!(body_aero) + + # Check that different panels have different deformations + for (idx, i) in enumerate(test_indices) + deformed_te = body_aero.panels[i].TE_point_1 + deformed_le = body_aero.panels[i].LE_point_1 + + # Points should have changed + @test !isapprox(original_points[idx].TE, deformed_te, atol=1e-2) + @test !isapprox(original_points[idx].LE, deformed_le, atol=1e-2) + end + + # Check that the deformation is applied correctly + # First panel should have smaller theta, last panel should have larger theta + @test body_aero.panels[1].delta < body_aero.panels[end].delta + + # Reset and verify + VortexStepMethod.deform!(wing, zeros(wing.n_panels), zeros(wing.n_panels)) + VortexStepMethod.reinit!(body_aero) + + for (idx, i) in enumerate(test_indices) + reset_te = body_aero.panels[i].TE_point_1 + reset_le = body_aero.panels[i].LE_point_1 + @test original_points[idx].TE ≈ reset_te atol=1e-4 + @test original_points[idx].LE ≈ reset_le atol=1e-4 + @test body_aero.panels[i].delta ≈ 0.0 atol=1e-4 + end + end + + @testset "Multiple Reinit Calls with NTuple aero_data" begin + # This test specifically checks the NTuple handling fix + simple_wing_file = test_data_path("yaml_geometry", "simple_wing.yaml") + wing = Wing(simple_wing_file; n_panels=4, n_groups=2) + + # Verify that sections have NTuple aero_data (for wings with simple polars) + # or other valid AeroData types + @test wing.sections[1].aero_data !== nothing + + # Perform multiple reinit! calls to ensure NTuple handling works + for _ in 1:5 + VortexStepMethod.reinit!(wing) + end + + # Wing should still be valid after multiple reinits + @test wing.sections[1].aero_data !== nothing + @test length(wing.sections) == 2 + end + + @testset "Deformation with BodyAerodynamics Reinit" begin + # Test that reinit! on BodyAerodynamics properly handles deformed wings + simple_wing_file = test_data_path("yaml_geometry", "simple_wing.yaml") + wing = Wing(simple_wing_file; n_panels=4, n_groups=2) + body_aero = BodyAerodynamics([wing]) + + # Apply deformation + theta_dist = fill(deg2rad(15.0), wing.n_panels) + delta_dist = fill(deg2rad(3.0), wing.n_panels) + VortexStepMethod.deform!(wing, theta_dist, delta_dist) + + # Store state after deformation + i = length(body_aero.panels) ÷ 2 + + # Multiple reinit calls should work without errors + for _ in 1:3 + VortexStepMethod.reinit!(body_aero; + va=zeros(3), + omega=zeros(3), + init_aero=true + ) + end + + # Panel should maintain deformation + @test body_aero.panels[i].delta ≈ deg2rad(3.0) atol=1e-6 + end + + @testset "Edge Cases" begin + simple_wing_file = test_data_path("yaml_geometry", "simple_wing.yaml") + wing = Wing(simple_wing_file; n_panels=2, n_groups=1) + body_aero = BodyAerodynamics([wing]) + + # Test zero deformation + VortexStepMethod.deform!(wing, zeros(wing.n_panels), zeros(wing.n_panels)) + VortexStepMethod.reinit!(body_aero) + @test all(p.delta ≈ 0.0 for p in body_aero.panels) + + # Test large deformation angles + theta_dist = fill(deg2rad(60.0), wing.n_panels) + delta_dist = fill(deg2rad(30.0), wing.n_panels) + + # Should not error even with large angles + VortexStepMethod.deform!(wing, theta_dist, delta_dist) + VortexStepMethod.reinit!(body_aero) + @test all(p.delta ≈ deg2rad(30.0) for p in body_aero.panels) + + # Test negative angles + theta_dist = fill(deg2rad(-20.0), wing.n_panels) + delta_dist = fill(deg2rad(-10.0), wing.n_panels) + VortexStepMethod.deform!(wing, theta_dist, delta_dist) + VortexStepMethod.reinit!(body_aero) + @test all(p.delta ≈ deg2rad(-10.0) for p in body_aero.panels) + end +end From 4ca52ac9ead629999d2260ce12f639ba9fd576db Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Thu, 6 Nov 2025 15:20:03 +0100 Subject: [PATCH 09/22] Option to not use data prefix --- src/settings.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/settings.jl b/src/settings.jl index fe538ed9..3448802d 100644 --- a/src/settings.jl +++ b/src/settings.jl @@ -40,9 +40,13 @@ end solver_settings::SolverSettings = SolverSettings() end -function VSMSettings(filename) +function VSMSettings(filename; data_prefix=true) # Uwe's suggested 3-line approach using StructMapping.jl (adapted) - data = YAML.load_file(joinpath("data", filename)) + if data_prefix + data = YAML.load_file(joinpath("data", filename)) + else + data = YAML.load_file(filename) + end # Use StructMapping for basic structure conversion # But handle special fields manually due to enum conversion needs @@ -110,4 +114,4 @@ function Base.show(io::IO, vsm_settings::VSMSettings) print(io, replace(repr(wing), "\n" => "\n ")) end print(io, replace(repr(vsm_settings.solver_settings), "\n" => "\n ")) -end \ No newline at end of file +end From 128f720fae4b8584642a08577a36cb58d6bf103b Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Fri, 7 Nov 2025 11:11:58 +0100 Subject: [PATCH 10/22] Add obj set --- src/settings.jl | 12 +++- src/yaml_geometry.jl | 71 ++++++++++++++++--- test/settings/test_settings.jl | 2 +- .../test_yaml_wing_deformation.jl | 12 ++-- 4 files changed, 80 insertions(+), 17 deletions(-) diff --git a/src/settings.jl b/src/settings.jl index 3448802d..ad3fb75f 100644 --- a/src/settings.jl +++ b/src/settings.jl @@ -8,8 +8,10 @@ end @with_kw mutable struct WingSettings name::String = "main_wing" geometry_file::String = "" # path to wing geometry YAML file + obj_file::String = "" # path to .obj geometry file + dat_file::String = "" # path to .dat airfoil file n_panels::Int64 = 40 - n_groups::Int64 = 40 + n_groups::Int64 = 40 spanwise_panel_distribution::PanelDistribution = LINEAR spanwise_direction::MVec3 = [0.0, 1.0, 0.0] remove_nan = true @@ -68,12 +70,18 @@ function VSMSettings(filename; data_prefix=true) if haskey(wing_data, "geometry_file") wing.geometry_file = wing_data["geometry_file"] end + if haskey(wing_data, "obj_file") + wing.obj_file = wing_data["obj_file"] + end + if haskey(wing_data, "dat_file") + wing.dat_file = wing_data["dat_file"] + end wing.n_panels = wing_data["n_panels"] wing.n_groups = wing_data["n_groups"] wing.spanwise_panel_distribution = eval(Symbol(wing_data["spanwise_panel_distribution"])) wing.spanwise_direction = MVec3(wing_data["spanwise_direction"]) wing.remove_nan = wing_data["remove_nan"] - + push!(vsm_settings.wings, wing) n_panels += wing.n_panels n_groups += wing.n_groups diff --git a/src/yaml_geometry.jl b/src/yaml_geometry.jl index 3ba618f5..8279d525 100644 --- a/src/yaml_geometry.jl +++ b/src/yaml_geometry.jl @@ -275,9 +275,13 @@ end Create a wing model from VSM settings configuration. -This constructor is a convenience wrapper that extracts wing configuration -from VSMSettings and creates a Wing using the YAML geometry file path and -parameters specified in the settings. +This constructor is a convenience wrapper that extracts wing configuration +from VSMSettings and creates a Wing using either: +- YAML geometry file (geometry_file field), or +- OBJ + DAT files (obj_file and dat_file fields) + +The constructor automatically determines which path to use based on which +fields are populated in the settings. # Arguments - `settings`: VSMSettings object containing wing configuration @@ -287,15 +291,64 @@ A fully initialized `Wing` instance ready for aerodynamic simulation. # Example ```julia -# Load settings and create wing in one step +# Using YAML geometry settings = VSMSettings("path/to/vsm_settings.yaml") wing = Wing(settings) + +# Settings can specify either: +# - geometry_file: "path/to/wing.yaml" # YAML-based +# - obj_file + dat_file # OBJ-based ``` """ function Wing(settings::VSMSettings) - Wing(settings.wings[1].geometry_file; - n_panels=settings.wings[1].n_panels, - n_groups=settings.wings[1].n_groups, - spanwise_distribution=settings.wings[1].spanwise_panel_distribution - ) + wing_settings = settings.wings[1] + + # Check which geometry format to use + has_yaml = !isempty(wing_settings.geometry_file) + has_obj = !isempty(wing_settings.obj_file) + has_dat = !isempty(wing_settings.dat_file) + + if has_yaml && (has_obj || has_dat) + throw(ArgumentError( + "Cannot specify both geometry_file and obj_file/dat_file" + )) + end + + if has_obj && !has_dat + throw(ArgumentError( + "obj_file requires dat_file to be specified" + )) + end + + if has_dat && !has_obj + throw(ArgumentError( + "dat_file requires obj_file to be specified" + )) + end + + if has_yaml + # Use YAML geometry constructor + Wing(wing_settings.geometry_file; + n_panels=wing_settings.n_panels, + n_groups=wing_settings.n_groups, + spanwise_distribution=wing_settings.spanwise_panel_distribution, + remove_nan=wing_settings.remove_nan + ) + elseif has_obj && has_dat + # Use ObjWing constructor + ObjWing( + wing_settings.obj_file, + wing_settings.dat_file; + n_panels=wing_settings.n_panels, + n_groups=wing_settings.n_groups, + spanwise_distribution=wing_settings.spanwise_panel_distribution, + spanwise_direction=wing_settings.spanwise_direction, + remove_nan=wing_settings.remove_nan + ) + else + throw(ArgumentError( + "WingSettings must specify either geometry_file or " * + "both obj_file and dat_file" + )) + end end diff --git a/test/settings/test_settings.jl b/test/settings/test_settings.jl index c7412894..80126bcf 100644 --- a/test/settings/test_settings.jl +++ b/test/settings/test_settings.jl @@ -16,6 +16,6 @@ using Test @test vss.wings isa Vector{WingSettings} @test length(vss.wings) == 2 io = IOBuffer(repr(vss)) - @test countlines(io) == 40 # Updated to match new output format + @test countlines(io) == 44 # Updated to match new output format end nothing diff --git a/test/yaml_geometry/test_yaml_wing_deformation.jl b/test/yaml_geometry/test_yaml_wing_deformation.jl index da931177..2631f2cc 100644 --- a/test/yaml_geometry/test_yaml_wing_deformation.jl +++ b/test/yaml_geometry/test_yaml_wing_deformation.jl @@ -29,8 +29,8 @@ using Test @test !isapprox(original_te_point, deformed_te_point, atol=1e-2) @test deformed_te_point[3] < original_te_point[3] # TE should move down with positive twist - # LE point should also change due to twist - @test !isapprox(original_le_point, deformed_le_point, atol=1e-2) + # LE point stays fixed (deform! rotates TE around LE) + @test isapprox(original_le_point, deformed_le_point, atol=1e-4) # Check delta is set correctly @test body_aero.panels[i].delta ≈ deg2rad(5.0) @@ -78,9 +78,10 @@ using Test deformed_te = body_aero.panels[i].TE_point_1 deformed_le = body_aero.panels[i].LE_point_1 - # Points should have changed + # TE points should have changed due to deformation @test !isapprox(original_points[idx].TE, deformed_te, atol=1e-2) - @test !isapprox(original_points[idx].LE, deformed_le, atol=1e-2) + # LE points stay fixed (deform! rotates TE around LE) + @test isapprox(original_points[idx].LE, deformed_le, atol=1e-4) end # Check that the deformation is applied correctly @@ -116,7 +117,8 @@ using Test # Wing should still be valid after multiple reinits @test wing.sections[1].aero_data !== nothing - @test length(wing.sections) == 2 + # Note: For deformation support, wing.sections now points to refined_sections + @test length(wing.sections) == wing.n_panels + 1 end @testset "Deformation with BodyAerodynamics Reinit" begin From 7a7e2038095fe8daeb4ca8d11dbaf0e03b61a19a Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Fri, 7 Nov 2025 13:09:45 +0100 Subject: [PATCH 11/22] Old tests not failing --- src/wing_geometry.jl | 32 +++++++++++++++++-- .../test_yaml_wing_deformation.jl | 11 +++++-- 2 files changed, 39 insertions(+), 4 deletions(-) diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index a1efc422..4ce34568 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -450,7 +450,7 @@ function deform!(wing::Wing) local_y .= normalize(section1.LE_point - section2.LE_point) chord .= section1.TE_point .- section1.LE_point normal .= chord × local_y - @. wing.sections[i].TE_point = section1.LE_point + cos(wing.theta_dist[i]) * chord - sin(wing.theta_dist[i]) * normal + @. wing.refined_sections[i].TE_point = section1.LE_point + cos(wing.theta_dist[i]) * chord - sin(wing.theta_dist[i]) * normal end return nothing end @@ -528,6 +528,28 @@ function flip_created_coord_in_pairs_if_needed!(coord::Matrix{Float64}) end +""" + update_non_deformed_sections!(wing::AbstractWing) + +Create non_deformed_sections to match refined_sections. +This enables deformation support for all wings (YAML and OBJ). +Should be called after refined_sections are populated for the FIRST time only. +Once populated, non_deformed_sections serves as the undeformed reference geometry. +""" +function update_non_deformed_sections!(wing::AbstractWing) + n_sections = wing.n_panels + 1 + + # Only populate non_deformed_sections if it's empty (initial setup) + # Once populated, it serves as the undeformed reference and should not be overwritten + if isempty(wing.non_deformed_sections) + wing.non_deformed_sections = [Section() for _ in 1:n_sections] + for i in 1:n_sections + reinit!(wing.non_deformed_sections[i], wing.refined_sections[i]) + end + end + return nothing +end + """ refine_aerodynamic_mesh!(wing::AbstractWing) @@ -542,6 +564,7 @@ function refine_aerodynamic_mesh!(wing::AbstractWing) if length(wing.refined_sections) == 0 if wing.spanwise_distribution == UNCHANGED || length(wing.sections) == n_sections wing.refined_sections = wing.sections + update_non_deformed_sections!(wing) return nothing else wing.refined_sections = Section[Section() for _ in 1:wing.n_panels+1] @@ -573,6 +596,7 @@ function refine_aerodynamic_mesh!(wing::AbstractWing) reinit!(wing.refined_sections[i], wing.sections[i]) end compute_refined_panel_mapping!(wing) + update_non_deformed_sections!(wing) return nothing end @@ -583,6 +607,7 @@ function refine_aerodynamic_mesh!(wing::AbstractWing) reinit!(wing.refined_sections[1], LE[1,:], TE[1,:], aero_model[1], aero_data[1]) reinit!(wing.refined_sections[2], LE[end,:], TE[end,:], aero_model[end], aero_data[end]) compute_refined_panel_mapping!(wing) + update_non_deformed_sections!(wing) return nothing end @@ -619,6 +644,9 @@ function refine_aerodynamic_mesh!(wing::AbstractWing) end end + # Create/update non_deformed_sections to match refined_sections + update_non_deformed_sections!(wing) + return nothing end @@ -1116,4 +1144,4 @@ function Base.getproperty(w::AbstractWing, s::Symbol) else return getfield(w, s) end -end \ No newline at end of file +end diff --git a/test/yaml_geometry/test_yaml_wing_deformation.jl b/test/yaml_geometry/test_yaml_wing_deformation.jl index 2631f2cc..75d64642 100644 --- a/test/yaml_geometry/test_yaml_wing_deformation.jl +++ b/test/yaml_geometry/test_yaml_wing_deformation.jl @@ -2,6 +2,12 @@ using VortexStepMethod using LinearAlgebra using Test +# Load TestSupport if not already loaded (for standalone execution) +if !@isdefined(TestSupport) + include(joinpath(@__DIR__, "..", "TestSupport.jl")) + using .TestSupport +end + @testset "YAML Wing Deformation Tests" begin @testset "Simple Wing Deformation" begin # Load existing simple_wing.yaml @@ -117,8 +123,9 @@ using Test # Wing should still be valid after multiple reinits @test wing.sections[1].aero_data !== nothing - # Note: For deformation support, wing.sections now points to refined_sections - @test length(wing.sections) == wing.n_panels + 1 + # Verify refined_sections and non_deformed_sections are properly populated + @test length(wing.refined_sections) == wing.n_panels + 1 + @test length(wing.non_deformed_sections) == wing.n_panels + 1 end @testset "Deformation with BodyAerodynamics Reinit" begin From 2ef5c3ed581cbeec08ca1e614c71d0386e61e1ba Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Fri, 7 Nov 2025 15:03:13 +0100 Subject: [PATCH 12/22] Passing tests --- src/body_aerodynamics.jl | 9 ++++++--- src/wing_geometry.jl | 16 ++++++++++++++-- .../test_yaml_wing_deformation.jl | 18 ++++++++++-------- 3 files changed, 30 insertions(+), 13 deletions(-) diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 9e09f54c..103745dc 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -114,7 +114,7 @@ function Base.setproperty!(obj::BodyAerodynamics, sym::Symbol, val) end """ - reinit!(body_aero::BodyAerodynamics; init_aero, va, omega) + reinit!(body_aero::BodyAerodynamics; init_aero, va, omega, refine_mesh) Initialize a BodyAerodynamics struct in-place by setting up panels and coefficients. @@ -125,6 +125,8 @@ Initialize a BodyAerodynamics struct in-place by setting up panels and coefficie - `init_aero::Bool`: Wether to initialize the aero data or not - `va=[15.0, 0.0, 0.0]`: Apparent wind vector - `omega=zeros(3)`: Turn rate in kite body frame x y and z +- `refine_mesh=true`: Whether to refine wing meshes. Set to `false` after + `deform!()` to preserve deformed geometry. # Returns nothing @@ -132,12 +134,13 @@ nothing function reinit!(body_aero::BodyAerodynamics; init_aero=true, va=[15.0, 0.0, 0.0], - omega=zeros(MVec3) + omega=zeros(MVec3), + refine_mesh=true ) idx = 1 vec = @MVector zeros(3) for wing in body_aero.wings - reinit!(wing) + reinit!(wing; refine_mesh=refine_mesh) panel_props = wing.panel_props # Create panels diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index 4ce34568..79554cb4 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -301,8 +301,20 @@ function Wing(n_panels::Int; ) end -function reinit!(wing::AbstractWing) - refine_aerodynamic_mesh!(wing) +""" + reinit!(wing::AbstractWing; refine_mesh=true) + +Reinitialize wing geometry and panel properties. + +# Keyword Arguments +- `refine_mesh::Bool=true`: Whether to refine the mesh. Set to `false` after + `deform!()` to preserve deformed geometry while updating panel properties. +""" +function reinit!(wing::AbstractWing; refine_mesh=true) + # Refine mesh unless explicitly disabled (e.g., to preserve deformation) + if refine_mesh + refine_aerodynamic_mesh!(wing) + end # Calculate panel properties update_panel_properties!( diff --git a/test/yaml_geometry/test_yaml_wing_deformation.jl b/test/yaml_geometry/test_yaml_wing_deformation.jl index 75d64642..580e3dc9 100644 --- a/test/yaml_geometry/test_yaml_wing_deformation.jl +++ b/test/yaml_geometry/test_yaml_wing_deformation.jl @@ -25,7 +25,7 @@ end delta_dist = fill(deg2rad(5.0), wing.n_panels) # 5 degrees trailing edge deflection VortexStepMethod.deform!(wing, theta_dist, delta_dist) - VortexStepMethod.reinit!(body_aero) + VortexStepMethod.reinit!(body_aero; refine_mesh=false) # Check if TE point changed after deformation deformed_te_point = copy(body_aero.panels[i].TE_point_1) @@ -46,7 +46,7 @@ end zero_delta_dist = zeros(wing.n_panels) VortexStepMethod.deform!(wing, zero_theta_dist, zero_delta_dist) - VortexStepMethod.reinit!(body_aero) + VortexStepMethod.reinit!(body_aero; refine_mesh=false) # Check if TE point returned to original position reset_te_point = copy(body_aero.panels[i].TE_point_1) @@ -77,7 +77,7 @@ end delta_dist = [deg2rad(-5.0 + 10.0 * i / wing.n_panels) for i in 1:wing.n_panels] # Varying deflection VortexStepMethod.deform!(wing, theta_dist, delta_dist) - VortexStepMethod.reinit!(body_aero) + VortexStepMethod.reinit!(body_aero; refine_mesh=false) # Check that different panels have different deformations for (idx, i) in enumerate(test_indices) @@ -96,7 +96,7 @@ end # Reset and verify VortexStepMethod.deform!(wing, zeros(wing.n_panels), zeros(wing.n_panels)) - VortexStepMethod.reinit!(body_aero) + VortexStepMethod.reinit!(body_aero; refine_mesh=false) for (idx, i) in enumerate(test_indices) reset_te = body_aero.panels[i].TE_point_1 @@ -138,6 +138,7 @@ end theta_dist = fill(deg2rad(15.0), wing.n_panels) delta_dist = fill(deg2rad(3.0), wing.n_panels) VortexStepMethod.deform!(wing, theta_dist, delta_dist) + VortexStepMethod.reinit!(body_aero; refine_mesh=false) # Store state after deformation i = length(body_aero.panels) ÷ 2 @@ -147,7 +148,8 @@ end VortexStepMethod.reinit!(body_aero; va=zeros(3), omega=zeros(3), - init_aero=true + init_aero=true, + refine_mesh=false ) end @@ -162,7 +164,7 @@ end # Test zero deformation VortexStepMethod.deform!(wing, zeros(wing.n_panels), zeros(wing.n_panels)) - VortexStepMethod.reinit!(body_aero) + VortexStepMethod.reinit!(body_aero; refine_mesh=false) @test all(p.delta ≈ 0.0 for p in body_aero.panels) # Test large deformation angles @@ -171,14 +173,14 @@ end # Should not error even with large angles VortexStepMethod.deform!(wing, theta_dist, delta_dist) - VortexStepMethod.reinit!(body_aero) + VortexStepMethod.reinit!(body_aero; refine_mesh=false) @test all(p.delta ≈ deg2rad(30.0) for p in body_aero.panels) # Test negative angles theta_dist = fill(deg2rad(-20.0), wing.n_panels) delta_dist = fill(deg2rad(-10.0), wing.n_panels) VortexStepMethod.deform!(wing, theta_dist, delta_dist) - VortexStepMethod.reinit!(body_aero) + VortexStepMethod.reinit!(body_aero; refine_mesh=false) @test all(p.delta ≈ deg2rad(-10.0) for p in body_aero.panels) end end From dad4e72a9ae82d6456b077379b9f0d9bc9a32214 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Fri, 7 Nov 2025 15:50:13 +0100 Subject: [PATCH 13/22] Add cl cd cm group array --- src/solver.jl | 38 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 36 insertions(+), 2 deletions(-) diff --git a/src/solver.jl b/src/solver.jl index 7f4daec3..742d921b 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -51,6 +51,9 @@ Struct for storing the solution of the [solve!](@ref) function. Must contain all moment_coeff_dist::MVector{P, Float64} = zeros(P) group_moment_dist::MVector{G, Float64} = zeros(G) group_moment_coeff_dist::MVector{G, Float64} = zeros(G) + cl_group_array::MVector{G, Float64} = zeros(G) + cd_group_array::MVector{G, Float64} = zeros(G) + cm_group_array::MVector{G, Float64} = zeros(G) solver_status::SolverStatus = FAILURE end @@ -296,8 +299,14 @@ function solve!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution= if length(solver.sol.group_moment_dist) > 0 group_moment_dist = solver.sol.group_moment_dist group_moment_coeff_dist = solver.sol.group_moment_coeff_dist + cl_group_array = solver.sol.cl_group_array + cd_group_array = solver.sol.cd_group_array + cm_group_array = solver.sol.cm_group_array group_moment_dist .= 0.0 group_moment_coeff_dist .= 0.0 + cl_group_array .= 0.0 + cd_group_array .= 0.0 + cm_group_array .= 0.0 panel_idx = 1 group_idx = 1 for wing in body_aero.wings @@ -306,21 +315,46 @@ function solve!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution= # Original method: divide panels into equally-sized sequential groups panels_per_group = wing.n_panels ÷ wing.n_groups for _ in 1:wing.n_groups + panel_count = 0 for _ in 1:panels_per_group group_moment_dist[group_idx] += moment_dist[panel_idx] group_moment_coeff_dist[group_idx] += moment_coeff_dist[panel_idx] + cl_group_array[group_idx] += solver.sol.cl_array[panel_idx] + cd_group_array[group_idx] += solver.sol.cd_array[panel_idx] + cm_group_array[group_idx] += solver.sol.cm_array[panel_idx] panel_idx += 1 + panel_count += 1 end + # Average the coefficients over panels in the group + cl_group_array[group_idx] /= panel_count + cd_group_array[group_idx] /= panel_count + cm_group_array[group_idx] /= panel_count group_idx += 1 end elseif wing.grouping_method == REFINE # REFINE method: group refined panels by their original unrefined section + # First pass: accumulate values + group_panel_counts = zeros(Int, wing.n_groups) for local_panel_idx in 1:wing.n_panels original_section_idx = wing.refined_panel_mapping[local_panel_idx] - group_moment_dist[group_idx + original_section_idx - 1] += moment_dist[panel_idx] - group_moment_coeff_dist[group_idx + original_section_idx - 1] += moment_coeff_dist[panel_idx] + target_group_idx = group_idx + original_section_idx - 1 + group_moment_dist[target_group_idx] += moment_dist[panel_idx] + group_moment_coeff_dist[target_group_idx] += moment_coeff_dist[panel_idx] + cl_group_array[target_group_idx] += solver.sol.cl_array[panel_idx] + cd_group_array[target_group_idx] += solver.sol.cd_array[panel_idx] + cm_group_array[target_group_idx] += solver.sol.cm_array[panel_idx] + group_panel_counts[original_section_idx] += 1 panel_idx += 1 end + # Second pass: average coefficients + for i in 1:wing.n_groups + target_group_idx = group_idx + i - 1 + if group_panel_counts[i] > 0 + cl_group_array[target_group_idx] /= group_panel_counts[i] + cd_group_array[target_group_idx] /= group_panel_counts[i] + cm_group_array[target_group_idx] /= group_panel_counts[i] + end + end group_idx += wing.n_groups end else From 1438b5572c521b247a564cfb70d8205f3b5b2c92 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Fri, 7 Nov 2025 16:23:17 +0100 Subject: [PATCH 14/22] Added tests and settings --- data/ram_air_kite/vsm_settings.yaml | 4 + src/settings.jl | 4 + test/runtests.jl | 1 + test/solver/test_group_coefficients.jl | 184 +++++++++++++++++++++++++ 4 files changed, 193 insertions(+) create mode 100644 test/solver/test_group_coefficients.jl diff --git a/data/ram_air_kite/vsm_settings.yaml b/data/ram_air_kite/vsm_settings.yaml index 5ae66dc9..8f726240 100644 --- a/data/ram_air_kite/vsm_settings.yaml +++ b/data/ram_air_kite/vsm_settings.yaml @@ -10,6 +10,9 @@ PanelDistribution: InitialGammaDistribution: ELLIPTIC: Elliptic distribution ZEROS: Constant distribution +PanelGroupingMethod: + EQUAL_SIZE: Divide panels into equally-sized sequential groups + REFINE: Group refined panels by their original unrefined section wings: - name: main_wing @@ -17,6 +20,7 @@ wings: n_groups: 40 spanwise_panel_distribution: LINEAR spanwise_direction: [0.0, 1.0, 0.0] + grouping_method: EQUAL_SIZE remove_nan: true solver_settings: n_panels: 40 diff --git a/src/settings.jl b/src/settings.jl index ad3fb75f..5916da52 100644 --- a/src/settings.jl +++ b/src/settings.jl @@ -14,6 +14,7 @@ end n_groups::Int64 = 40 spanwise_panel_distribution::PanelDistribution = LINEAR spanwise_direction::MVec3 = [0.0, 1.0, 0.0] + grouping_method::PanelGroupingMethod = EQUAL_SIZE remove_nan = true end @@ -80,6 +81,9 @@ function VSMSettings(filename; data_prefix=true) wing.n_groups = wing_data["n_groups"] wing.spanwise_panel_distribution = eval(Symbol(wing_data["spanwise_panel_distribution"])) wing.spanwise_direction = MVec3(wing_data["spanwise_direction"]) + if haskey(wing_data, "grouping_method") + wing.grouping_method = eval(Symbol(wing_data["grouping_method"])) + end wing.remove_nan = wing_data["remove_nan"] push!(vsm_settings.wings, wing) diff --git a/test/runtests.jl b/test/runtests.jl index 6c6f92cd..520eeab3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -28,6 +28,7 @@ end::Bool include("ram_geometry/test_kite_geometry.jl") include("settings/test_settings.jl") include("solver/test_solver.jl") + include("solver/test_group_coefficients.jl") include("VortexStepMethod/test_VortexStepMethod.jl") include("wake/test_wake.jl") include("wing_geometry/test_wing_geometry.jl") diff --git a/test/solver/test_group_coefficients.jl b/test/solver/test_group_coefficients.jl new file mode 100644 index 00000000..7eae60ef --- /dev/null +++ b/test/solver/test_group_coefficients.jl @@ -0,0 +1,184 @@ +using VortexStepMethod +using LinearAlgebra +using Test + +# Load test support +include("../TestSupport.jl") +using .TestSupport + +@testset "Group Coefficient Arrays Tests" begin + @testset "Group coefficients with EQUAL_SIZE method" begin + # Create a simple wing with groups + n_panels = 20 + n_groups = 4 + + # Create a test wing settings file + settings_file = create_temp_wing_settings("solver", "solver_test_wing.yaml"; alpha=5.0, beta=0.0, wind_speed=10.0) + + try + # Modify settings to use specific panel/group configuration + settings = VSMSettings(settings_file) + settings.wings[1].n_panels = n_panels + settings.wings[1].n_groups = n_groups + settings.wings[1].grouping_method = EQUAL_SIZE + settings.solver_settings.n_panels = n_panels + settings.solver_settings.n_groups = n_groups + + # Create wing and solver + wing = Wing(settings) + body_aero = BodyAerodynamics([wing]) + solver = Solver(body_aero, settings) + + # Set conditions and solve + va = [10.0, 0.0, 0.0] + set_va!(body_aero, va) + sol = solve!(solver, body_aero) + + # Test 1: Group arrays exist and have correct size + @test length(sol.cl_group_array) == n_groups + @test length(sol.cd_group_array) == n_groups + @test length(sol.cm_group_array) == n_groups + + # Test 2: Group arrays are not all zeros (solver computed them) + @test !all(sol.cl_group_array .== 0.0) + @test !all(sol.cd_group_array .== 0.0) + + # Test 3: Verify group coefficients are averages of panel coefficients + panels_per_group = n_panels ÷ n_groups + for group_idx in 1:n_groups + panel_start = (group_idx - 1) * panels_per_group + 1 + panel_end = group_idx * panels_per_group + + # Calculate expected average from panel coefficients + expected_cl = sum(sol.cl_array[panel_start:panel_end]) / panels_per_group + expected_cd = sum(sol.cd_array[panel_start:panel_end]) / panels_per_group + expected_cm = sum(sol.cm_array[panel_start:panel_end]) / panels_per_group + + # Check if group coefficients match expected averages + # Handle NaN values that can occur in INVISCID models + if isnan(expected_cl) + @test isnan(sol.cl_group_array[group_idx]) + else + @test isapprox(sol.cl_group_array[group_idx], expected_cl, rtol=1e-10) + end + if isnan(expected_cd) + @test isnan(sol.cd_group_array[group_idx]) + else + @test isapprox(sol.cd_group_array[group_idx], expected_cd, rtol=1e-10) + end + if isnan(expected_cm) + @test isnan(sol.cm_group_array[group_idx]) + else + @test isapprox(sol.cm_group_array[group_idx], expected_cm, rtol=1e-10) + end + end + + # Test 4: Verify physical consistency (lift coefficients should be positive at positive AoA) + # Skip test if values are NaN + if !any(isnan.(sol.cl_group_array)) + @test all(sol.cl_group_array .> 0.0) + end + + finally + rm(settings_file; force=true) + end + end + + @testset "Group coefficients with n_groups=0 (no grouping)" begin + # Create a wing with no groups + n_panels = 20 + n_groups = 0 + + settings_file = create_temp_wing_settings("solver", "solver_test_wing.yaml"; alpha=5.0, beta=0.0, wind_speed=10.0) + + try + settings = VSMSettings(settings_file) + settings.wings[1].n_panels = n_panels + settings.wings[1].n_groups = n_groups + settings.solver_settings.n_panels = n_panels + settings.solver_settings.n_groups = n_groups + + wing = Wing(settings) + body_aero = BodyAerodynamics([wing]) + solver = Solver(body_aero, settings) + + va = [10.0, 0.0, 0.0] + set_va!(body_aero, va) + sol = solve!(solver, body_aero) + + # Test: Group arrays should be empty when n_groups=0 + @test length(sol.cl_group_array) == 0 + @test length(sol.cd_group_array) == 0 + @test length(sol.cm_group_array) == 0 + + finally + rm(settings_file; force=true) + end + end + + @testset "Group coefficients with different group sizes" begin + # Test with various panel/group combinations + test_cases = [ + (n_panels=40, n_groups=8), + (n_panels=30, n_groups=5), + (n_panels=24, n_groups=6), + ] + + for (n_panels, n_groups) in test_cases + settings_file = create_temp_wing_settings("solver", "solver_test_wing.yaml"; alpha=5.0, beta=0.0, wind_speed=10.0) + + try + settings = VSMSettings(settings_file) + settings.wings[1].n_panels = n_panels + settings.wings[1].n_groups = n_groups + settings.wings[1].grouping_method = EQUAL_SIZE + settings.solver_settings.n_panels = n_panels + settings.solver_settings.n_groups = n_groups + + wing = Wing(settings) + body_aero = BodyAerodynamics([wing]) + solver = Solver(body_aero, settings) + + va = [10.0, 0.0, 0.0] + set_va!(body_aero, va) + sol = solve!(solver, body_aero) + + # Verify arrays have correct size + @test length(sol.cl_group_array) == n_groups + @test length(sol.cd_group_array) == n_groups + @test length(sol.cm_group_array) == n_groups + + # Verify group coefficients are computed correctly + panels_per_group = n_panels ÷ n_groups + for group_idx in 1:n_groups + panel_start = (group_idx - 1) * panels_per_group + 1 + panel_end = group_idx * panels_per_group + + expected_cl = sum(sol.cl_array[panel_start:panel_end]) / panels_per_group + expected_cd = sum(sol.cd_array[panel_start:panel_end]) / panels_per_group + expected_cm = sum(sol.cm_array[panel_start:panel_end]) / panels_per_group + + # Handle NaN for all coefficients + if isnan(expected_cl) + @test isnan(sol.cl_group_array[group_idx]) + else + @test isapprox(sol.cl_group_array[group_idx], expected_cl, rtol=1e-10) + end + if isnan(expected_cd) + @test isnan(sol.cd_group_array[group_idx]) + else + @test isapprox(sol.cd_group_array[group_idx], expected_cd, rtol=1e-10) + end + if isnan(expected_cm) + @test isnan(sol.cm_group_array[group_idx]) + else + @test isapprox(sol.cm_group_array[group_idx], expected_cm, rtol=1e-10) + end + end + + finally + rm(settings_file; force=true) + end + end + end +end From 80fc35760d0c3c95923f74b25679f5a85eb292d3 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sat, 8 Nov 2025 13:59:55 +0100 Subject: [PATCH 15/22] Improved obj file reading --- src/obj_geometry.jl | 103 ++++++++++++++---- src/wing_geometry.jl | 3 +- test/Project.toml | 1 + test/TestSupport.jl | 5 - test/ram_geometry/test_kite_geometry.jl | 58 ++++++++-- test/runtests.jl | 3 +- test/settings/test_settings.jl | 2 +- test/solver/test_group_coefficients.jl | 4 - .../test_yaml_wing_deformation.jl | 6 - 9 files changed, 136 insertions(+), 49 deletions(-) delete mode 100644 test/TestSupport.jl diff --git a/src/obj_geometry.jl b/src/obj_geometry.jl index 1dae4cbc..dcf2a2e3 100644 --- a/src/obj_geometry.jl +++ b/src/obj_geometry.jl @@ -112,7 +112,8 @@ Create interpolation functions for leading/trailing edges and area. - Where le_interp and te_interp are tuples themselves, containing the x, y and z interpolations """ function create_interpolations(vertices, circle_center_z, radius, gamma_tip, R=I(3); interp_steps=40) - gamma_range = range(-gamma_tip+1e-6, gamma_tip-1e-6, interp_steps) + gamma_range = range(-gamma_tip+gamma_tip/interp_steps*2, + gamma_tip-gamma_tip/interp_steps*2, interp_steps) stepsize = gamma_range.step.hi vz_centered = [v[3] - circle_center_z for v in vertices] @@ -122,33 +123,91 @@ function create_interpolations(vertices, circle_center_z, radius, gamma_tip, R=I leading_edges = zeros(3, length(gamma_range)) areas = zeros(length(gamma_range)) + n_slices = length(gamma_range) for (j, gamma) in enumerate(gamma_range) trailing_edges[1, j] = -Inf leading_edges[1, j] = Inf - for (i, v) in enumerate(vertices) - # Rotate y coordinate to check box containment - # rotated_y = v[2] * cos(-gamma) - vz_centered[i] * sin(-gamma) - gamma_v = atan(-v[2], vz_centered[i]) - if gamma ≤ 0 && gamma - stepsize ≤ gamma_v ≤ gamma - if v[1] > trailing_edges[1, j] - trailing_edges[:, j] .= v - te_gammas[j] = gamma_v - end - if v[1] < leading_edges[1, j] - leading_edges[:, j] .= v - le_gammas[j] = gamma_v + + # Determine if this is a tip slice and get search parameters + is_first_tip = (j == 1) + is_last_tip = (j == n_slices) + + if is_first_tip || is_last_tip + # Tip slices: use directional search within adjacent slice region + gamma_search = is_first_tip ? gamma_range[1] : gamma_range[end] + max_te_score = -Inf + max_le_score = -Inf + + for (i, v) in enumerate(vertices) + gamma_v = atan(-v[2], vz_centered[i]) + + # Check if vertex is in the adjacent slice region + in_range = if gamma_search ≤ 0 + gamma_search - stepsize ≤ gamma_v ≤ gamma_search + else + gamma_search ≤ gamma_v ≤ gamma_search + stepsize end - elseif gamma > 0 && gamma ≤ gamma_v ≤ gamma + stepsize - if v[1] > trailing_edges[1, j] - trailing_edges[:, j] .= v - te_gammas[j] = gamma_v + + if in_range + if is_first_tip + # TE: furthest in [X, Y, -Z] direction + te_score = v[1] + v[2] - v[3] + if te_score > max_te_score + trailing_edges[:, j] .= v + te_gammas[j] = gamma_v + max_te_score = te_score + end + # LE: furthest in [-X, Y, -Z] direction + le_score = -v[1] + v[2] - v[3] + if le_score > max_le_score + leading_edges[:, j] .= v + le_gammas[j] = gamma_v + max_le_score = le_score + end + else # is_last_tip + # TE: furthest in [X, -Y, -Z] direction + te_score = v[1] - v[2] - v[3] + if te_score > max_te_score + trailing_edges[:, j] .= v + te_gammas[j] = gamma_v + max_te_score = te_score + end + # LE: furthest in [-X, -Y, -Z] direction + le_score = -v[1] - v[2] - v[3] + if le_score > max_le_score + leading_edges[:, j] .= v + le_gammas[j] = gamma_v + max_le_score = le_score + end + end end - if v[1] < leading_edges[1, j] - leading_edges[:, j] .= v - le_gammas[j] = gamma_v + end + else + # Interior slices: use standard min/max x-coordinate search + for (i, v) in enumerate(vertices) + gamma_v = atan(-v[2], vz_centered[i]) + if gamma ≤ 0 && gamma - stepsize ≤ gamma_v ≤ gamma + if v[1] > trailing_edges[1, j] + trailing_edges[:, j] .= v + te_gammas[j] = gamma_v + end + if v[1] < leading_edges[1, j] + leading_edges[:, j] .= v + le_gammas[j] = gamma_v + end + elseif gamma > 0 && gamma ≤ gamma_v ≤ gamma + stepsize + if v[1] > trailing_edges[1, j] + trailing_edges[:, j] .= v + te_gammas[j] = gamma_v + end + if v[1] < leading_edges[1, j] + leading_edges[:, j] .= v + le_gammas[j] = gamma_v + end end end end + area = norm(leading_edges[:, j] - trailing_edges[:, j]) * stepsize * radius last_area = j > 1 ? areas[j-1] : 0.0 areas[j] = last_area + area @@ -159,9 +218,9 @@ function create_interpolations(vertices, circle_center_z, radius, gamma_tip, R=I trailing_edges[:, j] .= R * trailing_edges[:, j] end - le_interp = ntuple(i -> linear_interpolation(te_gammas, leading_edges[i, :], + le_interp = ntuple(i -> linear_interpolation(le_gammas, leading_edges[i, :], extrapolation_bc=Line()), 3) - te_interp = ntuple(i -> linear_interpolation(le_gammas, trailing_edges[i, :], + te_interp = ntuple(i -> linear_interpolation(te_gammas, trailing_edges[i, :], extrapolation_bc=Line()), 3) area_interp = linear_interpolation(gamma_range, areas, extrapolation_bc=Line()) diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index 79554cb4..6ac62f48 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -462,7 +462,8 @@ function deform!(wing::Wing) local_y .= normalize(section1.LE_point - section2.LE_point) chord .= section1.TE_point .- section1.LE_point normal .= chord × local_y - @. wing.refined_sections[i].TE_point = section1.LE_point + cos(wing.theta_dist[i]) * chord - sin(wing.theta_dist[i]) * normal + @. wing.refined_sections[i].TE_point = section1.LE_point + + cos(wing.theta_dist[i]) * chord - sin(wing.theta_dist[i]) * normal end return nothing end diff --git a/test/Project.toml b/test/Project.toml index 9ede8c7a..7b362c6d 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -13,6 +13,7 @@ Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +VortexStepMethod = "ed3cd733-9f0f-46a9-93e0-89b8d4998dd9" YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6" [compat] diff --git a/test/TestSupport.jl b/test/TestSupport.jl deleted file mode 100644 index 4b428502..00000000 --- a/test/TestSupport.jl +++ /dev/null @@ -1,5 +0,0 @@ -module TestSupport -export suppress_warnings, test_data_path, create_temp_wing_settings, - get_standard_wing_file, get_complete_settings_file -include("test_data_utils.jl") -end \ No newline at end of file diff --git a/test/ram_geometry/test_kite_geometry.jl b/test/ram_geometry/test_kite_geometry.jl index 47dee1d4..861bc880 100644 --- a/test/ram_geometry/test_kite_geometry.jl +++ b/test/ram_geometry/test_kite_geometry.jl @@ -189,18 +189,18 @@ using Serialization # Create an ObjWing for testing wing = ObjWing(test_obj_path, test_dat_path; remove_nan=true) body_aero = BodyAerodynamics([wing]) - + # Store original TE point for comparison i = length(body_aero.panels) ÷ 2 original_te_point = copy(body_aero.panels[i].TE_point_1) - + # Apply deformation with non-zero angles theta_dist = fill(deg2rad(30.0), wing.n_panels) # 10 degrees twist delta_dist = fill(deg2rad(5.0), wing.n_panels) # 5 degrees trailing edge deflection - + VortexStepMethod.deform!(wing, theta_dist, delta_dist) VortexStepMethod.reinit!(body_aero) - + # Check if TE point changed after deformation deformed_te_point = copy(body_aero.panels[i].TE_point_1) @test !isapprox(original_te_point, deformed_te_point, atol=1e-2) @@ -208,19 +208,61 @@ using Serialization @test deformed_te_point[2] ≈ original_te_point[2] atol=1e-2 # right hand rule @test deformed_te_point[1] < original_te_point[1] # right hand rule @test body_aero.panels[i].delta ≈ deg2rad(5.0) - + # Reset deformation with zero angles zero_theta_dist = zeros(wing.n_panels) zero_delta_dist = zeros(wing.n_panels) - + VortexStepMethod.deform!(wing, zero_theta_dist, zero_delta_dist) VortexStepMethod.reinit!(body_aero) - + # Check if TE point returned to original position reset_te_point = copy(body_aero.panels[i].TE_point_1) @test original_te_point ≈ reset_te_point atol=1e-4 end + + @testset "First and Last Section Deformation with group_deform!" begin + # Create an ObjWing with a small number of panels and groups + wing = ObjWing(test_obj_path, test_dat_path; + n_panels=4, n_groups=2, remove_nan=true) + + # Store original TE points from all refined_sections + # Wing has n_panels+1 sections (5 sections for 4 panels) + n_sections = wing.n_panels + 1 + original_te_points = [copy(wing.refined_sections[i].TE_point) + for i in 1:n_sections] + + @show wing.refined_sections[1].LE_point + @show wing.refined_sections[1].TE_point + @show wing.refined_sections[end].LE_point + @show wing.refined_sections[end].TE_point + + # Apply group_deform! with non-zero angles (2 groups, each controlling 2 panels) + theta_angles = [deg2rad(15.0), deg2rad(20.0)] + delta_angles = [deg2rad(5.0), deg2rad(10.0)] + + VortexStepMethod.group_deform!(wing, theta_angles, delta_angles; smooth=false) + + # Check that all sections' TE points have been deformed + for i in 1:n_sections + deformed_te = wing.refined_sections[i].TE_point + original_te = original_te_points[i] + + if i == 1 + # First section should be deformed + @test !isapprox(original_te, deformed_te, atol=1e-6) + @info "Section 1 (first): original=$original_te, deformed=$deformed_te" + elseif i == n_sections + # Last section (n_panels+1) should be deformed + @test !isapprox(original_te, deformed_te, atol=1e-6) + @info "Section $(n_sections) (last): original=$original_te, deformed=$deformed_te" + else + # Intermediate sections should also be deformed + @test !isapprox(original_te, deformed_te, atol=1e-6) + end + end + end rm(test_obj_path) rm(test_dat_path) -end \ No newline at end of file +end diff --git a/test/runtests.jl b/test/runtests.jl index 520eeab3..f89bea03 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,8 +2,7 @@ using Test, VortexStepMethod # Make paths robust (avoid cd("..")) cd(@__DIR__) # ensure we're in test/ no matter how tests are launched -include("TestSupport.jl") -using .TestSupport +include("test_data_utils.jl") println("Running tests...") diff --git a/test/settings/test_settings.jl b/test/settings/test_settings.jl index 37b3d920..e3200766 100644 --- a/test/settings/test_settings.jl +++ b/test/settings/test_settings.jl @@ -11,6 +11,6 @@ using Test @test vss.wings isa Vector{WingSettings} @test length(vss.wings) == 2 io = IOBuffer(repr(vss)) - @test countlines(io) == 44 # Updated to match new output format + @test countlines(io) == 46 # Updated to match new output format end nothing diff --git a/test/solver/test_group_coefficients.jl b/test/solver/test_group_coefficients.jl index 7eae60ef..ede46b2a 100644 --- a/test/solver/test_group_coefficients.jl +++ b/test/solver/test_group_coefficients.jl @@ -2,10 +2,6 @@ using VortexStepMethod using LinearAlgebra using Test -# Load test support -include("../TestSupport.jl") -using .TestSupport - @testset "Group Coefficient Arrays Tests" begin @testset "Group coefficients with EQUAL_SIZE method" begin # Create a simple wing with groups diff --git a/test/yaml_geometry/test_yaml_wing_deformation.jl b/test/yaml_geometry/test_yaml_wing_deformation.jl index 580e3dc9..578528b2 100644 --- a/test/yaml_geometry/test_yaml_wing_deformation.jl +++ b/test/yaml_geometry/test_yaml_wing_deformation.jl @@ -2,12 +2,6 @@ using VortexStepMethod using LinearAlgebra using Test -# Load TestSupport if not already loaded (for standalone execution) -if !@isdefined(TestSupport) - include(joinpath(@__DIR__, "..", "TestSupport.jl")) - using .TestSupport -end - @testset "YAML Wing Deformation Tests" begin @testset "Simple Wing Deformation" begin # Load existing simple_wing.yaml From b3ef2b0a1b597d5e5d65136ec15c9896a53e4b6c Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sat, 8 Nov 2025 14:11:12 +0100 Subject: [PATCH 16/22] Fixed deform --- src/wing_geometry.jl | 37 ++++++++++++++---- test/ram_geometry/test_kite_geometry.jl | 5 --- test/runtests.jl | 52 +++++++++++++++++-------- 3 files changed, 65 insertions(+), 29 deletions(-) diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index 6ac62f48..c212e8ab 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -456,14 +456,37 @@ function deform!(wing::Wing) chord = zeros(MVec3) normal = zeros(MVec3) - for i in 1:wing.n_panels - section1 = wing.non_deformed_sections[i] - section2 = wing.non_deformed_sections[i+1] - local_y .= normalize(section1.LE_point - section2.LE_point) - chord .= section1.TE_point .- section1.LE_point + # Process all sections (n_panels + 1) + # Each section gets angle(s) from adjacent panel(s) + for i in 1:(wing.n_panels + 1) + section = wing.non_deformed_sections[i] + + # Determine the angle for this section + if i == 1 + # First section: use angle from first panel + theta = wing.theta_dist[1] + elseif i == wing.n_panels + 1 + # Last section: use angle from last panel + theta = wing.theta_dist[wing.n_panels] + else + # Middle sections: average of adjacent panels + theta = 0.5 * (wing.theta_dist[i-1] + wing.theta_dist[i]) + end + + # Compute local coordinate system + if i < wing.n_panels + 1 + section2 = wing.non_deformed_sections[i+1] + local_y .= normalize(section.LE_point - section2.LE_point) + else + # For last section, use same local_y as previous + section_prev = wing.non_deformed_sections[i-1] + local_y .= normalize(section_prev.LE_point - section.LE_point) + end + + chord .= section.TE_point .- section.LE_point normal .= chord × local_y - @. wing.refined_sections[i].TE_point = section1.LE_point + - cos(wing.theta_dist[i]) * chord - sin(wing.theta_dist[i]) * normal + @. wing.refined_sections[i].TE_point = section.LE_point + + cos(theta) * chord - sin(theta) * normal end return nothing end diff --git a/test/ram_geometry/test_kite_geometry.jl b/test/ram_geometry/test_kite_geometry.jl index 861bc880..4eb07823 100644 --- a/test/ram_geometry/test_kite_geometry.jl +++ b/test/ram_geometry/test_kite_geometry.jl @@ -232,11 +232,6 @@ using Serialization original_te_points = [copy(wing.refined_sections[i].TE_point) for i in 1:n_sections] - @show wing.refined_sections[1].LE_point - @show wing.refined_sections[1].TE_point - @show wing.refined_sections[end].LE_point - @show wing.refined_sections[end].TE_point - # Apply group_deform! with non-zero angles (2 groups, each controlling 2 panels) theta_angles = [deg2rad(15.0), deg2rad(20.0)] delta_angles = [deg2rad(5.0), deg2rad(10.0)] diff --git a/test/runtests.jl b/test/runtests.jl index f89bea03..7a2abdd5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,7 +4,25 @@ using Test, VortexStepMethod cd(@__DIR__) # ensure we're in test/ no matter how tests are launched include("test_data_utils.jl") +# Support selective test execution via ]test test_args=["pattern"] +const test_patterns = isempty(ARGS) ? String[] : ARGS + println("Running tests...") +if !isempty(test_patterns) + println("Filtering tests matching: ", test_patterns) +end + +# Helper to check if a test file matches any pattern +function should_run_test(test_path::String) + isempty(test_patterns) && return true + for pattern in test_patterns + # Match directory (e.g., "solver") or specific file (e.g., "test_group_coefficients") + if occursin(pattern, test_path) + return true + end + end + return false +end # keep your env check as-is... const build_is_production_build_env_name = "BUILD_IS_PRODUCTION_BUILD" @@ -14,25 +32,25 @@ const build_is_production_build = let v = get(ENV, build_is_production_build_env end::Bool @testset verbose = true "Testing VortexStepMethod..." begin - if build_is_production_build + if build_is_production_build && should_run_test("bench") include("bench.jl") end - include("body_aerodynamics/test_body_aerodynamics.jl") - include("body_aerodynamics/test_results.jl") - include("filament/test_bound_filament.jl") - include("filament/test_semi_infinite_filament.jl") - include("panel/test_panel.jl") - include("plotting/test_plotting.jl") - include("polars/test_polars.jl") - include("ram_geometry/test_kite_geometry.jl") - include("settings/test_settings.jl") - include("solver/test_solver.jl") - include("solver/test_group_coefficients.jl") - include("VortexStepMethod/test_VortexStepMethod.jl") - include("wake/test_wake.jl") - include("wing_geometry/test_wing_geometry.jl") - include("yaml_geometry/test_yaml_geometry.jl") - include("Aqua.jl") + should_run_test("body_aerodynamics/test_body_aerodynamics.jl") && include("body_aerodynamics/test_body_aerodynamics.jl") + should_run_test("body_aerodynamics/test_results.jl") && include("body_aerodynamics/test_results.jl") + should_run_test("filament/test_bound_filament.jl") && include("filament/test_bound_filament.jl") + should_run_test("filament/test_semi_infinite_filament.jl") && include("filament/test_semi_infinite_filament.jl") + should_run_test("panel/test_panel.jl") && include("panel/test_panel.jl") + should_run_test("plotting/test_plotting.jl") && include("plotting/test_plotting.jl") + should_run_test("polars/test_polars.jl") && include("polars/test_polars.jl") + should_run_test("ram_geometry/test_kite_geometry.jl") && include("ram_geometry/test_kite_geometry.jl") + should_run_test("settings/test_settings.jl") && include("settings/test_settings.jl") + should_run_test("solver/test_solver.jl") && include("solver/test_solver.jl") + should_run_test("solver/test_group_coefficients.jl") && include("solver/test_group_coefficients.jl") + should_run_test("VortexStepMethod/test_VortexStepMethod.jl") && include("VortexStepMethod/test_VortexStepMethod.jl") + should_run_test("wake/test_wake.jl") && include("wake/test_wake.jl") + should_run_test("wing_geometry/test_wing_geometry.jl") && include("wing_geometry/test_wing_geometry.jl") + should_run_test("yaml_geometry/test_yaml_geometry.jl") && include("yaml_geometry/test_yaml_geometry.jl") + should_run_test("Aqua.jl") && include("Aqua.jl") end nothing From 8e28721446031a55ec1b6cfd2e523d165d37c065 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sat, 8 Nov 2025 14:19:28 +0100 Subject: [PATCH 17/22] Test on 1.10 and 1.11 --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 6407e6c9..48fdc83d 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -23,7 +23,7 @@ jobs: matrix: version: - '1.10' - - '1' + - '1.11' os: - ubuntu-latest build_is_production_build: From a6a8337f85ab814ed86a85dcd0bd084cd217e842 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sat, 8 Nov 2025 17:00:37 +0100 Subject: [PATCH 18/22] Disable mapping kwarg --- src/body_aerodynamics.jl | 9 ++++++--- src/wing_geometry.jl | 22 ++++++++++++++-------- 2 files changed, 20 insertions(+), 11 deletions(-) diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 103745dc..a05a289b 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -114,7 +114,7 @@ function Base.setproperty!(obj::BodyAerodynamics, sym::Symbol, val) end """ - reinit!(body_aero::BodyAerodynamics; init_aero, va, omega, refine_mesh) + reinit!(body_aero::BodyAerodynamics; init_aero, va, omega, refine_mesh, recompute_mapping) Initialize a BodyAerodynamics struct in-place by setting up panels and coefficients. @@ -127,6 +127,8 @@ Initialize a BodyAerodynamics struct in-place by setting up panels and coefficie - `omega=zeros(3)`: Turn rate in kite body frame x y and z - `refine_mesh=true`: Whether to refine wing meshes. Set to `false` after `deform!()` to preserve deformed geometry. +- `recompute_mapping=true`: Whether to recompute the refined panel mapping. + Set to `false` to skip mapping computation when it hasn't changed. # Returns nothing @@ -135,12 +137,13 @@ function reinit!(body_aero::BodyAerodynamics; init_aero=true, va=[15.0, 0.0, 0.0], omega=zeros(MVec3), - refine_mesh=true + refine_mesh=true, + recompute_mapping=true ) idx = 1 vec = @MVector zeros(3) for wing in body_aero.wings - reinit!(wing; refine_mesh=refine_mesh) + reinit!(wing; refine_mesh, recompute_mapping) panel_props = wing.panel_props # Create panels diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index c212e8ab..a5bcdbce 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -302,18 +302,20 @@ function Wing(n_panels::Int; end """ - reinit!(wing::AbstractWing; refine_mesh=true) + reinit!(wing::AbstractWing; refine_mesh=true, recompute_mapping=true) Reinitialize wing geometry and panel properties. # Keyword Arguments - `refine_mesh::Bool=true`: Whether to refine the mesh. Set to `false` after `deform!()` to preserve deformed geometry while updating panel properties. +- `recompute_mapping::Bool=true`: Whether to recompute the refined panel mapping. + Set to `false` to skip mapping computation when it hasn't changed. """ -function reinit!(wing::AbstractWing; refine_mesh=true) +function reinit!(wing::AbstractWing; refine_mesh=true, recompute_mapping=true) # Refine mesh unless explicitly disabled (e.g., to preserve deformation) if refine_mesh - refine_aerodynamic_mesh!(wing) + refine_aerodynamic_mesh!(wing; recompute_mapping) end # Calculate panel properties @@ -587,14 +589,18 @@ function update_non_deformed_sections!(wing::AbstractWing) end """ - refine_aerodynamic_mesh!(wing::AbstractWing) + refine_aerodynamic_mesh!(wing::AbstractWing; recompute_mapping=true) Refine the aerodynamic mesh of the wing based on spanwise panel distribution. +# Keyword Arguments +- `recompute_mapping::Bool=true`: Whether to recompute the refined panel mapping. + Set to `false` to skip mapping computation when it hasn't changed. + Returns: Vector{Section}: List of refined sections """ -function refine_aerodynamic_mesh!(wing::AbstractWing) +function refine_aerodynamic_mesh!(wing::AbstractWing; recompute_mapping=true) sort!(wing.sections, by=s -> s.LE_point[2], rev=true) n_sections = wing.n_panels + 1 if length(wing.refined_sections) == 0 @@ -631,7 +637,7 @@ function refine_aerodynamic_mesh!(wing::AbstractWing) for i in eachindex(wing.sections) reinit!(wing.refined_sections[i], wing.sections[i]) end - compute_refined_panel_mapping!(wing) + recompute_mapping && compute_refined_panel_mapping!(wing) update_non_deformed_sections!(wing) return nothing end @@ -642,7 +648,7 @@ function refine_aerodynamic_mesh!(wing::AbstractWing) if n_sections == 2 reinit!(wing.refined_sections[1], LE[1,:], TE[1,:], aero_model[1], aero_data[1]) reinit!(wing.refined_sections[2], LE[end,:], TE[end,:], aero_model[end], aero_data[end]) - compute_refined_panel_mapping!(wing) + recompute_mapping && compute_refined_panel_mapping!(wing) update_non_deformed_sections!(wing) return nothing end @@ -666,7 +672,7 @@ function refine_aerodynamic_mesh!(wing::AbstractWing) end # Compute panel mapping by finding closest unrefined panel for each refined panel - compute_refined_panel_mapping!(wing) + recompute_mapping && compute_refined_panel_mapping!(wing) # Validate REFINE grouping method if wing.grouping_method == REFINE && wing.n_groups > 0 From 3502599eb0f7c2f1d1a54473ae519d5c5de78ed5 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sat, 8 Nov 2025 17:33:57 +0100 Subject: [PATCH 19/22] Add kwarg to sort sections --- src/body_aerodynamics.jl | 9 ++++++--- src/wing_geometry.jl | 17 +++++++++++------ 2 files changed, 17 insertions(+), 9 deletions(-) diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index a05a289b..5d5ca33f 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -114,7 +114,7 @@ function Base.setproperty!(obj::BodyAerodynamics, sym::Symbol, val) end """ - reinit!(body_aero::BodyAerodynamics; init_aero, va, omega, refine_mesh, recompute_mapping) + reinit!(body_aero::BodyAerodynamics; init_aero, va, omega, refine_mesh, recompute_mapping, sort_sections) Initialize a BodyAerodynamics struct in-place by setting up panels and coefficients. @@ -129,6 +129,8 @@ Initialize a BodyAerodynamics struct in-place by setting up panels and coefficie `deform!()` to preserve deformed geometry. - `recompute_mapping=true`: Whether to recompute the refined panel mapping. Set to `false` to skip mapping computation when it hasn't changed. +- `sort_sections=true`: Whether to sort sections by spanwise position. + Set to `false` for REFINE wings where section order is determined by structural connectivity. # Returns nothing @@ -138,12 +140,13 @@ function reinit!(body_aero::BodyAerodynamics; va=[15.0, 0.0, 0.0], omega=zeros(MVec3), refine_mesh=true, - recompute_mapping=true + recompute_mapping=true, + sort_sections=true ) idx = 1 vec = @MVector zeros(3) for wing in body_aero.wings - reinit!(wing; refine_mesh, recompute_mapping) + reinit!(wing; refine_mesh, recompute_mapping, sort_sections) panel_props = wing.panel_props # Create panels diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index a5bcdbce..83889eff 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -302,7 +302,7 @@ function Wing(n_panels::Int; end """ - reinit!(wing::AbstractWing; refine_mesh=true, recompute_mapping=true) + reinit!(wing::AbstractWing; refine_mesh=true, recompute_mapping=true, sort_sections=true) Reinitialize wing geometry and panel properties. @@ -311,11 +311,13 @@ Reinitialize wing geometry and panel properties. `deform!()` to preserve deformed geometry while updating panel properties. - `recompute_mapping::Bool=true`: Whether to recompute the refined panel mapping. Set to `false` to skip mapping computation when it hasn't changed. +- `sort_sections::Bool=true`: Whether to sort sections by spanwise position. + Set to `false` for REFINE wings where section order is determined by structural connectivity. """ -function reinit!(wing::AbstractWing; refine_mesh=true, recompute_mapping=true) +function reinit!(wing::AbstractWing; refine_mesh=true, recompute_mapping=true, sort_sections=true) # Refine mesh unless explicitly disabled (e.g., to preserve deformation) if refine_mesh - refine_aerodynamic_mesh!(wing; recompute_mapping) + refine_aerodynamic_mesh!(wing; recompute_mapping, sort_sections) end # Calculate panel properties @@ -589,19 +591,22 @@ function update_non_deformed_sections!(wing::AbstractWing) end """ - refine_aerodynamic_mesh!(wing::AbstractWing; recompute_mapping=true) + refine_aerodynamic_mesh!(wing::AbstractWing; recompute_mapping=true, sort_sections=true) Refine the aerodynamic mesh of the wing based on spanwise panel distribution. # Keyword Arguments - `recompute_mapping::Bool=true`: Whether to recompute the refined panel mapping. Set to `false` to skip mapping computation when it hasn't changed. +- `sort_sections::Bool=true`: Whether to sort sections by spanwise position (y-coordinate). + Set to `false` for REFINE wings where section order is determined by structural connectivity. Returns: Vector{Section}: List of refined sections """ -function refine_aerodynamic_mesh!(wing::AbstractWing; recompute_mapping=true) - sort!(wing.sections, by=s -> s.LE_point[2], rev=true) +function refine_aerodynamic_mesh!(wing::AbstractWing; recompute_mapping=true, sort_sections=true) + # Only sort sections if requested (skip for REFINE wings with fixed structural order) + sort_sections && sort!(wing.sections, by=s -> s.LE_point[2], rev=true) n_sections = wing.n_panels + 1 if length(wing.refined_sections) == 0 if wing.spanwise_distribution == UNCHANGED || length(wing.sections) == n_sections From f55597f729d9755cd44cdd017ef5de2e7f397467 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sat, 8 Nov 2025 18:36:12 +0100 Subject: [PATCH 20/22] Add groups --- src/body_aerodynamics.jl | 12 +++++++-- src/solver.jl | 53 ++++++++++++++++++++++++++++++++++++++-- 2 files changed, 61 insertions(+), 4 deletions(-) diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 5d5ca33f..b96ff57d 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -24,6 +24,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru @with_kw mutable struct BodyAerodynamics{P} panels::Vector{Panel} wings::Vector{Wing} + groups::Vector{Panel} = Panel[] _va::MVec3 = zeros(MVec3) omega::MVec3 = zeros(MVec3) gamma_distribution::MVector{P, Float64} = zeros(P) @@ -73,6 +74,7 @@ function BodyAerodynamics( ) where T <: AbstractWing # Initialize panels panels = Panel[] + n_groups = 0 for wing in wings for section in wing.sections section.LE_point .-= kite_body_origin @@ -80,7 +82,7 @@ function BodyAerodynamics( end if wing.spanwise_distribution == UNCHANGED wing.refined_sections = wing.sections - !(wing.n_panels == length(wing.sections) - 1) && + !(wing.n_panels == length(wing.sections) - 1) && throw(ArgumentError("(wing.n_panels = $(wing.n_panels)) != (length(wing.sections) - 1 = $(length(wing.sections) - 1))")) else wing.refined_sections = Section[Section() for _ in 1:wing.n_panels+1] @@ -91,9 +93,15 @@ function BodyAerodynamics( panel = Panel() push!(panels, panel) end + + # Count total groups + n_groups += wing.n_groups end - body_aero = BodyAerodynamics{length(panels)}(; panels, wings) + # Initialize groups (unrefined panel representatives) + groups = [Panel() for _ in 1:n_groups] + + body_aero = BodyAerodynamics{length(panels)}(; panels, wings, groups) reinit!(body_aero; va, omega) return body_aero end diff --git a/src/solver.jl b/src/solver.jl index 742d921b..ed6f524b 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -316,43 +316,92 @@ function solve!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution= panels_per_group = wing.n_panels ÷ wing.n_groups for _ in 1:wing.n_groups panel_count = 0 + group_panel = body_aero.groups[group_idx] + # Zero out accumulated fields + group_panel.x_airf .= 0.0 + group_panel.y_airf .= 0.0 + group_panel.z_airf .= 0.0 + group_panel.va .= 0.0 + group_panel.chord = 0.0 + group_panel.width = 0.0 for _ in 1:panels_per_group + panel = body_aero.panels[panel_idx] group_moment_dist[group_idx] += moment_dist[panel_idx] group_moment_coeff_dist[group_idx] += moment_coeff_dist[panel_idx] cl_group_array[group_idx] += solver.sol.cl_array[panel_idx] cd_group_array[group_idx] += solver.sol.cd_array[panel_idx] cm_group_array[group_idx] += solver.sol.cm_array[panel_idx] + # Accumulate geometry for averaging + group_panel.x_airf .+= panel.x_airf + group_panel.y_airf .+= panel.y_airf + group_panel.z_airf .+= panel.z_airf + group_panel.va .+= panel.va + group_panel.chord += panel.chord + group_panel.width += panel.width panel_idx += 1 panel_count += 1 end - # Average the coefficients over panels in the group + # Average the coefficients and geometry over panels in the group cl_group_array[group_idx] /= panel_count cd_group_array[group_idx] /= panel_count cm_group_array[group_idx] /= panel_count + group_panel.x_airf ./= panel_count + group_panel.y_airf ./= panel_count + group_panel.z_airf ./= panel_count + group_panel.va ./= panel_count + group_panel.chord /= panel_count + group_panel.width /= panel_count group_idx += 1 end elseif wing.grouping_method == REFINE # REFINE method: group refined panels by their original unrefined section + # Initialize group panels + for i in 1:wing.n_groups + target_group_idx = group_idx + i - 1 + group_panel = body_aero.groups[target_group_idx] + group_panel.x_airf .= 0.0 + group_panel.y_airf .= 0.0 + group_panel.z_airf .= 0.0 + group_panel.va .= 0.0 + group_panel.chord = 0.0 + group_panel.width = 0.0 + end # First pass: accumulate values group_panel_counts = zeros(Int, wing.n_groups) for local_panel_idx in 1:wing.n_panels + panel = body_aero.panels[panel_idx] original_section_idx = wing.refined_panel_mapping[local_panel_idx] target_group_idx = group_idx + original_section_idx - 1 + group_panel = body_aero.groups[target_group_idx] group_moment_dist[target_group_idx] += moment_dist[panel_idx] group_moment_coeff_dist[target_group_idx] += moment_coeff_dist[panel_idx] cl_group_array[target_group_idx] += solver.sol.cl_array[panel_idx] cd_group_array[target_group_idx] += solver.sol.cd_array[panel_idx] cm_group_array[target_group_idx] += solver.sol.cm_array[panel_idx] + # Accumulate geometry + group_panel.x_airf .+= panel.x_airf + group_panel.y_airf .+= panel.y_airf + group_panel.z_airf .+= panel.z_airf + group_panel.va .+= panel.va + group_panel.chord += panel.chord + group_panel.width += panel.width group_panel_counts[original_section_idx] += 1 panel_idx += 1 end - # Second pass: average coefficients + # Second pass: average coefficients and geometry for i in 1:wing.n_groups target_group_idx = group_idx + i - 1 if group_panel_counts[i] > 0 + group_panel = body_aero.groups[target_group_idx] cl_group_array[target_group_idx] /= group_panel_counts[i] cd_group_array[target_group_idx] /= group_panel_counts[i] cm_group_array[target_group_idx] /= group_panel_counts[i] + group_panel.x_airf ./= group_panel_counts[i] + group_panel.y_airf ./= group_panel_counts[i] + group_panel.z_airf ./= group_panel_counts[i] + group_panel.va ./= group_panel_counts[i] + group_panel.chord /= group_panel_counts[i] + group_panel.width /= group_panel_counts[i] end end group_idx += wing.n_groups From 872bdfd7e4df92ff75c4493edd8128b38bbb6a6c Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sun, 9 Nov 2025 14:00:50 +0100 Subject: [PATCH 21/22] Correct width --- src/solver.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/solver.jl b/src/solver.jl index ed6f524b..cb7cc6be 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -54,6 +54,7 @@ Struct for storing the solution of the [solve!](@ref) function. Must contain all cl_group_array::MVector{G, Float64} = zeros(G) cd_group_array::MVector{G, Float64} = zeros(G) cm_group_array::MVector{G, Float64} = zeros(G) + alpha_group_array::MVector{G, Float64} = zeros(G) solver_status::SolverStatus = FAILURE end @@ -302,11 +303,13 @@ function solve!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution= cl_group_array = solver.sol.cl_group_array cd_group_array = solver.sol.cd_group_array cm_group_array = solver.sol.cm_group_array + alpha_group_array = solver.sol.alpha_group_array group_moment_dist .= 0.0 group_moment_coeff_dist .= 0.0 cl_group_array .= 0.0 cd_group_array .= 0.0 cm_group_array .= 0.0 + alpha_group_array .= 0.0 panel_idx = 1 group_idx = 1 for wing in body_aero.wings @@ -331,6 +334,7 @@ function solve!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution= cl_group_array[group_idx] += solver.sol.cl_array[panel_idx] cd_group_array[group_idx] += solver.sol.cd_array[panel_idx] cm_group_array[group_idx] += solver.sol.cm_array[panel_idx] + alpha_group_array[group_idx] += solver.sol.alpha_array[panel_idx] # Accumulate geometry for averaging group_panel.x_airf .+= panel.x_airf group_panel.y_airf .+= panel.y_airf @@ -345,12 +349,12 @@ function solve!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution= cl_group_array[group_idx] /= panel_count cd_group_array[group_idx] /= panel_count cm_group_array[group_idx] /= panel_count + alpha_group_array[group_idx] /= panel_count group_panel.x_airf ./= panel_count group_panel.y_airf ./= panel_count group_panel.z_airf ./= panel_count group_panel.va ./= panel_count group_panel.chord /= panel_count - group_panel.width /= panel_count group_idx += 1 end elseif wing.grouping_method == REFINE @@ -378,6 +382,7 @@ function solve!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution= cl_group_array[target_group_idx] += solver.sol.cl_array[panel_idx] cd_group_array[target_group_idx] += solver.sol.cd_array[panel_idx] cm_group_array[target_group_idx] += solver.sol.cm_array[panel_idx] + alpha_group_array[target_group_idx] += solver.sol.alpha_array[panel_idx] # Accumulate geometry group_panel.x_airf .+= panel.x_airf group_panel.y_airf .+= panel.y_airf @@ -396,6 +401,7 @@ function solve!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution= cl_group_array[target_group_idx] /= group_panel_counts[i] cd_group_array[target_group_idx] /= group_panel_counts[i] cm_group_array[target_group_idx] /= group_panel_counts[i] + alpha_group_array[target_group_idx] /= group_panel_counts[i] group_panel.x_airf ./= group_panel_counts[i] group_panel.y_airf ./= group_panel_counts[i] group_panel.z_airf ./= group_panel_counts[i] From a1df71207653e22e5b716a6824659d2ffefc9dd4 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Mon, 10 Nov 2025 16:15:32 +0100 Subject: [PATCH 22/22] Add makie plotting --- examples/Project.toml | 1 + examples/ram_air_kite.jl | 4 +- ext/VortexStepMethodControlPlotsExt.jl | 920 ++++++++++++++++++++++++- ext/VortexStepMethodMakieExt.jl | 690 ++++++++++++++++++- src/panel.jl | 2 +- src/plotting.jl | 920 ------------------------- 6 files changed, 1612 insertions(+), 925 deletions(-) delete mode 100644 src/plotting.jl diff --git a/examples/Project.toml b/examples/Project.toml index b94307b9..fea9b2fd 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -2,6 +2,7 @@ CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" ControlPlots = "23c2ee80-7a9e-4350-b264-8e670f12517c" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" VortexStepMethod = "ed3cd733-9f0f-46a9-93e0-89b8d4998dd9" diff --git a/examples/ram_air_kite.jl b/examples/ram_air_kite.jl index 13e60c62..bbf1480b 100644 --- a/examples/ram_air_kite.jl +++ b/examples/ram_air_kite.jl @@ -1,4 +1,4 @@ -using ControlPlots +using GLMakie using VortexStepMethod using LinearAlgebra @@ -118,4 +118,4 @@ PLOT && plot_polars( is_show=true, use_tex=USE_TEX ) -nothing \ No newline at end of file +nothing diff --git a/ext/VortexStepMethodControlPlotsExt.jl b/ext/VortexStepMethodControlPlotsExt.jl index ce9494be..b2d09e9f 100644 --- a/ext/VortexStepMethodControlPlotsExt.jl +++ b/ext/VortexStepMethodControlPlotsExt.jl @@ -5,6 +5,924 @@ import VortexStepMethod: calculate_filaments_for_plotting export plot_wing, plot_circulation_distribution, plot_geometry, plot_distribution, plot_polars, save_plot, show_plot, plot_polar_data -include("../src/plotting.jl") +""" + set_plot_style(titel_size=16; use_tex=false) + +Set the default style for plots using LaTeX. +`` +# Arguments: +- `titel_size`: size of the plot title in points (default: 16) +- `ùse_tex`: if the external `pdflatex` command shall be used +""" +function set_plot_style(titel_size=16; use_tex=false) + rcParams = plt.PyDict(plt.matplotlib."rcParams") + rcParams["text.usetex"] = use_tex + rcParams["font.family"] = "serif" + if use_tex + rcParams["font.serif"] = ["Computer Modern Roman"] + end + rcParams["axes.titlesize"] = titel_size + rcParams["axes.labelsize"] = 12 + rcParams["axes.linewidth"] = 1 + rcParams["lines.linewidth"] = 1 + rcParams["lines.markersize"] = 6 + rcParams["xtick.labelsize"] = 10 + rcParams["ytick.labelsize"] = 10 + rcParams["legend.fontsize"] = 10 + rcParams["figure.titlesize"] = 16 + if use_tex + rcParams["pgf.texsystem"] = "pdflatex" # Use pdflatex + end + rcParams["pgf.rcfonts"] = false + rcParams["figure.figsize"] = (10, 6) # Default figure size +end + + +""" + save_plot(fig, save_path, title; data_type=".pdf") + +Save a plot to a file. + +# Arguments +- `fig`: Plots figure object +- `save_path`: Path to save the plot +- `title`: Title of the plot + +# Keyword arguments +- `data_type`: File extension (default: ".pdf") +""" +function VortexStepMethod.save_plot(fig, save_path, title; data_type=".pdf") + isnothing(save_path) && throw(ArgumentError("save_path should be provided")) + + !isdir(save_path) && mkpath(save_path) + full_path = joinpath(save_path, title * data_type) + + @debug "Attempting to save figure to: $full_path" + @debug "Current working directory: $(pwd())" + + try + fig.savefig(full_path) + @debug "Figure saved as $data_type" + + if isfile(full_path) + @debug "File successfully saved to $full_path" + @debug "File size: $(filesize(full_path)) bytes" + else + @info "File does not exist after save attempt: $full_path" + end + catch e + @error "Error saving figure: $e" + @error "Error type: $(typeof(e))" + rethrow(e) + end +end + +""" + show_plot(fig; dpi=130) + +Display a plot at specified DPI. + +# Arguments +- `fig`: Plots figure object + +# Keyword arguments +- `dpi`: Dots per inch for the figure (default: 130) +""" +function VortexStepMethod.show_plot(fig; dpi=130) + plt.display(fig) +end + +""" + plot_line_segment!(ax, segment, color, label; width=3) + +Plot a line segment in 3D with arrow. + +# Arguments +- `ax`: Plot axis +- `segment`: Array of two points defining the segment +- `color`: Color of the segment +- `label`: Label for the legend + +# Keyword Arguments +- `width`: Line width (default: 3) +""" +function plot_line_segment!(ax, segment, color, label; width=3) + ax.plot( + [segment[1][1], segment[2][1]], + [segment[1][2], segment[2][2]], + [segment[1][3], segment[2][3]], + color=color, label=label, linewidth=width + ) + + dir = segment[2] - segment[1] + ax.quiver( + [segment[1][1]], [segment[1][2]], [segment[1][3]], + [dir[1]], [dir[2]], [dir[3]], + color=color + ) +end + +""" + set_axes_equal!(ax; zoom=1.8) + +Set 3D plot axes to equal scale. + +# Arguments +- ax: 3D plot axis + +# Keyword arguments +zoom: zoom factor (default: 1.8) +""" +function set_axes_equal!(ax; zoom=1.8) + x_lims = ax.get_xlim3d() ./ zoom + y_lims = ax.get_ylim3d() ./ zoom + z_lims = ax.get_zlim3d() ./ zoom + + x_range = abs(x_lims[2] - x_lims[1]) + y_range = abs(y_lims[2] - y_lims[1]) + z_range = abs(z_lims[2] - z_lims[1]) + + max_range = max(x_range, y_range, z_range) + + x_mid = mean(x_lims) + y_mid = mean(y_lims) + z_mid = mean(z_lims) + + ax.set_xlim3d((x_mid - max_range / 2, x_mid + max_range / 2)) + ax.set_ylim3d((y_mid - max_range / 2, y_mid + max_range / 2)) + ax.set_zlim3d((z_mid - max_range / 2, z_mid + max_range / 2)) +end + +""" + create_geometry_plot(body_aero::BodyAerodynamics, title, view_elevation, view_azimuth; + zoom=1.8, use_tex=false) + +Create a 3D plot of wing geometry including panels and filaments. + +# Arguments +- body_aero: struct of type BodyAerodynamics +- title: plot title +- view_elevation: initial view elevation angle [°] +- view_azimuth: initial view azimuth angle [°] + +# Keyword arguments +- zoom: zoom factor (default: 1.8) +""" +function create_geometry_plot(body_aero::BodyAerodynamics, title, view_elevation, view_azimuth; + zoom=1.8, use_tex=false) + set_plot_style(28; use_tex) + + panels = body_aero.panels + va = isa(body_aero.va, Tuple) ? body_aero.va[1] : body_aero.va + + # Extract geometric data + corner_points = [panel.corner_points for panel in panels] + control_points = [panel.control_point for panel in panels] + aero_centers = [panel.aero_center for panel in panels] + + # Create plot + fig = plt.figure(figsize=(14, 14)) + ax = fig.add_subplot(111, projection="3d") + ax.set_title(title) + + # Plot panels + legend_used = Dict{String,Bool}() + for (i, panel) in enumerate(panels) + # Plot panel edges and surfaces + corners = Matrix{Float64}(panel.corner_points) + x_corners = corners[1, :] + y_corners = corners[2, :] + z_corners = corners[3, :] + + push!(x_corners, x_corners[1]) + push!(y_corners, y_corners[1]) + push!(z_corners, z_corners[1]) + + ax.plot(x_corners, + y_corners, + z_corners, + color=:grey, + linewidth=1, + label=i == 1 ? "Panel Edges" : "") + + # Plot control points and aerodynamic centers + ax.scatter([control_points[i][1]], [control_points[i][2]], [control_points[i][3]], + color=:green, label=i == 1 ? "Control Points" : "") + ax.scatter([aero_centers[i][1]], [aero_centers[i][2]], [aero_centers[i][3]], + color=:blue, label=i == 1 ? "Aerodynamic Centers" : "") + + # Plot filaments + filaments = calculate_filaments_for_plotting(panel) + legends = ["Bound Vortex", "side1", "side2", "wake_1", "wake_2"] + + for (filament, legend) in zip(filaments, legends) + x1, x2, color = filament + @debug "Legend: $legend" + show_legend = !get(legend_used, legend, false) + plot_line_segment!(ax, [x1, x2], color, show_legend ? legend : "") + legend_used[legend] = true + end + end + + # Plot velocity vector + max_chord = maximum(panel.chord for panel in panels) + va_mag = norm(va) + va_vector_begin = -2 * max_chord * va / va_mag + va_vector_end = va_vector_begin + 1.5 * va / va_mag + plot_line_segment!(ax, [va_vector_begin, va_vector_end], :lightblue, "va") + + # Add legends for the first occurrence of each label + handles, labels = ax.get_legend_handles_labels() + # by_label = Dict(zip(labels, handles)) + # ax.legend(values(by_label), keys(by_label), bbox_to_anchor=(0, 0, 1.1, 1)) + + # Set labels and make axes equal + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.set_zlabel("z") + set_axes_equal!(ax; zoom) + + # Set the initial view + ax.view_init(elev=view_elevation, azim=view_azimuth) + + # Ensure the figure is fully rendered + # fig.canvas.draw() + plt.tight_layout(rect=(0, 0, 1, 0.97)) + + return fig +end + +""" + plot_geometry(body_aero::BodyAerodynamics, title; + data_type=".pdf", save_path=nothing, + is_save=false, is_show=false, + view_elevation=15, view_azimuth=-120, use_tex=false) + +Plot wing geometry from different viewpoints and optionally save/show plots. + +# Arguments: +- `body_aero`: the [BodyAerodynamics](@ref) to plot +- title: plot title + +# Keyword arguments: +- `data_type``: string with the file type postfix (default: ".pdf") +- `save_path`: path for saving the graphic (default: `nothing`) +- `is_save`: boolean value, indicates if the graphic shall be saved (default: `false`) +- `is_show`: boolean value, indicates if the graphic shall be displayed (default: `false`) +- `view_elevation`: initial view elevation angle (default: 15) [°] +- `view_azimuth`: initial view azimuth angle (default: -120) [°] +- `use_tex`: if the external `pdflatex` command shall be used (default: false) + +""" +function VortexStepMethod.plot_geometry(body_aero::BodyAerodynamics, title; + data_type=".pdf", + save_path=nothing, + is_save=false, + is_show=false, + view_elevation=15, + view_azimuth=-120, + use_tex=false) + + if is_save + plt.ioff() + # Angled view + fig = create_geometry_plot(body_aero, "$(title)_angled_view", 15, -120; use_tex) + save_plot(fig, save_path, "$(title)_angled_view", data_type=data_type) + + # Top view + fig = create_geometry_plot(body_aero, "$(title)_top_view", 90, 0; use_tex) + save_plot(fig, save_path, "$(title)_top_view", data_type=data_type) + + # Front view + fig = create_geometry_plot(body_aero, "$(title)_front_view", 0, 0; use_tex) + save_plot(fig, save_path, "$(title)_front_view", data_type=data_type) + + # Side view + fig = create_geometry_plot(body_aero, "$(title)_side_view", 0, -90; use_tex) + save_plot(fig, save_path, "$(title)_side_view", data_type=data_type) + end + + if is_show + plt.ion() + fig = create_geometry_plot(body_aero, title, view_elevation, view_azimuth; use_tex) + plt.display(fig) + else + fig = create_geometry_plot(body_aero, title, view_elevation, view_azimuth; use_tex) + end + fig +end + +""" + plot_distribution(y_coordinates_list, results_list, label_list; + title="spanwise_distribution", data_type=".pdf", + save_path=nothing, is_save=false, is_show=true, use_tex=false) + +Plot spanwise distributions of aerodynamic properties. + +# Arguments +- `y_coordinates_list`: List of spanwise coordinates +- `results_list`: List of result dictionaries +- `label_list`: List of labels for different results + +# Keyword arguments +- `title`: Plot title (default: "spanwise_distribution") +- `data_type`: File extension for saving (default: ".pdf") +- `save_path`: Path to save plots (default: nothing) +- `is_save`: Whether to save plots (default: false) +- `is_show`: Whether to display plots (default: true) +- `use_tex`: if the external `pdflatex` command shall be used +""" +function VortexStepMethod.plot_distribution(y_coordinates_list, results_list, label_list; + title="spanwise_distribution", + data_type=".pdf", + save_path=nothing, + is_save=false, + is_show=true, + use_tex=false) + + length(results_list) == length(label_list) || throw(ArgumentError( + "Number of results ($(length(results_list))) must match number of labels ($(length(label_list)))" + )) + + # Set the plot style + set_plot_style(; use_tex) + + # Initializing plot + fig, axs = plt.subplots(3, 3, figsize=(16, 10)) + fig.suptitle(title, fontsize=16) + + # CL plot + for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list) + value = "$(round(result_i["cl"], digits=2))" + if label_i == "LLT" + label = label_i * L" $~C_\mathrm{L}$: " * value + else + label = label_i * L" $C_\mathrm{L}$: " * value + end + axs[1, 1].plot( + y_coordinates_i, + result_i["cl_distribution"], + label=label + ) + end + axs[1, 1].set_title(L"$C_\mathrm{L}$ Distribution", size=16) + axs[1, 1].set_xlabel(L"Spanwise Position $y/b$") + axs[1, 1].set_ylabel(L"Lift Coefficient $C_\mathrm{L}$") + axs[1, 1].legend() + + # CD plot + for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list) + value = "$(round(result_i["cl"], digits=2))" + if label_i == "LLT" + label = label_i * L" $~C_\mathrm{D}$: " * value + else + label = label_i * L" $C_\mathrm{D}$: " * value + end + axs[1, 2].plot( + y_coordinates_i, + result_i["cd_distribution"], + label=label + ) + end + axs[1, 2].set_title(L"$C_\mathrm{D}$ Distribution", size=16) + axs[1, 2].set_xlabel(L"Spanwise Position $y/b$") + axs[1, 2].set_ylabel(L"Drag Coefficient $C_\mathrm{D}$") + axs[1, 2].legend() + + # Gamma Distribution + for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list) + axs[1, 3].plot( + y_coordinates_i, + result_i["gamma_distribution"], + label=label_i + ) + end + axs[1, 3].set_title(L"\Gamma~Distribution", size=16) + axs[1, 3].set_xlabel(L"Spanwise Position $y/b$") + axs[1, 3].set_ylabel(L"Circulation~\Gamma") + axs[1, 3].legend() + + # Geometric Alpha + for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list) + axs[2, 1].plot( + y_coordinates_i, + result_i["alpha_geometric"], + label=label_i + ) + end + axs[2, 1].set_title(L"$\alpha$ Geometric", size=16) + axs[2, 1].set_xlabel(L"Spanwise Position $y/b$") + axs[2, 1].set_ylabel(L"Angle of Attack $\alpha$ (deg)") + axs[2, 1].legend() + + # Calculated/ Corrected Alpha + for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list) + axs[2, 2].plot( + y_coordinates_i, + result_i["alpha_at_ac"], + label=label_i + ) + end + axs[2, 2].set_title(L"$\alpha$ result (corrected to aerodynamic center)", size=16) + axs[2, 2].set_xlabel(L"Spanwise Position $y/b$") + axs[2, 2].set_ylabel(L"Angle of Attack $\alpha$ (deg)") + axs[2, 2].legend() + + # Uncorrected Alpha plot + for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list) + axs[2, 3].plot( + y_coordinates_i, + result_i["alpha_uncorrected"], + label=label_i + ) + end + axs[2, 3].set_title(L"$\alpha$ Uncorrected (if VSM, at the control point)", size=16) + axs[2, 3].set_xlabel(L"Spanwise Position $y/b$") + axs[2, 3].set_ylabel(L"Angle of Attack $\alpha$ (deg)") + axs[2, 3].legend() + + # Force Components + for (idx, component) in enumerate(["x", "y", "z"]) + axs[3, idx].set_title("Force in $component direction", size=16) + axs[3, idx].set_xlabel(L"Spanwise Position $y/b$") + axs[3, idx].set_ylabel(raw"$F_\mathrm" * "{$component}" * raw"$") + for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) + # Extract force components for the current direction (idx) + forces = results["F_distribution"][idx, :] + # Verify dimensions match + if length(y_coords) != length(forces) + @warn "Dimension mismatch in force plotting" length(y_coords) length(forces) component + continue # Skip this component instead of throwing error + end + space = "" + if label == "LLT" + space = "~" + end + axs[3, idx].plot( + y_coords, + forces, + label="$label" * space * raw"$~\Sigma~F_\mathrm" * "{$component}:~" * + raw"$" * "$(round(results["F$component"], digits=2)) N" + ) + axs[3, idx].legend() + end + end + + fig.tight_layout() + + # Save and show plot + if is_save + save_plot(fig, save_path, title, data_type=data_type) + end + + if is_show + show_plot(fig) + end + + return fig +end + +""" + generate_polar_data(solver, body_aero::BodyAerodynamics, angle_range; + angle_type="angle_of_attack", angle_of_attack=0.0, + side_slip=0.0, v_a=10.0, use_latex=false) + +Generate polar data for aerodynamic analysis over a range of angles. + +# Arguments +- `solver`: Aerodynamic solver object +- `body_aero`: Wing aerodynamics struct +- `angle_range`: Range of angles to analyze + +# Keyword arguments +- `angle_type`: Type of angle variation ("angle_of_attack" or "side_slip") +- `angle_of_attack`: Initial angle of attack [rad] +- `side_slip`: Initial side slip angle in [rad] +- `v_a`: norm of apparent wind speed [m/s] + +# Returns +- Tuple of polar data array and Reynolds number +""" +function generate_polar_data( + solver, + body_aero::BodyAerodynamics, + angle_range; + angle_type="angle_of_attack", + angle_of_attack=0.0, + side_slip=0.0, + v_a=10.0, + use_latex=false +) + n_panels = length(body_aero.panels) + n_angles = length(angle_range) + + # Initialize arrays + cl = zeros(n_angles) + cd = zeros(n_angles) + cs = zeros(n_angles) + gamma_distribution = zeros(n_angles, n_panels) + cl_distribution = zeros(n_angles, n_panels) + cd_distribution = zeros(n_angles, n_panels) + cs_distribution = zeros(n_angles, n_panels) + reynolds_number = zeros(n_angles) + + # Previous gamma for initialization + gamma = nothing + + for (i, angle_i) in enumerate(angle_range) + # Set angle based on type + if angle_type == "angle_of_attack" + α = deg2rad(angle_i) + β = side_slip + elseif angle_type == raw"side_slip" + α = angle_of_attack + β = deg2rad(angle_i) + else + throw(ArgumentError("angle_type must be 'angle_of_attack' or 'side_slip'")) + end + + # Update inflow conditions + set_va!( + body_aero, + [ + cos(α) * cos(β), + sin(β), + sin(α) + ] * v_a + ) + + # Solve and store results + results = solve(solver, body_aero, gamma_distribution[i, :]) + + cl[i] = results["cl"] + cd[i] = results["cd"] + cs[i] = results["cs"] + gamma_distribution[i, :] = results["gamma_distribution"] + cl_distribution[i, :] = results["cl_distribution"] + cd_distribution[i, :] = results["cd_distribution"] + cs_distribution[i, :] = results["cs_distribution"] + reynolds_number[i] = results["Rey"] + + # Store gamma for next iteration + gamma = gamma_distribution[i, :] + end + + polar_data = [ + angle_range, + cl, + cd, + cs, + gamma_distribution, + cl_distribution, + cd_distribution, + cs_distribution, + reynolds_number + ] + + return polar_data, reynolds_number[1] +end + +""" + plot_polars(solver_list, body_aero_list, label_list; + literature_path_list=String[], + angle_range=range(0, 20, 2), angle_type="angle_of_attack", + angle_of_attack=0.0, side_slip=0.0, v_a=10.0, + title="polar", data_type=".pdf", save_path=nothing, + is_save=true, is_show=true, use_tex=false) + +Plot polar data comparing different solvers and configurations. + +# Arguments +- `solver_list`: List of aerodynamic solvers +- `body_aero_list`: List of wing aerodynamics objects +- `label_list`: List of labels for each configuration + +# Keyword arguments +- `literature_path_list`: Optional paths to literature data files +- `angle_range`: Range of angles to analyze [°] +- `angle_type`: "`angle_of_attack`" or "`side_slip`"; (default: `angle_of_attack`) +- `angle_of_attack:` AoA to be used for plotting the polars (default: 0.0) [rad] +- `side_slip`: side slip angle (default: 0.0) [rad] +- v_a: norm of apparent wind speed (default: 10.0) [m/s] +- title: plot title +- `data_type`: File extension for saving (default: ".pdf") +- `save_path`: Path to save plots (default: nothing) +- `is_save`: Whether to save plots (default: true) +- `is_show`: Whether to display plots (default: true) +- `use_tex`: if the external `pdflatex` command shall be used (default: false) +""" +function VortexStepMethod.plot_polars( + solver_list, + body_aero_list, + label_list; + literature_path_list=String[], + angle_range=range(0, 20, 2), + angle_type="angle_of_attack", + angle_of_attack=0.0, + side_slip=0.0, + v_a=10.0, + title="polar", + data_type=".pdf", + save_path=nothing, + is_save=true, + is_show=true, + use_tex=false +) + # Validate inputs + total_cases = length(body_aero_list) + length(literature_path_list) + if total_cases != length(label_list) || length(solver_list) != length(body_aero_list) + throw(ArgumentError("Mismatch in number of solvers ($(length(solver_list))), " * + "cases ($total_cases), and labels ($(length(label_list)))")) + end + main_title = replace(title, " " => "_") + set_plot_style(; use_tex) + + # Generate polar data + polar_data_list = [] + for (i, (solver, body_aero)) in enumerate(zip(solver_list, body_aero_list)) + polar_data, rey = generate_polar_data( + solver, body_aero, angle_range; + angle_type, + angle_of_attack, + side_slip, + v_a + ) + push!(polar_data_list, polar_data) + # Update label with Reynolds number + label_list[i] = "$(label_list[i]) Re = $(round(Int64, rey*1e-5))e5" + end + # Load literature data if provided + if !isempty(literature_path_list) + for path in literature_path_list + data = readdlm(path, ',') + header = lowercase.(string.(data[1, :])) + # Find column indices for alpha, CL, CD, CS (case-insensitive, allow common variants) + alpha_idx = findfirst(x -> occursin("alpha", x), header) + cl_idx = findfirst(x -> occursin("cl", x), header) + cd_idx = findfirst(x -> occursin("cd", x), header) + cs_idx = findfirst(x -> occursin("cs", x), header) + # Fallback: if CS not found, fill with zeros + cs_col = cs_idx === nothing ? zeros(size(data, 1)-1) : data[2:end, cs_idx] + # Push as [alpha, CL, CD, CS] + push!(polar_data_list, [ + data[2:end, alpha_idx], + data[2:end, cl_idx], + data[2:end, cd_idx], + cs_col + ]) + end + end + + # Initializing plot + fig, axs = plt.subplots(2, 2, figsize=(14, 14)) + + # Number of computational results (excluding literature) + n_solvers = length(solver_list) + for (i, (polar_data, label)) in enumerate(zip(polar_data_list, label_list)) + if i < n_solvers + linestyle = "-" + marker = "*" + markersize = 7 + else + linestyle = "-" + marker = "." + markersize = 5 + end + if contains(label, "LLT") + label = replace(label, "e5" => raw"\cdot10^5") + label = replace(label, " " => raw"~") + label = replace(label, "LLT" => raw"\mathrm{LLT}{~\,}") + label = raw"$" * label * raw"$" + else + label = replace(label, "e5" => raw"\cdot10^5") + label = replace(label, " " => "~") + label = replace(label, "VSM" => raw"\mathrm{VSM}") + label = raw"$" * label * raw"$" + end + axs[1, 1].plot( + polar_data[1], + polar_data[2], + label=label, + linestyle=linestyle, + marker=marker, + markersize=markersize, + ) + # Limit y-range if CL > 10 + if maximum(polar_data[2]) > 10 + axs[1, 1].set_ylim([-0.5, 2]) + end + title = raw"$C_\mathrm{L}" * raw"$" * " vs $angle_type [°]" + axs[1, 1].set_title(title) + axs[1, 1].set_xlabel("$angle_type [°]") + axs[1, 1].set_ylabel(L"$C_\mathrm{L}$") + axs[1, 1].legend() + end + + for (i, (polar_data, label)) in enumerate(zip(polar_data_list, label_list)) + if i < n_solvers + linestyle = "-" + marker = "*" + markersize = 7 + else + linestyle = "-" + marker = "." + markersize = 5 + end + if contains(label, "LLT") + label = replace(label, "e5" => raw"\cdot10^5") + label = replace(label, " " => raw"~") + label = replace(label, "LLT" => raw"\mathrm{LLT}{~\,}") + label = raw"$" * label * raw"$" + else + label = replace(label, "e5" => raw"\cdot10^5") + label = replace(label, " " => "~") + label = replace(label, "VSM" => raw"\mathrm{VSM}") + label = raw"$" * label * raw"$" + end + axs[1, 2].plot( + polar_data[1], + polar_data[3], + label=label, + linestyle=linestyle, + marker=marker, + markersize=markersize, + ) + # Limit y-range if CL > 10 + if maximum(polar_data[2]) > 10 + axs[1, 2].set_ylim([-0.5, 2]) + end + title = raw"$C_\mathrm{D}" * raw"$" * " vs $angle_type [°]" + axs[1, 2].set_title(title) + axs[1, 2].set_xlabel("$angle_type [°]") + axs[1, 2].set_ylabel(L"$C_\mathrm{D}$") + axs[1, 2].legend() + end + + + for (i, (polar_data, label)) in enumerate(zip(polar_data_list, label_list)) + if i < n_solvers + linestyle = "-" + marker = "*" + markersize = 7 + else + linestyle = "-" + marker = "." + markersize = 5 + end + if contains(label, "LLT") + label = replace(label, "e5" => raw"\cdot10^5") + label = replace(label, " " => raw"~") + label = replace(label, "LLT" => raw"\mathrm{LLT}{~\,}") + label = raw"$" * label * raw"$" + else + label = replace(label, "e5" => raw"\cdot10^5") + label = replace(label, " " => "~") + label = replace(label, "VSM" => raw"\mathrm{VSM}") + label = raw"$" * label * raw"$" + end + axs[2, 1].plot( + polar_data[1], + polar_data[4], + label=label, + linestyle=linestyle, + marker=marker, + markersize=markersize, + ) + # Limit y-range if CL > 10 + if maximum(polar_data[2]) > 10 + axs[2, 1].set_ylim([-0.5, 2]) + end + title = raw"$C_\mathrm{S}" * raw"$" * " vs $angle_type [°]" + axs[2, 1].set_title(title) + axs[2, 1].set_xlabel("$angle_type [°]") + axs[2, 1].set_ylabel(L"$C_\mathrm{S}$") + axs[2, 1].legend() + end + + for (i, (polar_data, label)) in enumerate(zip(polar_data_list, label_list)) + if i < n_solvers + linestyle = "-" + marker = "*" + markersize = 7 + else + linestyle = "-" + marker = "." + markersize = 5 + end + if contains(label, "LLT") + label = replace(label, "e5" => raw"\cdot10^5") + label = replace(label, " " => raw"~") + label = replace(label, "LLT" => raw"\mathrm{LLT}{~\,}") + label = raw"$" * label * raw"$" + else + label = replace(label, "e5" => raw"\cdot10^5") + label = replace(label, " " => "~") + label = replace(label, "VSM" => raw"\mathrm{VSM}") + label = raw"$" * label * raw"$" + end + axs[2, 2].plot( + polar_data[3], + polar_data[2], + label=label, + linestyle=linestyle, + marker=marker, + markersize=markersize, + ) + # Limit y-range if CL > 10 + if maximum(polar_data[2]) > 10 || maximum(polar_data[3]) > 10 + axs[2, 2].set_ylim([-0.5, 2]) + axs[2, 2].set_xlim([-0.5, 2]) + end + title = raw"$C_\mathrm{L}" * raw"$" * " vs " * raw"$C_\mathrm{D}" * raw"$" + axs[2, 2].set_title(title) + axs[2, 2].set_xlabel(L"$C_\mathrm{D}$") + axs[2, 2].set_ylabel(L"$C_\mathrm{L}$") + axs[2, 2].legend() + end + + fig.tight_layout(h_pad=3.5, rect=(0.01, 0.01, 0.99, 0.99)) + + # Save and show plot + if is_save && !isnothing(save_path) + save_plot(fig, save_path, main_title; data_type) + end + + if is_show + show_plot(fig) + end + + return fig +end + +""" + plot_polar_data(body_aero::BodyAerodynamics; alphas=collect(deg2rad.(-5:0.3:25)), delta_tes=collect(deg2rad.(-5:0.3:25))) + +Plot polar data (Cl, Cd, Cm) as 3D surfaces against alpha and delta_te angles. delta_te is the trailing edge deflection angle +relative to the 2d airfoil or panel chord line. + +# Arguments +- `body_aero`: Wing aerodynamics struct + +# Keyword arguments +- `alphas`: Range of angle of attack values in radians (default: -5° to 25° in 0.3° steps) +- `delta_tes`: Range of trailing edge angles in radians (default: -5° to 25° in 0.3° steps) +- `is_show`: Whether to display plots (default: true) +- `use_tex`: if the external `pdflatex` command shall be used +""" +function VortexStepMethod.plot_polar_data(body_aero::BodyAerodynamics; + alphas=collect(deg2rad.(-5:0.3:25)), + delta_tes = collect(deg2rad.(-5:0.3:25)), + is_show = true, + use_tex = false + ) + if body_aero.panels[1].aero_model == POLAR_MATRICES + set_plot_style() + + # Create figure with subplots + fig = plt.figure(figsize=(15, 6)) + + # Get interpolation functions and labels + interp_data = [ + (body_aero.panels[1].cl_interp, L"$C_l$"), + (body_aero.panels[1].cd_interp, L"$C_d$"), + (body_aero.panels[1].cm_interp, L"$C_m$") + ] + + # Create each subplot + for (idx, (interp, label)) in enumerate(interp_data) + ax = fig.add_subplot(1, 3, idx, projection="3d") + + # Create interpolation matrix + interp_matrix = zeros(length(alphas), length(delta_tes)) + interp_matrix .= [interp(alpha, delta_te) for alpha in alphas, delta_te in delta_tes] + X = collect(delta_tes) .+ zeros(length(alphas))' + Y = collect(alphas)' .+ zeros(length(delta_tes)) + + # Plot surface + ax.plot_wireframe(X, Y, interp_matrix, + edgecolor="blue", + lw=0.5, + rstride=5, + cstride=5, + alpha=0.6) + + # Set labels and title + ax.set_xlabel(L"$\delta$ [rad]") + ax.set_ylabel(L"$\alpha$ [rad]") + ax.set_zlabel(label) + ax.set_title(label * L" vs $\alpha$ and $\delta$") + ax.grid(true) + end + + # Adjust layout and display + plt.tight_layout(rect=(0.01, 0.01, 0.99, 0.99)) + if is_show + show_plot(fig) + end + return fig + else + throw(ArgumentError("Plotting polar data for $(body_aero.panels[1].aero_model) is not implemented.")) + end +end end \ No newline at end of file diff --git a/ext/VortexStepMethodMakieExt.jl b/ext/VortexStepMethodMakieExt.jl index f72ea3d8..69819d2a 100644 --- a/ext/VortexStepMethodMakieExt.jl +++ b/ext/VortexStepMethodMakieExt.jl @@ -1,5 +1,8 @@ module VortexStepMethodMakieExt -using Makie, VortexStepMethod +using Makie, VortexStepMethod, LinearAlgebra, Statistics, DelimitedFiles +import VortexStepMethod: calculate_filaments_for_plotting + +export plot_geometry, plot_distribution, plot_polars, save_plot, show_plot, plot_polar_data """ plot!(ax, panel::VortexStepMethod.Panel; kwargs...) @@ -62,4 +65,689 @@ function Makie.plot(body_aero::VortexStepMethod.BodyAerodynamics; size = (1200, return fig end +""" + save_plot(fig, save_path, title; data_type=".pdf") + +Save a Makie figure to a file. + +# Arguments +- `fig`: Makie Figure object +- `save_path`: Path to save the plot +- `title`: Title of the plot + +# Keyword arguments +- `data_type`: File extension (default: ".pdf") +""" +function VortexStepMethod.save_plot(fig, save_path, title; data_type=".pdf") + isnothing(save_path) && throw(ArgumentError("save_path should be provided")) + + !isdir(save_path) && mkpath(save_path) + full_path = joinpath(save_path, title * data_type) + + @debug "Attempting to save figure to: $full_path" + @debug "Current working directory: $(pwd())" + + try + save(full_path, fig) + @debug "Figure saved as $data_type" + + if isfile(full_path) + @debug "File successfully saved to $full_path" + @debug "File size: $(filesize(full_path)) bytes" + else + @info "File does not exist after save attempt: $full_path" + end + catch e + @error "Error saving figure: $e" + @error "Error type: $(typeof(e))" + rethrow(e) + end +end + +""" + show_plot(fig; dpi=130) + +Display a Makie figure. + +# Arguments +- `fig`: Makie Figure object + +# Keyword arguments +- `dpi`: Dots per inch for the figure (default: 130) - currently unused in Makie +""" +function VortexStepMethod.show_plot(fig; dpi=130) + display(fig) +end + +""" + plot_line_segment_makie!(ax, segment, color, label; width=3) + +Plot a line segment in 3D with arrow using Makie. + +# Arguments +- `ax`: Makie Axis3 +- `segment`: Array of two points defining the segment +- `color`: Color of the segment +- `label`: Label for the legend + +# Keyword Arguments +- `width`: Line width (default: 3) +""" +function plot_line_segment_makie!(ax, segment, color, label; width=3) + # Plot line + lines!(ax, [Point3f(segment[1]), Point3f(segment[2])]; + color=color, linewidth=width, label=label) + + # Plot arrow + dir = segment[2] - segment[1] + arrows!(ax, [Point3f(segment[1])], [Point3f(dir)]; + color=color, arrowsize=0.1) +end + +""" + set_axes_equal_makie!(ax, panels; zoom=1.8) + +Set 3D Makie axis to equal scale based on panel data. + +# Arguments +- `ax`: Makie Axis3 +- `panels`: Array of panels +- `zoom`: zoom factor (default: 1.8) +""" +function set_axes_equal_makie!(ax, panels; zoom=1.8) + # Calculate bounds from all panels + all_x = Float64[] + all_y = Float64[] + all_z = Float64[] + + for panel in panels + for i in 1:4 + push!(all_x, panel.corner_points[1, i]) + push!(all_y, panel.corner_points[2, i]) + push!(all_z, panel.corner_points[3, i]) + end + end + + x_range = (maximum(all_x) - minimum(all_x)) / zoom + y_range = (maximum(all_y) - minimum(all_y)) / zoom + z_range = (maximum(all_z) - minimum(all_z)) / zoom + + max_range = max(x_range, y_range, z_range) + + x_mid = mean([maximum(all_x), minimum(all_x)]) + y_mid = mean([maximum(all_y), minimum(all_y)]) + z_mid = mean([maximum(all_z), minimum(all_z)]) + + limits!(ax, + x_mid - max_range/2, x_mid + max_range/2, + y_mid - max_range/2, y_mid + max_range/2, + z_mid - max_range/2, z_mid + max_range/2) +end + +""" + create_geometry_plot_makie(body_aero::BodyAerodynamics, title, + view_elevation, view_azimuth; zoom=1.8) + +Create a 3D Makie plot of wing geometry including panels and filaments. + +# Arguments +- `body_aero`: struct of type BodyAerodynamics +- `title`: plot title +- `view_elevation`: initial view elevation angle [°] +- `view_azimuth`: initial view azimuth angle [°] + +# Keyword arguments +- `zoom`: zoom factor (default: 1.8) +""" +function create_geometry_plot_makie(body_aero::BodyAerodynamics, title, + view_elevation, view_azimuth; zoom=1.8) + panels = body_aero.panels + va = isa(body_aero.va, Tuple) ? body_aero.va[1] : body_aero.va + + # Create figure + fig = Figure(size=(1400, 1400)) + ax = Axis3(fig[1, 1]; + title=title, + xlabel="x", ylabel="y", zlabel="z", + aspect=:data, + azimuth=deg2rad(view_azimuth), + elevation=deg2rad(view_elevation)) + + # Plot panels + legend_used = Dict{String,Bool}() + for (i, panel) in enumerate(panels) + # Panel edges + corners = [Point3f(panel.corner_points[:, j]) for j in 1:4] + push!(corners, corners[1]) + lines!(ax, corners; color=:grey, linewidth=1, + label = i == 1 ? "Panel Edges" : nothing) + + # Control points + scatter!(ax, [Point3f(panel.control_point)]; + color=:green, markersize=10, + label = i == 1 ? "Control Points" : nothing) + + # Aerodynamic centers + scatter!(ax, [Point3f(panel.aero_center)]; + color=:blue, markersize=10, + label = i == 1 ? "Aerodynamic Centers" : nothing) + + # Plot filaments + filaments = calculate_filaments_for_plotting(panel) + legends = ["Bound Vortex", "side1", "side2", "wake_1", "wake_2"] + + for (filament, legend) in zip(filaments, legends) + x1, x2, color = filament + show_legend = !get(legend_used, legend, false) + plot_line_segment_makie!(ax, [x1, x2], color, + show_legend ? legend : nothing) + legend_used[legend] = true + end + end + + # Plot velocity vector + max_chord = maximum(panel.chord for panel in panels) + va_mag = norm(va) + va_vector_begin = -2 * max_chord * va / va_mag + va_vector_end = va_vector_begin + 1.5 * va / va_mag + plot_line_segment_makie!(ax, [va_vector_begin, va_vector_end], :lightblue, "va") + + # Set equal axes + set_axes_equal_makie!(ax, panels; zoom) + + # Add legend + axislegend(ax; position=:lt) + + return fig +end + +""" + plot_geometry(body_aero::BodyAerodynamics, title; + data_type=".pdf", save_path=nothing, + is_save=false, is_show=false, + view_elevation=15, view_azimuth=-120, use_tex=false) + +Plot wing geometry from different viewpoints using Makie. + +# Arguments: +- `body_aero`: the BodyAerodynamics to plot +- `title`: plot title + +# Keyword arguments: +- `data_type`: File extension (default: ".pdf") +- `save_path`: Path for saving (default: nothing) +- `is_save`: Whether to save (default: false) +- `is_show`: Whether to display (default: false) +- `view_elevation`: View elevation angle [°] (default: 15) +- `view_azimuth`: View azimuth angle [°] (default: -120) +- `use_tex`: Ignored for Makie (default: false) +""" +function VortexStepMethod.plot_geometry(body_aero::BodyAerodynamics, title; + data_type=".pdf", + save_path=nothing, + is_save=false, + is_show=false, + view_elevation=15, + view_azimuth=-120, + use_tex=false) + + if is_save + # Angled view + fig = create_geometry_plot_makie(body_aero, "$(title)_angled_view", 15, -120) + save_plot(fig, save_path, "$(title)_angled_view", data_type=data_type) + + # Top view + fig = create_geometry_plot_makie(body_aero, "$(title)_top_view", 90, 0) + save_plot(fig, save_path, "$(title)_top_view", data_type=data_type) + + # Front view + fig = create_geometry_plot_makie(body_aero, "$(title)_front_view", 0, 0) + save_plot(fig, save_path, "$(title)_front_view", data_type=data_type) + + # Side view + fig = create_geometry_plot_makie(body_aero, "$(title)_side_view", 0, -90) + save_plot(fig, save_path, "$(title)_side_view", data_type=data_type) + end + + fig = create_geometry_plot_makie(body_aero, title, view_elevation, view_azimuth) + + if is_show + display(fig) + end + + return fig +end + +""" + plot_distribution(y_coordinates_list, results_list, label_list; + title="spanwise_distribution", data_type=".pdf", + save_path=nothing, is_save=false, is_show=true, use_tex=false) + +Plot spanwise distributions of aerodynamic properties using Makie. + +# Arguments +- `y_coordinates_list`: List of spanwise coordinates +- `results_list`: List of result dictionaries +- `label_list`: List of labels for different results + +# Keyword arguments +- `title`: Plot title (default: "spanwise_distribution") +- `data_type`: File extension (default: ".pdf") +- `save_path`: Path to save plots (default: nothing) +- `is_save`: Whether to save (default: false) +- `is_show`: Whether to display (default: true) +- `use_tex`: Ignored for Makie (default: false) +""" +function VortexStepMethod.plot_distribution(y_coordinates_list, results_list, label_list; + title="spanwise_distribution", + data_type=".pdf", + save_path=nothing, + is_save=false, + is_show=true, + use_tex=false) + + length(results_list) == length(label_list) || throw(ArgumentError( + "Number of results ($(length(results_list))) must match labels ($(length(label_list)))" + )) + + # Create figure with 3x3 grid + fig = Figure(size=(1600, 1000)) + 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 Γ") + + # 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)") + + # 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") + + # Plot CL + for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) + value = round(results["cl"], digits=2) + lines!(ax_cl, Vector(y_coords), Vector(results["cl_distribution"]), + label="$label CL: $value") + end + axislegend(ax_cl, position=:lt) + + # Plot CD + for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) + value = round(results["cd"], digits=2) + lines!(ax_cd, Vector(y_coords), Vector(results["cd_distribution"]), + label="$label CD: $value") + end + axislegend(ax_cd, position=:lt) + + # Plot Gamma + for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) + lines!(ax_gamma, Vector(y_coords), Vector(results["gamma_distribution"]), + label=label) + end + axislegend(ax_gamma, position=:lt) + + # Plot alpha geometric + for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) + lines!(ax_alpha_geo, Vector(y_coords), Vector(results["alpha_geometric"]), + label=label) + end + axislegend(ax_alpha_geo, position=:lt) + + # Plot alpha at ac + for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) + lines!(ax_alpha_ac, Vector(y_coords), Vector(results["alpha_at_ac"]), + label=label) + end + axislegend(ax_alpha_ac, position=:lt) + + # Plot alpha uncorrected + for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) + lines!(ax_alpha_unc, Vector(y_coords), Vector(results["alpha_uncorrected"]), + label=label) + end + axislegend(ax_alpha_unc, position=:lt) + + # Plot force components + force_axes = [ax_fx, ax_fy, ax_fz] + components = ["x", "y", "z"] + for (idx, (ax, comp)) in enumerate(zip(force_axes, components)) + for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) + forces = results["F_distribution"][idx, :] + if length(y_coords) != length(forces) + @warn "Dimension mismatch" length(y_coords) length(forces) comp + continue + end + total_force = round(results["F$comp"], digits=2) + lines!(ax, Vector(y_coords), Vector(forces), + label="$label ΣF$comp: $total_force N") + end + axislegend(ax, position=:lt) + end + + # Save and show + if is_save + save_plot(fig, save_path, title, data_type=data_type) + end + + if is_show + display(fig) + end + + return fig +end + +""" + generate_polar_data(solver, body_aero::BodyAerodynamics, angle_range; + angle_type="angle_of_attack", angle_of_attack=0.0, + side_slip=0.0, v_a=10.0) + +Generate polar data for aerodynamic analysis over a range of angles. + +# Arguments +- `solver`: Aerodynamic solver object +- `body_aero`: Wing aerodynamics struct +- `angle_range`: Range of angles to analyze + +# Keyword arguments +- `angle_type`: Type of angle variation ("angle_of_attack" or "side_slip") +- `angle_of_attack`: Initial angle of attack [rad] +- `side_slip`: Initial side slip angle [rad] +- `v_a`: norm of apparent wind speed [m/s] + +# Returns +- Tuple of polar data array and Reynolds number +""" +function generate_polar_data_makie( + solver, + body_aero::BodyAerodynamics, + angle_range; + angle_type="angle_of_attack", + angle_of_attack=0.0, + side_slip=0.0, + v_a=10.0 +) + n_panels = length(body_aero.panels) + n_angles = length(angle_range) + + # Initialize arrays + cl = zeros(n_angles) + cd = zeros(n_angles) + cs = zeros(n_angles) + gamma_distribution = zeros(n_angles, n_panels) + reynolds_number = zeros(n_angles) + + for (i, angle_i) in enumerate(angle_range) + # Set angle based on type + if angle_type == "angle_of_attack" + α = deg2rad(angle_i) + β = side_slip + elseif angle_type == "side_slip" + α = angle_of_attack + β = deg2rad(angle_i) + else + throw(ArgumentError("angle_type must be 'angle_of_attack' or 'side_slip'")) + end + + # Update inflow conditions + set_va!( + body_aero, + [ + cos(α) * cos(β), + sin(β), + sin(α) + ] * v_a + ) + + # Solve and store results + results = solve(solver, body_aero, gamma_distribution[i, :]) + + cl[i] = results["cl"] + cd[i] = results["cd"] + cs[i] = results["cs"] + gamma_distribution[i, :] = results["gamma_distribution"] + reynolds_number[i] = results["Rey"] + end + + polar_data = [angle_range, cl, cd, cs] + return polar_data, reynolds_number[1] +end + +""" + plot_polars(solver_list, body_aero_list, label_list; + literature_path_list=String[], + angle_range=range(0, 20, 2), angle_type="angle_of_attack", + angle_of_attack=0.0, side_slip=0.0, v_a=10.0, + title="polar", data_type=".pdf", save_path=nothing, + is_save=true, is_show=true, use_tex=false) + +Plot polar data comparing different solvers using Makie. + +# Arguments +- `solver_list`: List of aerodynamic solvers +- `body_aero_list`: List of wing aerodynamics objects +- `label_list`: List of labels for each configuration + +# Keyword arguments +- `literature_path_list`: Optional paths to literature data files +- `angle_range`: Range of angles [°] +- `angle_type`: "angle_of_attack" or "side_slip" (default: angle_of_attack) +- `angle_of_attack`: AoA [rad] (default: 0.0) +- `side_slip`: Side slip angle [rad] (default: 0.0) +- `v_a`: Wind speed [m/s] (default: 10.0) +- `title`: Plot title +- `data_type`: File extension (default: ".pdf") +- `save_path`: Path to save (default: nothing) +- `is_save`: Whether to save (default: true) +- `is_show`: Whether to display (default: true) +- `use_tex`: Ignored for Makie (default: false) +""" +function VortexStepMethod.plot_polars( + solver_list, + body_aero_list, + label_list; + literature_path_list=String[], + angle_range=range(0, 20, 2), + angle_type="angle_of_attack", + angle_of_attack=0.0, + side_slip=0.0, + v_a=10.0, + title="polar", + data_type=".pdf", + save_path=nothing, + is_save=true, + is_show=true, + use_tex=false +) + # Validate inputs + total_cases = length(body_aero_list) + length(literature_path_list) + if total_cases != length(label_list) || length(solver_list) != length(body_aero_list) + throw(ArgumentError("Mismatch in solvers/cases/labels")) + end + main_title = replace(title, " " => "_") + + # Generate polar data + polar_data_list = [] + labels_with_re = copy(label_list) + for (i, (solver, body_aero)) in enumerate(zip(solver_list, body_aero_list)) + polar_data, rey = generate_polar_data_makie( + solver, body_aero, angle_range; + angle_type, angle_of_attack, side_slip, v_a + ) + push!(polar_data_list, polar_data) + labels_with_re[i] = "$(label_list[i]) Re = $(round(Int64, rey*1e-5))e5" + end + + # Load literature data if provided + if !isempty(literature_path_list) + for path in literature_path_list + data = readdlm(path, ',') + header = lowercase.(string.(data[1, :])) + alpha_idx = findfirst(x -> occursin("alpha", x), header) + cl_idx = findfirst(x -> occursin("cl", x), header) + cd_idx = findfirst(x -> occursin("cd", x), header) + cs_idx = findfirst(x -> occursin("cs", x), header) + cs_col = cs_idx === nothing ? zeros(size(data, 1)-1) : data[2:end, cs_idx] + push!(polar_data_list, [ + data[2:end, alpha_idx], + data[2:end, cl_idx], + data[2:end, cd_idx], + cs_col + ]) + end + end + + # Create figure with 2x2 grid + fig = Figure(size=(1400, 1400)) + + ax_cl = Axis(fig[1, 1], title="CL vs $angle_type [°]", + xlabel="$angle_type [°]", ylabel="CL") + ax_cd = Axis(fig[1, 2], title="CD vs $angle_type [°]", + xlabel="$angle_type [°]", ylabel="CD") + ax_cs = Axis(fig[2, 1], title="CS vs $angle_type [°]", + xlabel="$angle_type [°]", ylabel="CS") + ax_polar = Axis(fig[2, 2], title="CL vs CD", + xlabel="CD", ylabel="CL") + + # Number of computational results + n_solvers = length(solver_list) + + # Plot CL vs angle + for (i, (polar_data, label)) in enumerate(zip(polar_data_list, labels_with_re)) + marker = i <= n_solvers ? :star5 : :circle + markersize = i <= n_solvers ? 12 : 8 + scatterlines!(ax_cl, polar_data[1], polar_data[2]; + label=label, marker=marker, markersize=markersize) + if maximum(polar_data[2]) > 10 + ylims!(ax_cl, -0.5, 2) + end + end + axislegend(ax_cl, position=:lt) + + # Plot CD vs angle + for (i, (polar_data, label)) in enumerate(zip(polar_data_list, labels_with_re)) + marker = i <= n_solvers ? :star5 : :circle + markersize = i <= n_solvers ? 12 : 8 + scatterlines!(ax_cd, polar_data[1], polar_data[3]; + label=label, marker=marker, markersize=markersize) + if maximum(polar_data[2]) > 10 + ylims!(ax_cd, -0.5, 2) + end + end + axislegend(ax_cd, position=:lt) + + # Plot CS vs angle + for (i, (polar_data, label)) in enumerate(zip(polar_data_list, labels_with_re)) + marker = i <= n_solvers ? :star5 : :circle + markersize = i <= n_solvers ? 12 : 8 + scatterlines!(ax_cs, polar_data[1], polar_data[4]; + label=label, marker=marker, markersize=markersize) + if maximum(polar_data[2]) > 10 + ylims!(ax_cs, -0.5, 2) + end + end + axislegend(ax_cs, position=:lt) + + # Plot CL vs CD + for (i, (polar_data, label)) in enumerate(zip(polar_data_list, labels_with_re)) + marker = i <= n_solvers ? :star5 : :circle + markersize = i <= n_solvers ? 12 : 8 + scatterlines!(ax_polar, polar_data[3], polar_data[2]; + label=label, marker=marker, markersize=markersize) + if maximum(polar_data[2]) > 10 || maximum(polar_data[3]) > 10 + ylims!(ax_polar, -0.5, 2) + xlims!(ax_polar, -0.5, 2) + end + end + axislegend(ax_polar, position=:lt) + + # Save and show + if is_save && !isnothing(save_path) + save_plot(fig, save_path, main_title; data_type) + end + + if is_show + display(fig) + end + + return fig +end + +""" + plot_polar_data(body_aero::BodyAerodynamics; + alphas=collect(deg2rad.(-5:0.3:25)), + delta_tes=collect(deg2rad.(-5:0.3:25)), + is_show=true, use_tex=false) + +Plot polar data (Cl, Cd, Cm) as 3D surfaces using Makie. + +# Arguments +- `body_aero`: Wing aerodynamics struct + +# Keyword arguments +- `alphas`: Range of AoA values [rad] (default: -5° to 25° in 0.3° steps) +- `delta_tes`: Range of trailing edge angles [rad] (default: -5° to 25° in 0.3° steps) +- `is_show`: Whether to display (default: true) +- `use_tex`: Ignored for Makie (default: false) +""" +function VortexStepMethod.plot_polar_data(body_aero::BodyAerodynamics; + alphas=collect(deg2rad.(-5:0.3:25)), + delta_tes=collect(deg2rad.(-5:0.3:25)), + is_show=true, + use_tex=false) + + if body_aero.panels[1].aero_model == POLAR_MATRICES + # Create figure with 3 subplots + fig = Figure(size=(1500, 600)) + + # Get interpolation functions + interp_data = [ + (body_aero.panels[1].cl_interp, "Cl"), + (body_aero.panels[1].cd_interp, "Cd"), + (body_aero.panels[1].cm_interp, "Cm") + ] + + # Create each subplot + for (idx, (interp, label)) in enumerate(interp_data) + ax = Axis3(fig[1, idx]; + title="$label vs α and δ", + xlabel="δ [rad]", + ylabel="α [rad]", + zlabel=label, + azimuth=1.275*π) + + # Create interpolation matrix + interp_matrix = [interp(alpha, delta_te) + for alpha in alphas, delta_te in delta_tes] + + # Create wireframe + wireframe!(ax, delta_tes, alphas, interp_matrix; + color=:blue, linewidth=0.5, transparency=true) + end + + if is_show + display(fig) + end + return fig + else + throw(ArgumentError( + "Plotting polar data for $(body_aero.panels[1].aero_model) not implemented." + )) + end +end + end diff --git a/src/panel.jl b/src/panel.jl index 23886f4d..ba2e2448 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -535,4 +535,4 @@ function calculate_velocity_induced_bound_2D!( cross_square .= cross_.^2 U_2D .= (cross_ ./ sum(cross_square) ./ 2π) .* norm(r0) return nothing -end \ No newline at end of file +end diff --git a/src/plotting.jl b/src/plotting.jl deleted file mode 100644 index 3c602edb..00000000 --- a/src/plotting.jl +++ /dev/null @@ -1,920 +0,0 @@ - -""" - set_plot_style(titel_size=16; use_tex=false) - -Set the default style for plots using LaTeX. -`` -# Arguments: -- `titel_size`: size of the plot title in points (default: 16) -- `ùse_tex`: if the external `pdflatex` command shall be used -""" -function set_plot_style(titel_size=16; use_tex=false) - rcParams = plt.PyDict(plt.matplotlib."rcParams") - rcParams["text.usetex"] = use_tex - rcParams["font.family"] = "serif" - if use_tex - rcParams["font.serif"] = ["Computer Modern Roman"] - end - rcParams["axes.titlesize"] = titel_size - rcParams["axes.labelsize"] = 12 - rcParams["axes.linewidth"] = 1 - rcParams["lines.linewidth"] = 1 - rcParams["lines.markersize"] = 6 - rcParams["xtick.labelsize"] = 10 - rcParams["ytick.labelsize"] = 10 - rcParams["legend.fontsize"] = 10 - rcParams["figure.titlesize"] = 16 - if use_tex - rcParams["pgf.texsystem"] = "pdflatex" # Use pdflatex - end - rcParams["pgf.rcfonts"] = false - rcParams["figure.figsize"] = (10, 6) # Default figure size -end - - -""" - save_plot(fig, save_path, title; data_type=".pdf") - -Save a plot to a file. - -# Arguments -- `fig`: Plots figure object -- `save_path`: Path to save the plot -- `title`: Title of the plot - -# Keyword arguments -- `data_type`: File extension (default: ".pdf") -""" -function VortexStepMethod.save_plot(fig, save_path, title; data_type=".pdf") - isnothing(save_path) && throw(ArgumentError("save_path should be provided")) - - !isdir(save_path) && mkpath(save_path) - full_path = joinpath(save_path, title * data_type) - - @debug "Attempting to save figure to: $full_path" - @debug "Current working directory: $(pwd())" - - try - fig.savefig(full_path) - @debug "Figure saved as $data_type" - - if isfile(full_path) - @debug "File successfully saved to $full_path" - @debug "File size: $(filesize(full_path)) bytes" - else - @info "File does not exist after save attempt: $full_path" - end - catch e - @error "Error saving figure: $e" - @error "Error type: $(typeof(e))" - rethrow(e) - end -end - -""" - show_plot(fig; dpi=130) - -Display a plot at specified DPI. - -# Arguments -- `fig`: Plots figure object - -# Keyword arguments -- `dpi`: Dots per inch for the figure (default: 130) -""" -function VortexStepMethod.show_plot(fig; dpi=130) - plt.display(fig) -end - -""" - plot_line_segment!(ax, segment, color, label; width=3) - -Plot a line segment in 3D with arrow. - -# Arguments -- `ax`: Plot axis -- `segment`: Array of two points defining the segment -- `color`: Color of the segment -- `label`: Label for the legend - -# Keyword Arguments -- `width`: Line width (default: 3) -""" -function plot_line_segment!(ax, segment, color, label; width=3) - ax.plot( - [segment[1][1], segment[2][1]], - [segment[1][2], segment[2][2]], - [segment[1][3], segment[2][3]], - color=color, label=label, linewidth=width - ) - - dir = segment[2] - segment[1] - ax.quiver( - [segment[1][1]], [segment[1][2]], [segment[1][3]], - [dir[1]], [dir[2]], [dir[3]], - color=color - ) -end - -""" - set_axes_equal!(ax; zoom=1.8) - -Set 3D plot axes to equal scale. - -# Arguments -- ax: 3D plot axis - -# Keyword arguments -zoom: zoom factor (default: 1.8) -""" -function set_axes_equal!(ax; zoom=1.8) - x_lims = ax.get_xlim3d() ./ zoom - y_lims = ax.get_ylim3d() ./ zoom - z_lims = ax.get_zlim3d() ./ zoom - - x_range = abs(x_lims[2] - x_lims[1]) - y_range = abs(y_lims[2] - y_lims[1]) - z_range = abs(z_lims[2] - z_lims[1]) - - max_range = max(x_range, y_range, z_range) - - x_mid = mean(x_lims) - y_mid = mean(y_lims) - z_mid = mean(z_lims) - - ax.set_xlim3d((x_mid - max_range / 2, x_mid + max_range / 2)) - ax.set_ylim3d((y_mid - max_range / 2, y_mid + max_range / 2)) - ax.set_zlim3d((z_mid - max_range / 2, z_mid + max_range / 2)) -end - -""" - create_geometry_plot(body_aero::BodyAerodynamics, title, view_elevation, view_azimuth; - zoom=1.8, use_tex=false) - -Create a 3D plot of wing geometry including panels and filaments. - -# Arguments -- body_aero: struct of type BodyAerodynamics -- title: plot title -- view_elevation: initial view elevation angle [°] -- view_azimuth: initial view azimuth angle [°] - -# Keyword arguments -- zoom: zoom factor (default: 1.8) -""" -function create_geometry_plot(body_aero::BodyAerodynamics, title, view_elevation, view_azimuth; - zoom=1.8, use_tex=false) - set_plot_style(28; use_tex) - - panels = body_aero.panels - va = isa(body_aero.va, Tuple) ? body_aero.va[1] : body_aero.va - - # Extract geometric data - corner_points = [panel.corner_points for panel in panels] - control_points = [panel.control_point for panel in panels] - aero_centers = [panel.aero_center for panel in panels] - - # Create plot - fig = plt.figure(figsize=(14, 14)) - ax = fig.add_subplot(111, projection="3d") - ax.set_title(title) - - # Plot panels - legend_used = Dict{String,Bool}() - for (i, panel) in enumerate(panels) - # Plot panel edges and surfaces - corners = Matrix{Float64}(panel.corner_points) - x_corners = corners[1, :] - y_corners = corners[2, :] - z_corners = corners[3, :] - - push!(x_corners, x_corners[1]) - push!(y_corners, y_corners[1]) - push!(z_corners, z_corners[1]) - - ax.plot(x_corners, - y_corners, - z_corners, - color=:grey, - linewidth=1, - label=i == 1 ? "Panel Edges" : "") - - # Plot control points and aerodynamic centers - ax.scatter([control_points[i][1]], [control_points[i][2]], [control_points[i][3]], - color=:green, label=i == 1 ? "Control Points" : "") - ax.scatter([aero_centers[i][1]], [aero_centers[i][2]], [aero_centers[i][3]], - color=:blue, label=i == 1 ? "Aerodynamic Centers" : "") - - # Plot filaments - filaments = calculate_filaments_for_plotting(panel) - legends = ["Bound Vortex", "side1", "side2", "wake_1", "wake_2"] - - for (filament, legend) in zip(filaments, legends) - x1, x2, color = filament - @debug "Legend: $legend" - show_legend = !get(legend_used, legend, false) - plot_line_segment!(ax, [x1, x2], color, show_legend ? legend : "") - legend_used[legend] = true - end - end - - # Plot velocity vector - max_chord = maximum(panel.chord for panel in panels) - va_mag = norm(va) - va_vector_begin = -2 * max_chord * va / va_mag - va_vector_end = va_vector_begin + 1.5 * va / va_mag - plot_line_segment!(ax, [va_vector_begin, va_vector_end], :lightblue, "va") - - # Add legends for the first occurrence of each label - handles, labels = ax.get_legend_handles_labels() - # by_label = Dict(zip(labels, handles)) - # ax.legend(values(by_label), keys(by_label), bbox_to_anchor=(0, 0, 1.1, 1)) - - # Set labels and make axes equal - ax.set_xlabel("x") - ax.set_ylabel("y") - ax.set_zlabel("z") - set_axes_equal!(ax; zoom) - - # Set the initial view - ax.view_init(elev=view_elevation, azim=view_azimuth) - - # Ensure the figure is fully rendered - # fig.canvas.draw() - plt.tight_layout(rect=(0, 0, 1, 0.97)) - - return fig -end - -""" - plot_geometry(body_aero::BodyAerodynamics, title; - data_type=".pdf", save_path=nothing, - is_save=false, is_show=false, - view_elevation=15, view_azimuth=-120, use_tex=false) - -Plot wing geometry from different viewpoints and optionally save/show plots. - -# Arguments: -- `body_aero`: the [BodyAerodynamics](@ref) to plot -- title: plot title - -# Keyword arguments: -- `data_type``: string with the file type postfix (default: ".pdf") -- `save_path`: path for saving the graphic (default: `nothing`) -- `is_save`: boolean value, indicates if the graphic shall be saved (default: `false`) -- `is_show`: boolean value, indicates if the graphic shall be displayed (default: `false`) -- `view_elevation`: initial view elevation angle (default: 15) [°] -- `view_azimuth`: initial view azimuth angle (default: -120) [°] -- `use_tex`: if the external `pdflatex` command shall be used (default: false) - -""" -function VortexStepMethod.plot_geometry(body_aero::BodyAerodynamics, title; - data_type=".pdf", - save_path=nothing, - is_save=false, - is_show=false, - view_elevation=15, - view_azimuth=-120, - use_tex=false) - - if is_save - plt.ioff() - # Angled view - fig = create_geometry_plot(body_aero, "$(title)_angled_view", 15, -120; use_tex) - save_plot(fig, save_path, "$(title)_angled_view", data_type=data_type) - - # Top view - fig = create_geometry_plot(body_aero, "$(title)_top_view", 90, 0; use_tex) - save_plot(fig, save_path, "$(title)_top_view", data_type=data_type) - - # Front view - fig = create_geometry_plot(body_aero, "$(title)_front_view", 0, 0; use_tex) - save_plot(fig, save_path, "$(title)_front_view", data_type=data_type) - - # Side view - fig = create_geometry_plot(body_aero, "$(title)_side_view", 0, -90; use_tex) - save_plot(fig, save_path, "$(title)_side_view", data_type=data_type) - end - - if is_show - plt.ion() - fig = create_geometry_plot(body_aero, title, view_elevation, view_azimuth; use_tex) - plt.display(fig) - else - fig = create_geometry_plot(body_aero, title, view_elevation, view_azimuth; use_tex) - end - fig -end - -""" - plot_distribution(y_coordinates_list, results_list, label_list; - title="spanwise_distribution", data_type=".pdf", - save_path=nothing, is_save=false, is_show=true, use_tex=false) - -Plot spanwise distributions of aerodynamic properties. - -# Arguments -- `y_coordinates_list`: List of spanwise coordinates -- `results_list`: List of result dictionaries -- `label_list`: List of labels for different results - -# Keyword arguments -- `title`: Plot title (default: "spanwise_distribution") -- `data_type`: File extension for saving (default: ".pdf") -- `save_path`: Path to save plots (default: nothing) -- `is_save`: Whether to save plots (default: false) -- `is_show`: Whether to display plots (default: true) -- `use_tex`: if the external `pdflatex` command shall be used -""" -function VortexStepMethod.plot_distribution(y_coordinates_list, results_list, label_list; - title="spanwise_distribution", - data_type=".pdf", - save_path=nothing, - is_save=false, - is_show=true, - use_tex=false) - - length(results_list) == length(label_list) || throw(ArgumentError( - "Number of results ($(length(results_list))) must match number of labels ($(length(label_list)))" - )) - - # Set the plot style - set_plot_style(; use_tex) - - # Initializing plot - fig, axs = plt.subplots(3, 3, figsize=(16, 10)) - fig.suptitle(title, fontsize=16) - - # CL plot - for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list) - value = "$(round(result_i["cl"], digits=2))" - if label_i == "LLT" - label = label_i * L" $~C_\mathrm{L}$: " * value - else - label = label_i * L" $C_\mathrm{L}$: " * value - end - axs[1, 1].plot( - y_coordinates_i, - result_i["cl_distribution"], - label=label - ) - end - axs[1, 1].set_title(L"$C_\mathrm{L}$ Distribution", size=16) - axs[1, 1].set_xlabel(L"Spanwise Position $y/b$") - axs[1, 1].set_ylabel(L"Lift Coefficient $C_\mathrm{L}$") - axs[1, 1].legend() - - # CD plot - for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list) - value = "$(round(result_i["cl"], digits=2))" - if label_i == "LLT" - label = label_i * L" $~C_\mathrm{D}$: " * value - else - label = label_i * L" $C_\mathrm{D}$: " * value - end - axs[1, 2].plot( - y_coordinates_i, - result_i["cd_distribution"], - label=label - ) - end - axs[1, 2].set_title(L"$C_\mathrm{D}$ Distribution", size=16) - axs[1, 2].set_xlabel(L"Spanwise Position $y/b$") - axs[1, 2].set_ylabel(L"Drag Coefficient $C_\mathrm{D}$") - axs[1, 2].legend() - - # Gamma Distribution - for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list) - axs[1, 3].plot( - y_coordinates_i, - result_i["gamma_distribution"], - label=label_i - ) - end - axs[1, 3].set_title(L"\Gamma~Distribution", size=16) - axs[1, 3].set_xlabel(L"Spanwise Position $y/b$") - axs[1, 3].set_ylabel(L"Circulation~\Gamma") - axs[1, 3].legend() - - # Geometric Alpha - for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list) - axs[2, 1].plot( - y_coordinates_i, - result_i["alpha_geometric"], - label=label_i - ) - end - axs[2, 1].set_title(L"$\alpha$ Geometric", size=16) - axs[2, 1].set_xlabel(L"Spanwise Position $y/b$") - axs[2, 1].set_ylabel(L"Angle of Attack $\alpha$ (deg)") - axs[2, 1].legend() - - # Calculated/ Corrected Alpha - for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list) - axs[2, 2].plot( - y_coordinates_i, - result_i["alpha_at_ac"], - label=label_i - ) - end - axs[2, 2].set_title(L"$\alpha$ result (corrected to aerodynamic center)", size=16) - axs[2, 2].set_xlabel(L"Spanwise Position $y/b$") - axs[2, 2].set_ylabel(L"Angle of Attack $\alpha$ (deg)") - axs[2, 2].legend() - - # Uncorrected Alpha plot - for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list) - axs[2, 3].plot( - y_coordinates_i, - result_i["alpha_uncorrected"], - label=label_i - ) - end - axs[2, 3].set_title(L"$\alpha$ Uncorrected (if VSM, at the control point)", size=16) - axs[2, 3].set_xlabel(L"Spanwise Position $y/b$") - axs[2, 3].set_ylabel(L"Angle of Attack $\alpha$ (deg)") - axs[2, 3].legend() - - # Force Components - for (idx, component) in enumerate(["x", "y", "z"]) - axs[3, idx].set_title("Force in $component direction", size=16) - axs[3, idx].set_xlabel(L"Spanwise Position $y/b$") - axs[3, idx].set_ylabel(raw"$F_\mathrm" * "{$component}" * raw"$") - for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) - # Extract force components for the current direction (idx) - forces = results["F_distribution"][idx, :] - # Verify dimensions match - if length(y_coords) != length(forces) - @warn "Dimension mismatch in force plotting" length(y_coords) length(forces) component - continue # Skip this component instead of throwing error - end - space = "" - if label == "LLT" - space = "~" - end - axs[3, idx].plot( - y_coords, - forces, - label="$label" * space * raw"$~\Sigma~F_\mathrm" * "{$component}:~" * - raw"$" * "$(round(results["F$component"], digits=2)) N" - ) - axs[3, idx].legend() - end - end - - fig.tight_layout() - - # Save and show plot - if is_save - save_plot(fig, save_path, title, data_type=data_type) - end - - if is_show - show_plot(fig) - end - - return fig -end - -""" - generate_polar_data(solver, body_aero::BodyAerodynamics, angle_range; - angle_type="angle_of_attack", angle_of_attack=0.0, - side_slip=0.0, v_a=10.0, use_latex=false) - -Generate polar data for aerodynamic analysis over a range of angles. - -# Arguments -- `solver`: Aerodynamic solver object -- `body_aero`: Wing aerodynamics struct -- `angle_range`: Range of angles to analyze - -# Keyword arguments -- `angle_type`: Type of angle variation ("angle_of_attack" or "side_slip") -- `angle_of_attack`: Initial angle of attack [rad] -- `side_slip`: Initial side slip angle in [rad] -- `v_a`: norm of apparent wind speed [m/s] - -# Returns -- Tuple of polar data array and Reynolds number -""" -function generate_polar_data( - solver, - body_aero::BodyAerodynamics, - angle_range; - angle_type="angle_of_attack", - angle_of_attack=0.0, - side_slip=0.0, - v_a=10.0, - use_latex=false -) - n_panels = length(body_aero.panels) - n_angles = length(angle_range) - - # Initialize arrays - cl = zeros(n_angles) - cd = zeros(n_angles) - cs = zeros(n_angles) - gamma_distribution = zeros(n_angles, n_panels) - cl_distribution = zeros(n_angles, n_panels) - cd_distribution = zeros(n_angles, n_panels) - cs_distribution = zeros(n_angles, n_panels) - reynolds_number = zeros(n_angles) - - # Previous gamma for initialization - gamma = nothing - - for (i, angle_i) in enumerate(angle_range) - # Set angle based on type - if angle_type == "angle_of_attack" - α = deg2rad(angle_i) - β = side_slip - elseif angle_type == raw"side_slip" - α = angle_of_attack - β = deg2rad(angle_i) - else - throw(ArgumentError("angle_type must be 'angle_of_attack' or 'side_slip'")) - end - - # Update inflow conditions - set_va!( - body_aero, - [ - cos(α) * cos(β), - sin(β), - sin(α) - ] * v_a - ) - - # Solve and store results - results = solve(solver, body_aero, gamma_distribution[i, :]) - - cl[i] = results["cl"] - cd[i] = results["cd"] - cs[i] = results["cs"] - gamma_distribution[i, :] = results["gamma_distribution"] - cl_distribution[i, :] = results["cl_distribution"] - cd_distribution[i, :] = results["cd_distribution"] - cs_distribution[i, :] = results["cs_distribution"] - reynolds_number[i] = results["Rey"] - - # Store gamma for next iteration - gamma = gamma_distribution[i, :] - end - - polar_data = [ - angle_range, - cl, - cd, - cs, - gamma_distribution, - cl_distribution, - cd_distribution, - cs_distribution, - reynolds_number - ] - - return polar_data, reynolds_number[1] -end - -""" - plot_polars(solver_list, body_aero_list, label_list; - literature_path_list=String[], - angle_range=range(0, 20, 2), angle_type="angle_of_attack", - angle_of_attack=0.0, side_slip=0.0, v_a=10.0, - title="polar", data_type=".pdf", save_path=nothing, - is_save=true, is_show=true, use_tex=false) - -Plot polar data comparing different solvers and configurations. - -# Arguments -- `solver_list`: List of aerodynamic solvers -- `body_aero_list`: List of wing aerodynamics objects -- `label_list`: List of labels for each configuration - -# Keyword arguments -- `literature_path_list`: Optional paths to literature data files -- `angle_range`: Range of angles to analyze [°] -- `angle_type`: "`angle_of_attack`" or "`side_slip`"; (default: `angle_of_attack`) -- `angle_of_attack:` AoA to be used for plotting the polars (default: 0.0) [rad] -- `side_slip`: side slip angle (default: 0.0) [rad] -- v_a: norm of apparent wind speed (default: 10.0) [m/s] -- title: plot title -- `data_type`: File extension for saving (default: ".pdf") -- `save_path`: Path to save plots (default: nothing) -- `is_save`: Whether to save plots (default: true) -- `is_show`: Whether to display plots (default: true) -- `use_tex`: if the external `pdflatex` command shall be used (default: false) -""" -function VortexStepMethod.plot_polars( - solver_list, - body_aero_list, - label_list; - literature_path_list=String[], - angle_range=range(0, 20, 2), - angle_type="angle_of_attack", - angle_of_attack=0.0, - side_slip=0.0, - v_a=10.0, - title="polar", - data_type=".pdf", - save_path=nothing, - is_save=true, - is_show=true, - use_tex=false -) - # Validate inputs - total_cases = length(body_aero_list) + length(literature_path_list) - if total_cases != length(label_list) || length(solver_list) != length(body_aero_list) - throw(ArgumentError("Mismatch in number of solvers ($(length(solver_list))), " * - "cases ($total_cases), and labels ($(length(label_list)))")) - end - main_title = replace(title, " " => "_") - set_plot_style(; use_tex) - - # Generate polar data - polar_data_list = [] - for (i, (solver, body_aero)) in enumerate(zip(solver_list, body_aero_list)) - polar_data, rey = generate_polar_data( - solver, body_aero, angle_range; - angle_type, - angle_of_attack, - side_slip, - v_a - ) - push!(polar_data_list, polar_data) - # Update label with Reynolds number - label_list[i] = "$(label_list[i]) Re = $(round(Int64, rey*1e-5))e5" - end - # Load literature data if provided - if !isempty(literature_path_list) - for path in literature_path_list - data = readdlm(path, ',') - header = lowercase.(string.(data[1, :])) - # Find column indices for alpha, CL, CD, CS (case-insensitive, allow common variants) - alpha_idx = findfirst(x -> occursin("alpha", x), header) - cl_idx = findfirst(x -> occursin("cl", x), header) - cd_idx = findfirst(x -> occursin("cd", x), header) - cs_idx = findfirst(x -> occursin("cs", x), header) - # Fallback: if CS not found, fill with zeros - cs_col = cs_idx === nothing ? zeros(size(data, 1)-1) : data[2:end, cs_idx] - # Push as [alpha, CL, CD, CS] - push!(polar_data_list, [ - data[2:end, alpha_idx], - data[2:end, cl_idx], - data[2:end, cd_idx], - cs_col - ]) - end - end - - # Initializing plot - fig, axs = plt.subplots(2, 2, figsize=(14, 14)) - - # Number of computational results (excluding literature) - n_solvers = length(solver_list) - for (i, (polar_data, label)) in enumerate(zip(polar_data_list, label_list)) - if i < n_solvers - linestyle = "-" - marker = "*" - markersize = 7 - else - linestyle = "-" - marker = "." - markersize = 5 - end - if contains(label, "LLT") - label = replace(label, "e5" => raw"\cdot10^5") - label = replace(label, " " => raw"~") - label = replace(label, "LLT" => raw"\mathrm{LLT}{~\,}") - label = raw"$" * label * raw"$" - else - label = replace(label, "e5" => raw"\cdot10^5") - label = replace(label, " " => "~") - label = replace(label, "VSM" => raw"\mathrm{VSM}") - label = raw"$" * label * raw"$" - end - axs[1, 1].plot( - polar_data[1], - polar_data[2], - label=label, - linestyle=linestyle, - marker=marker, - markersize=markersize, - ) - # Limit y-range if CL > 10 - if maximum(polar_data[2]) > 10 - axs[1, 1].set_ylim([-0.5, 2]) - end - title = raw"$C_\mathrm{L}" * raw"$" * " vs $angle_type [°]" - axs[1, 1].set_title(title) - axs[1, 1].set_xlabel("$angle_type [°]") - axs[1, 1].set_ylabel(L"$C_\mathrm{L}$") - axs[1, 1].legend() - end - - for (i, (polar_data, label)) in enumerate(zip(polar_data_list, label_list)) - if i < n_solvers - linestyle = "-" - marker = "*" - markersize = 7 - else - linestyle = "-" - marker = "." - markersize = 5 - end - if contains(label, "LLT") - label = replace(label, "e5" => raw"\cdot10^5") - label = replace(label, " " => raw"~") - label = replace(label, "LLT" => raw"\mathrm{LLT}{~\,}") - label = raw"$" * label * raw"$" - else - label = replace(label, "e5" => raw"\cdot10^5") - label = replace(label, " " => "~") - label = replace(label, "VSM" => raw"\mathrm{VSM}") - label = raw"$" * label * raw"$" - end - axs[1, 2].plot( - polar_data[1], - polar_data[3], - label=label, - linestyle=linestyle, - marker=marker, - markersize=markersize, - ) - # Limit y-range if CL > 10 - if maximum(polar_data[2]) > 10 - axs[1, 2].set_ylim([-0.5, 2]) - end - title = raw"$C_\mathrm{D}" * raw"$" * " vs $angle_type [°]" - axs[1, 2].set_title(title) - axs[1, 2].set_xlabel("$angle_type [°]") - axs[1, 2].set_ylabel(L"$C_\mathrm{D}$") - axs[1, 2].legend() - end - - - for (i, (polar_data, label)) in enumerate(zip(polar_data_list, label_list)) - if i < n_solvers - linestyle = "-" - marker = "*" - markersize = 7 - else - linestyle = "-" - marker = "." - markersize = 5 - end - if contains(label, "LLT") - label = replace(label, "e5" => raw"\cdot10^5") - label = replace(label, " " => raw"~") - label = replace(label, "LLT" => raw"\mathrm{LLT}{~\,}") - label = raw"$" * label * raw"$" - else - label = replace(label, "e5" => raw"\cdot10^5") - label = replace(label, " " => "~") - label = replace(label, "VSM" => raw"\mathrm{VSM}") - label = raw"$" * label * raw"$" - end - axs[2, 1].plot( - polar_data[1], - polar_data[4], - label=label, - linestyle=linestyle, - marker=marker, - markersize=markersize, - ) - # Limit y-range if CL > 10 - if maximum(polar_data[2]) > 10 - axs[2, 1].set_ylim([-0.5, 2]) - end - title = raw"$C_\mathrm{S}" * raw"$" * " vs $angle_type [°]" - axs[2, 1].set_title(title) - axs[2, 1].set_xlabel("$angle_type [°]") - axs[2, 1].set_ylabel(L"$C_\mathrm{S}$") - axs[2, 1].legend() - end - - for (i, (polar_data, label)) in enumerate(zip(polar_data_list, label_list)) - if i < n_solvers - linestyle = "-" - marker = "*" - markersize = 7 - else - linestyle = "-" - marker = "." - markersize = 5 - end - if contains(label, "LLT") - label = replace(label, "e5" => raw"\cdot10^5") - label = replace(label, " " => raw"~") - label = replace(label, "LLT" => raw"\mathrm{LLT}{~\,}") - label = raw"$" * label * raw"$" - else - label = replace(label, "e5" => raw"\cdot10^5") - label = replace(label, " " => "~") - label = replace(label, "VSM" => raw"\mathrm{VSM}") - label = raw"$" * label * raw"$" - end - axs[2, 2].plot( - polar_data[3], - polar_data[2], - label=label, - linestyle=linestyle, - marker=marker, - markersize=markersize, - ) - # Limit y-range if CL > 10 - if maximum(polar_data[2]) > 10 || maximum(polar_data[3]) > 10 - axs[2, 2].set_ylim([-0.5, 2]) - axs[2, 2].set_xlim([-0.5, 2]) - end - title = raw"$C_\mathrm{L}" * raw"$" * " vs " * raw"$C_\mathrm{D}" * raw"$" - axs[2, 2].set_title(title) - axs[2, 2].set_xlabel(L"$C_\mathrm{D}$") - axs[2, 2].set_ylabel(L"$C_\mathrm{L}$") - axs[2, 2].legend() - end - - fig.tight_layout(h_pad=3.5, rect=(0.01, 0.01, 0.99, 0.99)) - - # Save and show plot - if is_save && !isnothing(save_path) - save_plot(fig, save_path, main_title; data_type) - end - - if is_show - show_plot(fig) - end - - return fig -end - -""" - plot_polar_data(body_aero::BodyAerodynamics; alphas=collect(deg2rad.(-5:0.3:25)), delta_tes=collect(deg2rad.(-5:0.3:25))) - -Plot polar data (Cl, Cd, Cm) as 3D surfaces against alpha and delta_te angles. delta_te is the trailing edge deflection angle -relative to the 2d airfoil or panel chord line. - -# Arguments -- `body_aero`: Wing aerodynamics struct - -# Keyword arguments -- `alphas`: Range of angle of attack values in radians (default: -5° to 25° in 0.3° steps) -- `delta_tes`: Range of trailing edge angles in radians (default: -5° to 25° in 0.3° steps) -- `is_show`: Whether to display plots (default: true) -- `use_tex`: if the external `pdflatex` command shall be used -""" -function VortexStepMethod.plot_polar_data(body_aero::BodyAerodynamics; - alphas=collect(deg2rad.(-5:0.3:25)), - delta_tes = collect(deg2rad.(-5:0.3:25)), - is_show = true, - use_tex = false - ) - if body_aero.panels[1].aero_model == POLAR_MATRICES - set_plot_style() - - # Create figure with subplots - fig = plt.figure(figsize=(15, 6)) - - # Get interpolation functions and labels - interp_data = [ - (body_aero.panels[1].cl_interp, L"$C_l$"), - (body_aero.panels[1].cd_interp, L"$C_d$"), - (body_aero.panels[1].cm_interp, L"$C_m$") - ] - - # Create each subplot - for (idx, (interp, label)) in enumerate(interp_data) - ax = fig.add_subplot(1, 3, idx, projection="3d") - - # Create interpolation matrix - interp_matrix = zeros(length(alphas), length(delta_tes)) - interp_matrix .= [interp(alpha, delta_te) for alpha in alphas, delta_te in delta_tes] - X = collect(delta_tes) .+ zeros(length(alphas))' - Y = collect(alphas)' .+ zeros(length(delta_tes)) - - # Plot surface - ax.plot_wireframe(X, Y, interp_matrix, - edgecolor="blue", - lw=0.5, - rstride=5, - cstride=5, - alpha=0.6) - - # Set labels and title - ax.set_xlabel(L"$\delta$ [rad]") - ax.set_ylabel(L"$\alpha$ [rad]") - ax.set_zlabel(label) - ax.set_title(label * L" vs $\alpha$ and $\delta$") - ax.grid(true) - end - - # Adjust layout and display - plt.tight_layout(rect=(0.01, 0.01, 0.99, 0.99)) - if is_show - show_plot(fig) - end - return fig - else - throw(ArgumentError("Plotting polar data for $(body_aero.panels[1].aero_model) is not implemented.")) - end -end