divergence tested and efie seems to work

examples/efie_gwp2.jl

@@ -0,0 +1,31 @@
+using CompScienceMeshes
+using BEAST
+
+Γ = readmesh(joinpath(dirname(pathof(BEAST)),"../examples/sphere2.in"))
+# Γ = meshsphere(radius=1.0, h=0.1)
+@show length(Γ)
+# X = raviartthomas(Γ)
+X = BEAST.gwpdiv(Γ; order=2)
+
+κ, η = 1.0, 1.0
+t = Maxwell3D.singlelayer(wavenumber=κ)
+E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
+e = (n × E) × n
+
+BEAST.@defaultquadstrat (t,X,X) BEAST.DoubleNumSauterQstrat(7,8,6,6,6,6)
+BEAST.@defaultquadstrat (e,X) BEAST.SingleNumQStrat(12)
+
+@hilbertspace j
+@hilbertspace k
+# efie = @discretise t[k,j]==e[k]  j∈X k∈X
+# u, ch = BEAST.gmres_ch(efie; restart=1500)
+
+A = assemble(t[k,j], X, X)
+b = assemble(e[k], X)
+Ai = BEAST.GMRESSolver(A; reltol=1e-8, restart=1500)
+u = Ai * b
+
+include("utils/postproc.jl")
+include("utils/plotresults.jl")
+
+# Id = BEAST.Identity()

examples/mfie_gwp2.jl

@@ -0,0 +1,30 @@
+using CompScienceMeshes, BEAST
+
+# Γ = readmesh(joinpath(dirname(pathof(BEAST)),"../examples/sphere2.in"))
+Γ = meshsphere(radius=1.0, h=0.1)
+# X, Y = raviartthomas(Γ), buffachristiansen(Γ)
+# X = brezzidouglasmarini(Γ)
+# Y = brezzidouglasmarini(Γ)
+# X = raviartthomas(Γ)
+# Y = raviartthomas(Γ)
+X = BEAST.gwpdiv(Γ; order=2)
+Y = BEAST.gwpdiv(Γ; order=2)
+
+ϵ, μ, ω = 1.0, 1.0, 1.0; κ = ω * √(ϵ*μ)
+# κ = 3.0
+NK, Id = BEAST.DoubleLayerRotatedMW3D(1.0, im*κ), Identity()
+E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
+H = -1/(im*μ*ω)*curl(E)
+h = n × H
+
+BEAST.@defaultquadstrat (NK,X,X) BEAST.DoubleNumSauterQstrat(7,8,6,6,6,6)
+BEAST.@defaultquadstrat (Id,X,X) BEAST.SingleNumQStrat(7)
+BEAST.@defaultquadstrat (h,X) BEAST.SingleNumQStrat(12)
+
+@hilbertspace j
+@hilbertspace m
+mfie = @discretise (NK-0.5Id)[m,j] == h[m]  j∈X m∈Y
+u = gmres(mfie)
+
+include("utils/postproc.jl")
+include("utils/plotresults.jl")

src/BEAST.jl

@@ -156,13 +156,17 @@ include("bases/local/ndlcdlocal.jl")
 include("bases/local/bdm3dlocal.jl")
 include("bases/local/rtqlocal.jl")
 include("bases/local/gwplocal.jl")
+include("bases/local/gwpdivlocal.jl")
+
 
 include("bases/basis.jl")
 include("bases/lincomb.jl")
 include("bases/trace.jl")
 include("bases/restrict.jl")
 include("bases/divergence.jl")
 include("bases/global/globaldofs.jl")
+include("bases/global/gwpglobal.jl")
+include("bases/global/gwpdivglobal.jl")
 
 include("bases/lagrange.jl")
 include("bases/rtspace.jl")

src/bases/global/gwpdivglobal.jl

@@ -0,0 +1,39 @@
+struct GWPDivSpace{T,M,P} <: Space{T}
+    geo::M
+    fns::Vector{Vector{Shape{T}}}
+    pos::P
+    degree::Int
+end
+
+function refspace(s::GWPDivSpace{T}) where {T} GWPDivRefSpace{T,s.degree}() end
+function subset(s::S,I) where {S<:GWPDivSpace} S(s.geo, s.fns[I], s.pos[I], s.degree) end
+
+function gwpdiv(mesh, edges=nothing; order)
+
+    T = coordtype(mesh)
+    S = Shape{T}
+
+    space = BEAST.gwpcurl(mesh, edges; order)
+    fns = Vector{Vector{S}}(undef, length(space.fns))
+    for (i,fn) in enumerate(space.fns)
+        fns[i] = [S(s.cellid, s.refid, -s.coeff) for s in fn]
+    end
+
+    return GWPDivSpace(space.geo, fns, space.pos, order)
+end
+
+@testitem "GWPcurl global: numfunctions" begin
+    using CompScienceMeshes
+
+    h = 0.5
+    mesh = meshrectangle(1.0, 1.0, 0.5)
+    edges = setminus(skeleton(mesh,1), boundary(mesh))
+
+    order = 2
+    gwp = BEAST.gwpdiv(mesh, edges; order=order)
+
+    ne = order+1
+    nf = order * (order+1)
+    Nt = length(edges)*ne + length(mesh)*nf
+    @test numfunctions(gwp) == Nt
+end

