Implement local version of density interpolation

Project.toml

@@ -6,10 +6,14 @@ version = "0.1.0"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 ElementaryPDESolutions = "88a69b33-da0f-4502-8c8c-d680cf4d883b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"
+Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

ext/IntiMeshesExt.jl

@@ -75,10 +75,10 @@ function Meshes.viz!(el::Inti.ReferenceInterpolant, args...; kwargs...)
 end
 
 function Meshes.viz(els::AbstractVector{<:Inti.ReferenceInterpolant}, args...; kwargs...)
-    return viz([to_meshes(el) for el in els])
+    return viz([to_meshes(el) for el in els], args...; kwargs...)
 end
 function Meshes.viz!(els::AbstractVector{<:Inti.ReferenceInterpolant}, args...; kwargs...)
-    return viz!([to_meshes(el) for el in els])
+    return viz!([to_meshes(el) for el in els], args...; kwargs...)
 end
 
 function Meshes.viz(msh::Inti.AbstractMesh, args...; kwargs...)

src/Inti.jl

@@ -14,6 +14,7 @@ using LinearMaps
 using NearestNeighbors
 using SparseArrays
 using StaticArrays
+using QuadGK
 using SpecialFunctions
 using Printf
 

src/adaptive.jl

@@ -42,7 +42,7 @@ function adaptive_correction(iop::IntegralOperator; tol, maxdist = nothing, maxs
     else
         maxdist
     end
-    dict_near = near_interaction_list(X, Y; tol = maxdist)
+    dict_near = elements_to_near_targets(X, Y; tol = maxdist)
     T = eltype(iop)
     msh = mesh(Y)
     correction = (I = Int[], J = Int[], V = T[])
@@ -166,14 +166,14 @@ function adaptive_integration_singular(f, τ̂::ReferenceLine, x̂ₛ; kwargs...
     end
 end
 
-function near_interaction_list(QX::Quadrature, QY::Quadrature; tol)
+function elements_to_near_targets(QX::Quadrature, QY::Quadrature; tol)
     X = [coords(q) for q in QX]
     msh = mesh(QY)
-    return near_interaction_list(X, msh; tol)
+    return elements_to_near_targets(X, msh; tol)
 end
-function near_interaction_list(X, QY::Quadrature; tol)
+function elements_to_near_targets(X, QY::Quadrature; tol)
     msh = mesh(QY)
-    return near_interaction_list(X, msh; tol)
+    return elements_to_near_targets(X, msh; tol)
 end
 
 """

src/bdim.jl

@@ -177,3 +177,152 @@ function bdim_correction(
     δD = sparse(Is, Js, Ds, num_trgs, n)
     return δS, δD
 end
+
+function local_bdim_correction(
+    pde,
+    target,
+    source::Quadrature;
+    green_multiplier::Vector{<:Real},
+    parameters = DimParameters(),
+    derivative::Bool = false,
+    maxdist = Inf,
+    kneighbor = 1,
+)
+    imat_cond = imat_norm = res_norm = rhs_norm = -Inf
+    T = default_kernel_eltype(pde) # Float64
+    N = ambient_dimension(source)
+    m, n = length(target), length(source)
+    msh = source.mesh
+    qnodes = source.qnodes
+    neighbors = topological_neighbors(msh, kneighbor)
+    dict_near = etype_to_nearest_points(target, source; maxdist)
+    # find first an appropriate set of source points to center the monopoles
+    qmax = sum(size(mat, 1) for mat in values(source.etype2qtags)) # max number of qnodes per el
+    ns   = ceil(Int, parameters.sources_oversample_factor * qmax)
+    # compute a bounding box for source points
+    low_corner = reduce((p, q) -> min.(coords(p), coords(q)), source)
+    high_corner = reduce((p, q) -> max.(coords(p), coords(q)), source)
+    xc = (low_corner + high_corner) / 2
+    R = parameters.sources_radius_multiplier * norm(high_corner - low_corner) / 2
+    xs = if N === 2
+        uniform_points_circle(ns, R, xc)
+    elseif N === 3
+        fibonnaci_points_sphere(ns, R, xc)
+    else
+        error("only 2D and 3D supported")
+    end
+    # figure out if we are dealing with a scalar or vector PDE
+    σ = if T <: Number
+        1
+    else
+        @assert allequal(size(T))
+        size(T, 1)
+    end
+    # compute traces of monopoles on the source mesh
+    G    = SingleLayerKernel(pde, T)
+    γ₁G  = DoubleLayerKernel(pde, T)
+    γ₁ₓG = AdjointDoubleLayerKernel(pde, T)
+    γ₀B  = Matrix{T}(undef, length(source), ns)
+    γ₁B  = Matrix{T}(undef, length(source), ns)
+    for k in 1:ns
+        for j in 1:length(source)
+            γ₀B[j, k] = G(source[j], xs[k])
+            γ₁B[j, k] = γ₁ₓG(source[j], xs[k])
+        end
+    end
+    Is, Js, Ss, Ds = Int[], Int[], T[], T[]
+    for (E, qtags) in source.etype2qtags
+        els = elements(msh, E)
+        near_list = dict_near[E]
+        nq, ne = size(qtags)
+        @assert length(near_list) == ne
+        # preallocate a local matrix to store interpolant values resulting
+        # weights. To benefit from Lapack, we must convert everything to
+        # matrices of scalars, so when `T` is an `SMatrix` we are careful to
+        # convert between the `Matrix{<:SMatrix}` and `Matrix{<:Number}` formats
+        # by viewing the elements of type `T` as `σ × σ` matrices of
+        # `eltype(T)`.
+        M_  = Matrix{eltype(T)}(undef, 2 * nq * σ, ns * σ)
+        W_  = Matrix{eltype(T)}(undef, 2 * nq * σ, σ)
+        W   = T <: Number ? W_ : Matrix{T}(undef, 2 * nq, 1)
+        Θi_ = Matrix{eltype(T)}(undef, σ, ns * σ)
+        Θi  = T <: Number ? Θi_ : Matrix{T}(undef, 1, ns)
+        K   = derivative ? γ₁ₓG : G
+        # for each element, we will solve Mᵀ W = Θiᵀ, where W is a vector of
+        # size 2nq, and Θi is a row vector of length(ns)
+        for n in 1:ne
+            # if there is nothing near, skip immediately to next element
+            isempty(near_list[n]) && continue
+            el = els[n]
+            # copy the monopoles/dipoles for the current element
+            jglob = @view qtags[:, n]
+            M0 = @view γ₀B[jglob, :]
+            M1 = @view γ₁B[jglob, :]
+            _copyto!(view(M_, 1:(nq*σ), :), M0)
+            _copyto!(view(M_, (nq*σ+1):2*nq*σ, :), M1)
+            F_ = qr!(transpose(M_))
+            @debug (imat_cond = max(cond(M_), imat_cond)) maxlog = 0
+            @debug (imat_norm = max(norm(M_), imat_norm)) maxlog = 0
+            # quadrature for auxiliary surface. In global dim, this is the same
+            # as the source quadrature, and independent of element. In local
+            # dim, this is constructed for each element using its neighbors.
+            function translate(q::QuadratureNode, x, s)
+                return QuadratureNode(coords(q) + x, weight(q), s * normal(q))
+            end
+            nei = neighbors[(E, n)]
+            qtags_nei = Int[]
+            for (E, m) in nei
+                append!(qtags_nei, source.etype2qtags[E][:, m])
+            end
+            qnodes_nei = source.qnodes[qtags_nei]
+            jac        = jacobian(el, 0.5)
+            ν          = -_normal(jac)
+            h          = sum(qnodes[i].weight for i in jglob)*4
+            qnodes_op  = map(q -> translate(q, h * ν, -1), qnodes_nei)
+            bindx      = boundary1d(nei, msh)
+            l, r       = nodes(msh)[-bindx[1]], nodes(msh)[bindx[2]]
+            Q, W       = gauss(4nq, 0, h)
+            qnodes_l   = [QuadratureNode(l.+q.*ν, w, SVector(-ν[2], ν[1])) for (q, w) in zip(Q, W)]
+            qnodes_r   = [QuadratureNode(r.+q.*ν, w, SVector(ν[2], -ν[1])) for (q, w) in zip(Q, W)]
+            qnodes_aux = append!(qnodes_nei, qnodes_op, qnodes_l, qnodes_r)
+            # qnodes_aux = source.qnodes # this is the global dim
+            for i in near_list[n]
+                # integrate the monopoles/dipoles over the auxiliary surface
+                # with target x: Θₖ <-- S[γ₁Bₖ](x) - D[γ₀Bₖ](x) + μ * Bₖ(x)
+                x = target[i]
+                μ = green_multiplier[i] # - 1/2
+                for k in 1:ns
+                    Θi[k] = μ * K(x, xs[k])
+                end
+                for q in qnodes_aux
+                    SK = G(x, q)
+                    DK = γ₁G(x, q)
+                    for k in 1:ns
+                        Θi[k] += (SK * γ₁ₓG(q, xs[k]) - DK * G(q, xs[k])) * weight(q)
+                    end
+                end
+                Θi_ = _copyto!(Θi_, Θi)
+                @debug (rhs_norm = max(rhs_norm, norm(Θi))) maxlog = 0
+                W_ = ldiv!(W_, F_, transpose(Θi_))
+                @debug (res_norm = max(norm(Matrix(F_) * W_ - transpose(Θi_)), res_norm)) maxlog =
+                    0
+                W = T <: Number ? W_ : _copyto!(W, W_)
+                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
+    @debug """Condition properties of bdim correction:
+    |-- max interp. matrix cond.: $imat_cond
+    |-- max interp. matrix norm : $imat_norm
+    |-- max residual error:       $res_norm
+    |-- max norm of source term:  $rhs_norm
+    """
+    δS = sparse(Is, Js, Ss, m, n)
+    δD = sparse(Is, Js, Ds, m, n)
+    return δS, δD
+end

src/mesh.jl

@@ -51,11 +51,11 @@ for a given mesh).
 function elements end
 
 """
-    element_types(msh::AbstractMesh)
+    elements(msh::AbstractMesh,E::DataType)
 
-Return the element types present in the `msh`.
+Return an iterator for all elements of type `E` on a mesh `msh`.
 """
-function element_types end
+elements(msh::AbstractMesh, E::DataType) = ElementIterator{E,typeof(msh)}(msh)
 
 """
     struct LagrangeMesh{N,T} <: AbstractMesh{N,T}
@@ -515,7 +515,7 @@ function connectivity(msh::SubMesh, E::DataType)
 end
 
 """
-    near_interaction_list(X,Y::AbstractMesh; tol)
+    elements_to_near_targets(X,Y::AbstractMesh; tol)
 
 For each element `el` of type `E` in `Y`, return the indices of the points in
 `X` which are closer than `tol` to the `center` of `el`.
@@ -526,7 +526,7 @@ fifth element of type `E`.
 
 If `tol` is a `Dict`, then `tol[E]` is the tolerance for elements of type `E`.
 """
-function near_interaction_list(
+function elements_to_near_targets(
     X::AbstractVector{<:SVector{N}},
     Y::AbstractMesh{N};
     tol,
@@ -538,13 +538,95 @@ function near_interaction_list(
     for E in element_types(Y)
         els = elements(Y, E)
         tol_ = isa(tol, Number) ? tol : tol[E]
-        idxs = _near_interaction_list(balltree, els, tol_)
+        idxs = _elements_to_near_targets(balltree, els, tol_)
         dict[E] = idxs
     end
     return dict
 end
 
-@noinline function _near_interaction_list(balltree, els, tol)
+@noinline function _elements_to_near_targets(balltree, els, tol)
     centers = map(center, els)
     return inrange(balltree, centers, tol)
 end
+
+"""
+    target_to_near_elements(X::AbstractVector{<:SVector{N}}, Y::AbstractMesh{N};
+    tol)
+
+For each target `x` in `X`, return a vector of tuples `(E, i)` where `E` is the
+type of the element in `Y` and `i` is the index of the element in `Y` such that
+`x` is closer than `tol` to the center of the element.
+"""
+function target_to_near_elements(
+    X::AbstractVector{<:SVector{N}},
+    Y::AbstractMesh{N};
+    tol,
+) where {N}
+    @assert isa(tol, Number) || isa(tol, Dict) "tol must be a number or a dictionary mapping element types to numbers"
+    dict = Dict{Int,Vector{Tuple{DataType,Int}}}()
+    balltree = BallTree(X)
+    for E in element_types(Y)
+        els = elements(Y, E)
+        tol_ = isa(tol, Number) ? tol : tol[E]
+        idxs = _target_to_near_elements(balltree, els, tol_)
+        for (i, idx) in enumerate(idxs)
+            dict[i] = get!(dict, i, Vector{Tuple{DataType,Int}}())
+            for j in idx
+                push!(dict[i], (E, j))
+            end
+        end
+    end
+    return dict
+end
+
+"""
+    topological_neighbors(msh::LagrangeMesh, k=1)
+
+Return the `k` neighbors of each element in `msh`. The one-neighbors are the
+elements that share a common vertex with the element, `k` neighbors are the
+one-neighbors of the `k-1` neighbors.
+
+This function returns a dictionary where the key is an `(Eᵢ,i)` tuple denoting
+the element `i` of type `E` in the mesh, and the value is a set of tuples
+`(Eⱼ,j)` where `Eⱼ` is the type of the element and `j` its index.
+"""
+function topological_neighbors(msh::AbstractMesh, k = 1)
+    # dictionary mapping a node index to all elements containing it. Note
+    # that the elements are stored as a tuple (type, index)
+    T = Tuple{DataType,Int}
+    node2els = Dict{Int,Vector{T}}()
+    for E in element_types(msh)
+        mat = connectivity(msh, E)::Matrix{Int} # connectivity matrix
+        np, Nel = size(mat)
+        for n in 1:Nel
+            for i in 1:np
+                idx = mat[i, n]
+                els = get!(node2els, idx, Vector{T}())
+                push!(els, (E, n))
+            end
+        end
+    end
+    # now revert the map to get the neighbors
+    one_neighbors = Dict{T,Set{T}}()
+    for (_, els) in node2els
+        for el in els
+            nei = get!(one_neighbors, el, Set{T}())
+            for el′ in els
+                push!(nei, el′)
+            end
+        end
+    end
+    # Recursively compute the neighbors from the one-neighbors
+    k_neighbors = deepcopy(one_neighbors)
+    while k > 1
+        # update neighborhood of each element
+        for el in keys(one_neighbors)
+            knn = k_neighbors[el]
+            for el′ in copy(knn)
+                union!(knn, one_neighbors[el′])
+            end
+        end
+        k -= 1
+    end
+    return k_neighbors
+end