Merge branch 'main' of https://github.com/Albatross-Kite-Transport/VortexStepMethod.jl

Author: 1-Bart-1

src/body_aerodynamics.jl

@@ -126,7 +126,8 @@ function init!(body_aero::BodyAerodynamics;
                 panel_props.x_airf[i],
                 panel_props.y_airf[i],
                 panel_props.z_airf[i],
-                beta
+                beta;
+                remove_nan=wing.remove_nan
             )
             idx += 1
         end

src/kite_geometry.jl

@@ -239,6 +239,56 @@ function calculate_inertia_tensor(vertices, faces, mass, com)
     return (mass / total_area) * I / 3
 end
 
+"""
+    interpolate_matrix_nans!(matrix::Matrix{Float64})
+
+Replace NaN values in a matrix by interpolating from nearest non-NaN neighbors.
+Uses an expanding search radius until valid neighbors are found.
+
+# Arguments
+- `matrix`: Matrix containing NaN values to be interpolated
+"""
+function interpolate_matrix_nans!(matrix::Matrix{Float64})
+    rows, cols = size(matrix)
+    nans_found = 0
+    while any(isnan, matrix)
+        for i in 1:rows, j in 1:cols
+            if isnan(matrix[i,j])
+                # Search in expanding radius until we find valid neighbors
+                radius = 1
+                values = Float64[]
+                weights = Float64[]
+                
+                while isempty(values) && radius < max(rows, cols)
+                    # Check all points at current Manhattan distance
+                    for di in -radius:radius, dj in -radius:radius
+                        if abs(di) + abs(dj) == radius  # Points exactly at distance 'radius'
+                            ni, nj = i + di, j + dj
+                            if 1 ≤ ni ≤ rows && 1 ≤ nj ≤ cols && !isnan(matrix[ni,nj])
+                                # Weight by inverse distance
+                                dist = sqrt(di^2 + dj^2)
+                                push!(values, matrix[ni,nj])
+                                push!(weights, 1/dist)
+                            end
+                        end
+                    end
+                    radius += 1
+                end
+                
+                if !isempty(values)
+                    # Calculate weighted average of found values
+                    matrix[i,j] = sum(values .* weights) / sum(weights)
+                    nans_found += 1
+                else
+                    throw(ArgumentError("Could not remove NaN"))
+                end
+            end
+        end
+    end
+    @info "Removed $nans_found NaNs from the matrix."
+    return matrix
+end
+
 
 """
     KiteWing
@@ -252,6 +302,7 @@ Represents a curved wing that inherits from Wing with additional geometric prope
   - `spanwise_direction::MVec3`: Wing span direction vector
   - `sections::Vector{Section}`: List of wing sections, see: [Section](@ref)
   -  refined_sections::Vector{Section}
+  - `remove_nan::Bool`: Wether to remove the NaNs from interpolations or not
 - Additional fields:
   - `center_of_mass::Vector{Float64}`: Center of mass coordinates
   - `circle_center_z::Vector{Float64}`: Center of circle coordinates
@@ -269,6 +320,7 @@ mutable struct KiteWing <: AbstractWing
     spanwise_direction::MVec3
     sections::Vector{Section}
     refined_sections::Vector{Section}
+    remove_nan::Bool
 
     # Additional fields for KiteWing
     non_deformed_sections::Vector{Section}
@@ -288,7 +340,7 @@ end
 """
     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=UNCHANGED, 
-             spanwise_direction=[0.0, 1.0, 0.0])
+             spanwise_direction=[0.0, 1.0, 0.0], remove_nan::Bool=true)
 
 Constructor for a [KiteWing](@ref) that allows to use an `.obj` and a `.dat` file as input.
 
@@ -306,10 +358,11 @@ Constructor for a [KiteWing](@ref) that allows to use an `.obj` and a `.dat` fil
 - `n_sections`=n_panels+1: Number of sections (there is a section on each side of each panel.)
 - `spanwise_panel_distribution`=UNCHANGED: see: [PanelDistribution](@ref)
 - `spanwise_direction`=[0.0, 1.0, 0.0]
