From f078c77e32118f36476cd565b5513ac031f13f15 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sat, 25 Oct 2025 22:36:47 +0200 Subject: [PATCH 1/6] 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 2/6] 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 3/6] 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 4/6] 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 2b0301504a3ac6a2c14054df13f28bc6ee5713e8 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Sat, 1 Nov 2025 23:31:49 +0100 Subject: [PATCH 5/6] 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 6/6] 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