From 901b946912047ee71adff6569be0a53e4b80a31e Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Thu, 21 May 2026 14:30:49 +0200 Subject: [PATCH 1/7] simplified quadrature structure and simd implementation --- src/helmholtz2d/hh2dops.jl | 16 +- src/helmholtz3d/nitsche.jl | 34 +- src/integralop.jl | 6 +- src/maxwell/nitsche.jl | 44 +- src/quadrature/doublenumqstrat.jl | 6 +- src/quadrature/doublenumsauterqstrat.jl | 110 ++-- .../doublenumwiltonbogaertqstrat.jl | 41 +- src/quadrature/doublenumwiltonsauterqstrat.jl | 60 +-- src/quadrature/dynamic_quadstrat.jl | 0 src/quadrature/quadstrats.jl | 49 +- .../selfsauterdnumotherwiseqstrat.jl | 21 +- src/quadrature/simddoublenum.jl | 469 ++++++++++++++++++ src/quadrature/strategies/quadstrat.jl | 2 +- src/volumeintegral/vieops.jl | 244 ++++----- src/volumeintegral/vsieops.jl | 170 +++---- test/test_assemble_ndKnd_nonconf.jl | 2 +- test/test_nonconf_quadrules.jl | 4 +- test/test_vie.jl | 2 +- 18 files changed, 905 insertions(+), 375 deletions(-) create mode 100644 src/quadrature/dynamic_quadstrat.jl create mode 100644 src/quadrature/simddoublenum.jl diff --git a/src/helmholtz2d/hh2dops.jl b/src/helmholtz2d/hh2dops.jl index 7f14d566..7602dc68 100644 --- a/src/helmholtz2d/hh2dops.jl +++ b/src/helmholtz2d/hh2dops.jl @@ -236,12 +236,12 @@ defaultquadstrat(op::HelmholtzOperator2D, tfs, bfs) = DoubleNumSauterQstrat(30,3 function quaddata(op::HelmholtzOperator2D, test_local_space::RefSpace, trial_local_space::RefSpace, - test_charts, trial_charts, qs::DoubleNumSauterQstrat) + test_charts, trial_charts, qs::SauterQStrat) T = coordtype(test_charts[1]) - tqd = quadpoints(test_local_space, test_charts, (qs.outer_rule,)) - bqd = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,)) + # tqd = quadpoints(test_local_space, test_charts, (qs.outer_rule,)) + # bqd = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,)) leg = ( convert.(NTuple{2,T},_legendre(qs.sauter_schwab_common_vert,0,1)), @@ -255,13 +255,14 @@ function quaddata(op::HelmholtzOperator2D, convert.(NTuple{2,T},BEAST.SauterSchwabQuadrature1D._MRWrules(qs.sauter_schwab_common_face,0,1)), ) - return (tpoints=tqd, bpoints=bqd, gausslegendre=leg, marokhlinwandura=mrw) + return (gausslegendre=leg, marokhlinwandura=mrw, + nestedqd = quaddata(op, test_local_space, trial_local_space, test_charts, trial_charts, qs.nested_strat)) end function quadrule(op::HelmholtzOperator2D, g::LagrangeRefSpace, f::LagrangeRefSpace, i, τ::CompScienceMeshes.Simplex{<:Any, 1}, j, σ::CompScienceMeshes.Simplex{<:Any, 1}, - qd, qs::DoubleNumSauterQstrat) + qd, qs::SauterQStrat) hits = _numhits(τ, σ) @assert hits <= 2 @@ -269,8 +270,5 @@ function quadrule(op::HelmholtzOperator2D, g::LagrangeRefSpace, f::LagrangeRefSp hits == 2 && return BEAST.SauterSchwabQuadrature1D.CommonEdge(qd.marokhlinwandura[2], qd.gausslegendre[2]) hits == 1 && return BEAST.SauterSchwabQuadrature1D.CommonVertex(qd.marokhlinwandura[1], qd.gausslegendre[1]) - return DoubleQuadRule( - qd.tpoints[1,i], - qd.bpoints[1,j], - ) + return quadrule(op, g, f, i, τ, j, σ, qd.nestedqd, qs.nested_strat) end \ No newline at end of file diff --git a/src/helmholtz3d/nitsche.jl b/src/helmholtz3d/nitsche.jl index 285e8c89..d022e544 100644 --- a/src/helmholtz3d/nitsche.jl +++ b/src/helmholtz3d/nitsche.jl @@ -4,28 +4,28 @@ mutable struct NitscheHH3{T} <: MaxwellOperator3D{T,T} gamma::T end -defaultquadstrat(::NitscheHH3, ::LagrangeRefSpace, ::LagrangeRefSpace) = DoubleNumWiltonSauterQStrat(10,8,10,8,3,3,3,3) +defaultquadstrat(::NitscheHH3, ::LagrangeRefSpace, ::LagrangeRefSpace) = DoubleNumQStrat(10,8) -function quaddata(operator::NitscheHH3, - localtestbasis::LagrangeRefSpace, - localtrialbasis::LagrangeRefSpace, - testelements, trialelements, qs::DoubleNumWiltonSauterQStrat) +# function quaddata(operator::NitscheHH3, +# localtestbasis::LagrangeRefSpace, +# localtrialbasis::LagrangeRefSpace, +# testelements, trialelements, qs::SauterQStrat) - tqd = quadpoints(localtestbasis, testelements, (qs.outer_rule_far,)) - bqd = quadpoints(x -> localtrialbasis(x), trialelements, (qs.inner_rule_far,)) +# # tqd = quadpoints(localtestbasis, testelements, (qs.outer_rule_far,)) +# # bqd = quadpoints(localtrialbasis, trialelements, (qs.inner_rule_far,)) - #return QuadData(tqd, bqd) - return (tpoints=tqd, bpoints=bqd) -end +# #return QuadData(tqd, bqd) +# return (nestedqd = quaddata(operator, localtestbasis, localtrialbasis, testelements, trialelements, qs.nested_strat),) +# end -function quadrule(op::NitscheHH3, g::LagrangeRefSpace, f::LagrangeRefSpace, i, τ, j, σ, qd, - qs::DoubleNumWiltonSauterQStrat) +# function quadrule(op::NitscheHH3, g::LagrangeRefSpace, f::LagrangeRefSpace, i, τ, j, σ, qd, +# qs::DoubleNumWiltonSauterQStrat) - DoubleQuadRule( - qd.tpoints[1,i], - qd.bpoints[1,j] - ) -end +# DoubleQuadRule( +# qd.tpoints[1,i], +# qd.bpoints[1,j] +# ) +# end struct KernelValsMaxwell3D{T,U,P,Q} diff --git a/src/integralop.jl b/src/integralop.jl index 8e9fba70..98d25385 100644 --- a/src/integralop.jl +++ b/src/integralop.jl @@ -141,7 +141,8 @@ end @test BlockArrays.blocksize(M) == (2,2) @test BlockArrays.blocksizes(M) == [(n1,n1) (n1,n2); (n2,n1) (n2,n2)] end - +# quadrule(biop, tshapes, bshapes, p, tcell, q, bcell, qd, qdcache, quadstrat) = +# quadrule(biop, tshapes, bshapes, p, tcell, q, bcell, qd, quadstrat) function assemblechunk_body!(biop, test_space, trial_space, test_elements, test_element_ptrs, test_assembly_data, active_test_els, @@ -154,6 +155,7 @@ function assemblechunk_body!(biop, test_space, trial_space, @tasks for p in eachindex(active_test_els) @set scheduler = scheduler @local begin + qd = quadcache(biop, refspace(test_space), refspace(trial_space), test_elements, trial_elements, qd, quadstrat) zlocal = zeros(scalartype(biop, test_space, trial_space), num_tshapes, num_bshapes) tadjq = Vector{eltype(trial_assembly_data.data)}(undef, size(trial_assembly_data.data,1)) end @@ -167,7 +169,7 @@ function assemblechunk_body!(biop, test_space, trial_space, fill!(zlocal, 0) qrule = quadrule(biop, refspace(test_space), refspace(trial_space), - P, tcell, Q, bcell, qd, quadstrat) + P, tcell, Q, bcell, qd,quadstrat) momintegrals!(zlocal, biop, test_space, tptr, tcell, trial_space, bptr, bcell, qrule) for j in 1 : num_bshapes diff --git a/src/maxwell/nitsche.jl b/src/maxwell/nitsche.jl index e159c2a8..95ead5fd 100644 --- a/src/maxwell/nitsche.jl +++ b/src/maxwell/nitsche.jl @@ -12,28 +12,28 @@ mutable struct SingleLayerTrace{T} <: MaxwellOperator3D{T,T} gamma::T end -defaultquadstrat(::SingleLayerTrace, ::LagrangeRefSpace, ::LagrangeRefSpace) = DoubleNumWiltonSauterQStrat(10,8,10,8,3,3,3,3) -function quaddata(operator::SingleLayerTrace, - localtestbasis::LagrangeRefSpace, - localtrialbasis::LagrangeRefSpace, - testelements, trialelements, qs::DoubleNumWiltonSauterQStrat) - - tqd = quadpoints(localtestbasis, testelements, (qs.outer_rule_far,)) - bqd = quadpoints(localtrialbasis, trialelements, (qs.outer_rule_near,)) - - #return QuadData(tqd, bqd) - return (tpoints=tqd, bpoints=bqd) -end - -# Use numerical quadrature for now -# Note: basis integral is over triangle, test over line -function quadrule(op::SingleLayerTrace, g::LagrangeRefSpace, f::LagrangeRefSpace, i, τ, j, σ, qd, - qs::DoubleNumWiltonSauterQStrat) +defaultquadstrat(::SingleLayerTrace, ::LagrangeRefSpace, ::LagrangeRefSpace) = DoubleNumQStrat(10,8) +# function quaddata(operator::SingleLayerTrace, +# localtestbasis::LagrangeRefSpace, +# localtrialbasis::LagrangeRefSpace, +# testelements, trialelements, qs::DoubleNumWiltonSauterQStrat) + +# tqd = quadpoints(localtestbasis, testelements, (qs.outer_rule_far,)) +# bqd = quadpoints(localtrialbasis, trialelements, (qs.outer_rule_near,)) + +# #return QuadData(tqd, bqd) +# return (tpoints=tqd, bpoints=bqd) +# end + +# # Use numerical quadrature for now +# # Note: basis integral is over triangle, test over line +# function quadrule(op::SingleLayerTrace, g::LagrangeRefSpace, f::LagrangeRefSpace, i, τ, j, σ, qd, +# qs::DoubleNumWiltonSauterQStrat) - DoubleQuadRule( - qd.tpoints[1,i], - qd.bpoints[1,j] - ) -end +# DoubleQuadRule( +# qd.tpoints[1,i], +# qd.bpoints[1,j] +# ) +# end integrand(op::SingleLayerTrace, kernel, g, τ, f, σ) = f[1]*g[1]*kernel.green diff --git a/src/quadrature/doublenumqstrat.jl b/src/quadrature/doublenumqstrat.jl index 5282d619..be5ad5c5 100644 --- a/src/quadrature/doublenumqstrat.jl +++ b/src/quadrature/doublenumqstrat.jl @@ -13,7 +13,7 @@ function quaddata(operator::IntegralOperator, test_quad_data = quadpoints(local_test_basis, test_elements, (qs.outer_rule,)) trial_quad_data = quadpoints(local_trial_basis, trial_elements, (qs.inner_rule,)) - return test_quad_data, trial_quad_data + return (tpoints = test_quad_data, bpoints = trial_quad_data) end @@ -22,8 +22,8 @@ function quadrule(operator::IntegralOperator, test_id, test_element, trial_id, trial_element, quad_data, qs::DoubleNumQStrat) - test_quad_rules = quad_data[1] - trial_quad_rules = quad_data[2] + test_quad_rules = quad_data.tpoints + trial_quad_rules = quad_data.bpoints DoubleQuadRule( test_quad_rules[1,test_id], diff --git a/src/quadrature/doublenumsauterqstrat.jl b/src/quadrature/doublenumsauterqstrat.jl index af965e4b..9d5255f3 100644 --- a/src/quadrature/doublenumsauterqstrat.jl +++ b/src/quadrature/doublenumsauterqstrat.jl @@ -1,21 +1,44 @@ -struct DoubleNumSauterQstrat{R,S} <: AbstractQuadStrat - outer_rule::R - inner_rule::R +# struct DoubleNumSauterQstrat{R,S} <: AbstractQuadStrat +# outer_rule::R +# inner_rule::R +# sauter_schwab_common_tetr::S +# sauter_schwab_common_face::S +# sauter_schwab_common_edge::S +# sauter_schwab_common_vert::S +# end +struct SauterQStrat{NestedStrat,S} <: AbstractQuadStrat + nested_strat::NestedStrat sauter_schwab_common_tetr::S sauter_schwab_common_face::S sauter_schwab_common_edge::S sauter_schwab_common_vert::S end +DoubleNumSauterQstrat(a,b,c,d,e,f) = SauterQStrat(DoubleNumQStrat(a,b),c,d,e,f) +DoubleNumSauterQStrat(a,b,c,d,e,f) = SauterQStrat(DoubleNumQStrat(a,b),c,d,e,f) +# function quaddata(op::IntegralOperator, +# test_local_space::RefSpace, trial_local_space::RefSpace, +# test_charts, trial_charts, qs::DoubleNumSauterQstrat) +# T = coordtype(test_charts[1]) + +# tqd = quadpoints(test_local_space, test_charts, (qs.outer_rule,)) +# bqd = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,)) + +# leg = ( +# convert.(NTuple{2,T},_legendre(qs.sauter_schwab_common_vert,0,1)), +# convert.(NTuple{2,T},_legendre(qs.sauter_schwab_common_edge,0,1)), +# convert.(NTuple{2,T},_legendre(qs.sauter_schwab_common_face,0,1)), +# convert.(NTuple{2,T},_legendre(qs.sauter_schwab_common_tetr,0,1)), +# ) + +# return (tpoints=tqd, bpoints=bqd, gausslegendre=leg) +# end function quaddata(op::IntegralOperator, test_local_space::RefSpace, trial_local_space::RefSpace, - test_charts, trial_charts, qs::DoubleNumSauterQstrat) + test_charts, trial_charts, qs::SauterQStrat) T = coordtype(test_charts[1]) - tqd = quadpoints(test_local_space, test_charts, (qs.outer_rule,)) - bqd = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,)) - leg = ( convert.(NTuple{2,T},_legendre(qs.sauter_schwab_common_vert,0,1)), convert.(NTuple{2,T},_legendre(qs.sauter_schwab_common_edge,0,1)), @@ -23,66 +46,83 @@ function quaddata(op::IntegralOperator, convert.(NTuple{2,T},_legendre(qs.sauter_schwab_common_tetr,0,1)), ) - return (tpoints=tqd, bpoints=bqd, gausslegendre=leg) + return (gausslegendre=leg, + nestedqd = quaddata(op, test_local_space, trial_local_space, test_charts, trial_charts, qs.nested_strat)) end +# function quaddata(op::IntegralOperator, +# test_local_space::RefSpace, trial_local_space::RefSpace, +# test_charts::Vector{<:CompScienceMeshes.Simplex{<:Any, 3}}, trial_charts::Vector{<:CompScienceMeshes.Simplex{<:Any, 3}}, qs::DoubleNumSauterQstrat) + +# #The combinations of rules (6,7) and (5,7 are) BAAAADDDD +# # they result in many near singularity evaluations with any +# # resemblence of accuracy going down the drain! Simply don't! +# # (same for (5,7) btw...). +# t_qp = quadpoints(test_local_space, test_charts, (qs.outer_rule,)) +# b_qp = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,)) + +# sing_qp = (SauterSchwab3D._legendre(qs.sauter_schwab_common_vert,0,1), +# SauterSchwab3D._shunnham2D(qs.sauter_schwab_common_edge), +# SauterSchwab3D._shunnham3D(qs.sauter_schwab_common_face), +# SauterSchwab3D._shunnham4D(qs.sauter_schwab_common_tetr),) +# return (tpoints=t_qp, bpoints=b_qp, sing_qp=sing_qp) +# end function quaddata(op::IntegralOperator, test_local_space::RefSpace, trial_local_space::RefSpace, - test_charts::Vector{<:CompScienceMeshes.Simplex{<:Any, 3}}, trial_charts::Vector{<:CompScienceMeshes.Simplex{<:Any, 3}}, qs::DoubleNumSauterQstrat) + test_charts::Vector{<:CompScienceMeshes.Simplex{<:Any, 3}}, trial_charts::Vector{<:CompScienceMeshes.Simplex{<:Any, 3}}, qs::SauterQStrat) #The combinations of rules (6,7) and (5,7 are) BAAAADDDD # they result in many near singularity evaluations with any # resemblence of accuracy going down the drain! Simply don't! # (same for (5,7) btw...). - t_qp = quadpoints(test_local_space, test_charts, (qs.outer_rule,)) - b_qp = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,)) + # t_qp = quadpoints(test_local_space, test_charts, (qs.outer_rule,)) + # b_qp = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,)) sing_qp = (SauterSchwab3D._legendre(qs.sauter_schwab_common_vert,0,1), SauterSchwab3D._shunnham2D(qs.sauter_schwab_common_edge), SauterSchwab3D._shunnham3D(qs.sauter_schwab_common_face), SauterSchwab3D._shunnham4D(qs.sauter_schwab_common_tetr),) - return (tpoints=t_qp, bpoints=b_qp, sing_qp=sing_qp) + return (sing_qp=sing_qp, + nestedqd = quaddata(op, test_local_space, trial_local_space, test_charts, trial_charts, qs.nested_strat)) end function quaddata(op::IntegralOperator, test_local_space::RefSpace, trial_local_space::RefSpace, - test_charts::Vector{<:CompScienceMeshes.Simplex{<:Any, 2}}, trial_charts::Vector{<:CompScienceMeshes.Simplex{<:Any, 3}}, qs::DoubleNumSauterQstrat) + test_charts::Vector{<:CompScienceMeshes.Simplex{<:Any, 2}}, trial_charts::Vector{<:CompScienceMeshes.Simplex{<:Any, 3}}, qs::SauterQStrat) #The combinations of rules (6,7) and (5,7 are) BAAAADDDD # they result in many near singularity evaluations with any # resemblence of accuracy going down the drain! Simply don't! # (same for (5,7) btw...). - t_qp = quadpoints(test_local_space, test_charts, (qs.outer_rule,)) - b_qp = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,)) + sing_qp = (SauterSchwab3D._legendre(qs.sauter_schwab_common_vert,0,1), SauterSchwab3D._shunnham2D(qs.sauter_schwab_common_edge), SauterSchwab3D._shunnham3D(qs.sauter_schwab_common_face), SauterSchwab3D._shunnham4D(qs.sauter_schwab_common_tetr),) - return (tpoints=t_qp, bpoints=b_qp, sing_qp=sing_qp) + return (sing_qp=sing_qp, + nestedqd = quaddata(op, test_local_space, trial_local_space, test_charts, trial_charts, qs.nested_strat)) end function quaddata(op::IntegralOperator, test_local_space::RefSpace, trial_local_space::RefSpace, - test_charts::Vector{<:CompScienceMeshes.Simplex{<:Any, 3}}, trial_charts::Vector{<:CompScienceMeshes.Simplex{<:Any, 2}}, qs::DoubleNumSauterQstrat) + test_charts::Vector{<:CompScienceMeshes.Simplex{<:Any, 3}}, trial_charts::Vector{<:CompScienceMeshes.Simplex{<:Any, 2}}, qs::SauterQStrat) #The combinations of rules (6,7) and (5,7 are) BAAAADDDD # they result in many near singularity evaluations with any # resemblence of accuracy going down the drain! Simply don't! # (same for (5,7) btw...). - t_qp = quadpoints(test_local_space, test_charts, (qs.outer_rule,)) - b_qp = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,)) - sing_qp = (SauterSchwab3D._legendre(qs.sauter_schwab_common_vert,0,1), SauterSchwab3D._shunnham2D(qs.sauter_schwab_common_edge), SauterSchwab3D._shunnham3D(qs.sauter_schwab_common_face), SauterSchwab3D._shunnham4D(qs.sauter_schwab_common_tetr),) - return (tpoints=t_qp, bpoints=b_qp, sing_qp=sing_qp) + return (sing_qp=sing_qp, + nestedqd = quaddata(op, test_local_space, trial_local_space, test_charts, trial_charts, qs.nested_strat)) end function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace, i, τ::CompScienceMeshes.Simplex{<:Any, 2}, j, σ::CompScienceMeshes.Simplex{<:Any, 2}, - qd, qs::DoubleNumSauterQstrat) + qd, qs::SauterQStrat) hits = _numhits(τ, σ) @assert hits <= 3 @@ -91,9 +131,7 @@ function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace, hits == 2 && return SauterSchwabQuadrature.CommonEdge(qd.gausslegendre[2]) hits == 1 && return SauterSchwabQuadrature.CommonVertex(qd.gausslegendre[1]) - return DoubleQuadRule( - qd.tpoints[1,i], - qd.bpoints[1,j],) + return quadrule(op, g, f, i, τ, j, σ, qd.nestedqd, qs.nested_strat) end struct _TransposedStrat{A} @@ -102,21 +140,21 @@ end function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace, i, τ::CompScienceMeshes.Simplex{<:Any, 3}, - j, σ::CompScienceMeshes.Simplex{<:Any, 3}, - qd, qs::DoubleNumSauterQstrat) + j, σ::CompScienceMeshes.Simplex{<:Any, 3}, + qd, qs::SauterQStrat) qr_volume(op, g, f, i, τ, j, σ, qd, qs) end function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace, i, τ::CompScienceMeshes.Simplex{<:Any, 3}, j, σ::CompScienceMeshes.Simplex{<:Any, 2}, - qd, qs::DoubleNumSauterQstrat) + qd, qs::SauterQStrat) qr_boundary(op, g, f, i, τ, j, σ, qd, qs) end function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace, i, τ::CompScienceMeshes.Simplex{<:Any, 2}, j, σ::CompScienceMeshes.Simplex{<:Any, 3}, - qd, qs::DoubleNumSauterQstrat) + qd, qs::SauterQStrat) _TransposedStrat(qr_boundary(op, g, f, i, τ, j, σ, qd, qs)) end @@ -156,9 +194,7 @@ function qr_volume(op::IntegralOperator, g::RefSpace, f::RefSpace, i, τ, j, σ, - return DoubleQuadRule( - qd[1][1,i], - qd[2][1,j]) + return quadrule(op, g, f, i, τ, j, σ, qd.nestedqd, qs.nested_strat) end @@ -198,16 +234,14 @@ function qr_boundary(op::IntegralOperator, g::RefSpace, f::RefSpace, i, τ, j, hits == 1 && return SauterSchwab3D.CommonVertex5D_S(SauterSchwab3D.Singularity5DPoint(idx_t,idx_s),(qd.sing_qp[3],qd.sing_qp[2])) - return DoubleQuadRule( - qd[1][1,i], - qd[2][1,j]) + return quadrule(op, g, f, i, τ, j, σ, qd.nestedqd, qs.nested_strat) end function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace, i, τ::CompScienceMeshes.Quadrilateral, j, σ::CompScienceMeshes.Quadrilateral, - qd, qs::DoubleNumSauterQstrat) + qd, qs::SauterQStrat) hits = _numhits(τ, σ) @assert hits != 3 @@ -217,9 +251,7 @@ function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace, hits == 2 && return SauterSchwabQuadrature.CommonEdgeQuad(qd.gausslegendre[2]) hits == 1 && return SauterSchwabQuadrature.CommonVertexQuad(qd.gausslegendre[1]) - return DoubleQuadRule( - qd.tpoints[1,i], - qd.bpoints[1,j],) + return quadrule(op, g, f, i, τ, j, σ, qd.nestedqd, qs.nested_strat) end diff --git a/src/quadrature/doublenumwiltonbogaertqstrat.jl b/src/quadrature/doublenumwiltonbogaertqstrat.jl index f1aff891..a9df3159 100644 --- a/src/quadrature/doublenumwiltonbogaertqstrat.jl +++ b/src/quadrature/doublenumwiltonbogaertqstrat.jl @@ -1,24 +1,28 @@ -struct DoubleNumWiltonBogaertQStrat{R} <: AbstractQuadStrat - outer_rule_far::R - inner_rule_far::R - outer_rule_near::R - inner_rule_near::R +# struct DoubleNumWiltonBogaertQStrat{R} <: AbstractQuadStrat +# outer_rule_far::R +# inner_rule_far::R +# outer_rule_near::R +# inner_rule_near::R +# end +DoubleNumWiltonBogaertQStrat(a,b,c,d) = BogaertQStrat(WiltonQStrat(DoubleNumQStrat(a,b),c,d)) +struct BogaertQStrat{NestedStrat} <: AbstractQuadStrat + nested_strat::NestedStrat end function quaddata(op::IntegralOperator, test_local_space::RefSpace, trial_local_space::RefSpace, - test_charts, trial_charts, qs::DoubleNumWiltonBogaertQStrat) + test_charts, trial_charts, qs::BogaertQStrat) - T = coordtype(test_charts[1]) + # T = coordtype(test_charts[1]) - tqd = quadpoints(test_local_space, test_charts, (qs.outer_rule_far,qs.outer_rule_near)) - bqd = quadpoints(trial_local_space, trial_charts, (qs.inner_rule_far,qs.inner_rule_near)) + # tqd = quadpoints(test_local_space, test_charts, (qs.outer_rule_far,qs.outer_rule_near)) + # bqd = quadpoints(trial_local_space, trial_charts, (qs.inner_rule_far,qs.inner_rule_near)) - return (tpoints=tqd, bpoints=bqd) + return (nestedqd = quaddata(op, test_local_space, trial_local_space, test_charts, trial_charts, qs.nested_strat),) end function quadrule(op::IntegralOperator, g::RTRefSpace, f::RTRefSpace, i, τ, j, σ, qd, - qs::DoubleNumWiltonBogaertQStrat) + qs::BogaertQStrat) dtol = 1.0e3 * eps(eltype(eltype(τ.vertices))) xtol = 0.2 @@ -41,17 +45,6 @@ function quadrule(op::IntegralOperator, g::RTRefSpace, f::RTRefSpace, i, τ, j, hits == 3 && return BogaertSelfPatchStrategy(5) hits == 2 && return BogaertEdgePatchStrategy(8, 4) hits == 1 && return BogaertPointPatchStrategy(2, 3) - rmin = xmin/k - xmin < xtol && return WiltonSERule( - qd.tpoints[1,i], - DoubleQuadRule( - qd.tpoints[2,i], - qd.bpoints[2,j], - ), - ) - return DoubleQuadRule( - qd.tpoints[1,i], - qd.bpoints[1,j], - ) - + return quadrule(op, g, f, i, τ, j, σ, qd.nestedqd, qs.nested_strat) + end \ No newline at end of file diff --git a/src/quadrature/doublenumwiltonsauterqstrat.jl b/src/quadrature/doublenumwiltonsauterqstrat.jl index 6853ce55..f053ea13 100644 --- a/src/quadrature/doublenumwiltonsauterqstrat.jl +++ b/src/quadrature/doublenumwiltonsauterqstrat.jl @@ -1,36 +1,42 @@ -struct DoubleNumWiltonSauterQStrat{R,S} <: AbstractQuadStrat - outer_rule_far::R - inner_rule_far::R +# struct DoubleNumWiltonSauterQStrat{R,S} <: AbstractQuadStrat +# outer_rule_far::R +# inner_rule_far::R +# outer_rule_near::R +# inner_rule_near::R +# sauter_schwab_common_tetr::S +# sauter_schwab_common_face::S +# sauter_schwab_common_edge::S +# sauter_schwab_common_vert::S +# end +struct WiltonQStrat{NestedStrat, R} <: AbstractQuadStrat + nested_strat::NestedStrat outer_rule_near::R inner_rule_near::R - sauter_schwab_common_tetr::S - sauter_schwab_common_face::S - sauter_schwab_common_edge::S - sauter_schwab_common_vert::S end - +DoubleNumWiltonSauterQStrat(a,b,c,d,e,f,g,h) = SauterQStrat(WiltonQStrat(DoubleNumQStrat(a,b),c,d),e,f,g,h) function quaddata(op::IntegralOperator, test_local_space, trial_local_space, - test_charts, trial_charts, qs::DoubleNumWiltonSauterQStrat) + test_charts, trial_charts, qs::WiltonQStrat) - T = coordtype(test_charts[1]) + # T = coordtype(test_charts[1]) # test_local_space = refspace(test_space) # trial_local_space = refspace(trial_space) - tqd = quadpoints(test_local_space, test_charts, (qs.outer_rule_far,qs.outer_rule_near)) - bqd = quadpoints(trial_local_space, trial_charts, (qs.inner_rule_far,qs.inner_rule_near)) + tqd = quadpoints(test_local_space, test_charts, (qs.outer_rule_near,)) + bqd = quadpoints(trial_local_space, trial_charts, (qs.inner_rule_near,)) - leg = ( - convert.(NTuple{2,T},_legendre(qs.sauter_schwab_common_vert,0,1)), - convert.(NTuple{2,T},_legendre(qs.sauter_schwab_common_edge,0,1)), - convert.(NTuple{2,T},_legendre(qs.sauter_schwab_common_face,0,1)),) + # leg = ( + # convert.(NTuple{2,T},_legendre(qs.sauter_schwab_common_vert,0,1)), + # convert.(NTuple{2,T},_legendre(qs.sauter_schwab_common_edge,0,1)), + # convert.(NTuple{2,T},_legendre(qs.sauter_schwab_common_face,0,1)),) - return (tpoints=tqd, bpoints=bqd, gausslegendre=leg) + return (tpoints=tqd, bpoints=bqd, + nestedqd = quaddata(op, test_local_space, trial_local_space, test_charts, trial_charts, qs.nested_strat),) end function quadrule(op::IntegralOperator, g, f, i, τ, j, σ, qd, - qs::DoubleNumWiltonSauterQStrat) + qs::WiltonQStrat) T = eltype(eltype(τ.vertices)) hits = 0 @@ -46,24 +52,22 @@ function quadrule(op::IntegralOperator, g, f, i, τ, j, σ, qd, end end - @assert hits <= 3 + # @assert hits <= 3 - hits == 3 && return SauterSchwabQuadrature.CommonFace(qd.gausslegendre[3]) - hits == 2 && return SauterSchwabQuadrature.CommonEdge(qd.gausslegendre[2]) - hits == 1 && return SauterSchwabQuadrature.CommonVertex(qd.gausslegendre[1]) + # hits == 3 && return SauterSchwabQuadrature.CommonFace(qd.gausslegendre[3]) + # hits == 2 && return SauterSchwabQuadrature.CommonEdge(qd.gausslegendre[2]) + # hits == 1 && return SauterSchwabQuadrature.CommonVertex(qd.gausslegendre[1]) h2 = volume(σ) xtol2 = 0.2 * 0.2 k2 = abs2(gamma(op)) if max(dmin2*k2, dmin2/16h2) < xtol2 return WiltonSERule( - qd.tpoints[2,i], + qd.tpoints[1,i], DoubleQuadRule( - qd.tpoints[2,i], - qd.bpoints[2,j],),) + qd.tpoints[1,i], + qd.bpoints[1,j],),) end - return DoubleQuadRule( - qd.tpoints[1,i], - qd.bpoints[1,j],) + return quadrule(op, g, f, i, τ, j, σ, qd.nestedqd, qs.nested_strat) end \ No newline at end of file diff --git a/src/quadrature/dynamic_quadstrat.jl b/src/quadrature/dynamic_quadstrat.jl new file mode 100644 index 00000000..e69de29b diff --git a/src/quadrature/quadstrats.jl b/src/quadrature/quadstrats.jl index ec9f2e58..da197baa 100644 --- a/src/quadrature/quadstrats.jl +++ b/src/quadrature/quadstrats.jl @@ -5,14 +5,14 @@ using InteractiveUtils -struct SauterSchwab3DQStrat{R,S} <: AbstractQuadStrat - outer_rule::R - inner_rule::R - sauter_schwab_1D::S - sauter_schwab_2D::S - sauter_schwab_3D::S - sauter_schwab_4D::S -end +# struct SauterSchwab3DQStrat{R,S} <: AbstractQuadStrat +# outer_rule::R +# inner_rule::R +# sauter_schwab_1D::S +# sauter_schwab_2D::S +# sauter_schwab_3D::S +# sauter_schwab_4D::S +# end struct OuterNumInnerAnalyticQStrat{R} <: AbstractQuadStrat outer_rule::R @@ -101,4 +101,35 @@ points and weights, geometric quantities, etc. The type of the returned quadrature rule will help in deciding which method of `momintegrals` to dispatch to. """ -function quadrule end \ No newline at end of file +function quadrule end + +function quadcache(biop, test_shapes, trial_shapes, test_elements, trial_elements, qd, quadstrat::NestedQuadStrat) + qdnew = (qd..., cache = nothing, nestedqd = quadcache(biop, test_shapes, trial_shapes, test_elements, trial_elements, qd.nestedqd, quadstrat.nested_strat)) + + + if :nested_strat in fieldnames(typeof(quadstrat)) + return (nestedcache = quadcache(biop, test_shapes, trial_shapes, test_elements, trial_elements, qd.nestedqd, quadstrat.nested_strat),) + end + return () +end + +@generated function quadcache(biop, test_shapes, trial_shapes, test_elements, trial_elements, qd::NamedTuple{Keys}, quadstrat::NestedQuadStrat) where {Keys} + newkeys = [] + newvalsexp = [] + for key in Keys + if key == :nestedqd + newkeys.push(key) + newvalsexp.push(:(quadcache($biop, $test_shapes, $trial_shapes, $test_elements, $trial_elements, $qd.nestedqd, $quadstrat.nested_strat))) + else + newkeys.push(key) + newvalsexp.push(:(qd.$key)) + end + end + ex = quote + return NamedTuple{($(newkeys...,))}((($(newvalsexp...,)),)) + end + return ex +end +function quadcache(biop, test_shapes, trial_shapes, test_elements, trial_elements, qd, quadstrat) + return qd +end diff --git a/src/quadrature/selfsauterdnumotherwiseqstrat.jl b/src/quadrature/selfsauterdnumotherwiseqstrat.jl index 1cbca3a4..ebb9474d 100644 --- a/src/quadrature/selfsauterdnumotherwiseqstrat.jl +++ b/src/quadrature/selfsauterdnumotherwiseqstrat.jl @@ -3,26 +3,29 @@ struct SelfSauterOtherwiseDNumQStrat{R,S} <: AbstractQuadStrat inner_rule::R sauter_schwab_common::S end - +struct SelfSauterQStrat{NestedStrat,S} <: AbstractQuadStrat + nested_strat::NestedStrat + sauter_schwab_common::S +end function quaddata(op::IntegralOperator, test_local_space::RefSpace, trial_local_space::RefSpace, - test_charts, trial_charts, qs::SelfSauterOtherwiseDNumQStrat) + test_charts, trial_charts, qs::SelfSauterQStrat) - T = coordtype(test_charts[1]) + # T = coordtype(test_charts[1]) - tqd = quadpoints(test_local_space, test_charts, (qs.outer_rule,)) - bqd = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,)) + # tqd = quadpoints(test_local_space, test_charts, (qs.outer_rule,)) + # bqd = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,)) leg = ( convert.(NTuple{2,T},_legendre(qs.sauter_schwab_common,0,1)), ) - return (tpoints=tqd, bpoints=bqd, gausslegendre=leg) + return (gausslegendre=leg, nestedqd = quaddata(op, test_local_space, trial_local_space, test_charts, trial_charts, qs.nested_strat)) end function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace, i, τ, j, σ, qd, - qs::SelfSauterOtherwiseDNumQStrat) + qs::SelfSauterQStrat) T = eltype(eltype(τ.vertices)) hits = 0 @@ -42,7 +45,5 @@ function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace, i, τ, j, σ, hits == 3 && return SauterSchwabQuadrature.CommonFace(qd.gausslegendre[1]) - return DoubleQuadRule( - qd.tpoints[1,i], - qd.bpoints[1,j],) + return quadrule(op, g, f, i, τ, j, σ, qd.nestedqd, qs.nested_strat) end \ No newline at end of file diff --git a/src/quadrature/simddoublenum.jl b/src/quadrature/simddoublenum.jl new file mode 100644 index 00000000..df37ff22 --- /dev/null +++ b/src/quadrature/simddoublenum.jl @@ -0,0 +1,469 @@ +using LoopVectorization + +abstract type Operations end +struct Cross <: Operations end +struct Dot <: Operations end +struct Times <: Operations end +struct STimesS <: Operations end +struct VTimesS <: Operations end +struct STimesV <: Operations end +function (::Cross)(a,b) + return cross(a,b) +end +function (::Dot)(a,b) + return transpose(a)*b +end +function (::Union{Times,STimesS,VTimesS,STimesV})(a,b) + return a*b +end + +function generate_cacheSIMDDoubleNum(operator::CompDoubleKern{T,Op1,Op2,Kern}, test_quad_data, trial_quad_data) where {T,Op1,Kern,Op2} + a = length(test_quad_data[1,1].weights) + b = length(trial_quad_data[1,1].weights) + return [zeros(a,b) for _ in 1:cachesizeSIMDDoubleNum(Kern)] +end +function generate_cacheSIMDDoubleNum(operator::MWSingleLayer3D, test_quad_data, trial_quad_data) + a = length(test_quad_data[1,1].weights) + b = length(trial_quad_data[1,1].weights) + return [zeros(a,b) for _ in 1:10] +end +cachesizeSIMDDoubleNum(::Type{<:HH3DGreen}) = 2 +cachesizeSIMDDoubleNum(::Type{<:HH3DGradGreen}) = 6 +function quadpoints_simd(f, els, rules, derivative_needed) + if derivative_needed + return map(Iterators.product(rules, els)) do (rule, el) + pws = quadpoints(el, rule) + (weights = [w for (p,w) in pws], + points = [cartesian(p)[i] for i in 1:length(cartesian(pws[1][1])), (p,w) in pws], + values = [f(p)[j][1][i] for i in 1:length(f(pws[1][1])[1][1]), (p,w) in pws,j in 1:length(f(pws[1][1]))], + derivatives =[f(p)[j][2][i] for i in 1:length(f(pws[1][1])[1][2]), (p,w) in pws,j in 1:length(f(pws[1][1]))]) + end + else + return map(Iterators.product(rules, els)) do (rule, el) + pws = quadpoints(el, rule) + (weights = [w for (p,w) in pws], + points = [cartesian(p)[i] for i in 1:length(cartesian(pws[1][1])), (p,w) in pws], + values = [f(p)[j][1][i] for i in 1:length(f(pws[1][1])[1][1]), (p,w) in pws,j in 1:length(f(pws[1][1]))]) + end + end + +end + +requires_derivative(::MWSingleLayer3D) = true +requires_derivative(::CompDoubleKern) = false + + +struct SIMDDoubleQuadRule{P,Q} + outer_quad_points::P + inner_quad_points::Q + g_cache::Vector{Matrix{Float64}} +end + +""" +momintegrals!(biop, tshs, bshs, tcell, bcell, interactions, strat) + +Function for the computation of moment integrals using simple double quadrature. +""" + + +function momintegrals!(biop::MWSingleLayer3D, + tshs, bshs, tcell, bcell, z, strat::SIMDDoubleQuadRule) + + igd = Integrand(biop, tshs, bshs, tcell, bcell) + + womps = strat.outer_quad_points + wimps = strat.inner_quad_points + + αim = imag(igd.operator.α) + αre = real(igd.operator.α) + βim = imag(igd.operator.β) + βre = real(igd.operator.β) + γim = imag(igd.operator.gamma) + γre = real(igd.operator.gamma) + @turbo for wimpid in eachindex(wimps.weights) + bgeoxi = wimps.points[1,wimpid] + bgeoyi = wimps.points[2,wimpid] + bgeozi = wimps.points[3,wimpid] + for wompid in eachindex(womps.weights) + tgeoxi = womps.points[1,wompid] + tgeoyi = womps.points[2,wompid] + tgeozi = womps.points[3,wompid] + rx = tgeoxi-bgeoxi; ry = tgeoyi-bgeoyi; rz = tgeozi-bgeozi + R = sqrt(rx*rx + ry*ry + rz*rz) + iR = 1 / R + e = exp(-γre*R) + s, c = sincos(γim*R) + scale = e * (i4pi * iR) + greenre = c * scale + greenim = -s * scale + + strat.g_cache[1][wompid, wimpid] = greenre + strat.g_cache[2][wompid, wimpid] = greenim + end + end + for l in 1:3 + for k in 1:3 + zre = 0.0 + zim = 0.0 + @turbo for wimpid in eachindex(wimps.weights) + bvalsxi = wimps.values[1,wimpid,l] + bvalsyi = wimps.values[2,wimpid,l] + bvalszi = wimps.values[3,wimpid,l] + bdiv = wimps.derivatives[wimpid,l] + jy = wimps.weights[wimpid] + + for wompid in eachindex(womps.weights) + tvalsxi = womps.values[1,wompid,k] + tvalsyi = womps.values[2,wompid,k] + tvalszi = womps.values[3,wompid,k] + tdiv = womps.derivatives[wompid,k] + jx = womps.weights[wompid] + j = jx * jy + greenre = strat.g_cache[1][wompid, wimpid] + greenim = strat.g_cache[2][wompid, wimpid] + dotprod = tvalsxi*bvalsxi + tvalsyi*bvalsyi + tvalszi*bvalszi + tdivbdiv = tdiv * bdiv + A_re = αre*dotprod + βre*tdivbdiv + A_im = αim*dotprod + βim*tdivbdiv + + zre += j * (greenre*A_re - greenim*A_im) + zim += j * (greenre*A_im + greenim*A_re) + end + end + + z[k,l] += zre + im*zim + + + end + end + + return z +end + +TransposedStrat(a::SIMDDoubleQuadRule) = a + + +@generated function momintegrals!(biop::CompDoubleKern{T,Op1,Op2,Kern}, + tshs::TSpace, bshs::BSpace, tcell, bcell, z, strat::SIMDDoubleQuadRule) where {T,Op1,Kern,Op2,TSpace,BSpace} + M = numfunctions(tshs, domain(tcell)) + N = numfunctions(bshs, domain(bcell)) + Core.println("M = $M, N = $N") + # load variables + vars = _load_variables(Kern) + + kernexp = _kernexp(Kern) + + assign_green_realexp = assign_green_real(Kern) + assign_green_imagexp = assign_green_imag(Kern) + + + + get_greenrealexp = get_greenreal(Kern) + get_greenimagexp = get_greenimag(Kern) + + + op1 = operationexp(Op1,Kern,bshs) + op2 = operationexp(Op2,tshs,Kern) + trialspacetovarmapexp = trialspacetovarmap(BSpace) + testspacetovarmapexp = testspacetovarmap(TSpace) + loadvaluestestexp = loadvalues_test(TSpace) + loadvaluestrialexp = loadvalues_trial(BSpace) + ex = quote + womps = strat.outer_quad_points + wimps = strat.inner_quad_points + $vars + @turbo for wimpid in eachindex(wimps.weights) + bgeoxi = wimps.points[1,wimpid] + bgeoyi = wimps.points[2,wimpid] + bgeozi = wimps.points[3,wimpid] + for wompid in eachindex(womps.weights) + tgeoxi = womps.points[1,wompid] + tgeoyi = womps.points[2,wompid] + tgeozi = womps.points[3,wompid] + + $kernexp # evaluates the kernel between the given points + $assign_green_realexp # assigns the real part of the kernel to the cache + $assign_green_imagexp # assigns the imaginary part of the kernel to the cache + + end + end + for li in 1:$N + for ki in 1:$M + zre = 0.0 + zim = 0.0 + @turbo for wimpid in eachindex(wimps.weights) + $loadvaluestrialexp # loads the trial space values at the given quadrature point into variables + jy = wimps.weights[wimpid] + + for wompid in eachindex(womps.weights) + $loadvaluestestexp # loads the test space values at the given quadrature point into variables + jx = womps.weights[wompid] + j = jx * jy + + $get_greenrealexp # retrieves the real part of the kernel from the cache + $trialspacetovarmapexp # maps the trial space values to the working variables used in the next operation + $op2 # evaluates the second operation (e.g. multiplying by the trial space values) + $testspacetovarmapexp # maps the test space values to the working variables used in the next operation + $op1 # evaluates the first operation (e.g. multiplying by the test space values) + zre += j * r + + $get_greenimagexp # retrieves the imaginary part of the kernel from the cache + $trialspacetovarmapexp # maps the trial space values to the working variables used in the next operation + $op2 # evaluates the second operation (e.g. multiplying by the trial space values) + $testspacetovarmapexp # maps the test space values to the working variables used in the next operation + $op1 # evaluates the first operation (e.g. multiplying by the test space values) + zim += j * r + + end + end + z[ki,li] += zre + im*zim + end + + end + + return z + end + return ex +end +valuetype(::Type{RTRefSpace{T}}) where {T} = SVector{3,T} +CompScienceMeshes.domain(a::Type{CompScienceMeshes.Simplex{A,B,C,D,T}}) where {A,B,C,D,T} = CompScienceMeshes.ReferenceSimplex{B} + +function _load_variables(::Type{U}) where {U <: Union{HH3DGreen, HH3DGradGreen}} + quote + γim = imag(biop.kernel.gamma) + γre = real(biop.kernel.gamma) + end +end +function loadvalues_test(::Type{U}) where {U} + if valuetype(U) <: Number + return quote + tvals = womps.values[1, wompid,ki] + end + end + if valuetype(U) <: SVector{3} + return quote + tvalsxi = womps.values[1,wompid,ki] + tvalsyi = womps.values[2,wompid,ki] + tvalszi = womps.values[3,wompid,ki] + end + end +end +function loadvalues_trial(::Type{U}) where {U} + if valuetype(U) <: Number + return quote + bvals = wimps.values[1,wimpid,li] + end + end + if valuetype(U) <: SVector{3} + return quote + bvalsxi = wimps.values[1,wimpid,li] + bvalsyi = wimps.values[2,wimpid,li] + bvalszi = wimps.values[3,wimpid,li] + end + end +end +function _kernexp(::Type{<:HH3DGreen}) + quote + rx = tgeoxi-bgeoxi; ry = tgeoyi-bgeoyi; rz = tgeozi-bgeozi + R = sqrt(rx*rx + ry*ry + rz*rz) + iR = 1 / R + e = exp(-γre*R) + s, c = sincos(γim*R) + scale = e * (i4pi * iR) + greenre = c * scale + greenim = -s * scale + end +end +function _kernexp(::Type{<:HH3DGradGreen}) + quote + rx = tgeoxi-bgeoxi; ry = tgeoyi-bgeoyi; rz = tgeozi-bgeozi + R = sqrt(rx*rx + ry*ry + rz*rz) + iR = 1 / R + e = exp(-γre*R) + s, c = sincos(γim*R) + scale = e * (i4pi * iR) + greenre = c*scale + greenim = -s*scale + + gradgreenre = iR*((-γre-iR)*greenre + γim*greenim) + gradgreenim = iR*((-γre-iR)*greenim - γim*greenre) + + + greenrex = gradgreenre * rx + greenrey = gradgreenre * ry + greenrez = gradgreenre * rz + greenimx = gradgreenim * rx + greenimy = gradgreenim * ry + greenimz = gradgreenim * rz + end +end +operationexp(::Type{Cross}, Kern, bshs) = _cross_exp() +operationexp(::Type{Dot}, Kern, bshs) = _dot_exp() +operationexp(::Type{STimesS}, Kern, bshs) = _scalar_times_scalar_exp() +operationexp(::Type{VTimesS}, Kern, bshs) = _vector_times_scalar_exp() +operationexp(::Type{STimesV}, Kern, bshs) = _scalar_times_vector_exp() + +function trialspacetovarmap(a::Type{T}) where {T <: RefSpace} + + if valuetype(T) <: SVector{3} + ex = quote + rx = bvalsxi; ry = bvalsyi; rz = bvalszi + end + elseif valuetype(T) <: Number + ex = quote + r = bvals + end + else + # k = _size(valuetype(T)) + # Core.println(k) + error("unsupported trial space value type") + end + return ex +end +function testspacetovarmap(a::Type{T}) where {T <: RefSpace} + + if valuetype(T) <: SVector{3} + ex = quote + lx = tvalsxi; ly = tvalsyi; lz = tvalszi + end + elseif valuetype(T) <: Number + ex = quote + l = tvals + end + else + error("unsupported trial space value type") + end + return ex +end +""" + valuetype(::Type{<:Space}) + +return the type of the values of the given space evaluated in a meshpoint +""" + +valuetype(::Type{_LocalBasisDot{T}}) where {T} = T +function valuetype(::Type{_LocalBasisTimes{T,U,V}}) where {T,U,V} + + valuetype(U) <: Number && valuetype(V) <: Number && return promote_type(U,V) + return SVector{3,T} +end +valuetype(::Type{_LocalBasisCross{T,U,V}}) where {T,U,V} = SVector{3,T} +valuetype(::Type{NormalVector}) = SVector{3,T} + +# _size(a::Type{<:SVector}) = size(a) +# _size(a::Type{<:Number}) = 1 +function _scalar_times_scalar_exp() + quote + sol = l*r + r = sol + end +end +function _scalar_times_vector_exp() + quote + solx = l*rx + soly = l*ry + solz = l*rz + rx = solx + ry = soly + rz = solz + end +end +function _vector_times_scalar_exp() + quote + solx = lx*r + soly = ly*r + solz = lz*r + rx = solx + ry = soly + rz = solz + end +end +function _dot_exp() + quote + r = rx*lx + ry*ly + rz*lz + end +end +function _cross_exp() + quote + solx = ry*lz - rz*ly + soly = rz*lx - rx*lz + solz = rx*ly - ry*lx + rx = solx + ry = soly + rz = solz + end +end + +function get_greenreal(kern::Type{<:HH3DGreen}) + quote + l = strat.g_cache[1][wompid, wimpid] + end +end +function get_greenimag(kern::Type{<:HH3DGreen}) + quote + l = strat.g_cache[2][wompid, wimpid] + end +end +function get_greenreal(kern::Type{<:HH3DGradGreen}) + quote + lx = strat.g_cache[1][wompid, wimpid] + ly = strat.g_cache[2][wompid, wimpid] + lz = strat.g_cache[3][wompid, wimpid] + end +end +function get_greenimag(kern::Type{<:HH3DGradGreen}) + quote + lx = strat.g_cache[4][wompid, wimpid] + ly = strat.g_cache[5][wompid, wimpid] + lz = strat.g_cache[6][wompid, wimpid] + end +end +function assign_green_real(kern::Type{<:HH3DGreen}) + quote + strat.g_cache[1][wompid, wimpid] = greenre + end +end +function assign_green_imag(kern::Type{<:HH3DGreen}) + quote + strat.g_cache[2][wompid, wimpid] = greenim + end +end +function assign_green_real(kern::Type{<:HH3DGradGreen}) + quote + strat.g_cache[1][wompid, wimpid] = greenrex + strat.g_cache[2][wompid, wimpid] = greenrey + strat.g_cache[3][wompid, wimpid] = greenrez + end +end +function assign_green_imag(kern::Type{<:HH3DGradGreen}) + quote + strat.g_cache[4][wompid, wimpid] = greenimx + strat.g_cache[5][wompid, wimpid] = greenimy + strat.g_cache[6][wompid, wimpid] = greenimz + end +end + +struct RefSIMDDoubleNumRule{R} <: RefQuadStrat + outer_rule::R + inner_rule::R +end + +function quadrule(operator::IntegralOperator, + local_test_basis, local_trial_basis, + test_id, test_element, trial_id, trial_element, qd) + qr = SIMDDoubleQuadRule(qd.tpoints[1, test_id], qd.bpoints[1, trial_id], cache.cache) +end +function quaddata(operator::IntegralOperator, + local_test_basis, local_trial_basis, + test_elements, trial_elements, rule::RefSIMDDoubleNumRule) + derivative_needed = requires_derivative(operator) + (tpoints = quadpoints_simd(local_test_basis, test_elements, (rule.outer_rule,),derivative_needed), + bpoints = quadpoints_simd(local_trial_basis, trial_elements, (rule.inner_rule,),derivative_needed)) +end + +function quadcache(operator::IntegralOperator, + local_test_basis, local_trial_basis, + test_elements, trial_elements, qd, rule::RefSIMDDoubleNumRule) + g_cache = generate_cacheSIMDDoubleNum(operator, qd.tpoints, qd.bpoints) + return (tpoints = qd.tpoints, bpoints = qd.bpoints, cache = g_cache,) +end \ No newline at end of file diff --git a/src/quadrature/strategies/quadstrat.jl b/src/quadrature/strategies/quadstrat.jl index 38c88ca7..3a835093 100644 --- a/src/quadrature/strategies/quadstrat.jl +++ b/src/quadrature/strategies/quadstrat.jl @@ -1,5 +1,5 @@ abstract type AbstractQuadStrat end - +abstract type NestedQuadStrat <: AbstractQuadStrat end function (qs::AbstractQuadStrat)(a, X, Y) qs end diff --git a/src/volumeintegral/vieops.jl b/src/volumeintegral/vieops.jl index 6be6c16d..47ed1ebf 100644 --- a/src/volumeintegral/vieops.jl +++ b/src/volumeintegral/vieops.jl @@ -453,133 +453,133 @@ end -defaultquadstrat(op::VIEOperator, tfs, bfs) = SauterSchwab3DQStrat(3,3,3,3,3,3) +# defaultquadstrat(op::VIEOperator, tfs, bfs) = SauterSchwab3DQStrat(3,3,3,3,3,3) -function quaddata(op::VIEOperator, - test_local_space::RefSpace, trial_local_space::RefSpace, - test_charts, trial_charts, qs::SauterSchwab3DQStrat) +# function quaddata(op::VIEOperator, +# test_local_space::RefSpace, trial_local_space::RefSpace, +# test_charts, trial_charts, qs::SauterSchwab3DQStrat) - #The combinations of rules (6,7) and (5,7 are) BAAAADDDD - # they result in many near singularity evaluations with any - # resemblence of accuracy going down the drain! Simply don't! - # (same for (5,7) btw...). - t_qp = quadpoints(test_local_space, test_charts, (qs.outer_rule,)) - b_qp = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,)) +# #The combinations of rules (6,7) and (5,7 are) BAAAADDDD +# # they result in many near singularity evaluations with any +# # resemblence of accuracy going down the drain! Simply don't! +# # (same for (5,7) btw...). +# t_qp = quadpoints(test_local_space, test_charts, (qs.outer_rule,)) +# b_qp = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,)) - sing_qp = (SauterSchwab3D._legendre(qs.sauter_schwab_1D,0,1), - SauterSchwab3D._shunnham2D(qs.sauter_schwab_2D), - SauterSchwab3D._shunnham3D(qs.sauter_schwab_3D), - SauterSchwab3D._shunnham4D(qs.sauter_schwab_4D),) - - - return (tpoints=t_qp, bpoints=b_qp, sing_qp=sing_qp) -end - - -function _hits(τ, σ) - T = coordtype(τ) - hits = 0 - dtol = 1.0e3 * eps(T) - idx_t = Int64[] - idx_s = Int64[] - sizehint!(idx_t,4) - sizehint!(idx_s,4) - - for (i,t) in enumerate(τ.vertices) - for (j,s) in enumerate(σ.vertices) - d2 = LinearAlgebra.norm_sqr(t-s) - if d2 < dtol - push!(idx_t,i) - push!(idx_s,j) - hits +=1 - break - end - end - end - - return hits, idx_t, idx_s -end - - -""" -tetrahedron-tetrahedron quadrule for the 6D integral ∫∫∫_Ω ∫∫∫_Ω -""" -function quadrule(op::VolumeOperator, g::RefSpace, f::RefSpace, - i, τ::CompScienceMeshes.Simplex{<:Any, 3, <:Any, 4}, - j, σ::CompScienceMeshes.Simplex{<:Any, 3, <:Any, 4}, - qd, qs::SauterSchwab3DQStrat) - - hits, idx_t, idx_s = _hits(τ, σ) - - @assert hits <= 4 - - hits == 4 && return SauterSchwab3D.CommonVolume6D_S(SauterSchwab3D.Singularity6DVolume(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[4])) - hits == 3 && return SauterSchwab3D.CommonFace6D_S(SauterSchwab3D.Singularity6DFace(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3])) - hits == 2 && return SauterSchwab3D.CommonEdge6D_S(SauterSchwab3D.Singularity6DEdge(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3],qd.sing_qp[4])) - hits == 1 && return SauterSchwab3D.CommonVertex6D_S(SauterSchwab3D.Singularity6DPoint(idx_t,idx_s),qd.sing_qp[3]) - #= - #Classic tensor-product quadrature rules - hits == 4 && return SauterSchwab3D.CommonVolume6D(SauterSchwab3D.Singularity6DVolume(idx_t,idx_s),qd.sing_qp[1]) - hits == 3 && return SauterSchwab3D.CommonFace6D(SauterSchwab3D.Singularity6DFace(idx_t,idx_s),qd.sing_qp[1]) - hits == 2 && return SauterSchwab3D.CommonEdge6D(SauterSchwab3D.Singularity6DEdge(idx_t,idx_s),qd.sing_qp[1]) - hits == 1 && return SauterSchwab3D.CommonVertex6D(SauterSchwab3D.Singularity6DPoint(idx_t,idx_s),qd.sing_qp[1]) - =# - - return DoubleQuadRule(qd[1][1,i], qd[2][1,j]) -end - - -""" -triangle-tetrahedron quadrule for the 5D integral ∫∫_Γ ∫∫∫_Ω -""" -function quadrule(op::BoundaryOperator, g::RefSpace, f::RefSpace, - i, τ::CompScienceMeshes.Simplex{<:Any, 2, <:Any, 3}, - j, σ::CompScienceMeshes.Simplex{<:Any, 3, <:Any, 4}, - qd, qs::SauterSchwab3DQStrat) - - hits, idx_t, idx_s = _hits(τ, σ) +# sing_qp = (SauterSchwab3D._legendre(qs.sauter_schwab_1D,0,1), +# SauterSchwab3D._shunnham2D(qs.sauter_schwab_2D), +# SauterSchwab3D._shunnham3D(qs.sauter_schwab_3D), +# SauterSchwab3D._shunnham4D(qs.sauter_schwab_4D),) + + +# return (tpoints=t_qp, bpoints=b_qp, sing_qp=sing_qp) +# end + + +# function _hits(τ, σ) +# T = coordtype(τ) +# hits = 0 +# dtol = 1.0e3 * eps(T) +# idx_t = Int64[] +# idx_s = Int64[] +# sizehint!(idx_t,4) +# sizehint!(idx_s,4) + +# for (i,t) in enumerate(τ.vertices) +# for (j,s) in enumerate(σ.vertices) +# d2 = LinearAlgebra.norm_sqr(t-s) +# if d2 < dtol +# push!(idx_t,i) +# push!(idx_s,j) +# hits +=1 +# break +# end +# end +# end + +# return hits, idx_t, idx_s +# end + + +# """ +# tetrahedron-tetrahedron quadrule for the 6D integral ∫∫∫_Ω ∫∫∫_Ω +# """ +# function quadrule(op::VolumeOperator, g::RefSpace, f::RefSpace, +# i, τ::CompScienceMeshes.Simplex{<:Any, 3, <:Any, 4}, +# j, σ::CompScienceMeshes.Simplex{<:Any, 3, <:Any, 4}, +# qd, qs::SauterSchwab3DQStrat) + +# hits, idx_t, idx_s = _hits(τ, σ) + +# @assert hits <= 4 + +# hits == 4 && return SauterSchwab3D.CommonVolume6D_S(SauterSchwab3D.Singularity6DVolume(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[4])) +# hits == 3 && return SauterSchwab3D.CommonFace6D_S(SauterSchwab3D.Singularity6DFace(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3])) +# hits == 2 && return SauterSchwab3D.CommonEdge6D_S(SauterSchwab3D.Singularity6DEdge(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3],qd.sing_qp[4])) +# hits == 1 && return SauterSchwab3D.CommonVertex6D_S(SauterSchwab3D.Singularity6DPoint(idx_t,idx_s),qd.sing_qp[3]) +# #= +# #Classic tensor-product quadrature rules +# hits == 4 && return SauterSchwab3D.CommonVolume6D(SauterSchwab3D.Singularity6DVolume(idx_t,idx_s),qd.sing_qp[1]) +# hits == 3 && return SauterSchwab3D.CommonFace6D(SauterSchwab3D.Singularity6DFace(idx_t,idx_s),qd.sing_qp[1]) +# hits == 2 && return SauterSchwab3D.CommonEdge6D(SauterSchwab3D.Singularity6DEdge(idx_t,idx_s),qd.sing_qp[1]) +# hits == 1 && return SauterSchwab3D.CommonVertex6D(SauterSchwab3D.Singularity6DPoint(idx_t,idx_s),qd.sing_qp[1]) +# =# + +# return DoubleQuadRule(qd[1][1,i], qd[2][1,j]) +# end + + +# """ +# triangle-tetrahedron quadrule for the 5D integral ∫∫_Γ ∫∫∫_Ω +# """ +# function quadrule(op::BoundaryOperator, g::RefSpace, f::RefSpace, +# i, τ::CompScienceMeshes.Simplex{<:Any, 2, <:Any, 3}, +# j, σ::CompScienceMeshes.Simplex{<:Any, 3, <:Any, 4}, +# qd, qs::SauterSchwab3DQStrat) + +# hits, idx_t, idx_s = _hits(τ, σ) - @assert hits <= 3 - - hits == 3 && return SauterSchwab3D.CommonFace5D_S(SauterSchwab3D.Singularity5DFace(idx_s,idx_t),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3])) - hits == 2 && return SauterSchwab3D.CommonEdge5D_S(SauterSchwab3D.Singularity5DEdge(idx_s,idx_t),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3])) - hits == 1 && return SauterSchwab3D.CommonVertex5D_S(SauterSchwab3D.Singularity5DPoint(idx_s,idx_t),(qd.sing_qp[3],qd.sing_qp[2])) - #= - #Classic tensor-product quadrature rules - hits == 3 && return SauterSchwab3D.CommonFace5D(SauterSchwab3D.Singularity5DFace(idx_s,idx_t),qd.sing_qp[1]) - hits == 2 && return SauterSchwab3D.CommonEdge5D(SauterSchwab3D.Singularity5DEdge(idx_s,idx_t),qd.sing_qp[1]) - hits == 1 && return SauterSchwab3D.CommonVertex5D(SauterSchwab3D.Singularity5DPoint(idx_s,idx_t),qd.sing_qp[1]) - =# - - return DoubleQuadRule(qd[1][1,i], qd[2][1,j]) -end - - -""" -tetrahedron-triangle quadrule for the 5D integral ∫∫∫_Ω ∫∫_Γ -""" -function quadrule(op::BoundaryOperator, g::RefSpace, f::RefSpace, - i, τ::CompScienceMeshes.Simplex{<:Any, 3, <:Any, 4}, - j, σ::CompScienceMeshes.Simplex{<:Any, 2, <:Any, 3}, - qd, qs::SauterSchwab3DQStrat) - - hits, idx_t, idx_s = _hits(τ, σ) - - @assert hits <= 3 - - hits == 3 && return SauterSchwab3D.CommonFace5D_S(SauterSchwab3D.Singularity5DFace(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3])) - hits == 2 && return SauterSchwab3D.CommonEdge5D_S(SauterSchwab3D.Singularity5DEdge(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3])) - hits == 1 && return SauterSchwab3D.CommonVertex5D_S(SauterSchwab3D.Singularity5DPoint(idx_t,idx_s),(qd.sing_qp[3],qd.sing_qp[2])) - #= - #Classic tensor-product quadrature rules - hits == 3 && return SauterSchwab3D.CommonFace5D(SauterSchwab3D.Singularity5DFace(idx_t,idx_s),qd.sing_qp[1]) - hits == 2 && return SauterSchwab3D.CommonEdge5D(SauterSchwab3D.Singularity5DEdge(idx_t,idx_s),qd.sing_qp[1]) - hits == 1 && return SauterSchwab3D.CommonVertex5D(SauterSchwab3D.Singularity5DPoint(idx_t,idx_s),qd.sing_qp[1]) - =# - - return DoubleQuadRule(qd[1][1,i], qd[2][1,j]) -end +# @assert hits <= 3 + +# hits == 3 && return SauterSchwab3D.CommonFace5D_S(SauterSchwab3D.Singularity5DFace(idx_s,idx_t),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3])) +# hits == 2 && return SauterSchwab3D.CommonEdge5D_S(SauterSchwab3D.Singularity5DEdge(idx_s,idx_t),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3])) +# hits == 1 && return SauterSchwab3D.CommonVertex5D_S(SauterSchwab3D.Singularity5DPoint(idx_s,idx_t),(qd.sing_qp[3],qd.sing_qp[2])) +# #= +# #Classic tensor-product quadrature rules +# hits == 3 && return SauterSchwab3D.CommonFace5D(SauterSchwab3D.Singularity5DFace(idx_s,idx_t),qd.sing_qp[1]) +# hits == 2 && return SauterSchwab3D.CommonEdge5D(SauterSchwab3D.Singularity5DEdge(idx_s,idx_t),qd.sing_qp[1]) +# hits == 1 && return SauterSchwab3D.CommonVertex5D(SauterSchwab3D.Singularity5DPoint(idx_s,idx_t),qd.sing_qp[1]) +# =# + +# return DoubleQuadRule(qd[1][1,i], qd[2][1,j]) +# end + + +# """ +# tetrahedron-triangle quadrule for the 5D integral ∫∫∫_Ω ∫∫_Γ +# """ +# function quadrule(op::BoundaryOperator, g::RefSpace, f::RefSpace, +# i, τ::CompScienceMeshes.Simplex{<:Any, 3, <:Any, 4}, +# j, σ::CompScienceMeshes.Simplex{<:Any, 2, <:Any, 3}, +# qd, qs::SauterSchwab3DQStrat) + +# hits, idx_t, idx_s = _hits(τ, σ) + +# @assert hits <= 3 + +# hits == 3 && return SauterSchwab3D.CommonFace5D_S(SauterSchwab3D.Singularity5DFace(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3])) +# hits == 2 && return SauterSchwab3D.CommonEdge5D_S(SauterSchwab3D.Singularity5DEdge(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3])) +# hits == 1 && return SauterSchwab3D.CommonVertex5D_S(SauterSchwab3D.Singularity5DPoint(idx_t,idx_s),(qd.sing_qp[3],qd.sing_qp[2])) +# #= +# #Classic tensor-product quadrature rules +# hits == 3 && return SauterSchwab3D.CommonFace5D(SauterSchwab3D.Singularity5DFace(idx_t,idx_s),qd.sing_qp[1]) +# hits == 2 && return SauterSchwab3D.CommonEdge5D(SauterSchwab3D.Singularity5DEdge(idx_t,idx_s),qd.sing_qp[1]) +# hits == 1 && return SauterSchwab3D.CommonVertex5D(SauterSchwab3D.Singularity5DPoint(idx_t,idx_s),qd.sing_qp[1]) +# =# + +# return DoubleQuadRule(qd[1][1,i], qd[2][1,j]) +# end diff --git a/src/volumeintegral/vsieops.jl b/src/volumeintegral/vsieops.jl index 7338e7d3..817671dc 100644 --- a/src/volumeintegral/vsieops.jl +++ b/src/volumeintegral/vsieops.jl @@ -228,111 +228,111 @@ function integrand(viop::VSIEDoubleLayerT, kerneldata, tvals, tgeo, bvals, bgeo) end =# -defaultquadstrat(op::VSIEOperator, tfs, bfs) = SauterSchwab3DQStrat(3,3,3,3,3,3) +# defaultquadstrat(op::VSIEOperator, tfs, bfs) = SauterSchwab3DQStrat(3,3,3,3,3,3) -function quaddata(op::VSIEOperator, - test_local_space::RefSpace, trial_local_space::RefSpace, - test_charts, trial_charts, qs::SauterSchwab3DQStrat) +# function quaddata(op::VSIEOperator, +# test_local_space::RefSpace, trial_local_space::RefSpace, +# test_charts, trial_charts, qs::SauterSchwab3DQStrat) - #The combinations of rules (6,7) and (5,7 are) BAAAADDDD - # they result in many near singularity evaluations with any - # resemblence of accuracy going down the drain! Simply don't! - # (same for (5,7) btw...). - t_qp = quadpoints(test_local_space, test_charts, (qs.outer_rule)) - b_qp = quadpoints(trial_local_space, trial_charts, (qs.inner_rule)) +# #The combinations of rules (6,7) and (5,7 are) BAAAADDDD +# # they result in many near singularity evaluations with any +# # resemblence of accuracy going down the drain! Simply don't! +# # (same for (5,7) btw...). +# t_qp = quadpoints(test_local_space, test_charts, (qs.outer_rule)) +# b_qp = quadpoints(trial_local_space, trial_charts, (qs.inner_rule)) - sing_qp = (SauterSchwab3D._legendre(qs.sauter_schwab_1D,0,1), - SauterSchwab3D._shunnham2D(qs.sauter_schwab_2D), - SauterSchwab3D._shunnham3D(qs.sauter_schwab_3D), - SauterSchwab3D._shunnham4D(qs.sauter_schwab_4D),) +# sing_qp = (SauterSchwab3D._legendre(qs.sauter_schwab_1D,0,1), +# SauterSchwab3D._shunnham2D(qs.sauter_schwab_2D), +# SauterSchwab3D._shunnham3D(qs.sauter_schwab_3D), +# SauterSchwab3D._shunnham4D(qs.sauter_schwab_4D),) - return (tpoints=t_qp, bpoints=b_qp, sing_qp=sing_qp) -end +# return (tpoints=t_qp, bpoints=b_qp, sing_qp=sing_qp) +# end -quadrule(op::VolumeSurfaceOperator, g::RefSpace, f::RefSpace, i, τ, j, σ, qd, qs) = qr_volume(op, g, f, i, τ, j, σ, qd, qs) +# quadrule(op::VolumeSurfaceOperator, g::RefSpace, f::RefSpace, i, τ, j, σ, qd, qs) = qr_volume(op, g, f, i, τ, j, σ, qd, qs) -function qr_volume(op::VolumeSurfaceOperator, g::RefSpace, f::RefSpace, i, τ, j, σ, qd, - qs::SauterSchwab3DQStrat) +# function qr_volume(op::VolumeSurfaceOperator, g::RefSpace, f::RefSpace, i, τ, j, σ, qd, +# qs::SauterSchwab3DQStrat) - @assert (length(τ.vertices)==3 && length(σ.vertices)==4) "Expected simplex wrong" - - dtol = 1.0e3 * eps(eltype(eltype(τ.vertices))) - - hits = 0 - idx_t = Int64[] - idx_s = Int64[] - sizehint!(idx_t,4) - sizehint!(idx_s,4) - dmin2 = floatmax(eltype(eltype(τ.vertices))) - D = dimension(τ)+dimension(σ) - for (i,t) in enumerate(τ.vertices) - for (j,s) in enumerate(σ.vertices) - d2 = LinearAlgebra.norm_sqr(t-s) - dmin2 = min(dmin2, d2) - if d2 < dtol - push!(idx_t,j) - push!(idx_s,i) - hits +=1 - break - end - end - end - - #singData = SauterSchwab3D.Singularity{D,hits}(idx_t, idx_s ) +# @assert (length(τ.vertices)==3 && length(σ.vertices)==4) "Expected simplex wrong" + +# dtol = 1.0e3 * eps(eltype(eltype(τ.vertices))) + +# hits = 0 +# idx_t = Int64[] +# idx_s = Int64[] +# sizehint!(idx_t,4) +# sizehint!(idx_s,4) +# dmin2 = floatmax(eltype(eltype(τ.vertices))) +# D = dimension(τ)+dimension(σ) +# for (i,t) in enumerate(τ.vertices) +# for (j,s) in enumerate(σ.vertices) +# d2 = LinearAlgebra.norm_sqr(t-s) +# dmin2 = min(dmin2, d2) +# if d2 < dtol +# push!(idx_t,j) +# push!(idx_s,i) +# hits +=1 +# break +# end +# end +# end + +# #singData = SauterSchwab3D.Singularity{D,hits}(idx_t, idx_s ) - hits == 3 && return SauterSchwab3D.CommonFace5D_S(SauterSchwab3D.Singularity5DFace(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3])) - hits == 2 && return SauterSchwab3D.CommonEdge5D_S(SauterSchwab3D.Singularity5DEdge(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3])) - hits == 1 && return SauterSchwab3D.CommonVertex5D_S(SauterSchwab3D.Singularity5DPoint(idx_t,idx_s),(qd.sing_qp[3],qd.sing_qp[2])) - #hits == 0 && return SauterSchwab3D.PositiveDistance5D_S(SauterSchwab3D.Singularity5DPositiveDistance(),(qd.sing_qp[3],qd.sing_qp[2])) +# hits == 3 && return SauterSchwab3D.CommonFace5D_S(SauterSchwab3D.Singularity5DFace(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3])) +# hits == 2 && return SauterSchwab3D.CommonEdge5D_S(SauterSchwab3D.Singularity5DEdge(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3])) +# hits == 1 && return SauterSchwab3D.CommonVertex5D_S(SauterSchwab3D.Singularity5DPoint(idx_t,idx_s),(qd.sing_qp[3],qd.sing_qp[2])) +# #hits == 0 && return SauterSchwab3D.PositiveDistance5D_S(SauterSchwab3D.Singularity5DPositiveDistance(),(qd.sing_qp[3],qd.sing_qp[2])) - return DoubleQuadRule( - qd[1][1,i], - qd[2][1,j]) +# return DoubleQuadRule( +# qd[1][1,i], +# qd[2][1,j]) -end +# end -quadrule(op::BoundarySurfaceOperator, g::RefSpace, f::RefSpace, i, τ, j, σ, qd, qs) = qr_boundary(op, g, f, i, τ, j, σ, qd, qs) +# quadrule(op::BoundarySurfaceOperator, g::RefSpace, f::RefSpace, i, τ, j, σ, qd, qs) = qr_boundary(op, g, f, i, τ, j, σ, qd, qs) -function qr_boundary(op::BoundarySurfaceOperator, g::RefSpace, f::RefSpace, i, τ, j, σ, qd, - qs::SauterSchwab3DQStrat) +# function qr_boundary(op::BoundarySurfaceOperator, g::RefSpace, f::RefSpace, i, τ, j, σ, qd, +# qs::SauterSchwab3DQStrat) - dtol = 1.0e3 * eps(eltype(eltype(τ.vertices))) - - hits = 0 - idx_t = Int64[] - idx_s = Int64[] - sizehint!(idx_t,4) - sizehint!(idx_s,4) - dmin2 = floatmax(eltype(eltype(τ.vertices))) - D = dimension(τ)+dimension(σ) - for (i,t) in enumerate(τ.vertices) - for (j,s) in enumerate(σ.vertices) - d2 = LinearAlgebra.norm_sqr(t-s) - dmin2 = min(dmin2, d2) - if d2 < dtol - push!(idx_t,i) - push!(idx_s,j) - hits +=1 - break - end - end - end - - #singData = SauterSchwab3D.Singularity{D,hits}(idx_t, idx_s ) +# dtol = 1.0e3 * eps(eltype(eltype(τ.vertices))) + +# hits = 0 +# idx_t = Int64[] +# idx_s = Int64[] +# sizehint!(idx_t,4) +# sizehint!(idx_s,4) +# dmin2 = floatmax(eltype(eltype(τ.vertices))) +# D = dimension(τ)+dimension(σ) +# for (i,t) in enumerate(τ.vertices) +# for (j,s) in enumerate(σ.vertices) +# d2 = LinearAlgebra.norm_sqr(t-s) +# dmin2 = min(dmin2, d2) +# if d2 < dtol +# push!(idx_t,i) +# push!(idx_s,j) +# hits +=1 +# break +# end +# end +# end + +# #singData = SauterSchwab3D.Singularity{D,hits}(idx_t, idx_s ) - hits == 3 && return SauterSchwab3D.CommonFace4D_S(SauterSchwab3D.Singularity4DFace(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[3])) - hits == 2 && return SauterSchwab3D.CommonEdge4D_S(SauterSchwab3D.Singularity4DEdge(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2])) - hits == 1 && return SauterSchwab3D.CommonVertex4D_S(SauterSchwab3D.Singularity4DPoint(idx_t,idx_s),(qd.sing_qp[2])) +# hits == 3 && return SauterSchwab3D.CommonFace4D_S(SauterSchwab3D.Singularity4DFace(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[3])) +# hits == 2 && return SauterSchwab3D.CommonEdge4D_S(SauterSchwab3D.Singularity4DEdge(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2])) +# hits == 1 && return SauterSchwab3D.CommonVertex4D_S(SauterSchwab3D.Singularity4DPoint(idx_t,idx_s),(qd.sing_qp[2])) - return DoubleQuadRule( - qd[1][1,i], - qd[2][1,j]) +# return DoubleQuadRule( +# qd[1][1,i], +# qd[2][1,j]) -end +# end diff --git a/test/test_assemble_ndKnd_nonconf.jl b/test/test_assemble_ndKnd_nonconf.jl index d115d3e4..b6bc24ae 100644 --- a/test/test_assemble_ndKnd_nonconf.jl +++ b/test/test_assemble_ndKnd_nonconf.jl @@ -13,7 +13,7 @@ # RT = -n x Nd nxY = n × Y - qs = BEAST.DoubleNumSauterQstrat{Int64, Int64}(2, 3, 4, 4, 4, 4) + qs = BEAST.DoubleNumSauterQstrat(2, 3, 4, 4, 4, 4) A = assemble(K, Y, X; quadstrat=qs) B = assemble(nxK, nxY, X; quadstrat=qs) diff --git a/test/test_nonconf_quadrules.jl b/test/test_nonconf_quadrules.jl index a94db2a6..a1bb7e48 100644 --- a/test/test_nonconf_quadrules.jl +++ b/test/test_nonconf_quadrules.jl @@ -10,7 +10,7 @@ E = Maxwell3D.planewave(;direction=point(0,0,1), polarization=point(1,0,0), wavenumber=0.0) e = n × E - cstrat = BEAST.DoubleNumWiltonSauterQStrat{Int64, Int64}(2, 3, 6, 7, 5, 5, 4, 3) + cstrat = BEAST.DoubleNumWiltonSauterQStrat(2, 3, 6, 7, 5, 5, 4, 3) # cstrat = BEAST.DoubleNumWiltonSauterQStrat{Int64, Int64}(12, 13, 12, 13, 7, 7, 7, 7) nstrat = BEAST.NonConformingIntegralOpQStrat(cstrat) @@ -51,7 +51,7 @@ end m1 = meshcuboid(1.0, 1.0, 1.0, h1; generator=:gmsh) m2 = meshcuboid(1.0, 1.0, 1.0, h2; generator=:gmsh) - cstrat = BEAST.DoubleNumWiltonSauterQStrat{Int64, Int64}(2, 3, 6, 7, 5, 5, 4, 3) + cstrat = BEAST.DoubleNumWiltonSauterQStrat(2, 3, 6, 7, 5, 5, 4, 3) nstrat = BEAST.NonConformingIntegralOpQStrat(cstrat) S = Maxwell3D.singlelayer(alpha=1.0, beta=1.0, gamma=1.0) diff --git a/test/test_vie.jl b/test/test_vie.jl index 9c5dc414..199777ca 100644 --- a/test/test_vie.jl +++ b/test/test_vie.jl @@ -14,7 +14,7 @@ This test does not include the Operators of the EVIE, DVIE and one operator of t @testset "Volume Integral Equations" begin - qs = BEAST.SauterSchwab3DQStrat(3,3,3,3,3,3) + qs = BEAST.DoubleNumSauterQstrat(3,3,3,3,3,3) r = 1.0 h = 0.4 From ddccb8c166cc69a974c7bacbdddfada1773bc2e0 Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Thu, 21 May 2026 16:44:03 +0200 Subject: [PATCH 2/7] update simd --- Project.toml | 6 +- src/BEAST.jl | 1 + src/bases/local/bdm3dlocal.jl | 1 + src/bases/local/bdmlocal.jl | 1 + src/bases/local/gwpdivlocal.jl | 4 ++ src/bases/local/gwplocal.jl | 4 ++ src/bases/local/laglocal.jl | 8 +++ src/bases/local/localcomposedbasis.jl | 2 + src/bases/local/ncrossbdmlocal.jl | 1 + src/bases/local/nd2local.jl | 3 +- src/bases/local/ndlcclocal.jl | 1 + src/bases/local/ndlcdlocal.jl | 1 + src/bases/local/ndlocal.jl | 1 + src/bases/local/rt2local.jl | 1 + src/bases/local/rtlocal.jl | 1 + src/bases/local/rtqlocal.jl | 1 + src/quadrature/doublenumsauterqstrat.jl | 2 +- .../doublenumwiltonbogaertqstrat.jl | 2 +- src/quadrature/doublenumwiltonsauterqstrat.jl | 2 +- src/quadrature/quadstrats.jl | 10 +-- src/quadrature/simddoublenum.jl | 63 ++++++++++++------- 21 files changed, 82 insertions(+), 34 deletions(-) diff --git a/Project.toml b/Project.toml index 67f335a0..bdbe1a75 100644 --- a/Project.toml +++ b/Project.toml @@ -23,6 +23,7 @@ Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" LiftedMaps = "d22a30c1-52ac-4762-a8c9-5838452405e0" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e" +LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" NestedUnitRanges = "032820ab-dc03-4b49-91f4-7d58d4da98b3" OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" @@ -61,10 +62,11 @@ IterativeSolvers = "0.9" Krylov = "0.10.1" LiftedMaps = "0.5.1" LinearMaps = "3.7 - 3.9, 3.11.2" +LoopVectorization = "0.12.170" NestedUnitRanges = "0.2.3" OhMyThreads = "0.8.3" -ProgressMeter = "1.11.0" PlotlyBase = "0.8.21" +ProgressMeter = "1.11.0" Requires = "1" SauterSchwab3D = "0.2" SauterSchwabQuadrature = "2.4.0" @@ -74,5 +76,3 @@ Suppressor = "0.2.8" TestItems = "0.1.1, 1" WiltonInts84 = "0.2.8" julia = "1.10" - - diff --git a/src/BEAST.jl b/src/BEAST.jl index bde48f83..f9b86b60 100644 --- a/src/BEAST.jl +++ b/src/BEAST.jl @@ -340,6 +340,7 @@ include("composedoperators/displacementmesh.jl") include("composedoperators/potentials.jl") include("composedoperators/trace.jl") include("composedoperators/analytic_excitation.jl") +include("quadrature/simddoublenum.jl") const x̂ = point(1, 0, 0) const ŷ = point(0, 1, 0) diff --git a/src/bases/local/bdm3dlocal.jl b/src/bases/local/bdm3dlocal.jl index e8a58d63..69e9b6a7 100644 --- a/src/bases/local/bdm3dlocal.jl +++ b/src/bases/local/bdm3dlocal.jl @@ -1,6 +1,7 @@ struct BDM3DRefSpace{T} <: RefSpace{T} end numfunctions(f::BDM3DRefSpace, ch::CompScienceMeshes.ReferenceSimplex{3}) = 12 +numfunctions(f::Type{<:BDM3DRefSpace}, ch::Type{<:CompScienceMeshes.ReferenceSimplex{3}}) = 12 function (f::BDM3DRefSpace)(p) diff --git a/src/bases/local/bdmlocal.jl b/src/bases/local/bdmlocal.jl index d1a87df1..f2816176 100644 --- a/src/bases/local/bdmlocal.jl +++ b/src/bases/local/bdmlocal.jl @@ -26,6 +26,7 @@ end divergence(ref::BDMRefSpace, sh, el) = [Shape(sh.cellid, 1, sh.coeff/(2*volume(el)))] numfunctions(x::BDMRefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 6 +numfunctions(x::Type{<:BDMRefSpace}, dom::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) = 6 const _vert_perms_bdm = [ (1,2,3), diff --git a/src/bases/local/gwpdivlocal.jl b/src/bases/local/gwpdivlocal.jl index f7d87a82..ac05440c 100644 --- a/src/bases/local/gwpdivlocal.jl +++ b/src/bases/local/gwpdivlocal.jl @@ -6,6 +6,10 @@ function numfunctions(x::GWPDivRefSpace{<:Any,D}, dom::CompScienceMeshes.ReferenceSimplex{2}) where {D} (D+1)*(D+3) end +function numfunctions(x::Type{<:GWPDivRefSpace{<:Any,D}}, + dom::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) where {D} + (D+1)*(D+3) +end function dimtype(x::GWPDivRefSpace{<:Any,D}, dom::CompScienceMeshes.ReferenceSimplex{2}) where {D} Val{(D+1)*(D+3)} diff --git a/src/bases/local/gwplocal.jl b/src/bases/local/gwplocal.jl index 9d9cdd1b..5ec3bd5b 100644 --- a/src/bases/local/gwplocal.jl +++ b/src/bases/local/gwplocal.jl @@ -6,6 +6,10 @@ function numfunctions(x::GWPCurlRefSpace{<:Any,D}, dom::CompScienceMeshes.ReferenceSimplex{2}) where {D} (D+1)*(D+3) end +function numfunctions(x::Type{GWPCurlRefSpace{<:Any,D}}, + dom::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) where {D} + (D+1)*(D+3) +end function dimtype(x::GWPCurlRefSpace{<:Any,D}, dom::CompScienceMeshes.ReferenceSimplex{2}) where {D} Val{(D+1)*(D+3)} diff --git a/src/bases/local/laglocal.jl b/src/bases/local/laglocal.jl index 2b140c88..18fee9de 100644 --- a/src/bases/local/laglocal.jl +++ b/src/bases/local/laglocal.jl @@ -9,6 +9,14 @@ numfunctions(s::LagrangeRefSpace{T,1,3}, ch::CompScienceMeshes.ReferenceSimplex{ numfunctions(s::LagrangeRefSpace{T,2,3}, ch::CompScienceMeshes.ReferenceSimplex{2}) where {T} = 6 numfunctions(s::LagrangeRefSpace{T,Dg}, ch::CompScienceMeshes.ReferenceSimplex{D}) where {T,Dg,D} = binomial(D+Dg,Dg) + + +numfunctions(s::Type{<:LagrangeRefSpace{T,D,2}}, ch::Type{<:CompScienceMeshes.ReferenceSimplex{1}}) where {T,D} = D+1 +numfunctions(s::LagrangeRefSpace{T,0,3}, ch::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) where {T} = 1 +numfunctions(s::LagrangeRefSpace{T,1,3}, ch::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) where {T} = 3 +numfunctions(s::LagrangeRefSpace{T,2,3}, ch::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) where {T} = 6 +numfunctions(s::LagrangeRefSpace{T,Dg}, ch::Type{<:CompScienceMeshes.ReferenceSimplex{D}}) where {T,Dg,D} = binomial(D+Dg,Dg) + # valuetype(ref::LagrangeRefSpace{T}, charttype) where {T} = # SVector{numfunctions(ref), Tuple{T,T}} valuetype(ref::LagrangeRefSpace{T}, charttype) where {T} = T diff --git a/src/bases/local/localcomposedbasis.jl b/src/bases/local/localcomposedbasis.jl index 297ffeb9..d7b16763 100644 --- a/src/bases/local/localcomposedbasis.jl +++ b/src/bases/local/localcomposedbasis.jl @@ -1,6 +1,8 @@ abstract type _LocalBasisOperations{T} <: RefSpace{T} end numfunctions(a::_LocalBasisOperations) = coalesce(numfunctions(a.el1) , numfunctions(a.el2)) numfunctions(a::_LocalBasisOperations,simp) = coalesce(numfunctions(a.el1,simp) , numfunctions(a.el2,simp)) +numfunctions(a::Type{<:_LocalBasisOperations}) = coalesce(numfunctions(a.el1) , numfunctions(a.el2)) +numfunctions(a::Type{<:_LocalBasisOperations},simp) = coalesce(numfunctions(a.el1,simp) , numfunctions(a.el2,simp)) # struct TriangleSupport end # struct TetraSupport end diff --git a/src/bases/local/ncrossbdmlocal.jl b/src/bases/local/ncrossbdmlocal.jl index cb332499..85b22d6c 100644 --- a/src/bases/local/ncrossbdmlocal.jl +++ b/src/bases/local/ncrossbdmlocal.jl @@ -20,3 +20,4 @@ function (f::NCrossBDMRefSpace{T})(p) where T end numfunctions(x::NCrossBDMRefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 6 +numfunctions(x::Type{<:NCrossBDMRefSpace}, dom::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) = 6 diff --git a/src/bases/local/nd2local.jl b/src/bases/local/nd2local.jl index 31dd94c0..8ba6063f 100644 --- a/src/bases/local/nd2local.jl +++ b/src/bases/local/nd2local.jl @@ -27,4 +27,5 @@ function (ϕ::ND2RefSpace)(nbd) end -numfunctions(x::ND2RefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 8 \ No newline at end of file +numfunctions(x::ND2RefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 8 +numfunctions(x::Type{<:ND2RefSpace}, dom::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) = 8 \ No newline at end of file diff --git a/src/bases/local/ndlcclocal.jl b/src/bases/local/ndlcclocal.jl index 37edb79c..7fc830b4 100644 --- a/src/bases/local/ndlcclocal.jl +++ b/src/bases/local/ndlcclocal.jl @@ -59,6 +59,7 @@ function (ϕ::NDLCCRefSpace)(ndlc) end numfunctions(x::NDLCCRefSpace, dom::CompScienceMeshes.ReferenceSimplex{3}) = 6 +numfunctions(x::Type{<:NDLCCRefSpace}, dom::Type{<:CompScienceMeshes.ReferenceSimplex{3}}) = 6 #check orientation function curl(ref::NDLCCRefSpace, sh, el) diff --git a/src/bases/local/ndlcdlocal.jl b/src/bases/local/ndlcdlocal.jl index 8b924607..40bfab9a 100644 --- a/src/bases/local/ndlcdlocal.jl +++ b/src/bases/local/ndlcdlocal.jl @@ -31,6 +31,7 @@ function (ϕ::NDLCDRefSpace)(ndlc) end numfunctions(x::NDLCDRefSpace, dom::CompScienceMeshes.ReferenceSimplex{3}) = 4 +numfunctions(x::Type{<:NDLCDRefSpace}, dom::Type{<:CompScienceMeshes.ReferenceSimplex{3}}) = 4 function ntrace(x::NDLCDRefSpace, el, q, fc) t = zeros(scalartype(x),1,4) diff --git a/src/bases/local/ndlocal.jl b/src/bases/local/ndlocal.jl index 3b3919f0..90303965 100644 --- a/src/bases/local/ndlocal.jl +++ b/src/bases/local/ndlocal.jl @@ -33,6 +33,7 @@ function (ϕ::NDRefSpace)(nbd) end numfunctions(x::NDRefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 3 +numfunctions(x::Type{<:NDRefSpace}, dom::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) = 3 function restrict(ϕ::NDRefSpace{T}, dom1, dom2) where T diff --git a/src/bases/local/rt2local.jl b/src/bases/local/rt2local.jl index bbecede9..66c29767 100644 --- a/src/bases/local/rt2local.jl +++ b/src/bases/local/rt2local.jl @@ -34,6 +34,7 @@ function (f::RT2RefSpace)(p) end numfunctions(x::RT2RefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 8 +numfunctions(x::Type{<:RT2RefSpace}, dom::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) = 8 function interpolate(fields, interpolant::BEAST.RT2RefSpace, chart) diff --git a/src/bases/local/rtlocal.jl b/src/bases/local/rtlocal.jl index a6ecf89a..71720895 100644 --- a/src/bases/local/rtlocal.jl +++ b/src/bases/local/rtlocal.jl @@ -27,6 +27,7 @@ function (ϕ::RTRefSpace)(mp) end numfunctions(x::RTRefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 3 +numfunctions(x::Type{<:RTRefSpace}, dom::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) = 3 divergence(ref::RTRefSpace, sh, el) = [Shape(sh.cellid, 1, sh.coeff/volume(el))] diff --git a/src/bases/local/rtqlocal.jl b/src/bases/local/rtqlocal.jl index cf6e13ea..1a9c354d 100644 --- a/src/bases/local/rtqlocal.jl +++ b/src/bases/local/rtqlocal.jl @@ -23,6 +23,7 @@ function (ϕ::RTQuadRefSpace{T})(p) where {T} end function numfunctions(ϕ::RTQuadRefSpace, dom::CompScienceMeshes.RefQuadrilateral) 4 end +function numfunctions(ϕ::Type{<:RTQuadRefSpace}, dom::Type{<:CompScienceMeshes.RefQuadrilateral}) 4 end function interpolate(fields, interpolant::RTQuadRefSpace{T}, chart) where {T} diff --git a/src/quadrature/doublenumsauterqstrat.jl b/src/quadrature/doublenumsauterqstrat.jl index 9d5255f3..e4196eac 100644 --- a/src/quadrature/doublenumsauterqstrat.jl +++ b/src/quadrature/doublenumsauterqstrat.jl @@ -6,7 +6,7 @@ # sauter_schwab_common_edge::S # sauter_schwab_common_vert::S # end -struct SauterQStrat{NestedStrat,S} <: AbstractQuadStrat +struct SauterQStrat{NestedStrat,S} <: NestedQuadStrat nested_strat::NestedStrat sauter_schwab_common_tetr::S sauter_schwab_common_face::S diff --git a/src/quadrature/doublenumwiltonbogaertqstrat.jl b/src/quadrature/doublenumwiltonbogaertqstrat.jl index a9df3159..ea1387d1 100644 --- a/src/quadrature/doublenumwiltonbogaertqstrat.jl +++ b/src/quadrature/doublenumwiltonbogaertqstrat.jl @@ -5,7 +5,7 @@ # inner_rule_near::R # end DoubleNumWiltonBogaertQStrat(a,b,c,d) = BogaertQStrat(WiltonQStrat(DoubleNumQStrat(a,b),c,d)) -struct BogaertQStrat{NestedStrat} <: AbstractQuadStrat +struct BogaertQStrat{NestedStrat} <: NestedQuadStrat nested_strat::NestedStrat end diff --git a/src/quadrature/doublenumwiltonsauterqstrat.jl b/src/quadrature/doublenumwiltonsauterqstrat.jl index f053ea13..47052fa0 100644 --- a/src/quadrature/doublenumwiltonsauterqstrat.jl +++ b/src/quadrature/doublenumwiltonsauterqstrat.jl @@ -8,7 +8,7 @@ # sauter_schwab_common_edge::S # sauter_schwab_common_vert::S # end -struct WiltonQStrat{NestedStrat, R} <: AbstractQuadStrat +struct WiltonQStrat{NestedStrat, R} <: NestedQuadStrat nested_strat::NestedStrat outer_rule_near::R inner_rule_near::R diff --git a/src/quadrature/quadstrats.jl b/src/quadrature/quadstrats.jl index da197baa..06debf0c 100644 --- a/src/quadrature/quadstrats.jl +++ b/src/quadrature/quadstrats.jl @@ -118,15 +118,15 @@ end newvalsexp = [] for key in Keys if key == :nestedqd - newkeys.push(key) - newvalsexp.push(:(quadcache($biop, $test_shapes, $trial_shapes, $test_elements, $trial_elements, $qd.nestedqd, $quadstrat.nested_strat))) + push!(newkeys,key) + push!(newvalsexp,:(quadcache(biop, test_shapes, trial_shapes, test_elements, trial_elements, qd.nestedqd, quadstrat.nested_strat))) else - newkeys.push(key) - newvalsexp.push(:(qd.$key)) + push!(newkeys,key) + push!(newvalsexp,:(qd.$key)) end end ex = quote - return NamedTuple{($(newkeys...,))}((($(newvalsexp...,)),)) + return NamedTuple{($(newkeys...,))}(($(newvalsexp...),)) end return ex end diff --git a/src/quadrature/simddoublenum.jl b/src/quadrature/simddoublenum.jl index df37ff22..edf21b82 100644 --- a/src/quadrature/simddoublenum.jl +++ b/src/quadrature/simddoublenum.jl @@ -59,6 +59,12 @@ struct SIMDDoubleQuadRule{P,Q} g_cache::Vector{Matrix{Float64}} end + +struct SIMDDoubleNumQStrat{R} <: AbstractQuadStrat + outer_rule::R + inner_rule::R +end + """ momintegrals!(biop, tshs, bshs, tcell, bcell, interactions, strat) @@ -147,7 +153,7 @@ TransposedStrat(a::SIMDDoubleQuadRule) = a tshs::TSpace, bshs::BSpace, tcell, bcell, z, strat::SIMDDoubleQuadRule) where {T,Op1,Kern,Op2,TSpace,BSpace} M = numfunctions(tshs, domain(tcell)) N = numfunctions(bshs, domain(bcell)) - Core.println("M = $M, N = $N") + # Core.println("M = $M, N = $N") # load variables vars = _load_variables(Kern) @@ -162,8 +168,8 @@ TransposedStrat(a::SIMDDoubleQuadRule) = a get_greenimagexp = get_greenimag(Kern) - op1 = operationexp(Op1,Kern,bshs) - op2 = operationexp(Op2,tshs,Kern) + op1 = operationexp(Op1) + op2 = operationexp(Op2) trialspacetovarmapexp = trialspacetovarmap(BSpace) testspacetovarmapexp = testspacetovarmap(TSpace) loadvaluestestexp = loadvalues_test(TSpace) @@ -191,7 +197,7 @@ TransposedStrat(a::SIMDDoubleQuadRule) = a for ki in 1:$M zre = 0.0 zim = 0.0 - @turbo for wimpid in eachindex(wimps.weights) + for wimpid in eachindex(wimps.weights) $loadvaluestrialexp # loads the trial space values at the given quadrature point into variables jy = wimps.weights[wimpid] @@ -223,6 +229,7 @@ TransposedStrat(a::SIMDDoubleQuadRule) = a return z end + # display(ex) return ex end valuetype(::Type{RTRefSpace{T}}) where {T} = SVector{3,T} @@ -297,11 +304,11 @@ function _kernexp(::Type{<:HH3DGradGreen}) greenimz = gradgreenim * rz end end -operationexp(::Type{Cross}, Kern, bshs) = _cross_exp() -operationexp(::Type{Dot}, Kern, bshs) = _dot_exp() -operationexp(::Type{STimesS}, Kern, bshs) = _scalar_times_scalar_exp() -operationexp(::Type{VTimesS}, Kern, bshs) = _vector_times_scalar_exp() -operationexp(::Type{STimesV}, Kern, bshs) = _scalar_times_vector_exp() +operationexp(::Type{Cross}) = _cross_exp() +operationexp(::Type{Dot}) = _dot_exp() +operationexp(::Type{STimesS}) = _scalar_times_scalar_exp() +operationexp(::Type{VTimesS}) = _vector_times_scalar_exp() +operationexp(::Type{STimesV}) = _scalar_times_vector_exp() function trialspacetovarmap(a::Type{T}) where {T <: RefSpace} @@ -385,9 +392,9 @@ function _dot_exp() end function _cross_exp() quote - solx = ry*lz - rz*ly - soly = rz*lx - rx*lz - solz = rx*ly - ry*lx + solx = -ry*lz + rz*ly + soly = -rz*lx + rx*lz + solz = -rx*ly + ry*lx rx = solx ry = soly rz = solz @@ -443,27 +450,39 @@ function assign_green_imag(kern::Type{<:HH3DGradGreen}) end end -struct RefSIMDDoubleNumRule{R} <: RefQuadStrat - outer_rule::R - inner_rule::R -end function quadrule(operator::IntegralOperator, local_test_basis, local_trial_basis, - test_id, test_element, trial_id, trial_element, qd) - qr = SIMDDoubleQuadRule(qd.tpoints[1, test_id], qd.bpoints[1, trial_id], cache.cache) + test_id, test_element, trial_id, trial_element, qd, qs::SIMDDoubleNumQStrat) + qr = SIMDDoubleQuadRule(qd.tpoints[1, test_id], qd.bpoints[1, trial_id], qd.cache) end function quaddata(operator::IntegralOperator, local_test_basis, local_trial_basis, - test_elements, trial_elements, rule::RefSIMDDoubleNumRule) + test_elements, trial_elements, qs::SIMDDoubleNumQStrat) derivative_needed = requires_derivative(operator) - (tpoints = quadpoints_simd(local_test_basis, test_elements, (rule.outer_rule,),derivative_needed), - bpoints = quadpoints_simd(local_trial_basis, trial_elements, (rule.inner_rule,),derivative_needed)) + (tpoints = quadpoints_simd(local_test_basis, test_elements, (qs.outer_rule,),derivative_needed), + bpoints = quadpoints_simd(local_trial_basis, trial_elements, (qs.inner_rule,),derivative_needed)) end function quadcache(operator::IntegralOperator, local_test_basis, local_trial_basis, - test_elements, trial_elements, qd, rule::RefSIMDDoubleNumRule) + test_elements, trial_elements, qd, qs::SIMDDoubleNumQStrat) g_cache = generate_cacheSIMDDoubleNum(operator, qd.tpoints, qd.bpoints) return (tpoints = qd.tpoints, bpoints = qd.bpoints, cache = g_cache,) +end + +@testitem "SIMD DoubleNumSauter Quadrature" begin + using CompScienceMeshes + + Γ = meshcuboid(1.0,1.0,1.0,0.3) + X = raviartthomas(Γ); + + ops = [BEAST.CompDoubleInt(x->x, BEAST.Dot(),BEAST.HH3DGreen(1.0*im),BEAST.STimesV(), x->x), + BEAST.CompDoubleInt(x->x, BEAST.Dot(),BEAST.HH3DGradGreen(1.0*im),BEAST.Cross(), x->x), + Maxwell3D.singlelayer(wavenumber = 1.0)] + for op2 in ops + M3 = assemble(op2,X,X,quadstrat = BEAST.DoubleNumSauterQstrat(6,6,6,6,6,6)); + M4 = assemble(op2,X,X,quadstrat = BEAST.SauterQStrat(BEAST.SIMDDoubleNumQStrat(6,6),6,6,6,6)); + @test isapprox(M3,M4, atol = 1e-5) + end end \ No newline at end of file From e4dc11f0c5571cac109bf04de36285cb9c30cc10 Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Thu, 21 May 2026 16:58:24 +0200 Subject: [PATCH 3/7] added doublelayer simd --- src/quadrature/simddoublenum.jl | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/src/quadrature/simddoublenum.jl b/src/quadrature/simddoublenum.jl index edf21b82..06449826 100644 --- a/src/quadrature/simddoublenum.jl +++ b/src/quadrature/simddoublenum.jl @@ -25,7 +25,12 @@ end function generate_cacheSIMDDoubleNum(operator::MWSingleLayer3D, test_quad_data, trial_quad_data) a = length(test_quad_data[1,1].weights) b = length(trial_quad_data[1,1].weights) - return [zeros(a,b) for _ in 1:10] + return [zeros(a,b) for _ in 1:2] +end +function generate_cacheSIMDDoubleNum(operator::MWDoubleLayer3D, test_quad_data, trial_quad_data) + a = length(test_quad_data[1,1].weights) + b = length(trial_quad_data[1,1].weights) + return [zeros(a,b) for _ in 1:6] end cachesizeSIMDDoubleNum(::Type{<:HH3DGreen}) = 2 cachesizeSIMDDoubleNum(::Type{<:HH3DGradGreen}) = 6 @@ -50,7 +55,7 @@ function quadpoints_simd(f, els, rules, derivative_needed) end requires_derivative(::MWSingleLayer3D) = true -requires_derivative(::CompDoubleKern) = false +requires_derivative(a) = false struct SIMDDoubleQuadRule{P,Q} @@ -70,7 +75,10 @@ momintegrals!(biop, tshs, bshs, tcell, bcell, interactions, strat) Function for the computation of moment integrals using simple double quadrature. """ - +function momintegrals!(biop::MWDoubleLayer3D, + tshs, bshs, tcell, bcell, z, strat::SIMDDoubleQuadRule) + momintegrals!(CompDoubleKern(Dot(),BEAST.HH3DGradGreen(biop.gamma),Cross()), tshs, bshs, tcell, bcell, z, strat::SIMDDoubleQuadRule) +end function momintegrals!(biop::MWSingleLayer3D, tshs, bshs, tcell, bcell, z, strat::SIMDDoubleQuadRule) @@ -479,7 +487,8 @@ end ops = [BEAST.CompDoubleInt(x->x, BEAST.Dot(),BEAST.HH3DGreen(1.0*im),BEAST.STimesV(), x->x), BEAST.CompDoubleInt(x->x, BEAST.Dot(),BEAST.HH3DGradGreen(1.0*im),BEAST.Cross(), x->x), - Maxwell3D.singlelayer(wavenumber = 1.0)] + Maxwell3D.singlelayer(wavenumber = 1.0), + Maxwell3D.doublelayer(wavenumber = 1.0)] for op2 in ops M3 = assemble(op2,X,X,quadstrat = BEAST.DoubleNumSauterQstrat(6,6,6,6,6,6)); M4 = assemble(op2,X,X,quadstrat = BEAST.SauterQStrat(BEAST.SIMDDoubleNumQStrat(6,6),6,6,6,6)); From 4a4446e930b724af214cc6c2c7069d3fefa74835 Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Wed, 27 May 2026 14:08:54 +0200 Subject: [PATCH 4/7] added documentation --- docs/make.jl | 2 ++ docs/src/manual/SIMDquadstrat.md | 22 ++++++++++++++++++++++ docs/src/manual/nestedquadstrat.md | 13 +++++++++++++ docs/src/manual/quadstrat.md | 2 ++ src/quadrature/quadstrats.jl | 25 ++++++++++--------------- src/quadrature/simddoublenum.jl | 3 +++ 6 files changed, 52 insertions(+), 15 deletions(-) create mode 100644 docs/src/manual/SIMDquadstrat.md create mode 100644 docs/src/manual/nestedquadstrat.md diff --git a/docs/make.jl b/docs/make.jl index 093b4867..c0b14b99 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -27,6 +27,8 @@ makedocs(; "Setting the Quadrature Strategy" => "manual/quadstrat.md", "Custom Quadrature Rules" => "manual/quadrule.md", "Custom Operators" => "manual/customop.md", + "Nested Quadrature Strategies" => "manual/nestedquadstrat.md", + "SIMDDoubleNum Quadrature Strategy" => "manual/simddoublenum.md", "Application Examples"=>Any[ "Time-Harmonic"=>Any[ "EFIE"=>"manual/examplesTH/efie.md", diff --git a/docs/src/manual/SIMDquadstrat.md b/docs/src/manual/SIMDquadstrat.md new file mode 100644 index 00000000..7c42397d --- /dev/null +++ b/docs/src/manual/SIMDquadstrat.md @@ -0,0 +1,22 @@ +# SIMDDoubleNum Quadstrat +This quadrature strategy uses the SIMD properties of the hardware to accelerate the double-num interactions. It should be used in combination with a quadrature strategy for the near interacting triangels, for example with the sacuterschwab quadrature. +The SIMD quadtarture is implemented for the Maxwell3D.singlelayer operator assembled in a raviart-thomas basis on triangles and for the composed operator structures at any basis embedded in a 3D space. +The Maxwell3D.doublelayer operator is assembled by casting the structure to the correpsonding compsoed operator. +Example +```julia + op = Maxwell3D.singlelayer(wavenumber=1.0) + assemble(op,X,X,quadstrat = BEAST.SIMDDoubleNumSauterQstrat(6,6,6,6,6,6)) + op = Maxwell3D.doublelayer(wavenumber=1.0) + assemble(op,X,X,quadstrat = BEAST.SIMDDoubleNumSauterQstrat(6,6,6,6,6,6)) + op = Maxwell3D.doublelayer(wavenumber=1.0) + assemble(op,X,X,quadstrat = BEAST.SIMDDoubleNumSauterQstrat(6,6,6,6,6,6)) + op = BEAST.CompDoubleInt(x->x, BEAST.Dot(),BEAST.HH3DGreen(1.0*im),BEAST.STimesV(), x->x) + assemble(op,X,X,quadstrat = BEAST.SIMDDoubleNumSauterQstrat(6,6,6,6,6,6)) +``` + +Note: The operators used in the operator fields of the CompDoubleInt should be of the specific types shown below +vector scalar multiplication: BEAST.VTimesS +scalar vector multiplication: BEAST.STimesV +scalar scalar multiplication: BEAST.STimesS +dot product: BEAST.Dot +cross product: BEAST.Cross diff --git a/docs/src/manual/nestedquadstrat.md b/docs/src/manual/nestedquadstrat.md new file mode 100644 index 00000000..a2a66327 --- /dev/null +++ b/docs/src/manual/nestedquadstrat.md @@ -0,0 +1,13 @@ +# Nested Quadrature Strategies + +Some quadrature strategies are only useble for a subest of all simplex-simplex interactions. If two simplicces do fullfill the right conditions this quadrature strategie returns a specified quadrature rule. If not, it calls the nested quadrature strategy on the simplex-simplex pair. This nested quadrature strategy is stored in the nested_strat field that a nested quadrature strategy must have. +An example of a nested quadstrat is given by +```julia +struct SauterQStrat{NestedStrat,S} <: NestedQuadStrat + nested_strat::NestedStrat + sauter_schwab_common_tetr::S + sauter_schwab_common_face::S + sauter_schwab_common_edge::S + sauter_schwab_common_vert::S +end +``` \ No newline at end of file diff --git a/docs/src/manual/quadstrat.md b/docs/src/manual/quadstrat.md index 1b7322c4..233d83af 100644 --- a/docs/src/manual/quadstrat.md +++ b/docs/src/manual/quadstrat.md @@ -118,4 +118,6 @@ Z3 = assemble(𝑇, RT, RT; quadstrat=myquadstrat2) Best practice is too return `BEAST.defaultquadstrat(op, testnfs, trialfns)` by default to ensure that all operators are supported, also those for which no explicit overwrite is specified. +# Quadrature cache +The `quadcache` function is called for each thread seperately and ats writable working memory to the quadrature data. diff --git a/src/quadrature/quadstrats.jl b/src/quadrature/quadstrats.jl index 06debf0c..56f9b58a 100644 --- a/src/quadrature/quadstrats.jl +++ b/src/quadrature/quadstrats.jl @@ -103,33 +103,28 @@ The type of the returned quadrature rule will help in deciding which method of """ function quadrule end -function quadcache(biop, test_shapes, trial_shapes, test_elements, trial_elements, qd, quadstrat::NestedQuadStrat) - qdnew = (qd..., cache = nothing, nestedqd = quadcache(biop, test_shapes, trial_shapes, test_elements, trial_elements, qd.nestedqd, quadstrat.nested_strat)) - - - if :nested_strat in fieldnames(typeof(quadstrat)) - return (nestedcache = quadcache(biop, test_shapes, trial_shapes, test_elements, trial_elements, qd.nestedqd, quadstrat.nested_strat),) - end - return () +""" + quadcache(operator, test_refspace, trial_refspace, test_elements, trial_elements, quad_data, quadstrat) + +This function is called for each thread separately and adds temporary writable memory +to the quaddata named tuple if necessary. If not, it leaves the quaddata as is. +""" +function quadcache(biop, test_shapes, trial_shapes, test_elements, trial_elements, qd, quadstrat) + return qd end @generated function quadcache(biop, test_shapes, trial_shapes, test_elements, trial_elements, qd::NamedTuple{Keys}, quadstrat::NestedQuadStrat) where {Keys} - newkeys = [] newvalsexp = [] for key in Keys if key == :nestedqd - push!(newkeys,key) push!(newvalsexp,:(quadcache(biop, test_shapes, trial_shapes, test_elements, trial_elements, qd.nestedqd, quadstrat.nested_strat))) else - push!(newkeys,key) push!(newvalsexp,:(qd.$key)) end end ex = quote - return NamedTuple{($(newkeys...,))}(($(newvalsexp...),)) + return NamedTuple{($(Keys...,))}(($(newvalsexp...),)) end return ex end -function quadcache(biop, test_shapes, trial_shapes, test_elements, trial_elements, qd, quadstrat) - return qd -end + diff --git a/src/quadrature/simddoublenum.jl b/src/quadrature/simddoublenum.jl index 06449826..5b985b11 100644 --- a/src/quadrature/simddoublenum.jl +++ b/src/quadrature/simddoublenum.jl @@ -479,6 +479,9 @@ function quadcache(operator::IntegralOperator, return (tpoints = qd.tpoints, bpoints = qd.bpoints, cache = g_cache,) end +# frequently used quadstrats +SIMDDoubleNumSauterQStrat(a,b,c,d,e,f) = SauterQStrat(SIMDDoubleNumQStrat(a,b),c,d,e,f) + @testitem "SIMD DoubleNumSauter Quadrature" begin using CompScienceMeshes From ce4943580ab301efd8b2cb77f8c93862a7332fd4 Mon Sep 17 00:00:00 2001 From: Kristof Cools Date: Tue, 26 May 2026 12:13:56 +0200 Subject: [PATCH 5/7] assemble on bilforms reuses blocks --- CHANGELOG.md | 6 +++++- examples/pmchwt.jl | 23 ++++++++++++++++++----- src/solvers/solver.jl | 23 ++++++++++++++++------- test/test_solvers_solver.jl | 27 +++++++++++++++++++++++++++ 4 files changed, 66 insertions(+), 13 deletions(-) create mode 100644 test/test_solvers_solver.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index 6cac0180..bd4b85f4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,8 +1,12 @@ -# Changelog +# New in 2.10.0 - `solve` and `solve!` on `GMRESSolver` allow for the defaults preconditioners to be overwritten - `materialize` keyword to assemble for bilinear forms allows for fine-grained control over the assembly of operators (e.g. dense representations vs H-matrix representations) - Restart for GMRES defaults to the maximum number of iterations +- `assemble` of bilinear forms checks for recurring blocks and reuses them. + +# New in 2.9.0 + - Cell coloring based lock-free multi-threaded assembly for frequency domain integral operators - Higher order Lagrange elements (cx and c0) on segments - All LinearMaps can be cast into bilinear forms. Assembly of the latter trivially returns the underlying LinearMap. This enables for example the construction inverses that explicitly exploit block diagonal structure. diff --git a/examples/pmchwt.jl b/examples/pmchwt.jl index 49f32942..3b281ab4 100644 --- a/examples/pmchwt.jl +++ b/examples/pmchwt.jl @@ -59,12 +59,25 @@ h = (n × H) × n @hilbertspace k l α, α′ = 1/η, 1/η′ -pmchwt = @discretise( - (η*T+η′*T′)[k,j] - (K+K′)[k,m] + - (K+K′)[l,j] + (α*T+α′*T′)[l,m] == -e[k] - h[l], - j∈X, m∈X, k∈X, l∈X) +# pmchwt = @discretise( +# (η*T+η′*T′)[k,j] - (K+K′)[k,m] + +# (K+K′)[l,j] + (α*T+α′*T′)[l,m] == -e[k] - h[l], +# j∈X, m∈X, k∈X, l∈X) -u = solve(pmchwt) +a = ( + η*T[k,j] + η′*T′[k,j] - K[k,m] - K′[k,m] + + K[l,j] + K′[l,j] + α*T[l,m] + α′*T′[l,m]) +l = -e[k] - h[l] + +𝕏 = X × X +A = assemble(a, 𝕏, 𝕏; threading=:cellcoloring) +b = assemble(l, 𝕏) + +A⁻¹ = BEAST.GMRESSolver(A, reltol=1e-5, maxiter=1000) +u = A⁻¹ * b + +# error() +# u = solve(pmchwt) #preconditioner #= diff --git a/src/solvers/solver.jl b/src/solvers/solver.jl index d59dd10a..873e5e9d 100644 --- a/src/solvers/solver.jl +++ b/src/solvers/solver.jl @@ -223,7 +223,8 @@ lift(a::ConvolutionOperators.AbstractConvOp ,I,J,U,V) = ConvolutionOperators.LiftedConvOp(a, U, V, I, J) -function assemble(bf::BilForm, X::DirectProductSpace, Y::DirectProductSpace; +function assemble(bf::BilForm, X::DirectProductSpace, Y::DirectProductSpace, + archive=Dict(); materialize=BEAST.assemble, kwargs...) T = Int32 @@ -272,18 +273,26 @@ function assemble(bf::BilForm, X::DirectProductSpace, Y::DirectProductSpace; y = op[end](op[1:end-1]..., y) end + term.coeff == 0 && continue a = term.kernel - # @show typeof(a) - # @show typeof(x) - # @show typeof(y) - z = (a isa BilForm) ? - assemble(a, x, y; materialize, kwargs...) : - materialize(a, x, y; kwargs...) + z = if a isa BilForm + assemble(a, x, y, archive; materialize, kwargs...) + else + get!(archive, (a,x,y)) do + materialize(a, x, y; kwargs...) + end + # materialize(a, x, y; kwargs...) + end + # z = (a isa BilForm) ? + # assemble(a, x, y; materialize, kwargs...) : + # materialize(a, x, y; kwargs...) Smap = term.coeff * lift(z, Block(term.test_id), Block(term.trial_id), U, V) T = promote_type(T, eltype(Smap)) push!(lincombv, Smap) end + @show length(archive) + if spaceTimeBasis return sum(lincombv) else diff --git a/test/test_solvers_solver.jl b/test/test_solvers_solver.jl new file mode 100644 index 00000000..13c774f1 --- /dev/null +++ b/test/test_solvers_solver.jl @@ -0,0 +1,27 @@ +@testitem "assemble BilForm: archive" begin + + using CompScienceMeshes + using LinearAlgebra + + fn = joinpath(pkgdir(BEAST), "test", "assets", "sphere45.in") + m = CompScienceMeshes.readmesh(fn) + + X = raviartthomas(m) + 𝕏 = X × X + + I = BEAST.Identity() + T = Maxwell3D.singlelayer(wavenumber=1.0) + + @hilbertspace m j + @hilbertspace k l + + a = ( + I[k,m] + 2*T[k,j] + + 3im*T[l,m] - I[l,j] + ) + + A = assemble(a, 𝕏, 𝕏; threading=:cellcoloring) + # this test is brittle but not sure how to do this otherwise... + @test A.maps[1].lmap.A.lmap === A.maps[4].lmap.A.lmap + @test A.maps[2].lmap.A.lmap === A.maps[3].lmap.A.lmap +end \ No newline at end of file From c14b6f166c1ac27824cb34429525e1cfe64a6a21 Mon Sep 17 00:00:00 2001 From: Kristof Cools Date: Tue, 26 May 2026 12:40:34 +0200 Subject: [PATCH 6/7] Use pinv to invert hypersingular HH2D at wavenumber zero --- test/test_hh2d_nearfield.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_hh2d_nearfield.jl b/test/test_hh2d_nearfield.jl index 0f81382e..2ab62c75 100644 --- a/test/test_hh2d_nearfield.jl +++ b/test/test_hh2d_nearfield.jl @@ -267,7 +267,7 @@ end ρ_EDPSL = M_EDPSL \ (-gD0) ρ_EDPDL = M_EDPDL \ (-gD1) - ρ_ENPDL = M_ENPDL \ gN + ρ_ENPDL = pinv(M_ENPDL) * gN ρ_ENPSL = M_ENPSL \ (-gN) testcircle = meshcircle(1.2 * r, 1.2 * 0.6 * r) From c0694692c29999bff1d9e4cadfc35fd5193872fa Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Wed, 27 May 2026 14:23:58 +0200 Subject: [PATCH 7/7] update --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index c0b14b99..3e582105 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -28,7 +28,7 @@ makedocs(; "Custom Quadrature Rules" => "manual/quadrule.md", "Custom Operators" => "manual/customop.md", "Nested Quadrature Strategies" => "manual/nestedquadstrat.md", - "SIMDDoubleNum Quadrature Strategy" => "manual/simddoublenum.md", + "SIMDDoubleNum Quadrature Strategy" => "manual/SIMDquadstrat.md", "Application Examples"=>Any[ "Time-Harmonic"=>Any[ "EFIE"=>"manual/examplesTH/efie.md",