+- `remove_nan::Bool`: Wether to remove the NaNs from interpolations or not
 """
 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=UNCHANGED, 
-                  spanwise_direction=[0.0, 1.0, 0.0])
+                  spanwise_direction=[0.0, 1.0, 0.0], remove_nan=true)
 
     !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"))
 
@@ -346,6 +399,11 @@ function KiteWing(obj_path, dat_path; alpha=0.0, crease_frac=0.75, wind_vel=10.,
 
     @info "Loading polars and kite info from $polar_path and $info_path"
     (alpha_range, beta_range, cl_matrix::Matrix, cd_matrix::Matrix, cm_matrix::Matrix) = deserialize(polar_path)
+    if remove_nan
+        interpolate_matrix_nans!(cl_matrix)
+        interpolate_matrix_nans!(cd_matrix)
+        interpolate_matrix_nans!(cm_matrix)
+    end
 
     (center_of_mass, inertia_tensor, circle_center_z, radius, gamma_tip, 
         le_interp, te_interp, area_interp) = deserialize(info_path)
@@ -359,7 +417,7 @@ function KiteWing(obj_path, dat_path; alpha=0.0, crease_frac=0.75, wind_vel=10.,
         push!(sections, Section(LE_point, TE_point, POLAR_MATRICES, aero_data))
     end
 
-    KiteWing(n_panels, spanwise_panel_distribution, spanwise_direction, sections, sections, sections,
+    KiteWing(n_panels, spanwise_panel_distribution, spanwise_direction, sections, sections, remove_nan, sections,
         mass, center_of_mass, circle_center_z, gamma_tip, inertia_tensor, radius,
         le_interp, te_interp, area_interp, zeros(n_panels), zeros(n_panels))
 end

src/panel.jl

@@ -1,6 +1,8 @@
 # static types for interpolations
 const I1 = Interpolations.FilledExtrapolation{Float64, 1, Interpolations.GriddedInterpolation{Float64, 1, Vector{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}}}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Float64}
-const I2 = Interpolations.FilledExtrapolation{Float64, 2, Interpolations.GriddedInterpolation{Float64, 2, Matrix{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}, Vector{Float64}}}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Float64}
+const I2 = Interpolations.Extrapolation{Float64, 1, Interpolations.GriddedInterpolation{Float64, 1, Vector{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}}}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Interpolations.Line{Nothing}}
+const I3 = Interpolations.FilledExtrapolation{Float64, 2, Interpolations.GriddedInterpolation{Float64, 2, Matrix{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}, Vector{Float64}}}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Float64}
+const I4 = Interpolations.Extrapolation{Float64, 2, Interpolations.GriddedInterpolation{Float64, 2, Matrix{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}, Vector{Float64}}}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Interpolations.Line{Nothing}}
 
 """
     @with_kw mutable struct Panel
@@ -50,9 +52,9 @@ Represents a panel in a vortex step method simulation. All points and vectors ar
     cl_coeffs::Vector{Float64} = zeros(Float64, 3)
     cd_coeffs::Vector{Float64} = zeros(Float64, 3)
     cm_coeffs::Vector{Float64} = zeros(Float64, 3)
-    cl_interp::Union{Nothing, I1, I2} = nothing
-    cd_interp::Union{Nothing, I1, I2} = nothing
-    cm_interp::Union{Nothing, I1, I2} = nothing
+    cl_interp::Union{Nothing, I1, I2, I3, I4} = nothing
+    cd_interp::Union{Nothing, I1, I2, I3, I4} = nothing
+    cm_interp::Union{Nothing, I1, I2, I3, I4} = nothing
     aero_center::MVec3 = zeros(MVec3)
     control_point::MVec3 = zeros(MVec3)
     bound_point_1::MVec3 = zeros(MVec3)
@@ -113,7 +115,8 @@ end
 function init_aero!(
     panel::Panel,
     section_1::Section,
-    section_2::Section,
+    section_2::Section;
+    remove_nan = true
 )
     # Validate aero model consistency
     panel.aero_model = section_1.aero_model
