continuous higher order scalar spaces

examples/hh3d_dirichlet_order3.jl

@@ -14,13 +14,13 @@ X = BEAST.lagrangecx(Γ; order=3)
 a = Helmholtz3D.singlelayer(wavenumber=κ)
 # b = Helmholtz3D.doublelayer_transposed(gamma=κ*im) +0.5Identity()
 
-BEAST.@defaultquadstrat (a,X,X) BEAST.DoubleNumSauterQstrat(7,8,6,6,6,6)
-BEAST.@defaultquadstrat (f,X) BEAST.SingleNumQStrat(12)
-
 uⁱ = Helmholtz3D.planewave(wavenumber=κ, direction=z)
 f = strace(uⁱ,Γ)
 # g = ∂n(uⁱ)
 
+BEAST.@defaultquadstrat (a,X,X) BEAST.DoubleNumSauterQstrat(7,8,6,6,6,6)
+BEAST.@defaultquadstrat (f,X) BEAST.SingleNumQStrat(12)
+
 @hilbertspace u
 @hilbertspace v
 
@@ -32,7 +32,7 @@ b = assemble(f[v], X)
 x1 = AbstractMatrix(A) \ b
 # x1 = BEAST.GMRESSolver(A; reltol=1e-10) * b
 
-x1 = gmres(eq1; tol=1e-6)
+# x1 = gmres(eq1; tol=1e-6)
 # x2 = gmres(eq2)
 
 fcr1, geo1 = facecurrents(x1, X)

examples/hh3d_neumann order3.jl

@@ -0,0 +1,56 @@
+using CompScienceMeshes, BEAST
+using LinearAlgebra, Pkg
+
+# Pkg.activate(@__DIR__)
+Γ = readmesh(joinpath(dirname(pathof(BEAST)),"../examples/sphere2.in"))
+# Γ = meshrectangle(1.0, 1.0, 0.05, 3)
+X = BEAST.lagrangec0(Γ; order=2)
+@show numfunctions(X)
+
+κ = 2π; γ = im*κ
+a = Helmholtz3D.hypersingular(gamma=γ)
+# b = Helmholtz3D.doublelayer(gamma=γ) - 0.5Identity()
+
+uⁱ = Helmholtz3D.planewave(wavenumber=κ, direction=ẑ)
+# f = strace(uⁱ,Γ)
+g = ∂n(uⁱ)
+
+BEAST.@defaultquadstrat (a,X,X) BEAST.DoubleNumSauterQstrat(7,8,6,6,6,6)
+BEAST.@defaultquadstrat (g,X) BEAST.SingleNumQStrat(12)
+
+@hilbertspace u
+@hilbertspace v
+
+A = assemble(a[v,u], X, X)
+b = assemble(g[v], X)
+x1 = AbstractMatrix(A) \ b
+
+# eq1 = @discretise a[v,u] == g[v] u∈X v∈X
+# eq2 = @discretise b[v,u] == f[v] u∈X v∈X
+
+# x1 = gmres(eq1)
+# x2 = gmres(eq2)
+
+# @assert norm(x1-x2)/norm(x1+x2) <  0.5e-2
+
+fcr1, geo1 = facecurrents(x1, X)
+# fcr2, geo2 = facecurrents(x2, X)
+
+using Plots
+Plots.plot(title="Comparse 1st and 2nd kind eqs.")
+Plots.plot!(norm.(fcr1),c=:blue,label="1st")
+# Plots.scatter!(norm.(fcr2),c=:red,label="2nd")
+
+import Plotly
+plt1 = Plotly.plot(patch(Γ, norm.(fcr1)))
+# plt2 = Plotly.plot(patch(Γ, norm.(fcr2)))
+# display([plt1 plt2])
+
+# ys = range(-2,2,length=200)
+# zs = range(-2,2,length=200)
+# pts = [point(0.5, y, z) for y in ys, z in zs]
+
+# nf = BEAST.HH3DDoubleLayerNear(wavenumber=κ)
+# near = BEAST.potential(nf, pts, x1, X, type=ComplexF64)
+# inc = uⁱ.(pts)
+# tot = near + inc

src/BEAST.jl

@@ -147,6 +147,7 @@ include("bases/lincomb.jl")
 include("bases/trace.jl")
 include("bases/restrict.jl")
 include("bases/divergence.jl")
+include("bases/global/globaldofs.jl")
 
 include("bases/local/laglocal.jl")
 include("bases/local/rtlocal.jl")

src/bases/global/globaldofs.jl

