add files from flowfmm branch

src/fmm.jl

@@ -0,0 +1,120 @@
+Base.getindex(sys::AbstractPanels, i, ::FastMultipole.Body) = sys.panels[i], sys.potential[i], sys.velocity[i], sys.normals[i]
+
+Base.getindex(sys::AbstractPanels, i, ::FastMultipole.Position) = sys.panels[i].control_point
+
+Base.getindex(sys::AbstractPanels, i, ::FastMultipole.Radius) = sys.panels[i].radius
+
+Base.getindex(sys::AbstractPanels, i, ::FastMultipole.ScalarPotential) = sys.potential[i]
+
+Base.getindex(sys::AbstractPanels, i, ::FastMultipole.Velocity) = sys.velocity[i]
+
+Base.getindex(sys::AbstractPanels, i, ::FastMultipole.VelocityGradient) = zero(SMatrix{3,3,Float64,3})
+
+Base.getindex(sys::AbstractPanels{ConstantNormalDoublet(),<:Any}, i, ::FastMultipole.Strength) = sys.panels[i].strength[1]
+Base.getindex(sys::AbstractPanels{ConstantSource(),<:Any}, i, ::FastMultipole.Strength) = sys.panels[i].strength[1]
+Base.getindex(sys::AbstractPanels, i, ::FastMultipole.Strength) = sys.panels[i].strength
+
+Base.getindex(sys::AbstractPanels, i, ::FastMultipole.Vertex, i_vertex) = sys.panels[i].vertices[i_vertex]
+
+Base.getindex(sys::AbstractPanels, i, ::FastMultipole.Normal) = sys.panels[i].normal
+
+function Base.setindex!(sys::AbstractPanels, val, i, ::FastMultipole.Body)
+    panel, potential, velocity, normal = val
+    sys.panels[i] = panel
+    sys.potential[i] = potential
+    sys.velocity[i] = velocity
+    sys.normals[i] = normal
+end
+
+Base.setindex!(sys::AbstractPanels, val, i, ::FastMultipole.ScalarPotential) = sys.potential[i] = val
+
+Base.setindex!(::AbstractPanels, val, i, ::FastMultipole.VectorPotential) = nothing
+
+Base.setindex!(sys::AbstractPanels, val, i, ::FastMultipole.Velocity) = sys.velocity[i] = val
+
+Base.setindex!(::AbstractPanels, val, i, ::FastMultipole.VelocityGradient) = nothing
+
+function Base.setindex!(system::AbstractPanels{<:Any,TF,NK,NS}, val, i, ::FastMultipole.Strength) where {TF,NK,NS}
+    panel = system.panels[i]
+    vertices = panel.vertices
+    control_point = panel.control_point
+    normal = panel.normal
+    radius = panel.radius
+    strength = SVector{NK,TF}(val)
+    system.panels[i] = Panel{TF,NK,NS}(vertices, control_point, normal, strength, radius)
+end
+
+FastMultipole.get_n_bodies(sys::AbstractPanels) = length(sys.panels)
+
+FastMultipole.buffer_element(sys::AbstractPanels) = deepcopy(sys.panels[1]), zero(eltype(sys.potential)), zero(eltype(sys.velocity)), zero(eltype(sys.normals))
+
+Base.eltype(::AbstractPanels{<:Any,TF,<:Any,<:Any}) where TF = TF
+
+FastMultipole.B2M!(system::AbstractPanels{ConstantSource(),<:Any,1,4}, branch, bodies_index, harmonics, expansion_order) = FastMultipole.B2M!_quadpanel(system, branch, bodies_index, harmonics, expansion_order, FastMultipole.UniformSourcePanel())
+
+FastMultipole.B2M!(system::AbstractPanels{ConstantSource(),<:Any,1,3}, branch, bodies_index, harmonics, expansion_order) = FastMultipole.B2M!_tripanel(system, branch, bodies_index, harmonics, expansion_order, FastMultipole.UniformSourcePanel())
+
+FastMultipole.B2M!(system::AbstractPanels{ConstantNormalDoublet(),<:Any,1,4}, branch, bodies_index, harmonics, expansion_order) = FastMultipole.B2M!_quadpanel(system, branch, bodies_index, harmonics, expansion_order, FastMultipole.UniformNormalDipolePanel())
+
+FastMultipole.B2M!(system::AbstractPanels{ConstantNormalDoublet(),<:Any,1,3}, branch, bodies_index, harmonics, expansion_order) = FastMultipole.B2M!_tripanel(system, branch, bodies_index, harmonics, expansion_order, FastMultipole.UniformNormalDipolePanel())
+
+@inline function convolve_kernel!(target_system, target_index, panel, kernel::AbstractRotatedKernel, derivatives_switch::DerivativesSwitch{PS,<:Any,VS,GS}) where {PS,VS,GS}
+    # rotate into source panel frame
+    Rprime, Rxprime, Ryprime, Rzprime = rotate_to_panel(panel)
+
+    # iterate over targets
+    for i_target in target_index
+        potential, velocity, gradient = _induced(target_system[i_target, FastMultipole.POSITION], panel, kernel, Rprime, Rxprime, Ryprime, Rzprime, derivatives_switch)
+        if PS
+            if FastMultipole.DEBUG_TOGGLE[]
+                println("before:")
+                @show target_system[i_target, FastMultipole.SCALAR_POTENTIAL]
+            end
+            target_system[i_target, FastMultipole.SCALAR_POTENTIAL] += potential
+            if FastMultipole.DEBUG_TOGGLE[]
+                println("after:")
+                @show target_system[i_target, FastMultipole.SCALAR_POTENTIAL]
+            end
+        end
+        # target_system[i_target, FastMultipole.SCALAR_POTENTIAL] = target_system[i_target, FastMultipole.SCALAR_POTENTIAL] + potential
+        if VS
+            target_system[i_target, FastMultipole.VELOCITY] += velocity
+        end
+        # target_system[i_target, FastMultipole.VELOCITY] = target_system[i_target, FastMultipole.VELOCITY] + velocity
+        if GS
+            target_system[i_target, FastMultipole.VELOCITY_GRADIENT] += gradient
+        end
+        # target_system[i_target, FastMultipole.VELOCITY_GRADIENT] = target_system[i_target, FastMultipole.VELOCITY_GRADIENT] + gradient
+    end
+
+    if FastMultipole.DEBUG_TOGGLE[]
+        println("after after:")
+        @show target_system[1, FastMultipole.SCALAR_POTENTIAL]
+    end
+end
+
+@inline function convolve_kernel!(target_system, target_index, panel, kernel::AbstractUnrotatedKernel, derivatives_switch::DerivativesSwitch{PS,<:Any,VS,GS}) where {PS,VS,GS}
+    # iterate over targets
+    for i_target in target_index
+        potential, velocity, gradient = _induced(target_system[i_target, FastMultipole.POSITION], panel, kernel, derivatives_switch::DerivativesSwitch{PS,<:Any,VS,GS}) where {PS,VS,GS}
+        if PS
+            target_system[i_target, FastMultipole.SCALAR_POTENTIAL] = target_system[i_target, FastMultipole.SCALAR_POTENTIAL] + potential
+        end
+        if VS
+            target_system[i_target, FastMultipole.VELOCITY] = target_system[i_target, FastMultipole.VELOCITY] + velocity
+        end
+        if GS
+            target_system[i_target, FastMultipole.VELOCITY_GRADIENT] = target_system[i_target, FastMultipole.VELOCITY_GRADIENT] + gradient
+        end
+    end
+end
+
+@inline function FastMultipole.direct!(target_system, target_index, derivatives_switch, source_system::AbstractPanels{kernel,<:Any,<:Any,<:Any}, source_index) where kernel
+    if FastMultipole.DEBUG_TOGGLE[]
+        @show source_index target_index
+    end
+    for i_panel in source_index
+        convolve_kernel!(target_system, target_index, source_system.panels[i_panel], kernel, derivatives_switch)
+    end
+    return nothing
+end