@@ -132,6 +135,12 @@ function init_aero!(
             throw(ArgumentError("Polar data must have same shape"))
         end
 
+        if remove_nan
+            extrapolation_bc = Line()
+        else
+            extrapolation_bc = NaN
+        end
+
         if panel.aero_model == POLAR_VECTORS
             !all(isapprox.(aero_1[1], aero_2[1])) && @error "Make sure you use the same alpha range for all your interpolations."
 
@@ -142,9 +151,9 @@ function init_aero!(
             )
             alphas = Vector{Float64}(aero_1[1])
 
-            panel.cl_interp = linear_interpolation(alphas, polar_data[1]; extrapolation_bc=NaN)
-            panel.cd_interp = linear_interpolation(alphas, polar_data[2]; extrapolation_bc=NaN)
-            panel.cm_interp = linear_interpolation(alphas, polar_data[3]; extrapolation_bc=NaN)
+            panel.cl_interp = linear_interpolation(alphas, polar_data[1]; extrapolation_bc)
+            panel.cd_interp = linear_interpolation(alphas, polar_data[2]; extrapolation_bc)
+            panel.cm_interp = linear_interpolation(alphas, polar_data[3]; extrapolation_bc)
 
         elseif panel.aero_model == POLAR_MATRICES
             !all(isapprox.(aero_1[1], aero_2[1])) && @error "Make sure you use the same alpha range for all your interpolations."
@@ -158,9 +167,9 @@ function init_aero!(
             alphas = Vector{Float64}(aero_1[1])
             betas = Vector{Float64}(aero_1[2])
 
-            panel.cl_interp = linear_interpolation((alphas, betas), polar_data[1]; extrapolation_bc=NaN)
-            panel.cd_interp = linear_interpolation((alphas, betas), polar_data[2]; extrapolation_bc=NaN)
-            panel.cm_interp = linear_interpolation((alphas, betas), polar_data[3]; extrapolation_bc=NaN)
+            panel.cl_interp = linear_interpolation((alphas, betas), polar_data[1]; extrapolation_bc)
+            panel.cd_interp = linear_interpolation((alphas, betas), polar_data[2]; extrapolation_bc)
+            panel.cm_interp = linear_interpolation((alphas, betas), polar_data[3]; extrapolation_bc)
         else
             throw(ArgumentError("Polar data in wrong format: $aero_1"))
         end
@@ -182,11 +191,12 @@ function init!(
     y_airf::PosVector,
     z_airf::PosVector,
     beta;
-    init_aero = true
+    init_aero = true,
+    remove_nan = true
 )
     init_pos!(panel, section_1, section_2, aero_center, control_point, bound_point_1, bound_point_2,
         x_airf, y_airf, z_airf, beta)
-    init_aero && init_aero!(panel, section_1, section_2)
+    init_aero && init_aero!(panel, section_1, section_2; remove_nan)
     return nothing
 end
 

src/wing_geometry.jl

@@ -69,6 +69,7 @@ Represents a wing composed of multiple sections with aerodynamic properties.
 - `spanwise_direction::MVec3`: Wing span direction vector
 - `sections::Vector{Section}`: Vector of wing sections, see: [Section](@ref)
 - `refined_sections::Vector{Section}`: Vector of refined wing sections, see: [Section](@ref)
+- `remove_nan::Bool`: Wether to remove the NaNs from interpolations or not
 
 """
 mutable struct Wing <: AbstractWing
@@ -77,12 +78,14 @@ mutable struct Wing <: AbstractWing
     spanwise_direction::MVec3
     sections::Vector{Section}
     refined_sections::Vector{Section}
+    remove_nan::Bool
 end
 
 """
     Wing(n_panels::Int;
          spanwise_panel_distribution::PanelDistribution=LINEAR,
-         spanwise_direction::PosVector=MVec3([0.0, 1.0, 0.0]))
+         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.
@@ -91,11 +94,13 @@ and refined sections as empty arrays.
 - `n_panels::Int64`: Number of panels in aerodynamic mesh
 - `spanwise_panel_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
 """
 function Wing(n_panels::Int;
     spanwise_panel_distribution::PanelDistribution=LINEAR,
-    spanwise_direction::PosVector=MVec3([0.0, 1.0, 0.0]))
-    Wing(n_panels, spanwise_panel_distribution, spanwise_direction, Section[], Section[])
+    spanwise_direction::PosVector=MVec3([0.0, 1.0, 0.0]),
+    remove_nan=true)
+    Wing(n_panels, spanwise_panel_distribution, spanwise_direction, Section[], Section[], remove_nan)
 end
 
 function init!(wing::AbstractWing; aero_center_location::Float64=0.25, control_point_location::Float64=0.75)
@@ -111,6 +116,36 @@ function init!(wing::AbstractWing; aero_center_location::Float64=0.25, control_p
     return panel_props
 end
 
+
+"""
+    remove_vector_nans(aero_data)
+
+Remove the indices from aero_data where a NaN is found.
+"""
+function remove_vector_nans(aero_data)
+    alpha_range, cl_vector, cd_vector, cm_vector = aero_data
+    alpha_range = collect(alpha_range)
+    nan_indices = Set{Int}()
+    for vec in (cl_vector, cd_vector, cm_vector)
+        union!(nan_indices, findall(isnan, vec))
+    end
+    if isempty(nan_indices)
+        return aero_data
+    end
+    # Convert to sorted array for consistent removal
+    nan_indices = sort(collect(nan_indices))
+    # Create mask for valid indices
+    valid_mask = trues(length(alpha_range))
+    valid_mask[nan_indices] .= false
+    # Remove NaN values from all vectors
+    clean_alpha = alpha_range[valid_mask]
+    clean_cl = cl_vector[valid_mask]
+    clean_cd = cd_vector[valid_mask]
+    clean_cm = cm_vector[valid_mask]
+    @info "Removed $(length(nan_indices)) indices containing NaNs from aero_data."
+    return (clean_alpha, clean_cl, clean_cd, clean_cm)
+end
+
 """
     add_section!(wing::Wing, LE_point::PosVector, TE_point::PosVector, 
                  aero_model, aero_data::AeroData=nothing)
@@ -126,6 +161,11 @@ Add a new section to the wing.
 """
 function add_section!(wing::Wing, LE_point::Vector{Float64}, 
                      TE_point::Vector{Float64}, aero_model::AeroModel, aero_data::AeroData=nothing)
+    if aero_model === POLAR_VECTORS && wing.remove_nan
+        aero_data = remove_vector_nans(aero_data)
+    elseif aero_model === POLAR_MATRICES && wing.remove_nan
+        interpolate_matrix_nans!.(aero_data[3:5])
+    end
     push!(wing.sections, Section(LE_point, TE_point, aero_model, aero_data))
     return nothing
 end

test/test_kite_geometry.jl

@@ -99,6 +99,9 @@ using Serialization
                 cm_matrix[i,j] = -0.1*sin(deg2rad(alphas[i]))
             end
         end
+        cl_matrix[end] = NaN
+        cd_matrix[end] = NaN
+        cm_matrix[end] = NaN
 
         # Serialize polar data
         polar_path = test_dat_path[1:end-4] * "_polar.bin"
@@ -140,7 +143,7 @@ using Serialization
     end
 
     @testset "KiteWing Construction" begin
-        wing = KiteWing(test_obj_path, test_dat_path)
+        wing = KiteWing(test_obj_path, test_dat_path; remove_nan=true)
 
         @test wing.n_panels == 54  # Default value
         @test wing.spanwise_panel_distribution == UNCHANGED
@@ -150,6 +153,14 @@ using Serialization
         @test length(wing.center_of_mass) == 3
         @test wing.radius ≈ r rtol=1e-2
         @test wing.gamma_tip ≈ π/4 rtol=1e-2
+        @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])
+
+        wing = KiteWing(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
 
     rm(test_obj_path)

test/test_wing_geometry.jl

@@ -39,6 +39,59 @@ end
         @test section.aero_model === INVISCID
     end
 
+    @testset "Vector polar data" begin
+        wing = Wing(2; remove_nan=true)
+        
+        # Create vector data with NaNs
+        alpha = range(-10, 10, 21)
+        cl = cos.(alpha)
+        cd = sin.(alpha)
+        cm = zeros(21)
+        
+        # Insert NaNs
+        cl[5] = NaN
+        cd[10] = NaN
+        cm[15] = NaN
+        
+        aero_data = (collect(alpha), cl, cd, cm)
+        add_section!(wing, [0.0, 0.0, 0.0], [1.0, 0.0, 0.0], POLAR_VECTORS, aero_data)
+        
+        # Check if NaNs were removed consistently
+        cleaned_data = wing.sections[1].aero_data
+        @test length(cleaned_data[1]) == 18  # 21 - 3 NaN positions
+        @test !any(isnan, cleaned_data[2])   # cl
+        @test !any(isnan, cleaned_data[3])   # cd
+        @test !any(isnan, cleaned_data[4])   # cm
+        @test all(diff(cleaned_data[1]) .> 0) # alpha still monotonic
+    end
+
+    @testset "Matrix polar data" begin
+        wing = Wing(2)
+        
+        # Create matrix data with NaNs
+        alpha = collect(range(-10, 10, 21))
+        beta = collect(range(-5, 5, 11))
+        cl = [cos(a) * cos(b) for a in alpha, b in beta]
+        cd = [sin(a) * sin(b) for a in alpha, b in beta]
+        cm = zeros(21, 11)
+        
+        # Insert NaNs at various positions
+        cl[5,3] = NaN
+        cd[10,5] = NaN
+        cm[15,7] = NaN
+        
+        aero_data = (alpha, beta, cl, cd, cm)
+        add_section!(wing, [0.0, 0.0, 0.0], [1.0, 0.0, 0.0], POLAR_MATRICES, aero_data)
+        
+        # Check if NaNs were removed consistently
+        cleaned_data = wing.sections[1].aero_data
+        @test !any(isnan, cleaned_data[3])   # cl
+        @test !any(isnan, cleaned_data[4])   # cd
+        @test !any(isnan, cleaned_data[5])   # cm
+        @test all(diff(cleaned_data[1]) .> 0) # alpha still monotonic
+        @test all(diff(cleaned_data[2]) .> 0) # beta still monotonic
+    end
+
     @testset "Robustness left to right" begin
         example_wing = Wing(10)
         # Test correct order