Add submerge key-word to panelize

I realized that it's pretty easy to find the z=0 point on a given strip and force strip below this point. This will let people trim ship hulls very easily. Also fixed up a bad test and the Wigley hull README example.

Author: weymouth

.gitignore

@@ -1 +1,2 @@
 /Manifest.toml
+.vscode/settings.json

README.md

@@ -278,9 +278,10 @@ See the docs for `NKPanelSystem`, `kelvin`, and the cited references for details
 
 As a final application, let's simulate the classic Wigley hull with `NKPanels`
 ```julia
-wigley(hᵤ;B=1/8,D=1/16,hᵥ=0.5hᵤ/D) = measure.(
-    (u,v)->SA[u-0.5,-2B*u*(1-u)*(v)*(2-v),D*(v-1)],
-    0.5hᵤ:hᵤ:1,(0.5hᵥ:hᵥ:1)',hᵤ,hᵥ)
+function wigley(hᵤ,hᵥ=hᵤ;B=1/8,D=1/6,kwargs...)
+    S(u,v) = SA[u-0.5,-2B*u*(1-u)*(v)*(2-v),D*(v-1)]
+    panelize(S;hᵤ,hᵥ,kwargs...)
+end
 NKsys = NKPanelSystem(wigley(0.025);ℓ=1/2π,sym_axes=2,contour=true) |> directsolve!
 viz(NKsys)
 ```

examples/NDLPwriteup.jl

@@ -143,9 +143,9 @@ The NDLP segmentation consistently produces a max deviation which is tight to th
 
 # ╔═╡ b38c2b1c-a48c-4f1f-954b-95faaba680d7
 begin
-	function metrics(method, r, u₀, u₁, Δs, dₙ)
+	function metrics(method, r, u₀, u₁, N′, dₙ)
 		L = seg_len(r, u₀, u₁)            # total length of curve `r`
