numfunctions now requires both the refspace and refdomain

src/BEAST.jl

@@ -141,14 +141,9 @@ include("utils/combinatorics.jl")
 include("utils/linearspace.jl")
 include("utils/zeromap.jl")
 include("utils/rank1map.jl")
+include("utils/lagpolys.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/localbasis.jl")
 include("bases/local/laglocal.jl")
 include("bases/local/rtlocal.jl")
 include("bases/local/rt2local.jl")
@@ -160,6 +155,14 @@ include("bases/local/ndlcclocal.jl")
 include("bases/local/ndlcdlocal.jl")
 include("bases/local/bdm3dlocal.jl")
 include("bases/local/rtqlocal.jl")
+include("bases/local/gwplocal.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/lagrange.jl")
 include("bases/rtspace.jl")

src/bases/basis.jl

@@ -1,5 +1,4 @@
-abstract type RefSpace{T,D} end
-abstract type DivRefSpace{T,D} <: RefSpace{T,D} end
+
 
 abstract type AbstractSpace end
 abstract type Space{T} <: AbstractSpace end
@@ -25,12 +24,12 @@ Returns the ReferenceSpace of local shape functions on which the basis is built.
 function refspace end
 
 
-"""
-    numfunctions(r)
+# """
+#     numfunctions(r)
 
