From aff0878afce8f4c538ea4f4a567d4241eda2beb1 Mon Sep 17 00:00:00 2001 From: lvchien Date: Mon, 8 Jun 2026 09:22:05 +0200 Subject: [PATCH 1/8] Muller hypersingular kernel --- src/maxwell/mwops.jl | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/maxwell/mwops.jl b/src/maxwell/mwops.jl index fc869ca66..b8c8f684c 100644 --- a/src/maxwell/mwops.jl +++ b/src/maxwell/mwops.jl @@ -14,6 +14,15 @@ struct MWSingleLayer3D{T,U} <: MaxwellOperator3D{T,U} β::U end +struct MWMuellerHyperSingular{T,U} <: MaxwellOperator3D{T,U} + gamma::T + β::U +end + +MWMuellerHyperSingular(gamma) = MWMuellerHyperSingular(gamma, -1/(gamma)) + +defaultquadstrat(op::MWMuellerHyperSingular, tfs::RTRefSpace, bfs::RTRefSpace) = DoubleNumSauterQstrat(6,7,5,5,4,3) + gamma(op::MWSingleLayer3D{Val{0}, U}) where {U} = zero(U) scalartype(op::MWSingleLayer3D{T,U}) where {T,U} = promote_type(T,U) @@ -137,6 +146,23 @@ function (igd::Integrand{<:MWSingleLayer3DReg})(x,y,f,g) end end +function (igd::Integrand{<:MWMuellerHyperSingular})(x,y,f,g) + β = igd.operator.β + γ = igd.operator.gamma + + r = cartesian(x) - cartesian(y) + R = norm(r) + γR = γ*R + + green = (expm1(-γR) + γR) / (4pi*R) + + βG = β * green + + _integrands(f,g) do fi,gj + βG * dot(fi.divergence, gj.divergence) + end +end + function (igd::Integrand{<:MWDoubleLayer3D})(x,y,f,g) From d5e6fa1674d696d66f6b65874c59d5855fac877b Mon Sep 17 00:00:00 2001 From: lvchien Date: Thu, 2 Jul 2026 15:35:31 +0200 Subject: [PATCH 2/8] Mueller hypersingular kernel extraction --- src/maxwell/mwops.jl | 28 +++++++++++++++++++++++++--- 1 file changed, 25 insertions(+), 3 deletions(-) diff --git a/src/maxwell/mwops.jl b/src/maxwell/mwops.jl index b8c8f684c..efd1a7f3c 100644 --- a/src/maxwell/mwops.jl +++ b/src/maxwell/mwops.jl @@ -19,6 +19,7 @@ struct MWMuellerHyperSingular{T,U} <: MaxwellOperator3D{T,U} β::U end +# Maxwell hypersingular operator, with its hypersingularity removed, can be assembled using basis functions other than div-conforming ones MWMuellerHyperSingular(gamma) = MWMuellerHyperSingular(gamma, -1/(gamma)) defaultquadstrat(op::MWMuellerHyperSingular, tfs::RTRefSpace, bfs::RTRefSpace) = DoubleNumSauterQstrat(6,7,5,5,4,3) @@ -146,16 +147,17 @@ function (igd::Integrand{<:MWSingleLayer3DReg})(x,y,f,g) end end -function (igd::Integrand{<:MWMuellerHyperSingular})(x,y,f,g) +# MWMuellerHyperSingular operator, discretized and tested by RT. Integration-by-parts is performed +function (igd::Integrand{<:MWMuellerHyperSingular, <:RTRefSpace, <:RTRefSpace})(x,y,f,g) β = igd.operator.β γ = igd.operator.gamma r = cartesian(x) - cartesian(y) R = norm(r) + iR = 1/R γR = γ*R - green = (expm1(-γR) + γR) / (4pi*R) - + green = expm1(-γR) * (i4pi * iR) + γ * i4pi βG = β * green _integrands(f,g) do fi,gj @@ -163,6 +165,26 @@ function (igd::Integrand{<:MWMuellerHyperSingular})(x,y,f,g) end end +# MWMuellerHyperSingular operator, discretized by RT and tested by other functions. Integration-by-parts is not performed. Instead, grad of regularized Green's function is used +function (igd::Integrand{<:MWMuellerHyperSingular, <:RefSpace, <:RTRefSpace})(x,y,f,g) + β = igd.operator.β + γ = igd.operator.gamma + + r = cartesian(x) - cartesian(y) + R = norm(r) + iR = 1/R + γR = γ*R + + gradgreen = -(γ * exp(-γR) * iR + expm1(-γR) * iR^2) * (i4pi * iR) * r + + # Minus sign is because integration-by-part is not performed + βgG = -β * gradgreen + + _integrands(f,g) do fi,gj + dot(dot(fi.value, βgG), gj.divergence) + end +end + function (igd::Integrand{<:MWDoubleLayer3D})(x,y,f,g) From 54e5a414bdd333f42af1eb87c87d0ca547842d13 Mon Sep 17 00:00:00 2001 From: lvchien Date: Thu, 2 Jul 2026 16:32:13 +0200 Subject: [PATCH 3/8] Correct the complex dot product --- src/maxwell/mwops.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/maxwell/mwops.jl b/src/maxwell/mwops.jl index efd1a7f3c..0a4b317e8 100644 --- a/src/maxwell/mwops.jl +++ b/src/maxwell/mwops.jl @@ -181,7 +181,7 @@ function (igd::Integrand{<:MWMuellerHyperSingular, <:RefSpace, <:RTRefSpace})(x, βgG = -β * gradgreen _integrands(f,g) do fi,gj - dot(dot(fi.value, βgG), gj.divergence) + dot(fi.value, βgG*gj.divergence) end end From ad9a382f715edc465120b4582639834fb27794b9 Mon Sep 17 00:00:00 2001 From: lvchien Date: Fri, 7 Aug 2026 13:56:13 +0200 Subject: [PATCH 4/8] Change name to MWStaticExtractedHyperSingular --- src/maxwell/mwops.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/maxwell/mwops.jl b/src/maxwell/mwops.jl index 0a4b317e8..e6d284c43 100644 --- a/src/maxwell/mwops.jl +++ b/src/maxwell/mwops.jl @@ -14,15 +14,15 @@ struct MWSingleLayer3D{T,U} <: MaxwellOperator3D{T,U} β::U end -struct MWMuellerHyperSingular{T,U} <: MaxwellOperator3D{T,U} +struct MWStaticExtractedHyperSingular{T,U} <: MaxwellOperator3D{T,U} gamma::T β::U end -# Maxwell hypersingular operator, with its hypersingularity removed, can be assembled using basis functions other than div-conforming ones -MWMuellerHyperSingular(gamma) = MWMuellerHyperSingular(gamma, -1/(gamma)) +# Maxwell hypersingular operator, with its static kernel removed, can be assembled using basis functions other than div-conforming ones +MWStaticExtractedHyperSingular(gamma) = MWStaticExtractedHyperSingular(gamma, -1/(gamma)) -defaultquadstrat(op::MWMuellerHyperSingular, tfs::RTRefSpace, bfs::RTRefSpace) = DoubleNumSauterQstrat(6,7,5,5,4,3) +defaultquadstrat(op::MWStaticExtractedHyperSingular, tfs::RTRefSpace, bfs::RTRefSpace) = DoubleNumSauterQstrat(6,7,5,5,4,3) gamma(op::MWSingleLayer3D{Val{0}, U}) where {U} = zero(U) @@ -147,8 +147,8 @@ function (igd::Integrand{<:MWSingleLayer3DReg})(x,y,f,g) end end -# MWMuellerHyperSingular operator, discretized and tested by RT. Integration-by-parts is performed -function (igd::Integrand{<:MWMuellerHyperSingular, <:RTRefSpace, <:RTRefSpace})(x,y,f,g) +# MWStaticExtractedHyperSingular operator, discretized and tested by RT. Integration-by-parts is performed +function (igd::Integrand{<:MWStaticExtractedHyperSingular, <:RTRefSpace, <:RTRefSpace})(x,y,f,g) β = igd.operator.β γ = igd.operator.gamma @@ -165,8 +165,8 @@ function (igd::Integrand{<:MWMuellerHyperSingular, <:RTRefSpace, <:RTRefSpace})( end end -# MWMuellerHyperSingular operator, discretized by RT and tested by other functions. Integration-by-parts is not performed. Instead, grad of regularized Green's function is used -function (igd::Integrand{<:MWMuellerHyperSingular, <:RefSpace, <:RTRefSpace})(x,y,f,g) +# MWStaticExtractedHyperSingular operator, discretized by RT and tested by other functions. Integration-by-parts is not performed. Instead, grad of regularized Green's function is used +function (igd::Integrand{<:MWStaticExtractedHyperSingular, <:RefSpace, <:RTRefSpace})(x,y,f,g) β = igd.operator.β γ = igd.operator.gamma @@ -175,7 +175,7 @@ function (igd::Integrand{<:MWMuellerHyperSingular, <:RefSpace, <:RTRefSpace})(x, iR = 1/R γR = γ*R - gradgreen = -(γ * exp(-γR) * iR + expm1(-γR) * iR^2) * (i4pi * iR) * r + gradgreen = -(expm1(-γR) * (1 + γR) + γR) * (i4pi * iR^3) * r # Minus sign is because integration-by-part is not performed βgG = -β * gradgreen From cb84d4d551f81d480e3c3bd333f594a070ae2d86 Mon Sep 17 00:00:00 2001 From: lvchien Date: Mon, 28 Sep 2026 15:35:06 +0200 Subject: [PATCH 5/8] correct assemble TDFunctional and tests --- src/bases/timebasis.jl | 10 ++-- test/test_assemble_tdexcitation.jl | 83 ++++++++++++++++++++++++++++++ 2 files changed, 88 insertions(+), 5 deletions(-) create mode 100644 test/test_assemble_tdexcitation.jl diff --git a/src/bases/timebasis.jl b/src/bases/timebasis.jl index 5b087ece9..f1350faff 100644 --- a/src/bases/timebasis.jl +++ b/src/bases/timebasis.jl @@ -204,13 +204,13 @@ function assemblydata(tbf::TimeBasisFunction) els = [ simplex(point((i-1)*Δt),point(i*Δt)) for i in 1:num_cells ] for k in 1 : numfunctions(tbf) - tk = (k-1) * Δt + tk = k * Δt for i in 1 : numintervals(tbf) # Focus on interval [(i-2)Δt,(i-1)Δt] p = tbf.polys[i] q = substitute(p,t-tk) - c = k + i - 2 + c = k + i - 1 1 <= c <= num_cells || continue for d = 0 : degree(q) r = d + 1 @@ -301,10 +301,10 @@ function assemblydata(tbf::TimeBasisDelta) num_funcs = zeros(Int, num_cells, num_refs) data = fill((0,z), max_num_funcs, num_refs, num_cells) - els = [ simplex(point((i-0)*Δt),point((i+1)*Δt)) for i in 1:num_cells ] + els = [ simplex(point((i-1)*Δt),point(i*Δt)) for i in 1:num_cells ] - for k in 1 : numfunctions(tbf)-1 - data[1,1,k] = (k+1,w) + for k in 1 : numfunctions(tbf) + data[1,1,k] = (k,w) end return els, AssemblyData(data) diff --git a/test/test_assemble_tdexcitation.jl b/test/test_assemble_tdexcitation.jl new file mode 100644 index 000000000..24182ba23 --- /dev/null +++ b/test/test_assemble_tdexcitation.jl @@ -0,0 +1,83 @@ +@testitem "testing TDFunctional" begin + + using CompScienceMeshes + + Γ = meshcuboid(1.0, 1.0, 1.0, 2.0; generator=:gmsh) + X = raviartthomas(Γ) + + struct ConstFunctional{T} <: Functional{T} + constant::T + end + + function(f::ConstFunctional)(r) + f.constant * [1.0, 0.0, 0.0] + end + + struct FuncXGaussian{T} <: TDFunctional{T} + functional::ConstFunctional{T} + gaussian::BEAST.Gaussian{T} + end + + function(f::FuncXGaussian)(r,t) + r = cartesian(r) + t = cartesian(t)[1] + + f.functional(r) * f.gaussian(t) + end + + Δt, Nt = 0.5, 20 + duration = 4 * Δt + delay = 10 * Δt + + func = ConstFunctional(1.0) + gaussian = creategaussian(duration, delay) + constxgaussian = FuncXGaussian(func, gaussian) + + + ### Testing with Dirac delta + δ = timebasisdelta(Δt, Nt) + A = assemble(constxgaussian, X⊗δ) + # Analytically temporal testing + Ats = assemble(func, X) * gaussian.(Δt*[1:1:Nt;])' + + @test maximum(abs.(A - Ats)) < 1e-12 + + + ### Testing with pulse functions + p = timebasiscxd0(Δt, Nt) + B = assemble(constxgaussian, X⊗p) + igaussian = integrate(gaussian) + Bts = assemble(func, X) * (igaussian.(Δt*[1:1:Nt;]) - igaussian.(Δt*[0:1:Nt-1;]))' + + @test maximum(abs.(B - Bts)) < 1e-12 + + + ### Testing with hat functions + import SpecialFunctions: erf + function integratewithh(g::BEAST.Gaussian, ts, Δt) + A = g.scaling + t0 = g.delay + w = g.width + + r = zeros(length(ts)) + + for i in ts + a = (i-1)*Δt + x = i*Δt + b = (i+1)*Δt + + ua = 4*(a-t0)/w + ux = 4*(x-t0)/w + ub = 4*(b-t0)/w + + r[i] = A/(2*Δt) * ((t0 - a)*(erf(ux) - erf(ua)) + (b - t0)*(erf(ub) - erf(ux))) + A*w/(8*√(π)*Δt)*(exp(-ua^2) + exp(-ub^2) - 2*exp(-ux^2)) + end + return r + end + + h = timebasisc0d1(Δt, Nt) + C = assemble(constxgaussian, X⊗h) + Cts = assemble(func, X) * integratewithh(gaussian, [1:1:Nt;], Δt)' + + @test maximum(abs.(C - Cts)) < 1e-12 +end From 1bf6a0bd93a83a2fd9ec24c08cd7f82a1a6392dd Mon Sep 17 00:00:00 2001 From: lvchien Date: Mon, 28 Sep 2026 15:35:36 +0200 Subject: [PATCH 6/8] correct fouriertransform and inversefouriertransform --- src/utils/specialfns.jl | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/utils/specialfns.jl b/src/utils/specialfns.jl index 18cda9708..fc4d40b89 100644 --- a/src/utils/specialfns.jl +++ b/src/utils/specialfns.jl @@ -28,16 +28,23 @@ end function fouriertransform(a::Array, dt, t0, dim=1) n = size(a,dim) dω = 2π / (n*dt) - b = fftshift(fft(a, dim), dim) * dt / sqrt(2π) ω0 = -dω * div(n,2) + ω = ω0 .+ dω .* (0:n-1) + + b = fftshift(fft(a, dim), dim) * dt / sqrt(2π) + b .*= reshape(exp.(-im .* ω .* t0), ntuple(j -> j == dim ? n : 1, ndims(a))) + b, dω, ω0 end function inversefouriertransform(a::Array, dω, ω0, dim=1) n = size(a,dim) dt = 2π/ (n*dω) - b = ifft(a,dim) * sqrt(2π) / dt + ω = ω0 .+ dω .* (0:n-1) t0 = -dt * div(n,2) + + b = ifft(ifftshift(a .* reshape(exp.(im .* ω .* t0), ntuple(j -> j == dim ? n : 1, ndims(a))), dim), dim) * sqrt(2π) / dt + b, dt, t0 end From f77483c9fbb0b68f6b49b7a244f6b39403d2f5a5 Mon Sep 17 00:00:00 2001 From: lvchien Date: Mon, 28 Sep 2026 19:11:14 +0200 Subject: [PATCH 7/8] correct the test for timebasisdelta --- test/test_tdrhs_scaling.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_tdrhs_scaling.jl b/test/test_tdrhs_scaling.jl index 5f7ea258b..211c1e694 100644 --- a/test/test_tdrhs_scaling.jl +++ b/test/test_tdrhs_scaling.jl @@ -44,12 +44,12 @@ derive(gaussian2).(taxis2)[120]/sol2 @test all(derive(gaussian1).(taxis1)/sol1 .≈ derive(gaussian2).(taxis2)/sol2) timeels1, timead1 = BEAST.assemblydata(U1) -@test timeels1[100][1][1] ≈ Δt1*100 -@test timeels1[100][2][1] ≈ Δt1*101 +@test timeels1[100][1][1] ≈ Δt1*99 +@test timeels1[100][2][1] ≈ Δt1*100 timeels2, timead2 = BEAST.assemblydata(U2) -@test timeels2[100][1][1] ≈ Δt2*100 -@test timeels2[100][2][1] ≈ Δt2*101 +@test timeels2[100][1][1] ≈ Δt2*99 +@test timeels2[100][2][1] ≈ Δt2*100 b1[120]/sol1 b2[120]/sol2 From 8b0a1a9890c14f620bb9e63626d849aad2bd2435 Mon Sep 17 00:00:00 2001 From: lvchien Date: Mon, 28 Sep 2026 19:43:38 +0200 Subject: [PATCH 8/8] fix the test --- test/test_assemble_tdexcitation.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_assemble_tdexcitation.jl b/test/test_assemble_tdexcitation.jl index 24182ba23..80e657a6b 100644 --- a/test/test_assemble_tdexcitation.jl +++ b/test/test_assemble_tdexcitation.jl @@ -5,7 +5,7 @@ Γ = meshcuboid(1.0, 1.0, 1.0, 2.0; generator=:gmsh) X = raviartthomas(Γ) - struct ConstFunctional{T} <: Functional{T} + struct ConstFunctional{T} <: BEAST.Functional{T} constant::T end @@ -13,7 +13,7 @@ f.constant * [1.0, 0.0, 0.0] end - struct FuncXGaussian{T} <: TDFunctional{T} + struct FuncXGaussian{T} <: BEAST.TDFunctional{T} functional::ConstFunctional{T} gaussian::BEAST.Gaussian{T} end