add non-allocating update pos functions

Author: 1-Bart-1

Project.toml

@@ -13,6 +13,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 Measures = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -34,13 +35,13 @@ Timers = "0.1"
 julia = "1.10, 1.11"
 
 [extras]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Xfoil = "19641d66-a62d-11e8-2441-8f57a969a9c4"
-BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 
 [targets]
 test = ["Test", "DataFrames", "CSV", "Distributed", "Documenter", "Xfoil", "BenchmarkTools"]

examples/update_pos.jl

@@ -0,0 +1,23 @@
+using VortexStepMethod
+using VortexStepMethod: BoundFilament, MVec3, SemiInfiniteFilament, update_pos!, Panel, Section
+
+
+x1 = zeros(MVec3)
+direction = zeros(MVec3)
+vel_mag = zero(Float64)
+filament_direction = zero(Int64)
+@time filament = SemiInfiniteFilament(x1, direction, vel_mag, filament_direction)
+@time update_pos!(filament, x1, direction, vel_mag, filament_direction) 
+
+bound_point_1 = zeros(MVec3)
+bound_point_2 = zeros(MVec3)
+@time filament = BoundFilament(bound_point_1, bound_point_2)
+@time update_pos!(filament, bound_point_1, bound_point_2)
+
+LE_point = zeros(MVec3)
+TE_point = zeros(MVec3)
+aero_input = :inviscid
+@time section = Section(LE_point, TE_point, aero_input)
+@time update_pos!(section, LE_point, TE_point)
+
+

src/VortexStepMethod.jl

@@ -11,7 +11,7 @@ using Measures
 using LaTeXStrings
 using NonlinearSolve
 using Interpolations
-using Interpolations: linear_interpolation, Line, Extrapolation, FilledExtrapolation
+using Interpolations: Extrapolation
 using Serialization
 using SharedArrays
 
@@ -39,14 +39,14 @@ const MVec3    = MVector{3, Float64}
 
 Position vector, either a `MVec3` or a `Vector` for use in function signatures.
 """
-const PosVector=Union{MVec3, Vector, SizedVector{3, Float64, Vector{Float64}}}
+const PosVector=Union{MVec3, Vector}
 
 """
    const VelVector=Union{MVec3, Vector}
 
 Velocity vector, either a `MVec3` or a `Vector` for use in function signatures.
 """
-const VelVector=Union{MVec3, Vector, SizedVector{3, Float64, Vector{Float64}}}
+const VelVector=Union{MVec3, Vector}
 
 """
    Model `VSM` `LLT`

src/body_aerodynamics.jl

@@ -104,6 +104,39 @@ function Base.setproperty!(obj::BodyAerodynamics, sym::Symbol, val)
     end
 end
 
+function update_pos!(body_aero::BodyAerodynamics, le_pos::AbstractMatrix, te_pos::AbstractMatrix)
+    panels = Panel[]
+    for wing in body_aero.wings
+        section_list = refine_aerodynamic_mesh(wing)
+        n_panels_per_wing = length(section_list) - 1
+        
+        # Calculate panel properties
+        panel_props = calculate_panel_properties(
+            section_list,
+            n_panels_per_wing,
+            aero_center_location,
+            control_point_location
+        )
+        
+        # Create panels
+        for panel in body_aero.panels
+
+            update_pos!(
+                panel,
+                section_list[i],
+                section_list[i+1],
+                panel_props.aero_centers[i],
+                panel_props.control_points[i],
+                panel_props.bound_points_1[i],
+                panel_props.bound_points_2[i],
+                panel_props.x_airf[i],
+                panel_props.y_airf[i],
+                panel_props.z_airf[i]
+            )
+        end
+    end
+end
+
 """
     PanelProperties
 
@@ -650,20 +683,13 @@ Set velocity array and update wake filaments.
 - `va`: Either a velocity vector or a list of velocity vectors
 - `omega`: Turn rate around x y and z axis
 """