-Return the number of functions in a `Space` or `RefSpace`.
-"""
-numfunctions(rs::RefSpace{T,D}) where {T,D} = D
+# Return the number of functions in a `Space` or `RefSpace`.
+# """
+# numfunctions(rs::RefSpace{T,D}) where {T,D} = D
 
 
 """
@@ -208,7 +207,10 @@ function assemblydata(basis::Space; onlyactives=true)
     num_cells = numcells(geo)
 
     num_bfs  = numfunctions(basis)
-    num_refs = numfunctions(refspace(basis))
+
+    ch = chart(geo, first(geo))
+    dom = domain(ch)
+    num_refs = numfunctions(refspace(basis), dom)
 
     # Determine the number of functions defined over a given cell
     # and per local shape function.

src/bases/global/globaldofs.jl

@@ -20,16 +20,7 @@ function trace(edge_ch, face_ch, localspace::RefSpace)
     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
@@ -38,7 +29,7 @@ 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]
+    [_lagpoly(nodes, i, s) for i in 2:order]
 end
 
 struct _LagrangeGlobalFaceDoFs
@@ -54,12 +45,12 @@ end
 #     idx = 1
 #     degree = dof.degree
 #     for i in 0:degree
-#         prodi = _lagrangepolynomial(nodes, i+1, s1, i)
+#         prodi = _lagpoly(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)
+#             prodj = _lagpoly(nodes, j+1, s2, j)
+#             prodk = _lagpoly(nodes, k+1, s3, k)
 #             r[idx] = prodi * prodj * prodk
 #     end end
 #     return r

src/bases/local/bdm3dlocal.jl

@@ -1,4 +1,4 @@
-struct BDM3DRefSpace{T} <: RefSpace{T,12} end
+struct BDM3DRefSpace{T} <: RefSpace{T} end
 
 function (f::BDM3DRefSpace)(p)
 

src/bases/local/bdmlocal.jl

@@ -1,4 +1,4 @@
-struct BDMRefSpace{T} <: RefSpace{T,6} end
+struct BDMRefSpace{T} <: RefSpace{T} end
 
 function (f::BDMRefSpace)(p)
 

src/bases/local/gwplocal.jl

@@ -0,0 +1,76 @@
+struct GWPCurlRefSpace{T,Degree} <: RefSpace{T} end
+
+function (ϕ::GWPCurlRefSpace{T,Degree})(p) where {T,Degree}
+    dom = domain(chart(p))
+    u = parametric(p)
+    vals = ϕ(dom, u)
+    pushforwardcurl(vals, p)
+end
+
+function (ϕ::GWPCurlRefSpace{T,Deg})(dom::CompScienceMeshes.ReferenceSimplex{Dim},
+    u) where {T,Deg,Dim}
+
+    d = Deg
+    u, v = u
+    w = 1-u-v
+
+    s = range(zero(T), one(T), length=d+3)    
+    # rwg1 = point(T, u-1, v)
+    # rwg2 = point(T, u, v-1)
+    # rwg3 = point(T, u, v)
+    nd1 = point(T, -v, u-1)
+    nd2 = point(T, -v+1, u)
+    nd3 = point(T, -v, u)
+
+    P = SVector{2,T}
+    vals = P[]
+    dffu = T[]
+    dffv = T[]
+
+    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)
+            push!(vals, Rᵢ*Rⱼ*Rₖ*nd1)
+
+            # du = _sylpoly_shift_diff(s, i+1, u)*Rⱼ*Rₖ
+    end
+
+    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)
+        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)
+        push!(vals, Rᵢ*Rⱼ*Rₖ*nd3)
+    end
+
+    for i in 1:d+1
+        for j in 1:d+1
+            k = (d+2)-i-j
+            k <= 0 && continue
+            Rsᵢ = _sylpoly_shift(s, i+1, u)
+            Rsⱼ = _sylpoly_shift(s, j+1, v)
+            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)
+    end end
+
+    NF = length(vals)
+    SVector{NF}([(value=f, curl=zero(T)) for f in vals])
+end
+
+

src/bases/local/laglocal.jl

@@ -1,13 +1,13 @@
 # T: coeff type
 # Degree: degree
 # Dim1: dimension of the support + 1
-struct LagrangeRefSpace{T,Degree,Dim1,NF} <: RefSpace{T,NF} end
+struct LagrangeRefSpace{T,Degree,Dim1,NF} <: RefSpace{T} end
 
-numfunctions(s::LagrangeRefSpace{T,D,2}) where {T,D} = D+1
-numfunctions(s::LagrangeRefSpace{T,0,3}) where {T} = 1
-numfunctions(s::LagrangeRefSpace{T,1,3}) where {T} = 3
-numfunctions(s::LagrangeRefSpace{T,2,3}) where {T} = 6
-numfunctions(s::LagrangeRefSpace{T,Dg,D1,NF}) where {T,Dg,D1,NF} = NF
+numfunctions(s::LagrangeRefSpace{T,D,2}, ch::CompScienceMeshes.ReferenceSimplex{1}) where {T,D} = D+1
+numfunctions(s::LagrangeRefSpace{T,0,3}, ch::CompScienceMeshes.ReferenceSimplex{2}) where {T} = 1
+numfunctions(s::LagrangeRefSpace{T,1,3}, ch::CompScienceMeshes.ReferenceSimplex{2}) where {T} = 3
+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)
 
 valuetype(ref::LagrangeRefSpace{T}, charttype) where {T} =
         SVector{numfunctions(ref), Tuple{T,T}}
@@ -110,7 +110,8 @@ function strace(x::LagrangeRefSpace, cell, localid, face)
     vals1 = x(P1)
     vals2 = x(P2)
 
-    for j in 1:numfunctions(x)
+    num_shapes = numfunctions(x, domain(cell))
+    for j in 1:num_shapes
         Q[1,j] = vals1[j].value
         Q[2,j] = vals2[j].value
     end
@@ -137,12 +138,13 @@ end
 
 
 function restrict(refs::LagrangeRefSpace{T,0}, dom1, dom2) where T
-    Q = Matrix{T}(I, numfunctions(refs), numfunctions(refs))
+    n = numfunctions(refs, domain(dom1))
+    Q = Matrix{T}(I, n, n)
 end
 
 function restrict(f::LagrangeRefSpace{T,1}, dom1, dom2) where T
 
-    D = numfunctions(f)
+    D = numfunctions(f, domain(dom1))
     Q = zeros(T, D, D)
 
     # for each point of the new domain

src/bases/local/ncrossbdmlocal.jl

@@ -1,4 +1,4 @@
-struct NCrossBDMRefSpace{T} <: RefSpace{T,6} end
+struct NCrossBDMRefSpace{T} <: RefSpace{T} end
 
 function (f::NCrossBDMRefSpace{T})(p) where T
 

src/bases/local/nd2local.jl

@@ -1,5 +1,5 @@
 
-mutable struct ND2RefSpace{T} <: RefSpace{T,8} end
+mutable struct ND2RefSpace{T} <: RefSpace{T} end
 
 function (ϕ::ND2RefSpace)(nbd)
 
@@ -26,3 +26,5 @@ function (ϕ::ND2RefSpace)(nbd)
     ))
 
 end
+
+numfunctions(x::ND2RefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 8

src/bases/local/ndlcclocal.jl

@@ -1,4 +1,4 @@
-struct NDLCCRefSpace{T} <: RefSpace{T,6} end
+struct NDLCCRefSpace{T} <: RefSpace{T} end
 
 function valuetype(ref::NDLCCRefSpace{T}, charttype::Type) where {T}
     SVector{universedimension(charttype),T}
@@ -58,6 +58,8 @@ function (ϕ::NDLCCRefSpace)(ndlc)
 ##        ))
 end
 
+numfunctions(x::NDLCCRefSpace, dom::CompScienceMeshes.ReferenceSimplex{3}) = 6
+
 #check orientation
 function curl(ref::NDLCCRefSpace, sh, el)
     a = [4,2,3,4,1,2]##[2,1,1,1,2,4]#[2,1,4,4,3,2]#[4,2,3,4,1,2]
@@ -71,7 +73,7 @@ end
 function restrict(ϕ::NDLCCRefSpace{T}, dom1, dom2) where {T}
     # dom2 is the smaller of the domains
 
-    K = numfunctions(ϕ)
+    K = numfunctions(ϕ, domain(dom1))
     D = dimension(dom1)
 
     @assert K == 6

src/bases/local/ndlcdlocal.jl

@@ -1,4 +1,4 @@
-struct NDLCDRefSpace{T} <: RefSpace{T,4} end
+struct NDLCDRefSpace{T} <: RefSpace{T} end
 
 function valuetype(ref::NDLCDRefSpace{T}, charttype::Type) where {T}
     SVector{universedimension(charttype),T}
@@ -30,6 +30,8 @@ function (ϕ::NDLCDRefSpace)(ndlc)
     ))
 end
 
+numfunctions(x::NDLCDRefSpace, dom::CompScienceMeshes.ReferenceSimplex{3}) = 4
+
 function ntrace(x::NDLCDRefSpace, el, q, fc)
     t = zeros(scalartype(x),1,4)
     t[q] = 1 / volume(fc)
@@ -82,7 +84,7 @@ divergence(ref::NDLCDRefSpace, sh, el) = Shape(sh.cellid, 1, sh.coeff/volume(el)
 function restrict(ϕ::NDLCDRefSpace{T}, dom1, dom2) where {T}
     # dom2 is the smaller of the domains
 
-    K = numfunctions(ϕ)
+    K = numfunctions(ϕ, domain(dom1))
     D = dimension(dom1)
 
     @assert K == 4

src/bases/local/ndlocal.jl

@@ -6,7 +6,7 @@ This is not the edge starting at vertex `r`. The downside of this local numberin
 scheme is that it cannot be extended to cells that are not simplices because
 there is no well defined concept of adjacent-ness.
 """
