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