Extend vdim to support vector-valued kernels

.gitignore

@@ -1,19 +1,5 @@
-*.jl.*.cov
-*.jl.cov
-*.jl.mem
-*.msh
-*.pos
-/Manifest.toml
+Manifest.toml
 /docs/Manifest.toml
-/docs/build/
 /test/Manifest.toml
-.vscode
-/docs/src/pluto-examples/*.md
-/sandbox/*
-/talks/*
+/docs/build/
 /test/vtk_output/*
-
-.markdownlint.json
-Inti.code-workspace
-docs/src/examples/heat.gif
-docs/src/.DS_Store

Project.toml

@@ -10,6 +10,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 HAdaptiveIntegration = "eaa5ad34-b243-48e9-b09c-54bc0655cecf"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
@@ -46,16 +47,19 @@ Bessels = "0.2"
 DataStructures = "0.18 - 0.19"
 ElementaryPDESolutions = "0.2 - 0.3"
 FMM2D = "0.2"
-FMM3D = "1"
+# FMM3D <2.1.0 has a different scaling convention for the kernels, not
+# including the 1/(4pi) prefactor
+FMM3D = "2.1.0"
 FMMLIB2D = "0.3"
 ForwardDiff = "0.10, 1"
 Gmsh = "0.3"
 HAdaptiveIntegration = "0.2"
 HMatrices = "0.2"
 LinearAlgebra = "1"
 LinearMaps = "3"
+Logging = "1"
 Makie = "0.21 - 0.24"
-Meshes = "0.42 - 0.55"
+Meshes = "0.42 - 0.56"
 NearestNeighbors = "0.4"
 OrderedCollections = "1"
 Pkg = "1"

ext/IntiFMM2DExt.jl

@@ -3,6 +3,7 @@ module IntiFMM2DExt
 import Inti
 import FMM2D
 import LinearMaps
+using StaticArrays
 
 function __init__()
     return @info "Loading Inti.jl FMM2D extension"
@@ -146,6 +147,148 @@ function Inti._assemble_fmm2d(iop::Inti.IntegralOperator; rtol = sqrt(eps()))
                 return copyto!(y, sum(xnormals .* out.gradtarg; dims = 1) |> vec)
             end
         end
+    elseif K isa Inti.GradientSingleLayerKernel{<:SVector{2}, <:Inti.Laplace{2}}
+        charges = Vector{Float64}(undef, n)
+        return LinearMaps.LinearMap{SVector{2, Float64}}(m, n) do y, x
+            # multiply by weights and constant
+            @. charges = -1 / (2 * π) * weights * x
+            if same_surface
+                out =
+                    FMM2D.rfmm2d(; sources = sources, charges = charges, eps = rtol, pg = 2)
+                return copyto!(y, reinterpret(SVector{2, Float64}, vec(out.grad)))
+            else
+                out = FMM2D.rfmm2d(;
+                    sources = sources,
+                    charges = charges,
+                    targets = targets,
+                    eps = rtol,
+                    pgt = 2,
+                )
+                return copyto!(y, reinterpret(SVector{2, Float64}, vec(out.gradtarg)))
+            end
+        end
+    elseif K isa Inti.GradientDoubleLayerKernel{<:SVector{2}, <:Inti.Laplace{2}}
+        normals = Matrix{Float64}(undef, 2, n)
+        for j in 1:n
+            normals[:, j] = Inti.normal(iop.source[j])
+        end
+        dipvecs = similar(normals)
+        return LinearMaps.LinearMap{SVector{2, Float64}}(m, n) do y, x
+            # multiply by weights and constant
+            for j in 1:n
+                dipvecs[:, j] = -1 / (2 * π) * view(normals, :, j) * weights[j]
+            end
+            dipstr = Vector{Float64}(undef, n)
+            for j in 1:n
+                dipstr[j] = x[j]
+            end
+            if same_surface
+                out = FMM2D.rfmm2d(;
+                    sources = sources,
+                    dipstr = dipstr,
+                    dipvecs = dipvecs,
+                    eps = rtol,
+                    pg = 2,
+                )
+                return copyto!(y, reinterpret(SVector{2, Float64}, vec(out.grad)))
+            else
+                out = FMM2D.rfmm2d(;
+                    sources = sources,
+                    targets = targets,
+                    dipstr = dipstr,
+                    dipvecs = dipvecs,
+                    eps = rtol,
+                    pgt = 2,
+                )
+                return copyto!(y, reinterpret(SVector{2, Float64}, vec(out.gradtarg)))
+            end
+        end
+    elseif K isa Inti.SourceGradientSingleLayerKernel{<:Any, <:Inti.Laplace{2}}
+        # ∇yG(x,y)⋅g : dipoles with vector strengths g (scalar output). The W operator
+        # W[g] = -∫∇yG⋅g applies the leading minus when this map is assembled.
+        dipvecs = Matrix{Float64}(undef, 2, n)
+        dipstr = ones(Float64, n)
+        return LinearMaps.LinearMap{Float64}(m, n) do y, x
+            for j in 1:n
+                dipvecs[:, j] = -(1 / (2 * π)) * x[j] * weights[j]
+            end
+            if same_surface
+                out = FMM2D.rfmm2d(;
+                    sources = sources,
+                    dipstr = dipstr,
+                    dipvecs = dipvecs,
+                    eps = rtol,
+                    pg = 1,
+                )
+                return copyto!(y, out.pot)
+            else
+                out = FMM2D.rfmm2d(;
+                    sources = sources,
+                    targets = targets,
+                    dipstr = dipstr,
+                    dipvecs = dipvecs,
+                    eps = rtol,
+                    pgt = 1,
+                )
+                return copyto!(y, out.pottarg)
+            end
+        end
+    elseif K isa Inti.HessianKernel{<:Any, <:Inti.Laplace{2}}
+        # X forward (PV part) = +∫∇ₓ∇ₓG⋅g = -∇ₓ(∫∇yG⋅g). The bracket is the
+        # `SourceGradientSingleLayer` dipole field above; its target-gradient (`gradtarg`,
+        # pgt = 2) is ∫∇ₓ∇yG⋅g = -X_forward. Negating the dipole strengths relative to the
+        # W branch folds in the leading minus, so `gradtarg` is the `SVector` output.
+        if K.charge_dipole == :charge
+            charges = Vector{Float64}(undef, n)
+            return LinearMaps.LinearMap{SMatrix{2, 2, Float64, 4}}(m, n) do y, x
+                @. charges = -1 / (2 * π) * weights * x
+                if same_surface
+                    out = FMM2D.rfmm2d(; sources = sources, charges = charges, eps = rtol, pg = 3)
+                    H = out.hess        # (3, n): ∂xx, ∂xy, ∂yy
+                else
+                    out = FMM2D.rfmm2d(;
+                        sources = sources,
+                        charges = charges,
+                        targets = targets,
+                        eps = rtol,
+                        pgt = 3,
+                    )
+                    H = out.hesstarg
+                end
+                @inbounds for i in 1:m
+                    y[i] = SMatrix{2, 2, Float64, 4}(H[1, i], H[2, i], H[2, i], H[3, i])
+                end
+                return y
+            end
+        else
+            dipvecs = Matrix{Float64}(undef, 2, n)
+            dipstr = ones(Float64, n)
+            return LinearMaps.LinearMap{SVector{2, Float64}}(m, n) do y, x
+                for j in 1:n
+                    dipvecs[:, j] = (1 / (2 * π)) * x[j] * weights[j]
+                end
+                if same_surface
+                    out = FMM2D.rfmm2d(;
+                        sources = sources,
+                        dipstr = dipstr,
+                        dipvecs = dipvecs,
+                        eps = rtol,
+                        pg = 2,
+                    )
+                    return copyto!(y, reinterpret(SVector{2, Float64}, vec(out.grad)))
+                else
+                    out = FMM2D.rfmm2d(;
+                        sources = sources,
+                        targets = targets,
+                        dipstr = dipstr,
+                        dipvecs = dipvecs,
+                        eps = rtol,
+                        pgt = 2,
+                    )
+                    return copyto!(y, reinterpret(SVector{2, Float64}, vec(out.gradtarg)))
+                end
+            end
+        end
         # Helmholtz
     elseif K isa Inti.SingleLayerKernel{ComplexF64, <:Inti.Helmholtz{2}}
         charges = Vector{ComplexF64}(undef, n)
@@ -287,6 +430,178 @@ function Inti._assemble_fmm2d(iop::Inti.IntegralOperator; rtol = sqrt(eps()))
                 return copyto!(y, sum(xnormals .* out.gradtarg; dims = 1) |> vec)
             end
         end
+    elseif K isa Inti.GradientSingleLayerKernel{<:SVector{2}, <:Inti.Helmholtz{2}}
+        charges = Vector{ComplexF64}(undef, n)
+        zk = ComplexF64(K.op.k)
+        return LinearMaps.LinearMap{SVector{2, ComplexF64}}(m, n) do y, x
+            # multiply by weights and constant
+            @. charges = weights * x
+            if same_surface
+                out = FMM2D.hfmm2d(;
+                    zk = zk,
+                    sources = sources,
+                    charges = charges,
+                    eps = rtol,
+                    pg = 2,
+                )
+                return copyto!(y, reinterpret(SVector{2, ComplexF64}, vec(out.grad)))
+            else
+                out = FMM2D.hfmm2d(;
+                    zk = zk,
+                    sources = sources,
+                    charges = charges,
+                    targets = targets,
+                    eps = rtol,
+                    pgt = 2,
+                )
+                return copyto!(y, reinterpret(SVector{2, ComplexF64}, vec(out.gradtarg)))
+            end
+        end
+    elseif K isa Inti.GradientDoubleLayerKernel{<:SVector{2}, <:Inti.Helmholtz{2}}
+        normals = Matrix{Float64}(undef, 2, n)
+        for j in 1:n
+            normals[:, j] = Inti.normal(iop.source[j])
+        end
+        dipvecs = similar(normals, Float64)
+        dipstrs = Vector{ComplexF64}(undef, n)
+        zk = ComplexF64(K.op.k)
+        return LinearMaps.LinearMap{SVector{2, ComplexF64}}(m, n) do y, x
+            # multiply by weights and constant
+            for j in 1:n
+                dipvecs[:, j] = view(normals, :, j) * weights[j]
+            end
+            for j in 1:n
+                dipstrs[j] = x[j]
+            end
+            if same_surface
+                out = FMM2D.hfmm2d(;
+                    zk = zk,
+                    sources = sources,
+                    dipstr = dipstrs,
+                    dipvecs = dipvecs,
+                    eps = rtol,
+                    pg = 2,
+                )
+                return copyto!(y, reinterpret(SVector{2, ComplexF64}, vec(out.grad)))
+            else
+                out = FMM2D.hfmm2d(;
+                    zk = zk,
+                    sources = sources,
+                    targets = targets,
+                    dipstr = dipstrs,
+                    dipvecs = dipvecs,
+                    eps = rtol,
+                    pgt = 2,
+                )
+                return copyto!(y, reinterpret(SVector{2, ComplexF64}, vec(out.gradtarg)))
+            end
+        end
+
+    elseif K isa Inti.SourceGradientSingleLayerKernel{<:Any, <:Inti.Helmholtz{2}}
+        # ∫∇yG⋅g = Σ wⱼ gⱼ⋅∇yG : dipoles with vector strengths g. The W operator
+        # W[g] = -∫∇yG⋅g applies the leading minus when this map is assembled.
+        # FMM2D's `hfmm2d` requires real dipole directions (with a complex
+        # strength), so the complex density is split into real/imaginary parts:
+        #   Re(g): dipvec=Re(gⱼ), dipstr=wⱼ ;  Im(g): dipvec=Im(gⱼ), dipstr=i wⱼ.
+        dipvecs = Matrix{Float64}(undef, 2, n)
+        dipstr = Vector{ComplexF64}(undef, n)
+        zk = ComplexF64(K.op.k)
+        return LinearMaps.LinearMap{ComplexF64}(m, n) do y, x
+            fill!(y, 0)
+            for (getpart, strconst) in ((real, one(ComplexF64)), (imag, 1im))
+                for j in 1:n
+                    dipvecs[:, j] = getpart.(x[j])
+                    dipstr[j] = strconst * weights[j]
+                end
+                if same_surface
+                    out = FMM2D.hfmm2d(;
+                        zk = zk,
+                        sources = sources,
+                        dipstr = dipstr,
+                        dipvecs = dipvecs,
+                        eps = rtol,
+                        pg = 1,
+                    )
+                    y .+= out.pot
+                else
+                    out = FMM2D.hfmm2d(;
+                        zk = zk,
+                        sources = sources,
+                        targets = targets,
+                        dipstr = dipstr,
+                        dipvecs = dipvecs,
+                        eps = rtol,
+                        pgt = 1,
+                    )
+                    y .+= out.pottarg
+                end
+            end
+            return y
+        end
+    elseif K isa Inti.HessianKernel{<:Any, <:Inti.Helmholtz{2}}
+        # Charge→Hessian realization of the Hessian single-layer *volume* operator used by the
+        # `X = ∇W` VDIM correction in 2D: a scalar density `ρ` maps to `∫∇ₓ∇ₓG(x,y)ρ(y)dy`
+        # (a 2×2 `SMatrix` per target), i.e. the Hessian of the single-layer potential of
+        # charges `ρ`. `rfmm2d` returns the 3 unique second derivatives per point as `(3,·)` in
+        # the order ∂xx, ∂xy, ∂yy;
+        if K.charge_dipole == :charge
+            #`hess`/`hesstarg` is (3,·) = ∂xx, ∂xy, ∂yy.
+            charges = Vector{ComplexF64}(undef, n)
+            zk = ComplexF64(K.op.k)
+            return LinearMaps.LinearMap{SMatrix{2, 2, ComplexF64, 4}}(m, n) do y, x
+                @. charges = weights * x
+                if same_surface
+                    out = FMM2D.hfmm2d(; zk = zk, sources = sources, charges = charges, eps = rtol, pg = 3)
+                    H = out.hess
+                else
+                    out = FMM2D.hfmm2d(;
+                        zk = zk,
+                        sources = sources,
+                        charges = charges,
+                        targets = targets,
+                        eps = rtol,
+                        pgt = 3,
+                    )
+                    H = out.hesstarg
+                end
+                @inbounds for i in 1:m
+                    y[i] = SMatrix{2, 2, ComplexF64, 4}(H[1, i], H[2, i], H[2, i], H[3, i])
+                end
+                return y
+            end
+        # X forward = +∫∇ₓ∇ₓG⋅g = -∇ₓ(∫∇yG⋅g). The bracket is the ∇yG-dipole field (W
+        # forward); its target-gradient is ∫∇ₓ∇yG⋅g = -X_forward, so the dipole strengths
+        # are negated and `grad`/`gradtarg` is the output. As for the W forward, hfmm2d
+        # needs real dipvecs with complex strengths, so the complex density is split into
+        # real/imag parts (dipstr = wⱼ resp. i·wⱼ).
+        else
+            dipvecs = Matrix{Float64}(undef, 2, n)
+            dipstr = Vector{ComplexF64}(undef, n)
+            zk = ComplexF64(K.op.k)
+            return LinearMaps.LinearMap{SVector{2, ComplexF64}}(m, n) do y, x
+                fill!(y, zero(SVector{2, ComplexF64}))
+                for (getpart, strconst) in ((real, one(ComplexF64)), (imag, 1im))
+                    for j in 1:n
+                        dipvecs[:, j] = -getpart.(x[j])
+                        dipstr[j] = strconst * weights[j]
+                    end
+                    if same_surface
+                        out = FMM2D.hfmm2d(;
+                            zk = zk, sources = sources, dipstr = dipstr, dipvecs = dipvecs,
+                            eps = rtol, pg = 2,
+                        )
+                        y .+= reinterpret(SVector{2, ComplexF64}, vec(out.grad))
+                    else
+                        out = FMM2D.hfmm2d(;
+                            zk = zk, sources = sources, targets = targets, dipstr = dipstr,
+                            dipvecs = dipvecs, eps = rtol, pgt = 2,
+                        )
+                        y .+= reinterpret(SVector{2, ComplexF64}, vec(out.gradtarg))
+                    end
+                end
+                return y
+            end
+        end
     else
         error("integral operator not supported by Inti's FMM2D wrapper")
     end