From b1dc8dc9f2d148af509f5d499244027cca9f9f46 Mon Sep 17 00:00:00 2001 From: "Simon B. Adrian" Date: Sat, 8 Nov 2025 16:40:57 +0100 Subject: [PATCH 1/2] Add 1D Higher Order Lagrange Functions (Cx and C0) Polynomials are implemented hard-coded up to order 10 for maximal speed. TODO: The Mie unit test should be adapted once curvilinear meshes are available. So far, we only confirm that higher-order functions do not break the performance. However, since the mesh is not refined the geometrical error starts to dominate quickly. --- src/BEAST.jl | 1 + src/bases/lagrange.jl | 171 ++++++++++++++++++- src/bases/local/laglocal.jl | 183 ++++++++++++++++++++ src/interpolation.jl | 19 +++ test/runtests.jl | 2 + test/test_hh2d_mie_higher_order.jl | 262 +++++++++++++++++++++++++++++ test/test_lagrange.jl | 84 +++++++++ 7 files changed, 719 insertions(+), 3 deletions(-) create mode 100644 test/test_hh2d_mie_higher_order.jl create mode 100644 test/test_lagrange.jl diff --git a/src/BEAST.jl b/src/BEAST.jl index 88f38061..4a517e81 100644 --- a/src/BEAST.jl +++ b/src/BEAST.jl @@ -34,6 +34,7 @@ export planewave export RefSpace, numfunctions, coordtype, scalartype, assemblydata, geometry, refspace, valuetype export lagrangecxd0, lagrangec0d1, duallagrangec0d1, lagrangec0d2, unitfunctioncxd0, unitfunctionc0d1 export duallagrangecxd0 +export lagrangec0, lagrangecx export lagdimension export restrict export raviartthomas, raowiltonglisson, positions diff --git a/src/bases/lagrange.jl b/src/bases/lagrange.jl index b044dce0..cb39257b 100644 --- a/src/bases/lagrange.jl +++ b/src/bases/lagrange.jl @@ -177,7 +177,10 @@ end """ lagrangec0d1(mesh[, bnd]) -Construct the basis of continuous, piecewise linear basis functions subordinate to mesh `mesh`. Basis functions are constructed at vertices in the interionr of the mesh and on the closure of 'bnd'. In particular, leaving out the second argument creates a finite element space subject to homogeneous Dirichlet boundary conditions. +Construct the basis of continuous, piecewise linear basis functions subordinate +to mesh `mesh`. Basis functions are constructed at vertices in the interior of +the mesh and on the closure of 'bnd'. In particular, leaving out the second argument +creates a finite element space subject to homogeneous Dirichlet boundary conditions. """ function lagrangec0d1_dirichlet(mesh) @@ -198,6 +201,24 @@ function lagrangec0d1_dirichlet(mesh) lagrangec0d1(mesh, vertexlist, Val{dimension(mesh)+1}) end +function lagrangec0_dirichlet(mesh; order=1) + + verts = skeleton(mesh, 0) + detached = trues(numvertices(mesh)) + for v in cells(verts) + detached[v] = false + end + + bnd = boundary(mesh) + bndverts = skeleton(bnd, 0) + notonbnd = trues(numvertices(mesh)) + for v in cells(bndverts) + notonbnd[v] = false + end + + vertexlist = findall(notonbnd .& .!detached) + lagrangec0(mesh, vertexlist, Val{dimension(mesh)+1}; order=order) +end function interior_and_junction_vertices(mesh, jct) verts = skeleton(mesh, 0) @@ -310,7 +331,8 @@ end """ lagrangec0d1(mesh; dirichlet=[true|false]) -> basis -Build lagrangec0d1 elements, including (dirichlet=false) or excluding (dirichlet=true) those attached to boundary vertices. +Build lagrangec0d1 elements, including (dirichlet=false) +or excluding (dirichlet=true) those attached to boundary vertices. """ function lagrangec0d1(mesh; dirichlet::Bool=true) if dirichlet == false @@ -321,6 +343,17 @@ function lagrangec0d1(mesh; dirichlet::Bool=true) end end +function lagrangec0(mesh; order=1, dirichlet::Bool=true) + + @assert order>0 + if dirichlet == false + # return lagrangec0d1(mesh, boundary(mesh)) + return lagrangec0(mesh, skeleton(mesh,0); order=order) + else + return lagrangec0_dirichlet(mesh; order=order) + end +end + function lagrangec0d1(mesh, jct) vertexlist = interior_and_junction_vertices(mesh, jct) lagrangec0d1(mesh, vertexlist, Val{dimension(mesh)+1}) @@ -413,8 +446,65 @@ function lagrangec0d1(mesh, vertexlist, ::Type{Val{2}}) LagrangeBasis{1,0,NF}(geometry, fns, pos) end +function lagrangec0(mesh, vertexlist, ::Type{Val{2}}; order=1) + + if order > 10 + error("Maximal 1D Lagrange polynomial order supported is 10") + end + + T = coordtype(mesh) + U = universedimension(mesh) + P = vertextype(mesh) + S = Shape{T} + + geometry = mesh + + cellids, ncells = vertextocellmap(mesh) + + # create the local shapes + numverts = numvertices(mesh) + + fns = Vector{Vector{S}}() + pos = Vector{P}() + + sizehint!(fns, length(vertexlist)) + sizehint!(pos, length(vertexlist)) + for v in vertexlist + + numshapes = ncells[v] + numshapes == 0 && continue # skip detached vertices + + shapes = Vector{S}(undef,numshapes) + for s in 1: numshapes + c = cellids[v,s] + cell = mesh.faces[c] + if cell[1] == v + shapes[s] = Shape(c, 1, T(1.0)) + elseif cell[2] == v + shapes[s] = Shape(c, 2, T(1.0)) + else + error("Junctions not supported") + end + end + + push!(fns, shapes) + push!(pos, mesh.vertices[v]) + end -function lagrangec0d1(mesh, nodes::CompScienceMeshes.AbstractMesh{U,1} where {U}) + NF = order + 1 + + for (c,cell) in enumerate(mesh) + ch = chart(mesh,cell) + for r in 3:NF + push!(fns, S[S(c,r,T(1.0))]) + push!(pos, cartesian(neighborhood(ch, (r-2)/(NF-1)))) + end + end + + LagrangeBasis{order,0,NF}(geometry, fns, pos) +end + +function lagrangec0d1(mesh, nodes::CompScienceMeshes.AbstractMesh{<:Any,1}) Conn = connectivity(nodes, mesh, abs) rows = rowvals(Conn) @@ -438,9 +528,52 @@ function lagrangec0d1(mesh, nodes::CompScienceMeshes.AbstractMesh{U,1} where {U} end NF = dimension(mesh) + 1 + LagrangeBasis{1,0,NF}(mesh, fns, pos) end +function lagrangec0(mesh, nodes::CompScienceMeshes.AbstractMesh{<:Any,1}; order=1) + + if order > 10 + error("Maximal 1D Lagrange polynomial order supported is 10") + end + + Conn = connectivity(nodes, mesh, abs) + rows = rowvals(Conn) + vals = nonzeros(Conn) + + T = coordtype(mesh) + P = vertextype(mesh) + S = Shape{T} + + fns = Vector{Vector{S}}() + pos = Vector{P}() + + u = T(1.0) + + for (i,node) in enumerate(nodes) + fn = Vector{S}() + for k in nzrange(Conn,i) + cellid = rows[k] + refid = vals[k] + push!(fn, Shape(cellid, refid, u)) + end + push!(fns,fn) + push!(pos,cartesian(center(chart(nodes,node)))) + end + + NF = order + 1 + + for (c, cell) in enumerate(mesh) + ch = chart(mesh, cell) + for r in 3:NF + push!(fns, S[S(c, r, u)]) + push!(pos, cartesian(neighborhood(ch, SVector(u - (r-2)/(NF-1))))) + end + end + + LagrangeBasis{order,0,NF}(mesh, fns, pos) +end function lagrangec0d2(mesh::CompScienceMeshes.AbstractMesh{U,3}, nodes::CompScienceMeshes.AbstractMesh{U,1}, @@ -1027,7 +1160,39 @@ function dual0forms_body(mesh::CompScienceMeshes.AbstractMesh{<:Any,3}, refd, bn LagrangeBasis{1,0,3}(refd, bfs, pos) end +function lagrangecx(mesh::CompScienceMeshes.AbstractMesh{<:Any,2}; order) + + T = coordtype(mesh) + + P = vertextype(mesh) + S = Shape{T} + + fns = Vector{Vector{S}}() + pos = Vector{P}() + u = one(T) + for (c,cell) in enumerate(mesh) + ch = chart(mesh,cell) + + push!(fns, S[S(c,1,u)]) + push!(pos, cartesian(neighborhood(ch, T(1.0)))) + + push!(fns, S[S(c,2,u)]) + push!(pos, cartesian(neighborhood(ch, T(0.0)))) + end + + NF = order + 1 + + for (c,cell) in enumerate(mesh) + ch = chart(mesh,cell) + for r in 3:NF + push!(fns, S[S(c,r,u)]) + push!(pos, cartesian(neighborhood(ch, u - (r-2)/(NF-1)))) + end + end + + return LagrangeBasis{order,-1,NF}(mesh, fns, pos) +end function lagrangecx(mesh::CompScienceMeshes.AbstractMesh{<:Any,3}; order) diff --git a/src/bases/local/laglocal.jl b/src/bases/local/laglocal.jl index e590f338..2b140c88 100644 --- a/src/bases/local/laglocal.jl +++ b/src/bases/local/laglocal.jl @@ -19,6 +19,7 @@ valuetype(ref::LagrangeRefSpace{T}, charttype) where {T} = T # Evaluate linear Lagrange elements on a segment +# The derivative denotes the tangential derivative function (f::LagrangeRefSpace{T,1,2})(mp) where T u = mp.bary[1] j = jacobian(mp) @@ -27,6 +28,188 @@ function (f::LagrangeRefSpace{T,1,2})(mp) where T (value=1-u, derivative= 1/j)) end +function (f::LagrangeRefSpace{T,2,2})(mp) where T + u = mp.bary[1] + j = jacobian(mp) + SVector( + (value=-u*(1 - 2*u), derivative=-(4*u - 1)/j), + (value=(1 - 2*u)*(1 - u), derivative=-(4*u - 3)/j), + (value=4*u*(1 - u), derivative=-(4 - 8*u)/j) + ) +end + +function (f::LagrangeRefSpace{T,3,2})(mp) where T + u = mp.bary[1] + j = jacobian(mp) + SVector( + (value=u*(1 - 3*u)*(2 - 3*u)/2, derivative=-(27*u^2/2 - 9*u + 1)/j), + (value=(1 - 3*u)*(1 - u)*(2 - 3*u)/2, derivative=(27*u^2/2 - 18*u + 11/2)/j), + (value=-9*u*(1 - 3*u)*(1 - u)/2, derivative=(81*u^2/2 - 36*u + 9/2)/j), + (value=9*u*(1 - u)*(2 - 3*u)/2, derivative=-(81*u^2/2 - 45*u + 9)/j) + ) +end + +function (f::LagrangeRefSpace{T,4,2})(mp) where T + u = mp.bary[1] + j = jacobian(mp) + SVector( + (value=-u*(1 - 4*u)*(1 - 2*u)*(3 - 4*u)/3, + derivative=-(128*u^3/3 - 48*u^2 + 44*u/3 - 1)/j), + (value=(1 - 4*u)*(1 - 2*u)*(1 - u)*(3 - 4*u)/3, + derivative=-(128*u^3/3 - 80*u^2 + 140*u/3 - 25/3)/j), + (value=16*u*(1 - 4*u)*(1 - 2*u)*(1 - u)/3, + derivative=-( -512*u^3/3 + 224*u^2 - 224*u/3 + 16/3)/j), + (value=-4*u*(1 - 4*u)*(1 - u)*(3 - 4*u), + derivative=-(256*u^3 - 384*u^2 + 152*u - 12)/j), + (value=16*u*(1 - 2*u)*(1 - u)*(3 - 4*u)/3, + derivative=-( -512*u^3/3 + 288*u^2 - 416*u/3 + 16)/j) + ) +end + +function (f::LagrangeRefSpace{T,5,2})(mp) where T + u = mp.bary[1] + j = jacobian(mp) + SVector( + (value=u*(1 - 5*u)*(2 - 5*u)*(3 - 5*u)*(4 - 5*u)/24, + derivative=-(3125*u^4/24 - 625*u^3/3 + 875*u^2/8 - 125*u/6 + 1)/j), + (value=(1 - 5*u)*(1 - u)*(2 - 5*u)*(3 - 5*u)*(4 - 5*u)/24, + derivative=(3125*u^4/24 - 625*u^3/2 + 2125*u^2/8 - 375*u/4 + 137/12)/j), + (value=-25*u*(1 - 5*u)*(1 - u)*(2 - 5*u)*(3 - 5*u)/24, + derivative=(15625*u^4/24 - 6875*u^3/6 + 5125*u^2/8 - 1525*u/12 + 25/4)/j), + (value=25*u*(1 - 5*u)*(1 - u)*(2 - 5*u)*(4 - 5*u)/12, + derivative=-(15625*u^4/12 - 2500*u^3 + 6125*u^2/4 - 325*u + 50/3)/j), + (value=-25*u*(1 - 5*u)*(1 - u)*(3 - 5*u)*(4 - 5*u)/12, + derivative=(15625*u^4/12 - 8125*u^3/3 + 7375*u^2/4 - 2675*u/6 + 25)/j), + (value=25*u*(1 - u)*(2 - 5*u)*(3 - 5*u)*(4 - 5*u)/24, + derivative=-(15625*u^4/24 - 4375*u^3/3 + 8875*u^2/8 - 1925*u/6 + 25)/j) + ) +end + +function (f::LagrangeRefSpace{T,6,2})(mp) where T + u = mp.bary[1] + j = jacobian(mp) + SVector( + (value=-u*(1 - 6*u)*(1 - 3*u)*(1 - 2*u)*(2 - 3*u)*(5 - 6*u)/10, + derivative=-(1944*u^5/5 - 810*u^4 + 612*u^3 - 405*u^2/2 + 137*u/5 - 1)/j), + (value=(1 - 6*u)*(1 - 3*u)*(1 - 2*u)*(1 - u)*(2 - 3*u)*(5 - 6*u)/10, + derivative=(1944*u^5/5 - 1134*u^4 + 1260*u^3 - 1323*u^2/2 + 812*u/5 - 147/10)/j), + (value=18*u*(1 - 6*u)*(1 - 3*u)*(1 - 2*u)*(1 - u)*(2 - 3*u)/5, + derivative=-( -11664*u^5/5 + 5184*u^4 - 4104*u^3 + 1404*u^2 - 972*u/5 + 36/5)/j), + (value=-9*u*(1 - 6*u)*(1 - 3*u)*(1 - 2*u)*(1 - u)*(5 - 6*u)/2, + derivative=-(5832*u^5 - 13770*u^4 + 11556*u^3 - 8289*u^2/2 + 594*u - 45/2)/j), + (value=4*u*(1 - 6*u)*(1 - 3*u)*(1 - u)*(2 - 3*u)*(5 - 6*u), + derivative=-( -7776*u^5 + 19440*u^4 - 17424*u^3 + 6696*u^2 - 1016*u + 40)/j), + (value=-9*u*(1 - 6*u)*(1 - 2*u)*(1 - u)*(2 - 3*u)*(5 - 6*u)/2, + derivative=-(5832*u^5 - 15390*u^4 + 14796*u^3 - 12447*u^2/2 + 1053*u - 45)/j), + (value=18*u*(1 - 3*u)*(1 - 2*u)*(1 - u)*(2 - 3*u)*(5 - 6*u)/5, + derivative=-( -11664*u^5/5 + 6480*u^4 - 6696*u^3 + 3132*u^2 - 3132*u/5 + 36)/j) + ) +end + +function (f::LagrangeRefSpace{T,7,2})(mp) where T + u = mp.bary[1] + j = jacobian(mp) + SVector( + (value=u*(7*u - 6)*(7*u - 5)*(7*u - 4)*(7*u - 3)*(7*u - 2)*(7*u - 1)/720, + derivative=-(823543*u^6/720 - 117649*u^5/40 + 420175*u^4/144 - 16807*u^3/12 + 9947*u^2/30 - 343*u/10 + 1)/j), + (value=-(u - 1)*(7*u - 6)*(7*u - 5)*(7*u - 4)*(7*u - 3)*(7*u - 2)*(7*u - 1)/720, + derivative=(823543*u^6/720 - 117649*u^5/30 + 386561*u^4/72 - 33614*u^3/9 + 331681*u^2/240 - 22981*u/90 + 363/20)/j), + (value=-49*u*(u - 1)*(7*u - 5)*(7*u - 4)*(7*u - 3)*(7*u - 2)*(7*u - 1)/720, + derivative=(5764801*u^6/720 - 1294139*u^5/60 + 1596665*u^4/72 - 98441*u^3/9 + 634207*u^2/240 - 49931*u/180 + 49/6)/j), + (value=49*u*(u - 1)*(7*u - 6)*(7*u - 4)*(7*u - 3)*(7*u - 2)*(7*u - 1)/240, + derivative=-(5764801*u^6/240 - 2705927*u^5/40 + 1159683*u^4/16 - 444185*u^3/12 + 45962*u^2/5 - 9849*u/10 + 147/5)/j), + (value=-49*u*(u - 1)*(7*u - 6)*(7*u - 5)*(7*u - 3)*(7*u - 2)*(7*u - 1)/144, + derivative=(5764801*u^6/144 - 117649*u^5 + 9495955*u^4/72 - 211288*u^3/3 + 872935*u^2/48 - 2009*u + 245/4)/j), + (value=49*u*(u - 1)*(7*u - 6)*(7*u - 5)*(7*u - 4)*(7*u - 2)*(7*u - 1)/144, + derivative=-(5764801*u^6/144 - 2941225*u^5/24 + 20756645*u^4/144 - 2926819*u^3/36 + 133427*u^2/6 - 46501*u/18 + 245/3)/j), + (value=-49*u*(u - 1)*(7*u - 6)*(7*u - 5)*(7*u - 4)*(7*u - 3)*(7*u - 1)/240, + derivative=(5764801*u^6/240 - 1529437*u^5/20 + 756315*u^4/8 - 170471*u^3/3 + 1347647*u^2/80 - 43071*u/20 + 147/2)/j), + (value=49*u*(u - 1)*(7*u - 6)*(7*u - 5)*(7*u - 4)*(7*u - 3)*(7*u - 2)/720, + derivative=-(5764801*u^6/720 - 1058841*u^5/40 + 4958065*u^4/144 - 88837*u^3/4 + 109417*u^2/15 - 10927*u/10 + 49)/j) + ) +end + +function (f::LagrangeRefSpace{T,8,2})(mp) where T + u = mp.bary[1] + j = jacobian(mp) + SVector( + (value=u*(2*u - 1)*(4*u - 3)*(4*u - 1)*(8*u - 7)*(8*u - 5)*(8*u - 3)*(8*u - 1)/315, + derivative=-(1048576*u^7/315 - 458752*u^6/45 + 188416*u^5/15 - 71680*u^4/9 + 123776*u^3/45 - 7504*u^2/15 + 1452*u/35 - 1)/j), + (value=(u - 1)*(2*u - 1)*(4*u - 3)*(4*u - 1)*(8*u - 7)*(8*u - 5)*(8*u - 3)*(8*u - 1)/315, + derivative=-(1048576*u^7/315 - 65536*u^6/5 + 106496*u^5/5 - 18432*u^4 + 136832*u^3/15 - 12816*u^2/5 + 118124*u/315 - 761/35)/j), + (value=-64*u*(u - 1)*(2*u - 1)*(4*u - 3)*(4*u - 1)*(8*u - 5)*(8*u - 3)*(8*u - 1)/315, + derivative=(8388608*u^7/315 - 3801088*u^6/45 + 1605632*u^5/15 - 624640*u^4/9 + 1097728*u^3/45 - 67456*u^2/15 + 13184*u/35 - 64/7)/j), + (value=16*u*(u - 1)*(2*u - 1)*(4*u - 1)*(8*u - 7)*(8*u - 5)*(8*u - 3)*(8*u - 1)/45, + derivative=-(4194304*u^7/45 - 917504*u^6/3 + 1998848*u^5/5 - 266240*u^4 + 1435136*u^3/15 - 17952*u^2 + 68576*u/45 - 112/3)/j), + (value=-64*u*(u - 1)*(2*u - 1)*(4*u - 3)*(4*u - 1)*(8*u - 7)*(8*u - 3)*(8*u - 1)/45, + derivative=(8388608*u^7/45 - 28442624*u^6/45 + 12812288*u^5/15 - 5285888*u^4/9 + 9773056*u^3/45 - 626048*u^2/15 + 18048*u/5 - 448/5)/j), + (value=4*u*(u - 1)*(4*u - 3)*(4*u - 1)*(8*u - 7)*(8*u - 5)*(8*u - 3)*(8*u - 1)/9, + derivative=-(2097152*u^7/9 - 7340032*u^6/9 + 3424256*u^5/3 - 7331840*u^4/9 + 2814208*u^3/9 - 186496*u^2/3 + 5528*u - 140)/j), + (value=-64*u*(u - 1)*(2*u - 1)*(4*u - 3)*(4*u - 1)*(8*u - 7)*(8*u - 5)*(8*u - 1)/45, + derivative=(8388608*u^7/45 - 10092544*u^6/15 + 4882432*u^5/5 - 727040*u^4 + 4390912*u^3/15 - 306048*u^2/5 + 256384*u/45 - 448/3)/j), + (value=16*u*(u - 1)*(2*u - 1)*(4*u - 3)*(8*u - 7)*(8*u - 5)*(8*u - 3)*(8*u - 1)/45, + derivative=-(4194304*u^7/45 - 15597568*u^6/45 + 7831552*u^5/15 - 3665920*u^4/9 + 7827968*u^3/45 - 587296*u^2/15 + 19872*u/5 - 112)/j), + (value=-64*u*(u - 1)*(2*u - 1)*(4*u - 3)*(4*u - 1)*(8*u - 7)*(8*u - 5)*(8*u - 3)/315, + derivative=(8388608*u^7/315 - 917504*u^6/9 + 2392064*u^5/15 - 1177600*u^4/9 + 2695168*u^3/45 - 44672*u^2/3 + 61568*u/35 - 64)/j) + ) +end + +function (f::LagrangeRefSpace{T,9,2})(mp) where T + u = mp.bary[1] + j = jacobian(mp) + SVector( + (value=u*(3*u - 2)*(3*u - 1)*(9*u - 8)*(9*u - 7)*(9*u - 5)*(9*u - 4)*(9*u - 2)*(9*u - 1)/4480, + derivative=-(43046721*u^8/4480 - 4782969*u^7/140 + 16120377*u^6/320 - 1594323*u^5/40 + 2337903*u^4/128 - 194643*u^3/40 + 797337*u^2/1120 - 6849*u/140 + 1)/j), + (value=-(u - 1)*(3*u - 2)*(3*u - 1)*(9*u - 8)*(9*u - 7)*(9*u - 5)*(9*u - 4)*(9*u - 2)*(9*u - 1)/4480, + derivative=(43046721*u^8/4480 - 4782969*u^7/112 + 5137263*u^6/64 - 2657205*u^5/32 + 6589431*u^4/128 - 623295*u^3/32 + 122121*u^2/28 - 58635*u/112 + 7129/280)/j), + (value=-81*u*(u - 1)*(3*u - 2)*(3*u - 1)*(9*u - 7)*(9*u - 5)*(9*u - 4)*(9*u - 2)*(9*u - 1)/4480, + derivative=(387420489*u^8/4480 - 176969853*u^7/560 + 152523567*u^6/320 - 61470009*u^5/160 + 22878207*u^4/128 - 7712091*u^3/160 + 3986901*u^2/560 - 275967*u/560 + 81/8)/j), + (value=81*u*(u - 1)*(3*u - 2)*(3*u - 1)*(9*u - 8)*(9*u - 5)*(9*u - 4)*(9*u - 2)*(9*u - 1)/1120, + derivative=-(387420489*u^8/1120 - 90876411*u^7/70 + 80247591*u^6/40 - 66075831*u^5/40 + 25043337*u^4/32 - 2142531*u^3/10 + 8968887*u^2/280 - 78327*u/35 + 324/7)/j), + (value=-9*u*(u - 1)*(3*u - 1)*(9*u - 8)*(9*u - 7)*(9*u - 5)*(9*u - 4)*(9*u - 2)*(9*u - 1)/160, + derivative=(129140163*u^8/160 - 62178597*u^7/20 + 197164611*u^6/40 - 166341033*u^5/40 + 64448703*u^4/32 - 22480173*u^3/40 + 6828867*u^2/80 - 60381*u/10 + 126)/j), + (value=81*u*(u - 1)*(3*u - 2)*(3*u - 1)*(9*u - 8)*(9*u - 7)*(9*u - 4)*(9*u - 2)*(9*u - 1)/320, + derivative=-(387420489*u^8/320 - 4782969*u^7 + 249245829*u^6/32 - 54029835*u^5/8 + 215023653*u^4/64 - 3844017*u^3/4 + 2386017*u^2/16 - 21465*u/2 + 1134/5)/j), + (value=-81*u*(u - 1)*(3*u - 2)*(3*u - 1)*(9*u - 8)*(9*u - 7)*(9*u - 5)*(9*u - 2)*(9*u - 1)/320, + derivative=(387420489*u^8/320 - 196101729*u^7/40 + 1313190711*u^6/160 - 586888011*u^5/80 + 241241409*u^4/64 - 89119521*u^3/80 + 14257053*u^2/80 - 526419*u/40 + 567/2)/j), + (value=9*u*(u - 1)*(3*u - 2)*(9*u - 8)*(9*u - 7)*(9*u - 5)*(9*u - 4)*(9*u - 2)*(9*u - 1)/160, + derivative=-(129140163*u^8/160 - 33480783*u^7/10 + 115322697*u^6/20 - 213107841*u^5/40 + 91020753*u^4/32 - 8776431*u^3/10 + 5878089*u^2/40 - 56601*u/5 + 252)/j), + (value=-81*u*(u - 1)*(3*u - 2)*(3*u - 1)*(9*u - 8)*(9*u - 7)*(9*u - 5)*(9*u - 4)*(9*u - 1)/1120, + derivative=(387420489*u^8/1120 - 205667667*u^7/140 + 26040609*u^6/10 - 99733761*u^5/40 + 44529507*u^4/32 - 18152829*u^3/40 + 45570519*u^2/560 - 475389*u/70 + 162)/j), + (value=81*u*(u - 1)*(3*u - 2)*(3*u - 1)*(9*u - 8)*(9*u - 7)*(9*u - 5)*(9*u - 4)*(9*u - 2)/4480, + derivative=-(387420489*u^8/4480 - 52612659*u^7/140 + 219485133*u^6/320 - 13640319*u^5/20 + 51221727*u^4/128 - 5589243*u^3/40 + 30921993*u^2/1120 - 373329*u/140 + 81)/j) + ) +end + +function (f::LagrangeRefSpace{T,10,2})(mp) where T + u = mp.bary[1] + j = jacobian(mp) + SVector( + (value=u*(2*u - 1)*(5*u - 4)*(5*u - 3)*(5*u - 2)*(5*u - 1)*(10*u - 9)*(10*u - 7)*(10*u - 3)*(10*u - 1)/4536, + derivative=-(15625000*u^9/567 - 781250*u^8/7 + 36250000*u^7/189 - 546875*u^6/3 + 1883125*u^5/18 - 296875*u^4/8 + 4523000*u^3/567 - 162875*u^2/168 + 7129*u/126 - 1)/j), + (value=(u - 1)*(2*u - 1)*(5*u - 4)*(5*u - 3)*(5*u - 2)*(5*u - 1)*(10*u - 9)*(10*u - 7)*(10*u - 3)*(10*u - 1)/4536, + derivative=-(15625000*u^9/567 - 8593750*u^8/63 + 55000000*u^7/189 - 9453125*u^6/27 + 4695625*u^5/18 - 26846875*u^4/216 + 42711625*u^3/1134 - 10511875*u^2/1512 + 177133*u/252 - 7381/252)/j), + (value=-25*u*(u - 1)*(2*u - 1)*(5*u - 4)*(5*u - 3)*(5*u - 2)*(5*u - 1)*(10*u - 7)*(10*u - 3)*(10*u - 1)/1134, + derivative=(156250000*u^9/567 - 71875000*u^8/63 + 377500000*u^7/189 - 52062500*u^6/27 + 10090625*u^5/9 - 21709375*u^4/54 + 49435250*u^3/567 - 4033825*u^2/378 + 13150*u/21 - 100/9)/j), + (value=25*u*(u - 1)*(2*u - 1)*(5*u - 3)*(5*u - 2)*(5*u - 1)*(10*u - 9)*(10*u - 7)*(10*u - 3)*(10*u - 1)/504, + derivative=-(78125000*u^9/63 - 36718750*u^8/7 + 590000000*u^7/63 - 82796875*u^6/9 + 32584375*u^5/6 - 142028125*u^4/72 + 54486625*u^3/126 - 2990025*u^2/56 + 88325*u/28 - 225/4)/j), + (value=-50*u*(u - 1)*(2*u - 1)*(5*u - 4)*(5*u - 3)*(5*u - 2)*(5*u - 1)*(10*u - 9)*(10*u - 3)*(10*u - 1)/189, + derivative=(625000000*u^9/189 - 100000000*u^8/7 + 1640000000*u^7/63 - 234625000*u^6/9 + 15662500*u^5 - 52006250*u^4/9 + 242639000*u^3/189 - 1121950*u^2/7 + 200600*u/21 - 1200/7)/j), + (value=25*u*(u - 1)*(2*u - 1)*(5*u - 4)*(5*u - 2)*(5*u - 1)*(10*u - 9)*(10*u - 7)*(10*u - 3)*(10*u - 1)/108, + derivative=-(156250000*u^9/27 - 76562500*u^8/3 + 47500000*u^7 - 437281250*u^6/9 + 89384375*u^5/3 - 403334375*u^4/36 + 68357750*u^3/27 - 11544725*u^2/36 + 174025*u/9 - 350)/j), + (value=-u*(u - 1)*(5*u - 4)*(5*u - 3)*(5*u - 2)*(5*u - 1)*(10*u - 9)*(10*u - 7)*(10*u - 3)*(10*u - 1)/9, + derivative=(62500000*u^9/9 - 31250000*u^8 + 535000000*u^7/9 - 560000000*u^6/9 + 117216250*u^5/3 - 135371875*u^4/9 + 31274500*u^3/9 - 448875*u^2 + 27508*u - 504)/j), + (value=25*u*(u - 1)*(2*u - 1)*(5*u - 4)*(5*u - 3)*(5*u - 1)*(10*u - 9)*(10*u - 7)*(10*u - 3)*(10*u - 1)/108, + derivative=-(156250000*u^9/27 - 26562500*u^8 + 155000000*u^7/3 - 498968750*u^6/9 + 107321875*u^5/3 - 510353125*u^4/36 + 91073375*u^3/27 - 1792225*u^2/4 + 168775*u/6 - 525)/j), + (value=-50*u*(u - 1)*(2*u - 1)*(5*u - 4)*(5*u - 3)*(5*u - 2)*(5*u - 1)*(10*u - 9)*(10*u - 7)*(10*u - 1)/189, + derivative=(625000000*u^9/189 - 325000000*u^8/21 + 1940000000*u^7/63 - 305375000*u^6/9 + 67737500*u^5/3 - 83431250*u^4/9 + 433739000*u^3/189 - 20028950*u^2/63 + 1308200*u/63 - 400)/j), + (value=25*u*(u - 1)*(2*u - 1)*(5*u - 4)*(5*u - 3)*(5*u - 2)*(10*u - 9)*(10*u - 7)*(10*u - 3)*(10*u - 1)/504, + derivative=-(78125000*u^9/63 - 41406250*u^8/7 + 758750000*u^7/63 - 122828125*u^6/9 + 56396875*u^5/6 - 289909375*u^4/72 + 66191750*u^3/63 - 8694225*u^2/56 + 153025*u/14 - 225)/j), + (value=-25*u*(u - 1)*(2*u - 1)*(5*u - 4)*(5*u - 3)*(5*u - 2)*(5*u - 1)*(10*u - 9)*(10*u - 7)*(10*u - 3)/1134, + derivative=(156250000*u^9/567 - 9375000*u^8/7 + 527500000*u^7/189 - 29312500*u^6/9 + 20965625*u^5/9 - 18878125*u^4/18 + 165985250*u^3/567 - 1997825*u^2/42 + 243050*u/63 - 100)/j), + ) +end + # Evaluete linear lagrange elements on a triangle function (f::LagrangeRefSpace{T,1,3})(t) where T u,v,w, = barycentric(t) diff --git a/src/interpolation.jl b/src/interpolation.jl index a20a69e9..ad520572 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -35,6 +35,25 @@ function DofInterpolate(basis::LagrangeBasis, field) return res end + +function DofInterpolate(basis::LagrangeBasis{D,C,M}, field) where {D,C,U,M<:CompScienceMeshes.Mesh{U,2}} + + T = promote_type(scalartype(basis), scalartype(field)) + + num_bfs = numfunctions(basis) + + res = Vector{T}(undef, num_bfs) + + for b in 1 : num_bfs + bfs = basis.fns[b] + + res[b] = field(basis.pos[b]) + end + + return res +end + + ### Piecewise constant elements require separate treatment ### TODO: Probably a rewrite is advisable to also take into account ### dual elements properly. diff --git a/test/runtests.jl b/test/runtests.jl index 372faaf4..ce0cd25b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,6 +18,7 @@ include("test_fourier.jl") include("test_specials.jl") include("test_basis.jl") +include("test_lagrange.jl") include("test_directproduct.jl") include("test_raviartthomas.jl") include("test_rt.jl") @@ -56,6 +57,7 @@ include("test_dipole.jl") include("test_sauterschwabints1D.jl") include("test_hh2d_exec.jl") include("test_hh2d_mie.jl") +include("test_hh2d_mie_higher_order.jl") include("test_hh2d_nearfield.jl") include("test_wiltonints.jl") diff --git a/test/test_hh2d_mie_higher_order.jl b/test/test_hh2d_mie_higher_order.jl new file mode 100644 index 00000000..83b9ef22 --- /dev/null +++ b/test/test_hh2d_mie_higher_order.jl @@ -0,0 +1,262 @@ +using BEAST +using CompScienceMeshes +using StaticArrays +using LinearAlgebra +using Test +using SpecialFunctions + + +# We check only +for order = [2; 3] + + ε0 = 8.854187821e-12 + μ0 = 4π*1e-7 + c0 = 1/sqrt(ε0*μ0) + η0 = sqrt(μ0/ε0) + + f = 1e9 # 1 GHz + ω = 2π*f + λ = c0/f + k = 2π/λ + h = λ/10 + a = 1.0 # radius of the scatterer + + circle = CompScienceMeshes.meshcircle(a, h) + X0 = lagrangecx(circle, order=order) + X1 = lagrangec0(circle, order=order) + + # Computing the fields on a circle of radius rc + rc = 500.0*a + pts = meshcircle(rc, 0.3 * rc).vertices + + ## Analytical solutions + + dbesselj(n,x) = (besselj(n-1, x) - besselj(n+1, x)) / 2 + dhankelh2(n,x) = (hankelh2(n-1, x) - hankelh2(n+1, x)) / 2 + + cart2polar(x,y) = SVector(sqrt(x^2 + y^2), atan(y, x)) + + + + +## Analytical Solutions + + # TM planewave solution (based on Sec 6.4.2, Jin, Theory and Computation of Electromagnetic Fields]) + function TM_pec_planewave(E0, k, a, ρ, φ) + n = 0 + + # Jin (6.4.11) + a_n(n) = -(1.0im)^(-n)*(besselj(n, k*a) / hankelh2(n, k*a)) + valEz(n) = a_n(n) * hankelh2(n, k*ρ)*exp(im*n*φ) + valHφ(n) = 1/(im * ω * μ0) * k * a_n(n) * dhankelh2(n,k*ρ) * exp(im*n*φ) + valHρ(n) = -1/(im * ω * μ0) * (im * n / ρ) * a_n(n) * hankelh2(n, k*ρ) * exp(im*n*φ) + + retEz = valEz(0) + retHφ = valHφ(0) + retHρ = valHρ(0) + + while true + n += 1 + bufEz = valEz(n) + valEz(-n) + if abs(bufEz) >= abs(retEz) * eps(eltype(real(retEz))) * 1e3 + retEz += bufEz + else + break + end + end + + n=0 + while true + n += 1 + bufHφ = valHφ(n) + valHφ(-n) + if abs(bufHφ) >= abs(retHφ) * eps(eltype(real(retHφ))) * 1e3 + retHφ += bufHφ + else + break + end + end + + n=0 + while true + n += 1 + bufHρ = valHρ(n) + valHρ(-n) + if abs(bufHρ) >= abs(retHρ) * eps(eltype(real(retHρ))) * 1e3 + retHρ += bufHρ + else + break + end + end + + return E0 * retEz, E0*(retHρ*SVector(cos(φ), sin(φ)) + retHφ*SVector(-sin(φ), cos(φ))) + end + + function TM_pec_planewave_E(E0, k, a, pts) + return [TM_pec_planewave(E0, k, a, cart2polar(p[1], p[2])...)[1] for p in pts] + end + + function TM_pec_planewave_H(E0, k, a, pts) + return [TM_pec_planewave(E0, k, a, cart2polar(p[1], p[2])...)[2] for p in pts] + end + + + # TE planewave solution + function TE_pec_planewave(H0, k, a, ρ, φ) + n = 0 + + # Jin (6.4.19) + b(n) = -(1.0im)^(-n)*dbesselj(n,k*a)/dhankelh2(n,k*a) + + # Hz is not needed, but maybe we use it later for MFIE + valHz(n) = b(n) * hankelh2(n, k*ρ)*exp(im*n*φ) + valEφ(n) = -1/(im * ω * ε0) * k * b(n) * dhankelh2(n, k*ρ)*exp(im*n*φ) + valEρ(n) = 1/(im * ω * ε0) * (im * n / ρ) * b(n) * hankelh2(n, k*ρ)*exp(im*n*φ) + + retHz = valHz(0) + retEφ = valEφ(0) + retEρ = valEρ(0) + + while true + n += 1 + bufHz = valHz(n) + valHz(-n) + if abs(bufHz) >= abs(retHz) * eps(eltype(real(retHz))) * 1e3 + retHz += bufHz + else + break + end + end + + n = 0 + while true + n += 1 + bufEφ = valEφ(n) + valEφ(-n) + if abs(bufEφ) >= abs(retEφ) * eps(eltype(real(retEφ))) * 1e3 + retEφ += bufEφ + else + break + end + end + + n = 0 + while true + n += 1 + bufEρ = valEρ(n) + valEρ(-n) + if abs(bufEρ) >= abs(retEρ) * eps(eltype(real(retEρ))) * 1e3 + retEρ += bufEρ + else + break + end + end + + return H0*retHz, H0*(retEρ*SVector(cos(φ), sin(φ)) + retEφ*SVector(-sin(φ), cos(φ))) + end + + function TE_pec_planewave_H(H0, k, a, pts) + return [TE_pec_planewave(H0, k, a, cart2polar(p[1],p[2])...)[1] for p in pts] + end + + function TE_pec_planewave_E(H0, k, a, pts) + return [TE_pec_planewave(H0, k, a, cart2polar(p[1],p[2])...)[2] for p in pts] + end + + + + +## TM Scattering + + # TM-EFIE system matrix + 𝒮 = Helmholtz2D.singlelayer(; alpha=im*k*η0, wavenumber=k) + M_TMEFIE = assemble(𝒮, X0, X0) + + # TM-MFIE system matrix + 𝒟ᵀ = Helmholtz2D.doublelayer_transposed(; wavenumber=k) + Dᵀ = assemble(𝒟ᵀ, X0, X0) + + I0 = assemble(BEAST.Identity(), X0, X0) + M_TMMFIE = +0.5*I0 + Dᵀ + + # 1. Excitation: Plane wave (The incoming planewave is along the x-axis) + E0 = 1.0 # Amplitude + Ez_pw_inc = Helmholtz2D.planewave(; amplitude=E0, wavenumber=k, direction=SVector(1.0, 0.0)) + + ez_pw_inc = assemble(DirichletTrace(Ez_pw_inc), X0) + j_TMEFIE_pw = M_TMEFIE \ ez_pw_inc + + # TM-EFIE: We compute the scattered Ez component (scalar) + Ez_pw_sca_num = -potential(HH2DSingleLayerNear(𝒮), pts, j_TMEFIE_pw, X0; type=ComplexF64) + Ez_pw_sca_ana = TM_pec_planewave_E(E0, k, a, pts) + + # Will not improve until we have curvilinear elements + # But must not become worse either + # TODO: Once curvilinear elements are support: refine order of functions and + # mesh in a lock step and update these values + Ez_pw_sca_bounds = [0.002; 0.002; 0.0014] + @test norm(Ez_pw_sca_ana - Ez_pw_sca_num) ./ norm(Ez_pw_sca_ana) < Ez_pw_sca_bounds[order] + + Ht_pw_sca_num = -potential(HH2DDoubleLayerTransposedNear(𝒟ᵀ), pts, j_TMEFIE_pw, X0; type=SVector{2,ComplexF64}) + Ht_pw_sca_ana = TM_pec_planewave_H(E0, k, a, pts) + + Ht_pw_sca_bounds = [0.002; 0.002; 0.0015] + @test norm(Ref(ẑ) .× Ht_pw_sca_num - Ht_pw_sca_ana) ./ norm(Ht_pw_sca_ana) < Ht_pw_sca_bounds[order] + + # TM-MFIE + Ht_pw_inc = - 1.0 / (im * ω * μ0) * curl(Ez_pw_inc) + ht_pw_inc = -assemble(TangentTrace(Ht_pw_inc), X0) + + j_TMMFIE_pw = M_TMMFIE \ ht_pw_inc + + j_pw_sca_bounds = [0.05; 0.001; 0.0009] + @test norm(j_TMMFIE_pw - j_TMEFIE_pw) / norm(j_TMEFIE_pw) < j_pw_sca_bounds[order] + + + +## TE Scattering + + I = 1.0 + ρp = 2.0 + φp = 0.0 + r_inc = 1.0 + + pts_inc = meshcircle(r_inc, 0.6*r_inc).vertices + + # Scattered Field + # TE-MFIE + 𝒟 = Helmholtz2D.doublelayer(; wavenumber=k) + D = assemble(𝒟, X1, X1) + I1 = assemble(BEAST.Identity(), X1, X1) + M_TEMFIE = 0.5I1 - D + + # TE-EFIE + 𝒩 = Helmholtz2D.hypersingular(;alpha=im*ω*μ0, beta=1.0/(im*ω*ε0), wavenumber=k) + M_TEEFIE = assemble(𝒩, X1, X1) + + # 1. Excitation: Planewave + H0 = 1.0 + Hz_pw_inc = Helmholtz2D.planewave(; amplitude=H0, wavenumber=k, direction=SVector(1.0, 0.0)) + hz_pw_inc = assemble(DirichletTrace(Hz_pw_inc), X1) + j_TEMFIE_pw = M_TEMFIE \ hz_pw_inc + + Hz_pw_sca_num = potential(HH2DDoubleLayerNear(𝒟), pts, j_TEMFIE_pw, X1; type=ComplexF64) + Hz_pw_sca_ana = TE_pec_planewave_H(H0, k, a, pts) + + Hz_pw_sca_bounds = [0.0015; 0.0015; 0.0015] + @test norm(Hz_pw_sca_num - Hz_pw_sca_ana) /norm(Hz_pw_sca_ana) < Hz_pw_sca_bounds[order] + + Et_pw_inc = 1 / (im * ω * ε0) * curl(Hz_pw_inc) + @test Et_pw_inc.direction == Hz_pw_inc.direction + @test Et_pw_inc.polarization == SVector(0.0, 1.0) + @test Et_pw_inc.amplitude ≈ Et_pw_inc.gamma / (im * ω * ε0) atol=1e-15 + et_pw_inc = assemble(TangentTrace(Et_pw_inc), X1) + + j_TEEFIE_pw = M_TEEFIE \ et_pw_inc + + j_pw_sca_bounds = [0.007; 0.0002; 1.19e-5] + @test norm(j_TEEFIE_pw - j_TEMFIE_pw)/norm(j_TEMFIE_pw) < j_pw_sca_bounds[order] + + Et_pw_sca_num = -potential(HH2DHyperSingularNear(𝒩), pts, j_TEEFIE_pw, X1; type=SVector{2, ComplexF64}) + Et_pw_sca_ana = TE_pec_planewave_E(H0, k, a, pts) + + # We compute the scattered Ez component (scalar) + Et_pw_sca_bounds = [0.003; 0.0016; 0.0016] + @test norm(Ref(ẑ) .× Et_pw_sca_num - Et_pw_sca_ana) / norm(Et_pw_sca_ana) < Et_pw_sca_bounds[order] + +end \ No newline at end of file diff --git a/test/test_lagrange.jl b/test/test_lagrange.jl new file mode 100644 index 00000000..6ba9b7d2 --- /dev/null +++ b/test/test_lagrange.jl @@ -0,0 +1,84 @@ +using Test +using BEAST +using CompScienceMeshes + +@testset "Higher order 1D Lagrange Elements" begin + line = meshsegment(2.0, 0.5) + numsegments = CompScienceMeshes.numcells(line) + g = BEAST.ScalarTrace{Float64}(x -> x[1]^5) + glbf = BEAST.GlobalFunction(g, line, Vector(1:numcells(line))) + @test BEAST.Lp_integrate(glbf; p=1) ≈ 10.66666666666 + + + @testset "C0" begin + + for order in 1:10 + C0 = lagrangec0(line, order=order, dirichlet=false) + @test numfunctions(C0) == (order-1)*numsegments + numsegments + 1 + end + + # We integrate x^5 from 0 to 2 + # Analytic outcome: 32/3 + + # Order 1 + C0 = lagrangec0(line, order=1, dirichlet=false) + coeffs = DofInterpolate(C0, g) + gfc0 = BEAST.FEMFunction(coeffs, C0) + @test BEAST.Lp_integrate(gfc0; p=1) ≈ 12.3125 + + # Order 2 + C0 = lagrangec0(line, order=2, dirichlet=false) + coeffs = DofInterpolate(C0, g) + gfc0 = BEAST.FEMFunction(coeffs, C0) + @test BEAST.Lp_integrate(gfc0; p=1) ≈ 10.672838520122093 + + # Order 3 + C0 = lagrangec0(line, order=3, dirichlet=false) + coeffs = DofInterpolate(C0, g) + gfc0 = BEAST.FEMFunction(coeffs, C0) + @test BEAST.Lp_integrate(gfc0; p=1) ≈ 10.6689814814888 + + @test C0.pos[6][1] ≈ 0.166666666666666 + + # Order 4 + C0 = lagrangec0(line, order=4, dirichlet=false) + coeffs = DofInterpolate(C0, g) + gfc0 = BEAST.FEMFunction(coeffs, C0) + @test BEAST.Lp_integrate(gfc0; p=1) ≈ 10.66668405 + + # Order 5 + C0 = lagrangec0(line, order=5, dirichlet=false) + coeffs = DofInterpolate(C0, g) + gfc0 = BEAST.FEMFunction(coeffs, C0) + @test BEAST.Lp_integrate(gfc0; p=1) ≈ 10.6666666666 + + @test C0.pos[6][1] ≈ 0.1 + end + ## + + @testset "CX" begin + + for order in 1:10 + CX = lagrangecx(line, order=order) + @test numfunctions(CX) == (order+1)*numsegments + end + + # Order 1 + CX = lagrangecx(line, order=1) + coeffs = DofInterpolate(CX, g) + gfc0 = BEAST.FEMFunction(coeffs, CX) + @test BEAST.Lp_integrate(gfc0; p=1) ≈ 12.3125 + + CX = lagrangecx(line, order=2) + coeffs = DofInterpolate(CX, g) + gfc0 = BEAST.FEMFunction(coeffs, CX) + @test BEAST.Lp_integrate(gfc0; p=1) ≈ 10.672838520122093 + + CX = lagrangecx(line, order=5) + coeffs = DofInterpolate(CX, g) + gfc0 = BEAST.FEMFunction(coeffs, CX) + @test BEAST.Lp_integrate(gfc0; p=1) ≈ 10.6666666666 + + @test CX.pos[9] ≈ point(0.1, 0.0) + end +end \ No newline at end of file From 5ab9e1261504794fcc592ce33d22c1ddf1a50ebb Mon Sep 17 00:00:00 2001 From: "Simon B. Adrian" Date: Sat, 8 Nov 2025 17:49:25 +0100 Subject: [PATCH 2/2] Remove historical type instabilities Greatly improves the speed. I expect further improvements are possible, but this was the most significant type instability that I have found. --- src/helmholtz2d/hh2dops.jl | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/helmholtz2d/hh2dops.jl b/src/helmholtz2d/hh2dops.jl index f9486089..a467664e 100644 --- a/src/helmholtz2d/hh2dops.jl +++ b/src/helmholtz2d/hh2dops.jl @@ -102,14 +102,13 @@ struct HH2DHyperSingularFDBIO{T, K} <: HelmholtzOperator2D{T, K} end end -mutable struct KernelValsHelmholtz2D - gamma - wavenumber - vect - dist - green - gradgreen - txty +mutable struct KernelValsHelmholtz2D{K1,K2,F} + gamma::K1 + vect::SVector{2,F} + dist::F + green::K2 + gradgreen::SVector{2,K2} + txty::F end # Kernel values for 2D Laplace BIE operators @@ -148,7 +147,7 @@ function kernelvals(biop::HelmholtzOperator2D{T, K}, tgeo, bgeo) where {T, K <: txty = dot(normal(tgeo), normal(bgeo)) - KernelValsHelmholtz2D(gamma, k, r, R, green, gradgreen, txty) + KernelValsHelmholtz2D(gamma, r, R, green, gradgreen, txty) end shapevals(op::HelmholtzOperator2D, ϕ, ts) = shapevals(ValDiff(), ϕ, ts)