Fixed LCS affine transformation

src/AeroMDAO.jl

@@ -78,9 +78,9 @@ export DoubletLine3D, velocity
 #==========================================================================================#
 
 include("Geometry/PanelGeometry/PanelGeometry.jl")
-import .PanelGeometry: AbstractPanel, AbstractPanel2D, Panel2D, WakePanel2D, AbstractPanel3D, Panel3D, panel_length, transform_panel, transform_panel_points, panel_angle, tangent_vector, normal_vector, panel_location, panel_area, panel_coordinates, transform, midpoint, panel_points, wake_panel, wake_panels, reverse_panel, panel_velocity, panel_scalar, trailing_edge_panel, get_surface_values, panel_vector, distance, average_chord, average_width, wetted_area, make_panels, local_coordinate_system, trailing_edge_info, panel_coordinates
+import .PanelGeometry: AbstractPanel, AbstractPanel2D, Panel2D, WakePanel2D, AbstractPanel3D, Panel3D, panel_length, transform_panel, transform_panel_points, panel_angle, tangent_vector, normal_vector, panel_location, panel_area, panel_coordinates, transform, midpoint, panel_points, wake_panel, wake_panels, reverse_panel, panel_velocity, panel_scalar, trailing_edge_panel, get_surface_values, panel_vector, distance, average_chord, average_width, wetted_area, make_panels, local_coordinate_system, get_transformation, trailing_edge_info, panel_coordinates
 
-export AbstractPanel, AbstractPanel2D, Panel2D, WakePanel2D, AbstractPanel3D, Panel3D, transform, normal_vector, midpoint, panel_location, tangent_vector, panel_points, distance, wake_panel, wake_panels, panel_area, reverse_panel, panel_length, transform_panel, panel_angle, panel_vector, panel_velocity, panel_scalar, trailing_edge_panel, get_surface_values, average_chord, average_width, wetted_area, make_panels, local_coordinate_system, trailing_edge_info, panel_coordinates
+export AbstractPanel, AbstractPanel2D, Panel2D, WakePanel2D, AbstractPanel3D, Panel3D, transform, normal_vector, midpoint, panel_location, tangent_vector, panel_points, distance, wake_panel, wake_panels, panel_area, reverse_panel, panel_length, transform_panel, panel_angle, panel_vector, panel_velocity, panel_scalar, trailing_edge_panel, get_surface_values, average_chord, average_width, wetted_area, make_panels, local_coordinate_system, get_transformation, trailing_edge_info, panel_coordinates
 
 ## Aircraft geometry
 #==========================================================================================#
@@ -191,6 +191,6 @@ export AerostructWing, make_beam_mesh, transform_stiffy, permute_stiffy, build_b
 
 include("Tools/plot_tools.jl")
 
-export plot_panels, plot_streamlines, plot_planform, plot_surface
+export plot_panel, plot_panels, plot_streamlines, plot_planform, plot_surface
 
 end

src/Geometry/PanelGeometry/3d_panels.jl

