integral equation solvers

Author: maltezfaria

Project.toml

@@ -4,10 +4,14 @@ version = "0.2.5"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 

src/IO/plotsIO.jl

@@ -82,6 +82,18 @@ end
     end
 end
 
+@recipe function f(Q::NystromMesh)
+    @series begin
+        mesh(Q)
+    end
+    @series begin
+        aspect_ratio --> :equal
+        label := "Quadrature nodes"
+        seriestype := :scatter
+        PlotPoints(), collect(qcoords(Q))
+    end
+end
+
 @recipe function f(msh::GenericMesh)
     return msh, domain(msh)
 end

src/IntegralEquations/IntegralEquations.jl

@@ -0,0 +1,11 @@
+#=
+    Functionality related to defining and solving integral equations, with a
+    strong focus on boundary integral equations
+=#
+
+include("nystrommesh.jl")
+include("kernels.jl")
+include("integralpotential.jl")
+include("integraloperator.jl")
+include("density.jl")
+include("dim.jl")

src/IntegralEquations/adaptive.jl

@@ -0,0 +1,105 @@
+#=
+Adaptive quadrature technique to compute a sparse correction to integral
+operators when the kernel is weakly singular. The technique consists of:
+    1. Identifying a priorily the entries of the integral operator which need
+    correction. Those are given by entries whose target point and source element
+    lie closer than a certain threshold
+    2. For each element in the source surface, and for each target point close
+    to that element, compute the matrix entries by performing an adaptive
+    integration
+=#
+
+"""
+    near_interaction_list_adaptive(X,Y::AbstractMesh; atol)
+
+For each element `el` of type `E` in `Y`, return the indices of the points in `X` which
+are closer than `atol` to the `center` of `el`.
+
+This function returns a dictionary where e.g. `Dict[E][5] --> Vector{Int}` gives
+the indices of points in `X` which are closer than atol to the center of the
+fifth element of type `E`.
+"""
+function near_interaction_list(X, M::AbstractMesh; atol)
+    dict = Dict{DataType,Vector{Vector{Int}}}()
+    for E in keys(M)
+        pts = [center(el) for el in M[E]]
+        kdtree = KDTree(pts)
+        push!(dict, E => inrange(kdtree, X, atol))
+    end
+    return dict
+end
+
+function adaptive_correction(iop::IntegralOperator)
+    X, Y = target_surface(iop), source_surface(iop)
+    # maximum element size
+    lmax = -Inf
+    dict_near = near_interaction_list(dofs(X), Y; atol=0.1)
+    return
+end
+
+function assemble_gk(iop; compress=Matrix)
+    X, Y = target_surface(iop), source_surface(iop)
+    @timeit_debug "assemble dense part" begin
+        out = compress(iop)
+    end
+    @timeit_debug "compute near interaction list" begin
+        dict_near = near_interaction_list(dofs(X), Y; atol=0.1)
+    end
+    @timeit_debug "compute singular correction" begin
+        correction = singular_weights_gk(iop, dict_near)
+    end
+    return axpy!(1, correction, out) # out <-- out + correction
+end
+
+function singular_weights_gk(iop::IntegralOperator, dict_near)
+    X, Y = target_surface(iop), source_surface(iop)
+    T = eltype(iop)
+    Is = Int[]
+    Js = Int[]
+    Vs = T[]
+    # loop over all integration elements by types
+    for E in keys(Y)
+        iter = ElementIterator(Y, E)
+        qreg = Y.etype2qrule[E]
+        lag_basis = lagrange_basis(qnodes(qreg))
+        @timeit_debug "singular weights" begin
+            _singular_weights_gk!(Is, Js, Vs, iop, iter, dict_near, lag_basis, qreg)
+        end
+    end
+    Sp = sparse(Is, Js, Vs, size(iop)...)
+    return Sp
+end
+
+function _singular_weights_gk!(Is, Js, Vs, iop, iter, dict_near, lag_basis, qreg)
+    yi = qnodes(qreg)
+    K = kernel(iop)
+    X = target_surface(iop)
+    Y = source_surface(iop)
+    E = eltype(iter)
+    list_near = dict_near[E]
+    # Over the lines below we loop over all element, find which target nodes
+    # need to be regularized, and compute the regularized weights on the source
+    # points for that target point.
+    for (n, el) in enumerate(iter)               # loop over elements
+        jglob = Y.elt2dof[E][:, n]
+        for (i, jloc) in list_near[n]            # loop over near targets
+            ys = yi[jloc]
+            xdof = dofs(X)[i]
+            for (m, l) in enumerate(lag_basis)   # loop over lagrange basis
+                val, _ = quadgk(0, ys[1], 1; atol=1e-14) do ŷ
+                    y = el(ŷ)
+                    jac = jacobian(el, ŷ)
+                    μ = integration_measure(jac)
+                    ydof = NystromDOF(y, -1.0, jac, -1, -1)
+                    return K(xdof, ydof) * μ * l(ŷ)
+                end
+                # push data to sparse matrix
+                push!(Is, i)
+                j = jglob[m]
+                push!(Js, j)
+                push!(Vs, val - iop[i, j])
+            end
+        end
+    end
+    return Is, Js, Vs
+end

