From 75746dfa1b5a786a556f9ee757cca89a294ddb0e Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Wed, 16 Sep 2026 12:25:17 +0200 Subject: [PATCH 1/2] Delete the inertia functions ObjAdapter duplicated from SymbolicAWEModels center_to_com!, calculate_inertia_tensor and calc_inertia_y_rotation are removed from ObjAdapter with their export, tests and docs entry. read_faces stays. Also drops the test scaffolding only they used: the serialized _info.bin, the .obj the alignment test read back, and the Serialization test dependency. Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 8 ++ docs/src/private_functions.md | 2 - src/obj_adapter/ObjAdapter.jl | 2 +- src/obj_adapter/obj_geometry.jl | 144 ------------------------ test/Project.toml | 2 - test/obj_adapter/test_obj_adapter.jl | 5 - test/ram_geometry/test_kite_geometry.jl | 61 +--------- 7 files changed, 11 insertions(+), 213 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 532b0816..d5ea457b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,13 @@ # Changelog +## Unreleased + +### Changed + +- BREAKING: `ObjAdapter` no longer exports `center_to_com!`, `calculate_inertia_tensor` + or `calc_inertia_y_rotation`; they are gone. Mesh mass properties are computed by + SymbolicAWEModels, which reads the mesh with `read_faces`. + ## VortexStepMethod v5.1.1 2026-09-12 ### Fixed diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index d7807c42..93dd8d3e 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -244,8 +244,6 @@ densify_contour create_interpolations find_circle_center_and_radius march_edges -calculate_inertia_tensor -center_to_com! airfoils_from_yaml write_geometry_yaml resolve_aero_geometry diff --git a/src/obj_adapter/ObjAdapter.jl b/src/obj_adapter/ObjAdapter.jl index 492ca228..b60a6821 100644 --- a/src/obj_adapter/ObjAdapter.jl +++ b/src/obj_adapter/ObjAdapter.jl @@ -45,7 +45,7 @@ function plot_slices_3d end export obj_to_yaml, resolve_aero_geometry export write_geometry_yaml, write_yaml, airfoils_from_yaml export perpendicular_sections -export read_faces, center_to_com!, calculate_inertia_tensor, calc_inertia_y_rotation +export read_faces export plot_airfoil_fit, plot_airfoils, plot_slices_3d end diff --git a/src/obj_adapter/obj_geometry.jl b/src/obj_adapter/obj_geometry.jl index 0a5ac73a..09ee3675 100644 --- a/src/obj_adapter/obj_geometry.jl +++ b/src/obj_adapter/obj_geometry.jl @@ -226,147 +226,3 @@ function create_interpolations(vertices, circle_center_z, radius, gamma_tip, R=I return (le_interp, te_interp, area_interp) end - -""" - center_to_com!(vertices, faces) - -Calculate center of mass of a mesh and translate vertices so that COM is at origin. - -# Arguments -- `vertices`: Vector of 3D point coordinates -- `faces`: Vector of vertex indices for each face (can be triangular or non-triangular) - -# Returns -- Vector representing the original center of mass before translation - -# Notes -- Non-triangular faces are automatically triangulated into triangles -- Assumes uniform surface density -""" -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" - com[2] = 0.0 - for v in vertices - v .+= com - end - return com -end - -""" - calculate_inertia_tensor(vertices, faces, mass, com) - -Calculate the inertia tensor for a triangulated surface mesh, assuming a thin shell with uniform -surface density. - -# Arguments -- `vertices`: Vector of 3D point coordinates representing mesh vertices -- `faces`: Vector of triangle indices, each defining a face of the mesh -- `mass`: Total mass of the shell in kg -- `com`: Center of mass coordinates [x,y,z] - -# Method -Uses the thin shell approximation where: -1. Mass is distributed uniformly over the surface area -2. Each triangle contributes to the inertia based on its area and position -3. For each triangle vertex p, contribution to diagonal terms is: area * (sum(p²) - p_i²) -4. For off-diagonal terms: area * (-`p_i` * `p_j`) -5. Final tensor is scaled by mass/(3*total_area) to get correct units - -# Returns -- 3×3 matrix representing the inertia tensor in kg⋅m² -""" -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 - # Diagonal terms - I[i,i] += area * (sum(p.^2) - p[i]^2) - else - # Off-diagonal terms - I[i,j] -= area * (p[i] * p[j]) - end - end - end - end - end - - # Scale by mass/total_area to get actual inertia tensor - return (mass / total_area) * I / 3 -end - -function calc_inertia_y_rotation(I_b_tensor) - # Function for nonlinear solver - off-diagonal element should be zero - function eq!(F, theta, _) - # Rotation matrix around y-axis - R_y = [ - cos(theta[1]) 0 sin(theta[1]); - 0 1 0; - -sin(theta[1]) 0 cos(theta[1]) - ] - # 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] - 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; - -sin(theta_opt) 0 cos(theta_opt) - ] - # Calculate diagonalized inertia tensor - I_diag = R_b_p * I_b_tensor * R_b_p' - @assert isapprox(I_diag[1,3], 0.0, atol=1e-5) - return I_diag, R_b_p -end - - diff --git a/test/Project.toml b/test/Project.toml index 5e1ffca3..7e0423ef 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -14,7 +14,6 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" MakieControlPlots = "6d616b69-6563-4f6e-8472-6f6c706c6f74" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" @@ -36,7 +35,6 @@ LinearAlgebra = "1" Logging = "1" MakieControlPlots = "0.1.5" Random = "1.10.0" -Serialization = "1" StaticArrays = "1" Statistics = "1" Test = "1" diff --git a/test/obj_adapter/test_obj_adapter.jl b/test/obj_adapter/test_obj_adapter.jl index e2f0ae91..50732486 100644 --- a/test/obj_adapter/test_obj_adapter.jl +++ b/test/obj_adapter/test_obj_adapter.jl @@ -113,11 +113,6 @@ obj_path = normpath(joinpath(@__DIR__, "..", "..", tables) == "" end - @testset "center_to_com! rejects non-triangular faces" begin - verts = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 1.0, 0.0]] - @test_throws ArgumentError center_to_com!(verts, [[1, 2, 3, 4]]; prn=false) - end - @testset "write_yaml emits nested and scalar values" begin dir = mktempdir() nested = joinpath(dir, "nested.yaml") diff --git a/test/ram_geometry/test_kite_geometry.jl b/test/ram_geometry/test_kite_geometry.jl index 4c9fdc8d..09fbdb1c 100644 --- a/test/ram_geometry/test_kite_geometry.jl +++ b/test/ram_geometry/test_kite_geometry.jl @@ -4,10 +4,9 @@ using VortexStepMethod using VortexStepMethod: read_aero_matrix using VortexStepMethod.AirfoilAero: write_aero_matrix using VortexStepMethod.ObjAdapter: create_interpolations, find_circle_center_and_radius, - calculate_inertia_tensor, center_to_com!, read_faces, calc_inertia_y_rotation + read_faces using LinearAlgebra using Interpolations -using Serialization @testset "Kite Geometry Tests" begin work_dir = mktempdir() @@ -34,30 +33,6 @@ using Serialization @test faces[1] == [1, 2, 3] end - @testset "Center of Mass Calculation" begin - vertices = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]] - faces = [[1, 2, 3]] - - com = center_to_com!(vertices, faces) - expected_com = [-1/3, 0.0, -1/3] - - @test isapprox(com, expected_com, rtol=1e-5) - end - - @testset "Inertia Tensor Calculation" begin - vertices = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]] - faces = [[1, 2, 3]] - mass = 1.0 - com = [1/3, 1/3, 0.0] - - I = calculate_inertia_tensor(vertices, faces, mass, com) - - # Test properties of inertia tensor - @test size(I) == (3,3) - @test isapprox(I, I', rtol=1e-10) # Symmetric - @test all(diag(I) .≥ 0) # Non-negative diagonal - end - @testset "Circle Fitting" begin # Create simple curved wing vertices r = 5.0 @@ -120,17 +95,6 @@ using Serialization write_aero_matrix(cd_polar_path, cd_matrix, deg2rad.(alphas), deg2rad.(d_trailing_edge_angles), "C_d") write_aero_matrix(cm_polar_path, cm_matrix, deg2rad.(alphas), deg2rad.(d_trailing_edge_angles), "C_m") - # Create and serialize obj file - faces = [[i, i+1, i+2] for i in 1:3:length(vertices)-2] - open(test_obj_path, "w") do io - for v in vertices - println(io, "v $(v[1]) $(v[2]) $(v[3])") - end - for f in faces - println(io, "f $(f[1]) $(f[2]) $(f[3])") - end - end - # Test reading back the matrices cl_read, alphas_read, deltas_read = read_aero_matrix(cl_polar_path) # write_aero_matrix stores coefficients rounded to 4 decimals @@ -139,34 +103,13 @@ using Serialization @test alphas_read ≈ deg2rad.(alphas) @test deltas_read ≈ deg2rad.(d_trailing_edge_angles) - # Create info file - info_path = test_obj_path[1:end-4] * "_info.bin" le_interp, te_interp, area_interp = create_interpolations(vertices, z_center, r, π/4, I(3)) - center_of_mass = center_to_com!(vertices, faces) - inertia_tensor = calculate_inertia_tensor(vertices, faces, 1.0, zeros(3)) - - serialize(info_path, (inertia_tensor, center_of_mass, I(3), r, π/4, - le_interp, te_interp, area_interp)) - + # Test interpolation at middle point @test isapprox([le_interp[i](0.0) for i in 1:3], [0.0, 0.0, r+z_center], atol=0.03) @test isapprox([te_interp[i](0.0) for i in 1:3], [1.0, 0.0, r+z_center], atol=0.03) end - @testset "Alignment to principal frame" begin - vertices, faces = read_faces(test_obj_path) - center_of_mass = center_to_com!(vertices, faces) - inertia_tensor_b = calculate_inertia_tensor(vertices, faces, 1.0, zeros(3)) - inertia_tensor_p, R_b_p = calc_inertia_y_rotation(inertia_tensor_b) - for v in vertices - v .= R_b_p * v - end - inertia_tensor_b2 = calculate_inertia_tensor(vertices, faces, 1.0, zeros(3)) - inertia_tensor_p2, R_b_p2 = calc_inertia_y_rotation(inertia_tensor_b2) - @test inertia_tensor_p ≈ inertia_tensor_p2 - @test R_b_p2 ≈ I(3) - end - @testset "Converted-wing construction and deformation" begin # TODO: redesign. These previously tested ObjWing internals (radius, # gamma_tip, UNCHANGED distribution, obj deform\!) that were dropped when From e344c8e08e8e9bcd818eac9594b37a6e4d824794 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Thu, 17 Sep 2026 00:01:19 +0200 Subject: [PATCH 2/2] Deprecate the ObjAdapter inertia functions instead of deleting them Per review, center_to_com!, calculate_inertia_tensor and calc_inertia_y_rotation come back with their export, tests and docs entry, and each now emits a deprecation warning naming SymbolicAWEModels. Their tests assert the warning with @test_deprecated. The changelog entry is a deprecation rather than BREAKING; removal waits for the next breaking release. The dead _info.bin serialization in test_kite_geometry.jl and the Serialization test dependency stay removed. Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 6 +- docs/src/private_functions.md | 3 + src/obj_adapter/ObjAdapter.jl | 2 +- src/obj_adapter/obj_geometry.jl | 162 ++++++++++++++++++++++++ test/obj_adapter/test_obj_adapter.jl | 5 + test/ram_geometry/test_kite_geometry.jl | 51 +++++++- 6 files changed, 224 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c6b93860..3f30a36a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,9 +4,9 @@ ### Changed -- BREAKING: `ObjAdapter` no longer exports `center_to_com!`, `calculate_inertia_tensor` - or `calc_inertia_y_rotation`; they are gone. Mesh mass properties are computed by - SymbolicAWEModels, which reads the mesh with `read_faces`. +- `ObjAdapter.center_to_com!`, `calculate_inertia_tensor` and `calc_inertia_y_rotation` + are deprecated and will be removed in the next breaking release. Mesh mass properties + are computed by SymbolicAWEModels, which reads the mesh with `read_faces`. - Requires Julia 1.12 or 1.13; 1.10 and 1.11 keep resolving v5.1.1. ## VortexStepMethod v5.1.1 2026-09-12 diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 93dd8d3e..415c743f 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -244,6 +244,9 @@ densify_contour create_interpolations find_circle_center_and_radius march_edges +calculate_inertia_tensor +center_to_com! +calc_inertia_y_rotation airfoils_from_yaml write_geometry_yaml resolve_aero_geometry diff --git a/src/obj_adapter/ObjAdapter.jl b/src/obj_adapter/ObjAdapter.jl index b60a6821..492ca228 100644 --- a/src/obj_adapter/ObjAdapter.jl +++ b/src/obj_adapter/ObjAdapter.jl @@ -45,7 +45,7 @@ function plot_slices_3d end export obj_to_yaml, resolve_aero_geometry export write_geometry_yaml, write_yaml, airfoils_from_yaml export perpendicular_sections -export read_faces +export read_faces, center_to_com!, calculate_inertia_tensor, calc_inertia_y_rotation export plot_airfoil_fit, plot_airfoils, plot_slices_3d end diff --git a/src/obj_adapter/obj_geometry.jl b/src/obj_adapter/obj_geometry.jl index 09ee3675..e427d6f2 100644 --- a/src/obj_adapter/obj_geometry.jl +++ b/src/obj_adapter/obj_geometry.jl @@ -226,3 +226,165 @@ function create_interpolations(vertices, circle_center_z, radius, gamma_tip, R=I return (le_interp, te_interp, area_interp) end + +function depwarn_inertia(name::Symbol) + Base.depwarn("`$name` is deprecated and will be removed in the next breaking release; " * + "SymbolicAWEModels computes mesh mass properties.", name) +end + +""" + center_to_com!(vertices, faces) + +Calculate center of mass of a mesh and translate vertices so that COM is at origin. + +# Arguments +- `vertices`: Vector of 3D point coordinates +- `faces`: Vector of vertex indices for each face (can be triangular or non-triangular) + +# Returns +- Vector representing the original center of mass before translation + +# Notes +- Non-triangular faces are automatically triangulated into triangles +- Assumes uniform surface density + +Deprecated; removed in the next breaking release. +""" +function center_to_com!(vertices, faces; prn=true) + depwarn_inertia(:center_to_com!) + 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" + com[2] = 0.0 + for v in vertices + v .+= com + end + return com +end + +""" + calculate_inertia_tensor(vertices, faces, mass, com) + +Calculate the inertia tensor for a triangulated surface mesh, assuming a thin shell with uniform +surface density. + +# Arguments +- `vertices`: Vector of 3D point coordinates representing mesh vertices +- `faces`: Vector of triangle indices, each defining a face of the mesh +- `mass`: Total mass of the shell in kg +- `com`: Center of mass coordinates [x,y,z] + +# Method +Uses the thin shell approximation where: +1. Mass is distributed uniformly over the surface area +2. Each triangle contributes to the inertia based on its area and position +3. For each triangle vertex p, contribution to diagonal terms is: area * (sum(p²) - p_i²) +4. For off-diagonal terms: area * (-`p_i` * `p_j`) +5. Final tensor is scaled by mass/(3*total_area) to get correct units + +# Returns +- 3×3 matrix representing the inertia tensor in kg⋅m² + +Deprecated; removed in the next breaking release. +""" +function calculate_inertia_tensor(vertices, faces, mass, com) + depwarn_inertia(:calculate_inertia_tensor) + # 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 + # Diagonal terms + I[i,i] += area * (sum(p.^2) - p[i]^2) + else + # Off-diagonal terms + I[i,j] -= area * (p[i] * p[j]) + end + end + end + end + end + + # Scale by mass/total_area to get actual inertia tensor + return (mass / total_area) * I / 3 +end + +""" + calc_inertia_y_rotation(I_b_tensor) + +Rotate the inertia tensor `I_b_tensor` about the y-axis until its xz product of inertia +vanishes. Returns the rotated tensor and the rotation matrix `R_b_p`. + +Deprecated; removed in the next breaking release. +""" +function calc_inertia_y_rotation(I_b_tensor) + depwarn_inertia(:calc_inertia_y_rotation) + # Function for nonlinear solver - off-diagonal element should be zero + function eq!(F, theta, _) + # Rotation matrix around y-axis + R_y = [ + cos(theta[1]) 0 sin(theta[1]); + 0 1 0; + -sin(theta[1]) 0 cos(theta[1]) + ] + # 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] + 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; + -sin(theta_opt) 0 cos(theta_opt) + ] + # Calculate diagonalized inertia tensor + I_diag = R_b_p * I_b_tensor * R_b_p' + @assert isapprox(I_diag[1,3], 0.0, atol=1e-5) + return I_diag, R_b_p +end diff --git a/test/obj_adapter/test_obj_adapter.jl b/test/obj_adapter/test_obj_adapter.jl index 50732486..e2f0ae91 100644 --- a/test/obj_adapter/test_obj_adapter.jl +++ b/test/obj_adapter/test_obj_adapter.jl @@ -113,6 +113,11 @@ obj_path = normpath(joinpath(@__DIR__, "..", "..", tables) == "" end + @testset "center_to_com! rejects non-triangular faces" begin + verts = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 1.0, 0.0]] + @test_throws ArgumentError center_to_com!(verts, [[1, 2, 3, 4]]; prn=false) + end + @testset "write_yaml emits nested and scalar values" begin dir = mktempdir() nested = joinpath(dir, "nested.yaml") diff --git a/test/ram_geometry/test_kite_geometry.jl b/test/ram_geometry/test_kite_geometry.jl index 09fbdb1c..8fe5d479 100644 --- a/test/ram_geometry/test_kite_geometry.jl +++ b/test/ram_geometry/test_kite_geometry.jl @@ -4,7 +4,7 @@ using VortexStepMethod using VortexStepMethod: read_aero_matrix using VortexStepMethod.AirfoilAero: write_aero_matrix using VortexStepMethod.ObjAdapter: create_interpolations, find_circle_center_and_radius, - read_faces + calculate_inertia_tensor, center_to_com!, read_faces, calc_inertia_y_rotation using LinearAlgebra using Interpolations @@ -33,6 +33,30 @@ using Interpolations @test faces[1] == [1, 2, 3] end + @testset "Center of Mass Calculation" begin + vertices = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]] + faces = [[1, 2, 3]] + + com = @test_deprecated center_to_com!(vertices, faces) + expected_com = [-1/3, 0.0, -1/3] + + @test isapprox(com, expected_com, rtol=1e-5) + end + + @testset "Inertia Tensor Calculation" begin + vertices = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]] + faces = [[1, 2, 3]] + mass = 1.0 + com = [1/3, 1/3, 0.0] + + I = @test_deprecated calculate_inertia_tensor(vertices, faces, mass, com) + + # Test properties of inertia tensor + @test size(I) == (3,3) + @test isapprox(I, I', rtol=1e-10) # Symmetric + @test all(diag(I) .≥ 0) # Non-negative diagonal + end + @testset "Circle Fitting" begin # Create simple curved wing vertices r = 5.0 @@ -95,6 +119,17 @@ using Interpolations write_aero_matrix(cd_polar_path, cd_matrix, deg2rad.(alphas), deg2rad.(d_trailing_edge_angles), "C_d") write_aero_matrix(cm_polar_path, cm_matrix, deg2rad.(alphas), deg2rad.(d_trailing_edge_angles), "C_m") + # Create obj file + faces = [[i, i+1, i+2] for i in 1:3:length(vertices)-2] + open(test_obj_path, "w") do io + for v in vertices + println(io, "v $(v[1]) $(v[2]) $(v[3])") + end + for f in faces + println(io, "f $(f[1]) $(f[2]) $(f[3])") + end + end + # Test reading back the matrices cl_read, alphas_read, deltas_read = read_aero_matrix(cl_polar_path) # write_aero_matrix stores coefficients rounded to 4 decimals @@ -110,6 +145,20 @@ using Interpolations @test isapprox([te_interp[i](0.0) for i in 1:3], [1.0, 0.0, r+z_center], atol=0.03) end + @testset "Alignment to principal frame" begin + vertices, faces = read_faces(test_obj_path) + center_of_mass = center_to_com!(vertices, faces) + inertia_tensor_b = calculate_inertia_tensor(vertices, faces, 1.0, zeros(3)) + inertia_tensor_p, R_b_p = @test_deprecated calc_inertia_y_rotation(inertia_tensor_b) + for v in vertices + v .= R_b_p * v + end + inertia_tensor_b2 = calculate_inertia_tensor(vertices, faces, 1.0, zeros(3)) + inertia_tensor_p2, R_b_p2 = calc_inertia_y_rotation(inertia_tensor_b2) + @test inertia_tensor_p ≈ inertia_tensor_p2 + @test R_b_p2 ≈ I(3) + end + @testset "Converted-wing construction and deformation" begin # TODO: redesign. These previously tested ObjWing internals (radius, # gamma_tip, UNCHANGED distribution, obj deform\!) that were dropped when