-		Δs *= L                           # scale by L to compare across curves
+		Δs = L/N′                         # segment length limit for `r`
 		u = method(r, u₀, u₁, Δs, dₙ)     # sample points using `method`
 		dl,δ = seg_len(r,u),seg_dev(r,u)  # segment lengths and deviations
 		return (;δ∞ = maximum(δ)/(Δs*dₙ), # scaled max deviation, should => 1!
@@ -189,7 +189,7 @@ begin
 	end
 	function NDLP(r, u₀, u₁, Δs, dₙ) # one-shot sampling! 🤓
 		speed(u) = max(arcspeed(r)(u),√(Δs*aₙ(r,u)/8dₙ))
-		rtol = 1e-6Δs # only needed since convergence study lets Δs→0
+		rtol = 1e-6Δs # only needed since sensitivity study lets Δs→0
 		S,s⁻¹ = ∫speed(speed, u₀, u₁; rtol)
 		return s⁻¹.(range(0, S, round(Int,S/Δs)+1))
 	end
@@ -220,21 +220,21 @@ end; plot!(aspect_ratio=:equal,xlabel="x",ylabel="y")
 
 # ╔═╡ f7a8cdf1-2cdb-4d6b-a10f-3f2d980c836e
 let
-	invΔs = 33; dₙ = 1/100; # held constant for all methods/test_curves
-	println("Metrics using Δs = L/$invΔs and dₙ=$dₙ")
+	N′ = 33; dₙ = 1/100; # held constant for all methods/test_curves
+	println("Metrics using Δs = L/$N′ and dₙ=$dₙ")
 	for method in (subdivision,κ_weighted,NDLP)
 		println("\nMethod: $method")
 		map(test_curves) do (name, r, range, _)
-			(;name,metrics(method,r,range...,1/invΔs,dₙ)...)
+			(;name,metrics(method,r,range...,N′,dₙ)...)
 		end |> Table |> display
 	end
 end
 
 # ╔═╡ 561a1373-57a4-4b67-be5a-9d94c9496360
 md"""
-## Convergence study
+## Sampling sensitivity study
 
-We also evaluate the convergence of the NDLP segmentation metrics as the $\Delta s$ and $d_n$ limits vary. We use the cubic spline fish as a representative curve.
+We also evaluate the sensitivity of the NDLP sampling as the $\Delta s$ and $d_n$ limits vary. We use the cubic spline fish as a representative curve.
 
 Holding $d_n=1\%$ constant and *reducing* $\Delta s$ shows two distance phases.
  - In the first phase, the deviation $\max(\delta)$ goes rapidly to the limit $d_n\Delta s$ and holds steady while the excess number of segments and total variation in the panel lengths drops to zero with $\Delta s$.
@@ -244,14 +244,14 @@ Holding $d_n=1\%$ constant and *reducing* $\Delta s$ shows two distance phases.
 # ╔═╡ dcb68886-2771-4b37-91b6-9c983e118728
 let
 	dₙ = 1e-2
-	convergeΔs = map(logrange(1e-1,1e-4,70)) do Δs
-		(;Δs,metrics(NDLP,fish_spline,0,1,Δs,dₙ)...)
+	N′sweep = map(logrange(10,10000,70)) do N′
+		(;N′,metrics(NDLP,fish_spline,0,1,N′,dₙ)...)
 	end |> Table
-	plot(convergeΔs.Δs,convergeΔs.δ∞,label="max(δ)/dₙΔs")
-	plot!(convergeΔs.Δs,convergeΔs.σ,label="NΔs/L-1")
-	plot!(convergeΔs.Δs,convergeΔs.Rₜᵥ,label="sum(R)/L")
-	plot!(convergeΔs.Δs,convergeΔs.R∞,label="max(R)/Δs")
-	plot!(ylabel="scaled metrics",xlabel="Δs",xscale=:log10,xflip=true)
+	plot(N′sweep.N′,N′sweep.δ∞,label="max(δ)/dₙΔs")
+	plot!(N′sweep.N′,N′sweep.σ,label="NΔs/L-1")
+	plot!(N′sweep.N′,N′sweep.Rₜᵥ,label="sum(R)/L")
+	plot!(N′sweep.N′,N′sweep.R∞,label="max(R)/Δs")
+	plot!(ylabel="scaled metrics",xlabel="L/Δs",xscale=:log10)
 	plot!(title="NDLP segmentation metrics with fixed dₙ=$dₙ")
 end
 
@@ -264,16 +264,16 @@ Holding $\Delta s=L/100$ and _increasing_ $d_n$ shows a similar trend.
 
 # ╔═╡ f67fc432-8031-4230-899b-b24fcb7b44f4
 let
-	invΔs = 100
-	convergedₙ = map(logrange(1,5e-4,100)) do dₙ
-		(;dₙ,metrics(NDLP,fish_spline,0,1,1/invΔs,dₙ)...)
+	N′ = 100
+	dₙsweep = map(logrange(1,5e-4,100)) do dₙ
+		(;dₙ,metrics(NDLP,fish_spline,0,1,N′,dₙ)...)
 	end |> Table
-	plot(convergedₙ.dₙ,convergedₙ.δ∞,label="max(δ)/dₙΔs")
-	plot!(convergedₙ.dₙ,convergedₙ.σ,label="NΔs/L-1")
-	plot!(convergedₙ.dₙ,convergedₙ.Rₜᵥ,label="sum(R)/L")
-	plot!(convergedₙ.dₙ,convergedₙ.R∞,label="max(R)/Δs")
-	plot!(ylabel="scaled metrics",xlabel="dₙ",xscale=:log10)
-	plot!(title="NDLP segmentation metrics with fixed Δs=L/$invΔs")
+	plot(dₙsweep.dₙ,dₙsweep.δ∞,label="max(δ)/dₙΔs")
+	plot!(dₙsweep.dₙ,dₙsweep.σ,label="NΔs/L-1")
+	plot!(dₙsweep.dₙ,dₙsweep.Rₜᵥ,label="sum(R)/L")
+	plot!(dₙsweep.dₙ,dₙsweep.R∞,label="max(R)/Δs")
+	plot!(ylabel="scaled metrics",xlabel="dₙ",xscale=:log10,ylims=(-0.04,1.8))
+	plot!(title="NDLP segmentation metrics with fixed Δs=L/$N′")
 end
 
 # ╔═╡ 6094d19e-e64a-4434-b6db-e5647fd71e78
@@ -309,7 +309,7 @@ TypedTables = "~1.4.6"
 PLUTO_MANIFEST_TOML_CONTENTS = """
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.11.5"
+julia_version = "1.11.2"
 manifest_format = "2.0"
 project_hash = "84beccf13d730d2ef8c64387b8e5b55843d0e9dd"
 
@@ -1082,7 +1082,7 @@ version = "0.3.27+1"
 [[deps.OpenLibm_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
-version = "0.8.5+0"
+version = "0.8.1+2"
 
 [[deps.OpenSSL]]
 deps = ["BitFlags", "Dates", "MozillaCACerts_jll", "NetworkOptions", "OpenSSL_jll", "Sockets"]

examples/readme.jl

@@ -0,0 +1,57 @@
+using NeumannKelvin
+h = 0.1 # spacing
+freesurf = measure.((u,v)->SA[u,v,0],2:-h:-4,(2:-h:-2)',h,h)
+S(θ₁,θ₂) = SA[0.5cos(θ₁),-0.1sin(θ₂)*sin(θ₁),0.1cos(θ₂)*sin(θ₁)-0.15]
+body = panelize(S,0,π,0,2π,hᵤ=h)
+sys = BodyPanelSystem(body) |> directsolve!
+
+bigbody = panelize(S,0,π,0,2π,hᵤ=1/200,N_max=Inf); #20k panels
+sys = BodyPanelSystem(bigbody,wrap=PanelTree) |> gmressolve!
+
+halfbody = panelize(S,0,π,0,π,hᵤ=1/200,N_max=Inf); # half the θ₂ range
+BodyPanelSystem(halfbody,wrap=PanelTree,sym_axes=2) |> gmressolve!
+
+directsolve!(BodyPanelSystem(body,U=SA[0,-1,0]),verbose=false) |> addedmass
+
+ℓ = 1/4; h = 0.3ℓ # Froude-length and spacing
+# Make the free-surface grid using y-symmetry and resolving ℓ
+freesurf = measure.((u,v)->SA[u,-v,0],2:-h:-4,(h/2:h:2)',h,h)
+halfbody = panelize(S,0,π,0,π,hᵤ=h)
+
+FSsys = FSPanelSystem(halfbody,freesurf;
+                  ℓ,sym_axes=2,θ²=16) |> gmressolve!
+steadyforce(FSsys)
+
+# Read mesh, place it, and filter out panels intersecting z=0
+using GeometryBasics,FileIO
+function affine(mesh, A, b)  # rotate, scale, and shift the mesh
+    position = [Point3f(A * p + b) for p in mesh.position]
+    GeometryBasics.Mesh(position, mesh.faces)
+end
+dolphin = let 
+    mesh = load("examples//LowPolyDolphin.stl")
+    mesh = affine(mesh, SA[0 -1 0;1 0 0;0 0 1]/65,SA[0.043,0,-0.09])
+    filter(p->abs(p.x[3])>0.01, panelize(mesh))
+end
+
+# The rest of the system matches the example above
+h = 0.04; freesurf = measure.((u,v)->SA[u,-v,0],2/3:-h:-4/3,(-2/3:h:2/3)',h,h,T=Float32); # Float32 to match the Mesh
+Meshsys = FSPanelSystem(dolphin,freesurf;ℓ=9f-2) |> gmressolve!
+
+using GLMakie # can also use Plots, or WGLMakie (for browsers)
+viz(Meshsys)
+
+ℓ = 1/4; h = 0.3ℓ # Froude-length and spacing
+halfbody = panelize(S,0,π,0,π,hᵤ=h)
+NKsys = NKPanelSystem(halfbody;ℓ,sym_axes=2) |> directsolve!
+steadyforce(NKsys)
+
+DBsys = BodyPanelSystem(halfbody;sym_axes=(2,3)) |> directsolve! |> steadyforce
+
+function wigley(hᵤ,hᵥ=hᵤ;B=1/8,D=1/6,kwargs...)
+    S(u,v) = SA[u-0.5,-2B*u*(1-u)*(v)*(2-v),D*(v-1)]
+    panelize(S;hᵤ,hᵥ,kwargs...)
+end
+NKsys = NKPanelSystem(wigley(0.025);
+          ℓ=1/2π,sym_axes=2,contour=true) |> directsolve!
+viz(NKsys)

src/panels.jl

@@ -1,17 +1,18 @@
 """
-    panelize(surface,u₀=0,u₁=1,v₀=0,v₁=1;hᵤ=1,hᵥ=hᵤ,devlimit=0.05,transpose=false,flip=false,N_max=1000,kwargs...)
+    panelize(surface,u₀=0,u₁=1,v₀=0,v₁=1;hᵤ=1,hᵥ=hᵤ,devlimit=0.05,
+            transpose=false,flip=false,submerge=false,N_max=1000,kwargs...)
 
 Panelize a parametric `surface` of `u∈[u₀,u₁]` and `v∈[v₀,v₁]`, returning a `Table` of panels.
 
 The surface is split into strips roughly `hᵤ` wide which are split into panels roughly `hᵥ` high which are then `measure`d.
-Use `transpose=true` to change the strip direction and `flip=true` to flip the normal direction.
-The parameter `devlimit` sets the max deviation of a flat panel from the surface by reducing panel size in regions of high curvature.
+The parameter `devlimit` sets the max deviation `δ/(hᵤ+hᵥ)` from the surface by reducing panel size in regions of high curvature.
+Use `transpose=true` to change the strip direction, `flip=true` to flip the normal direction, and submerge=`true` to trim surface for `z≥0`.
 This function throws an error if the adaptive routine gives more than `N_max` panels.
 """
 function panelize(surface,u₀=0.,u₁=1.,v₀=0.,v₁=1.;hᵤ=1.,hᵥ=hᵤ,devlimit=0.05,
-                  transpose=false,flip=false,N_max=1000,verbose=false,T=Float64,kwargs...)
+                  transpose=false,flip=false,N_max=1000,verbose=false,submerge=false,T=Float64,kwargs...)
     # Transpose arguments u,v -> v,u
-    transpose && return panelize((v,u)->surface(u,v),v₀,v₁,u₀,u₁,;hᵤ=hᵥ,hᵥ=hᵤ,devlimit,
+    transpose && return panelize((v,u)->surface(u,v),v₀,v₁,u₀,u₁,;hᵤ,hᵥ,devlimit,
                                  transpose=false,flip=!flip,N_max,verbose,T,kwargs...)
 
     # Check inputs and get output type
@@ -37,7 +38,8 @@ function panelize(surface,u₀=0.,u₁=1.,v₀=0.,v₁=1.;hᵤ=1.,hᵥ=hᵤ,devl
         u,du = 0.5uᵢ[i+1]+0.5uᵢ[i],uᵢ[i+1]-uᵢ[i] # Parametric center & width
 
         # Find equidistant points along strip center
-        S,s⁻¹ = arclength(v->surface(u,v),hᵥ,devlimit,v₀,v₁)
+        uv₀,uv₁ = submerge ? newlimits(v->surface(u,v),v₀,v₁) : (v₀,v₁)
+        S,s⁻¹ = arclength(v->surface(u,v),hᵥ,devlimit,uv₀,uv₁)
         verbose && @show i,S
         S ≤ 0.5hᵥ && return init               # not enough height
         ve = s⁻¹(range(0,S,round(Int,S/hᵥ)+1)) # panel endpoints
@@ -53,7 +55,11 @@ function panelize(surface,u₀=0.,u₁=1.,v₀=0.,v₁=1.;hᵤ=1.,hᵥ=hᵤ,devl
     length(panels) ≤ N_max && return Table(panels)
     throw(ArgumentError("length(panels)=$(length(panels))>$N_max. Increase hᵤ,hᵥ,devlimit and/or N_max."))
 end
-
+function newlimits(r,low,high)
+    z(u) = r(u)[3]+√eps(u)
+    zero = find_zero(z,(low,high))
+    return z(low)<0 ? (low,zero) : (zero,high)
+end
 using QuadGK,DataInterpolations
 """
     arclength(r, Δs, devlimit, low, high) -> S, u(s)

test/runtests.jl

@@ -68,9 +68,13 @@ end
     area_checks(panels.dA,0.18)
     @test sum(panels.dA) ≈ 4π^2*0.3
 
-    # Check sign is correct when flipped
-    panels_bad = panelize(torus,0,2pi,0,2pi,hᵤ=0.6,hᵥ=0.3,devlimit=Inf)
-    @test panels[1].n ⋅ panels_bad[1].n > 0.99
+    # Check sign is flipped
+    panels = panelize(sphere,0,pi,0,2pi,hᵤ=0.6,flip=true)
+    @test all(p->p.x'p.n<-0.9,panels)
+
+    # Check submerged
+    panels = panelize(spheroid,0,2pi/3,0,2pi,hᵤ=0.6,hᵥ=0.3,transpose=true,submerge=true)
+    @test all(x->x[3]<0,panels.verts)
 
     # Check inputs
     @test_throws ArgumentError panelize(spheroid,0,pi,0,2pi,hᵤ=0)