Implement local version of vdim

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/api.jl

@@ -358,6 +358,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}
@@ -608,6 +608,69 @@ function connectivity(msh::SubMesh, E::DataType)
     return map(t -> g2l[t], view(msh.parent.etype2mat[E], :, eltags))
 end
 
+"""
+    Domain(f::Function, msh::AbstractMesh)
+
+Call `Domain(f, ents)` on `ents = entities(msh).`
+"""
+Domain(f::Function, msh::AbstractMesh) = Domain(f, entities(msh))
+
+"""
+    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
+    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
+
 """
     near_interaction_list(X,Y::AbstractMesh; tol)
 
@@ -643,12 +706,9 @@ end
     return inrange(balltree, centers, tol)
 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
 
 function node2etags(msh)
     # dictionary mapping a node index to all elements containing it. Note

src/nystrom.jl

@@ -76,15 +76,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
 
@@ -117,7 +120,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]
@@ -119,6 +124,40 @@ 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}}(),
+    )
+    ori = ones(Int64, length(elementlist))
+    # loop element types and generate quadrature for each
+    _build_quadrature!(quad, elementlist, ori, 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)
@@ -134,14 +173,18 @@ end
         quad::Quadrature{N, T},
         els::AbstractVector{E},
         orientation::Vector{Int},
-        qrule::ReferenceQuadrature,
-    ) where {N, T, E}
+        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
     x̂ = map(x̂ -> T.(x̂), qcoords(qrule))
     ŵ = map(ŵ -> T.(ŵ), qweights(qrule))
     num_nodes = length(ŵ)
     M = geometric_dimension(domain(E))
     codim = N - M
     istart = length(quad.qnodes) + 1
+    sizehint!(quad.qnodes, length(els) * length(x̂))
     @assert length(els) == length(orientation)
     for (s, el) in zip(orientation, els)
         # and all qnodes for that element
@@ -151,7 +194,7 @@ end
             μ = _integration_measure(jac)
             w = μ * ŵi
             ν = codim == 1 ? T.(s * _normal(jac)) : nothing
-            qnode = QuadratureNode(T.(x), T.(w), ν)
+            qnode = QuadratureNode(T.((x - center) / scale), T.(w / scale^M), ν)
             push!(quad.qnodes, qnode)
         end
     end

src/reference_integration.jl

@@ -169,17 +169,20 @@ struct Gauss{D, N} <: ReferenceQuadrature{D}
     function Gauss(; domain, order)
         domain == :segment && (domain = ReferenceLine())
         domain == :triangle && (domain = ReferenceTriangle())
-        domain == :tetrahedron && (domain = ReferenceTetrahedron())
-        msg = "quadrature of order $order not available for $domain"
+        domain == :tetrehedron && (domain = ReferenceTetrahedron())
         if domain isa ReferenceLine
             # TODO: support Gauss-Legendre quadratures of arbitrary order
-            order == 13 || error(msg)
+            order == 13 || error("quadrature of order $order not available for $domain")
             n = 7
         elseif domain isa ReferenceTriangle
-            haskey(TRIANGLE_GAUSS_ORDER_TO_NPTS, order) || error(msg)
+            if !haskey(TRIANGLE_GAUSS_ORDER_TO_NPTS, order)
+                error("quadrature of order $order not available for $domain")
+            end
             n = TRIANGLE_GAUSS_ORDER_TO_NPTS[order]
         elseif domain isa ReferenceTetrahedron
-            haskey(TETRAHEDRON_GAUSS_ORDER_TO_NPTS, order) || error(msg)
+            if !haskey(TETRAHEDRON_GAUSS_ORDER_TO_NPTS, order)
+                error("quadrature of order $order not available for $domain")
+            end
             n = TETRAHEDRON_GAUSS_ORDER_TO_NPTS[order]
         else
             error(
@@ -272,12 +275,18 @@ struct VioreanuRokhlin{D, N} <: ReferenceQuadrature{D}
         domain == :triangle && (domain = ReferenceTriangle())
         domain == :tetrahedron && (domain = ReferenceTetrahedron())
         if domain isa ReferenceTriangle
-            msg = "VioreanuRokhlin quadrature of order $order not available for ReferenceTriangle"
-            haskey(TRIANGLE_VR_ORDER_TO_NPTS, order) || error(msg)
+            if !haskey(TRIANGLE_VR_ORDER_TO_NPTS, order)
+                error(
+                    "VioreanuRokhlin quadrature of order $order not available for ReferenceTriangle",
+                )
+            end
             n = TRIANGLE_VR_ORDER_TO_NPTS[order]
         elseif domain isa ReferenceTetrahedron
-            msg = "VioreanuRokhlin quadrature of order $order not available for ReferenceTetrahedron"
-            haskey(TETRAHEDRON_VR_ORDER_TO_NPTS, order) || error(msg)
+            if !haskey(TETRAHEDRON_VR_ORDER_TO_NPTS, order)
+                error(
+                    "VioreanuRokhlin quadrature of order $order not available for ReferenceTetrahedron",
+                )
+            end
             n = TETRAHEDRON_VR_ORDER_TO_NPTS[order]
         else
             error(