-mutable struct NDRefSpace{T} <: RefSpace{T,3} end
+mutable struct NDRefSpace{T} <: RefSpace{T} end
 
 function (ϕ::NDRefSpace)(nbd)
 
@@ -27,10 +27,11 @@ function (ϕ::NDRefSpace)(nbd)
 
 end
 
+numfunctions(x::NDRefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 3
 
 function restrict(ϕ::NDRefSpace{T}, dom1, dom2) where T
 
-    K = numfunctions(ϕ)
+    K = numfunctions(ϕ, domain(dom1))
     D = dimension(dom1)
 
     @assert K == 3

src/bases/local/rt2local.jl

@@ -1,4 +1,4 @@
-struct RT2RefSpace{T} <: RefSpace{T,8} end
+struct RT2RefSpace{T} <: RefSpace{T} end
 
 function (f::RT2RefSpace)(p)
 
@@ -33,6 +33,7 @@ function (f::RT2RefSpace)(p)
     )
 end
 
+numfunctions(x::RT2RefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 8
 
 function interpolate(fields, interpolant::BEAST.RT2RefSpace, chart)
 

src/bases/local/rtlocal.jl

@@ -1,4 +1,4 @@
-struct RTRefSpace{T} <: DivRefSpace{T,3} end
+struct RTRefSpace{T} <: DivRefSpace{T} end
 
 # valuetype(ref::RTRefSpace{T}, charttype) where {T} = SVector{3,Tuple{SVector{universedimension(charttype),T},T}}
 function valuetype(ref::RTRefSpace{T}, charttype::Type) where {T}
@@ -26,6 +26,8 @@ function (ϕ::RTRefSpace)(mp)
     ))
 end
 
+numfunctions(x::RTRefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 3
+
 divergence(ref::RTRefSpace, sh, el) = Shape(sh.cellid, 1, sh.coeff/volume(el))
 
 """

src/bases/local/rtqlocal.jl

@@ -1,4 +1,4 @@
-struct RTQuadRefSpace{T} <: DivRefSpace{T,4} end
+struct RTQuadRefSpace{T} <: DivRefSpace{T} end
 
 function (ϕ::RTQuadRefSpace{T})(p) where {T}