diff --git a/ext/IntiMakieExt.jl b/ext/IntiMakieExt.jl index e8a4d509..34f1eecb 100644 --- a/ext/IntiMakieExt.jl +++ b/ext/IntiMakieExt.jl @@ -76,10 +76,10 @@ function Meshes.viz!(el::Inti.ReferenceInterpolant, args...; kwargs...) end function Meshes.viz(els::AbstractVector{<:Inti.ReferenceInterpolant}, args...; kwargs...) - return viz([to_meshes(el) for el in els]) + return viz([to_meshes(el) for el in els], args...; kwargs...) end function Meshes.viz!(els::AbstractVector{<:Inti.ReferenceInterpolant}, args...; kwargs...) - return viz!([to_meshes(el) for el in els]) + return viz!([to_meshes(el) for el in els], args...; kwargs...) end function Meshes.viz(msh::Inti.AbstractMesh, args...; kwargs...) @@ -89,4 +89,45 @@ function Meshes.viz!(msh::Inti.AbstractMesh, args...; kwargs...) return viz!(to_meshes(msh), args...; kwargs...) end +function Inti.viz_elements(els, msh) + E = first(Inti.element_types(msh)) + Els = [Inti.elements(msh, E)[i] for (E, i) in els] + fig, _, _ = viz(Els) + viz!(msh; color = 0, showsegments = true, alpha = 0.3) + return display(fig) +end + +function Inti.viz_elements_bords(Ei, els, ell, bords, msh; quad = nothing) + E = first(Inti.element_types(msh)) + #ell = collect(Ei[(E, 1)])[1] + el = Inti.elements(msh, ell[1])[ell[2]] + fig, _, _ = viz(msh; color = 0, showsegments = true, alpha = 0.3) + viz!(el; color = 0, showsegments = true, alpha = 0.5) + for (E, i) in els + el = Inti.elements(msh, E)[i] + viz!(el; showsegments = true, alpha = 0.7) + end + viz!(bords; color = 4, showsegments = false, segmentsize = 5, segmentcolor = 4) + # if !isnothing(quad) + # xs = [qnode.coords[1] for qnode in quad.qnodes] + # ys = [qnode.coords[2] for qnode in quad.qnodes] + # zs = [qnode.coords[3] for qnode in quad.qnodes] + # nx = [Inti.normal(qnode)[1] for qnode in quad.qnodes] + # ny = [Inti.normal(qnode)[2] for qnode in quad.qnodes] + # nz = [Inti.normal(qnode)[3] for qnode in quad.qnodes] + # Mke.arrows!( + # xs, + # ys, + # zs, + # nx, + # ny, + # nz; + # lengthscale = 0.05, + # linewidth = 0.01, + # arrowsize = 0.01, + # ) + # end + return display(fig) +end + end # module diff --git a/src/api.jl b/src/api.jl index ed9862a8..9fa96df2 100644 --- a/src/api.jl +++ b/src/api.jl @@ -358,6 +358,25 @@ function volume_potential(; op, target, source::Quadrature, compression, correct correction.maxdist, interpolation_order, ) + elseif correction.method == :ldim + loc = target === source ? :inside : correction.target_location + μ = _green_multiplier(loc) + green_multiplier = fill(μ, length(target)) + shift = Val(true) + δV = local_vdim_correction( + op, + eltype(V), + target, + source, + correction.mesh, + correction.bdry_nodes; + green_multiplier, + correction.maxdist, + correction.interpolation_order, + correction.quadrature_order, + correction.meshsize, + shift, + ) else error("Unknown correction method. Available options: $CORRECTION_METHODS") end diff --git a/src/mesh.jl b/src/mesh.jl index 066072fc..b0f2ddd8 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 Mesh{N,T} <: AbstractMesh{N,T} @@ -608,6 +608,69 @@ function connectivity(msh::SubMesh, E::DataType) return map(t -> g2l[t], view(msh.parent.etype2mat[E], :, eltags)) end +""" + Domain(f::Function, msh::AbstractMesh) + +Call `Domain(f, ents)` on `ents = entities(msh).` +""" +Domain(f::Function, msh::AbstractMesh) = Domain(f, entities(msh)) + +""" + topological_neighbors(msh::LagrangeMesh, k=1) + +Return the `k` neighbors of each element in `msh`. The one-neighbors are the +elements that share a common vertex with the element, `k` neighbors are the +one-neighbors of the `k-1` neighbors. + +This function returns a dictionary where the key is an `(Eᵢ,i)` tuple denoting +the element `i` of type `E` in the mesh, and the value is a set of tuples +`(Eⱼ,j)` where `Eⱼ` is the type of the element and `j` its index. +""" +function topological_neighbors(msh::AbstractMesh, k = 1) + # dictionary mapping a node index to all elements containing it. Note + # that the elements are stored as a tuple (type, index) + T = Tuple{DataType, Int} + node2els = Dict{Int, Vector{T}}() + for E in element_types(msh) + mat = connectivity(msh, E)::Matrix{Int} # connectivity matrix + np, Nel = size(mat) + for n in 1:Nel + for i in 1:np + idx = mat[i, n] + els = get!(node2els, idx, Vector{T}()) + push!(els, (E, n)) + end + end + end + # now revert the map to get the neighbors + one_neighbors = Dict{T, Set{T}}() + for (_, els) in node2els + for el in els + nei = get!(one_neighbors, el, Set{T}()) + for el′ in els + push!(nei, el′) + end + end + end + # Recursively compute the neighbors from the one-neighbors + if k > 1 + k_neighbors = deepcopy(one_neighbors) + else + k_neighbors = one_neighbors + end + while k > 1 + # update neighborhood of each element + for el in keys(one_neighbors) + knn = k_neighbors[el] + for el′ in copy(knn) + union!(knn, one_neighbors[el′]) + end + end + k -= 1 + end + return k_neighbors +end + """ near_interaction_list(X,Y::AbstractMesh; tol) @@ -643,12 +706,9 @@ end return inrange(balltree, centers, tol) end -""" - Domain(f::Function, msh::AbstractMesh) +function viz_elements(args...; kwargs...) end -Call `Domain(f, ents)` on `ents = entities(msh).` -""" -Domain(f::Function, msh::AbstractMesh) = Domain(f, entities(msh)) +function viz_elements_bords(args...; kwargs...) end function node2etags(msh) # dictionary mapping a node index to all elements containing it. Note diff --git a/src/nystrom.jl b/src/nystrom.jl index 86af35fc..69c74bac 100644 --- a/src/nystrom.jl +++ b/src/nystrom.jl @@ -76,15 +76,18 @@ kernel(iop::IntegralOperator) = iop.kernel target(iop::IntegralOperator) = iop.target source(iop::IntegralOperator) = iop.source -function IntegralOperator(k, X, Y::Quadrature = X) +function IntegralOperator(k, X, Y = X) # check that all entities in the quadrature are of the same dimension - if !allequal(geometric_dimension(ent) for ent in entities(Y)) - msg = "entities in the target quadrature have different geometric dimensions" - throw(ArgumentError(msg)) + if Y isa Quadrature && !isnothing(Y.mesh) + if !allequal(geometric_dimension(ent) for ent in entities(Y)) + msg = "entities in the target quadrature have different geometric dimensions" + throw(ArgumentError(msg)) + end end T = return_type(k, eltype(X), eltype(Y)) - msg = """IntegralOperator of nonbits being created: $T""" - isbitstype(T) || (@warn msg) + # FIXME This cripples performance for local VDIM + #msg = """IntegralOperator of nonbits being created: $T""" + #isbitstype(T) || (@warn msg) return IntegralOperator{T, typeof(k), typeof(X), typeof(Y)}(k, X, Y) end @@ -117,7 +120,7 @@ end @noinline function _assemble_matrix!(out, K, X, Y::Quadrature, threads) @usethreads threads for j in 1:length(Y) for i in 1:length(X) - out[i, j] = K(X[i], Y[j]) * weight(Y[j]) + @inbounds out[i, j] = K(X[i], Y[j]) * weight(Y[j]) end end return out diff --git a/src/quadrature.jl b/src/quadrature.jl index af2d02b4..131e859b 100644 --- a/src/quadrature.jl +++ b/src/quadrature.jl @@ -56,6 +56,8 @@ function Base.show(io::IO, q::QuadratureNode) return print(io, "-- weight: $(q.weight)") end +const Maybe{T} = Union{T, Nothing} + """ struct Quadrature{N,T} <: AbstractVector{QuadratureNode{N,T}} @@ -63,7 +65,7 @@ A collection of [`QuadratureNode`](@ref)s used to integrate over an [`AbstractMesh`](@ref). """ struct Quadrature{N, T} <: AbstractVector{QuadratureNode{N, T}} - mesh::AbstractMesh{N, T} + mesh::Maybe{AbstractMesh{N, T}} etype2qrule::Dict{DataType, ReferenceQuadrature} qnodes::Vector{QuadratureNode{N, T}} etype2qtags::Dict{DataType, Matrix{Int}} @@ -75,7 +77,10 @@ Base.getindex(quad::Quadrature, i) = quad.qnodes[i] Base.setindex!(quad::Quadrature, q, i) = (quad.qnodes[i] = q) qnodes(quad::Quadrature) = quad.qnodes -mesh(quad::Quadrature) = quad.mesh +function mesh(quad::Quadrature) + isnothing(quad.mesh) && error("The Quadrature has no mesh!") + return quad.mesh +end etype2qtags(quad::Quadrature, E) = quad.etype2qtags[E] quadrature_rule(quad::Quadrature, E) = quad.etype2qrule[E] @@ -119,6 +124,40 @@ function Quadrature(msh::AbstractMesh{N, T}, etype2qrule::Dict) where {N, T} return quad end +# Quadrature constructor for list of volume elements for local vdim +function Quadrature( + ::Type{T}, + elementlist::AbstractVector{E}, + etype2qrule::Dict{DataType, Q}, + qrule::Q; + center::SVector{N, Float64} = zero(SVector{N, Float64}), + scale::Float64 = 1.0, + ) where {N, T, E, Q} + # initialize mesh with empty fields + quad = Quadrature{N, T}( + nothing, + etype2qrule, + QuadratureNode{N, T}[], + Dict{DataType, Matrix{Int}}(), + ) + ori = ones(Int64, length(elementlist)) + # loop element types and generate quadrature for each + _build_quadrature!(quad, elementlist, ori, qrule; center, scale) + + # check for entities with negative orientation and flip normal vectors if + # present + #for ent in entities(msh) + # if (sign(tag(ent)) < 0) && (N - geometric_dimension(ent) == 1) + # @debug "Flipping normals of $ent" + # tags = dom2qtags(quad, Domain(ent)) + # for i in tags + # quad[i] = flip_normal(quad[i]) + # end + # end + #end + return quad +end + function Quadrature(msh::AbstractMesh{N, T}, qrule::ReferenceQuadrature) where {N, T} etype2qrule = Dict(E => qrule for E in element_types(msh)) return Quadrature(msh, etype2qrule) @@ -134,14 +173,18 @@ end quad::Quadrature{N, T}, els::AbstractVector{E}, orientation::Vector{Int}, - qrule::ReferenceQuadrature, - ) where {N, T, E} + qrule::ReferenceQuadrature; + center::SVector{N, Float64} = zero(SVector{N, Float64}), + scale::Float64 = 1.0, + ) where {E, N, T} + x̂, ŵ = qrule() # nodes and weights on reference element x̂ = map(x̂ -> T.(x̂), qcoords(qrule)) ŵ = map(ŵ -> T.(ŵ), qweights(qrule)) num_nodes = length(ŵ) M = geometric_dimension(domain(E)) codim = N - M istart = length(quad.qnodes) + 1 + sizehint!(quad.qnodes, length(els) * length(x̂)) @assert length(els) == length(orientation) for (s, el) in zip(orientation, els) # and all qnodes for that element @@ -151,7 +194,7 @@ end μ = _integration_measure(jac) w = μ * ŵi ν = codim == 1 ? T.(s * _normal(jac)) : nothing - qnode = QuadratureNode(T.(x), T.(w), ν) + qnode = QuadratureNode(T.((x - center) / scale), T.(w / scale^M), ν) push!(quad.qnodes, qnode) end end diff --git a/src/reference_integration.jl b/src/reference_integration.jl index 2696a2e1..4c2a64a9 100644 --- a/src/reference_integration.jl +++ b/src/reference_integration.jl @@ -169,17 +169,20 @@ struct Gauss{D, N} <: ReferenceQuadrature{D} function Gauss(; domain, order) domain == :segment && (domain = ReferenceLine()) domain == :triangle && (domain = ReferenceTriangle()) - domain == :tetrahedron && (domain = ReferenceTetrahedron()) - msg = "quadrature of order $order not available for $domain" + domain == :tetrehedron && (domain = ReferenceTetrahedron()) if domain isa ReferenceLine # TODO: support Gauss-Legendre quadratures of arbitrary order - order == 13 || error(msg) + order == 13 || error("quadrature of order $order not available for $domain") n = 7 elseif domain isa ReferenceTriangle - haskey(TRIANGLE_GAUSS_ORDER_TO_NPTS, order) || error(msg) + if !haskey(TRIANGLE_GAUSS_ORDER_TO_NPTS, order) + error("quadrature of order $order not available for $domain") + end n = TRIANGLE_GAUSS_ORDER_TO_NPTS[order] elseif domain isa ReferenceTetrahedron - haskey(TETRAHEDRON_GAUSS_ORDER_TO_NPTS, order) || error(msg) + if !haskey(TETRAHEDRON_GAUSS_ORDER_TO_NPTS, order) + error("quadrature of order $order not available for $domain") + end n = TETRAHEDRON_GAUSS_ORDER_TO_NPTS[order] else error( @@ -272,12 +275,18 @@ struct VioreanuRokhlin{D, N} <: ReferenceQuadrature{D} domain == :triangle && (domain = ReferenceTriangle()) domain == :tetrahedron && (domain = ReferenceTetrahedron()) if domain isa ReferenceTriangle - msg = "VioreanuRokhlin quadrature of order $order not available for ReferenceTriangle" - haskey(TRIANGLE_VR_ORDER_TO_NPTS, order) || error(msg) + if !haskey(TRIANGLE_VR_ORDER_TO_NPTS, order) + error( + "VioreanuRokhlin quadrature of order $order not available for ReferenceTriangle", + ) + end n = TRIANGLE_VR_ORDER_TO_NPTS[order] elseif domain isa ReferenceTetrahedron - msg = "VioreanuRokhlin quadrature of order $order not available for ReferenceTetrahedron" - haskey(TETRAHEDRON_VR_ORDER_TO_NPTS, order) || error(msg) + if !haskey(TETRAHEDRON_VR_ORDER_TO_NPTS, order) + error( + "VioreanuRokhlin quadrature of order $order not available for ReferenceTetrahedron", + ) + end n = TETRAHEDRON_VR_ORDER_TO_NPTS[order] else error( diff --git a/src/reference_interpolation.jl b/src/reference_interpolation.jl index 58329639..18f9be87 100644 --- a/src/reference_interpolation.jl +++ b/src/reference_interpolation.jl @@ -314,6 +314,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. """ @@ -356,17 +357,18 @@ 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) - I = vertices_idxs(el) +function boundary_idxs(T::Type{<:LagrangeTriangle}) + I = vertices_idxs(T) return (I[1], I[2]), (I[2], I[3]), (I[3], I[1]) end -function boundary_idxs(el::LagrangeSquare) - I = vertices_idxs(el) +function boundary_idxs(T::Type{<:LagrangeSquare}) + I = vertices_idxs(T) return (I[1], I[2]), (I[2], I[3]), (I[3], I[4]), (I[4], I[1]) end @@ -562,3 +564,43 @@ function lagrange_basis(::Type{LagrangeElement{D, N, T}}) where {D, N, T} vals = svector(i -> svector(j -> i == j, N), N) return LagrangeElement{D}(vals) end + +function boundarynd(::Type{T}, els, msh) where {T} + bdi = Inti.boundary_idxs(T) + nedges = length(els) * length(bdi) + edgelist = Vector{SVector{length(bdi[1]), Int64}}(undef, nedges) + edgelist_unsrt = Vector{SVector{length(bdi[1]), Int64}}(undef, nedges) + bords = Vector{MVector{length(bdi[1]), Int64}}(undef, length(bdi)) + for i in 1:length(bdi) + bords[i] = MVector{length(bdi[1]), Int64}(undef) + end + j = 1 + for ii in els + for k in 1:length(bdi) + for jjj in 1:length(bdi[k]) + bords[k][jjj] = Inti.connectivity(msh, T)[bdi[k][jjj], ii] + end + end + for q in bords + edgelist_unsrt[j] = q[:] + edgelist[j] = sort!(q) + j += 1 + end + end + I = sortperm(edgelist) + uniqlist = Int64[] + sizehint!(uniqlist, length(els)) + i = 1 + while i <= length(edgelist) - 1 + if isequal(edgelist[I[i]], edgelist[I[i + 1]]) + i += 1 + else + push!(uniqlist, I[i]) + end + i += 1 + end + if !isequal(edgelist[I[end - 1]], edgelist[I[end]]) + push!(uniqlist, I[end]) + end + return edgelist_unsrt[uniqlist] +end diff --git a/src/utils.jl b/src/utils.jl index 355f8d22..0d018524 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -198,12 +198,12 @@ function _normalize_compression(compression, target, source) end function _normalize_correction(correction, target, source) - methods = (:dim, :adaptive, :none) + methods = (:dim, :ldim, :adaptive, :none) # check that method is valid correction.method ∈ methods || error("Unknown correction.method $(correction.method). Available options: $methods") # set default values if absent - if correction.method == :dim + if correction.method == :dim || correction.method == :ldim haskey(correction, :target_location) && target === source && correction.target_location != :on && @@ -219,6 +219,10 @@ function _normalize_correction(correction, target, source) haskey(correction, :maxdist) || target === source || @warn("missing maxdist field in correction: setting to Inf") + if correction.method == :ldim + haskey(correction, :mesh) || + error("missing mesh information needed for local dim") + end correction = merge( (maxdist = Inf, interpolation_order = nothing, center = nothing), correction, diff --git a/src/vdim.jl b/src/vdim.jl index 88243f30..bad5ebc0 100644 --- a/src/vdim.jl +++ b/src/vdim.jl @@ -27,6 +27,7 @@ See [anderson2024fast](@cite) for more details on the method. - `shift`: a boolean indicating whether the basis functions should be shifted and rescaled to each element. """ + function vdim_correction( op, target, @@ -128,6 +129,367 @@ function vdim_correction( return δV end +# function barrier for type stability purposes +function build_vander(vals_trg, pts, PFE_p, c, r) + tmp = Vector{Float64}(undef, length(c)) + for i in 1:length(pts) + tmp .= (pts[i].coords - c) / r + ElementaryPDESolutions.fast_evaluate!(view(vals_trg, :, i), tmp, PFE_p) + end + return vals_trg +end + +function _scaled_operator(op::AbstractDifferentialOperator, scale) + if op isa Helmholtz + return Helmholtz(; k = scale * op.k, dim = ambient_dimension(op)) + elseif op isa Laplace + return op + else + error("Unsupported operator for stabilized Local VDIM") + end +end + +function _lowfreq_operator(op::AbstractDifferentialOperator{N}) where {N} + if op isa Helmholtz + return Laplace(; dim = N) + elseif op isa Laplace + return Laplace(; dim = N) + else + error("Unsupported operator for stabilized Local VDIM") + end +end + +function local_vdim_correction( + op, + ::Type{Eltype}, + target, + source::Quadrature, + mesh::AbstractMesh, + bdry_nodes; + green_multiplier::Vector{<:Real}, + interpolation_order = nothing, + quadrature_order = nothing, + meshsize = 1.0, + maxdist = Inf, + center = nothing, + shift::Val{SHIFT} = Val(false), + ) where {SHIFT, Eltype} + # variables for debugging the condition properties of the method + vander_cond = vander_norm = rhs_norm = res_norm = shift_norm = -Inf + # figure out if we are dealing with a scalar or vector PDE + m, n = length(target), length(source) + N = ambient_dimension(op) + @assert ambient_dimension(source) == N "vdim only works for volume potentials" + m, n = length(target), length(source) + # a reasonable interpolation_order if not provided + isnothing(interpolation_order) && + (interpolation_order = maximum(order, values(source.etype2qrule))) + + # Helmholtz PDE operator in x̂ coordinates where x = scale * x̂ + s = meshsize + op_hat = _scaled_operator(op, s) + op_lowfreq = _lowfreq_operator(op) + PFE_p_lowfreq, PFE_P_lowfreq, multiindices_lowfreq, monomials_indices_lowfreq = + polynomial_solutions_local_vdim(op_lowfreq, interpolation_order + 2) + PFE_p, PFE_P, multiindices, monomials_indices = + polynomial_solutions_local_vdim(op_hat, interpolation_order) + + dict_near = etype_to_nearest_points(target, source; maxdist) + bdry_kdtree = KDTree(bdry_nodes) + # compute sparse correction + Is = Int[] + Js = Int[] + Vs = Eltype[] + for (E, qtags) in source.etype2qtags + els = elements(source.mesh, E) + near_list = dict_near[E] + nq, ne = size(qtags) + @assert length(near_list) == ne + sizehint!(Is, ne * nq * nq) + sizehint!(Js, ne * nq * nq) + sizehint!(Vs, ne * nq * nq) + num_basis = binomial(interpolation_order + N, N) + L̃ = Matrix{Float64}(undef, nq, num_basis) + vals_trg = Matrix{Float64}(undef, num_basis, nq) + + bdry_qorder = 2 * quadrature_order + if N == 3 + bdry_qrule = _qrule_for_reference_shape(ReferenceSimplex{2}(), bdry_qorder) + bdry_etype2qrule = Dict(ReferenceSimplex{2} => bdry_qrule) + else + bdry_qrule = _qrule_for_reference_shape(ReferenceHyperCube{1}(), bdry_qorder) + bdry_etype2qrule = Dict(ReferenceHyperCube{1} => bdry_qrule) + end + vol_qrule = VioreanuRokhlin(; domain = domain(E), order = quadrature_order) + vol_etype2qrule = Dict(E => vol_qrule) + + topo_neighs = 1 + neighbors = topological_neighbors(mesh, topo_neighs) + + for n in 1:ne + # indices of nodes in element `n` + isempty(near_list[n]) && continue + c, r, diam = translation_and_scaling(els[n]) + s = 1.0 + if false + #if r * op.k < 10^(-3) + lowfreq = true + Yvol, Ybdry, need_layer_corr = _local_vdim_construct_local_quadratures( + N, + mesh, + neighbors, + n, + c, + r, + diam, + bdry_kdtree, + bdry_etype2qrule, + vol_etype2qrule, + bdry_qrule, + vol_qrule, + ) + R = _lowfreq_vdim_auxiliary_quantities( + op, + op_lowfreq, + c, + r, + num_basis, + PFE_p_lowfreq, + PFE_P_lowfreq, + multiindices, + multiindices_lowfreq, + monomials_indices, + monomials_indices_lowfreq, + target[near_list[n]], + green_multiplier, + Yvol, + Ybdry, + diam, + need_layer_corr + ) + Yvol_s1, Ybdry_s1, need_layer_corr_s1 = _local_vdim_construct_local_quadratures( + N, + mesh, + neighbors, + n, + c, + 1.0, + diam, + bdry_kdtree, + bdry_etype2qrule, + vol_etype2qrule, + bdry_qrule, + vol_qrule, + ) + R_s1, b_s1 = _local_vdim_auxiliary_quantities( + op_hat, + c, + 1.0, + PFE_p, + PFE_P, + target[near_list[n]], + green_multiplier, + Yvol_s1, + Ybdry_s1, + diam, + need_layer_corr_s1 + ) + else + lowfreq = false + Yvol, Ybdry, need_layer_corr = _local_vdim_construct_local_quadratures( + N, + mesh, + neighbors, + n, + c, + s, + diam, + bdry_kdtree, + bdry_etype2qrule, + vol_etype2qrule, + bdry_qrule, + vol_qrule, + ) + #if isdefined(Main, :Infiltrator) + # Main.infiltrate(@__MODULE__, Base.@locals, @__FILE__, @__LINE__) + # end + R, b = _local_vdim_auxiliary_quantities( + op_hat, + c, + s, + PFE_p, + PFE_P, + target[near_list[n]], + green_multiplier, + Yvol, + Ybdry, + diam, + need_layer_corr + ) + Yvol_lowfreq, Ybdry_lowfreq, need_layer_corr_lowfreq = _local_vdim_construct_local_quadratures( + N, + mesh, + neighbors, + n, + c, + r, + diam, + bdry_kdtree, + bdry_etype2qrule, + vol_etype2qrule, + bdry_qrule, + vol_qrule, + ) + R_lowfreq = _lowfreq_vdim_auxiliary_quantities( + op, + op_lowfreq, + c, + r, + num_basis, + PFE_p_lowfreq, + PFE_P_lowfreq, + multiindices, + multiindices_lowfreq, + monomials_indices, + monomials_indices_lowfreq, + target[near_list[n]], + green_multiplier, + Yvol_lowfreq, + Ybdry_lowfreq, + diam, + need_layer_corr_lowfreq + ) + #Yvol_s1, Ybdry_s1, need_layer_corr_s1 = _local_vdim_construct_local_quadratures( + # N, + # mesh, + # neighbors, + # n, + # c, + # 1.0, + # diam, + # bdry_kdtree, + # bdry_etype2qrule, + # vol_etype2qrule, + # bdry_qrule, + # vol_qrule, + #) + #R_s1, b_s1 = _local_vdim_auxiliary_quantities( + # op_hat, + # c, + # 1.0, + # PFE_p, + # PFE_P, + # target[near_list[n]], + # green_multiplier, + # Yvol_s1, + # Ybdry_s1, + # diam, + # need_layer_corr_s1; + #) + #if isdefined(Main, :Infiltrator) + # Main.infiltrate(@__MODULE__, Base.@locals, @__FILE__, @__LINE__) + # end + end + jglob = @view qtags[:, n] + # compute translation and scaling + if SHIFT + # TODO copy this from (part of) the output of _local_vdim_auxiliary_quantities ? + L̃ .= transpose(build_vander(vals_trg, view(source, jglob), PFE_p, c, r)) + Linv = pinv(L̃) + if !lowfreq + S = s^2 * Diagonal((s / r) .^ (abs.(multiindices))) + #S_s1 = Diagonal((1 / r) .^ (abs.(multiindices))) + wei = transpose(Linv) * S * transpose(R) + #wei_s1 = transpose(Linv) * S_s1 * transpose(R_s1) + + #area = 0.0 + #for i in 1:length(Yvol) + # area += Yvol[i].weight + #end + #area_s1 = 0.0 + #for i in 1:length(Yvol_s1) + # area_s1 += Yvol_s1[i].weight + #end + #area_contour = 0.0 + #for i in 1:length(Ybdry) + # x = Ybdry[i].coords[1] + # y = Ybdry[i].coords[2] + # nx = Ybdry[i].normal[1] + # ny = Ybdry[i].normal[2] + # area_contour += (x/2 * nx + y/2*ny) * Ybdry[i].weight + #end + #area_contour_s1 = 0.0 + #for i in 1:length(Ybdry_s1) + # x = Ybdry_s1[i].coords[1] + # y = Ybdry_s1[i].coords[2] + # nx = Ybdry_s1[i].normal[1] + # ny = Ybdry_s1[i].normal[2] + # area_contour_s1 += (x/2 * nx + y/2*ny) * Ybdry_s1[i].weight + #end + + #Vint = 0.0 + #for i in 1:length(Yvol) + # x = Yvol[i].coords[1] + # y = Yvol[i].coords[2] + # Vint += (3*x^2 + 3*y^2) * Yvol[i].weight + #end + #Vint_s1 = 0.0 + #for i in 1:length(Yvol) + # x = Yvol_s1[i].coords[1] + # y = Yvol_s1[i].coords[2] + # Vint_s1 += (3*x^2 + 3*y^2) * Yvol_s1[i].weight + #end + #Vcontour = 0.0 + #for i in 1:length(Ybdry) + # x = Ybdry[i].coords[1] + # y = Ybdry[i].coords[2] + # nx = Ybdry[i].normal[1] + # ny = Ybdry[i].normal[2] + # Vcontour += (x^3 * nx + y^3*ny) * Ybdry[i].weight + #end + #Vcontour_s1 = 0.0 + #for i in 1:length(Ybdry_s1) + # x = Ybdry_s1[i].coords[1] + # y = Ybdry_s1[i].coords[2] + # nx = Ybdry_s1[i].normal[1] + # ny = Ybdry_s1[i].normal[2] + # Vcontour_s1 += (x^3 * nx + y^3*ny) * Ybdry_s1[i].weight + #end + if isdefined(Main, :Infiltrator) + Main.infiltrate(@__MODULE__, Base.@locals, @__FILE__, @__LINE__) + end + else + D = Vector{Float64}(undef, num_basis) + D .= r^2 + D = Diagonal(D) + wei = transpose(Linv) * D * transpose(R) + + S_s1 = Diagonal((1 / r) .^ (abs.(multiindices))) + wei_s1 = transpose(Linv) * S_s1 * transpose(R_s1) + if isdefined(Main, :Infiltrator) + Main.infiltrate(@__MODULE__, Base.@locals, @__FILE__, @__LINE__) + end + end + else + error("unsupported local VDIM without shifting") + end + # correct each target near the current element + append!(Is, repeat(near_list[n]; inner = nq)) + append!(Js, repeat(jglob; outer = length(near_list[n]))) + append!(Vs, wei) + end + end + @debug """Condition properties of vdim correction: + |-- max interp. matrix condition: $vander_cond + |-- max norm of source term: $rhs_norm + |-- max residual error: $res_norm + |-- max interp. matrix norm : $vander_norm + |-- max shift norm : $shift_norm + """ + δV = sparse(Is, Js, Vs, m, n) + return δV +end + function change_of_basis(multiindices, p, c, r) nbasis = length(multiindices) P = zeros(nbasis, nbasis) @@ -151,6 +513,7 @@ function translation_and_scaling(el::LagrangeTriangle) l1 = norm(vertices[1] - vertices[2]) l2 = norm(vertices[2] - vertices[3]) l3 = norm(vertices[3] - vertices[1]) + diam = max(l1, l2, l3) if ((l1^2 + l2^2 >= l3^2) && (l2^2 + l3^2 >= l1^2) && (l3^2 + l1^2 > l2^2)) acuteright = true else @@ -163,7 +526,7 @@ function translation_and_scaling(el::LagrangeTriangle) Cp = vertices[3] - vertices[1] Dp = 2 * (Bp[1] * Cp[2] - Bp[2] * Cp[1]) Upx = 1 / Dp * (Cp[2] * (Bp[1]^2 + Bp[2]^2) - Bp[2] * (Cp[1]^2 + Cp[2]^2)) - Upy = 1 / Dp * (Bp[1] * (Cp[1]^2 + Cp[2]^2) - Cp[2] * (Bp[1]^2 + Bp[2]^2)) + Upy = 1 / Dp * (Bp[1] * (Cp[1]^2 + Cp[2]^2) - Cp[1] * (Bp[1]^2 + Bp[2]^2)) Up = SVector{2}(Upx, Upy) r = norm(Up) c = Up + vertices[1] @@ -179,7 +542,7 @@ function translation_and_scaling(el::LagrangeTriangle) r = l3 / 2 end end - return c, r + return c, r, diam end function translation_and_scaling(el::ParametricElement{ReferenceSimplex{3}}) @@ -207,6 +570,7 @@ function translation_and_scaling(el::LagrangeTetrahedron) d = norm(vertices[3] - vertices[2]) e = norm(vertices[3] - vertices[1]) f = norm(vertices[2] - vertices[1]) + diam = max(a, b, c, d, e, f) f² = f^2 a² = a^2 b² = b^2 @@ -254,7 +618,297 @@ function translation_and_scaling(el::LagrangeTetrahedron) R = f / 2 end end - return center, R + return center, R, diam +end + +# function barrier for type stability purposes +_newbord_line(vtxs) = LagrangeLine(SVector{2}(vtxs)) + +# function barrier for type stability purposes +_newbord_tri(vtxs) = LagrangeElement{ReferenceSimplex{2}}(SVector{3}(vtxs)) + +function _local_vdim_construct_local_quadratures( + N, + mesh, + neighbors, + el, + center, + scale, + diam, + bdry_kdtree, + bdry_etype2qrule, + vol_etype2qrule, + bdry_qrule, + vol_qrule + ) + # construct the local region + Etype = first(element_types(mesh)) + el_neighs = neighbors[(Etype, el)] + + T = first(el_neighs)[1] + els_idxs = [i[2] for i in collect(el_neighs)] + els_list = mesh.etype2els[Etype][els_idxs] + + loc_bdry = boundarynd(T, els_idxs, mesh) + # TODO handle curved boundary of Γ?? + if N == 2 + bords = LagrangeElement{ReferenceHyperCube{N - 1}, 2, SVector{N, Float64}}[] + else + bords = LagrangeElement{ReferenceSimplex{N - 1}, 3, SVector{N, Float64}}[] + end + + for idxs in loc_bdry + vtxs = nodes(mesh)[idxs] + if N == 2 + bord = _newbord_line(vtxs) + else + bord = _newbord_tri(vtxs) + end + push!(bords, bord) + end + + # Check if we need to do near-singular layer potential evaluation + vertices = mesh.etype2els[Etype][el].vals[vertices_idxs(Etype)] + need_layer_corr = sum(inrangecount(bdry_kdtree, vertices, diam / 2)) > 0 + + # Now begin working in x̂ coordinates where x = scale * x̂ + + # build O(h) volume neighbors + Yvol = Quadrature(Float64, els_list, vol_etype2qrule, vol_qrule; center, scale) + Ybdry = Quadrature(Float64, bords, bdry_etype2qrule, bdry_qrule; center, scale) + + return Yvol, Ybdry, need_layer_corr +end + +function _local_vdim_auxiliary_quantities( + op::AbstractDifferentialOperator{N}, + center, + scale, + PFE_p, + PFE_P, + X, + μ, + Yvol, + Ybdry, + diam, + need_layer_corr + ) where {N} + # TODO handle derivative case + G = SingleLayerKernel(op) + dG = DoubleLayerKernel(op) + Xshift = [(q.coords - center) / scale for q in X] + Sop = IntegralOperator(G, Xshift, Ybdry) + Dop = IntegralOperator(dG, Xshift, Ybdry) + Vop = IntegralOperator(G, Xshift, Yvol) + Smat = assemble_matrix(Sop) + Dmat = assemble_matrix(Dop) + Vmat = assemble_matrix(Vop) + if need_layer_corr + μloc = _green_multiplier(:inside) + green_multiplier = fill(μloc, length(X)) + δS, δD = bdim_correction( + op, + Xshift, + Ybdry, + Smat, + Dmat; + green_multiplier, + maxdist = diam / scale, + derivative = false, + ) + + Smat += δS + Dmat += δD + end + + num_basis = length(PFE_P) + num_targets = length(X) + b = Matrix{Float64}(undef, length(Yvol), num_basis) + γ₁B = Matrix{Float64}(undef, length(Ybdry), num_basis) + γ₀B = Matrix{Float64}(undef, length(Ybdry), num_basis) + P = Matrix{Float64}(undef, length(X), num_basis) + grad = Array{Float64}(undef, num_basis, N, length(Ybdry)) + + for i in 1:length(Yvol) + ElementaryPDESolutions.fast_evaluate!(view(b, i, :), Yvol[i].coords, PFE_p) + end + for i in 1:length(X) + ElementaryPDESolutions.fast_evaluate!(view(P, i, :), Xshift[i], PFE_P) + end + for i in 1:length(Ybdry) + ElementaryPDESolutions.fast_evaluate_with_jacobian!( + view(γ₀B, i, :), + view(grad, :, :, i), + Ybdry[i].coords, + PFE_P, + ) + end + for i in 1:length(Ybdry) + for j in 1:num_basis + γ₁B[i, j] = 0 + for k in 1:N + γ₁B[i, j] += grad[j, k, i] * Ybdry[i].normal[k] #nrml_bdry_vec[i][k] + end + end + end + + Θ = zeros(eltype(Vop), num_targets, num_basis) + # Compute Θ <-- S * γ₁B - D * γ₀B - V * b + σ * B(x) using in-place matvec + for n in 1:num_basis + @views mul!(Θ[:, n], Smat, γ₁B[:, n]) + @views mul!(Θ[:, n], Dmat, γ₀B[:, n], -1, 1) + @views mul!(Θ[:, n], Vmat, b[:, n], -1, 1) + for i in 1:num_targets + Θ[i, n] += μ[i] * P[i, n] + end + end + return Θ, b +end + +function _lowfreq_vdim_auxiliary_quantities( + op::Laplace{2}, + op_lowfreq::Laplace{2}, + center, + scale, + num_basis, + PFE_p_lowfreq, + PFE_P_lowfreq, + multiindices, + multiindices_lowfreq, + monomials_indices, + monomials_indices_lowfreq, + X, + μ, + Yvol, + Ybdry, + diam, + need_layer_corr + ) + θ, b = _local_vdim_auxiliary_quantities( + op_lowfreq, + center, + scale, + PFE_p_lowfreq, + PFE_P_lowfreq, + X, + μ, + Yvol, + Ybdry, + diam, + need_layer_corr + ) + Xshift = [(q.coords - center) / scale for q in X] + num_targets = length(Xshift) + + Hmat = Matrix{eltype(θ)}(undef, num_targets, length(Yvol)) + for n in 1:num_targets + for j in 1:length(Yvol) + Hmat[n, j] = Yvol[j].weight + end + end + + R = Matrix{eltype(θ)}(undef, num_targets, num_basis) + for n in 1:num_basis + β = multiindices[n] + R[:, n] = θ[:, monomials_indices_lowfreq[β]] + R[:, n] += log(scale) * Hmat * b[:, monomials_indices_lowfreq[β]] + end + if isdefined(Main, :Infiltrator) + Main.infiltrate(@__MODULE__, Base.@locals, @__FILE__, @__LINE__) + end + return R +end + +function _lowfreq_vdim_auxiliary_quantities( + op::Helmholtz{2}, + op_lowfreq::Laplace{2}, + center, + scale, + num_basis, + PFE_p_lowfreq, + PFE_P_lowfreq, + multiindices, + multiindices_lowfreq, + monomials_indices, + monomials_indices_lowfreq, + X, + μ, + Yvol, + Ybdry, + diam, + need_layer_corr + ) + Xshift = [(q.coords - center) / scale for q in X] + # Laplace + θ, b = _local_vdim_auxiliary_quantities( + op_lowfreq, + center, + scale, + PFE_p_lowfreq, + PFE_P_lowfreq, + X, + μ, + Yvol, + Ybdry, + diam, + need_layer_corr, + ) + + num_targets = length(X) + + # Set up integral operators just to find the element type of the correction matrix R + G = SingleLayerKernel(op) + Vop = IntegralOperator(G, Xshift, Yvol) + R = zeros(eltype(Vop), num_targets, num_basis) + kr2 = (op.k * scale)^2 + γ = 0.5772156649015328606 + + Hmat = Matrix{ComplexF64}(undef, length(X), length(Yvol)) + #Hmat .= 0 + #@show scale * op.k + for i in 1:num_targets + for j in 1:length(Yvol) + z2 = + kr2 * ( + (Xshift[i][1] - Yvol[j].coords[1])^2 + + (Xshift[i][2] - Yvol[j].coords[2])^2 + ) + bessj0 = 0 + for k in 0:3 + bessj0 += (-1)^k * (1 / 4 * z2)^k / (factorial(k)^2) + end + Hmat[i, j] = + ( + (1 + 2 * im / pi * (γ + 1 / 2 * log(kr2 / 4))) * bessj0 + + 2 * im / pi * + (1 / 4 * z2 - 3 / 2 * (1 / 4 * z2)^2 / 4 + 11 / 3456 * z2^3) + ) * + Yvol[j].weight * + scale^2 + end + end + + for n in 1:num_basis + beta = multiindices[n] + beta10 = beta + Inti.MultiIndex((1, 0)) + beta01 = beta + Inti.MultiIndex((0, 1)) + beta20 = beta + Inti.MultiIndex((2, 0)) + beta02 = beta + Inti.MultiIndex((0, 2)) + for j in 1:num_targets + x1t = Xshift[j][1] + x2t = Xshift[j][2] + R[j, n] = + (1 - 1 / 4 * kr2 * (x1t^2 + x2t^2)) * θ[j, monomials_indices_lowfreq[beta]] + +1 / 2 * kr2 * x1t * θ[j, monomials_indices_lowfreq[beta10]] + +1 / 2 * kr2 * x2t * θ[j, monomials_indices_lowfreq[beta01]] + -1 / 4 * kr2 * θ[j, monomials_indices_lowfreq[beta20]] + -1 / 4 * kr2 * θ[j, monomials_indices_lowfreq[beta02]] + end + R[:, n] .*= 2 * im / pi * scale^2 + # Commenting this yields more accuracy + R[:, n] += Hmat * b[:, monomials_indices_lowfreq[beta]] + end + return R, b end function _vdim_auxiliary_quantities( @@ -308,6 +962,41 @@ function vdim_mesh_center(msh::AbstractMesh) end return xc / M end +""" + polynomial_solutions_local_vdim(op, order) + +For every monomial term `pₙ` of degree `order`, compute a polynomial `Pₙ` such +that `ℒ[Pₙ] = pₙ`, where `ℒ` is the differential operator `op`. +This function returns `{pₙ,Pₙ,γ₁Pₙ}`, where `γ₁Pₙ` is the generalized Neumann +trace of `Pₙ`. +""" +function polynomial_solutions_local_vdim(op::AbstractDifferentialOperator, order::Integer) + N = ambient_dimension(op) + # create empty arrays to store the monomials, solutions, and traces. For the + # neumann trace, we try to infer the concrete return type instead of simply + # having a vector of `Function`. + monomials = Vector{ElementaryPDESolutions.Polynomial{N, Float64}}() + poly_solutions = Vector{ElementaryPDESolutions.Polynomial{N, Float64}}() + multiindices = Vector{MultiIndex{N}}() + # iterate over N-tuples going from 0 to order + for I in Iterators.product(ntuple(i -> 0:order, N)...) + sum(I) > order && continue + # define the monomial basis functions, and the corresponding solutions. + # TODO: adapt this to vectorial case + p = ElementaryPDESolutions.Polynomial(I => 1 / factorial(MultiIndex(I))) + P = polynomial_solution(op, p) + push!(multiindices, MultiIndex(I)) + push!(monomials, p) + push!(poly_solutions, P) + end + monomials_indices = Dict(multiindices .=> 1:length(multiindices)) + + PFE_monomials = ElementaryPDESolutions.assemble_fastevaluator(monomials, Float64) + PFE_polysolutions = + ElementaryPDESolutions.assemble_fastevaluator(poly_solutions, Float64) + + return PFE_monomials, PFE_polysolutions, multiindices, monomials_indices +end """ polynomial_solutions_vdim(op, order[, center]) @@ -354,7 +1043,8 @@ function polynomial_solutions_vdim( return (q) -> f(coords(q) - center) end neumann_shift = map(neumann_traces) do f - return (q) -> f((coords = q.coords - center, normal = q.normal)) + # return (q) -> f((coords = q.coords - center, normal = q.normal)) + return (q) -> f(coords(q) - center, normal(q)) end return monomial_shift, dirchlet_shift, neumann_shift, multiindices # return monomials, dirchlet_traces, neumann_traces, multiindices @@ -387,7 +1077,8 @@ end function _normal_derivative(P::ElementaryPDESolutions.Polynomial{N, T}) where {N, T} ∇P = ElementaryPDESolutions.gradient(P) - return (q) -> dot(normal(q), ∇P(coords(q))) + # return (q) -> dot(normal(q), ∇P(coords(q))) + return (x, n) -> dot(n, ∇P(x)) end function (∇P::NTuple{N, <:ElementaryPDESolutions.Polynomial})(x) where {N} diff --git a/test/Project.toml b/test/Project.toml index 7559c3b7..7ec5c3f9 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,7 +1,9 @@ [deps] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" FMM2D = "2d63477d-9690-4b75-bcc1-c3461d43fecc" FMM3D = "1e13804c-f9b7-11ea-0ef0-29f3b1745df8" +FixedPolynomials = "3dd14ad9-0029-526e-86e9-8aa0bdd2ab0d" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb" diff --git a/test/boundary_test.jl b/test/boundary_test.jl new file mode 100644 index 00000000..0c47cc54 --- /dev/null +++ b/test/boundary_test.jl @@ -0,0 +1,84 @@ +using Test +using LinearAlgebra +using Inti +using Random +using Meshes +import 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) +# ell = collect(Ei[(E, 1)])[1] +# el = Inti.elements(msh, ell[1])[ell[2]] +# fig, _, _ = viz(msh; color = 0, showsegments = true,alpha=0.3) +# viz!(el; color = 0, showsegments = true,alpha=0.5) +# for (E, i) in els +# el = Inti.elements(msh, E)[i] +# viz!(el; showsegments = true, alpha=0.7) +# end +# viz!(bords;color=4,showsegments = false,segmentsize=5,segmentcolor=4) +# display(fig) +#end + +# el_in_set(el, set) = any(x->sort(x) == sort(el), set) + +I = 3 +test_els = union(copy(nei[(E, I)])) +els = Inti.elements(test_msh, E) +#test_els = union(copy(nei[(E,1)]), nei[(E,2)]) +#test_els = union(copy(nei[(E,1)]), nei[(E,2)], nei[(E,3)], nei[(E,4)]) +Inti.viz_elements(test_els, test_msh) + +components = Inti.connected_components(test_els, nei) + +test_els = copy(nei[(E, I)]) +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 + +els_idxs = [i[2] for i in collect(test_els)] +E = collect(test_els)[1][1] +els_list = test_msh.etype2els[E][els_idxs] +newquad = Inti.Quadrature(test_msh, els_list; qorder = 2) + +Inti.viz_elements_bords(nei, test_els, (E, I), bords, test_msh) diff --git a/test/convergence/vdim_potential_script.jl b/test/convergence/vdim_potential_script.jl index b7815324..a9043661 100644 --- a/test/convergence/vdim_potential_script.jl +++ b/test/convergence/vdim_potential_script.jl @@ -6,7 +6,7 @@ using Gmsh using LinearAlgebra using HMatrices using FMMLIB2D -using CairoMakie +using GLMakie function domain_and_mesh(; meshsize, meshorder = 1) Inti.clear_entities!() @@ -31,18 +31,14 @@ interpolation_order = 4 VR_qorder = Inti.Triangle_VR_interpolation_order_to_quadrature_order(interpolation_order) bdry_qorder = 2 * VR_qorder -tmesh = @elapsed begin - Ω, msh = domain_and_mesh(; meshsize) -end -@info "Mesh generation time: $tmesh" - +Ω, msh = domain_and_mesh(; meshsize) Γ = Inti.external_boundary(Ω) #Γₕ = msh[Γ] #Ωₕ = msh[Ω] ψ = (t) -> [cos(2 * π * t), sin(2 * π * t)] θ = 3 # smoothness order of curved elements -crvmsh = Inti.curve_mesh(msh, ψ, θ, 500 * Int(1 / meshsize)) +crvmsh = Inti.curve_mesh(msh, ψ, θ; patch_sample_num = 500 * Int(1 / meshsize)) Ωₕ = view(crvmsh, Ω) Γₕ = view(crvmsh, Γ) # @@ -54,7 +50,6 @@ tquad = @elapsed begin Q = Inti.VioreanuRokhlin(; domain = :triangle, order = VR_qorder) dict = Dict(E => Q for E in Inti.element_types(Ωₕ)) Ωₕ_quad = Inti.Quadrature(Ωₕ, dict) - # Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorders[1]) Γₕ_quad = Inti.Quadrature(Γₕ; qorder = bdry_qorder) end @info "Quadrature generation time: $tquad" diff --git a/test/lvdim_test.jl b/test/lvdim_test.jl new file mode 100644 index 00000000..c3c0a339 --- /dev/null +++ b/test/lvdim_test.jl @@ -0,0 +1,121 @@ +# # Testing local vdim + +using DynamicPolynomials +using FixedPolynomials +using Inti +using StaticArrays +using Gmsh +using LinearAlgebra +using HMatrices +using FMMLIB2D +using Meshes + +#meshsize = 0.001/8 +#meshsize = 0.125/8 +meshsize = 0.000125 +interpolation_order = 4 +VR_qorder = Inti.Triangle_VR_interpolation_order_to_quadrature_order(4) +bdry_qorder = 2 * VR_qorder + +function gmsh_disk(; name, meshsize, order = 1, center = (0, 0), paxis = (2, 1)) + return try + gmsh.initialize() + gmsh.option.setNumber("General.Terminal", 0) + gmsh.model.add("circle-mesh") + gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) + gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) + gmsh.model.occ.addDisk(center[1], center[2], 0, paxis[1], paxis[2]) + gmsh.model.occ.synchronize() + gmsh.model.mesh.generate(2) + gmsh.model.mesh.setOrder(order) + gmsh.write(name) + finally + gmsh.finalize() + end +end + +name = joinpath(@__DIR__, "disk.msh") +gmsh_disk(; meshsize, order = 2, name, paxis = (meshsize * 20, meshsize * 10)) + +Inti.clear_entities!() # empty the entity cache +msh = Inti.import_mesh(name; dim = 2) +Ω = Inti.Domain(e -> Inti.geometric_dimension(e) == 2, Inti.entities(msh)) +Γ = Inti.boundary(Ω) + +Ωₕ = msh[Ω] +Γₕ = msh[Γ] +Ωₕ_Sub = view(msh, Ω) +Γₕ_Sub = view(msh, Γ) + +tquad = @elapsed begin + # Use VDIM with the Vioreanu-Rokhlin quadrature rule for Ωₕ + Q = Inti.VioreanuRokhlin(; domain = :triangle, order = VR_qorder) + dict = Dict(E => Q for E in Inti.element_types(Ωₕ)) + Ωₕ_quad = Inti.Quadrature(Ωₕ, dict) + Ωₕ_Sub_quad = Inti.Quadrature(Ωₕ_Sub, dict) + # Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorders[1]) + Γₕ_quad = Inti.Quadrature(Γₕ; qorder = bdry_qorder) + Γₕ_Sub_quad = Inti.Quadrature(Γₕ_Sub; qorder = bdry_qorder) +end +@info "Quadrature generation time: $tquad" + +k0 = 1 +k = 0 +θ = (cos(π / 3), sin(π / 3)) +#u = (x) -> exp(im * k0 * dot(x, θ)) +#du = (x,n) -> im * k0 * dot(θ, n) * exp(im * k0 * dot(x, θ)) +u = (x) -> cos(k0 * dot(x, θ)) +du = (x, n) -> -k0 * dot(θ, n) * sin(k0 * dot(x, θ)) +f = (x) -> (k^2 - k0^2) * u(x) + +#s = 4 +#u = (x) -> 1 / (k^2 - k0^2) * exp(im * k0 * dot(x, θ)) + 1 / (k^2 - 4 * s) * exp(-s * norm(x)^2) +#du = (x, n) -> im * k0 * dot(θ, n) / (k^2 - k0^2) * exp(im * k0 * dot(x, θ)) - 2 * s / (k^2 - 4 * s) * dot(x, n) * exp(-s * norm(x)^2) +#f = (x) -> exp(im * k0 * dot(x, θ)) + 1 / (k^2 - 4 * s) * (4 * s^2 * norm(x)^2 - 4 * s + k^2) * exp(-s * norm(x)^2) + +u_d = map(q -> u(q.coords), Ωₕ_quad) +u_b = map(q -> u(q.coords), Γₕ_quad) +du_b = map(q -> du(q.coords, q.normal), Γₕ_quad) +f_d = map(q -> f(q.coords), Ωₕ_quad) + +op = k == 0 ? Inti.Laplace(; dim = 2) : Inti.Helmholtz(; dim = 2, k) + +## Boundary operators +tbnd = @elapsed begin + S_b2d, D_b2d = Inti.single_double_layer(; + op, + target = Ωₕ_quad, + source = Γₕ_quad, + compression = (method = :fmm, tol = 1.0e-14), + correction = (method = :dim, maxdist = 5 * meshsize, target_location = :inside), + ) +end +@info "Boundary operators time: $tbnd" + +## Volume potentials +#tvol = @elapsed begin +V_d2d = Inti.volume_potential(; + op, + target = Ωₕ_quad, + source = Ωₕ_quad, + compression = (method = :fmm, tol = 1.0e-14), + correction = ( + method = :ldim, + mesh = Ωₕ, + interpolation_order, + quadrature_order = VR_qorder, + bdry_nodes = Γₕ.nodes, + maxdist = 5 * meshsize, + meshsize = meshsize, + ), +) +#end +#@info "Volume potential time: $tvol" + +vref = -u_d - D_b2d * u_b + S_b2d * du_b +vapprox = V_d2d * f_d +er = vref - vapprox + +ndofs = length(er) + +@show ndofs, meshsize, norm(er, Inf) diff --git a/test/lvdim_test_3d.jl b/test/lvdim_test_3d.jl new file mode 100644 index 00000000..6ffdbd19 --- /dev/null +++ b/test/lvdim_test_3d.jl @@ -0,0 +1,112 @@ +# # High-order convergence of vdim + +using DynamicPolynomials +using FixedPolynomials +using Inti +using Meshes +using StaticArrays +using Gmsh +using LinearAlgebra +using HMatrices +using FMM3D +using GLMakie + +meshsize = 0.05 +interpolation_order = 4 +VR_qorder = Inti.Tetrahedron_VR_interpolation_order_to_quadrature_order(interpolation_order) +bdry_qorder = 2 * VR_qorder + +function gmsh_sphere(; order = 1, name, meshsize) + return try + gmsh.initialize() + gmsh.option.setNumber("General.Terminal", 0) + gmsh.model.add("sphere-mesh") + gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) + gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) + gmsh.model.occ.addSphere(0, 0, 0, 1) + gmsh.model.occ.synchronize() + gmsh.model.mesh.generate() + gmsh.model.mesh.setOrder(order) + gmsh.write(name) + finally + gmsh.finalize() + end +end + +name = joinpath(@__DIR__, "sphere.msh") +gmsh_sphere(; order = 1, name, meshsize) + +Inti.clear_entities!() # empty the entity cache +msh = Inti.import_mesh(name; dim = 3) +Ω = Inti.Domain(e -> Inti.geometric_dimension(e) == 3, Inti.entities(msh)) +Γ = Inti.boundary(Ω) + +Ωₕ = msh[Ω] +Γₕ = msh[Γ] +Ωₕ_Sub = view(msh, Ω) +Γₕ_Sub = view(msh, Γ) + +tquad = @elapsed begin + # Use VDIM with the Vioreanu-Rokhlin quadrature rule for Ωₕ + Q = Inti.VioreanuRokhlin(; domain = :tetrahedron, order = VR_qorder) + dict = Dict(E => Q for E in Inti.element_types(Ωₕ)) + Ωₕ_quad = Inti.Quadrature(Ωₕ, dict) + # Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorders[1]) + Γₕ_quad = Inti.Quadrature(Γₕ; qorder = bdry_qorder) +end +@info "Quadrature generation time: $tquad" + +k0 = π +k = 0 +θ = (sin(π / 3) * cos(π / 3), sin(π / 3) * sin(π / 3), cos(π / 3)) +#u = (x) -> exp(im * k0 * dot(x, θ)) +#du = (x,n) -> im * k0 * dot(θ, n) * exp(im * k0 * dot(x, θ)) +u = (x) -> cos(k0 * dot(x, θ)) +du = (x, n) -> -k0 * dot(θ, n) * sin(k0 * dot(x, θ)) +f = (x) -> (k^2 - k0^2) * u(x) + +u_d = map(q -> u(q.coords), Ωₕ_quad) +u_b = map(q -> u(q.coords), Γₕ_quad) +du_b = map(q -> du(q.coords, q.normal), Γₕ_quad) +f_d = map(q -> f(q.coords), Ωₕ_quad) + +pde = k == 0 ? Inti.Laplace(; dim = 3) : Inti.Helmholtz(; dim = 3, k) + +## Boundary operators +tbnd = @elapsed begin + S_b2d, D_b2d = Inti.single_double_layer(; + pde, + target = Ωₕ_quad, + source = Γₕ_quad, + compression = (method = :fmm, tol = 1.0e-8), + correction = (method = :dim, maxdist = 5 * meshsize, target_location = :inside), + ) +end +@info "Boundary operators time: $tbnd" + +## Volume potentials +tvol = @elapsed begin + V_d2d = Inti.volume_potential(; + pde, + target = Ωₕ_quad, + source = Ωₕ_quad, + compression = (method = :fmm, tol = 1.0e-8), + correction = ( + method = :ldim, + interpolation_order, + quadrature_order = VR_qorder, + mesh = Ωₕ, + bdry_nodes = Γₕ.nodes, + maxdist = 5 * meshsize, + ), + ) +end +@info "Volume potential time: $tvol" + +vref = -u_d - D_b2d * u_b + S_b2d * du_b +vapprox = V_d2d * f_d +er = vref - vapprox + +ndofs = length(er) + +@show ndofs, meshsize, norm(er, Inf) diff --git a/test/poly_test.jl b/test/poly_test.jl new file mode 100644 index 00000000..671e1faa --- /dev/null +++ b/test/poly_test.jl @@ -0,0 +1,77 @@ +using Inti +using StaticArrays +using LinearAlgebra +import DynamicPolynomials: @polyvar +using FixedPolynomials + +#Run with npts = 1000000 +function poly_test(npts) + #npts = 10000 + pde = Inti.Laplace(; dim = 2) + interpolation_order = 4 + p, P, γ₁P, multiindices = Inti.polynomial_solutions_vdim(pde, interpolation_order) + + @polyvar x y + + N = 2 + PolArray = Array{FixedPolynomials.Polynomial{Float64}}(undef, length(P)) + + tsetup = @elapsed begin + for (polind, ElemPolySolsPols) in enumerate(P) + pp = ElemPolySolsPols.f.order2coeff + exp_data = Matrix{Int64}(undef, N, length(pp)) + coeff_data = Vector{Float64}(undef, length(pp)) + for (i, pol) in enumerate(pp) + exp_data[:, i] = [q for q in pol[1]] + coeff_data[i] = pol[2] + end + PolArray[polind] = + FixedPolynomials.Polynomial(exp_data, coeff_data, [:x, :y]) + end + PolSystem = System(PolArray) + pts = Vector{Vector{Float64}}(undef, npts) + pts[1] = [0, 0] + cfg = JacobianConfig(PolSystem, pts[1]) + end + @info "FixedPolynomials.jl setup time: $tsetup" + + pts = Vector{Vector{Float64}}(undef, npts) + for i in 1:npts + pts[i] = rand(2) + end + cfg = JacobianConfig(PolSystem, pts[1]) + res1 = Matrix{Float64}(undef, length(PolArray), length(pts)) + res2 = Matrix{Float64}(undef, length(PolArray), length(pts)) + res3 = Vector{MVector{length(PolArray), Float64}}(undef, npts) + + cfg = JacobianConfig(PolSystem, pts[1]) + u = Vector{Float64}(undef, length(PolArray)) + tfixed = @elapsed begin + for i in 1:npts + evaluate!(view(res1, :, i), PolSystem, pts[i], cfg) + end + end + @info "FixedPolynomials.jl time: $tfixed" + evaluator = (xx) -> evaluate(PolSystem, xx, cfg) + tbroadcast = @elapsed begin + res3 .= evaluator.(pts) + end + @info "FixedPolynomials.jl w/ broadcast time: $tbroadcast" + tregular = @elapsed begin + for i in 1:length(P) + res2[i, :] .= P[i].f.(pts) + end + end + @info "ElementaryPDESolutions.jl time: $tregular" + + # Evaluate Jacobian + u = Matrix{Float64}(undef, length(P), npts) + U = Array{Float64}(undef, length(P), 2, npts) + tjacob = @elapsed begin + for i in 1:npts + evaluate_and_jacobian!(view(u, :, i), view(U, :, :, i), PolSystem, pts[i], cfg) + end + end + @info "FixedPolynomials.jl Jacobian+Eval time: $tjacob" + return res1, res2, res3, tfixed, tregular, tbroadcast, tjacob +end diff --git a/test/test_utils.jl b/test/test_utils.jl index 3dfe77fc..80dd3305 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(String[], false) 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()