Merge pull request #70 from IntegralEquations/stokesfmm

Stokes FMM3D support
IntegralEquations · Jun 25, 2024 · 399d4ab · 399d4ab
2 parents ec7f3af + cf83013
commit 399d4ab
Show file tree

Hide file tree

Showing 2 changed files with 37 additions and 2 deletions.
diff --git a/ext/IntiFMM3DExt.jl b/ext/IntiFMM3DExt.jl
@@ -3,6 +3,7 @@ module IntiFMM3DExt
 import Inti
 import FMM3D
 import LinearMaps
+using StaticArrays # For Stokes types
 
 function __init__()
     @info "Loading Inti.jl FMM3D extension"
@@ -131,6 +132,33 @@ function Inti._assemble_fmm3d(iop::Inti.IntegralOperator; rtol = sqrt(eps()))
             out = FMM3D.hfmm3d(rtol, zk, sources; dipvecs, targets, pgt = 2)
             return copyto!(y, sum(xnormals .* out.gradtarg; dims = 1) |> vec)
         end
+        # Stokes
+    elseif K isa Inti.SingleLayerKernel{SMatrix{3,3,Float64,9},<:Inti.Stokes{3,Float64}}
+        T = SVector{3,Float64}
+        stoklet = Matrix{Float64}(undef, 3, n)
+        return LinearMaps.LinearMap{SVector{3,Float64}}(m, n) do y, x
+            # multiply by weights and constant
+            stoklet[:] = 1 / (4 * π * K.pde.μ) .* reinterpret(Float64, weights .* x)
+            out = FMM3D.stfmm3d(rtol, sources; stoklet, targets, ppregt = 1)
+            return copyto!(y, reinterpret(T, out.pottarg))
+        end
+    elseif K isa Inti.DoubleLayerKernel{SMatrix{3,3,Float64,9},<:Inti.Stokes{3,Float64}}
+        T = SVector{3,Float64}
+        normals = Matrix{Float64}(undef, 3, n)
+        for j in 1:n
+            normals[:, j] = Inti.normal(iop.source[j])
+        end
+        strsvec = similar(normals, Float64)
+        strslet = similar(normals, Float64)
+        # multiply by weights and constant
+        for j in 1:n
+            strsvec[:, j] = -1 / (4 * π) * view(normals, :, j) .* weights[j]
+        end
+        return LinearMaps.LinearMap{SVector{3,Float64}}(m, n) do y, x
+            strslet[:] = reinterpret(Float64, x)
+            out = FMM3D.stfmm3d(rtol, sources; strslet, strsvec, targets, ppregt = 1)
+            return copyto!(y, reinterpret(T, out.pottarg))
+        end
     else
         error("integral operator not supported by Inti's FMM3D wrapper")
     end

diff --git a/test/fmm3d_test.jl b/test/fmm3d_test.jl
@@ -16,18 +16,25 @@ include("test_utils.jl")
 Γ₂ = Inti.external_boundary(Ω₂)
 Γ₂_quad = Inti.Quadrature(view(msh₂, Γ₂); qorder = 3)
 
-for pde in (Inti.Laplace(; dim = 3), Inti.Helmholtz(; dim = 3, k = 1.2))
+for pde in (
+    Inti.Laplace(; dim = 3),
+    Inti.Helmholtz(; dim = 3, k = 1.2),
+    Inti.Stokes(; dim = 3, μ = 0.5),
+)
     @testset "PDE: $pde" begin
         for K in (
             Inti.DoubleLayerKernel(pde),
             Inti.SingleLayerKernel(pde),
             Inti.AdjointDoubleLayerKernel(pde),
             Inti.HyperSingularKernel(pde),
         )
+            # TODO Stokes has only single and double layer implemented for now
+            (K isa Inti.AdjointDoubleLayerKernel && pde isa Inti.Stokes) && continue
+            (K isa Inti.HyperSingularKernel && pde isa Inti.Stokes) && continue
             for Γ_quad in (Γ₁_quad, Γ₂_quad)
                 iop = Inti.IntegralOperator(K, Γ₁_quad, Γ_quad)
                 iop_fmm = Inti.assemble_fmm(iop; rtol = 1e-8)
-                x = rand(eltype(iop), size(iop, 2))
+                x = rand(Inti.default_density_eltype(pde), size(iop, 2))
                 yapprox = iop_fmm * x
                 # test on a given index set
                 idx_test = rand(1:size(iop, 1), 10)