Skip to content
Open
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
@@ -1,5 +1,13 @@
# Changelog

## Unreleased

### Changed

- BREAKING: `ObjAdapter` no longer exports `center_to_com!`, `calculate_inertia_tensor`

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Keep the functions, mark them deprecated, and rm them in the next breaking release, because otherwise we have too many breaking changes.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done in e344c8e: the three functions are back with their export, tests and docs, each emits a deprecation warning, and the changelog entry is a deprecation instead of BREAKING.

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
Expand Down
2 changes: 0 additions & 2 deletions docs/src/private_functions.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/obj_adapter/ObjAdapter.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
144 changes: 0 additions & 144 deletions src/obj_adapter/obj_geometry.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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


2 changes: 0 additions & 2 deletions test/Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -36,7 +35,6 @@ LinearAlgebra = "1"
Logging = "1"
MakieControlPlots = "0.1.5"
Random = "1.10.0"
Serialization = "1"
StaticArrays = "1"
Statistics = "1"
Test = "1"
Expand Down
5 changes: 0 additions & 5 deletions test/obj_adapter/test_obj_adapter.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
61 changes: 2 additions & 59 deletions test/ram_geometry/test_kite_geometry.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
Loading