Merge pull request #24 from byuflowlab/fastmultipole-mvp

Fastmultipole MVP and multipole new features for unstructured grids

Author: EdoAlvarezR

.github/workflows/docs.yaml

@@ -22,7 +22,7 @@ jobs:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
         with:
-          version: '1.10'
+          version: '1.11'
       - uses: julia-actions/cache@v1
       - name: Override Conda Python
         run: if [ "$RUNNER_OS" = "Linux" ]; then julia -e 'import Pkg; Pkg.add("PyCall"); ENV["PYTHON"]="/usr/bin/python3"; Pkg.build("PyCall")'; fi

.github/workflows/test.yaml

@@ -7,7 +7,7 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        julia-version: ['1.8', '1.10']
+        julia-version: ['1.8', '1.10', '1.11']
         julia-arch: [x64]
         os: [ubuntu-latest, macOS-latest]
 
@@ -24,3 +24,5 @@ jobs:
         shell: bash
       # - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
+        with:
+          test_args: 'remote'

Project.toml

@@ -1,33 +1,38 @@
 name = "FLOWPanel"
 uuid = "6be8c882-484d-4309-b349-d23112750151"
 author = ["Eduardo Alvarez <Edo.AlvarezR@gmail.com>"]
-version = "1.1.1"
+version = "1.2.1"
 
 [deps]
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
-Requires = "ae029012-a4dd-5104-9daa-d747884805df"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
+# FastMultipole = "ce07d0d3-2b9f-49ba-89eb-12c800257c85"
 GeometricTools = "83792f5e-c6a1-11e8-2e0a-93511f02ae5f"
 ImplicitAD = "e7cbb90b-9b31-4eb2-a8c8-45099c074ee1"
+Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearOperators = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125"
+Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
+
+# [sources]
+# FastMultipole = {url = "https://github.com/byuflowlab/FastMultipole", rev = "d21d242d081ea8c9c19655097d83413b1f783d2d"}
 
 [compat]
-julia = "1.6"
 Dierckx = "0.5"
 GeometricTools = "2.2"
 ImplicitAD = "0.3"
 Krylov = "0.9"
 Requires = "1"
+julia = "1.6"
 
 [extras]
-Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Pkg", "Test", "Printf",
-        "ForwardDiff",
-        "CSV", "DataFrames"]
+test = ["Pkg", "Test", "Printf", "ForwardDiff", "CSV", "DataFrames"]

docs/src/examples/cessna-aero.md

@@ -37,8 +37,6 @@ import Meshes
 import GeoIO
 import Rotations: RotX, RotY, RotZ
 
-# import CUDA                               # Uncomment this to use GPU (if available)
-
 
 run_name        = "cessna"                  # Name of this run
 
@@ -53,20 +51,20 @@ rho             = 0.71461                   # (kg/m^3) air density at 17300 ft
 
 
 # ----------------- GEOMETRY DESCRIPTION ---------------------------------------
-meshfile        = joinpath(read_path, "cessna.msh")    # Gmsh file to read
+meshfile        = "cessna.msh"              # Gmsh file to read
 
 offset          = [0, 0, 0]                 # Offset to center the mesh
 rotation        = RotX(0*pi/180) * RotY(0*pi/180) * RotX(0*pi/180) # Rotation to align mesh
 scaling         = 0.3048                    # Factor to scale mesh dimensions (ft -> m)
 
-                                            # Gmsh files with trailing edges
-trailingedges   = [ #   ( Gmsh file,                 span direction )
-                        ("cessna-TE-leftwing.msh",      [0, 1, 0]),
-                        ("cessna-TE-rightwing.msh",     [0, 1, 0]),
-                        ("cessna-TE-leftelevator.msh",  [0, 1, 0]),
-                        ("cessna-TE-rightelevator.msh", [0, 1, 0]),
-                        ("cessna-TE-rudder.msh",        [0, 0, 1]),
-                    ]
+trailingedges   = [                         # Gmsh files with trailing edges
+#   (  Gmsh file,                    span sorting function,   junction criterion,     closed )
+    ( "cessna-TE-leftwing.msh",      X -> dot(X, [0, 1, 0]),  X -> abs(X[2]) - 0.67, false  ),
+    ( "cessna-TE-rightwing.msh",     X -> dot(X, [0, 1, 0]),  X -> abs(X[2]) - 0.67, false  ),
+    ( "cessna-TE-leftelevator.msh",  X -> dot(X, [0, 1, 0]),  X -> abs(X[2]) - 0.23, false  ),
+    ( "cessna-TE-rightelevator.msh", X -> dot(X, [0, 1, 0]),  X -> abs(X[2]) - 0.23, false  ),
+    ( "cessna-TE-rudder.msh",        X -> dot(X, [0, 0, 1]),  X -> X[3] - 0.65,      false  ),
+]
 
 flip            = true                      # Whether to flip control points against the direction of normals
                                             # NOTE: use `flip=true` if the normals
