From 3b2d94adf14786b5ba5805a016b639a21af24d59 Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Fri, 7 Jun 2024 08:38:47 +0200 Subject: [PATCH 01/13] WIP on local bdim --- src/bdim.jl | 147 ++++++++++++++++++++++++++++++++++++++++++++++ src/mesh.jl | 109 +++++++++++++++++++++++++++++++++- src/quadrature.jl | 2 + 3 files changed, 255 insertions(+), 3 deletions(-) diff --git a/src/bdim.jl b/src/bdim.jl index d6ffb105..82f961c2 100644 --- a/src/bdim.jl +++ b/src/bdim.jl @@ -177,3 +177,150 @@ function bdim_correction( δD = sparse(Is, Js, Ds, num_trgs, n) return δS, δD end + +function local_bdim_correction( + pde, + target, + source::Quadrature, + Sop, + Dop; + green_multiplier::Vector{<:Real}, + parameters = DimParameters(), + derivative::Bool = false, + maxdist = Inf, +) + imat_cond = imat_norm = res_norm = rhs_norm = -Inf + T = eltype(Sop) + 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) + 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, 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] + 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..71017d8c 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} @@ -548,3 +548,106 @@ end centers = map(center, els) return inrange(balltree, centers, tol) 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. +""" +function topological_neighbors(msh::LagrangeMesh, k = 1) + @assert 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 = msh.etype2mat[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 + #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 diff --git a/src/quadrature.jl b/src/quadrature.jl index 73fef9d0..0e84ecfb 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) From 62dd139029c8ae497f57b7a77464f864e99e8fae Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Fri, 7 Jun 2024 09:06:44 +0200 Subject: [PATCH 02/13] begin working towards `local_bdim` --- src/adaptive.jl | 10 ++--- src/bdim.jl | 39 ++++++++--------- src/mesh.jl | 107 ++++++++++++++++------------------------------ test/ldim_test.jl | 69 ++++++++++++++++++++++++++++++ 4 files changed, 130 insertions(+), 95 deletions(-) create mode 100644 test/ldim_test.jl 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 82f961c2..bf10c2e8 100644 --- a/src/bdim.jl +++ b/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 diff --git a/src/mesh.jl b/src/mesh.jl index 71017d8c..0cf511c1 100644 --- a/src/mesh.jl +++ b/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 diff --git a/test/ldim_test.jl b/test/ldim_test.jl new file mode 100644 index 00000000..845c0021 --- /dev/null +++ b/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) From cb315c6d61f98e99d544bf71bef2ce8beec90098 Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Fri, 7 Jun 2024 10:46:00 +0200 Subject: [PATCH 03/13] small changes --- src/reference_interpolation.jl | 1 + test/ldim_test.jl | 2 +- test/test_utils.jl | 3 ++- 3 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/reference_interpolation.jl b/src/reference_interpolation.jl index 112318dd..6d481619 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. """ diff --git a/test/ldim_test.jl b/test/ldim_test.jl index 845c0021..170618e8 100644 --- a/test/ldim_test.jl +++ b/test/ldim_test.jl @@ -12,7 +12,7 @@ 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) +Ω, msh = gmsh_disk(; center = [0.0, 0.0], rx = 1.0, ry = 1.0, meshsize = 0.2, order = 2) Γ = Inti.external_boundary(Ω) Γ_msh = msh[Γ] diff --git a/test/test_utils.jl b/test/test_utils.jl index 5112c107..3d74eff5 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,7 @@ 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() From 814bfd1c8de9071229c51b899f370a52dcd1c3d7 Mon Sep 17 00:00:00 2001 From: Dongchen He Date: Mon, 10 Jun 2024 21:31:15 +0800 Subject: [PATCH 04/13] sketch of a function returning boundary --- src/reference_interpolation.jl | 47 +++++++++++++++++++++++++++++---- test/ldim_test.jl | 48 +++++++++++++++++++++++++++++----- test/test_utils.jl | 3 ++- 3 files changed, 86 insertions(+), 12 deletions(-) diff --git a/src/reference_interpolation.jl b/src/reference_interpolation.jl index 112318dd..a94b2221 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,16 +285,52 @@ 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{P}}) where {P} + return 1, P 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() + E, _ = first(els) + 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) + end + end + return 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 #= diff --git a/test/ldim_test.jl b/test/ldim_test.jl index 845c0021..046a0341 100644 --- a/test/ldim_test.jl +++ b/test/ldim_test.jl @@ -12,24 +12,60 @@ 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) +Ω, 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[Ω] -nei = Inti.topological_neighbors(Ω_msh) |> collect -function viz_neighbors(i) +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, ax, pl = viz(el; color = 0, showsegments = true) + el = Inti.elements(msh, E)[idx] + fig, _, _ = viz(el; color = 0, showsegments = true) for (E, i) in v - el = Inti.elements(Ω_msh, E)[i] + 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) diff --git a/test/test_utils.jl b/test/test_utils.jl index 5112c107..3d74eff5 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,7 @@ 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() From 4934db59bb7736813085519b672a992aa5db7900 Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Tue, 11 Jun 2024 18:01:49 +0200 Subject: [PATCH 05/13] tweak `topological_neighbors` to work on `AbstractMesh` --- src/mesh.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesh.jl b/src/mesh.jl index 0cf511c1..94ba229d 100644 --- a/src/mesh.jl +++ b/src/mesh.jl @@ -590,14 +590,14 @@ 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) +function topological_neighbors(msh::AbstractMesh, k = 1) @assert 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 = msh.etype2mat[E]::Matrix{Int} # connectivity matrix + mat = connectivity(msh, E)::Matrix{Int} # connectivity matrix np, Nel = size(mat) for n in 1:Nel for i in 1:np From 78b74d98a6f1e8532f4c2b92b9d9903dfa40fcc9 Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Tue, 11 Jun 2024 18:27:49 +0200 Subject: [PATCH 06/13] implement `k` topological neighbors --- src/mesh.jl | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/src/mesh.jl b/src/mesh.jl index 94ba229d..d71460e4 100644 --- a/src/mesh.jl +++ b/src/mesh.jl @@ -591,7 +591,6 @@ 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) - @assert 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} @@ -617,6 +616,17 @@ function topological_neighbors(msh::AbstractMesh, k = 1) end end end - #TODO: for k > 1, recursively compute the neighbors from the one-neighbors - return one_neighbors + # 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 From 942af12ea9444c6556d91398c024fa3a1923a8d9 Mon Sep 17 00:00:00 2001 From: Dongchen He Date: Wed, 12 Jun 2024 17:44:51 +0800 Subject: [PATCH 07/13] sketch implement of local dim --- Project.toml | 3 ++ ext/IntiMeshesExt.jl | 4 +- src/Inti.jl | 1 + src/bdim.jl | 38 ++++++++++--------- src/reference_interpolation.jl | 7 ++-- test/boundary_test.jl | 69 ++++++++++++++++++++++++++++++++++ test/ldim_test.jl | 60 +++++------------------------ 7 files changed, 108 insertions(+), 74 deletions(-) create mode 100644 test/boundary_test.jl diff --git a/Project.toml b/Project.toml index 6d5b8ae4..5c679b85 100644 --- a/Project.toml +++ b/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" 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/bdim.jl b/src/bdim.jl index bf10c2e8..c2245ab0 100644 --- a/src/bdim.jl +++ b/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) diff --git a/src/reference_interpolation.jl b/src/reference_interpolation.jl index a94b2221..89e20cb2 100644 --- a/src/reference_interpolation.jl +++ b/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) diff --git a/test/boundary_test.jl b/test/boundary_test.jl new file mode 100644 index 00000000..7a4c2f69 --- /dev/null +++ b/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)) \ No newline at end of file diff --git a/test/ldim_test.jl b/test/ldim_test.jl index 046a0341..316e1b8a 100644 --- a/test/ldim_test.jl +++ b/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 From 6c88ace411838082d0fb28c7e3c0bc2f173983b8 Mon Sep 17 00:00:00 2001 From: Dongchen He Date: Thu, 13 Jun 2024 21:02:43 +0800 Subject: [PATCH 08/13] Implemented k neighbor ldim --- Project.toml | 1 + src/bdim.jl | 7 +- src/reference_interpolation.jl | 8 +-- test/boundary_test.jl | 11 +--- test/ldim_test.jl | 114 ++++++++++++++++++--------------- 5 files changed, 75 insertions(+), 66 deletions(-) diff --git a/Project.toml b/Project.toml index 5c679b85..56a1c10f 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ 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" diff --git a/src/bdim.jl b/src/bdim.jl index c2245ab0..74651766 100644 --- a/src/bdim.jl +++ b/src/bdim.jl @@ -186,6 +186,7 @@ function local_bdim_correction( parameters = DimParameters(), derivative::Bool = false, maxdist = Inf, + kneighbor = 1, ) imat_cond = imat_norm = res_norm = rhs_norm = -Inf T = default_kernel_eltype(pde) # Float64 @@ -193,7 +194,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, 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 @@ -276,11 +277,11 @@ function local_bdim_correction( qnodes_nei = source.qnodes[qtags_nei] jac = jacobian(el, 0.5) ν = -_normal(jac) - h = sum(qnodes[i].weight for i in jglob) + 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(3nq, 0, h) + 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) diff --git a/src/reference_interpolation.jl b/src/reference_interpolation.jl index 89e20cb2..398f031a 100644 --- a/src/reference_interpolation.jl +++ b/src/reference_interpolation.jl @@ -286,8 +286,8 @@ 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(::Type{<:LagrangeLine{P}}) where {P} - return 1, P +function boundary_idxs(::Type{<:LagrangeLine}) + return 1, 2 end function boundary_idxs(::Type{<:LagrangeTriangle{3}}) @@ -299,12 +299,12 @@ function boundary_idxs(::Type{<:LagrangeTriangle{6}}) end function boundary1d(els, msh) - res = Set() + 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]]] + for bord in (-vertices[bdi[1]], vertices[bdi[2]]) -bord in res ? delete!(res, -bord) : push!(res, bord) end end diff --git a/test/boundary_test.jl b/test/boundary_test.jl index 7a4c2f69..cb60dd13 100644 --- a/test/boundary_test.jl +++ b/test/boundary_test.jl @@ -12,14 +12,15 @@ 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_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) |> collect +nei = Inti.topological_neighbors(test_msh, k) |> collect function viz_neighbors(i, msh) k, v = nei[i] E, idx = k @@ -61,9 +62,3 @@ for idxs in BD 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)) \ No newline at end of file diff --git a/test/ldim_test.jl b/test/ldim_test.jl index 9b070d36..ad6dc3db 100644 --- a/test/ldim_test.jl +++ b/test/ldim_test.jl @@ -3,60 +3,72 @@ using LinearAlgebra using Inti using Random using Meshes -using GLMakie +using Plots 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) -Γ = 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) -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) +pde = Inti.Stokes(; dim = N, μ=1.2) + +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 From e2833947d64c40bcbccc322174247b6fd8f697c9 Mon Sep 17 00:00:00 2001 From: Dongchen He Date: Fri, 14 Jun 2024 08:59:15 +0800 Subject: [PATCH 09/13] function connected_components --- src/reference_interpolation.jl | 24 ++++++++++++++++++++++++ test/boundary_test.jl | 26 ++++++++++++++++++-------- 2 files changed, 42 insertions(+), 8 deletions(-) diff --git a/src/reference_interpolation.jl b/src/reference_interpolation.jl index 398f031a..f00eb358 100644 --- a/src/reference_interpolation.jl +++ b/src/reference_interpolation.jl @@ -332,6 +332,30 @@ function boundarynd(els, msh) 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. TODO: Eventually this could/should be automated. diff --git a/test/boundary_test.jl b/test/boundary_test.jl index cb60dd13..44b02ab7 100644 --- a/test/boundary_test.jl +++ b/test/boundary_test.jl @@ -19,8 +19,9 @@ Inti.clear_entities!() k = 2 Γ_msh = msh[Γ] Ω_msh = msh[Ω] -test_msh = Γ_msh -nei = Inti.topological_neighbors(test_msh, k) |> collect +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 @@ -33,8 +34,15 @@ function viz_neighbors(i, msh) 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 = el = Inti.elements(msh, Ei[1])[Ei[2]] + 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 @@ -47,12 +55,14 @@ 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]) +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.boundary1d(test_els, test_msh) -bords = [Inti.nodes(test_msh)[abs(i)] for i in BD] +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 From 12a0bd694bb44dc3405f603f694859d8e114b14b Mon Sep 17 00:00:00 2001 From: Dongchen He Date: Wed, 19 Jun 2024 21:24:43 +0800 Subject: [PATCH 10/13] function etype_to_near_elements and test --- src/quadrature.jl | 56 +++++++++++++++++++++++++++++++++++++++++++ test/ldim_test.jl | 2 +- test/nearlist_test.jl | 31 ++++++++++++++++++++++++ test/test_utils.jl | 24 +++++++++++++++++++ 4 files changed, 112 insertions(+), 1 deletion(-) create mode 100644 test/nearlist_test.jl diff --git a/src/quadrature.jl b/src/quadrature.jl index 0e84ecfb..f1384b00 100644 --- a/src/quadrature.jl +++ b/src/quadrature.jl @@ -278,6 +278,62 @@ function _etype_to_nearest_points(X, Y::Quadrature, maxdist) return etype2nearlist end +""" + etype_to_nearest_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, _ = size(Q) + etype2nearlist[E] = source2el = [Set{Int}() for _ in 1:length(X)] + for (l, x) in enumerate(X) + for q in inrange(tree, x, tol) + push!(source2el[l], ceil(Int, q/P)) + end + end + end + return etype2nearlist +end + +""" + etype_to_nearest_elements(Y::Quadrature; tol) + +Return `Nl = [[el_j in Y.mesh && dist(el_j, el_i) ≤ tol] for el_i in Y.mesh]` +""" +function etype_to_near_elements(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] = el2el = [Set{Int}() for _ in 1:N] + for n in 1:N + for i in 1:P + for q in inrange(tree, coords(qnodes(Y)[Q[i,n]]), tol) + push!(el2el[n], ceil(Int, q/P)) + end + end + end + end + return etype2nearlist +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/test/ldim_test.jl b/test/ldim_test.jl index ad6dc3db..dcc72e94 100644 --- a/test/ldim_test.jl +++ b/test/ldim_test.jl @@ -10,7 +10,7 @@ Random.seed!(1) N = 2 t = :interior -pde = Inti.Stokes(; dim = N, μ=1.2) +pde = Inti.Laplace(; dim = N) K = 5:5 H = [0.2*2.0^(-i) for i in 0:0] diff --git a/test/nearlist_test.jl b/test/nearlist_test.jl new file mode 100644 index 00000000..844be3de --- /dev/null +++ b/test/nearlist_test.jl @@ -0,0 +1,31 @@ +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.etype_to_near_elements(quad;tol=1) +fig, _, _ = viz(msh;showsegments = false,alpha=0.3) + +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 3d74eff5..ba5bcb63 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -23,6 +23,30 @@ function gmsh_disk(; center, rx, ry, meshsize, order = 1) 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() + end + Ω = Inti.Domain(Inti.entities(msh)) do e + return Inti.geometric_dimension(e) == 2 + end + return Ω, msh +end + function gmsh_ball(; center, radius, meshsize) msh = try gmsh.initialize() From 328d717bab2c1c9d998b8c2e19ecbc0e54cf6aa6 Mon Sep 17 00:00:00 2001 From: Dongchen He Date: Thu, 20 Jun 2024 16:50:51 +0800 Subject: [PATCH 11/13] etype_to_near_elements correct version --- src/quadrature.jl | 22 +++++++++++++++++++--- test/nearlist_test.jl | 2 +- 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/src/quadrature.jl b/src/quadrature.jl index f1384b00..04d2209e 100644 --- a/src/quadrature.jl +++ b/src/quadrature.jl @@ -288,11 +288,19 @@ function etype_to_near_elements(X,Y::Quadrature; tol) tree = BallTree(y) etype2nearlist = Dict{DataType,Vector{Set{Int}}}() for (E, Q) in Y.etype2qtags - P, _ = size(Q) + 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!(source2el[l], ceil(Int, q/P)) + 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 @@ -311,13 +319,21 @@ function etype_to_near_elements(Y::Quadrature; tol) for (E, Q) in Y.etype2qtags P, N = size(Q) etype2nearlist[E] = el2el = [Set{Int}() for _ in 1:N] + quad2el = [Set{Int}() for _ in 1:length(y)] for n in 1:N for i in 1:P for q in inrange(tree, coords(qnodes(Y)[Q[i,n]]), tol) - push!(el2el[n], ceil(Int, q/P)) + push!(quad2el[q], n) end end end + for n in 1:N + for i in 1:P + for m in quad2el[Q[i,n]] + push!(el2el[n], m) + end + end + end end return etype2nearlist end diff --git a/test/nearlist_test.jl b/test/nearlist_test.jl index 844be3de..0eced9c7 100644 --- a/test/nearlist_test.jl +++ b/test/nearlist_test.jl @@ -18,7 +18,7 @@ Inti.clear_entities!() Γ = Inti.external_boundary(Ω) quad = Inti.Quadrature(msh[Γ]; qorder = 3) -Nl = Inti.etype_to_near_elements(quad;tol=1) +Nl = Inti.etype_to_near_elements(quad;tol=0.01) fig, _, _ = viz(msh;showsegments = false,alpha=0.3) for (E, nl) in Nl From 56a6a5954b152cbcad640f9648452eadafdafc46 Mon Sep 17 00:00:00 2001 From: Dongchen He Date: Mon, 24 Jun 2024 18:36:39 +0800 Subject: [PATCH 12/13] (etype_to_)near_elements takes different types of elements into account --- src/quadrature.jl | 28 +++++++++++++++------------- test/nearlist_test.jl | 20 +++++++++++++------- 2 files changed, 28 insertions(+), 20 deletions(-) diff --git a/src/quadrature.jl b/src/quadrature.jl index 04d2209e..6ff97214 100644 --- a/src/quadrature.jl +++ b/src/quadrature.jl @@ -279,7 +279,7 @@ function _etype_to_nearest_points(X, Y::Quadrature, maxdist) end """ - etype_to_nearest_elements(X,Y::Quadrature; tol) + etype_to_near_elements(X,Y::Quadrature; tol) Return `Nl = [[el in Y.mesh && dist(x, el) ≤ tol] for x in X]` """ @@ -308,34 +308,36 @@ function etype_to_near_elements(X,Y::Quadrature; tol) end """ - etype_to_nearest_elements(Y::Quadrature; tol) + etype_to_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 etype_to_near_elements(Y::Quadrature; tol) +function near_elements(Y::Quadrature; tol) y = [coords(q) for q in Y] tree = BallTree(y) - etype2nearlist = Dict{DataType,Vector{Set{Int}}}() + 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) - etype2nearlist[E] = el2el = [Set{Int}() for _ in 1:N] - quad2el = [Set{Int}() for _ in 1:length(y)] 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], n) + 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 - for i in 1:P - for m in quad2el[Q[i,n]] - push!(el2el[n], m) - end - end + el2el[(E, n)] = union([quad2el[Q[i,n]] for i in 1:P]...) end end - return etype2nearlist + return el2el end # function _geometric_center_circum_radius(Y::Quadrature, E, Q, P, N) diff --git a/test/nearlist_test.jl b/test/nearlist_test.jl index 0eced9c7..b987533c 100644 --- a/test/nearlist_test.jl +++ b/test/nearlist_test.jl @@ -18,14 +18,20 @@ Inti.clear_entities!() Γ = Inti.external_boundary(Ω) quad = Inti.Quadrature(msh[Γ]; qorder = 3) -Nl = Inti.etype_to_near_elements(quad;tol=0.01) +Nl = Inti.near_elements(quad;tol=0.2) fig, _, _ = viz(msh;showsegments = false,alpha=0.3) -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 +E = first(keys(Nl))[1]; i = 1 +viz!(Inti.elements(msh[Γ], E)[i];color=:red) +for (E_, j) in Nl[(E, i)] + 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 From e609f932d399923e8b76c13aa8c981436a110c8b Mon Sep 17 00:00:00 2001 From: Dongchen He Date: Mon, 24 Jun 2024 18:47:17 +0800 Subject: [PATCH 13/13] near_components --- src/quadrature.jl | 13 ++++++++++++- test/nearlist_test.jl | 7 ++++--- 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/src/quadrature.jl b/src/quadrature.jl index 6ff97214..94268e6e 100644 --- a/src/quadrature.jl +++ b/src/quadrature.jl @@ -308,7 +308,7 @@ function etype_to_near_elements(X,Y::Quadrature; tol) end """ - etype_to_near_elements(Y::Quadrature; tol) + near_elements(Y::Quadrature; tol) Return `Nl = [[el_j in Y.mesh && dist(el_j, el_i) ≤ tol] for el_i in Y.mesh]` """ @@ -340,6 +340,17 @@ function near_elements(Y::Quadrature; tol) 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]]) diff --git a/test/nearlist_test.jl b/test/nearlist_test.jl index b987533c..435f3e09 100644 --- a/test/nearlist_test.jl +++ b/test/nearlist_test.jl @@ -18,12 +18,13 @@ Inti.clear_entities!() Γ = Inti.external_boundary(Ω) quad = Inti.Quadrature(msh[Γ]; qorder = 3) -Nl = Inti.near_elements(quad;tol=0.2) +# 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(Nl))[1]; i = 1 +E = first(keys(Ncl))[1]; i = 1 viz!(Inti.elements(msh[Γ], E)[i];color=:red) -for (E_, j) in Nl[(E, i)] +for (E_, j) in Ncl[(E, i)][2] viz!(Inti.elements(msh[Γ], E_)[j];color=:blue,alpha=0.3) end