diff --git a/CHANGELOG.md b/CHANGELOG.md index b69c9a1b..73cdbd1a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,9 @@ ### Changed +- BREAKING: `ObjAdapter.center_to_com!`, `calculate_inertia_tensor` and + `calc_inertia_y_rotation` are removed. 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. - The Makie `plot!` methods for a `Panel` or a `BodyAerodynamics` return a `Vector{Makie.AbstractPlot}` instead of a `Vector{Any}`; for a `BodyAerodynamics` diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index ad67330d..8580b5fd 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -250,8 +250,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 8697ad85..74fced87 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 64d0d1ec..f6f7c676 100644 --- a/test/ram_geometry/test_kite_geometry.jl +++ b/test/ram_geometry/test_kite_geometry.jl @@ -4,10 +4,8 @@ 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() @@ -36,30 +34,6 @@ using Serialization @test faces isa Vector{Vector{Int64}} 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 +94,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 +102,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