@@ -78,6 +76,12 @@ Sref            = 16.2344                   # (m^2) reference area
 Xac             = [2.815, 0, 0.4142]        # (m) aerodynamic center for
                                             # calculation of moments
 
+# Define function used for reading Gmsh files
+meshreader(file) = GeoIO.load(file).geometry
+
+# Format input for `generate_multibody(...)`
+meshfiles       = [ ("Airframe", meshfile, flip) ]
+
 # ----------------- SOLVER SETTINGS -------------------------------------------
 
 # Solver: direct linear solver for open bodies
@@ -88,71 +92,11 @@ bodytype        = pnl.RigidWakeBody{pnl.VortexRing, 2}
 
 
 # ----------------- GENERATE BODY ----------------------------------------------
-# Read Gmsh mesh
-msh = GeoIO.load(meshfile)
-msh = msh.geometry
-
-# Transform the original mesh: Translate, rotate, and scale
-msh = msh |> Meshes.Translate(offset...) |> Meshes.Rotate(rotation) |> Meshes.Scale(scaling)
-
-# Wrap Meshes object into a Grid object from GeometricTools
-grid = pnl.gt.GridTriangleSurface(msh)
-
-# Read all trailing edges
-sheddings = []
-
-for (trailingedgefile, spandir) in trailingedges
-
-    # Read Gmsh line of trailing edge
-    TEmsh = GeoIO.load(joinpath(read_path, trailingedgefile))
-    TEmsh = TEmsh.geometry
-
-    # Apply the same transformations of the mesh to the trailing edge
-    TEmsh = TEmsh |> Meshes.Translate(offset...) |> Meshes.Rotate(rotation) |> Meshes.Scale(scaling)
-
-    # Convert TE Meshes object into a matrix of points used to identify the trailing edge
-    trailingedge = pnl.gt.vertices2nodes(TEmsh.vertices)
-
-    # Sort TE points from "left" to "right" according to span direction
-    trailingedge = sortslices(trailingedge; dims=2, by = X -> pnl.dot(X, spandir))
-
-    # Function for identifying if a point is close to a junction
-    if trailingedgefile in ["cessna-TE-leftwing.msh", "cessna-TE-rightwing.msh"]
-
-        criterion = X -> abs(X[2]) <= 0.67
-
-    elseif trailingedgefile in ["cessna-TE-leftelevator.msh", "cessna-TE-rightelevator.msh"]
-
-        criterion = X -> abs(X[2]) <= 0.23
-
-    elseif trailingedgefile in ["cessna-TE-rudder.msh"]
-
-        criterion = X -> X[3] <= 0.65
-
-    else
-
-        criterion = X -> false
-
-    end
-
-    # Filter out any points that are close to junctions
-    tokeep = filter( i -> !criterion(trailingedge[:, i]), 1:size(trailingedge, 2) )
-    trailingedge = trailingedge[:, tokeep]
-
-    # Generate TE shedding matrix
-    shedding = pnl.calc_shedding(grid, trailingedge; tolerance=0.001*bref)
-
-    push!(sheddings, shedding)
-
-end
-
-# Combine all TE shedding matrices into one
-shedding = hcat(sheddings...)
-
-# Generate paneled body
-body = bodytype(grid, shedding; CPoffset=(-1)^flip * 1e-14)
-
-println("Number of panels:\t$(body.ncells)")
+body, elsprescribe = pnl.generate_multibody(bodytype, meshfiles, trailingedges, meshreader;
+                                tolerance=0.001*bref,
+                                read_path=read_path,
+                                offset=offset, rotation=rotation, scaling=scaling,
+                                )
 
 
 # ----------------- CALL SOLVER ------------------------------------------------
@@ -171,10 +115,7 @@ Dbs = repeat(Vinf/magVinf, 1, body.nsheddings)
 
 # Solve body (panel strengths) giving `Uinfs` as boundary conditions and
 # `Das` and `Dbs` as trailing edge rigid wake direction
-@time pnl.solve(body, Uinfs, Das, Dbs)
-
-# Uncomment this to use GPU instead (if available)
-# @time pnl.solve(body, Uinfs, Das, Dbs; GPUArray=CUDA.CuArray{Float32})
+@time pnl.solve(body, Uinfs, Das, Dbs; elprescribe=elsprescribe)
 
 # ----------------- POST PROCESSING ----------------------------------------
 println("Post processing...")

docs/src/generate_examples_cessna.jl

@@ -203,7 +203,7 @@ open(joinpath(output_path, output_name*"-aero.md"), "w") do fout
 
     println(fout, "```julia")
 
-    open(joinpath(pnl.examples_path, "cessna.jl"), "r") do fin
+    open(joinpath(pnl.examples_path, "cessna2.jl"), "r") do fin
         for l in eachline(fin)
             # if contains(l, "COMPARISON TO EXPERIMENTAL DATA")
             #     break