WIP on local bdim

Author: maltezfaria

src/bdim.jl

@@ -177,3 +177,150 @@ function bdim_correction(
     δD = sparse(Is, Js, Ds, num_trgs, n)
     return δS, δD
 end
+
+function local_bdim_correction(
+    pde,
+    target,
+    source::Quadrature,
+    Sop,
+    Dop;
+    green_multiplier::Vector{<:Real},
+    parameters = DimParameters(),
+    derivative::Bool = false,
+    maxdist = Inf,
+)
+    imat_cond = imat_norm = res_norm = rhs_norm = -Inf
+    T = eltype(Sop)
+    N = ambient_dimension(source)
+    @assert eltype(Dop) == T "eltype of S and D must match"
+    m, n = length(target), length(source)
+    msh = source.mesh
+    qnodes = source.qnodes
+    neighbors = topological_neighbors(msh)
+    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, n) in nei
+                append!(qtags_nei, source.etype2qtags[E][:, n])
+            end
+            qnodes_nei = source.qnodes[qtags_nei]
+            jac        = jacobian(el, 0.5)
+            ν          = _normal(jac)
+            h          = sum(qnodes[i].weight for i in jglob)
+            qnodes_op  = map(q -> translate(q, h * ν, -1), qnodes_nei)
+            a, b       = external_boundary()
+            # qnodes_aux = source.qnodes[jglob]
+            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]
+                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}
@@ -548,3 +548,106 @@ end
     centers = map(center, els)
     return inrange(balltree, centers, tol)
 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.
+"""
+function topological_neighbors(msh::LagrangeMesh, k = 1)
+    @assert 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 = msh.etype2mat[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
+    #TODO: for k > 1, recursively compute the neighbors from the one-neighbors
+    return one_neighbors
+end
+
+"""
+    element_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`.
+
+This function returns a dictionary where e.g. `dict[E][5] --> Vector{Int}` gives
+the indices of points in `X` which are closer than `tol` to the center of the
+fifth element of type `E`.
+
+If `tol` is a `Dict`, then `tol[E]` is the tolerance for elements of type `E`.
+"""
+function element_to_near_targets(
+    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"
+    # for each element type, build the list of targets close to a given element
+    dict = Dict{DataType,Vector{Vector{Int}}}()
+    balltree = BallTree(X)
+    for E in element_types(Y)
+        els = elements(Y, E)
+        tol_ = isa(tol, Number) ? tol : tol[E]
+        idxs = _element_to_near_targets(balltree, els, tol_)
+        dict[E] = idxs
+    end
+    return dict
+end
+
+@noinline function _element_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

src/quadrature.jl

@@ -45,6 +45,8 @@ flip_normal(q::QuadratureNode) = QuadratureNode(q.coords, q.weight, -q.normal)
 
 weight(q::QuadratureNode) = q.weight
 
+translate(q::QuadratureNode, x) = QuadratureNode(coords(q) + x, weight(q), normal(q))
+
 # useful for using either a quadrature node or a just a simple point in
 # `IntegralOperators`.
 coords(x::Union{SVector,Tuple}) = SVector(x)