src/fmm_FLOWPanel.jl

@@ -0,0 +1,92 @@
+"""
+  Three-dimensional panel method for high-Reynolds aerodynamics.
+
+  # AUTHORSHIP
+    * Created by    : Eduardo J. Alvarez
+    * Email         : [email protected]
+    * Date          : Jun 2018 originally as MyPanel.jl
+    * License       : MIT License
+"""
+module FLOWPanel
+
+export  solve, save, Uind!, phi!,
+        get_ndivscells, get_ndivsnodes,
+        get_cart2lin_cells, get_cart2lin_nodes,
+        get_field, get_fieldval, add_field,
+        calc_normals!, calc_normals,
+        calc_tangents!, calc_tangents,
+        calc_obliques!, calc_obliques,
+        calc_controlpoints!, calc_controlpoints,
+        calc_areas!, calc_areas
+
+# ------------ GENERIC MODULES -------------------------------------------------
+import Dierckx
+import LinearAlgebra as LA
+import LinearAlgebra: I
+import Krylov
+import Requires: @require
+
+# ------------ FLOW LAB MODULES ------------------------------------------------
+# GeometricTools from https://github.com/byuflowlab/GeometricTools.jl
+import GeometricTools
+const gt = GeometricTools
+import GeometricTools: Meshes
+import ImplicitAD as IAD
+import ImplicitAD: ForwardDiff as FD, ReverseDiff as RD
+
+# ------------ GLOBAL VARIABLES AND DATA STRUCTURES ----------------------------
+const module_path = splitdir(@__FILE__)[1]      # Path to this module
+                                                # Default path to data files
+const def_data_path = joinpath(module_path, "..", "docs", "resources", "data")
+                                            # Default path to airfoil data files
+const def_rfl_path = joinpath(def_data_path, "airfoils")
+                                                # Path to examples
+const examples_path = joinpath(module_path, "..", "examples")
+
+# const RType = Union{Float64,                    # Concrete real types
+#                     Int64,
+#                     # ForwardDiff.Dual{Nothing,Float64,3},
+#                     # ForwardDiff.Dual{Nothing,Int64,3}
+#                     }
+
+# Discretization parameter type
+const ndivstype = Union{Float64, gt.multidisctype, Nothing}
+
+# Identity matrix
+const Im = Array(1.0I, 3, 3)
+
+# Shedding matrix for a RigidWakeBody without shedding
+const noshedding = zeros(Int, 6, 0)
+
+# ------------ HEADERS ---------------------------------------------------------
+for header_name in ["elements", "linearsolver",
+                    "abstractbody", "nonliftingbody",
+                    "abstractliftingbody", "liftingbody",
+                    "multibody",
+                    "utils", "postprocess"
+                    ]
+  include("FLOWPanel_"*header_name*".jl")
+end
+
+# Conditionally load monitors if PyPlot is available
+function __init__()
+
+    try
+        @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" begin
+
+            import .PyPlot as plt
+            import .PyPlot: @L_str
+
+            for header_name in ["monitor"]
+              include("FLOWPanel_"*header_name*".jl")
+            end
+
+        end
+
+    catch e
+        @warn "PyPlot is not available; monitors will not be loaded"
+    end
+
+end
+
+end # END OF MODULE

