From 67f84f7872621e759d648aee95291e3f78d6e08a Mon Sep 17 00:00:00 2001 From: maltezfaria Date: Mon, 19 Jan 2026 10:00:35 +0100 Subject: [PATCH 01/20] simplify `.gitignore` - use `.git/info/exclude` for local ignores instead --- .gitignore | 18 ++---------------- 1 file changed, 2 insertions(+), 16 deletions(-) diff --git a/.gitignore b/.gitignore index 8ef96d1e..56a8913e 100644 --- a/.gitignore +++ b/.gitignore @@ -1,19 +1,5 @@ -*.jl.*.cov -*.jl.cov -*.jl.mem -*.msh -*.pos -/Manifest.toml +Manifest.toml /docs/Manifest.toml -/docs/build/ /test/Manifest.toml -.vscode -/docs/src/pluto-examples/*.md -/sandbox/* -/talks/* +/docs/build/ /test/vtk_output/* - -.markdownlint.json -Inti.code-workspace -docs/src/examples/heat.gif -docs/src/.DS_Store From 29993a377d3e499a79ea132a37d4fd5a9959612b Mon Sep 17 00:00:00 2001 From: maltezfaria Date: Mon, 19 Jan 2026 10:04:58 +0100 Subject: [PATCH 02/20] cleanup --- src/blockarray.jl | 2 +- test/Project.toml | 4 + test/curved_test_2d.jl | 224 ----------------------------------------- test/curved_test_3d.jl | 77 -------------- 4 files changed, 5 insertions(+), 302 deletions(-) delete mode 100644 test/curved_test_2d.jl delete mode 100644 test/curved_test_3d.jl diff --git a/src/blockarray.jl b/src/blockarray.jl index f080dac0..f70aa50d 100644 --- a/src/blockarray.jl +++ b/src/blockarray.jl @@ -84,7 +84,7 @@ function BlockArray(src::AbstractArray{T}) where {T <: SArray} return dest end -function Array(src::BlockArray{T, N}) where {T, N} +function Base.Array(src::BlockArray{T, N}) where {T, N} dest = Array{T, N}(undef, size(src)) dest .= src return dest diff --git a/test/Project.toml b/test/Project.toml index 7559c3b7..1c9efc78 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,5 +1,6 @@ [deps] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +ElementaryPDESolutions = "88a69b33-da0f-4502-8c8c-d680cf4d883b" FMM2D = "2d63477d-9690-4b75-bcc1-c3461d43fecc" FMM3D = "1e13804c-f9b7-11ea-0ef0-29f3b1745df8" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" @@ -16,3 +17,6 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" + +[sources] +Inti = {path = ".."} diff --git a/test/curved_test_2d.jl b/test/curved_test_2d.jl deleted file mode 100644 index ff42d19e..00000000 --- a/test/curved_test_2d.jl +++ /dev/null @@ -1,224 +0,0 @@ -using Inti -using Meshes -using StaticArrays -using GLMakie -using Gmsh -using LinearAlgebra -using ForwardDiff -using Test - -@testset begin - function domain_and_mesh(; meshsize, meshorder = 1) - Inti.clear_entities!() - gmsh.initialize() - gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) - gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) - gmsh.model.occ.addDisk(0, 0, 0, 1, 1) - gmsh.model.occ.synchronize() - gmsh.model.mesh.generate(2) - gmsh.model.mesh.setOrder(meshorder) - msh = Inti.import_mesh(; dim = 2) - Ω = Inti.Domain(Inti.entities(msh)) do ent - return Inti.geometric_dimension(ent) == 2 - end - gmsh.finalize() - return Ω, msh - end - - meshsize = 0.1 - - tmesh = @elapsed begin - Ω, msh = domain_and_mesh(; meshsize) - end - @info "Mesh generation time: $tmesh" - - Γ = Inti.external_boundary(Ω) - Ωₕ = view(msh, Ω) - Γₕ = view(msh, Γ) - - ψ = (t) -> [cos(2 * π * t), sin(2 * π * t)] - θ = 6 # smoothness order of curved elements - crvmsh = Inti.curve_mesh(msh, ψ, θ) - - Γₕ = crvmsh[Γ] - Ωₕ = crvmsh[Ω] - - qorder = 2 - Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) - Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) - @test isapprox(Inti.integrate(x -> 1, Ωₕ_quad), π, rtol = 1.0e-7) - @test isapprox(Inti.integrate(q -> q.coords[1]^4, Ωₕ_quad), π / 8, rtol = 1.0e-5) - - qorder = 5 - Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) - Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) - @test isapprox(Inti.integrate(x -> 1, Ωₕ_quad), π, rtol = 1.0e-11) - @test isapprox(Inti.integrate(q -> q.coords[1]^4, Ωₕ_quad), π / 8, rtol = 1.0e-10) - - qorder = 8 - Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) - Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) - @test isapprox(Inti.integrate(x -> 1, Ωₕ_quad), π, rtol = 1.0e-14) - @test isapprox(Inti.integrate(q -> q.coords[1]^4, Ωₕ_quad), π / 8, rtol = 1.0e-14) - - Fvol = (x) -> x[2]^2 - 2 * x[2] * x[1]^3 - F = (x) -> [x[1] * x[2]^2, x[1]^3 * x[2]^2] - #Fvol = (x) -> 1.0 - #F = (x) -> [1/2*x[1], 1/2*x[2]] - divvol = Inti.integrate(q -> Fvol(q.coords), Ωₕ_quad) - divline = Inti.integrate(q -> dot(F(q.coords), q.normal), Γₕ_quad) - @test isapprox(divline, divvol, rtol = 1.0e-13) -end - -@testset begin - meshsize = 2π / 4 / 8 - Inti.clear_entities!() - - gmsh.initialize() - gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) - gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) - - # Two kites - kite = Inti.gmsh_curve(0, 1; meshsize) do s - return SVector(0.25, 0.0) + - SVector(cos(2π * s) + 0.65 * cos(4π * s[1]) - 0.65, 1.5 * sin(2π * s)) - end - cl = gmsh.model.occ.addCurveLoop([kite]) - surf = gmsh.model.occ.addPlaneSurface([cl]) - kite_trans = Inti.gmsh_curve(0, 1; meshsize) do s - return SVector(4.5, 0.0) + - SVector(cos(2π * s) + 0.65 * cos(4π * s[1]) - 0.65, 1.5 * sin(2π * s)) - end - cl_trans = gmsh.model.occ.addCurveLoop([kite_trans]) - surf_trans = gmsh.model.occ.addPlaneSurface([cl_trans]) - gmsh.model.occ.synchronize() - - # Add tags for stable identification of the entities - gmsh.model.addPhysicalGroup(2, [surf], -1, "c1") - gmsh.model.addPhysicalGroup(2, [surf_trans], -1, "c2") - - gmsh.model.mesh.generate(2) - msh = Inti.import_mesh(; dim = 2) - Ω = Inti.Domain(Inti.entities(msh)) do ent - return Inti.geometric_dimension(ent) == 2 - end - gmsh.finalize() - - Γ = Inti.external_boundary(Ω) - Ωₕ = view(msh, Ω) - Γₕ = view(msh, Γ) - - # Two kites - ψ₁ = (t) -> [0.25 + cos(2π * t) + 0.65 * cos(4π * t) - 0.65, 1.5 * sin(2π * t)] - ψ₂ = (t) -> [4.5 + cos(2π * t) + 0.65 * cos(4π * t) - 0.65, 1.5 * sin(2π * t)] - entity_parametrizations = Dict{Inti.EntityKey, Function}() - for e in Inti.entities(Ω) - l = Inti.labels(e) - if "c1" in l - entity_parametrizations[e] = ψ₁ - elseif "c2" in l - entity_parametrizations[e] = ψ₂ - elseif "c3" in l - entity_parametrizations[e] = ψ₃ - end - end - θ = 6 # smoothness order of curved elements - crvmsh = Inti.curve_mesh(msh, entity_parametrizations, θ) - - Γₕ = crvmsh[Γ] - Ωₕ = crvmsh[Ω] - - qorder = 8 - Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) - Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) - - Fvol = (x) -> x[2]^2 - 2 * x[2] * x[1]^3 - F = (x) -> [x[1] * x[2]^2, x[1]^3 * x[2]^2] - #Fvol = (x) -> 1.0 - #F = (x) -> [1/2*x[1], 1/2*x[2]] - divvol = Inti.integrate(q -> Fvol(q.coords), Ωₕ_quad) - divline = Inti.integrate(q -> dot(F(q.coords), q.normal), Γₕ_quad) - @test isapprox(divline, divvol, rtol = 1.0e-11) -end - -@testset begin - meshsize = 0.075 - Inti.clear_entities!() - - gmsh.initialize() - gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) - gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) - - # Three circles - c1 = gmsh.model.occ.addDisk(0, 0, 0, 1, 1) - c2 = gmsh.model.occ.addDisk(0, 3.0, 0, 1, 1) - c3 = gmsh.model.occ.addDisk(0, 8.0, 0, 2, 2) - gmsh.model.occ.synchronize() - - # Add tags for stable identification of the entities - gmsh.model.addPhysicalGroup(2, [c1], -1, "c1") - gmsh.model.addPhysicalGroup(2, [c2], -1, "c2") - gmsh.model.addPhysicalGroup(2, [c3], -1, "c3") - - gmsh.model.mesh.generate(2) - msh = Inti.import_mesh(; dim = 2) - Ω = Inti.Domain(Inti.entities(msh)) do ent - return Inti.geometric_dimension(ent) == 2 - end - gmsh.finalize() - - Γ = Inti.external_boundary(Ω) - Ωₕ = view(msh, Ω) - Γₕ = view(msh, Γ) - - # Three circles - ψ₁ = (t) -> [cos(2 * π * t), sin(2 * π * t)] - ψ₂ = (t) -> [cos(2 * π * t), 3.0 + sin(2 * π * t)] - ψ₃ = (t) -> [2 * cos(2 * π * t), 8.0 + 2 * sin(2 * π * t)] - entity_parametrizations = Dict{Inti.EntityKey, Function}() - for e in Inti.entities(Ω) - l = Inti.labels(e) - if "c1" in l - entity_parametrizations[e] = ψ₁ - elseif "c2" in l - entity_parametrizations[e] = ψ₂ - elseif "c3" in l - entity_parametrizations[e] = ψ₃ - end - end - - θ = 6 # smoothness order of curved elements - crvmsh = Inti.curve_mesh(msh, entity_parametrizations, θ) - - Γₕ = crvmsh[Γ] - Ωₕ = crvmsh[Ω] - - qorder = 2 - Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) - Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) - @test isapprox(Inti.integrate(x -> 1, Ωₕ_quad), 6π, rtol = 1.0e-6) - - qorder = 5 - Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) - Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) - @test isapprox(Inti.integrate(x -> 1, Ωₕ_quad), 6π, rtol = 1.0e-11) - - qorder = 8 - Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) - Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) - @test isapprox(Inti.integrate(x -> 1, Ωₕ_quad), 6π, rtol = 1.0e-14) - - # The third circle has larger volume than others, so check that we can pick sub-domains correctly - Ω_sub = Inti.Domain(collect(keys(Ω))[3]) - Ωₕ_sub = crvmsh[Ω_sub] - Ωₕ_sub_quad = Inti.Quadrature(Ωₕ_sub; qorder = qorder) - @test isapprox(Inti.integrate(x -> 1, Ωₕ_sub_quad), 4π, rtol = 1.0e-14) - - Fvol = (x) -> x[2]^2 - 2 * x[2] * x[1]^3 - F = (x) -> [x[1] * x[2]^2, x[1]^3 * x[2]^2] - #Fvol = (x) -> 1.0 - #F = (x) -> [1/2*x[1], 1/2*x[2]] - divvol = Inti.integrate(q -> Fvol(q.coords), Ωₕ_quad) - divline = Inti.integrate(q -> dot(F(q.coords), q.normal), Γₕ_quad) - @test isapprox(divline, divvol, rtol = 1.0e-13) -end diff --git a/test/curved_test_3d.jl b/test/curved_test_3d.jl deleted file mode 100644 index e4e89f7c..00000000 --- a/test/curved_test_3d.jl +++ /dev/null @@ -1,77 +0,0 @@ -using Test -using Inti -using FMM3D -using Gmsh -using LinearAlgebra -using Random -using StaticArrays -using ForwardDiff - -include("test_utils.jl") - -@testset begin - # create a boundary and area meshes and quadrature only once - Inti.clear_entities!() - meshsize = 0.1 - r1 = 1.0 - r2 = 0.5 - Ω, msh = gmsh_torus(; center = [0.0, 0.0, 0.0], r1 = r1, r2 = r2, meshsize = meshsize) - Γ = Inti.external_boundary(Ω) - Γ_msh = view(msh, Γ) - - function face_element_on_torus(nodelist, R, r) - return all( - [ - (sqrt(node[1]^2 + node[2]^2) - R^2)^2 + node[3]^2 ≈ r^2 for node in nodelist - ] - ) - end - face_element_on_curved_surface = (nodelist) -> face_element_on_torus(nodelist, r1, r2) - - function ψ(v::AbstractVector) - return [ - (r1 + r2 * sin(v[1])) * cos(v[2]), - (r1 + r2 * sin(v[1])) * sin(v[2]), - r2 * cos(v[1]), - ] - end - θ = 6 # smoothness order of curved elements - crvmsh = Inti.curve_mesh( - msh, - ψ, - θ; - face_element_on_curved_surface = face_element_on_curved_surface, - ) - - Γₕ = crvmsh[Γ] - Ωₕ = crvmsh[Ω] - - truevol = 2 * π^2 * r2^2 * r1 - truesfcarea = 4 * π^2 * r1 * r2 - - qorder = 2 - Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) - Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) - @test isapprox(Inti.integrate(x -> 1, Ωₕ_quad), truevol, rtol = 1.0e-5) - @test isapprox(Inti.integrate(x -> 1, Γₕ_quad), truesfcarea, rtol = 1.0e-5) - - qorder = 5 - Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) - Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) - @test isapprox(Inti.integrate(x -> 1, Ωₕ_quad), truevol, rtol = 1.0e-7) - @test isapprox(Inti.integrate(x -> 1, Γₕ_quad), truesfcarea, rtol = 1.0e-7) - - qorder = 7 - Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) - Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) - @test isapprox(Inti.integrate(x -> 1, Ωₕ_quad), truevol, rtol = 1.0e-11) - @test isapprox(Inti.integrate(x -> 1, Γₕ_quad), truesfcarea, rtol = 1.0e-14) - - divF = (x) -> x[3] + x[3]^2 + x[2]^3 - F = (x) -> [x[1] * x[3], x[2] * x[3]^2, x[2]^3 * x[3]] - #divF = (x) -> 1.0 - #F = (x) -> 1/3*[x[1], x[2], x[3]] - divtest_vol = Inti.integrate(q -> divF(q.coords), Ωₕ_quad) - divtest_sfc = Inti.integrate(q -> dot(F(q.coords), q.normal), Γₕ_quad) - @test isapprox(divtest_vol, divtest_sfc, rtol = 1.0e-9) -end From 2bb11d2ac8fe79450006cff8fbbc1bf34e5d1aa2 Mon Sep 17 00:00:00 2001 From: maltezfaria Date: Mon, 19 Jan 2026 10:05:48 +0100 Subject: [PATCH 03/20] avoid use of `@eval` to help static analysis tools --- src/entities.jl | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/src/entities.jl b/src/entities.jl index 45e7496c..ede31e9e 100644 --- a/src/entities.jl +++ b/src/entities.jl @@ -20,16 +20,12 @@ Base.hash(ent::EntityKey, h::UInt) = hash((ent.dim, abs(ent.tag)), h) Base.:(==)(e1::EntityKey, e2::EntityKey) = e1.dim == e2.dim && abs(e1.tag) == abs(e2.tag) # defer some functions on EntityKey to the corresponding GeometricEntity -for f in ( - :labels, - :boundary, - :pushforward, - :ambient_dimension, - :hasparametrization, - :parametrization, - ) - @eval $f(k::EntityKey) = $f(global_get_entity(k)) -end +labels(k::EntityKey) = labels(global_get_entity(k)) +boundary(k::EntityKey) = boundary(global_get_entity(k)) +pushforward(k::EntityKey) = pushforward(global_get_entity(k)) +ambient_dimension(k::EntityKey) = ambient_dimension(global_get_entity(k)) +hasparametrization(k::EntityKey) = hasparametrization(global_get_entity(k)) +parametrization(k::EntityKey) = parametrization(global_get_entity(k)) function (k::EntityKey)(x) hasparametrization(k) || error("$k has no parametrization") From dee3b8d522a80cff905af877a59d7dea8bdc0147 Mon Sep 17 00:00:00 2001 From: maltezfaria Date: Mon, 19 Jan 2026 10:09:26 +0100 Subject: [PATCH 04/20] `range_dimension` of a kernel --- src/kernels.jl | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/kernels.jl b/src/kernels.jl index c0796ac4..bc66b26e 100644 --- a/src/kernels.jl +++ b/src/kernels.jl @@ -33,6 +33,17 @@ abstract type AbstractDifferentialOperator{N} end ambient_dimension(::AbstractDifferentialOperator{N}) where {N} = N +function range_dimension(op::AbstractDifferentialOperator) + T = default_density_eltype(op) + if T <: Number + return 1 + elseif hasmethod(length, Tuple{Type{T}}) + return length(T) + else + error("default_density_eltype($(typeof(op))) = $T does not define length(::Type{$T}); cannot determine range dimension") + end +end + # convenient constructor for e.g. SingleLayerKernel(op,Float64) or DoubleLayerKernel(op,ComplexF64) function (::Type{K})( op::Op, From 5607c919866cdc34edbfa8c46b0e30d9bcb669fd Mon Sep 17 00:00:00 2001 From: maltezfaria Date: Mon, 19 Jan 2026 11:32:59 +0100 Subject: [PATCH 05/20] `vdim` for vector valued (maybe working?) --- src/vdim.jl | 414 +++++++++++++++++++--------------------------------- 1 file changed, 149 insertions(+), 265 deletions(-) diff --git a/src/vdim.jl b/src/vdim.jl index 88243f30..ef699851 100644 --- a/src/vdim.jl +++ b/src/vdim.jl @@ -28,9 +28,9 @@ See [anderson2024fast](@cite) for more details on the method. and rescaled to each element. """ function vdim_correction( - op, + op::AbstractDifferentialOperator{N}, target, - source::Quadrature, + source::Quadrature{N}, boundary::Quadrature, Sop, Dop, @@ -38,81 +38,79 @@ function vdim_correction( green_multiplier::Vector{<:Real}, interpolation_order = nothing, maxdist = Inf, - center = nothing, - shift::Val{SHIFT} = Val(false), - ) where {SHIFT} + ) where {N} # variables for debugging the condition properties of the method vander_cond = vander_norm = rhs_norm = res_norm = shift_norm = -Inf T = eltype(Vop) + # determine type for dense matrices + Dense = T <: SMatrix ? BlockArray : Array @assert eltype(Dop) == eltype(Sop) == T "eltype of Sop, Dop, and Vop must match" # 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) + num_target, num_source = length(target), length(source) # a reasonable interpolation_order if not provided isnothing(interpolation_order) && (interpolation_order = maximum(order, values(source.etype2qrule))) # by default basis centered at origin - center = isnothing(center) ? zero(SVector{N, Float64}) : center - p, P, γ₁P, multiindices = polynomial_solutions_vdim(op, interpolation_order, center) + p, P, γ₁P = polynomial_solutions_vdim(op, interpolation_order, T) dict_near = etype_to_nearest_points(target, source; maxdist) - R = _vdim_auxiliary_quantities( - p, - P, - γ₁P, - target, - source, - boundary, - green_multiplier, - Sop, - Dop, - Vop, - ) + num_basis = length(p) + b = Dense{T}(undef, length(source), num_basis) + γ₀B = Dense{T}(undef, length(boundary), num_basis) + γ₁B = Dense{T}(undef, length(boundary), num_basis) + b = [f(q) for q in source, f in p] + for k in 1:num_basis, j in 1:length(boundary) + γ₀B[j, k] = P[k](boundary[j]) + γ₁B[j, k] = γ₁P[k](boundary[j]) + end + T = eltype(γ₀B) + Θ = zeros(T, num_target, num_basis) + # Compute Θ <-- S * γ₁B - D * γ₀B - V * b + σ * B(x) using in-place matvec + for n in 1:num_basis + @views mul!(Θ[:, n], Sop, γ₁B[:, n]) + @views mul!(Θ[:, n], Dop, γ₀B[:, n], -1, 1) + @views mul!(Θ[:, n], Vop, b[:, n], -1, 1) + for i in 1:num_target + Θ[i, n] += green_multiplier[i] * P[n](target[i]) + end + end # compute sparse correction Is = Int[] Js = Int[] Vs = eltype(Vop)[] for (E, qtags) in source.etype2qtags - els = elements(source.mesh, E) + # els = elements(source.mesh, E) near_list = dict_near[E] nq, ne = size(qtags) @assert length(near_list) == ne + # preallocate local arrays to store interpolant values and weights + # Similar to bdim, we need to handle the case where T is an SMatrix + # by converting between Matrix{<:SMatrix} and Matrix{<:Number} formats + L_arr = Dense{T}(undef, num_basis, nq) + b_arr = Dense{T}(undef, num_basis, 1) + wei_arr = Dense{T}(undef, nq, 1) + Ldata, bdata, weidata = parent(L_arr)::Matrix, parent(b_arr)::Matrix, parent(wei_arr)::Matrix for n in 1:ne # indices of nodes in element `n` isempty(near_list[n]) && continue jglob = @view qtags[:, n] - # compute translation and scaling - c, r = translation_and_scaling(els[n]) - if SHIFT - iszero(center) || error("SHIFT is not implemented for non-zero center") - L̃ = [f((q.coords - c) / r) for f in p, q in view(source, jglob)] - S = change_of_basis(multiindices, p, c, r) - F = svd(L̃) - @debug (vander_cond = max(vander_cond, cond(L̃))) maxlog = 0 - @debug (shift_norm = max(shift_norm, norm(S))) maxlog = 0 - @debug (vander_norm = max(vander_norm, norm(L̃))) maxlog = 0 - else - L = [f(q.coords) for f in p, q in view(source, jglob)] - F = svd(L) - @debug (vander_cond = max(vander_cond, cond(L))) maxlog = 0 - @debug (shift_norm = max(shift_norm, 1)) maxlog = 0 - @debug (vander_norm = max(vander_norm, norm(L))) maxlog = 0 + # Fill the interpolation matrix + for k in 1:nq, m in 1:num_basis + L_arr[m, k] = p[m](view(source, jglob)[k].coords) end + F = svd(Ldata) + @debug (vander_cond = max(vander_cond, cond(Ldata))) maxlog = 0 + @debug (shift_norm = max(shift_norm, 1)) maxlog = 0 + @debug (vander_norm = max(vander_norm, norm(Ldata))) maxlog = 0 # correct each target near the current element for i in near_list[n] - b = @views R[i, :] - wei = SHIFT ? F \ (S * b) : F \ b # weights for the current element and target i - rhs_norm = max(rhs_norm, norm(b)) - res_norm = if SHIFT - max(res_norm, norm(L̃ * wei - S * b)) - else - max(res_norm, norm(L * wei - b)) - end + b_arr .= @views Θ[i:i, :]' + @debug (rhs_norm = max(rhs_norm, norm(bdata))) maxlog = 0 + ldiv!(weidata, F, bdata) + @debug (res_norm = max(res_norm, norm(Ldata * weidata - bdata))) maxlog = 0 for k in 1:nq push!(Is, i) push!(Js, jglob[k]) - push!(Vs, wei[k]) + push!(Vs, wei_arr[k, 1]) end end end @@ -124,139 +122,10 @@ function vdim_correction( |-- max interp. matrix norm : $vander_norm |-- max shift norm : $shift_norm """ - δV = sparse(Is, Js, Vs, m, n) + δV = sparse(Is, Js, Vs, num_target, num_source) return δV end -function change_of_basis(multiindices, p, c, r) - nbasis = length(multiindices) - P = zeros(nbasis, nbasis) - for i in 1:nbasis - α = multiindices[i] - for j in 1:nbasis - β = multiindices[j] - β ≤ α || continue - # P[i, j] = prod((-c) .^ ((α - β).indices)) / r^abs(α) / factorial(α - # - β) - γ = α - β - p_γ = p[findfirst(x -> x == γ, multiindices)] # p_{\alpha - \beta} - P[i, j] = p_γ(-c) / r^abs(α) - end - end - return P -end - -function translation_and_scaling(el::LagrangeTriangle) - vertices = el.vals[1:3] - l1 = norm(vertices[1] - vertices[2]) - l2 = norm(vertices[2] - vertices[3]) - l3 = norm(vertices[3] - vertices[1]) - if ((l1^2 + l2^2 >= l3^2) && (l2^2 + l3^2 >= l1^2) && (l3^2 + l1^2 > l2^2)) - acuteright = true - else - acuteright = false - end - - if acuteright - # Compute the circumcenter and circumradius - Bp = vertices[2] - vertices[1] - 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)) - Up = SVector{2}(Upx, Upy) - r = norm(Up) - c = Up + vertices[1] - else - if (l1 >= l2) && (l1 >= l3) - c = (vertices[1] + vertices[2]) / 2 - r = l1 / 2 - elseif (l2 >= l1) && (l2 >= l3) - c = (vertices[2] + vertices[3]) / 2 - r = l2 / 2 - else - c = (vertices[1] + vertices[3]) / 2 - r = l3 / 2 - end - end - return c, r -end - -function translation_and_scaling(el::ParametricElement{ReferenceSimplex{3}}) - straight_nodes = - [el([0.0, 0.0, 0.0]), el([1.0, 0.0, 0.0]), el([0.0, 1.0, 0.0]), el([0.0, 0.0, 1.0])] - return translation_and_scaling( - LagrangeElement{ReferenceSimplex{3}, 4, SVector{3, Float64}}(straight_nodes), - ) -end - -function translation_and_scaling(el::ParametricElement{ReferenceSimplex{2}}) - straight_nodes = [el([1.0e-18, 1.0e-18]), el([1.0, 0.0]), el([0.0, 1.0])] - return translation_and_scaling( - LagrangeElement{ReferenceSimplex{2}, 3, SVector{2, Float64}}(straight_nodes), - ) -end - -function translation_and_scaling(el::LagrangeTetrahedron) - vertices = el.vals[1:4] - # Compute the circumcenter in barycentric coordinates - # formulas here are due to: https://math.stackexchange.com/questions/2863613/tetrahedron-centers - a = norm(vertices[4] - vertices[1]) - b = norm(vertices[2] - vertices[4]) - c = norm(vertices[3] - vertices[4]) - d = norm(vertices[3] - vertices[2]) - e = norm(vertices[3] - vertices[1]) - f = norm(vertices[2] - vertices[1]) - f² = f^2 - a² = a^2 - b² = b^2 - c² = c^2 - d² = d^2 - e² = e^2 - - ρ = - a² * d² * (-d² + e² + f²) + b² * e² * (d² - e² + f²) + c² * f² * (d² + e² - f²) - - 2 * d² * e² * f² - α = - a² * d² * (b² + c² - d²) + e² * b² * (-b² + c² + d²) + f² * c² * (b² - c² + d²) - - 2 * b² * c² * d² - β = - b² * e² * (a² + c² - e²) + d² * a² * (-a² + c² + e²) + f² * c² * (a² - c² + e²) - - 2 * a² * c² * e² - γ = - c² * f² * (a² + b² - f²) + d² * a² * (-a² + b² + f²) + e² * b² * (a² - b² + f²) - - 2 * a² * b² * f² - if (ρ >= 0 && α >= 0 && β >= 0 + γ >= 0) - # circumcenter lays inside `el` - center = - (α * vertices[1] + β * vertices[2] + γ * vertices[3] + ρ * vertices[4]) / - (ρ + α + β + γ) - # ref: https://math.stackexchange.com/questions/1087011/calculating-the-radius-of-the-circumscribed-sphere-of-an-arbitrary-tetrahedron - R = sqrt(1 / 2 * (β * f² + γ * e² + ρ * a²) / (ρ + α + β + γ)) - else - if (a >= b && a >= c && a >= d && a >= e && a >= f) - center = (vertices[1] + vertices[4]) / 2 - R = a / 2 - elseif (b >= a && b >= c && b >= d && b >= e && b >= f) - center = (vertices[2] + vertices[4]) / 2 - R = b / 2 - elseif (c >= a && c >= b && c >= d && c >= e && c >= f) - center = (vertices[3] + vertices[4]) / 2 - R = c / 2 - elseif (d >= a && d >= b && d >= c && d >= e && d >= f) - center = (vertices[3] + vertices[2]) / 2 - R = d / 2 - elseif (e >= a && e >= b && e >= c && e >= d && e >= f) - center = (vertices[3] + vertices[1]) / 2 - R = e / 2 - else - center = (vertices[2] + vertices[1]) / 2 - R = f / 2 - end - end - return center, R -end - function _vdim_auxiliary_quantities( p, P, @@ -274,7 +143,8 @@ function _vdim_auxiliary_quantities( b = [f(q) for q in Y, f in p] γ₀B = [f(q) for q in Γ, f in P] γ₁B = [f(q) for q in Γ, f in γ₁P] - Θ = zeros(eltype(Vop), num_targets, num_basis) + T = eltype(γ₀B) + Θ = zeros(T, 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], Sop, γ₁B[:, n]) @@ -288,113 +158,127 @@ function _vdim_auxiliary_quantities( end """ - vdim_mesh_center(msh) - -Point `x` which minimizes ∑ (x-xⱼ)²/r²ⱼ, where xⱼ and rⱼ are the circumcenter -and circumradius of the elements of `msh`, respectively. -""" -function vdim_mesh_center(msh::AbstractMesh) - N = ambient_dimension(msh) - M = 0.0 - xc = zero(SVector{N, Float64}) - for E in element_types(msh) - for el in elements(msh, E) - c, r = translation_and_scaling(el) - # w = 1/r^2 - w = 1 - M += w - xc += c * w - end - end - return xc / M -end + polynomial_solutions_vdim(op, order) -""" - polynomial_solutions_vdim(op, order[, center]) - -For every monomial term `pₙ` of degree `order`, compute a polynomial `Pₙ` such -that `ℒ[Pₙ] = pₙ`, where `ℒ` is the differential operator associated with `op`. -This function returns `{pₙ,Pₙ,γ₁Pₙ}`, where `γ₁Pₙ` is the generalized Neumann -trace of `Pₙ`. - -Passing a point `center` will shift the monomials and solutions accordingly. +For every monomial term `pₙ` of degree at most `order`, compute a polynomial `Pₙ` such that +`ℒ[Pₙ] = pₙ`, where `ℒ` is the differential operator associated with `op`. This function +returns `{pₙ,Pₙ,γ₁Pₙ, Iₙ}`, where `γ₁Pₙ` is the generalized Neumann trace of `Pₙ`, and `Iₙ` +is the multi-index associated with `pₙ`. """ function polynomial_solutions_vdim( - op::AbstractDifferentialOperator, + op::AbstractDifferentialOperator{N}, order::Integer, - center = nothing, - ) - N = ambient_dimension(op) - center = isnothing(center) ? zero(SVector{N, Float64}) : center - # 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}}() - dirchlet_traces = Vector{ElementaryPDESolutions.Polynomial{N, Float64}}() - T = return_type(neumann_trace, typeof(op), eltype(dirchlet_traces)) - neumann_traces = Vector{T}() - multiindices = Vector{MultiIndex{N}}() - # iterate over N-tuples going from 0 to order + ::Type{T} = default_kernel_eltype(op), + ) where {N, T} + # create empty arrays to store the monomials, solutions, and traces. + monomials = Vector{Polynomial{N, T}}() + solutions = Vector{Polynomial{N, T}}() 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) - γ₁P = neumann_trace(op, P) - push!(multiindices, MultiIndex(I)) - push!(monomials, p) - push!(dirchlet_traces, P) - push!(neumann_traces, γ₁P) - end - monomial_shift = map(monomials) do f - return (q) -> f(coords(q) - center) + monomial = Polynomial(I => one(T)) + sol = monomial_solution(op, monomial) + push!(monomials, monomial) + push!(solutions, sol) end - dirchlet_shift = map(dirchlet_traces) do f - return (q) -> f(coords(q) - center) + neumann_traces = map(solutions) do sol + neumann_trace(op, sol) end - neumann_shift = map(neumann_traces) do f - return (q) -> f((coords = q.coords - center, normal = q.normal)) - end - return monomial_shift, dirchlet_shift, neumann_shift, multiindices - # return monomials, dirchlet_traces, neumann_traces, multiindices + return monomials, solutions, neumann_traces end -# dispatch to the correct solver in ElementaryPDESolutions -function polynomial_solution(::Laplace, p::ElementaryPDESolutions.Polynomial) - P = ElementaryPDESolutions.solve_laplace(p) - return ElementaryPDESolutions.convert_coefs(P, Float64) -end +""" + monomial_solution(op::AbstractDifferentialOperator{N}, p::Polynomial{N}) + +Compute a polynomial solution `P` to the equation `ℒ[P] = p`, where + +- `ℒ` is the differential operator associated with `op`, +- `p(x) = x₁^α₁ * x₂^α₂ * ... * x_N^α_N * I` is the monomial of multi-index `α`, +- `I` is the `σ × σ` identity tensor, with `σ` the range dimension of `op` (i.e., number of + components of the solution). -function polynomial_solution(op::Helmholtz, p::ElementaryPDESolutions.Polynomial) - k = op.k - P = ElementaryPDESolutions.solve_helmholtz(p; k) - return ElementaryPDESolutions.convert_coefs(P, Float64) +Returns `(p, P, γ₁P)`, where `γ₁P` is the generalized Neumann trace of `P`. + +Both `p` and `P` are given as `ElementaryPDESolutions.Polynomial` objects, while +""" +monomial_solution(::Laplace{N}, p::Polynomial{N}) where {N} = ElementaryPDESolutions.solve_laplace(p) + +function neumann_trace(::Laplace{N}, p::Polynomial{N, T}) where {N, T} + ∇P = ElementaryPDESolutions.gradient(p) + γ₁P = (q) -> dot(normal(q), ∇P(coords(q))) + return γ₁P end -function polynomial_solution(op::Yukawa, p::ElementaryPDESolutions.Polynomial) - k = im * op.λ - P = ElementaryPDESolutions.solve_helmholtz(p; k) - return ElementaryPDESolutions.convert_coefs(P, Float64) +function monomial_solution(op::Elastostatic{N}, p::Polynomial{N, T}) where {N, T} + @assert T <: StaticMatrix + @assert size(T) == (N, N) + # extract exponent of p, and make sure it is a monomial + ord2coef = p.order2coeff + @assert length(ord2coef) == 1 "Input polynomial must be a monomial" + coef = first(values(ord2coef)) + val = first(keys(ord2coef)) + # extract material parameters + μ = op.μ + λ = op.λ + ν = λ / (2 * (λ + μ)) + # compute solution for each columns of the tensor + ptuple = ntuple(N) do n + p = ntuple(d -> Polynomial(val => coef[d, n]), N) + ElementaryPDESolutions.solve_elastostatic(p; μ, ν) + end + p = flatten_polynomial_ntuple(ptuple) + return p end function neumann_trace( - ::Union{Laplace, Helmholtz, Yukawa}, - P::ElementaryPDESolutions.Polynomial{N, T}, + op::Elastostatic{N}, + P::Polynomial{N, T}, ) where {N, T} - return _normal_derivative(P) + μ = op.μ + λ = op.λ + ∇P = ElementaryPDESolutions.gradient(P) + γ₁P = (q) -> begin + ν = normal(q) + x = coords(q) + M = ∇P(x) + cols = svector(N) do i + divu = tr(M[i]) + return λ * divu * ν + μ * (M[i] + M[i]') * ν + end + reduce(hcat, cols) + end + return γ₁P end -function _normal_derivative(P::ElementaryPDESolutions.Polynomial{N, T}) where {N, T} - ∇P = ElementaryPDESolutions.gradient(P) - return (q) -> dot(normal(q), ∇P(coords(q))) + +function flatten_polynomial_ntuple(P::NTuple{N, NTuple{N, <:Polynomial{DIM, T}}}) where {N, DIM, T <: Number} + V = SMatrix{N, N, T, N * N} + # collect all multi-indices + idxs = Set{NTuple{DIM, Int}}() + foreach(p -> union!(idxs, keys(p.order2coeff)), Iterators.flatten(P)) + # now loop over keys and build flattened coefficients + idx2coef = Dict{NTuple{DIM, Int}, V}() + for idx in idxs + coef_tuple = ntuple(N^2) do n + m = div(n - 1, N) + 1 + l = mod(n - 1, N) + 1 + p = P[m][l] + get(p.order2coeff, idx, zero(T)) + end + idx2coef[idx] = V(coef_tuple) + end + return Polynomial{DIM, V}(idx2coef) +end + +function (P::NTuple{N, <:Polynomial})(x) where {N} + return svector(n -> P[n](x), N) end -function (∇P::NTuple{N, <:ElementaryPDESolutions.Polynomial})(x) where {N} - return ntuple(n -> ∇P[n](x), N) +function (P::Polynomial)(q::QuadratureNode) + x = coords(q) + return P(x) end -function (P::ElementaryPDESolutions.Polynomial)(q::QuadratureNode) - x = q.coords.data +function (P::NTuple{N, <:Polynomial})(q::QuadratureNode) where {N} + x = coords(q) return P(x) end From 72ecf96e9bdae14b03650c995c7f6068504ecb89 Mon Sep 17 00:00:00 2001 From: maltezfaria Date: Mon, 19 Jan 2026 13:41:52 +0100 Subject: [PATCH 06/20] fix --- src/vdim.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/vdim.jl b/src/vdim.jl index ef699851..19ca80ea 100644 --- a/src/vdim.jl +++ b/src/vdim.jl @@ -95,7 +95,7 @@ function vdim_correction( jglob = @view qtags[:, n] # Fill the interpolation matrix for k in 1:nq, m in 1:num_basis - L_arr[m, k] = p[m](view(source, jglob)[k].coords) + L_arr[m, k] = p[m](view(source, jglob)[k]) end F = svd(Ldata) @debug (vander_cond = max(vander_cond, cond(Ldata))) maxlog = 0 @@ -239,10 +239,12 @@ function neumann_trace( γ₁P = (q) -> begin ν = normal(q) x = coords(q) - M = ∇P(x) - cols = svector(N) do i - divu = tr(M[i]) - return λ * divu * ν + μ * (M[i] + M[i]') * ν + M = ∇P(x) # M[j] = ∂P/∂xⱼ + cols = svector(N) do m + # Build gradient of m-th column: (∇uₘ)[:, j] = M[j][:, m] + gradu = hcat(ntuple(j -> M[j][:, m], N)...) + divu = tr(gradu) + return λ * divu * ν + μ * (gradu + gradu') * ν end reduce(hcat, cols) end From c81c4a00d90e5504a04fbfc0defef0c870ccde04 Mon Sep 17 00:00:00 2001 From: maltezfaria Date: Mon, 19 Jan 2026 13:42:01 +0100 Subject: [PATCH 07/20] add tests --- test/volume_potential_test.jl | 168 ++++++++++++++++++++++++++++++++++ 1 file changed, 168 insertions(+) create mode 100644 test/volume_potential_test.jl diff --git a/test/volume_potential_test.jl b/test/volume_potential_test.jl new file mode 100644 index 00000000..de7f9ce0 --- /dev/null +++ b/test/volume_potential_test.jl @@ -0,0 +1,168 @@ +#= + +Basic tests for the volume potentials based on + +S*γ₁u - D*γ₀u + V*f - σ*γ₀u ≈ 0 when Lu = f + +where S, D are the single and double layer boundary integral operators, V is the volume +potential operator, γ₀u and γ₁u are the Dirichlet and Neumann traces of u on the boundary, σ +is 1/2 or -1/2 depending on whether the target point is inside or outside the domain, and L +is the differential operator (Laplace, Helmholtz, Stokes, Elastostatic). + +For testing, we use known polynomial solutions of the PDEs to generate u and f = Lu. + +=# + +using Test +using LinearAlgebra +using Inti +using Random +using StaticArrays +using ForwardDiff + +Random.seed!(1) + +## Test parameters +rtol = 1.0e-2 # relative tolerance for volume potential tests +meshsize = 0.2 +meshorder = 1 +bdry_qorder = 8 +interpolation_order = 2 + +""" + test_volume_potential(op, Ω, Γ, meshsize; interpolation_order=2, bdry_qorder=8) + +Test the volume potential operator for a given PDE operator `op` on domain `Ω` with boundary `Γ`. +Verifies the identity: S*γ₁u - D*γ₀u + V*f - σ*γ₀u ≈ 0 for polynomial solutions. +""" +function test_volume_potential(op, Ω, Γ, msh; interpolation_order = 2, bdry_qorder = 8) + # TODO: make it work for 3D as well + Ωₕ = view(msh, Ω) + Γₕ = view(msh, Γ) + + # Setup quadrature + VR_qorder = Inti.Triangle_VR_interpolation_order_to_quadrature_order(interpolation_order) + 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 = bdry_qorder) + + # Build boundary operators + S_b2d, D_b2d = Inti.single_double_layer(; + op, + target = Ωₕ_quad, + source = Γₕ_quad, + compression = (method = :none,), + correction = (method = :dim, maxdist = 5 * meshsize, target_location = :inside), + ) + + # Build volume potential + V_d2d = Inti.volume_potential(; + op, + target = Ωₕ_quad, + source = Ωₕ_quad, + compression = (method = :none,), + correction = (method = :none,), + ) + + # Build VDIM correction + δV_d2d = Inti.vdim_correction( + op, Ωₕ_quad, Ωₕ_quad, Γₕ_quad, S_b2d, D_b2d, V_d2d; + green_multiplier = -ones(length(Ωₕ_quad)), + interpolation_order, + maxdist = Inf + ) + + # Get polynomial solutions for this operator + monomials, dir_traces, neumann_traces = Inti.polynomial_solutions_vdim(op, interpolation_order) + + # Test with each polynomial basis function + errors_uncorrected = Float64[] + errors_corrected = Float64[] + + for idx in 1:length(monomials) + # For vector-valued operators, multiply by a constant vector + if Inti.default_density_eltype(op) <: SVector + N = Inti.ambient_dimension(Ωₕ) + c = SVector(ntuple(i -> rand(), N)...) + f = (q) -> monomials[idx](q) * c + u = (q) -> dir_traces[idx](q) * c + t = (q) -> neumann_traces[idx](q) * c + else + # Scalar case + f = (q) -> monomials[idx](q) + u = (q) -> dir_traces[idx](q) + t = (q) -> neumann_traces[idx](q) + end + + # Evaluate on quadrature nodes + u_d = map(q -> u(q), Ωₕ_quad) + u_b = map(q -> u(q), Γₕ_quad) + du_b = map(q -> t(q), Γₕ_quad) + f_d = map(q -> f(q), Ωₕ_quad) + + # Compute reference solution: -u - D*u_b + S*du_b + # This comes from Green's representation: u = S*t - D*u + V*f + # So V*f = u - S*t + D*u, and we test -V*f = S*t - D*u - u + vref = -u_d - D_b2d * u_b + S_b2d * du_b + + # Compute uncorrected approximation + vapprox = V_d2d * f_d + er = vref - vapprox + push!(errors_uncorrected, norm(er, Inf)) + + # Compute corrected approximation + vapprox_corr = vapprox + δV_d2d * f_d + er_corr = vref - vapprox_corr + push!(errors_corrected, norm(er_corr, Inf)) + end + + return errors_uncorrected, errors_corrected +end + +@testset "Volume potential operators" begin + @testset "2D Laplace" begin + Inti.clear_entities!() + gmsh.initialize() + gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) + gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) + gmsh.model.occ.addDisk(0, 0, 0, 1, 1) + gmsh.model.occ.synchronize() + gmsh.model.mesh.generate(2) + gmsh.model.mesh.setOrder(meshorder) + msh = Inti.import_mesh(; dim = 2) + Ω = Inti.Domain(e -> Inti.geometric_dimension(e) == 2, Inti.entities(msh)) + gmsh.finalize() + + Γ = Inti.external_boundary(Ω) + op = Inti.Laplace(; dim = 2) + + err_uncorr, err_corr = test_volume_potential(op, Ω, Γ, msh; interpolation_order, bdry_qorder) + + @test maximum(err_corr) < rtol + @test maximum(err_corr) < maximum(err_uncorr) # Correction should improve accuracy + end + + @testset "2D Elastostatic" begin + Inti.clear_entities!() + gmsh.initialize() + gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) + gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) + gmsh.model.occ.addDisk(0, 0, 0, 1, 1) + gmsh.model.occ.synchronize() + gmsh.model.mesh.generate(2) + gmsh.model.mesh.setOrder(meshorder) + msh = Inti.import_mesh(; dim = 2) + Ω = Inti.Domain(e -> Inti.geometric_dimension(e) == 2, Inti.entities(msh)) + gmsh.finalize() + + Γ = Inti.external_boundary(Ω) + op = Inti.Elastostatic(; μ = 1.0, λ = 1.0, dim = 2) + + err_uncorr, err_corr = test_volume_potential(op, Ω, Γ, msh; interpolation_order, bdry_qorder) + @show err_uncorr + @show err_corr + @test maximum(err_corr) < rtol + @test maximum(err_corr) < maximum(err_uncorr) + end +end From 4161ea4f7264666b45b376b894834887ac6cda9e Mon Sep 17 00:00:00 2001 From: maltezfaria Date: Mon, 19 Jan 2026 14:12:02 +0100 Subject: [PATCH 08/20] add helmholtz --- src/vdim.jl | 11 ++++++++++- test/volume_potential_test.jl | 25 +++++++++++++++++++++++-- 2 files changed, 33 insertions(+), 3 deletions(-) diff --git a/src/vdim.jl b/src/vdim.jl index 19ca80ea..53b5f13e 100644 --- a/src/vdim.jl +++ b/src/vdim.jl @@ -208,6 +208,16 @@ function neumann_trace(::Laplace{N}, p::Polynomial{N, T}) where {N, T} return γ₁P end +function monomial_solution(op::Helmholtz{N}, p::Polynomial{N}) where {N} + return ElementaryPDESolutions.solve_helmholtz(p, op.k^2) +end + +function neumann_trace(::Helmholtz{N}, p::Polynomial{N, T}) where {N, T} + ∇P = ElementaryPDESolutions.gradient(p) + γ₁P = (q) -> dot(normal(q), ∇P(coords(q))) + return γ₁P +end + function monomial_solution(op::Elastostatic{N}, p::Polynomial{N, T}) where {N, T} @assert T <: StaticMatrix @assert size(T) == (N, N) @@ -251,7 +261,6 @@ function neumann_trace( return γ₁P end - function flatten_polynomial_ntuple(P::NTuple{N, NTuple{N, <:Polynomial{DIM, T}}}) where {N, DIM, T <: Number} V = SMatrix{N, N, T, N * N} # collect all multi-indices diff --git a/test/volume_potential_test.jl b/test/volume_potential_test.jl index de7f9ce0..2f9ca4e0 100644 --- a/test/volume_potential_test.jl +++ b/test/volume_potential_test.jl @@ -19,6 +19,7 @@ using Inti using Random using StaticArrays using ForwardDiff +using Gmsh Random.seed!(1) @@ -143,6 +144,28 @@ end @test maximum(err_corr) < maximum(err_uncorr) # Correction should improve accuracy end + @testset "2D Helmholtz" begin + Inti.clear_entities!() + gmsh.initialize() + gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) + gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) + gmsh.model.occ.addDisk(0, 0, 0, 1, 1) + gmsh.model.occ.synchronize() + gmsh.model.mesh.generate(2) + gmsh.model.mesh.setOrder(meshorder) + msh = Inti.import_mesh(; dim = 2) + Ω = Inti.Domain(e -> Inti.geometric_dimension(e) == 2, Inti.entities(msh)) + gmsh.finalize() + + Γ = Inti.external_boundary(Ω) + op = Inti.Helmholtz(; k = 1.0, dim = 2) + + err_uncorr, err_corr = test_volume_potential(op, Ω, Γ, msh; interpolation_order, bdry_qorder) + + @test maximum(err_corr) < rtol + @test maximum(err_corr) < maximum(err_uncorr) + end + @testset "2D Elastostatic" begin Inti.clear_entities!() gmsh.initialize() @@ -160,8 +183,6 @@ end op = Inti.Elastostatic(; μ = 1.0, λ = 1.0, dim = 2) err_uncorr, err_corr = test_volume_potential(op, Ω, Γ, msh; interpolation_order, bdry_qorder) - @show err_uncorr - @show err_corr @test maximum(err_corr) < rtol @test maximum(err_corr) < maximum(err_uncorr) end From 2cdb9d95fad0584b94cc0cbdccecbbeaf5c841df Mon Sep 17 00:00:00 2001 From: maltezfaria Date: Mon, 19 Jan 2026 14:14:27 +0100 Subject: [PATCH 09/20] refactor test --- test/vdim_elastostatic_test.jl | 87 -------------------------------- test/volume_potential_test.jl | 90 ++++++++++------------------------ 2 files changed, 27 insertions(+), 150 deletions(-) delete mode 100644 test/vdim_elastostatic_test.jl diff --git a/test/vdim_elastostatic_test.jl b/test/vdim_elastostatic_test.jl deleted file mode 100644 index 983ae082..00000000 --- a/test/vdim_elastostatic_test.jl +++ /dev/null @@ -1,87 +0,0 @@ -# # High-order convergence of vdim - -using Inti -using StaticArrays -using Gmsh -using LinearAlgebra -using HMatrices -using CairoMakie - -meshsize = 0.2 -meshorder = 1 -bdry_qorder = 8 -interpolation_order = 2 - -Inti.clear_entities!() -gmsh.initialize() -gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) -gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) -gmsh.model.occ.addDisk(0, 0, 0, 1, 1) -gmsh.model.occ.synchronize() -gmsh.model.mesh.generate(2) -gmsh.model.mesh.setOrder(meshorder) -msh = Inti.import_mesh(; dim = 2) -Ω = Inti.Domain(e -> Inti.geometric_dimension(e) == 2, Inti.entities(msh)) -gmsh.finalize() - -Γ = Inti.external_boundary(Ω) -Ωₕ = view(msh, Ω) -Γₕ = view(msh, Γ) -## - -VR_qorder = Inti.Triangle_VR_interpolation_order_to_quadrature_order(interpolation_order) - -# 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) -# Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorders[1]) -Γₕ_quad = Inti.Quadrature(Γₕ; qorder = bdry_qorder) - -## - -# build exact solution -μ = 1.0 -λ = 1.0 - -u = (x) -> SVector(1.0, 1.0) -du = (x, n) -> SVector(0.0, 0.0) -f = (x) -> SVector(0.0, 0.0) - -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 = Inti.Elastostatic(; λ = λ, μ = μ, dim = 2) - -# m, d, n = Inti._polynomial_solutions_vec(op, 1) -# x = @SVector rand(2) -# m[1](x) -# d[1](x) -# n[1]((; coords = x, normal = x)) - -## Boundary operators -S_b2d, D_b2d = Inti.single_double_layer(; - op, - target = Ωₕ_quad, - source = Γₕ_quad, - compression = (method = :none, tol = 1.0e-14), - correction = (method = :dim, maxdist = 5 * meshsize, target_location = :inside), -) - -## Volume potentials -V_d2d = Inti.volume_potential(; - op, - target = Ωₕ_quad, - source = Ωₕ_quad, - compression = (method = :none, tol = 1.0e-14), - correction = (method = :dim, interpolation_order), -) - -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), t diff --git a/test/volume_potential_test.jl b/test/volume_potential_test.jl index 2f9ca4e0..34d434f7 100644 --- a/test/volume_potential_test.jl +++ b/test/volume_potential_test.jl @@ -122,68 +122,32 @@ function test_volume_potential(op, Ω, Γ, msh; interpolation_order = 2, bdry_qo end @testset "Volume potential operators" begin - @testset "2D Laplace" begin - Inti.clear_entities!() - gmsh.initialize() - gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) - gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) - gmsh.model.occ.addDisk(0, 0, 0, 1, 1) - gmsh.model.occ.synchronize() - gmsh.model.mesh.generate(2) - gmsh.model.mesh.setOrder(meshorder) - msh = Inti.import_mesh(; dim = 2) - Ω = Inti.Domain(e -> Inti.geometric_dimension(e) == 2, Inti.entities(msh)) - gmsh.finalize() - - Γ = Inti.external_boundary(Ω) - op = Inti.Laplace(; dim = 2) - - err_uncorr, err_corr = test_volume_potential(op, Ω, Γ, msh; interpolation_order, bdry_qorder) - - @test maximum(err_corr) < rtol - @test maximum(err_corr) < maximum(err_uncorr) # Correction should improve accuracy - end - - @testset "2D Helmholtz" begin - Inti.clear_entities!() - gmsh.initialize() - gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) - gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) - gmsh.model.occ.addDisk(0, 0, 0, 1, 1) - gmsh.model.occ.synchronize() - gmsh.model.mesh.generate(2) - gmsh.model.mesh.setOrder(meshorder) - msh = Inti.import_mesh(; dim = 2) - Ω = Inti.Domain(e -> Inti.geometric_dimension(e) == 2, Inti.entities(msh)) - gmsh.finalize() - - Γ = Inti.external_boundary(Ω) - op = Inti.Helmholtz(; k = 1.0, dim = 2) - - err_uncorr, err_corr = test_volume_potential(op, Ω, Γ, msh; interpolation_order, bdry_qorder) - - @test maximum(err_corr) < rtol - @test maximum(err_corr) < maximum(err_uncorr) - end - - @testset "2D Elastostatic" begin - Inti.clear_entities!() - gmsh.initialize() - gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) - gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) - gmsh.model.occ.addDisk(0, 0, 0, 1, 1) - gmsh.model.occ.synchronize() - gmsh.model.mesh.generate(2) - gmsh.model.mesh.setOrder(meshorder) - msh = Inti.import_mesh(; dim = 2) - Ω = Inti.Domain(e -> Inti.geometric_dimension(e) == 2, Inti.entities(msh)) - gmsh.finalize() - - Γ = Inti.external_boundary(Ω) - op = Inti.Elastostatic(; μ = 1.0, λ = 1.0, dim = 2) - - err_uncorr, err_corr = test_volume_potential(op, Ω, Γ, msh; interpolation_order, bdry_qorder) - @test maximum(err_corr) < rtol - @test maximum(err_corr) < maximum(err_uncorr) + # Create 2D geometry once + Inti.clear_entities!() + gmsh.initialize() + gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) + gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) + gmsh.model.occ.addDisk(0, 0, 0, 1, 1) + gmsh.model.occ.synchronize() + gmsh.model.mesh.generate(2) + gmsh.model.mesh.setOrder(meshorder) + msh = Inti.import_mesh(; dim = 2) + Ω = Inti.Domain(e -> Inti.geometric_dimension(e) == 2, Inti.entities(msh)) + gmsh.finalize() + Γ = Inti.external_boundary(Ω) + + # Define operators to test + operators = [ + ("2D Laplace", Inti.Laplace(; dim = 2)), + ("2D Helmholtz", Inti.Helmholtz(; k = 1.0, dim = 2)), + ("2D Elastostatic", Inti.Elastostatic(; μ = 1.0, λ = 1.0, dim = 2)), + ] + + for (name, op) in operators + @testset "$name" begin + err_uncorr, err_corr = test_volume_potential(op, Ω, Γ, msh; interpolation_order, bdry_qorder) + @test maximum(err_corr) < rtol + @test maximum(err_corr) < maximum(err_uncorr) # Correction should improve accuracy + end end end From e1e09557996ce10b2b9bff617afc6aaba2ee19e6 Mon Sep 17 00:00:00 2001 From: maltezfaria Date: Mon, 19 Jan 2026 14:22:50 +0100 Subject: [PATCH 10/20] tests in 2D and 3D --- src/Inti.jl | 17 +-- test/curved_test_2d.jl | 224 ++++++++++++++++++++++++++++++++++ test/curved_test_3d.jl | 77 ++++++++++++ test/volume_potential_test.jl | 58 +++++++-- 4 files changed, 360 insertions(+), 16 deletions(-) create mode 100644 test/curved_test_2d.jl create mode 100644 test/curved_test_3d.jl diff --git a/src/Inti.jl b/src/Inti.jl index 7df4521f..0da21d84 100644 --- a/src/Inti.jl +++ b/src/Inti.jl @@ -7,24 +7,26 @@ module Inti const PROJECT_ROOT = pkgdir(Inti) +import Bessels +import ElementaryPDESolutions +import HAdaptiveIntegration +import SpecialFunctions + +import ElementaryPDESolutions: Polynomial + using DataStructures using ForwardDiff using LinearAlgebra using LinearMaps using NearestNeighbors using Pkg +using Printf using QuadGK +using Richardson using Scratch using SparseArrays using StaticArrays -using Printf using TOML -using Richardson - -import ElementaryPDESolutions -import SpecialFunctions -import Bessels # faster than SpecialFunctions for Bessel functions with real args -import HAdaptiveIntegration # helper functions include("utils.jl") @@ -32,6 +34,7 @@ include("blockarray.jl") # basic interpolation and integration include("reference_shapes.jl") + include("polynomials.jl") include("reference_interpolation.jl") include("quad_rules_tables.jl") diff --git a/test/curved_test_2d.jl b/test/curved_test_2d.jl new file mode 100644 index 00000000..ff42d19e --- /dev/null +++ b/test/curved_test_2d.jl @@ -0,0 +1,224 @@ +using Inti +using Meshes +using StaticArrays +using GLMakie +using Gmsh +using LinearAlgebra +using ForwardDiff +using Test + +@testset begin + function domain_and_mesh(; meshsize, meshorder = 1) + Inti.clear_entities!() + gmsh.initialize() + gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) + gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) + gmsh.model.occ.addDisk(0, 0, 0, 1, 1) + gmsh.model.occ.synchronize() + gmsh.model.mesh.generate(2) + gmsh.model.mesh.setOrder(meshorder) + msh = Inti.import_mesh(; dim = 2) + Ω = Inti.Domain(Inti.entities(msh)) do ent + return Inti.geometric_dimension(ent) == 2 + end + gmsh.finalize() + return Ω, msh + end + + meshsize = 0.1 + + tmesh = @elapsed begin + Ω, msh = domain_and_mesh(; meshsize) + end + @info "Mesh generation time: $tmesh" + + Γ = Inti.external_boundary(Ω) + Ωₕ = view(msh, Ω) + Γₕ = view(msh, Γ) + + ψ = (t) -> [cos(2 * π * t), sin(2 * π * t)] + θ = 6 # smoothness order of curved elements + crvmsh = Inti.curve_mesh(msh, ψ, θ) + + Γₕ = crvmsh[Γ] + Ωₕ = crvmsh[Ω] + + qorder = 2 + Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) + Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) + @test isapprox(Inti.integrate(x -> 1, Ωₕ_quad), π, rtol = 1.0e-7) + @test isapprox(Inti.integrate(q -> q.coords[1]^4, Ωₕ_quad), π / 8, rtol = 1.0e-5) + + qorder = 5 + Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) + Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) + @test isapprox(Inti.integrate(x -> 1, Ωₕ_quad), π, rtol = 1.0e-11) + @test isapprox(Inti.integrate(q -> q.coords[1]^4, Ωₕ_quad), π / 8, rtol = 1.0e-10) + + qorder = 8 + Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) + Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) + @test isapprox(Inti.integrate(x -> 1, Ωₕ_quad), π, rtol = 1.0e-14) + @test isapprox(Inti.integrate(q -> q.coords[1]^4, Ωₕ_quad), π / 8, rtol = 1.0e-14) + + Fvol = (x) -> x[2]^2 - 2 * x[2] * x[1]^3 + F = (x) -> [x[1] * x[2]^2, x[1]^3 * x[2]^2] + #Fvol = (x) -> 1.0 + #F = (x) -> [1/2*x[1], 1/2*x[2]] + divvol = Inti.integrate(q -> Fvol(q.coords), Ωₕ_quad) + divline = Inti.integrate(q -> dot(F(q.coords), q.normal), Γₕ_quad) + @test isapprox(divline, divvol, rtol = 1.0e-13) +end + +@testset begin + meshsize = 2π / 4 / 8 + Inti.clear_entities!() + + gmsh.initialize() + gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) + gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) + + # Two kites + kite = Inti.gmsh_curve(0, 1; meshsize) do s + return SVector(0.25, 0.0) + + SVector(cos(2π * s) + 0.65 * cos(4π * s[1]) - 0.65, 1.5 * sin(2π * s)) + end + cl = gmsh.model.occ.addCurveLoop([kite]) + surf = gmsh.model.occ.addPlaneSurface([cl]) + kite_trans = Inti.gmsh_curve(0, 1; meshsize) do s + return SVector(4.5, 0.0) + + SVector(cos(2π * s) + 0.65 * cos(4π * s[1]) - 0.65, 1.5 * sin(2π * s)) + end + cl_trans = gmsh.model.occ.addCurveLoop([kite_trans]) + surf_trans = gmsh.model.occ.addPlaneSurface([cl_trans]) + gmsh.model.occ.synchronize() + + # Add tags for stable identification of the entities + gmsh.model.addPhysicalGroup(2, [surf], -1, "c1") + gmsh.model.addPhysicalGroup(2, [surf_trans], -1, "c2") + + gmsh.model.mesh.generate(2) + msh = Inti.import_mesh(; dim = 2) + Ω = Inti.Domain(Inti.entities(msh)) do ent + return Inti.geometric_dimension(ent) == 2 + end + gmsh.finalize() + + Γ = Inti.external_boundary(Ω) + Ωₕ = view(msh, Ω) + Γₕ = view(msh, Γ) + + # Two kites + ψ₁ = (t) -> [0.25 + cos(2π * t) + 0.65 * cos(4π * t) - 0.65, 1.5 * sin(2π * t)] + ψ₂ = (t) -> [4.5 + cos(2π * t) + 0.65 * cos(4π * t) - 0.65, 1.5 * sin(2π * t)] + entity_parametrizations = Dict{Inti.EntityKey, Function}() + for e in Inti.entities(Ω) + l = Inti.labels(e) + if "c1" in l + entity_parametrizations[e] = ψ₁ + elseif "c2" in l + entity_parametrizations[e] = ψ₂ + elseif "c3" in l + entity_parametrizations[e] = ψ₃ + end + end + θ = 6 # smoothness order of curved elements + crvmsh = Inti.curve_mesh(msh, entity_parametrizations, θ) + + Γₕ = crvmsh[Γ] + Ωₕ = crvmsh[Ω] + + qorder = 8 + Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) + Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) + + Fvol = (x) -> x[2]^2 - 2 * x[2] * x[1]^3 + F = (x) -> [x[1] * x[2]^2, x[1]^3 * x[2]^2] + #Fvol = (x) -> 1.0 + #F = (x) -> [1/2*x[1], 1/2*x[2]] + divvol = Inti.integrate(q -> Fvol(q.coords), Ωₕ_quad) + divline = Inti.integrate(q -> dot(F(q.coords), q.normal), Γₕ_quad) + @test isapprox(divline, divvol, rtol = 1.0e-11) +end + +@testset begin + meshsize = 0.075 + Inti.clear_entities!() + + gmsh.initialize() + gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) + gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) + + # Three circles + c1 = gmsh.model.occ.addDisk(0, 0, 0, 1, 1) + c2 = gmsh.model.occ.addDisk(0, 3.0, 0, 1, 1) + c3 = gmsh.model.occ.addDisk(0, 8.0, 0, 2, 2) + gmsh.model.occ.synchronize() + + # Add tags for stable identification of the entities + gmsh.model.addPhysicalGroup(2, [c1], -1, "c1") + gmsh.model.addPhysicalGroup(2, [c2], -1, "c2") + gmsh.model.addPhysicalGroup(2, [c3], -1, "c3") + + gmsh.model.mesh.generate(2) + msh = Inti.import_mesh(; dim = 2) + Ω = Inti.Domain(Inti.entities(msh)) do ent + return Inti.geometric_dimension(ent) == 2 + end + gmsh.finalize() + + Γ = Inti.external_boundary(Ω) + Ωₕ = view(msh, Ω) + Γₕ = view(msh, Γ) + + # Three circles + ψ₁ = (t) -> [cos(2 * π * t), sin(2 * π * t)] + ψ₂ = (t) -> [cos(2 * π * t), 3.0 + sin(2 * π * t)] + ψ₃ = (t) -> [2 * cos(2 * π * t), 8.0 + 2 * sin(2 * π * t)] + entity_parametrizations = Dict{Inti.EntityKey, Function}() + for e in Inti.entities(Ω) + l = Inti.labels(e) + if "c1" in l + entity_parametrizations[e] = ψ₁ + elseif "c2" in l + entity_parametrizations[e] = ψ₂ + elseif "c3" in l + entity_parametrizations[e] = ψ₃ + end + end + + θ = 6 # smoothness order of curved elements + crvmsh = Inti.curve_mesh(msh, entity_parametrizations, θ) + + Γₕ = crvmsh[Γ] + Ωₕ = crvmsh[Ω] + + qorder = 2 + Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) + Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) + @test isapprox(Inti.integrate(x -> 1, Ωₕ_quad), 6π, rtol = 1.0e-6) + + qorder = 5 + Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) + Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) + @test isapprox(Inti.integrate(x -> 1, Ωₕ_quad), 6π, rtol = 1.0e-11) + + qorder = 8 + Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) + Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) + @test isapprox(Inti.integrate(x -> 1, Ωₕ_quad), 6π, rtol = 1.0e-14) + + # The third circle has larger volume than others, so check that we can pick sub-domains correctly + Ω_sub = Inti.Domain(collect(keys(Ω))[3]) + Ωₕ_sub = crvmsh[Ω_sub] + Ωₕ_sub_quad = Inti.Quadrature(Ωₕ_sub; qorder = qorder) + @test isapprox(Inti.integrate(x -> 1, Ωₕ_sub_quad), 4π, rtol = 1.0e-14) + + Fvol = (x) -> x[2]^2 - 2 * x[2] * x[1]^3 + F = (x) -> [x[1] * x[2]^2, x[1]^3 * x[2]^2] + #Fvol = (x) -> 1.0 + #F = (x) -> [1/2*x[1], 1/2*x[2]] + divvol = Inti.integrate(q -> Fvol(q.coords), Ωₕ_quad) + divline = Inti.integrate(q -> dot(F(q.coords), q.normal), Γₕ_quad) + @test isapprox(divline, divvol, rtol = 1.0e-13) +end diff --git a/test/curved_test_3d.jl b/test/curved_test_3d.jl new file mode 100644 index 00000000..e4e89f7c --- /dev/null +++ b/test/curved_test_3d.jl @@ -0,0 +1,77 @@ +using Test +using Inti +using FMM3D +using Gmsh +using LinearAlgebra +using Random +using StaticArrays +using ForwardDiff + +include("test_utils.jl") + +@testset begin + # create a boundary and area meshes and quadrature only once + Inti.clear_entities!() + meshsize = 0.1 + r1 = 1.0 + r2 = 0.5 + Ω, msh = gmsh_torus(; center = [0.0, 0.0, 0.0], r1 = r1, r2 = r2, meshsize = meshsize) + Γ = Inti.external_boundary(Ω) + Γ_msh = view(msh, Γ) + + function face_element_on_torus(nodelist, R, r) + return all( + [ + (sqrt(node[1]^2 + node[2]^2) - R^2)^2 + node[3]^2 ≈ r^2 for node in nodelist + ] + ) + end + face_element_on_curved_surface = (nodelist) -> face_element_on_torus(nodelist, r1, r2) + + function ψ(v::AbstractVector) + return [ + (r1 + r2 * sin(v[1])) * cos(v[2]), + (r1 + r2 * sin(v[1])) * sin(v[2]), + r2 * cos(v[1]), + ] + end + θ = 6 # smoothness order of curved elements + crvmsh = Inti.curve_mesh( + msh, + ψ, + θ; + face_element_on_curved_surface = face_element_on_curved_surface, + ) + + Γₕ = crvmsh[Γ] + Ωₕ = crvmsh[Ω] + + truevol = 2 * π^2 * r2^2 * r1 + truesfcarea = 4 * π^2 * r1 * r2 + + qorder = 2 + Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) + Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) + @test isapprox(Inti.integrate(x -> 1, Ωₕ_quad), truevol, rtol = 1.0e-5) + @test isapprox(Inti.integrate(x -> 1, Γₕ_quad), truesfcarea, rtol = 1.0e-5) + + qorder = 5 + Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) + Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) + @test isapprox(Inti.integrate(x -> 1, Ωₕ_quad), truevol, rtol = 1.0e-7) + @test isapprox(Inti.integrate(x -> 1, Γₕ_quad), truesfcarea, rtol = 1.0e-7) + + qorder = 7 + Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorder) + Γₕ_quad = Inti.Quadrature(Γₕ; qorder = qorder) + @test isapprox(Inti.integrate(x -> 1, Ωₕ_quad), truevol, rtol = 1.0e-11) + @test isapprox(Inti.integrate(x -> 1, Γₕ_quad), truesfcarea, rtol = 1.0e-14) + + divF = (x) -> x[3] + x[3]^2 + x[2]^3 + F = (x) -> [x[1] * x[3], x[2] * x[3]^2, x[2]^3 * x[3]] + #divF = (x) -> 1.0 + #F = (x) -> 1/3*[x[1], x[2], x[3]] + divtest_vol = Inti.integrate(q -> divF(q.coords), Ωₕ_quad) + divtest_sfc = Inti.integrate(q -> dot(F(q.coords), q.normal), Γₕ_quad) + @test isapprox(divtest_vol, divtest_sfc, rtol = 1.0e-9) +end diff --git a/test/volume_potential_test.jl b/test/volume_potential_test.jl index 34d434f7..6b81f1f4 100644 --- a/test/volume_potential_test.jl +++ b/test/volume_potential_test.jl @@ -25,9 +25,9 @@ Random.seed!(1) ## Test parameters rtol = 1.0e-2 # relative tolerance for volume potential tests -meshsize = 0.2 +meshsize = 0.4 meshorder = 1 -bdry_qorder = 8 +bdry_qorder = 5 interpolation_order = 2 """ @@ -37,13 +37,20 @@ Test the volume potential operator for a given PDE operator `op` on domain `Ω` Verifies the identity: S*γ₁u - D*γ₀u + V*f - σ*γ₀u ≈ 0 for polynomial solutions. """ function test_volume_potential(op, Ω, Γ, msh; interpolation_order = 2, bdry_qorder = 8) - # TODO: make it work for 3D as well Ωₕ = view(msh, Ω) Γₕ = view(msh, Γ) - # Setup quadrature - VR_qorder = Inti.Triangle_VR_interpolation_order_to_quadrature_order(interpolation_order) - Q = Inti.VioreanuRokhlin(; domain = :triangle, order = VR_qorder) + # Setup quadrature based on dimension + dim = Inti.geometric_dimension(Ω) + if dim == 2 + VR_qorder = Inti.Triangle_VR_interpolation_order_to_quadrature_order(interpolation_order) + Q = Inti.VioreanuRokhlin(; domain = :triangle, order = VR_qorder) + elseif dim == 3 + VR_qorder = Inti.Tetrahedron_VR_interpolation_order_to_quadrature_order(interpolation_order) + Q = Inti.VioreanuRokhlin(; domain = :tetrahedron, order = VR_qorder) + else + error("Unsupported dimension: $dim") + end dict = Dict(E => Q for E in Inti.element_types(Ωₕ)) Ωₕ_quad = Inti.Quadrature(Ωₕ, dict) Γₕ_quad = Inti.Quadrature(Γₕ; qorder = bdry_qorder) @@ -121,10 +128,11 @@ function test_volume_potential(op, Ω, Γ, msh; interpolation_order = 2, bdry_qo return errors_uncorrected, errors_corrected end -@testset "Volume potential operators" begin +@testset "Volume potential operators 2D" begin # Create 2D geometry once Inti.clear_entities!() gmsh.initialize() + gmsh.option.setNumber("General.Verbosity", 0) gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) gmsh.model.occ.addDisk(0, 0, 0, 1, 1) @@ -139,8 +147,40 @@ end # Define operators to test operators = [ ("2D Laplace", Inti.Laplace(; dim = 2)), - ("2D Helmholtz", Inti.Helmholtz(; k = 1.0, dim = 2)), - ("2D Elastostatic", Inti.Elastostatic(; μ = 1.0, λ = 1.0, dim = 2)), + ("2D Helmholtz", Inti.Helmholtz(; k = 0.7, dim = 2)), + ("2D Elastostatic", Inti.Elastostatic(; μ = 0.8, λ = 1.3, dim = 2)), + ] + + for (name, op) in operators + @testset "$name" begin + err_uncorr, err_corr = test_volume_potential(op, Ω, Γ, msh; interpolation_order, bdry_qorder) + @test maximum(err_corr) < rtol + @test maximum(err_corr) < maximum(err_uncorr) # Correction should improve accuracy + end + end +end + +@testset "Volume potential operators 3D" begin + # Create 3D geometry once + Inti.clear_entities!() + gmsh.initialize() + gmsh.option.setNumber("General.Verbosity", 0) + 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(3) + gmsh.model.mesh.setOrder(meshorder) + msh = Inti.import_mesh(; dim = 3) + Ω = Inti.Domain(e -> Inti.geometric_dimension(e) == 3, Inti.entities(msh)) + gmsh.finalize() + Γ = Inti.external_boundary(Ω) + + # Define operators to test + operators = [ + ("3D Laplace", Inti.Laplace(; dim = 3)), + ("3D Helmholtz", Inti.Helmholtz(; k = 1.2, dim = 3)), + ("3D Elastostatic", Inti.Elastostatic(; μ = 1.1, λ = 0.9, dim = 3)), ] for (name, op) in operators From 55627708c89ba1645face98447f04072424cd55a Mon Sep 17 00:00:00 2001 From: maltezfaria Date: Mon, 19 Jan 2026 15:12:32 +0100 Subject: [PATCH 11/20] tweak test parameters to reduce testing time --- src/vdim.jl | 2 +- test/aqua_test.jl | 2 +- test/green_identities_test.jl | 8 ++++---- test/kernels_test.jl | 6 +++--- test/runtests.jl | 1 + 5 files changed, 10 insertions(+), 9 deletions(-) diff --git a/src/vdim.jl b/src/vdim.jl index 53b5f13e..30393f53 100644 --- a/src/vdim.jl +++ b/src/vdim.jl @@ -261,7 +261,7 @@ function neumann_trace( return γ₁P end -function flatten_polynomial_ntuple(P::NTuple{N, NTuple{N, <:Polynomial{DIM, T}}}) where {N, DIM, T <: Number} +function flatten_polynomial_ntuple(P::NTuple{N, NTuple{N, Polynomial{DIM, T}}}) where {N, DIM, T <: Number} V = SMatrix{N, N, T, N * N} # collect all multi-indices idxs = Set{NTuple{DIM, Int}}() diff --git a/test/aqua_test.jl b/test/aqua_test.jl index 650e6f22..fad7a67c 100644 --- a/test/aqua_test.jl +++ b/test/aqua_test.jl @@ -6,7 +6,7 @@ using Aqua Aqua.test_all( Inti; ambiguities = false, # test only `Inti` for ambiguities later - unbound_args = true, + unbound_args = (; broken = true), # broken due to use of NTuple in some signatures undefined_exports = true, project_extras = true, stale_deps = true, diff --git a/test/green_identities_test.jl b/test/green_identities_test.jl index b72611e0..3b16001c 100644 --- a/test/green_identities_test.jl +++ b/test/green_identities_test.jl @@ -29,10 +29,10 @@ Random.seed!(1) rtol1 = 1.0e-2 # single and double layer rtol2 = 5.0e-2 # hypersingular (higher tolerance to avoid use of fine mesh + long unit tests) dims = (2, 3) -meshsize = 0.5 +meshsize = 0.2 types = (:interior, :exterior) -corrections = [(method = :dim,), (method = :adaptive, maxdist = 2 * meshsize, rtol = 1.0e-2)] +corrections = [(method = :dim,), (method = :adaptive, maxdist = 2 * meshsize, rtol = 1.0e-2, atol = 1.0e-2)] for correction in corrections @testset "Method = $(correction.method)" begin @@ -43,12 +43,12 @@ for correction in corrections Γ = Inti.parametric_curve(x -> SVector(cos(x), sin(x)), 0.0, 2π) |> Inti.Domain - quad = Inti.Quadrature(Γ; meshsize, qorder = 5) + quad = Inti.Quadrature(Γ; meshsize, qorder = 3) else # Ω, msh = gmsh_ball(; center = [0.0, 0.0, 0.0], radius = 1.0, meshsize = 0.2) Ω = Inti.GeometricEntity("ellipsoid") |> Inti.Domain Γ = Inti.external_boundary(Ω) - quad = Inti.Quadrature(Γ; meshsize, qorder = 5) + quad = Inti.Quadrature(Γ; meshsize, qorder = 3) end for t in types σ = t == :interior ? 1 / 2 : -1 / 2 diff --git a/test/kernels_test.jl b/test/kernels_test.jl index 80ddecb1..51123549 100644 --- a/test/kernels_test.jl +++ b/test/kernels_test.jl @@ -85,19 +85,19 @@ end @test isapprox( ForwardDiff.derivative(t -> G((coords = x + t * nx,), y), 0), dGdnx((coords = x, normal = nx), y); - atol = 1.0e-6, + atol = 1.0e-4, ) @test isapprox( ForwardDiff.derivative(t -> G(x, (coords = y + t * ny,)), 0), dGdny(x, (coords = y, normal = ny)); - atol = 1.0e-6, + atol = 1.0e-4, ) @test isapprox( ForwardDiff.derivative( t -> dGdny((coords = x + t * nx,), (coords = y, normal = ny)), 0, ), d2Gdnxy((coords = x, normal = nx), (coords = y, normal = ny)); - atol = 1.0e-6, + atol = 1.0e-4, ) # test quasi-periodicity period = 2π diff --git a/test/runtests.jl b/test/runtests.jl index 61b6439e..84a705d6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -31,6 +31,7 @@ using Aqua @testset verbose = true "Corrections (Green identities)" begin include("green_identities_test.jl") + include("volume_potential_test.jl") end @testset "Accelerated density interpolation" include("dim_test.jl") From 627afd13b32cb9f246c42aca6abe74e66d2fdf60 Mon Sep 17 00:00:00 2001 From: maltezfaria Date: Mon, 19 Jan 2026 18:44:25 +0100 Subject: [PATCH 12/20] update compat for `Meshes` --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 48fa461f..5c819320 100644 --- a/Project.toml +++ b/Project.toml @@ -54,7 +54,7 @@ HMatrices = "0.2" LinearAlgebra = "1" LinearMaps = "3" Makie = "0.21 - 0.24" -Meshes = "0.42 - 0.55" +Meshes = "0.42 - 0.56" NearestNeighbors = "0.4" Pkg = "1" Printf = "1" From 13c2ed6965c8a4e4d717f2ec66da9f00d5bd416e Mon Sep 17 00:00:00 2001 From: maltezfaria Date: Mon, 19 Jan 2026 18:48:34 +0100 Subject: [PATCH 13/20] try to fix doc generation isssues --- docs/make.jl | 6 +++--- src/vdim.jl | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 2247f492..782bf04c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -143,9 +143,9 @@ makedocs(; mathengine = Documenter.KaTeX( Dict( :delimiters => [ - Dict(:left => raw"$", :right => raw"$", display => false), - Dict(:left => raw"$$", :right => raw"$$", display => true), - Dict(:left => raw"\[", :right => raw"\]", display => true), + Dict(:left => raw"$", :right => raw"$", :display => false), + Dict(:left => raw"$$", :right => raw"$$", :display => true), + Dict(:left => raw"\[", :right => raw"\]", :display => true), ], :macros => Dict( "\\RR" => "\\mathbb{R}", diff --git a/src/vdim.jl b/src/vdim.jl index 30393f53..698505fc 100644 --- a/src/vdim.jl +++ b/src/vdim.jl @@ -162,7 +162,7 @@ end For every monomial term `pₙ` of degree at most `order`, compute a polynomial `Pₙ` such that `ℒ[Pₙ] = pₙ`, where `ℒ` is the differential operator associated with `op`. This function -returns `{pₙ,Pₙ,γ₁Pₙ, Iₙ}`, where `γ₁Pₙ` is the generalized Neumann trace of `Pₙ`, and `Iₙ` +returns `{pₙ,Pₙ,γ₁Pₙ}`, where `γ₁Pₙ` is the generalized Neumann trace of `Pₙ`, and `Iₙ` is the multi-index associated with `pₙ`. """ function polynomial_solutions_vdim( From 66ed2f12f884f691e62d811cb7616de9fe1dd1b2 Mon Sep 17 00:00:00 2001 From: maltezfaria Date: Mon, 19 Jan 2026 21:32:19 +0100 Subject: [PATCH 14/20] restructure polynomial solution and add Stokes --- src/vdim.jl | 184 +++++++++++++++++++--------------- test/volume_potential_test.jl | 18 ++-- 2 files changed, 111 insertions(+), 91 deletions(-) diff --git a/src/vdim.jl b/src/vdim.jl index 698505fc..72a19300 100644 --- a/src/vdim.jl +++ b/src/vdim.jl @@ -51,16 +51,16 @@ function vdim_correction( isnothing(interpolation_order) && (interpolation_order = maximum(order, values(source.etype2qrule))) # by default basis centered at origin - p, P, γ₁P = polynomial_solutions_vdim(op, interpolation_order, T) + basis = polynomial_solutions_vdim(op, interpolation_order, T) dict_near = etype_to_nearest_points(target, source; maxdist) - num_basis = length(p) + num_basis = length(basis) b = Dense{T}(undef, length(source), num_basis) γ₀B = Dense{T}(undef, length(boundary), num_basis) γ₁B = Dense{T}(undef, length(boundary), num_basis) - b = [f(q) for q in source, f in p] + b = [basis[m].source(q) for q in source, m in eachindex(basis)] for k in 1:num_basis, j in 1:length(boundary) - γ₀B[j, k] = P[k](boundary[j]) - γ₁B[j, k] = γ₁P[k](boundary[j]) + γ₀B[j, k] = basis[k].solution(boundary[j]) + γ₁B[j, k] = basis[k].neumann_trace(boundary[j]) end T = eltype(γ₀B) Θ = zeros(T, num_target, num_basis) @@ -70,7 +70,7 @@ function vdim_correction( @views mul!(Θ[:, n], Dop, γ₀B[:, n], -1, 1) @views mul!(Θ[:, n], Vop, b[:, n], -1, 1) for i in 1:num_target - Θ[i, n] += green_multiplier[i] * P[n](target[i]) + Θ[i, n] += green_multiplier[i] * basis[n].solution(target[i]) end end # compute sparse correction @@ -95,7 +95,7 @@ function vdim_correction( jglob = @view qtags[:, n] # Fill the interpolation matrix for k in 1:nq, m in 1:num_basis - L_arr[m, k] = p[m](view(source, jglob)[k]) + L_arr[m, k] = basis[m].source(view(source, jglob)[k]) end F = svd(Ldata) @debug (vander_cond = max(vander_cond, cond(Ldata))) maxlog = 0 @@ -127,9 +127,7 @@ function vdim_correction( end function _vdim_auxiliary_quantities( - p, - P, - γ₁P, + basis, X, Y::Quadrature, Γ::Quadrature, @@ -138,11 +136,11 @@ function _vdim_auxiliary_quantities( Dop, Vop, ) - num_basis = length(p) + num_basis = length(basis) num_targets = length(X) - b = [f(q) for q in Y, f in p] - γ₀B = [f(q) for q in Γ, f in P] - γ₁B = [f(q) for q in Γ, f in γ₁P] + b = [basis[m].source(q) for q in Y, m in eachindex(basis)] + γ₀B = [basis[m].solution(q) for q in Γ, m in eachindex(basis)] + γ₁B = [basis[m].neumann_trace(q) for q in Γ, m in eachindex(basis)] T = eltype(γ₀B) Θ = zeros(T, num_targets, num_basis) # Compute Θ <-- S * γ₁B - D * γ₀B - V * b + σ * B(x) using in-place matvec @@ -151,114 +149,134 @@ function _vdim_auxiliary_quantities( @views mul!(Θ[:, n], Dop, γ₀B[:, n], -1, 1) @views mul!(Θ[:, n], Vop, b[:, n], -1, 1) for i in 1:num_targets - Θ[i, n] += μ[i] * P[n](X[i]) + Θ[i, n] += μ[i] * basis[n].solution(X[i]) end end return Θ end """ - polynomial_solutions_vdim(op, order) + polynomial_solutions_vdim(op, order, [T]) -For every monomial term `pₙ` of degree at most `order`, compute a polynomial `Pₙ` such that -`ℒ[Pₙ] = pₙ`, where `ℒ` is the differential operator associated with `op`. This function -returns `{pₙ,Pₙ,γ₁Pₙ}`, where `γ₁Pₙ` is the generalized Neumann trace of `Pₙ`, and `Iₙ` -is the multi-index associated with `pₙ`. +Build a basis of polynomial solutions for the VDIM method. + +For every monomial `pₙ` of degree at most `order`, computes a polynomial solution `Pₙ` +satisfying `ℒ[Pₙ] = pₙ`, where `ℒ` is the differential operator of `op`. + +Returns a vector of named tuples `(source = pₙ, solution = Pₙ, neumann_trace = γ₁Pₙ)`. """ function polynomial_solutions_vdim( op::AbstractDifferentialOperator{N}, order::Integer, ::Type{T} = default_kernel_eltype(op), ) where {N, T} - # create empty arrays to store the monomials, solutions, and traces. - monomials = Vector{Polynomial{N, T}}() - solutions = Vector{Polynomial{N, T}}() - for I in Iterators.product(ntuple(i -> 0:order, N)...) - sum(I) > order && continue + indices = [I for I in Iterators.product(ntuple(i -> 0:order, N)...) if sum(I) <= order] + return map(indices) do I monomial = Polynomial(I => one(T)) - sol = monomial_solution(op, monomial) - push!(monomials, monomial) - push!(solutions, sol) + P, γ₁P = basis_from_monomial(op, monomial) + (source = monomial, solution = P, neumann_trace = γ₁P) end - neumann_traces = map(solutions) do sol - neumann_trace(op, sol) - end - return monomials, solutions, neumann_traces end """ - monomial_solution(op::AbstractDifferentialOperator{N}, p::Polynomial{N}) - -Compute a polynomial solution `P` to the equation `ℒ[P] = p`, where + basis_from_monomial(op, monomial) -> (solution, neumann_trace) -- `ℒ` is the differential operator associated with `op`, -- `p(x) = x₁^α₁ * x₂^α₂ * ... * x_N^α_N * I` is the monomial of multi-index `α`, -- `I` is the `σ × σ` identity tensor, with `σ` the range dimension of `op` (i.e., number of - components of the solution). +Compute a polynomial solution `P` to `ℒ[P] = monomial` and its Neumann trace `γ₁P`. -Returns `(p, P, γ₁P)`, where `γ₁P` is the generalized Neumann trace of `P`. - -Both `p` and `P` are given as `ElementaryPDESolutions.Polynomial` objects, while +Each operator implements this to handle its specific PDE structure, including any +auxiliary fields (e.g., pressure for Stokes) needed to compute the Neumann trace. """ -monomial_solution(::Laplace{N}, p::Polynomial{N}) where {N} = ElementaryPDESolutions.solve_laplace(p) +function basis_from_monomial end -function neumann_trace(::Laplace{N}, p::Polynomial{N, T}) where {N, T} - ∇P = ElementaryPDESolutions.gradient(p) - γ₁P = (q) -> dot(normal(q), ∇P(coords(q))) - return γ₁P -end +# Laplace -function monomial_solution(op::Helmholtz{N}, p::Polynomial{N}) where {N} - return ElementaryPDESolutions.solve_helmholtz(p, op.k^2) +function basis_from_monomial(::Laplace{N}, monomial::Polynomial{N, T}) where {N, T} + P = ElementaryPDESolutions.solve_laplace(monomial) + ∇P = ElementaryPDESolutions.gradient(P) + γ₁P = q -> dot(normal(q), ∇P(coords(q))) + return P, γ₁P end -function neumann_trace(::Helmholtz{N}, p::Polynomial{N, T}) where {N, T} - ∇P = ElementaryPDESolutions.gradient(p) - γ₁P = (q) -> dot(normal(q), ∇P(coords(q))) - return γ₁P +# Helmholtz + +function basis_from_monomial(op::Helmholtz{N}, monomial::Polynomial{N, T}) where {N, T} + P = ElementaryPDESolutions.solve_helmholtz(monomial, op.k^2) + ∇P = ElementaryPDESolutions.gradient(P) + γ₁P = q -> dot(normal(q), ∇P(coords(q))) + return P, γ₁P end -function monomial_solution(op::Elastostatic{N}, p::Polynomial{N, T}) where {N, T} - @assert T <: StaticMatrix - @assert size(T) == (N, N) - # extract exponent of p, and make sure it is a monomial - ord2coef = p.order2coeff - @assert length(ord2coef) == 1 "Input polynomial must be a monomial" +# Elastostatic + +function basis_from_monomial(op::Elastostatic{N}, monomial::Polynomial{N, T}) where {N, T} + @assert T <: StaticMatrix && size(T) == (N, N) + S = eltype(T) + ord2coef = monomial.order2coeff + @assert length(ord2coef) == 1 "Input must be a monomial" coef = first(values(ord2coef)) - val = first(keys(ord2coef)) - # extract material parameters - μ = op.μ - λ = op.λ + idx = first(keys(ord2coef)) + μ, λ = op.μ, op.λ ν = λ / (2 * (λ + μ)) - # compute solution for each columns of the tensor - ptuple = ntuple(N) do n - p = ntuple(d -> Polynomial(val => coef[d, n]), N) - ElementaryPDESolutions.solve_elastostatic(p; μ, ν) + # Solve for each column of the tensor + sol_tuple = ntuple(N) do n + p = ntuple(d -> Polynomial(idx => coef[d, n]), N) + u = ElementaryPDESolutions.solve_elastostatic(p; μ, ν) + ntuple(d -> convert(Polynomial{N, S}, u[d]), N) end - p = flatten_polynomial_ntuple(ptuple) - return p -end - -function neumann_trace( - op::Elastostatic{N}, - P::Polynomial{N, T}, - ) where {N, T} - μ = op.μ - λ = op.λ + P = flatten_polynomial_ntuple(sol_tuple) ∇P = ElementaryPDESolutions.gradient(P) - γ₁P = (q) -> begin - ν = normal(q) + # Neumann trace: traction vector + γ₁P = q -> begin + n = normal(q) x = coords(q) M = ∇P(x) # M[j] = ∂P/∂xⱼ cols = svector(N) do m - # Build gradient of m-th column: (∇uₘ)[:, j] = M[j][:, m] gradu = hcat(ntuple(j -> M[j][:, m], N)...) divu = tr(gradu) - return λ * divu * ν + μ * (gradu + gradu') * ν + λ * divu * n + μ * (gradu + gradu') * n + end + reduce(hcat, cols) + end + return P, γ₁P +end + +# Stokes + +function basis_from_monomial(op::Stokes{N}, monomial::Polynomial{N, T}) where {N, T} + @assert T <: StaticMatrix && size(T) == (N, N) + S = eltype(T) + ord2coef = monomial.order2coeff + @assert length(ord2coef) == 1 "Input must be a monomial" + coef = first(values(ord2coef)) + idx = first(keys(ord2coef)) + μ = op.μ + # Solve for each column: velocity and pressure + # Note: ElementaryPDESolutions uses μΔu - ∇p = f, while Inti uses -μΔu + ∇p = f, + # so we negate the source term + solutions = ntuple(N) do n + f = ntuple(d -> Polynomial(idx => -coef[d, n]), N) + u, p = ElementaryPDESolutions.solve_stokes(f; μ) + vel = ntuple(d -> convert(Polynomial{N, S}, u[d]), N) + pres = convert(Polynomial{N, S}, p) + (velocity = vel, pressure = pres) + end + velocities = ntuple(n -> solutions[n].velocity, N) + pressures = ntuple(n -> solutions[n].pressure, N) + U = flatten_polynomial_ntuple(velocities) + ∇U = ElementaryPDESolutions.gradient(U) + # Neumann trace: traction using velocity gradient and pressure + γ₁U = q -> begin + n = normal(q) + x = coords(q) + M = ∇U(x) # M[j] = ∂U/∂xⱼ + cols = svector(N) do m + gradu = hcat(ntuple(j -> M[j][:, m], N)...) + p_val = pressures[m](x) + -p_val * n + μ * (gradu + gradu') * n end reduce(hcat, cols) end - return γ₁P + return U, γ₁U end function flatten_polynomial_ntuple(P::NTuple{N, NTuple{N, Polynomial{DIM, T}}}) where {N, DIM, T <: Number} diff --git a/test/volume_potential_test.jl b/test/volume_potential_test.jl index 6b81f1f4..8c965ff6 100644 --- a/test/volume_potential_test.jl +++ b/test/volume_potential_test.jl @@ -82,25 +82,25 @@ function test_volume_potential(op, Ω, Γ, msh; interpolation_order = 2, bdry_qo ) # Get polynomial solutions for this operator - monomials, dir_traces, neumann_traces = Inti.polynomial_solutions_vdim(op, interpolation_order) + basis = Inti.polynomial_solutions_vdim(op, interpolation_order) # Test with each polynomial basis function errors_uncorrected = Float64[] errors_corrected = Float64[] - for idx in 1:length(monomials) + for idx in 1:length(basis) # For vector-valued operators, multiply by a constant vector if Inti.default_density_eltype(op) <: SVector N = Inti.ambient_dimension(Ωₕ) c = SVector(ntuple(i -> rand(), N)...) - f = (q) -> monomials[idx](q) * c - u = (q) -> dir_traces[idx](q) * c - t = (q) -> neumann_traces[idx](q) * c + f = (q) -> basis[idx].source(q) * c + u = (q) -> basis[idx].solution(q) * c + t = (q) -> basis[idx].neumann_trace(q) * c else # Scalar case - f = (q) -> monomials[idx](q) - u = (q) -> dir_traces[idx](q) - t = (q) -> neumann_traces[idx](q) + f = (q) -> basis[idx].source(q) + u = (q) -> basis[idx].solution(q) + t = (q) -> basis[idx].neumann_trace(q) end # Evaluate on quadrature nodes @@ -149,6 +149,7 @@ end ("2D Laplace", Inti.Laplace(; dim = 2)), ("2D Helmholtz", Inti.Helmholtz(; k = 0.7, dim = 2)), ("2D Elastostatic", Inti.Elastostatic(; μ = 0.8, λ = 1.3, dim = 2)), + ("2D Stokes", Inti.Stokes(; μ = 1.0, dim = 2)), ] for (name, op) in operators @@ -181,6 +182,7 @@ end ("3D Laplace", Inti.Laplace(; dim = 3)), ("3D Helmholtz", Inti.Helmholtz(; k = 1.2, dim = 3)), ("3D Elastostatic", Inti.Elastostatic(; μ = 1.1, λ = 0.9, dim = 3)), + ("3D Stokes", Inti.Stokes(; μ = 1.0, dim = 3)), ] for (name, op) in operators From 66f73383f870eff63bfc73160d0d4a433b146a2e Mon Sep 17 00:00:00 2001 From: maltezfaria Date: Wed, 21 Jan 2026 18:44:47 +0100 Subject: [PATCH 15/20] fix adjoint/transpose issues and make test more stringent --- src/vdim.jl | 4 ++-- test/volume_potential_test.jl | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/vdim.jl b/src/vdim.jl index 72a19300..0a9ae9c2 100644 --- a/src/vdim.jl +++ b/src/vdim.jl @@ -103,14 +103,14 @@ function vdim_correction( @debug (vander_norm = max(vander_norm, norm(Ldata))) maxlog = 0 # correct each target near the current element for i in near_list[n] - b_arr .= @views Θ[i:i, :]' + b_arr .= @views transpose(Θ[i:i, :]) @debug (rhs_norm = max(rhs_norm, norm(bdata))) maxlog = 0 ldiv!(weidata, F, bdata) @debug (res_norm = max(res_norm, norm(Ldata * weidata - bdata))) maxlog = 0 for k in 1:nq push!(Is, i) push!(Js, jglob[k]) - push!(Vs, wei_arr[k, 1]) + push!(Vs, transpose(wei_arr[k])) end end end diff --git a/test/volume_potential_test.jl b/test/volume_potential_test.jl index 8c965ff6..3118c2c0 100644 --- a/test/volume_potential_test.jl +++ b/test/volume_potential_test.jl @@ -24,7 +24,7 @@ using Gmsh Random.seed!(1) ## Test parameters -rtol = 1.0e-2 # relative tolerance for volume potential tests +rtol = 1.0e-10 # relative tolerance for volume potential tests meshsize = 0.4 meshorder = 1 bdry_qorder = 5 From d5c96cbd0b4fe20d1f8b17db368160a6db8bbf0a Mon Sep 17 00:00:00 2001 From: maltezfaria Date: Mon, 26 Jan 2026 21:28:36 +0100 Subject: [PATCH 16/20] fix vector valued `vdim` to work with FMM3D --- src/vdim.jl | 37 ++++++++++++++++++++++++++++++++----- test/fmm3d_test.jl | 44 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 76 insertions(+), 5 deletions(-) diff --git a/src/vdim.jl b/src/vdim.jl index 0a9ae9c2..c85dcb5c 100644 --- a/src/vdim.jl +++ b/src/vdim.jl @@ -57,18 +57,45 @@ function vdim_correction( b = Dense{T}(undef, length(source), num_basis) γ₀B = Dense{T}(undef, length(boundary), num_basis) γ₁B = Dense{T}(undef, length(boundary), num_basis) - b = [basis[m].source(q) for q in source, m in eachindex(basis)] + for k in 1:num_basis, j in 1:length(source) + b[j, k] = basis[k].source(source[j]) + end for k in 1:num_basis, j in 1:length(boundary) γ₀B[j, k] = basis[k].solution(boundary[j]) γ₁B[j, k] = basis[k].neumann_trace(boundary[j]) end T = eltype(γ₀B) - Θ = zeros(T, num_target, num_basis) + Θ = Dense{T}(undef, num_target, num_basis) + fill!(Θ, zero(T)) # Compute Θ <-- S * γ₁B - D * γ₀B - V * b + σ * B(x) using in-place matvec + if Dense <: Array || (Sop isa BlockArray && Dop isa BlockArray && Vop isa BlockArray) + for n in 1:num_basis + @views mul!(Θ[:, n], Sop, γ₁B[:, n]) + @views mul!(Θ[:, n], Dop, γ₀B[:, n], -1, 1) + @views mul!(Θ[:, n], Vop, b[:, n], -1, 1) + end + else + # For vector-valued problems with FMM (LinearMap operators), we need to + # perform multiplication column-by-column since FMM expects vector densities + # (SVector) not matrix densities (SMatrix). See bdim.jl for similar handling. + P, Q = size(T) + S = eltype(T) + Θ_data = parent(Θ)::Matrix + γ₀B_data = parent(γ₀B)::Matrix + γ₁B_data = parent(γ₁B)::Matrix + b_data = parent(b)::Matrix + for k in 1:size(Θ_data, 2) + y = reinterpret(SVector{P, S}, @view Θ_data[:, k]) + x = reinterpret(SVector{Q, S}, @view γ₁B_data[:, k]) + mul!(y, Sop, x) + x = reinterpret(SVector{Q, S}, @view γ₀B_data[:, k]) + mul!(y, Dop, x, -1, 1) + x = reinterpret(SVector{Q, S}, @view b_data[:, k]) + mul!(y, Vop, x, -1, 1) + end + end + # Add σ * B(x) term for n in 1:num_basis - @views mul!(Θ[:, n], Sop, γ₁B[:, n]) - @views mul!(Θ[:, n], Dop, γ₀B[:, n], -1, 1) - @views mul!(Θ[:, n], Vop, b[:, n], -1, 1) for i in 1:num_target Θ[i, n] += green_multiplier[i] * basis[n].solution(target[i]) end diff --git a/test/fmm3d_test.jl b/test/fmm3d_test.jl index fc805d87..dbb01be8 100644 --- a/test/fmm3d_test.jl +++ b/test/fmm3d_test.jl @@ -4,6 +4,7 @@ using FMM3D using Gmsh using LinearAlgebra using Random +using StaticArrays include("test_utils.jl") @@ -44,3 +45,46 @@ for op in ( end end end + +# Test VDIM correction with FMM for Stokes (vector-valued PDE) +@testset "VDIM + FMM for Stokes 3D" begin + op = Inti.Stokes(; μ = 1.0, dim = 3) + + # Create a coarse mesh for this test + Ω_coarse, msh_coarse = gmsh_ball(; center = [0.0, 0.0, 0.0], radius = 1.0, meshsize = 0.4) + Γ_coarse = Inti.external_boundary(Ω_coarse) + Γₕ_quad = Inti.Quadrature(view(msh_coarse, Γ_coarse); qorder = 4) + + # Create volume quadrature + Ωₕ = view(msh_coarse, Ω_coarse) + VR_qorder = Inti.Tetrahedron_VR_interpolation_order_to_quadrature_order(2) + Q = Inti.VioreanuRokhlin(; domain = :tetrahedron, order = VR_qorder) + Ωₕ_quad = Inti.Quadrature(Ωₕ, Dict(E => Q for E in Inti.element_types(Ωₕ))) + + # Build operators with FMM + S, D = Inti.single_double_layer(; + op, target = Ωₕ_quad, source = Γₕ_quad, + compression = (method = :fmm, tol = 1.0e-10), + correction = (method = :dim, maxdist = 0.5, target_location = :inside) + ) + V = Inti.volume_potential(; + op, target = Ωₕ_quad, source = Ωₕ_quad, + compression = (method = :fmm, tol = 1.0e-10), correction = (method = :none,) + ) + δV = Inti.vdim_correction( + op, Ωₕ_quad, Ωₕ_quad, Γₕ_quad, S, D, V; + green_multiplier = -ones(length(Ωₕ_quad)), interpolation_order = 2 + ) + + # Test Green's identity with polynomial solution + basis = Inti.polynomial_solutions_vdim(op, 2) + c = SVector(1.0, 2.0, 3.0) + u_d = [basis[1].solution(q) * c for q in Ωₕ_quad] + u_b = [basis[1].solution(q) * c for q in Γₕ_quad] + du_b = [basis[1].neumann_trace(q) * c for q in Γₕ_quad] + f_d = [basis[1].source(q) * c for q in Ωₕ_quad] + + vref = -u_d - D * u_b + S * du_b + vapprox = V * f_d + δV * f_d + @test norm(vref - vapprox, Inf) < 1.0e-10 +end From 3b6eeca83d35a14d1cde5287706164f75e265104 Mon Sep 17 00:00:00 2001 From: maltezfaria Date: Fri, 13 Mar 2026 09:51:10 +0100 Subject: [PATCH 17/20] sign mistake --- src/vdim.jl | 4 +--- test/volume_potential_test.jl | 3 +-- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/src/vdim.jl b/src/vdim.jl index c85dcb5c..be4bbbe5 100644 --- a/src/vdim.jl +++ b/src/vdim.jl @@ -278,10 +278,8 @@ function basis_from_monomial(op::Stokes{N}, monomial::Polynomial{N, T}) where {N idx = first(keys(ord2coef)) μ = op.μ # Solve for each column: velocity and pressure - # Note: ElementaryPDESolutions uses μΔu - ∇p = f, while Inti uses -μΔu + ∇p = f, - # so we negate the source term solutions = ntuple(N) do n - f = ntuple(d -> Polynomial(idx => -coef[d, n]), N) + f = ntuple(d -> Polynomial(idx => coef[d, n]), N) u, p = ElementaryPDESolutions.solve_stokes(f; μ) vel = ntuple(d -> convert(Polynomial{N, S}, u[d]), N) pres = convert(Polynomial{N, S}, p) diff --git a/test/volume_potential_test.jl b/test/volume_potential_test.jl index 3118c2c0..34fd1018 100644 --- a/test/volume_potential_test.jl +++ b/test/volume_potential_test.jl @@ -51,8 +51,7 @@ function test_volume_potential(op, Ω, Γ, msh; interpolation_order = 2, bdry_qo else error("Unsupported dimension: $dim") end - dict = Dict(E => Q for E in Inti.element_types(Ωₕ)) - Ωₕ_quad = Inti.Quadrature(Ωₕ, dict) + Ωₕ_quad = Inti.Quadrature(Ωₕ, Q) Γₕ_quad = Inti.Quadrature(Γₕ; qorder = bdry_qorder) # Build boundary operators From bab1690a1aed02e3b80daa4f16a817bf875877ae Mon Sep 17 00:00:00 2001 From: maltezfaria Date: Fri, 13 Mar 2026 21:13:34 +0100 Subject: [PATCH 18/20] consistent signs internally - sign mistake of volume potential of Green identity of `vdim` - sign mistake in the polynomial solution. Negation needed due to `-\Delta` convention --- src/vdim.jl | 45 +++++++---------------------------- test/volume_potential_test.jl | 2 +- 2 files changed, 10 insertions(+), 37 deletions(-) diff --git a/src/vdim.jl b/src/vdim.jl index be4bbbe5..b4d47abb 100644 --- a/src/vdim.jl +++ b/src/vdim.jl @@ -67,12 +67,12 @@ function vdim_correction( T = eltype(γ₀B) Θ = Dense{T}(undef, num_target, num_basis) fill!(Θ, zero(T)) - # Compute Θ <-- S * γ₁B - D * γ₀B - V * b + σ * B(x) using in-place matvec + # Compute Θ <-- S * γ₁B - D * γ₀B + V * b + σ * B(x) using in-place matvec if Dense <: Array || (Sop isa BlockArray && Dop isa BlockArray && Vop isa BlockArray) for n in 1:num_basis @views mul!(Θ[:, n], Sop, γ₁B[:, n]) @views mul!(Θ[:, n], Dop, γ₀B[:, n], -1, 1) - @views mul!(Θ[:, n], Vop, b[:, n], -1, 1) + @views mul!(Θ[:, n], Vop, b[:, n], 1, 1) end else # For vector-valued problems with FMM (LinearMap operators), we need to @@ -91,7 +91,7 @@ function vdim_correction( x = reinterpret(SVector{Q, S}, @view γ₀B_data[:, k]) mul!(y, Dop, x, -1, 1) x = reinterpret(SVector{Q, S}, @view b_data[:, k]) - mul!(y, Vop, x, -1, 1) + mul!(y, Vop, x, 1, 1) end end # Add σ * B(x) term @@ -137,7 +137,7 @@ function vdim_correction( for k in 1:nq push!(Is, i) push!(Js, jglob[k]) - push!(Vs, transpose(wei_arr[k])) + push!(Vs, -transpose(wei_arr[k])) end end end @@ -153,35 +153,6 @@ function vdim_correction( return δV end -function _vdim_auxiliary_quantities( - basis, - X, - Y::Quadrature, - Γ::Quadrature, - μ, - Sop, - Dop, - Vop, - ) - num_basis = length(basis) - num_targets = length(X) - b = [basis[m].source(q) for q in Y, m in eachindex(basis)] - γ₀B = [basis[m].solution(q) for q in Γ, m in eachindex(basis)] - γ₁B = [basis[m].neumann_trace(q) for q in Γ, m in eachindex(basis)] - T = eltype(γ₀B) - Θ = zeros(T, 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], Sop, γ₁B[:, n]) - @views mul!(Θ[:, n], Dop, γ₀B[:, n], -1, 1) - @views mul!(Θ[:, n], Vop, b[:, n], -1, 1) - for i in 1:num_targets - Θ[i, n] += μ[i] * basis[n].solution(X[i]) - end - end - return Θ -end - """ polynomial_solutions_vdim(op, order, [T]) @@ -215,10 +186,10 @@ auxiliary fields (e.g., pressure for Stokes) needed to compute the Neumann trace """ function basis_from_monomial end -# Laplace +# Laplace function basis_from_monomial(::Laplace{N}, monomial::Polynomial{N, T}) where {N, T} - P = ElementaryPDESolutions.solve_laplace(monomial) + P = ElementaryPDESolutions.solve_laplace(-monomial) ∇P = ElementaryPDESolutions.gradient(P) γ₁P = q -> dot(normal(q), ∇P(coords(q))) return P, γ₁P @@ -227,7 +198,7 @@ end # Helmholtz function basis_from_monomial(op::Helmholtz{N}, monomial::Polynomial{N, T}) where {N, T} - P = ElementaryPDESolutions.solve_helmholtz(monomial, op.k^2) + P = ElementaryPDESolutions.solve_helmholtz(-monomial, op.k^2) ∇P = ElementaryPDESolutions.gradient(P) γ₁P = q -> dot(normal(q), ∇P(coords(q))) return P, γ₁P @@ -236,6 +207,7 @@ end # Elastostatic function basis_from_monomial(op::Elastostatic{N}, monomial::Polynomial{N, T}) where {N, T} + monomial = -monomial @assert T <: StaticMatrix && size(T) == (N, N) S = eltype(T) ord2coef = monomial.order2coeff @@ -270,6 +242,7 @@ end # Stokes function basis_from_monomial(op::Stokes{N}, monomial::Polynomial{N, T}) where {N, T} + monomial = -monomial @assert T <: StaticMatrix && size(T) == (N, N) S = eltype(T) ord2coef = monomial.order2coeff diff --git a/test/volume_potential_test.jl b/test/volume_potential_test.jl index 34fd1018..c86c171f 100644 --- a/test/volume_potential_test.jl +++ b/test/volume_potential_test.jl @@ -111,7 +111,7 @@ function test_volume_potential(op, Ω, Γ, msh; interpolation_order = 2, bdry_qo # Compute reference solution: -u - D*u_b + S*du_b # This comes from Green's representation: u = S*t - D*u + V*f # So V*f = u - S*t + D*u, and we test -V*f = S*t - D*u - u - vref = -u_d - D_b2d * u_b + S_b2d * du_b + vref = u_d + D_b2d * u_b - S_b2d * du_b # Compute uncorrected approximation vapprox = V_d2d * f_d From 3a9a22e56c1e7fdf2a23492b8e346961dedc8318 Mon Sep 17 00:00:00 2001 From: maltezfaria Date: Fri, 13 Mar 2026 21:15:44 +0100 Subject: [PATCH 19/20] improve debugging messages --- src/Inti.jl | 1 + src/bdim.jl | 18 ++++++++++++------ src/utils.jl | 9 +++++++++ src/vdim.jl | 16 +++++++++++----- 4 files changed, 33 insertions(+), 11 deletions(-) diff --git a/src/Inti.jl b/src/Inti.jl index 14bcf43e..e6f99312 100644 --- a/src/Inti.jl +++ b/src/Inti.jl @@ -18,6 +18,7 @@ using DataStructures using ForwardDiff using LinearAlgebra using LinearMaps +using Logging using NearestNeighbors using Pkg using Printf diff --git a/src/bdim.jl b/src/bdim.jl index 213315ba..1f73db71 100644 --- a/src/bdim.jl +++ b/src/bdim.jl @@ -66,6 +66,8 @@ function bdim_correction( N = ambient_dimension(source) @assert eltype(Dop) == T "eltype of S and D must match" m, n = length(target), length(source) + # check if we are in debug mode to avoid expensive computations + do_debug = debug_mode() if isnothing(filter_target_params) dict_near = etype_to_nearest_points(target, source; maxdist) num_trgs = m @@ -161,17 +163,21 @@ function bdim_correction( # TODO: get ride of all this transposing mumble jumble by assembling # the matrix in the correct orientation in the first place F = qr!(transpose(Mdata)) - @debug (imat_cond = max(cond(Mdata), imat_cond)) maxlog = 0 - @debug (imat_norm = max(norm(Mdata), imat_norm)) maxlog = 0 + if do_debug + imat_cond = max(cond(Mdata), imat_cond) + imat_norm = max(norm(Mdata), imat_norm) + end for i in near_list[n] j = glob_loc_near_trgs[i] Θi .= Θ[j:j, :] - @debug (rhs_norm = max(rhs_norm, norm(Θidata))) maxlog = 0 + if do_debug + rhs_norm = max(rhs_norm, norm(Θidata)) + end ldiv!(Wdata, F, transpose(Θidata)) - @debug ( + if do_debug res_norm = max(norm(Matrix(F) * Wdata - transpose(Θidata)), res_norm) - ) maxlog = 0 - @debug (theta_norm = max(theta_norm, norm(Wdata))) maxlog = 0 + theta_norm = max(theta_norm, norm(Wdata)) + end for k in 1:nq push!(Is, i) push!(Js, jglob[k]) diff --git a/src/utils.jl b/src/utils.jl index 355f8d22..f672b9df 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -243,6 +243,15 @@ function notimplemented() return error("not (yet) implemented") end +""" + debug_mode() + +Return `true` if the current logger's minimum enabled level is `Debug` or lower. +""" +function debug_mode() + return Logging.min_enabled_level(Logging.current_logger()) <= Logging.Debug +end + """ MultiIndex{N} diff --git a/src/vdim.jl b/src/vdim.jl index b4d47abb..9209f899 100644 --- a/src/vdim.jl +++ b/src/vdim.jl @@ -50,6 +50,8 @@ function vdim_correction( # a reasonable interpolation_order if not provided isnothing(interpolation_order) && (interpolation_order = maximum(order, values(source.etype2qrule))) + # check if we are in debug mode to avoid expensive computations + do_debug = debug_mode() # by default basis centered at origin basis = polynomial_solutions_vdim(op, interpolation_order, T) dict_near = etype_to_nearest_points(target, source; maxdist) @@ -125,15 +127,19 @@ function vdim_correction( L_arr[m, k] = basis[m].source(view(source, jglob)[k]) end F = svd(Ldata) - @debug (vander_cond = max(vander_cond, cond(Ldata))) maxlog = 0 - @debug (shift_norm = max(shift_norm, 1)) maxlog = 0 - @debug (vander_norm = max(vander_norm, norm(Ldata))) maxlog = 0 + if do_debug + vander_cond = max(vander_cond, cond(Ldata)) + shift_norm = max(shift_norm, 1) + vander_norm = max(vander_norm, norm(Ldata)) + end # correct each target near the current element for i in near_list[n] b_arr .= @views transpose(Θ[i:i, :]) - @debug (rhs_norm = max(rhs_norm, norm(bdata))) maxlog = 0 ldiv!(weidata, F, bdata) - @debug (res_norm = max(res_norm, norm(Ldata * weidata - bdata))) maxlog = 0 + if do_debug + rhs_norm = max(rhs_norm, norm(bdata)) + res_norm = max(res_norm, norm(Ldata * weidata - bdata)) + end for k in 1:nq push!(Is, i) push!(Js, jglob[k]) From d2f2d15e1ba6bef58a8dd9d178fc8fb81140d206 Mon Sep 17 00:00:00 2001 From: "Thomas G. Anderson" Date: Sat, 28 Mar 2026 00:22:54 -0500 Subject: [PATCH 20/20] Properly set the LinearMaps wrapping of HMatrix objects for layer and volume potentials --- src/api.jl | 7 ++++--- src/vdim.jl | 7 ++++++- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/api.jl b/src/api.jl index 26cbccf1..72ae5e9f 100644 --- a/src/api.jl +++ b/src/api.jl @@ -176,8 +176,9 @@ function single_double_layer(; S = axpy!(true, δS, Smat) D = axpy!(true, δD, Dmat) else - S = LinearMap(Smat) + LinearMap(δS) - D = LinearMap(Dmat) + LinearMap(δD) + T = default_density_eltype(op) + S = LinearMap{T}(Smat) + LinearMap{T}(δS) + D = LinearMap{T}(Dmat) + LinearMap{T}(δD) end elseif compression.method == :fmm S = Smat + LinearMap(δS) @@ -365,7 +366,7 @@ function volume_potential(; op, target, source::Quadrature, compression, correct if compression.method ∈ (:hmatrix, :none) # TODO: in the hmatrix case, we may want to add the correction directly # to the HMatrix so that a direct solver can be later used - V = LinearMap(Vmat) + LinearMap(δV) + V = LinearMap{default_density_eltype(op)}(Vmat) + LinearMap{default_density_eltype(op)}(δV) # if target === source # V = axpy!(true, δV, Vmat) # else diff --git a/src/vdim.jl b/src/vdim.jl index 9209f899..4c43401b 100644 --- a/src/vdim.jl +++ b/src/vdim.jl @@ -44,7 +44,12 @@ function vdim_correction( T = eltype(Vop) # determine type for dense matrices Dense = T <: SMatrix ? BlockArray : Array - @assert eltype(Dop) == eltype(Sop) == T "eltype of Sop, Dop, and Vop must match" + # FIXME? In order to wrap HMatrices for vector-valued layer *potentials* we + # use LinearMap{default_density_eltype(op)} for the layer potential maps. + # However, this breaks the following line, as e.g. instead of + # SMatrix{3,3,Float64} the eltype will be SVector{3,Float64}. Is there a way + # to keep this assert working? + #@assert eltype(Dop) == eltype(Sop) == T "eltype of Sop, Dop, and Vop must match" # figure out if we are dealing with a scalar or vector PDE num_target, num_source = length(target), length(source) # a reasonable interpolation_order if not provided