src/IntegralEquations/density.jl

@@ -0,0 +1,31 @@
+"""
+    struct Density{T,S} <: AbstractVector{T}
+
+Discrete density with values `vals` on an `AbstractMesh`
+"""
+struct Density{V,S<:AbstractMesh} <: AbstractVector{V}
+    vals::Vector{V}
+    mesh::S
+end
+
+# AbstractArray interface
+Base.size(σ::Density) = size(vals(σ))
+Base.getindex(σ::Density, args...) = getindex(vals(σ), args...)
+Base.setindex!(σ::Density, args...) = setindex!(vals(σ), args...)
+Base.similar(σ::Density) = Density(similar(vals(σ)), mesh(σ))
+
+vals(σ::Density) = σ.vals
+mesh(σ::Density) = σ.mesh
+
+function Density(f::Function, X::NystromMesh)
+    vals = [f(dof) for dof in qnodes(X)]
+    return Density(vals, X)
+end
+
+function γ₀(f, X::NystromMesh)
+    return Density(x -> f(coords(x)), X)
+end
+
+function γ₁(f, X::NystromMesh)
+    return Density(dof -> f(coords(dof), normal(dof)), X)
+end

src/IntegralEquations/dim.jl

@@ -0,0 +1,152 @@
+Base.@kwdef struct DimParameters
+    sources_oversample_factor::Float64 = 3
+    sources_radius_multiplier::Float64 = 5
+end
+
+"""
+    dim_correction(pde,X,Y,S,D[;location,p,derivative])
+
+Given a `pde` and a (possibly innacurate) discretizations of its single and
+double-layer operators `S` and `D` with domain `Y` and range `X`, compute
+corrections `δS` and `δD` such that `S + δS` and `D + δD` are more accurate
+approximations of the underlying integral operators.
+
+The following optional keyword arguments are available:
+- `location`: whether the target `X` lies on, inside, or outside the domain `Y`.
+- `p::DimParameters`: parameters associated with the density interpolation
+  method
+- `derivative`: if true, compute the correction to the adjoint double-layer and
+  hypersingular operators instead. In this case, `S` and `D` should contain a
+  discretization of adjoint double-layer and hypersingular operators,
+  respectively.
+"""
+function dim_correction(pde, X, Y, S, D; location=:onsurface, p=DimParameters(),
+                        derivative=false)
+    T = eltype(S)
+    isvectorial = T <: SMatrix
+    msg = "unrecognized value for kw `location`: received $location.
+    Valid options are `:onsurface`, `:inside` and `:outside`."
+    σ = location === :onsurface ? -0.5 :
+        location === :inside ? -1 : location === :outside ? 0 : error(msg)
+    dict_near = nearest_point_to_element(X, Y)
+    # find first an appropriate set of source points to center the monopoles
+    qmax = sum(size(etype2qtags(Y, E), 1) for E in keys(Y))
+    ns = ceil(Int, p.sources_oversample_factor * qmax)
+    pts = qcoords(Y)
+    bbox = HyperRectangle(pts)
+    xc = center(bbox)
+    r = p.sources_radius_multiplier * radius(bbox)
+    if ambient_dimension(Y) == 2
+        xs = uniform_points_circle(ns, r, xc)
+    elseif ambient_dimension(Y) == 3
+        xs = fibonnaci_points_sphere(ns, r, xc)
+    else
+        notimplemented()
+    end
+    # compute traces of monopoles on the source mesh
+    Xnodes = qnodes(X)
+    Ynodes = qnodes(Y)
+    G = SingleLayerKernel(pde, T)
+    γ₁G = AdjointDoubleLayerKernel(pde, T)
+    γ₀B = Matrix{T}(undef, length(Ynodes), ns)
+    γ₁B = Matrix{T}(undef, length(Ynodes), ns)
+    for k in 1:ns
+        for j in 1:length(Ynodes)
+            γ₀B[j, k] = G(Ynodes[j], xs[k])
+            γ₁B[j, k] = γ₁G(Ynodes[j], xs[k])
+        end
+    end
+    # integrate the monopoles/dipoles over Y with target on X. This is the
+    # slowest step, and passing a custom S,D can accelerate this computation
+    Θ = S * γ₁B - D * γ₀B
+    for k in 1:ns
+        for i in 1:length(Xnodes)
+            if derivative
+                Θ[i, k] += σ * γ₁G(Xnodes[i], xs[k])
+            else
+                Θ[i, k] += σ * G(Xnodes[i], xs[k])
+            end
+        end
+    end
+    # finally compute the corrected weights as sparse matrices
+    Is, Js, Ss, Ds = Int[], Int[], T[], T[]
+    dict_near = nearest_point_to_element(X, Y)
+    for E in keys(Y)
+        qtags = etype2qtags(Y, E)
+        near_list = dict_near[E]
+        nq, ne = size(qtags)
+        @assert length(near_list) == ne
+        M = Matrix{T}(undef, 2nq, ns)
+        for n in 1:ne # loop over elements of type E
+            # if there is nothing near, skip immediately to next element
+            isempty(near_list[n]) && continue
+            # the weights for points in near_list[n] must be corrected when
+            # integrating over the current element
+            jglob = @view qtags[:, n]
+            M0 = @view γ₀B[jglob, :]
+            M1 = @view γ₁B[jglob, :]
+            copy!(view(M, 1:nq, :), M0)
+            copy!(view(M, (nq + 1):(2nq), :), M1)
+            F = qr!(blockmatrix_to_matrix(M))
+            for i in near_list[n]
+                Θi = @view Θ[i:i, :]
+                W_ = (blockmatrix_to_matrix(Θi) / F.R) * adjoint(F.Q)
+                W = matrix_to_blockmatrix(W_, T)
+                for k in 1:nq
+                    push!(Is, i)
+                    push!(Js, jglob[k])
+                    push!(Ss, -W[nq + k]) # single layer corresponds to α=0,β=-1
+                    push!(Ds, W[k]) # double layer corresponds to α=1,β=0
+                end
+            end
+        end
+    end
+    δS = sparse(Is, Js, Ss)
+    δD = sparse(Is, Js, Ds)
+    return δS, δD
+end
+
+"""
+    nearest_point_to_element(X::NystroMesh,Y::NystromMesh; tol)
+
+For each element `el`, return a list with the indices of all points in `X` for
+which `el` is the nearest element. Ignore indices for which the distance exceeds `tol`.
+"""
+function nearest_point_to_element(X, Y::NystromMesh; tol=Inf)
+    y = collect(qcoords(Y))
+    kdtree = KDTree(y)
+    dict = Dict(j => Int[] for j in 1:length(y))
+    for i in eachindex(X)
+        qtag, d = nn(kdtree, X[i])
+        d > tol || push!(dict[qtag], i)
+    end
+    # dict[j] now contains indices in X for which the j quadrature node in Y is
+    # the closest. Next we reverse the map
+    etype2nearlist = Dict{DataType,Vector{Vector{Int}}}()
+    for E in keys(Y)
+        tags = etype2qtags(Y, E)
+        nq, ne = size(tags)
+        etype2nearlist[E] = nearlist = [Int[] for _ in 1:ne]
+        for j in 1:ne # loop over each element of type E
+            for q in 1:nq # loop over qnodes in the element
+                qtag = tags[q, j]
+                append!(nearlist[j], dict[qtag])
+            end
+        end
+    end
+    return etype2nearlist
+end
+function nearest_point_to_element(X::NystromMesh, Y::NystromMesh)
+    if X === Y
+        # when both surfaces are the same, the "near points" of an element are
+        # simply its own quadrature points
+        dict = Dict{DataType,Vector{Vector{Int}}}()
+        for E in keys(Y)
+            idx_dofs = etype2qtags(Y, E)
+            dict[E] = map(i -> collect(i), eachcol(idx_dofs))
+        end
+    else
+        dict = nearest_point_to_element(collect(qcoords(X)), Y)
+    end
+    return dict
+end

