diff --git a/Project.toml b/Project.toml index 6d5b8ae4..56a1c10f 100644 --- a/Project.toml +++ b/Project.toml @@ -6,10 +6,14 @@ 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" +Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b" 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" diff --git a/ext/IntiMeshesExt.jl b/ext/IntiMeshesExt.jl index 43d577f7..16440ca8 100644 --- a/ext/IntiMeshesExt.jl +++ b/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...) diff --git a/src/Inti.jl b/src/Inti.jl index 55d79610..b2321b24 100644 --- a/src/Inti.jl +++ b/src/Inti.jl @@ -14,6 +14,7 @@ using LinearMaps using NearestNeighbors using SparseArrays using StaticArrays +using QuadGK using SpecialFunctions using Printf diff --git a/src/adaptive.jl b/src/adaptive.jl index 97ea7148..df6b4fdf 100644 --- a/src/adaptive.jl +++ b/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 """ diff --git a/src/bdim.jl b/src/bdim.jl index d6ffb105..74651766 100644 --- a/src/bdim.jl +++ b/src/bdim.jl @@ -177,3 +177,152 @@ function bdim_correction( δD = sparse(Is, Js, Ds, num_trgs, n) return δS, δD end + +function local_bdim_correction( + pde, + target, + source::Quadrature; + green_multiplier::Vector{<:Real}, + parameters = DimParameters(), + derivative::Bool = false, + maxdist = Inf, + kneighbor = 1, +) + imat_cond = imat_norm = res_norm = rhs_norm = -Inf + T = default_kernel_eltype(pde) # Float64 + N = ambient_dimension(source) + m, n = length(target), length(source) + msh = source.mesh + qnodes = source.qnodes + neighbors = topological_neighbors(msh, kneighbor) + 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, 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)*4 + 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(4nq, 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) + x = target[i] + μ = green_multiplier[i] # - 1/2 + 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 diff --git a/src/mesh.jl b/src/mesh.jl index 86e4307e..d71460e4 100644 --- a/src/mesh.jl +++ b/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} @@ -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,13 +538,95 @@ 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::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 + k_neighbors = deepcopy(one_neighbors) + 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 diff --git a/src/quadrature.jl b/src/quadrature.jl index 73fef9d0..94268e6e 100644 --- a/src/quadrature.jl +++ b/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) @@ -276,6 +278,91 @@ function _etype_to_nearest_points(X, Y::Quadrature, maxdist) return etype2nearlist end +""" + etype_to_near_elements(X,Y::Quadrature; tol) + +Return `Nl = [[el in Y.mesh && dist(x, el) ≤ tol] for x in X]` +""" +function etype_to_near_elements(X,Y::Quadrature; tol) + y = [coords(q) for q in Y] + tree = BallTree(y) + etype2nearlist = Dict{DataType,Vector{Set{Int}}}() + for (E, Q) in Y.etype2qtags + P, N = size(Q) + etype2nearlist[E] = source2el = [Set{Int}() for _ in 1:length(X)] + quad2source = [Set{Int}() for _ in 1:length(y)] + for (l, x) in enumerate(X) + for q in inrange(tree, x, tol) + push!(quad2source[q], l) + end + end + for n in 1:N + for i in 1:P + for l in quad2source[Q[i,n]] + push!(source2el[l], n) + end + end + end + end + return etype2nearlist +end + +""" + near_elements(Y::Quadrature; tol) + +Return `Nl = [[el_j in Y.mesh && dist(el_j, el_i) ≤ tol] for el_i in Y.mesh]` +""" +function near_elements(Y::Quadrature; tol) + y = [coords(q) for q in Y] + tree = BallTree(y) + el2el = Dict{Tuple{DataType,Int},Set{Tuple{DataType,Int}}}() + quad2el = [Set{Tuple{DataType,Int}}() for _ in 1:length(y)] + # for each element, loop over its qnodes, find q∈y close to one of its qnodes, add the element to quad2el[q] + # quad2el[q] is the set of elements whose qnodes are close to q + for (E, Q) in Y.etype2qtags + P, N = size(Q) + for n in 1:N + for i in 1:P + for q in inrange(tree, coords(qnodes(Y)[Q[i,n]]), tol) + push!(quad2el[q], (E, n)) + end + end + end + end + # for each element, the set of elements close to it is the + # union of the sets of elements close to its qnodes + for (E, Q) in Y.etype2qtags + P, N = size(Q) + for n in 1:N + el2el[(E, n)] = union([quad2el[Q[i,n]] for i in 1:P]...) + end + end + return el2el +end + +""" + near_components(Y::Quadrature; tol) + +Calculate the connected components of each near_elements +""" +function near_components(Y::Quadrature; tol) + topo_neighbor = topological_neighbors(Y.mesh) + dist_neighbor = near_elements(Y; tol) + return Dict(E => connected_components(nl, topo_neighbor) for (E, nl) in dist_neighbor) +end + +# function _geometric_center_circum_radius(Y::Quadrature, E, Q, P, N) +# C = [(sum(1:P) do i +# coords(qnodes(Y)[Q[i,n]]) .* weight(qnodes(Y)[Q[i,n]]) +# end) / +# (sum(1:P) do i +# weight(qnodes(Y)[Q[i,n]]) +# end) for n in 1:N] +# r = [maximum(i->norm(C[n]-nodes(mesh(Y))[i]), connectivity(mesh(Y), E)[:,n]) for n in 1:N] +# return C, r +# end + + """ quadrature_to_node_vals(Q::Quadrature, qvals::AbstractVector) diff --git a/src/reference_interpolation.jl b/src/reference_interpolation.jl index 112318dd..f00eb358 100644 --- a/src/reference_interpolation.jl +++ b/src/reference_interpolation.jl @@ -262,6 +262,7 @@ const LagrangeCube = LagrangeElement{ReferenceCube} """ vertices_idxs(el::LagrangeElement) + vertices_idxs(::Type{LagrangeElement}) The indices of the nodes in `el` that define the vertices of the element. """ @@ -284,17 +285,76 @@ vertices(el::LagrangeElement) = view(vals(el), vertices_idxs(el)) The indices of the nodes in `el` that define the boundary of the element. """ -function boundary_idxs(el::LagrangeLine) - return 1, length(vals(el)) + +function boundary_idxs(::Type{<:LagrangeLine}) + return 1, 2 end -function boundary_idxs(el::LagrangeTriangle{3}) +function boundary_idxs(::Type{<:LagrangeTriangle{3}}) return (1, 2), (2, 3), (3, 1) end -function boundary_idxs(el::LagrangeTriangle{6}) - return (1, 2), (2, 3), (3, 1) +function boundary_idxs(::Type{<:LagrangeTriangle{6}}) + return (1, 2, 4), (2, 3, 5), (3, 1, 6) +end + +function boundary1d(els, msh) + res = Set{Int}() + E, _ = first(els) + bdi = Inti.boundary_idxs(E) + for (E, i) in els + vertices = Inti.connectivity(msh, E)[:,i] + for bord in (-vertices[bdi[1]], vertices[bdi[2]]) + -bord in res ? delete!(res, -bord) : push!(res, bord) + end + end + return sort([res...]) +end + +function boundarynd(els, msh) + res = Set() + E, _ = first(els) + bdi = Inti.boundary_idxs(E) + for (E, i) in els + vertices = Inti.connectivity(msh, E)[:,i] + bords = [[vertices[i] for i in bi] for bi in bdi] + for new_bord in bords + flag = true + for old_bord in res + if sort(new_bord) == sort(old_bord) + delete!(res, old_bord) + flag = false + end + end + flag && push!(res, new_bord) + end + end + return res +end + +## +function _dfs!(comp, el, nei, els) + for el_nei in nei[el] + if el_nei in els + push!(comp, el_nei) + delete!(els, el_nei) + _dfs!(comp, el_nei, nei, els) + end + end +end + +function connected_components(els, nei) + components = Set{Tuple{DataType,Int}}[] + while !isempty(els) + el = pop!(els) + comp = Set{Tuple{DataType,Int}}() + push!(comp, el) + _dfs!(comp, el, nei, els) + push!(components, comp) + end + return components end +## #= Hardcode some basic elements. diff --git a/test/boundary_test.jl b/test/boundary_test.jl new file mode 100644 index 00000000..44b02ab7 --- /dev/null +++ b/test/boundary_test.jl @@ -0,0 +1,74 @@ +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 = 2) +# Ω, msh = gmsh_ball(; center = [0.0, 0.0, 0.0], radius=1.0, meshsize = 0.2) +Γ = Inti.external_boundary(Ω) + +k = 2 +Γ_msh = msh[Γ] +Ω_msh = msh[Ω] +test_msh = Ω_msh +nei = Inti.topological_neighbors(test_msh, 1) +E = first(Inti.element_types(test_msh)) +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(els, msh) + Els = [Inti.elements(msh, E)[i] for (E, i) in els] + fig, _, _ = viz(Els) + viz!(msh; color = 0, showsegments = true,alpha=0.3) + display(fig) +end + +function viz_elements_bords(Ei, els, bords, msh) + 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, J = 1, 3 +test_els = union(copy(nei[(E,1)]), nei[(E,2)], nei[(E,3)], nei[(E,4)]) +viz_elements(test_els, test_msh) + +components = Inti.connected_components(test_els, nei) + +BD = Inti.boundarynd(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) diff --git a/test/ldim_test.jl b/test/ldim_test.jl new file mode 100644 index 00000000..dcc72e94 --- /dev/null +++ b/test/ldim_test.jl @@ -0,0 +1,74 @@ +using Test +using LinearAlgebra +using Inti +using Random +using Meshes +using Plots + +include("test_utils.jl") +Random.seed!(1) + +N = 2 +t = :interior +pde = Inti.Laplace(; dim = N) + +K = 5:5 +H = [0.2*2.0^(-i) for i in 0:0] +fig = plot(xscale=:log,yscale=:log,xticks=H) +for k in K + err = [] + for h in H + Inti.clear_entities!() + Ω, msh = gmsh_disk(; center = [0.0, 0.0], rx = 1.0, ry = 1.0, meshsize = h, order = 2) + Γ = Inti.external_boundary(Ω) + + ## + + 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) + 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) + + # 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, kneighbor=k) + 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) + push!(err, e1) + end + plot!(fig, H, err;lw=2,marker=:o,label=k) +end +display(fig) \ No newline at end of file diff --git a/test/nearlist_test.jl b/test/nearlist_test.jl new file mode 100644 index 00000000..435f3e09 --- /dev/null +++ b/test/nearlist_test.jl @@ -0,0 +1,38 @@ +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) +h = 0.2 + +Inti.clear_entities!() +Ω, msh = gmsh_disks([([0.0,0.0],1.0,1.0), ([-2.1,0.0],1.0,1.0)];meshsize = h, order = 2) +Γ = Inti.external_boundary(Ω) +quad = Inti.Quadrature(msh[Γ]; qorder = 3) + +# Nl = Inti.near_elements(quad;tol=0.2) +Ncl = Inti.near_components(quad;tol=0.2) +fig, _, _ = viz(msh;showsegments = false,alpha=0.3) + +E = first(keys(Ncl))[1]; i = 1 +viz!(Inti.elements(msh[Γ], E)[i];color=:red) +for (E_, j) in Ncl[(E, i)][2] + viz!(Inti.elements(msh[Γ], E_)[j];color=:blue,alpha=0.3) +end + +# for (E, nl) in Nl +# i = 1 +# viz!(Inti.elements(msh[Γ], E)[i];color=:red) +# for j in nl[i] +# viz!(Inti.elements(msh[Γ], E)[j];color=:blue,alpha=0.3) +# end +# end +display(fig) \ No newline at end of file diff --git a/test/test_utils.jl b/test/test_utils.jl index 5112c107..ba5bcb63 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -1,7 +1,7 @@ using Inti using Gmsh -function gmsh_disk(; center, rx, ry, meshsize) +function gmsh_disk(; center, rx, ry, meshsize, order = 1) msh = try gmsh.initialize() gmsh.option.setNumber("General.Verbosity", 2) @@ -12,6 +12,31 @@ function gmsh_disk(; center, rx, ry, meshsize) gmsh.model.occ.addDisk(center[1], center[2], 0, rx, ry) gmsh.model.occ.synchronize() gmsh.model.mesh.generate(2) + gmsh.model.mesh.setOrder(order) + Inti.import_mesh(; dim = 2) + finally + gmsh.finalize() + end + Ω = Inti.Domain(Inti.entities(msh)) do e + return Inti.geometric_dimension(e) == 2 + end + return Ω, msh +end + +function gmsh_disks(disks;meshsize, order = 1) + msh = try + gmsh.initialize() + gmsh.option.setNumber("General.Verbosity", 2) + gmsh.model.add("disk") + # set max and min meshsize to meshsize + gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) + gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) + for (center, rx, ry) in disks + gmsh.model.occ.addDisk(center[1], center[2], 0, rx, ry) + end + gmsh.model.occ.synchronize() + gmsh.model.mesh.generate(2) + gmsh.model.mesh.setOrder(order) Inti.import_mesh(; dim = 2) finally gmsh.finalize()