-function set_va!(body_aero::BodyAerodynamics, va, omega=zeros(MVec3))
-    omega = MVec3(omega)
-
-    # Add length check for va
-    if va isa Vector{Float64} && length(va) != 3 && length(va) != length(body_aero.panels)
-        throw(ArgumentError("va must be length 3 or match number of panels"))
-    end
+function set_va!(body_aero::BodyAerodynamics, va::VelVector, omega=zeros(MVec3))
+    length(va) != 3 && throw(ArgumentError("va must be length 3"))
 
     # Calculate va_distribution based on input type
-    va_distribution = if length(va) == 3 && all(omega .== 0.0)
+    va_distribution = if all(omega .== 0.0)
         repeat(reshape(va, 1, 3), length(body_aero.panels))
-    elseif length(va) == length(body_aero.panels)
-        va
-    elseif !all(omega .== 0.0) && length(va) == 3
+    elseif !all(omega .== 0.0)
         va_dist = zeros(length(body_aero.panels), 3)
 
         for wing in body_aero.wings
@@ -677,8 +703,6 @@ function set_va!(body_aero::BodyAerodynamics, va, omega=zeros(MVec3))
             end
         end
         va_dist
-    else
-        throw(ArgumentError("Invalid va distribution: length(va)=$(length(va)) ≠ n_panels=$(length(body_aero.panels))"))
     end
 
     # Update panel velocities
@@ -691,3 +715,16 @@ function set_va!(body_aero::BodyAerodynamics, va, omega=zeros(MVec3))
     body_aero._va = va
     return nothing
 end
+
+function set_va!(body_aero::BodyAerodynamics, va_distribution::Vector{VelVector}, omega=zeros(MVec3))
+    length(va) != length(body_aero.panels) && throw(ArgumentError("Length of va distribution should be equal to number of panels."))
+    
+    for (i, panel) in enumerate(body_aero.panels)
+        panel.va = va_distribution[i]
+    end
+    
+    # Update wake elements
+    body_aero.panels = frozen_wake(va_distribution, body_aero.panels)
+    body_aero._va = va
+    return nothing
+end

src/filament.jl

@@ -19,7 +19,7 @@ Represents a bound vortex filament defined by two points.
 - length: Filament length
 - r0::MVec3: Vector from x1 to x2
 """
-struct BoundFilament <: Filament
+mutable struct BoundFilament <: Filament
     x1::MVec3       # First point
     x2::MVec3       # Second point
     length::Float64 # Filament length
@@ -30,6 +30,14 @@ struct BoundFilament <: Filament
     end
 end
 
+function update_pos!(filament::BoundFilament, x1::PosVector, x2::PosVector)
+    filament.x1 .= x1
+    filament.x2 .= x2
+    filament.length = norm(x2 - x1)
+    filament.r0 .= x2 .- x1
+    return nothing
+end
+
 """
     velocity_3D_bound_vortex(filament::BoundFilament, XVP::Vector{Float64}, 
                            gamma::Float64, core_radius_fraction::Float64)
@@ -151,13 +159,21 @@ end
 
 Represents a semi-infinite vortex filament.
 """
-struct SemiInfiniteFilament <: Filament
+mutable struct SemiInfiniteFilament <: Filament
     x1::MVec3                    # Starting point
     direction::MVec3             # Direction vector
     vel_mag::Float64             # Velocity magnitude
     filament_direction::Int64    # Direction indicator (-1 or 1)
 end
 
+function update_pos!(filament::SemiInfiniteFilament, x1::PosVector, direction::PosVector, vel_mag, filament_direction)
+    filament.x1 .= x1
+    filament.direction .= direction
+    filament.vel_mag = vel_mag
+    filament.filament_direction = filament_direction
+    return nothing
+end
+
 """
     velocity_3D_trailing_vortex_semiinfinite(filament::SemiInfiniteFilament, 
                                            Vf::Vector{Float64}, XVP::Vector{Float64},

src/kite_geometry.jl

@@ -205,7 +205,7 @@ mutable struct KiteWing <: AbstractWing
     area_interp::Extrapolation
 
     function KiteWing(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_panel_distribution=LINEAR, spanwise_direction=[0.0, 1.0, 0.0])
+            n_panels=54, n_sections=n_panels+1, spanwise_panel_distribution=UNCHANGED, spanwise_direction=[0.0, 1.0, 0.0])
 
         !isapprox(spanwise_direction, [0.0, 1.0, 0.0]) && @error "Spanwise direction has to be [0.0, 1.0, 0.0]"
 
