sketch implement of local dim

Project.toml

@@ -6,10 +6,13 @@ 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"
 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/bdim.jl

@@ -193,7 +193,7 @@ function local_bdim_correction(
     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
@@ -265,22 +265,26 @@ 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()
-            # qnodes_aux = source.qnodes[jglob]
-            qnodes_aux = source.qnodes # this is the global dim
+            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)
+            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(3nq, 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)

src/reference_interpolation.jl

@@ -304,12 +304,11 @@ function boundary1d(els, msh)
     bdi = Inti.boundary_idxs(E)
     for (E, i) in els
         vertices = Inti.connectivity(msh, E)[:,i]
-        bords = [vertices[bi] for bi in bdi]
-        for bord in bords
-            bord in res ? delete!(res, bord) : push!(res, bord)
+        for bord in [-vertices[bdi[1]], vertices[bdi[2]]]
+            -bord in res ? delete!(res, -bord) : push!(res, bord)
         end
     end
-    return res
+    return sort([res...])
 end
 
 function boundarynd(els, msh)

test/boundary_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, order = 1)
+# Ω, msh = gmsh_ball(; center = [0.0, 0.0, 0.0], radius=1.0, meshsize = 0.2)
+Γ = Inti.external_boundary(Ω)
+
+Γ_msh = msh[Γ]
+Ω_msh = msh[Ω]
+test_msh = Γ_msh
+nei = Inti.topological_neighbors(test_msh) |> collect
+function viz_neighbors(i, msh)
+    k, v = nei[i]
+    E, idx = k
+    el = Inti.elements(msh, E)[idx]
+    fig, _, _ = 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
+
+function viz_elements_bords(Ei, els, bords, msh)
+    el = el = Inti.elements(msh, Ei[1])[Ei[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)
+    display(fig)
+end
+
+# el_in_set(el, set) = any(x->sort(x) == sort(el), set)
+
+I = 10
+test_els = copy(nei[I][2])
+push!(test_els, nei[I][1])
+
+BD = Inti.boundary1d(test_els, test_msh)
+bords = [Inti.nodes(test_msh)[abs(i)] for i in BD]
+
+bords = Inti.LagrangeElement[]
+for idxs in BD
+    vtxs = Inti.nodes(Ω_msh)[idxs]
+    bord = Inti.LagrangeLine(vtxs...)
+    push!(bords, bord)
+end
+
+viz_elements_bords(nei[I][1], test_els, bords, test_msh)
+
+for bord in bords
+    viz!(bord;color=4)
+end
+viz(bords)
+viz(first(bords))

test/ldim_test.jl

@@ -16,59 +16,9 @@ Inti.clear_entities!()
 # Ω, msh = gmsh_ball(; center = [0.0, 0.0, 0.0], radius=1.0, meshsize = 0.2)
 Γ = Inti.external_boundary(Ω)
 
-Γ_msh = msh[Γ]
-Ω_msh = msh[Ω]
-test_msh = Γ_msh
-nei = Inti.topological_neighbors(test_msh) |> collect
-function viz_neighbors(i, msh)
-    k, v = nei[i]
-    E, idx = k
-    el = Inti.elements(msh, E)[idx]
-    fig, _, _ = 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
-
-function viz_elements_bords(els, bords, msh)
-    fig, _, _ = viz(msh; color = 0, showsegments = false,alpha=0.5)
-    for (E, i) in els
-        el = Inti.elements(msh, E)[i]
-        viz!(el; showsegments = false, alpha=0.7)
-    end
-    viz!(bords;color=4,showsegments = false)
-    display(fig)
-end
-
-el_in_set(el, set) = any(x->sort(x) == sort(el), set)
-
-
-I = 5
-test_els = copy(nei[I][2])
-push!(test_els, nei[I][1])
-
-BD = Inti.boundary1d(test_els, test_msh)
-
-bords = Inti.LagrangeElement[]
-for idxs in BD
-    vtxs = Inti.nodes(Ω_msh)[idxs]
-    bord = Inti.LagrangeLine(vtxs...)
-    push!(bords, bord)
-end
-
-viz_elements_bords(test_els, bords, test_msh)
-
-for bord in bords
-    viz!(bord;color=4)
-end
-viz(bords)
-viz(first(bords))
-
 ##
 
-quad = Inti.Quadrature(view(msh, Γ); qorder = 3)
+quad = Inti.Quadrature(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)
@@ -90,6 +40,14 @@ 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)
+
+# qnodes = Inti.local_bdim_correction(pde, quad, quad; green_multiplier)
+# X = [q.coords[1] for q in qnodes]; Y = [q.coords[2] for q in qnodes]
+# u = [q.normal[1] for q in qnodes]; v = [q.normal[2] for q in qnodes]
+# fig, _, _ = scatter(X, Y)
+# arrows!(X, Y, u, v, lengthscale=0.01)
+# display(fig)
+
 δS, δD = Inti.local_bdim_correction(pde, quad, quad; green_multiplier)
 Sdim = Smat + δS
 Ddim = Dmat + δD