From dd123feaaa85b24bcd9c2330b053b0d255dfbb4e Mon Sep 17 00:00:00 2001 From: Simon Adrian Date: Thu, 4 Dec 2025 16:48:33 +0100 Subject: [PATCH] Example implementation of 1D curvilinear elements Requires an updated CMS version! Some explanations: 1) Overall, the changes to the existing code are minimal. My understanding is that the way Sauter-Schwab works these days does not necessitate a swapping of the (control) vertices, instead we work with a reference element and basically remap the barycentric coordinates. I adapted the unit test for 1D Sauter-Schwab. I do not understand why I picked SauterSchwab as a reference for itself. Now I use the Gauss-Legendre. Of course, accuracy is low, but if something would go fundamentally wrong with the normal or the reordering of vertices I would expect it to be visible. 2) My 2D Mie series suffered from early convergence thus the rewrite. 3) The changes for lagrange.jl where necessary because the number of vertices of an element does not help to assess the type in the higher order case. (3 vertices can be a triangle or a quadratic segment) --- src/bases/lagrange.jl | 13 ++- src/helmholtz2d/hh2dops.jl | 4 +- src/identityop.jl | 2 +- src/quadrature/doublenumsauterqstrat.jl | 4 +- test/test_hh2d_mie_higher_order.jl | 109 ++++++++---------------- test/test_sauterschwabints1D.jl | 104 ++++++++++++++++++---- 6 files changed, 137 insertions(+), 99 deletions(-) diff --git a/src/bases/lagrange.jl b/src/bases/lagrange.jl index cb39257b..1104674c 100644 --- a/src/bases/lagrange.jl +++ b/src/bases/lagrange.jl @@ -1194,10 +1194,17 @@ function lagrangecx(mesh::CompScienceMeshes.AbstractMesh{<:Any,2}; order) return LagrangeBasis{order,-1,NF}(mesh, fns, pos) end -function lagrangecx(mesh::CompScienceMeshes.AbstractMesh{<:Any,3}; order) +function lagrangecx(mesh::CompScienceMeshes.AbstractMesh; order) T = coordtype(mesh) - NF = binomial(2+order, 2) + + if dimension(mesh) == 1 + NF = 1 + order + elseif dimension(mesh) == 2 + NF = binomial(2+order, 2) + else + error("We do not yet support tetrahedron") + end P = vertextype(mesh) S = Shape{T} @@ -1320,7 +1327,7 @@ end -function lagrangec0(mesh::CompScienceMeshes.AbstractMesh{<:Any,3}; order) +function lagrangec0(mesh::CompScienceMeshes.Mesh{<:Any,3}; order) T = coordtype(mesh) atol = sqrt(eps(T)) diff --git a/src/helmholtz2d/hh2dops.jl b/src/helmholtz2d/hh2dops.jl index a467664e..cab47dc5 100644 --- a/src/helmholtz2d/hh2dops.jl +++ b/src/helmholtz2d/hh2dops.jl @@ -230,8 +230,8 @@ function quaddata(op::HelmholtzOperator2D, end function quadrule(op::HelmholtzOperator2D, g::LagrangeRefSpace, f::LagrangeRefSpace, - i, τ::CompScienceMeshes.Simplex{<:Any, 1}, - j, σ::CompScienceMeshes.Simplex{<:Any, 1}, + i, τ::CompScienceMeshes.AbstractSimplex{<:Any, 1}, + j, σ::CompScienceMeshes.AbstractSimplex{<:Any, 1}, qd, qs::DoubleNumSauterQstrat) hits = _numhits(τ, σ) diff --git a/src/identityop.jl b/src/identityop.jl index ee6f8cc1..baa228ea 100644 --- a/src/identityop.jl +++ b/src/identityop.jl @@ -100,7 +100,7 @@ end defaultquadstrat(::LocalOperator, ::LagrangeRefSpace{T,D1,2}, ::LagrangeRefSpace{T,D2,2}) where {T,D1,D2} = SingleNumQStrat(6) function quaddata(op::LocalOperator, g::LagrangeRefSpace{T,Deg,2} where {T,Deg}, f::LagrangeRefSpace, tels::Vector, bels::Vector, qs::SingleNumQStrat) - U = typeof(tels[1].volume) + U = typeof(volume(tels[1])) # works for straight and curvilinear charts u, w = legendre(qs.quad_rule, 0.0, 1.0) qd = [(w[i],u[i]) for i in eachindex(w)] A = _alloc_workspace(Tuple{U,U}.(qd), g, f, tels, bels) diff --git a/src/quadrature/doublenumsauterqstrat.jl b/src/quadrature/doublenumsauterqstrat.jl index af965e4b..f28a42cc 100644 --- a/src/quadrature/doublenumsauterqstrat.jl +++ b/src/quadrature/doublenumsauterqstrat.jl @@ -228,8 +228,8 @@ function _numhits(τ, σ) hits = 0 dtol = 1.0e3 * eps(T) dmin2 = floatmax(T) - for t in vertices(τ) - for s in vertices(σ) + for t in nodes(τ) + for s in nodes(σ) d2 = LinearAlgebra.norm_sqr(t-s) d = norm(t-s) dmin2 = min(dmin2, d2) diff --git a/test/test_hh2d_mie_higher_order.jl b/test/test_hh2d_mie_higher_order.jl index 83b9ef22..2a385d57 100644 --- a/test/test_hh2d_mie_higher_order.jl +++ b/test/test_hh2d_mie_higher_order.jl @@ -5,9 +5,8 @@ using LinearAlgebra using Test using SpecialFunctions - -# We check only -for order = [2; 3] +# We check only order 5 (implies four internal control vertices) +for order = [5] ε0 = 8.854187821e-12 μ0 = 4π*1e-7 @@ -21,7 +20,7 @@ for order = [2; 3] h = λ/10 a = 1.0 # radius of the scatterer - circle = CompScienceMeshes.meshcircle(a, h) + circle = CompScienceMeshes.gmshcircle(a, h; order=order) X0 = lagrangecx(circle, order=order) X1 = lagrangec0(circle, order=order) @@ -31,20 +30,15 @@ for order = [2; 3] ## Analytical solutions - dbesselj(n,x) = (besselj(n-1, x) - besselj(n+1, x)) / 2 - dhankelh2(n,x) = (hankelh2(n-1, x) - hankelh2(n+1, x)) / 2 + dbesselj(n,x) = (besselj(n-1, x) - besselj(n+1, x)) / 2 + dhankelh2(n,x) = (hankelh2(n-1, x) - hankelh2(n+1, x)) / 2 cart2polar(x,y) = SVector(sqrt(x^2 + y^2), atan(y, x)) - - - ## Analytical Solutions # TM planewave solution (based on Sec 6.4.2, Jin, Theory and Computation of Electromagnetic Fields]) function TM_pec_planewave(E0, k, a, ρ, φ) - n = 0 - # Jin (6.4.11) a_n(n) = -(1.0im)^(-n)*(besselj(n, k*a) / hankelh2(n, k*a)) valEz(n) = a_n(n) * hankelh2(n, k*ρ)*exp(im*n*φ) @@ -55,33 +49,22 @@ for order = [2; 3] retHφ = valHφ(0) retHρ = valHρ(0) - while true + tol = eps(eltype(real(retEz))) * 1e1 + n = 0 + while n < 100 n += 1 - bufEz = valEz(n) + valEz(-n) - if abs(bufEz) >= abs(retEz) * eps(eltype(real(retEz))) * 1e3 - retEz += bufEz - else - break - end - end - n=0 - while true - n += 1 + bufEz = valEz(n) + valEz(-n) bufHφ = valHφ(n) + valHφ(-n) - if abs(bufHφ) >= abs(retHφ) * eps(eltype(real(retHφ))) * 1e3 - retHφ += bufHφ - else - break - end - end - - n=0 - while true - n += 1 bufHρ = valHρ(n) + valHρ(-n) - if abs(bufHρ) >= abs(retHρ) * eps(eltype(real(retHρ))) * 1e3 - retHρ += bufHρ + + if abs(bufEz) >= abs(retEz) * tol || + abs(bufHφ) >= abs(retHφ) * tol || + abs(bufHρ) >= abs(retHρ) * tol + + bufEz != NaN ? retEz += bufEz : ~ + bufHφ != NaN ? retHφ += bufHφ : ~ + bufHρ != NaN ? retHρ += bufHρ : ~ else break end @@ -98,10 +81,8 @@ for order = [2; 3] return [TM_pec_planewave(E0, k, a, cart2polar(p[1], p[2])...)[2] for p in pts] end - # TE planewave solution function TE_pec_planewave(H0, k, a, ρ, φ) - n = 0 # Jin (6.4.19) b(n) = -(1.0im)^(-n)*dbesselj(n,k*a)/dhankelh2(n,k*a) @@ -115,33 +96,22 @@ for order = [2; 3] retEφ = valEφ(0) retEρ = valEρ(0) - while true - n += 1 - bufHz = valHz(n) + valHz(-n) - if abs(bufHz) >= abs(retHz) * eps(eltype(real(retHz))) * 1e3 - retHz += bufHz - else - break - end - end - + tol = eps(eltype(real(retHz))) * 1e1 n = 0 - while true - n += 1 - bufEφ = valEφ(n) + valEφ(-n) - if abs(bufEφ) >= abs(retEφ) * eps(eltype(real(retEφ))) * 1e3 - retEφ += bufEφ - else - break - end - end - n = 0 - while true + while n < 100 n += 1 + bufHz = valHz(n) + valHz(-n) + bufEφ = valEφ(n) + valEφ(-n) bufEρ = valEρ(n) + valEρ(-n) - if abs(bufEρ) >= abs(retEρ) * eps(eltype(real(retEρ))) * 1e3 - retEρ += bufEρ + + if abs(bufHz) >= abs(retHz) * tol || + abs(bufEφ) >= abs(retEφ) * tol || + abs(bufEρ) >= abs(retEρ) * tol + + bufHz != NaN ? retHz += bufHz : ~ + bufEφ != NaN ? retEφ += bufEφ : ~ + bufEρ != NaN ? retEρ += bufEρ : ~ else break end @@ -158,9 +128,6 @@ for order = [2; 3] return [TE_pec_planewave(H0, k, a, cart2polar(p[1],p[2])...)[2] for p in pts] end - - - ## TM Scattering # TM-EFIE system matrix @@ -185,17 +152,13 @@ for order = [2; 3] Ez_pw_sca_num = -potential(HH2DSingleLayerNear(𝒮), pts, j_TMEFIE_pw, X0; type=ComplexF64) Ez_pw_sca_ana = TM_pec_planewave_E(E0, k, a, pts) - # Will not improve until we have curvilinear elements - # But must not become worse either - # TODO: Once curvilinear elements are support: refine order of functions and - # mesh in a lock step and update these values - Ez_pw_sca_bounds = [0.002; 0.002; 0.0014] + Ez_pw_sca_bounds = [0.002; 1.6e-8; 2.3e-9; 5.0e-13; 5.1e-13] @test norm(Ez_pw_sca_ana - Ez_pw_sca_num) ./ norm(Ez_pw_sca_ana) < Ez_pw_sca_bounds[order] Ht_pw_sca_num = -potential(HH2DDoubleLayerTransposedNear(𝒟ᵀ), pts, j_TMEFIE_pw, X0; type=SVector{2,ComplexF64}) Ht_pw_sca_ana = TM_pec_planewave_H(E0, k, a, pts) - Ht_pw_sca_bounds = [0.002; 0.002; 0.0015] + Ht_pw_sca_bounds = [0.002; 1.6e-8; 2.3e-9; 5e-13; 5.1e-13] @test norm(Ref(ẑ) .× Ht_pw_sca_num - Ht_pw_sca_ana) ./ norm(Ht_pw_sca_ana) < Ht_pw_sca_bounds[order] # TM-MFIE @@ -204,11 +167,9 @@ for order = [2; 3] j_TMMFIE_pw = M_TMMFIE \ ht_pw_inc - j_pw_sca_bounds = [0.05; 0.001; 0.0009] + j_pw_sca_bounds = [0.05; 0.0002; 5e-6; 4e-7; 6.3e-9] @test norm(j_TMMFIE_pw - j_TMEFIE_pw) / norm(j_TMEFIE_pw) < j_pw_sca_bounds[order] - - ## TE Scattering I = 1.0 @@ -238,7 +199,7 @@ for order = [2; 3] Hz_pw_sca_num = potential(HH2DDoubleLayerNear(𝒟), pts, j_TEMFIE_pw, X1; type=ComplexF64) Hz_pw_sca_ana = TE_pec_planewave_H(H0, k, a, pts) - Hz_pw_sca_bounds = [0.0015; 0.0015; 0.0015] + Hz_pw_sca_bounds = [0.0015; 7.2e-7; 2.4e-9; 1.8e-10; 3.2e-10] @test norm(Hz_pw_sca_num - Hz_pw_sca_ana) /norm(Hz_pw_sca_ana) < Hz_pw_sca_bounds[order] Et_pw_inc = 1 / (im * ω * ε0) * curl(Hz_pw_inc) @@ -249,14 +210,14 @@ for order = [2; 3] j_TEEFIE_pw = M_TEEFIE \ et_pw_inc - j_pw_sca_bounds = [0.007; 0.0002; 1.19e-5] + j_pw_sca_bounds = [0.007; 0.0002; 1e-5; 2.3e-7; 1.2e-8] @test norm(j_TEEFIE_pw - j_TEMFIE_pw)/norm(j_TEMFIE_pw) < j_pw_sca_bounds[order] Et_pw_sca_num = -potential(HH2DHyperSingularNear(𝒩), pts, j_TEEFIE_pw, X1; type=SVector{2, ComplexF64}) Et_pw_sca_ana = TE_pec_planewave_E(H0, k, a, pts) # We compute the scattered Ez component (scalar) - Et_pw_sca_bounds = [0.003; 0.0016; 0.0016] + Et_pw_sca_bounds = [0.003; 1.6e-5; 1.9e-8; 3.3e-11; 5.3e-13] @test norm(Ref(ẑ) .× Et_pw_sca_num - Et_pw_sca_ana) / norm(Et_pw_sca_ana) < Et_pw_sca_bounds[order] end \ No newline at end of file diff --git a/test/test_sauterschwabints1D.jl b/test/test_sauterschwabints1D.jl index 0cdb2427..faba1392 100644 --- a/test/test_sauterschwabints1D.jl +++ b/test/test_sauterschwabints1D.jl @@ -3,11 +3,26 @@ using LinearAlgebra using BEAST, CompScienceMeshes, StaticArrays - # Float32 not working since hankelh2 returns always F64 for T in [Float64] - T = Float64 - Γto = meshsegment(T(1.0), T(0.5)) + + v1 = SVector(0.0, 0.0)*0.5 + v2 = SVector(1.0, 0.0)*0.5 + v3 = SVector(0.25, 0.0)*0.5 + v4 = SVector(0.50, 0.0)*0.5 + v5 = SVector(0.75, 0.0)*0.5 + + shift = SVector(0.5, 0.0) + + v6 = v2 + shift + v7 = v3 + shift + v8 = v4 + shift + v9 = v5 + shift + + verts = [v1, v2, v3, v4, v5, v6, v7, v8, v9] + faces = [SVector(1, 2, 3, 4, 5), SVector(2, 6, 7, 8, 9)] + + Γto = CompScienceMeshes.CurvilinearMesh(verts, faces, 4) Γt = Γto Γs = Γt @@ -21,7 +36,7 @@ for T in [Float64] Helmholtz2D.singlelayer(; wavenumber=k) ] - refstrat = BEAST.DoubleNumSauterQstrat(3,3,0,4,30,30) + refstrat = BEAST.DoubleNumQStrat(1000, 1001) for op in ops Sref = assemble(op, Xt, Xs; quadstrat=refstrat) @@ -38,7 +53,7 @@ for T in [Float64] #glrel = norm(gl - ref) / norm(ref) ssrel = norm(ss - ref) / norm(ref) - @test ssrel < 1e-14 + @test ssrel < 1e-3 end Xt = lagrangec0d1(Γt) @@ -50,7 +65,7 @@ for T in [Float64] ] 𝒮 = - refstrat = BEAST.DoubleNumSauterQstrat(3,3,0,4,30,30) + refstrat = BEAST.DoubleNumQStrat(1000, 1001) for op in ops Sref = assemble(op, Xt, Xs; quadstrat=refstrat) @@ -67,19 +82,30 @@ for T in [Float64] #glrel = norm(gl - ref) / norm(ref) ssrel = norm(ss - ref) / norm(ref) - @test ssrel < 1e-14 + @test ssrel < 1e-3 end ## Test accuracy of vertex integrations - T = Float64 - Γto = meshsegment(T(1.0), T(1.0)) - Γt = Γto - Γs = translate(Γto, SVector(1.0, 0.0)) + v1 = SVector(0.0, 0.0) + v2 = SVector(1.0, 0.0) + v3 = SVector(0.25, 0.0) + v4 = SVector(0.50, 0.0) + v5 = SVector(0.75, 0.0) + v6 = SVector(0.0, 1.0) + v7 = SVector(0.0, 0.25) + v8 = SVector(0.0, 0.50) + v9 = SVector(0.0, 0.75) - Γs = Mesh([SVector(0.0, 0.0), SVector(0.0, 1.0)], [SVector(1, 2)]) - Xt = lagrangecxd0(Γt) - Xs = lagrangecxd0(Γs) + verts = [v1, v2, v3, v4, v5, v6, v7, v8, v9] + facest = [SVector(1, 2, 3, 4, 5)] + facess = [SVector(1, 6, 7, 8, 9)] + + Γt = CompScienceMeshes.CurvilinearMesh(verts, facest, 4) + Γs = CompScienceMeshes.CurvilinearMesh(verts, facess, 4) + + Xt = lagrangecx(Γt, order=3) + Xs = lagrangecx(Γs, order=3) λ = 10 k = 2π/λ @@ -89,9 +115,50 @@ for T in [Float64] Helmholtz2D.doublelayer(; wavenumber=k) Helmholtz2D.doublelayer_transposed(; wavenumber=k) ] - 𝒮 = - refstrat = BEAST.DoubleNumSauterQstrat(3,3,0,4,30,30) + refstrat = BEAST.DoubleNumQStrat(1000, 1001) + + for op in ops + Sref = assemble(op, Xt, Xs; quadstrat=refstrat) + ref = Sref[1, 1] + + n = 25 + + #Sgl = assemble(op, Xt, Xs; quadstrat=BEAST.DoubleNumQStrat(n, n+1)) + Sss = assemble(op, Xt, Xs; quadstrat=BEAST.DoubleNumSauterQstrat(3,3,0,4,n,n)) + + #gl = Sgl[1, 1] + ss = Sss[1, 1] + + #glrel = norm(gl - ref) / norm(ref) + ssrel = norm(ss - ref) / norm(ref) + + @test ssrel < 1e-4 + + Ssrel = norm(Sss - Sref)/norm(Sref) + + @test Ssrel < 1e-4 + end + + +## + facess = [SVector(6, 1, 9, 8, 7)] + + Γs = CompScienceMeshes.CurvilinearMesh(verts, facess, 4) + + Xt = lagrangecx(Γt, order=3) + Xs = lagrangecx(Γs, order=3) + + λ = 10 + k = 2π/λ + + ops = [ + Helmholtz2D.singlelayer(; wavenumber=k) + Helmholtz2D.doublelayer(; wavenumber=k) + Helmholtz2D.doublelayer_transposed(; wavenumber=k) + ] + + refstrat = BEAST.DoubleNumQStrat(1000, 1001) for op in ops Sref = assemble(op, Xt, Xs; quadstrat=refstrat) @@ -108,6 +175,9 @@ for T in [Float64] #glrel = norm(gl - ref) / norm(ref) ssrel = norm(ss - ref) / norm(ref) - @test ssrel < 1e-14 + @test ssrel < 1e-4 + + Ssrel = norm(Sss - Sref)/norm(Sref) + @test Ssrel < 1e-4 end end \ No newline at end of file