@@ -45,7 +45,9 @@ collocation_point(panel :: AbstractPanel3D, a = 0.5) = (p1(panel) + p2(panel) +
 +(p :: AbstractPanel3D, v) = Panel3D(p.p1 + v, p.p2 + v, p.p3 + v, p.p4 + v)
 -(p :: AbstractPanel3D, v) = Panel3D(p.p1 - v, p.p2 - v, p.p3 - v, p.p4 - v)
 ×(p :: AbstractPanel3D, v) = Panel3D(p.p1 × v, p.p2 × v, p.p3 × v, p.p4 × v)
+
 (T :: AffineMap)(p :: AbstractPanel3D) = typeof(p)(T(p.p1), T(p.p2), T(p.p3), T(p.p4))
+(T :: LinearMap)(p :: AbstractPanel3D) = typeof(p)(T(p.p1), T(p.p2), T(p.p3), T(p.p4))
 
 Base.length(:: Panel3D) = 1
 
@@ -81,20 +83,17 @@ function transform(panel :: Panel3D, rotation, translation)
     Panel3D(T(panel.p1), T(panel.p2), T(panel.p3), T(panel.p4))
 end
 
-(trans :: AffineMap)(p :: Panel3D) = Panel3D(trans(p.p1), trans(p.p2), trans(p.p3), trans(p.p4))
-
 """
     midpoint(panel :: AbstractPanel3D)
 
 Compute the midpoint of an `AbstractPanel3D`.
 """
 midpoint(panel :: AbstractPanel3D) = (p1(panel) + p2(panel) + p3(panel) + p4(panel)) / 4
 
-
 """
     normal_vector(panel :: Panel3D)
 
-Compute the unit normal vector of an `AbstractPanel3D` normalisation.
+Compute the normal vector of an `AbstractPanel3D`.
 """
 normal_vector(panel :: AbstractPanel3D) = let p31 = panel.p3 - panel.p1, p42 = panel.p4 - panel.p2; cross(p31, p42) end
 
@@ -112,7 +111,7 @@ transform_normal(panel :: Panel3D, h_l, g_l) = g_l * cross(h_l, normal_vector(pa
 
 Compute the area of a planar quadrilateral 3D panel.
 """
-panel_area(panel :: Panel3D) = 1/2 * norm(normal_vector(panel)) # (norm ∘ cross)(average_chord(panel), average_width(panel))
+panel_area(panel :: Panel3D) = 1/2 * norm(normal_vector(panel))
 
 """
     wetted_area(panels :: Array{Panel3D})
@@ -129,17 +128,17 @@ Reflect a Panel3D with respect to the ``x``-``z`` plane of its reference coordin
 """
 reflect_xz(panel :: Panel3D) = Panel3D((reflect_xz ∘ p1)(panel), (reflect_xz ∘ p2)(panel), (reflect_xz ∘ p3)(panel), (reflect_xz ∘ p4)(panel))
 
-# Determine axis system of panel
+# Determine panel's local axis system
 function local_coordinate_system(stream, normie) 
     s_hat = normalize(stream)
     n_hat = normalize(normie)
-    c_hat = s_hat × n_hat
+    l_hat = n_hat × s_hat
 
-    return [ c_hat s_hat n_hat ]
+    return [ s_hat l_hat n_hat ]
 end
 
 # Compute local axis coordinates
-local_coordinate_system(panel :: AbstractPanel3D) = local_coordinate_system((panel.p4 - panel.p1 + panel.p3 - panel.p2) / 2, normal_vector(panel))
+local_coordinate_system(panel :: AbstractPanel3D) = local_coordinate_system((panel.p2 + panel.p3) / 2 - midpoint(panel), normal_vector(panel))
 
 function transform_panel(panel :: AbstractPanel3D, point :: Point3D)
     T = get_transformation(panel)
@@ -151,7 +150,7 @@ end
 
 Generate an `AffineMap` to transform a point from global coordinates to an `AbstractPanel3D`'s local coordinate system.
 """
-get_transformation(p) = AffineMap(local_coordinate_system(p)', zeros(3))
+get_transformation(p, P = I(3)) = let T = P * local_coordinate_system(p)'; AffineMap(T, -T * midpoint(p)) end
 
 
 """

src/Tools/MathTools.jl

@@ -10,6 +10,8 @@ struct Point2D{T <: Real} <: FieldVector{2, T}
     y :: T
 end
 
+StaticArrays.similar_type(::Type{<:Point2D}, ::Type{T}, s::Size{(2,)}) where {T} = Point2D{T}
+
 x(p :: Point2D) = p.x
 y(p :: Point2D) = p.y
 
@@ -19,6 +21,8 @@ struct Point3D{T <: Real} <: FieldVector{3, T}
     z :: T
 end
 
+StaticArrays.similar_type(::Type{<:Point3D}, ::Type{T}, s::Size{(3,)}) where {T} = Point3D{T}
+
 x(p :: Point3D) = p.x
 y(p :: Point3D) = p.y
 z(p :: Point3D) = p.z

src/Tools/plot_tools.jl

@@ -1,5 +1,7 @@
-plot_panel(panel :: AbstractPanel3D) = combinedimsview(panel_coordinates(panel), (1))
-
+function plot_panel(panel :: AbstractPanel3D)
+    p = panel_coordinates(panel)
+    combinedimsview([ p; [p[1]]], (1))
+end
 """
     plot_panels(panels :: Array{Panel3D})
 

test/panel_test.jl

@@ -7,43 +7,76 @@ using AeroMDAO
 ##
 panel = Panel3D(
     Point3D( 1.0, -1.0,  0.0), #4
-    Point3D(-1.0, -1.0,  0.0), #3
-    Point3D(-1.0,  1.0,  0.0), #2
+    Point3D(-1.0, -0.9,  0.0), #3
+    Point3D(-1.0,  1.1,  0.0), #2
     Point3D( 1.0,  1.0,  0.0), #1
 )
 
-panel1 = Panel3D(
-    Point3D( 1.0, -1.0,  1.0), #4
-    Point3D(-1.0, -1.0,  1.0), #3
-    Point3D(-1.0,  1.0,  1.0), #2
-    Point3D( 1.0,  1.0,  1.0), #1
-)
-
-point1 = Point3D(0., 0., 10.)
-
-ϵ = 1.0e-5
+# point = Point3D(0.024698438432502412, -0.1209022472975474, -5.551115123125783e-17)
+point = Point3D(2,3,4)
+ϵ = 1.0e-7
 pertx = Point3D(ϵ, 0., 0.)
 perty = Point3D(0., ϵ, 0.)
 pertz = Point3D(0., 0., ϵ)
 
-(quadrilateral_doublet_potential(1, panel, point + pertx) - quadrilateral_doublet_potential(1, panel, point)) / ϵ
-(quadrilateral_doublet_potential(1, panel, point + perty) - quadrilateral_doublet_potential(1, panel, point)) / ϵ
-(quadrilateral_source_potential(1, panel, point + pertz) - quadrilateral_source_potential(1, panel, point)) / ϵ
+(quadrilateral_doublet_potential(panel, point + pertx) - quadrilateral_doublet_potential(panel, point)) / ϵ
+(quadrilateral_doublet_potential(panel, point + perty) - quadrilateral_doublet_potential(panel, point)) / ϵ
+(quadrilateral_doublet_potential(panel, point + pertz) - quadrilateral_doublet_potential(panel, point)) / ϵ
 
-quadrilateral_source_velocity(1, panel, point)
-quadrilateral_source_velocity_farfield(1, panel, point)
+quadrilateral_doublet_velocity(panel, point)
+# quadrilateral_source_velocity_farfield(1, panel, point)
 
-quadrilateral_source_velocity(1, panel, Point3D(0,0,1e-7))
-
-quadrilateral_doublet_potential(1, panel, Point3D(2,0,0))
+##
+plot(xlabel="x", ylabel="y", zlim=(-2,3))
+[plot!(p) for p in plot_panels([
+    panelview[3],
+    panelview[2],
+    panelview[end-2],
+    panelview[end-1],
+])]
+plot!(aspect_ratio=:equal)
 
 ##
+plot(xlabel="x", ylabel="y")
+[plot!(p) for p in plot_panels([
+    surf_pans[1:3];
+    surf_pans[end-2:end]
+])]
+plot!()
 
-plot!(xs(foil_pans[npansp+1]), ys(foil_pans[npansp+1]), zs(foil_pans[npansp+1]), dpi=300, xlabel="x", ylabel="y", zlabel="z")
-plot!(xs(foil_pans[npansp+2]), ys(foil_pans[npansp+2]), zs(foil_pans[npansp+2]), dpi=300, xlabel="x", ylabel="y", zlabel="z")
-plot!(xs(foil_pans[npansp+3]), ys(foil_pans[npansp+3]), zs(foil_pans[npansp+3]), dpi=300, xlabel="x", ylabel="y", zlabel="z")
+##
+plot(xlabel="x", ylabel="y", zlim=(-2,3))
+[plot!(p) for p in plot_panels(
+    panelview[:]
+)]
+plot!(aspect_ratio=:equal)
 
 ##
-plot(xs(thispan), ys(thispan), zs(thispan))
-scatter!([thispoint.x],[thispoint.y],[thispoint.z])
+plot(xlabel="x", ylabel="y", zlim=(-2,3))
+plot!(plot_panels([panel])[1], aspect_ratio=:equal)
+scatter!([point[1]], [point[2]], [point[3]])
+
+x = LinRange(-0.01,0.01,1001)
+y = (x -> quadrilateral_doublet_potential(panel, Point3D(2,3,4+x))).(x)
+plot(x,getindex.(y,3))
+
+av = point - coord[i]
+as = norm(av)
+bv = point - coord[i%4+1]
+bs = norm(bv)
+dv = coord[i%4+1] - coord[i]
+ds = norm(dv)
+al = av ⋅ l
+am = av ⋅ m
+bl = bv ⋅ l
+bm = bv ⋅ m
+dl = dv ⋅ l
+dm = dv ⋅ m
+
+A = dl * (am^2 + gn^2) - al * am * dm
+B = dl * (bm^2 + gn^2) - bl * bm * dm
+
+NUM = A * B + gn^2 * as * bs * dm^2
+DEN = dm * gn * (bs * A - as * B)
 
+atan(NUM/DEN)

test/runtests.jl

@@ -1,7 +1,7 @@
 using AeroMDAO
 using Test
 
-@testset "NACA-4 Doublet-Source Panel Method" begin
+@testset "NACA-4 Doublet-Source 2D Panel Method" begin
     # Define airfoil
     airfoil = (naca4)((0,0,1,2))
 
@@ -24,7 +24,7 @@ using Test
     @test sum(cls_2) ≈ 0.6007449 atol = 1e-6
 end
 
-@testset "Airfoil Processing and Doublet-Source Panel Method" begin
+@testset "Airfoil Processing and Doublet-Source 2D Panel Method" begin
     # Import and read airfoil coordinates
     coo_foil = naca4((2,4,1,2))
 
@@ -93,6 +93,22 @@ end
     @test wing_mac ≈ [0.4209310, 1.3343524, 0.0] atol = 1e-6
 end
 
+@testset "Geometry - 3D Panel" begin
+    panel = Panel3D(
+        [1.0, -1., 0.0], 
+        [0.0, -1., -0.5], 
+        [0.0, 0.0, -0.5], 
+        [1.0, 0.0, 0.0]
+    )
+
+    mp = midpoint(panel)
+    ε = 1.
+    p = SVector(mp.x + ε, mp.y + ε, mp.z + ε)
+    T = get_transformation(panel)
+
+    @test inv(T)(T(mp) + T(p)) ≈ p atol = 1e-6
+end
+
 @testset "Vortex Lattice Method (Incompressible) - NACA 0012 Tapered Wing" begin
     # Define wing
     wing = Wing(foils     = [ naca4((0,0,1,2)) for i ∈ 1:2 ],