Implement local version of vdim

Project.toml

@@ -7,6 +7,7 @@ Bessels = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 ElementaryPDESolutions = "88a69b33-da0f-4502-8c8c-d680cf4d883b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"

ext/IntiMakieExt.jl

@@ -76,10 +76,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...)
@@ -89,4 +89,45 @@ function Meshes.viz!(msh::Inti.AbstractMesh, args...; kwargs...)
     return viz!(to_meshes(msh), args...; kwargs...)
 end
 
+function Inti.viz_elements(els, msh)
+    E = first(Inti.element_types(msh))
+    Els = [Inti.elements(msh, E)[i] for (E, i) in els]
+    fig, _, _ = viz(Els)
+    viz!(msh; color = 0, showsegments = true, alpha = 0.3)
+    return display(fig)
+end
+
+function Inti.viz_elements_bords(Ei, els, ell, bords, msh; quad = nothing)
+    E = first(Inti.element_types(msh))
+    #ell = collect(Ei[(E, 1)])[1]
+    el = Inti.elements(msh, ell[1])[ell[2]]
+    fig, _, _ = viz(msh; color = 0, showsegments = true, alpha = 0.3)
+    viz!(el; color = 0, showsegments = true, alpha = 0.5)
+    for (E, i) in els
+        el = Inti.elements(msh, E)[i]
+        viz!(el; showsegments = true, alpha = 0.7)
+    end
+    viz!(bords; color = 4, showsegments = false, segmentsize = 5, segmentcolor = 4)
+    # if !isnothing(quad)
+    #     xs = [qnode.coords[1] for qnode in quad.qnodes]
+    #     ys = [qnode.coords[2] for qnode in quad.qnodes]
+    #     zs = [qnode.coords[3] for qnode in quad.qnodes]
+    #     nx = [Inti.normal(qnode)[1] for qnode in quad.qnodes]
+    #     ny = [Inti.normal(qnode)[2] for qnode in quad.qnodes]
+    #     nz = [Inti.normal(qnode)[3] for qnode in quad.qnodes]
+    #     Mke.arrows!(
+    #         xs,
+    #         ys,
+    #         zs,
+    #         nx,
+    #         ny,
+    #         nz;
+    #         lengthscale = 0.05,
+    #         linewidth = 0.01,
+    #         arrowsize = 0.01,
+    #     )
+    # end
+    return display(fig)
+end
+
 end # module

src/adaptive.jl

@@ -151,7 +151,7 @@ end
     adaptive_integration_singular(f, τ̂, x̂ₛ; kwargs...)
 
 Similar to [`adaptive_integration`](@ref), but indicates that `f` has an
-isolated (integrable) singularity at `x̂ₛ ∈ x̂ₛ`.
+isolated (integrable) singularity at `x̂ₛ ∈ τ̂`.
 
 The integration is performed by splitting `τ̂` so that `x̂ₛ` is a fixed vertex,
 guaranteeing that `f` is never evaluated at `x̂ₛ`. Aditionally, a suitable

src/api.jl

@@ -335,6 +335,25 @@ function volume_potential(; op, target, source::Quadrature, compression, correct
             correction.maxdist,
             interpolation_order,
         )
+    elseif correction.method == :ldim
+        loc = target === source ? :inside : correction.target_location
+        μ = _green_multiplier(loc)
+        green_multiplier = fill(μ, length(target))
+        shift = Val(true)
+        δV = local_vdim_correction(
+            op,
+            eltype(V),
+            target,
+            source,
+            correction.mesh,
+            correction.bdry_nodes;
+            green_multiplier,
+            correction.maxdist,
+            correction.interpolation_order,
+            correction.quadrature_order,
+            correction.meshsize,
+            shift,
+        )
     else
         error("Unknown correction method. Available options: $CORRECTION_METHODS")
     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 Mesh{N,T} <: AbstractMesh{N,T}
