Added benchmark and improved surface_integral

AK.mapreduce was really slowing down surface_integral. I replaced it with mapbody!, followed by a simple sum over the panels. It's much faster.

Author: weymouth

examples/benchmark.jl

@@ -0,0 +1,13 @@
+using NeumannKelvin,BenchmarkTools,Test
+
+S(θ,φ) = SA[sin(θ)*cos(φ), sin(θ)*sin(φ), cos(θ)]
+panels = panelize(S, 0, π/2, 0, π, hᵤ=1/64, N_max=Inf) # quite a few panels
+
+@btime extrema(cₚ(directsolve!(BodyPanelSystem($panels,sym_axes=(2,3)),verbose=false)))
+@btime extrema(cₚ(gmressolve!(BodyPanelSystem($panels,sym_axes=(2,3)),verbose=false)))
+@btime extrema(cₚ(gmressolve!(BodyPanelSystem($panels,sym_axes=(2,3),wrap=PanelTree),verbose=false)))
+@btime extrema(cₚ(gmressolve!(BodyPanelSystem($panels,sym_axes=(2,3),wrap=PanelTree,θ²=1),verbose=false)))
+
+sys = gmressolve!(BodyPanelSystem(panels,sym_axes=(2,3),wrap=PanelTree,θ²=1),verbose=false);
+@btime extrema(cₚ($sys))
+@btime steadyforce($sys)

src/panel_method.jl

@@ -79,10 +79,8 @@ See also: [`cₚ`](@ref)
 """
 steadyforce(sys;S=bodyarea(sys)) = -surface_integral(cₚ,sys)/S
 @inline function surface_integral(f,sys)
-    init = neutral = zero(eltype(sys.body.n))
-    AK.mapreduce(+, sys.body, AK.get_backend(sys.body.q); init, neutral) do p
-        f(p.x,sys) * p.n * p.dA
-    end
+    a = mapbody!(f,similar(sys.body.q),sys)
+    sum(i->a[i]*sys.body.n[i]*sys.body.dA[i],eachindex(a)) # sum quick
 end
 
 """