src/IntegralEquations/integraloperator.jl

@@ -0,0 +1,48 @@
+"""
+    struct IntegralOperator{T} <: AbstractMatrix{T}
+
+A discrete linear integral operator given by
+```math
+I[u](x) = \\int_{\\Gamma\\_s} K(x,y)u(y) ds_y, x \\in \\Gamma_{t}
+```
+where ``\\Gamma_s`` and ``\\Gamma_t`` are the source and target domains, respectively.
+"""
+struct IntegralOperator{T} <: AbstractMatrix{T}
+    kernel::AbstractKernel
+    target_mesh::AbstractMesh
+    source_mesh::AbstractMesh
+end
+
+kernel(iop::IntegralOperator) = iop.kernel
+target_mesh(iop::IntegralOperator) = iop.target_mesh
+source_mesh(iop::IntegralOperator) = iop.source_mesh
+
+function IntegralOperator(k, X, Y=X)
+    T = return_type(k)
+    msg = """IntegralOperator of nonbits being created"""
+    isbitstype(T) || (@warn msg)
+    return IntegralOperator{T}(k, X, Y)
+end
+
+function Base.size(iop::IntegralOperator)
+    X = target_mesh(iop)
+    Y = source_mesh(iop)
+    return (length(qnodes(X)), length(qnodes(Y)))
+end
+
+function Base.getindex(iop::IntegralOperator, i::Integer, j::Integer)
+    k = kernel(iop)
+    targets = qnodes(target_mesh(iop))
+    sources = qnodes(source_mesh(iop))
+    return k(targets[i], sources[j]) * weight(sources[j])
+end
+
+# convenience constructors
+SingleLayerOperator(op::AbstractPDE, X, Y=X) = IntegralOperator(SingleLayerKernel(op), X, Y)
+DoubleLayerOperator(op::AbstractPDE, X, Y=X) = IntegralOperator(DoubleLayerKernel(op), X, Y)
+function AdjointDoubleLayerOperator(op::AbstractPDE, X, Y=X)
+    return IntegralOperator(AdjointDoubleLayerKernel(op), X, Y)
+end
+function HyperSingularOperator(op::AbstractPDE, X, Y=X)
+    return IntegralOperator(HyperSingularKernel(op), X, Y)
+end