@@ -565,43 +565,68 @@ function connectivity(msh::SubMesh, E::DataType)
 end
 
 """
-    near_interaction_list(X,Y::AbstractMesh; tol)
+    Domain(f::Function, msh::AbstractMesh)
+
+Call `Domain(f, ents)` on `ents = entities(msh).`
+"""
+Domain(f::Function, msh::AbstractMesh) = Domain(f, entities(msh))
 
-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`.
+"""
+    topological_neighbors(msh::LagrangeMesh, k=1)
 
-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`.
+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.
 
-If `tol` is a `Dict`, then `tol[E]` is the tolerance for elements of type `E`.
+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 near_interaction_list(
-    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 = _near_interaction_list(balltree, els, tol_)
-        dict[E] = idxs
+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
-    return dict
-end
-
-@noinline function _near_interaction_list(balltree, els, tol)
-    centers = map(center, els)
-    return inrange(balltree, centers, tol)
+    # 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
+    if k > 1
+        k_neighbors = deepcopy(one_neighbors)
+    else
+        k_neighbors = one_neighbors
+    end
+    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
 
-"""
-    Domain(f::Function, msh::AbstractMesh)
+function viz_elements(args...; kwargs...) end
 
-Call `Domain(f, ents)` on `ents = entities(msh).`
-"""
-Domain(f::Function, msh::AbstractMesh) = Domain(f, entities(msh))
+function viz_elements_bords(args...; kwargs...) end

src/nystrom.jl

@@ -46,7 +46,7 @@ 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{V,K,T,S<:Quadrature} <: AbstractMatrix{V}
+struct IntegralOperator{V,K,T,S} <: AbstractMatrix{V}
     kernel::K
     # since the target can be as simple as a vector of points, leave it untyped
     target::T
@@ -58,15 +58,18 @@ kernel(iop::IntegralOperator) = iop.kernel
 target(iop::IntegralOperator) = iop.target
 source(iop::IntegralOperator) = iop.source
 
-function IntegralOperator(k, X, Y::Quadrature = X)
+function IntegralOperator(k, X, Y = X)
     # check that all entities in the quadrature are of the same dimension
-    if !allequal(geometric_dimension(ent) for ent in entities(Y))
-        msg = "entities in the target quadrature have different geometric dimensions"
-        throw(ArgumentError(msg))
+    if Y isa Quadrature && !isnothing(Y.mesh)
+        if !allequal(geometric_dimension(ent) for ent in entities(Y))
+            msg = "entities in the target quadrature have different geometric dimensions"
+            throw(ArgumentError(msg))
+        end
     end
     T = return_type(k, eltype(X), eltype(Y))
-    msg = """IntegralOperator of nonbits being created: $T"""
-    isbitstype(T) || (@warn msg)
+    # FIXME This cripples performance for local VDIM
+    #msg = """IntegralOperator of nonbits being created: $T"""
+    #isbitstype(T) || (@warn msg)
     return IntegralOperator{T,typeof(k),typeof(X),typeof(Y)}(k, X, Y)
 end
 
@@ -99,7 +102,7 @@ end
 @noinline function _assemble_matrix!(out, K, X, Y::Quadrature, threads)
     @usethreads threads for j in 1:length(Y)
         for i in 1:length(X)
-            out[i, j] = K(X[i], Y[j]) * weight(Y[j])
+            @inbounds out[i, j] = K(X[i], Y[j]) * weight(Y[j])
         end
     end
     return out

src/quadrature.jl

@@ -56,14 +56,16 @@ function Base.show(io::IO, q::QuadratureNode)
     return print(io, "-- weight: $(q.weight)")
 end
 
+const Maybe{T} = Union{T,Nothing}
+
 """
     struct Quadrature{N,T} <: AbstractVector{QuadratureNode{N,T}}
 
 A collection of [`QuadratureNode`](@ref)s used to integrate over an
 [`AbstractMesh`](@ref).
 """
 struct Quadrature{N,T} <: AbstractVector{QuadratureNode{N,T}}
-    mesh::AbstractMesh{N,T}
+    mesh::Maybe{AbstractMesh{N,T}}
     etype2qrule::Dict{DataType,ReferenceQuadrature}
     qnodes::Vector{QuadratureNode{N,T}}
     etype2qtags::Dict{DataType,Matrix{Int}}
@@ -75,7 +77,10 @@ Base.getindex(quad::Quadrature, i) = quad.qnodes[i]
 Base.setindex!(quad::Quadrature, q, i) = (quad.qnodes[i] = q)
 
 qnodes(quad::Quadrature) = quad.qnodes
-mesh(quad::Quadrature) = quad.mesh
+function mesh(quad::Quadrature)
+    isnothing(quad.mesh) && error("The Quadrature has no mesh!")
+    return quad.mesh
+end
 etype2qtags(quad::Quadrature, E) = quad.etype2qtags[E]
 
 quadrature_rule(quad::Quadrature, E) = quad.etype2qrule[E]
@@ -129,6 +134,39 @@ function Quadrature(msh::AbstractMesh{N,T}, etype2qrule::Dict) where {N,T}
     return quad
 end
 
+# Quadrature constructor for list of volume elements for local vdim
+function Quadrature(
+    ::Type{T},
+    elementlist::AbstractVector{E},
+    etype2qrule::Dict{DataType,Q},
+    qrule::Q;
+    center::SVector{N,Float64} = zero(SVector{N,Float64}),
+    scale::Float64 = 1.0,
+) where {N,T,E,Q}
+    # initialize mesh with empty fields
+    quad = Quadrature{N,T}(
+        nothing,
+        etype2qrule,
+        QuadratureNode{N,T}[],
+        Dict{DataType,Matrix{Int}}(),
+    )
+    # loop element types and generate quadrature for each
+    _build_quadrature!(quad, elementlist, qrule; center, scale)
+
+    # check for entities with negative orientation and flip normal vectors if
+    # present
+    #for ent in entities(msh)
+    #    if (sign(tag(ent)) < 0) && (N - geometric_dimension(ent) == 1)
+    #        @debug "Flipping normals of $ent"
+    #        tags = dom2qtags(quad, Domain(ent))
+    #        for i in tags
+    #            quad[i] = flip_normal(quad[i])
+    #        end
+    #    end
+    #end
+    return quad
+end
+
 function Quadrature(msh::AbstractMesh{N,T}, qrule::ReferenceQuadrature) where {N,T}
     etype2qrule = Dict(E => qrule for E in element_types(msh))
     return Quadrature(msh, etype2qrule)
@@ -141,16 +179,18 @@ function Quadrature(msh::AbstractMesh; qorder)
 end
 
 @noinline function _build_quadrature!(
-    quad,
+    quad::Quadrature{N,T},
     els::AbstractVector{E},
-    qrule::ReferenceQuadrature,
-) where {E}
-    N = ambient_dimension(quad)
+    qrule::ReferenceQuadrature;
+    center::SVector{N,Float64} = zero(SVector{N,Float64}),
+    scale::Float64 = 1.0,
+) where {E,N,T}
     x̂, ŵ = qrule() # nodes and weights on reference element
     num_nodes = length(ŵ)
     M = geometric_dimension(domain(E))
     codim = N - M
     istart = length(quad.qnodes) + 1
+    sizehint!(quad.qnodes, length(els) * length(x̂))
     for el in els
         # and all qnodes for that element
         for (x̂i, ŵi) in zip(x̂, ŵ)
@@ -159,7 +199,7 @@ end
             μ = _integration_measure(jac)
             w = μ * ŵi
             ν = codim == 1 ? _normal(jac) : nothing
-            qnode = QuadratureNode(x, w, ν)
+            qnode = QuadratureNode((x - center) / scale, w / scale^M, ν)
             push!(quad.qnodes, qnode)
         end
     end