Unify Wing and RamAirWing

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
 ```

examples/Project.toml

@@ -0,0 +1,3 @@
+[deps]
+ControlPlots = "23c2ee80-7a9e-4350-b264-8e670f12517c"
+VortexStepMethod = "ed3cd733-9f0f-46a9-93e0-89b8d4998dd9"

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])

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

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")

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)
@@ -129,20 +129,20 @@ 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
 
         # 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
@@ -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,31 +283,46 @@ 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
                 break
             end
             cl_old = cl
         end
-        
-        push!(stall_angles, panel_stall)
+
+        stall_angles[idx] = panel_stall
     end
-    
-    return stall_angles
+
+    return nothing
 end
 
 """