self-interpolation verified

src/bases/global/globaldofs.jl

@@ -1,4 +1,4 @@
-function trace(edge_ch, face_ch, localspace::RefSpace)
+function trace(edge_ch, face_ch, fields)
     T = coordtype(edge_ch)
     atol = sqrt(eps(T))
 
@@ -14,7 +14,7 @@ function trace(edge_ch, face_ch, localspace::RefSpace)
     function f(s)
         u = neighborhood(injection, s)
         p = neighborhood(face_ch, cartesian(u))
-        localspace(p)
+        fields(p)
     end
 
     return f

src/bases/local/gwplocal.jl

@@ -30,9 +30,9 @@ function (ϕ::GWPCurlRefSpace{T,Deg})(dom::CompScienceMeshes.ReferenceSimplex{Di
     i = 0
     for j in 1:d+1
         k = (d+2)-i-j
-            Rᵢ = _sylpoly_shift(s, i+1, u)
-            Rⱼ = _sylpoly(s, j+1, v)
-            Rₖ = _sylpoly(s, k+1, w)
+            Rᵢ = _sylpoly(s, i+1, u)
+            Rⱼ = _sylpoly_shift(s, j+1, v)
+            Rₖ = _sylpoly_shift(s, k+1, w)
             push!(vals, Rᵢ*Rⱼ*Rₖ*nd1)
 
             # du = _sylpoly_shift_diff(s, i+1, u)*Rⱼ*Rₖ
@@ -41,18 +41,18 @@ function (ϕ::GWPCurlRefSpace{T,Deg})(dom::CompScienceMeshes.ReferenceSimplex{Di
     for i in 1:d+1
         j = 0
         k = (d+2)-i-j
-        Rᵢ = _sylpoly(s, i+1, u)
-        Rⱼ = _sylpoly_shift(s, j+1, v)
-        Rₖ = _sylpoly(s, k+1, w)
+        Rᵢ = _sylpoly_shift(s, i+1, u)
+        Rⱼ = _sylpoly(s, j+1, v)
+        Rₖ = _sylpoly_shift(s, k+1, w)
         push!(vals, Rᵢ*Rⱼ*Rₖ*nd2)
     end
 
     for i in 1:d+1
         j = (d+2)-i
         k = 0
-        Rᵢ = _sylpoly(s, i+1, u)
-        Rⱼ = _sylpoly(s, j+1, v)
-        Rₖ = _sylpoly_shift(s, k+1, w)
+        Rᵢ = _sylpoly_shift(s, i+1, u)
+        Rⱼ = _sylpoly_shift(s, j+1, v)
+        Rₖ = _sylpoly(s, k+1, w)
         push!(vals, Rᵢ*Rⱼ*Rₖ*nd3)
     end
 
@@ -62,15 +62,83 @@ function (ϕ::GWPCurlRefSpace{T,Deg})(dom::CompScienceMeshes.ReferenceSimplex{Di
             k <= 0 && continue
             Rsᵢ = _sylpoly_shift(s, i+1, u)
             Rsⱼ = _sylpoly_shift(s, j+1, v)
+            Rsₖ = _sylpoly_shift(s, k+1, w)
             Rᵢ = _sylpoly(s, i+1, u)
             Rⱼ = _sylpoly(s, j+1, v)
             Rₖ = _sylpoly(s, k+1, w)
-            push!(vals, Rsᵢ*Rⱼ*Rₖ*nd1)
-            push!(vals, Rᵢ*Rsⱼ*Rₖ*nd2)
+            S1 = Rᵢ*Rsⱼ*Rsₖ*nd1
+            S2 = Rsᵢ*Rⱼ*Rsₖ*nd2
+            S3 = Rsᵢ*Rsⱼ*Rₖ*nd3
+            N1 = (d+2)/(d+2-i)
+            N2 = (d+2)/(d+2-j)
+            N3 = (d+2)/(d+2-k)
+            push!(vals, S2-S3)
+            push!(vals, S3-S1)
     end end
 
     NF = length(vals)
     SVector{NF}([(value=f, curl=zero(T)) for f in vals])
 end
 
 
+function interpolate(fields, interpolant::GWPCurlRefSpace{T,Degree}, chart) where {T,Degree}
+
+    d = Degree
+    dim = (d+1)*(d+3)
+
+    s = range(zero(T), one(T), length=d+3)
+
+    edges = faces(chart)
+
+    edge = edges[1]
+    fields_edge = trace(edge, chart, fields)
+    i = 0
+    Q1 = stack(1:d+1) do j
+        k = (d+2)-i-j
+        u_edge = s[j+1]
+        p_edge = neighborhood(edge, (u_edge,))
+        @show cartesian(p_edge)
+        t_edge = -tangents(p_edge, 1)
+        vals = fields_edge(u_edge)
+        [dot(t_edge, val) for val in vals]
+    end
+
+    edge = edges[2]
+    fields_edge = trace(edge, chart, fields)
+    Q2 = stack(1:d+1) do i
+        j = 0
+        k = (d+2)-i-j
+        u_edge = 1 - s[i+1]
+        p_edge = neighborhood(edge, (u_edge,))
+        t_edge = -tangents(p_edge, 1)
+        vals = fields_edge(u_edge)
+        [dot(t_edge, val) for val in vals]
+    end
+
+    edge = edges[3]
+    fields_edge = trace(edge, chart, fields)
+    Q3 = stack(1:d+1) do i
+        j = (d+2)-i
+        k = 0
+        u_edge = 1-s[j+1]
+        p_edge = neighborhood(edge, (u_edge,))
+        t_edge = -tangents(p_edge, 1)
+        vals = fields_edge(u_edge)
+        [dot(t_edge, val) for val in vals]
+    end
+
+    Q = hcat(Q1,Q2,Q3)
+    if 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)
+        for (i,j,k) in S
+            p_chart = neighborhood(chart, (s[i+1],s[j+1]))
+            t_i = tangents(p_chart, 1)
+            t_j = tangents(p_chart, 2)
+            vals = fields(p_chart)
+            q_i = [dot(t_i, val) for val in vals]
+            q_j = [dot(t_j, val) for val in vals]
+            Q = hcat(Q, q_i, q_j)
+        end
+    end
+    return Q
+end

test/test_gwp.jl

@@ -1,4 +1,4 @@
-@testitem "GWPRefSpace eval" begin
+@testitem "refspace: dimension" begin
     using CompScienceMeshes
 
     T = Float64
@@ -9,4 +9,33 @@
 
     u = (one(T)/3, one(T)/3)
     vals = ϕ(dom, u)
+
+    nf = (Degree+1)*(Degree+3)
+    @test length(vals) == nf
+end
+
+
+@testitem "refspace: self-interpolate" begin
+    using CompScienceMeshes
+
+    T, Dim, Nverts, Degree = Float64, 2, 3, 2
+    dom =  CompScienceMeshes.ReferenceSimplex{Dim,T,Nverts}()
+    ϕ = BEAST.GWPCurlRefSpace{T,Degree}()
+
+    fields(p) = [v.value for v in ϕ(dom, p)]
+    # fields(p) = [ϕ(dom,p)[1].value]
+    coeffs = BEAST.interpolate(fields, ϕ, dom)
+
+    display(round.(coeffs, digits=3))
+
+    # pts = [
+    #     point(T, 0.3, 0.1),
+    #     point(T, 0.1, 0.3),
+    #     point(T, 0.5, 0.0)]
+
+    # for u in pts
+    #     vals = ϕ(dom, u)
+    #     r = sum(c * v.value for (c,v) in zip(coeffs, vals))
+    #     @show r
+    # end
 end