@@ -0,0 +1,148 @@
+function trace(edge_ch, face_ch, localspace::RefSpace)
+    T = coordtype(edge_ch)
+    atol = sqrt(eps(T))
+
+    face_dom = domain(face_ch)
+
+    edge_verts = vertices(edge_ch)
+    face_verts = vertices(face_ch)
+
+    perm = [something(findfirst(w -> isapprox(v, w; atol), face_verts), 0) for v in edge_verts]
+    # @show perm
+    injection = simplex(vertices(face_dom)[perm])
+
+    function f(s)
+        u = neighborhood(injection, s)
+        p = neighborhood(face_ch, cartesian(u))
+        localspace(p)
+    end
+
+    return f
+end
+
+function _lagrangepolynomial(nodes, i, s, i1=length(nodes))
+    r = one(T)
+    si = nodes[i]
+    for j in 1:i1
+        j == i && continue
+        sj = nodes[j]
+        r *= (s - sj) / (si - sj)
+    end
+    return r
+end 
+
+struct _LagrangeGlobalEdgeDoFs
+    order::Int
+end
+function numfunctions(dof::_LagrangeGlobalEdgeDoFs) dof.order-1 end
+function (dof::_LagrangeGlobalEdgeDoFs)(s)
+    T = typeof(s)
+    nodes = range(zero(T), one(T), length=order+1)
+    [_lagrangepolynomial(nodes, i, s) for i in 2:order]
+end
+
+struct _LagrangeGlobalFaceDoFs
+    order::Int
+end
+# function numfunctions(dof::_LagrangeGlobalFaceDoFs) div((order-1)*(order-2),2) end
+# function (dof::_LagrangeGlobalFaceDoFs)(s)
+#     T = eltype(s)
+#     nodes = range(zero(T), one(T), length=order+1)
+#     r = zeros(T, numfunctions(dof))
+#     s1, s2 = s
+#     s3 = 1 - s1 - s2
+#     idx = 1
+#     degree = dof.degree
+#     for i in 0:degree
+#         prodi = _lagrangepolynomial(nodes, i+1, s1, i)
+#         for j in 0:degree
+#             k = degree - i - j
+#             k < 0 && continue
+#             prodj = _lagrangepolynomial(nodes, j+1, s2, j)
+#             prodk = _lagrangepolynomial(nodes, k+1, s3, k)
+#             r[idx] = prodi * prodj * prodk
+#     end end
+#     return r
+# end
+
+function globaldofs(edge_ch, face_ch, localspace, dof::_LagrangeGlobalEdgeDoFs)
+
+    T = coordtype(edge_ch)
+    f = trace(edge_ch, face_ch, localspace)
+    # r = zero(numfunctions(localspace), numfunctions(dof))
+    # for (s,w) in CompScienceMeshes.quadpoints(edge_dom, 2*order)
+    #     u = neighborhood(injection, s)
+    #     p = neighborhood(face_ch, cartesian(u))
+    #     r .+= w * [x.value for x in localspace(p)] * dof(s)'
+    # end
+
+    ds = one(T) / dof.order
+    return stack(range(ds, step=ds, length=dof.order-1)) do s
+        [x.value for x in f(s)]
+    end
+end
+
+
+function globaldofs(edge_ch, face_ch, localspace, dof::_LagrangeGlobalFaceDoFs)
+
+    T = coordtype(edge_ch)
+    f = trace(edge_ch, face_ch, localspace)
+    # r = zero(numfunctions(localspace), numfunctions(dof))
+    # for (s,w) in CompScienceMeshes.quadpoints(edge_dom, 2*order)
+    #     u = neighborhood(injection, s)
+    #     p = neighborhood(face_ch, cartesian(u))
+    #     r .+= w * [x.value for x in localspace(p)] * dof(s)'
+    # end
+
+    d = dof.order
+    S = ((i,j,d-i-j) for i in 0:d for j in 0:d if (i+j < d && i > 0 && j > 0))
+    return stack(S) do s
+        s = (s[1]/d,s[2]/d)
+        [x.value for x in f(s)]
+    end
+end
+
+struct _LagrangeGlobalNodesDoFs
+    order::Int
+end
+
+function globaldofs(edge_ch, face_ch, localspace, dof::_LagrangeGlobalNodesDoFs)
+    f = trace(edge_ch, face_ch, localspace)
+    return stack([()]) do s
+        [x.value for x in f(s)]
+    end
+end
+
+
+@testitem "globaldofs: interpolatory on edge" begin
+    using CompScienceMeshes
+
+    T = Float64
+    face_ch = simplex(
+        point(3,0,0),
+        point(2,0,-1),
+        point(0,0,1),
+    )
+    edge_ch = simplex(
+        point(0,0,1),
+        point(2,0,-1),
+    )
+
+    localspace = BEAST.LagrangeRefSpace{T,3,3,10}()
+    L = BEAST._LagrangeGlobalEdgeDoFs(3)
+    dofs = BEAST.globaldofs(edge_ch, face_ch, localspace, L)
+
+    @test size(dofs) == (10,2)
+    A = [
+        0 0
+        0 1
+        1 0
+        0 0
+        0 0
+        0 0
+        0 0
+        0 0
+        0 0
+        0 0]
+    @test dofs ≈ A
+end