src/bases/global/gwpglobal.jl

@@ -0,0 +1,179 @@
+struct GWPCurlSpace{T,M,P} <: Space{T}
+    geo::M
+    fns::Vector{Vector{Shape{T}}}
+    pos::P
+    degree::Int
+end
+
+function refspace(s::GWPCurlSpace{T}) where {T} GWPCurlRefSpace{T,s.degree}() end
+function subset(s::S,I) where {S<:GWPCurlSpace} S(s.geo, s.fns[I], s.pos[I], s.degree) end
+
+struct GWPGlobalEdgeDoFs order::Int end
+struct GWPGlobalFaceDoFs order::Int end
+
+function globaldofs(edge_ch, face_ch, localspace, dof::GWPGlobalEdgeDoFs)
+    f = trace(edge_ch, face_ch, localspace)
+    T = coordtype(edge_ch)
+    ds = one(T) / (dof.order+2)
+    return stack(range(ds, step=ds, length=dof.order+1)) do s
+        p = neighborhood(edge_ch, s)
+        t = tangents(p, 1)
+        [dot(t,x.value) for x in f(s)]
+    end
+end
+
+@testitem "GWP edge global dofs" 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),
+    )
+
+    degree = 3
+    localspace = BEAST.GWPCurlRefSpace{T,degree}()
+    dofs = BEAST.GWPGlobalEdgeDoFs(degree)
+    A = BEAST.globaldofs(edge_ch, face_ch, localspace, dofs)
+    display(round.(A, digits=3))
+end
+
+function globaldofs(edge_ch, face_ch, localspace, dof::GWPGlobalFaceDoFs)
+
+    d = dof.order
+    T = coordtype(edge_ch)
+    ds = one(T) / (d+2)
+
+    f = trace(edge_ch, face_ch, localspace)
+    Q = zeros(T, numfunctions(localspace, domain(face_ch)), d*(d+1))
+    S = ((i,j,d+2-i-j) for i in 1:d+1 for j in 1:d+1 if d+2-i-j > 0)
+    # return stack(S) do (i,j,k)
+    idx = 1
+    for (i,j,k) in S
+        u = (i*ds,j*ds)
+        p = neighborhood(edge_ch, u)
+        t1 = tangents(p,1)
+        t2 = tangents(p,2)
+        vals = f(u)
+        q1 = [dot(t1,x.value) for x in vals]
+        q2 = [dot(t2,x.value) for x in vals]
+        Q[:,idx] = q1; idx += 1
+        Q[:,idx] = q2; idx += 1
+    end
+
+    return Q
+end
+
+
+@testitem "GWP face global dofs" 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),
+    #     point(3,0,0)
+    # )
+    edge_ch = simplex(
+        point(3,0,0),
+        point(2,0,-1),
+        point(0,0,1),
+    )
+
+    degree = 3
+    localspace = BEAST.GWPCurlRefSpace{T,degree}()
+    dofs = BEAST.GWPGlobalFaceDoFs(degree)
+    A = BEAST.globaldofs(edge_ch, face_ch, localspace, dofs)
+    display(round.(A, digits=3))
+end
+
+function _addshapes!(fns, cell, gids, lids, β)
+    T = eltype(β)
+    atol = sqrt(eps(T))
+    for i in axes(β,1)
+        fn = fns[gids[i]]
+        for j in axes(β,2)
+            isapprox(β[i,j], 0; atol) && continue
+            push!(fn, Shape{T}(cell, lids[j], β[i,j]))
+        end
+    end
+end
+
+
+function gwpcurl(mesh, edges=nothing; order)
+
+    T = coordtype(mesh)
+    atol = sqrt(eps(T))
+
+    if edges == nothing
+        edges = setminus(skeleton(mesh, 1), boundary(mesh))
+    end
+
+    C12 = connectivity(mesh, edges, abs)
+    R12 = rowvals(C12)
+    V12 = nonzeros(C12)
+
+    ne = order+1
+    nf = order*(order+1)
+
+    Ne = ne * length(edges)
+    Nf = nf * length(mesh)
+    Nt = Ne + Nf
+
+    localspace = GWPCurlRefSpace{T,order}()
+    refdom = domain(chart(mesh, first(mesh)))
+    localdim = numfunctions(localspace, refdom)
+
+    fns = [Shape{T}[] for n in 1:Nt]
+    pos = fill(point(0,0,0), Nt)
+    for cell in mesh
+        cell_ch = chart(mesh, cell)
+        for k in nzrange(C12, cell)
+            e = R12[k]
+            i = V12[k]
+            gids = (e-1)*ne .+ (1:ne)
+            lids = localindices(localspace, refdom, Val{1}, i)
+            edge_ch = chart(edges, e)
+            vals = globaldofs(edge_ch, cell_ch, localspace, GWPGlobalEdgeDoFs(order))
+            α = vals[lids,:]
+            _addshapes!(fns, cell, gids, lids, inv(α'))
+        end
+
+        order < 2 && continue
+        face_ch = chart(mesh, cell)
+        gids = Ne + (cell-1)*nf .+ (1:nf)
+        lids = localindices(localspace, refdom, Val{2}, 1)
+        vals = globaldofs(face_ch, face_ch, localspace, GWPGlobalFaceDoFs(order))
+        α = vals[lids,:]
+        _addshapes!(fns, cell, gids, lids, inv(α'))
+    end
+
+    for e in edges pos[(e-1)*ne .+ (1:ne)] .= Ref(cartesian(center(chart(edges, e)))) end
+    for f in mesh pos[Ne+(f-1)*nf .+ (1: nf)] .= Ref(cartesian(center(chart(mesh, f)))) end
+
+    return GWPCurlSpace(mesh, fns, pos, order)
+end
+
+@testitem "GWPcurl global: numfunctions" begin
+    using CompScienceMeshes
+
+    h = 0.5
+    mesh = meshrectangle(1.0, 1.0, 0.5)
+    edges = setminus(skeleton(mesh,1), boundary(mesh))
+
+    order = 2
+    gwp = BEAST.gwpcurl(mesh, edges; order=order)
+
+    ne = order+1
+    nf = order * (order+1)
+    Nt = length(edges)*ne + length(mesh)*nf
+    @test numfunctions(gwp) == Nt
+end

src/bases/lagrange.jl

@@ -977,7 +977,8 @@ function lagrangec0(mesh::CompScienceMeshes.AbstractMesh{<:Any,3}; order)
     N = nV + nE + nF
 
     localspace = LagrangeRefSpace{T,order,3,binomial(2+order,2)}()
-    localdim = numfunctions(localspace)
+    dom = domain(chart(mesh, first(mesh)))
+    localdim = numfunctions(localspace, dom)
 
     d = order
     fns = [S[] for n in 1:(nV+nE+nF)]

src/bases/local/gwpdivlocal.jl

@@ -0,0 +1,69 @@
+struct GWPDivRefSpace{T,Degree} <: RefSpace{T} end
+
+function numfunctions(x::GWPDivRefSpace{<:Any,D},
+    dom::CompScienceMeshes.ReferenceSimplex{2}) where{D} (D+1)*(D+3) end
+
+function (ϕ::GWPDivRefSpace{T,Degree})(p) where {T,Degree}
+    ψ = GWPCurlRefSpace{T,Degree}()
+    dom = domain(chart(p))
+    u = parametric(p)
+    vals = ψ(dom, u)
+    vals = pushforwardcurl(vals, p)
+    n = normal(p)
+    map(vals) do v
+        (value = -n × v.value, divergence=v.curl)
+    end
+end
+
+function interpolate(fields, interpolant::GWPDivRefSpace{T,Degree}, chart) where {T,Degree}
+
+    function nxfields(p)
+        vals = fields(p)
+        n = normal(p)
+        return map(v -> n×v, vals)
+    end
+
+    return interpolate(nxfields, GWPCurlRefSpace{T,Degree}(), chart)
+end
+
+@testitem "divergence" begin
+    using CompScienceMeshes, LinearAlgebra
+
+    T = Float64
+    s = simplex(
+        point(3,0,0),
+        point(0,2,1),
+        point(-1,-1,-1),
+    )
+
+    order = 2
+    ϕ = BEAST.GWPDivRefSpace{T,order}()
+
+    u = T(0.2432); dx = sqrt(eps(T))
+    v = T(0.5786); dy = sqrt(eps(T))
+
+    # p00 = neighborhood(s, (u,v))
+    # p10 = neighborhood(s, (u+du,v))
+    # p01 = neighborhood(s, (u, v+dv))
+
+    p00 = neighborhood(s, (u,v))
+    t1 = normalize(tangents(p00,1))
+    t2 = normal(p00) × t1
+
+    p10 = neighborhood(s, carttobary(s, cartesian(p00) + dx*t1 + 0*t2))
+    p01 = neighborhood(s, carttobary(s, cartesian(p00) + 0*t1 + dy*t2))
+
+    ϕ00 = ϕ(p00)
+    ϕ10 = ϕ(p10)
+    ϕ01 = ϕ(p01)
+
+
+    for (f00, f10, f01) in zip(ϕ00, ϕ10, ϕ01)
+        df1dx = dot(t1, f10.value - f00.value) / dx
+        df2dy = dot(t2, f01.value - f00.value) / dy
+        curl_num = df1dx + df2dy
+        curl_ana = f00.divergence
+        # @show curl_num,  curl_ana, abs(curl_num - curl_ana)
+        @test curl_num ≈ curl_ana atol=sqrt(eps(T))*100
+    end
+end