begin working towards local_bdim

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

@@ -181,22 +181,19 @@ end
 function local_bdim_correction(
     pde,
     target,
-    source::Quadrature,
-    Sop,
-    Dop;
+    source::Quadrature;
     green_multiplier::Vector{<:Real},
     parameters = DimParameters(),
     derivative::Bool = false,
     maxdist = Inf,
 )
     imat_cond = imat_norm = res_norm = rhs_norm = -Inf
-    T = eltype(Sop)
+    T = default_kernel_eltype(pde) # Float64
     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)
+    # 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
@@ -268,27 +265,27 @@ function local_bdim_correction(
             # 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()
+            # 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]
+                μ = green_multiplier[i] # - 1/2
                 for k in 1:ns
                     Θi[k] = μ * K(x, xs[k])
                 end

src/mesh.jl

@@ -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,23 +538,57 @@ 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::LagrangeMesh, k = 1)
     @assert k == 1
@@ -586,68 +620,3 @@ function topological_neighbors(msh::LagrangeMesh, k = 1)
     #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

test/ldim_test.jl

@@ -0,0 +1,69 @@
+using Test
+using LinearAlgebra
+using Inti
+using Random
+using Meshes
+using GLMakie
+
+include("test_utils.jl")
+Random.seed!(1)
+
+N = 2
+t = :interior
+pde = Inti.Laplace(; dim = N)
+Inti.clear_entities!()
+Ω, msh = gmsh_disk(; center = [0.0, 0.0], rx = 1.0, ry = 1.0, meshsize = 0.2)
+Γ = Inti.external_boundary(Ω)
+
+Γ_msh = msh[Γ]
+Ω_msh = msh[Ω]
+nei = Inti.topological_neighbors(Ω_msh) |> collect
+function viz_neighbors(i)
+    k, v = nei[i]
+    E, idx = k
+    el = Inti.elements(Ω_msh, E)[idx]
+    fig, ax, pl = viz(el; color = 0, showsegments = true)
+    for (E, i) in v
+        el = Inti.elements(Ω_msh, E)[i]
+        viz!(el; color = 1 / 2, showsegments = true, alpha = 0.2)
+    end
+    return display(fig)
+end
+
+##
+
+quad = Inti.Quadrature(view(msh, Γ); qorder = 3)
+σ = t == :interior ? 1 / 2 : -1 / 2
+xs = t == :interior ? ntuple(i -> 3, N) : ntuple(i -> 0.1, N)
+T = Inti.default_density_eltype(pde)
+c = rand(T)
+u = (qnode) -> Inti.SingleLayerKernel(pde)(qnode, xs) * c
+dudn = (qnode) -> Inti.AdjointDoubleLayerKernel(pde)(qnode, xs) * c
+γ₀u = map(u, quad)
+γ₁u = map(dudn, quad)
+γ₀u_norm = norm(norm.(γ₀u, Inf), Inf)
+γ₁u_norm = norm(norm.(γ₁u, Inf), Inf)
+# single and double layer
+G = Inti.SingleLayerKernel(pde)
+S = Inti.IntegralOperator(G, quad)
+Smat = Inti.assemble_matrix(S)
+dG = Inti.DoubleLayerKernel(pde)
+D = Inti.IntegralOperator(dG, quad)
+Dmat = Inti.assemble_matrix(D)
+e0 = norm(Smat * γ₁u - Dmat * γ₀u - σ * γ₀u, Inf) / γ₀u_norm
+
+green_multiplier = fill(-0.5, length(quad))
+# δS, δD = Inti.bdim_correction(pde, quad, quad, Smat, Dmat; green_multiplier)
+δS, δD = Inti.local_bdim_correction(pde, quad, quad; green_multiplier)
+Sdim = Smat + δS
+Ddim = Dmat + δD
+# Sdim, Ddim = Inti.single_double_layer(;
+#     pde,
+#     target      = quad,
+#     source      = quad,
+#     compression = (method = :none,),
+#     correction  = (method = :ldim,),
+# )
+e1 = norm(Sdim * γ₁u - Ddim * γ₀u - σ * γ₀u, Inf) / γ₀u_norm
+@show norm(e0, Inf)
+@show norm(e1, Inf)