src/freestream.jl

@@ -0,0 +1,20 @@
+function apply_freestream!(panel_array::AbstractPanels, freestream::AbstractVector{<:Number})
+    # unpack panels
+    # (; panels, potential, velocity) = panel_array
+    panels = panel_array.panels
+    potential = panel_array.potential
+    velocity = panel_array.velocity
+
+    # potential
+    for (i,panel) in enumerate(panels)
+        x = panel.control_point
+        potential[i] = -dot(x, freestream)
+    end
+
+    # velocity
+    for i in eachindex(velocity)
+        velocity[i] += freestream
+    end
+
+    return nothing
+end

src/geometry.jl

@@ -0,0 +1,57 @@
+#####
+##### general geometric helper functions
+#####
+@inline function compute_centroid(p1,p2,p3)
+    return (p1+p2+p3) ./ 3
+end
+
+function compute_centroid_normal_area(points...)
+    n_points = length(points)
+    p1 = points[1]
+    centroid = @SVector zeros(eltype(p1),3)
+    normal = @SVector zeros(eltype(p1),3)
+    area = zero(eltype(p1))
+    for i_triangle in 2:n_points-1
+        v1 = points[i_triangle] - p1
+        v2 = points[i_triangle+1] - p1
+        this_cross = cross(v1,v2)
+        this_cross_norm = norm(this_cross)
+        this_area = this_cross_norm/2
+        centroid = centroid + this_area * compute_centroid(points[1], points[i_triangle], points[i_triangle+1])
+        normal = normal + this_cross * this_area
+        area += this_area
+    end
+    return centroid ./ area, normal ./ norm(normal), area
+end
+
+#####
+##### generators for common or test geometries
+#####
+function generate_cube(side_length, centroid, kernel)
+    p1x = centroid[1] - side_length/2
+    p1y = centroid[2] - side_length/2
+    p1z = centroid[3] - side_length/2
+    points = [
+        SVector{3}(p1x,p1y,p1z),
+        SVector{3}(p1x+side_length,p1y,p1z),
+        SVector{3}(p1x+side_length,p1y+side_length,p1z),
+        SVector{3}(p1x,p1y+side_length,p1z),
+        SVector{3}(p1x,p1y,p1z+side_length),
+        SVector{3}(p1x+side_length,p1y,p1z+side_length),
+        SVector{3}(p1x+side_length,p1y+side_length,p1z+side_length),
+        SVector{3}(p1x,p1y+side_length,p1z+side_length),
+    ]
+    cells = [
+        MeshCell(VTKCellTypes.VTK_QUAD,[1,4,3,2]),
+        MeshCell(VTKCellTypes.VTK_QUAD,[1,2,6,5]),
+        MeshCell(VTKCellTypes.VTK_QUAD,[5,6,7,8]),
+        MeshCell(VTKCellTypes.VTK_QUAD,[5,8,4,1]),
+        MeshCell(VTKCellTypes.VTK_QUAD,[4,8,7,3]),
+        MeshCell(VTKCellTypes.VTK_QUAD,[3,7,6,2]),
+    ]
+    return PanelArray(points, cells, kernel)
+end
+
+function generate_plane()
+
+end