@@ -272,4 +272,4 @@ function rotate_v_around_k(v, k, θ)
     k = normalize(k)
     v_rot = v * cos(θ) + (k × v) * sin(θ)  + k * (k ⋅ v) * (1 - cos(θ))
     return v_rot
-end
+end

src/panel.jl

@@ -84,9 +84,6 @@ mutable struct Panel
             throw(ArgumentError("Both sections must have the same aero_input, not $aero_model and $aero_model_2"))
         end
 
-        # Initialize aerodynamic properties
-        cl_coeffs, cd_coeffs, cm_coeffs = zeros(3), zeros(3), zeros(3)
-        
         # Calculate width
         width = norm(bound_point_2 - bound_point_1)
 
@@ -95,9 +92,12 @@ mutable struct Panel
             BoundFilament(bound_point_2, bound_point_1),
             BoundFilament(bound_point_1, TE_point_1),
             BoundFilament(TE_point_2, bound_point_2)
-        ]
-
+            ]
+            
         cl_interp, cd_interp, cm_interp = nothing, nothing, nothing
+        
+        # Initialize aerodynamic properties
+        cl_coeffs, cd_coeffs, cm_coeffs = zeros(3), zeros(3), zeros(3)
 
         if aero_model === :lei_airfoil_breukels
             cl_coeffs, cd_coeffs, cm_coeffs = compute_lei_coefficients(section_1, section_2)
@@ -171,12 +171,9 @@ function update_pos!(
     y_airf::PosVector,
     z_airf::PosVector
 )
-    width = norm(bound_point_2 - bound_point_1)
-    filaments = [
-        BoundFilament(bound_point_2, bound_point_1),
-        BoundFilament(bound_point_1, TE_point_1),
-        BoundFilament(TE_point_2, bound_point_2)
-    ]
+    update_filament!(filaments[1], bound_point_2, bound_point_1)
+    update_filament!(filaments[2], bound_point_1, TE_point_1)
+    update_filament!(filaments[3], TE_point_2, bound_point_2)
 
     panel.TE_point_1 .= section_1.TE_point
     panel.LE_point_1 .= section_1.LE_point
@@ -194,7 +191,7 @@ function update_pos!(
     panel.x_airf .= x_airf
     panel.y_airf .= y_airf
     panel.z_airf .= z_airf
-    panel.width = width
+    panel.width = norm(bound_point_2 - bound_point_1)
     panel.filaments .= filaments
     nothing
 end

src/wing_geometry.jl

@@ -24,14 +24,13 @@ struct Section{T}
 end
 
 """
-    update_section!(wing::Wing, LE_point::PosVector, TE_point::PosVector)
+    update_pos!(wing::Wing, LE_point::PosVector, TE_point::PosVector)
 
 Update the section leading edge and trailing edge positions.
 """
-function update_section!(section::Section, LE_point::AbstractVector, TE_point::AbstractVector)
-    size(LE_point) != (3,) && throw(ArgumentError("LE_point has wrong size"))
-    size(TE_point) != (3,) && throw(ArgumentError("TE_point has wrong size"))
+function update_pos!(section::Section, LE_point::AbstractVector, TE_point::AbstractVector)
     section.LE_point .= LE_point
+    section.TE_point .= TE_point
     return nothing
 end
 

test/test_kite_geometry.jl

@@ -141,7 +141,7 @@ using Serialization
         wing = KiteWing(test_obj_path, test_dat_path)
 
         @test wing.n_panels == 54  # Default value
-        @test wing.spanwise_panel_distribution == LINEAR
+        @test wing.spanwise_panel_distribution == UNCHANGED
         @test wing.spanwise_direction ≈ [0.0, 1.0, 0.0]
         @test length(wing.sections) > 0  # Should have sections now
         @test wing.mass ≈ 1.0