support to take the divergence of GWP space

examples/divgwp_expansion.jl

@@ -0,0 +1,66 @@
+using BEAST
+using CompScienceMeshes
+const CSM = CompScienceMeshes
+using  StaticArrays
+
+T = Float64
+D = 4
+NF = binomial(2+D, 2)
+gwp = BEAST.GWPDivRefSpace{T,D}()
+lgx = BEAST.LagrangeRefSpace{T,D,3,10}()
+
+function fields(p)
+    map(gwp(p)) do x
+        x.divergence
+    end
+end
+
+ch = CSM.simplex(
+    point(1,0,0),
+    point(0,1,0),
+    point(0,0,0))
+
+# p = CSM.center(ch)
+# v = fields(p)
+
+coeffs = BEAST.interpolate(fields, lgx, ch)
+Q = rationalize.(coeffs; tol=sqrt(eps(T)))
+
+using LinearAlgebra
+norm(Q - coeffs)
+
+Q3 = Rational{Int64}[89//3 223//81 -7//81 -2 223//81 -7//81 -2 -7//81 -2 -2; -59//2 236//27 -13//54 13 -83//9 -4//27 8 -1//18 3 -2; 13 -13//54 236//27 -59//2 8 -4//27 -83//9 3 -1//18 -2; -2 -7//81 223//81 89//3 -2 -7//81 223//81 -2 -7//81 -2; 89//3 223//81 -7//81 -2 223//81 -7//81 -2 -7//81 -2 -2; -59//2 -83//9 -1//18 -2 236//27 -4//27 3 -13//54 8 13; 13 8 3 -2 -13//54 -4//27 -1//18 236//27 -83//9 -59//2; -2 -2 -2 -2 -7//81 -7//81 -7//81 223//81 223//81 89//3; -2 -7//81 223//81 89//3 -2 -7//81 223//81 -2 -7//81 -2; -2 -1//18 -83//9 -59//2 3 -4//27 236//27 8 -13//54 13; -2 3 8 13 -1//18 -4//27 -13//54 -83//9 236//27 -59//2; -2 -2 -2 -2 -7//81 -7//81 -7//81 223//81 223//81 89//3; -50 890//81 55//81 35//3 -565//81 125//162 20//3 70//81 5//3 -10//3; 50 565//81 -70//81 10//3 -890//81 -125//162 -5//3 -55//81 -20//3 -35//3; 30 -355//27 355//27 -30 80//9 0 -80//9 -10//9 10//9 0; -40 440//27 85//27 -10 5//2 -125//27 205//18 35//3 95//9 -25//2; -35//3 -55//81 -890//81 50 -20//3 -125//162 565//81 -5//3 -70//81 10//3; 15 -5//27 485//27 0 5 0 -485//27 -5 5//27 -15; 40 -5//2 -35//3 25//2 -440//27 125//27 -95//9 -85//27 -205//18 10; -30 -80//9 10//9 0 355//27 0 -10//9 -355//27 80//9 30; -25//2 35//3 5//2 -40 95//9 -125//27 440//27 205//18 85//27 -10; 25//2 -95//9 -205//18 10 -35//3 125//27 -85//27 -5//2 -440//27 40; -15 -5 5 15 5//27 0 -5//27 -485//27 485//27 0; 35//3 20//3 5//3 -10//3 55//81 125//162 70//81 890//81 -565//81 -50]
+function BEAST.divergence(localspace::BEAST.GWPDivRefSpace, sh, ch)
+    BEAST.divergence(localspace, sh, ch, BEAST.dimtype(localspace, domain(ch)))
+end
+
+function BEAST.divergence(localspace::BEAST.GWPDivRefSpace{T,D}, cellid, ch, ::Type{Val{N}}) where {N,D}
+    function fields(p)
+        map(localspace(p)) do x
+            x.divergence
+        end
+    end
+    T = coordtype(ch)
+    Dim = 2
+    NFout = div((Dim+1)*(Dim+2), 2)
+    lag = BEAST.LagrangeRefSpace{T,D,Dim+1,NFout}()
+    coeffs = BEAST.interpolate(fields, lag, ch)
+    S = BEAST.Shape{T}
+    A = Vector{Vector{S}}(undef, size(coeffs,1))
+    for r in axes(coeffs,1)
+        A[r] = collect(S(cellid, c, coeffs[r,c]) for c in axes(coeffs,2))
+    end
+    return SVector{N}(A)
+end
+
+divfns = BEAST.divergence(gwp, 1, ch)
+
+p = neighborhood(ch, (0.2, 0.6))
+gwp_vals = gwp(p)
+val1 = gwp_vals[1].divergence
+
+lgx_vals = lgx(p)
+val2 = zero(T)
+for sh in divfns[1]
+    val2 += sh.coeff * lgx_vals[sh.refid].value
+end

src/bases/divergence.jl

@@ -17,10 +17,11 @@ function divergence(x::Space)
     fns = x.fns
     dvs = similar(fns)
     for (i,fn) in enumerate(fns)
-        dvs[i] = similar(fns[i])
+        dvs[i] = similar(fns[i], 0)
         for (j,sh) in enumerate(fn)
             el = els[sh.cellid]
-            dvs[i][j] = divergence(ref, sh, el)
+            append!(dvs[i], divergence(ref, sh, el))
+            # dvs[i][j] = divergence(ref, sh, el)
         end
     end
     divergence(x, geo, dvs)

src/bases/global/gwpdivglobal.jl

@@ -36,4 +36,62 @@ end
     nf = order * (order+1)
     Nt = length(edges)*ne + length(mesh)*nf
     @test numfunctions(gwp) == Nt
+end
+
+
+function divergence(X::GWPDivSpace, geo, fns)
+    ch = chart(geo, first(geo))
+    dim = dimension(ch)
+    dom = domain(ch)
+    # dim = dimension(dom)
+
+    P = X.degree
+    Cont = -1
+    NF = binomial(dim+P, dim)
+
+    LagrangeBasis{P,Cont,NF}(geo, fns, deepcopy(positions(X)))
+end
+
+
+@testitem "divergence - global" begin
+    using CompScienceMeshes
+    const CSM = CompScienceMeshes
+
+    T = Float64
+
+    m = CSM.meshrectangle(1.0, 1.0, 0.5, 3)
+    X = BEAST.gwpdiv(m; order=4)
+    divX = BEAST.divergence(X)
+
+    x = BEAST.refspace(X)
+    divx = BEAST.refspace(divX)
+
+    err = zero(T)
+    for i in eachindex(X.fns)
+        fn = X.fns[i]
+        for j in eachindex(fn)
+            cellid = X.fns[i][j].cellid
+            ch = chart(m, cellid)
+
+            u = (0.2341, 0.4312)
+            p = neighborhood(ch, u)
+
+            r1 = zero(T)
+            ϕp = x(p)
+            for sh in X.fns[i]
+                sh.cellid == cellid || continue
+                r1 += sh.coeff * ϕp[sh.refid].divergence
+            end
+
+            r2 = zero(T)
+            ϕp = divx(p)
+            for sh in divX.fns[i]
+                sh.cellid == cellid || continue
+                r2 += sh.coeff * ϕp[sh.refid].value
+            end
+            global err = max(err, abs(r1-r2))
+        end
+    end
+
+    @test err < sqrt(eps(T))
 end

src/bases/lagrange.jl

@@ -11,6 +11,7 @@ function lagdimension end
 # M: mesh type
 # T: field type
 # NF: number of local shape functions
+# P: point type
 mutable struct LagrangeBasis{D,C,M,T,NF,P} <: Space{T}
   geo::M
   fns::Vector{Vector{Shape{T}}}

src/bases/local/bdmlocal.jl

@@ -19,7 +19,7 @@ function (f::BDMRefSpace)(p)
         (value=v*tv/j,         divergence=d),]
 end
 
-divergence(ref::BDMRefSpace, sh, el) = Shape(sh.cellid, 1, sh.coeff/(2*volume(el)))
+divergence(ref::BDMRefSpace, sh, el) = [Shape(sh.cellid, 1, sh.coeff/(2*volume(el)))]
 
 const _vert_perms_bdm = [
     (1,2,3),

src/bases/local/gwpdivlocal.jl

@@ -2,6 +2,8 @@ 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 dimtype(x::GWPDivRefSpace{<:Any,D},
+    dom::CompScienceMeshes.ReferenceSimplex{2}) where{D} Val{(D+1)*(D+3)} end
 
 function (ϕ::GWPDivRefSpace{T,Degree})(p) where {T,Degree}
     ψ = GWPCurlRefSpace{T,Degree}()
@@ -26,7 +28,7 @@ function interpolate(fields, interpolant::GWPDivRefSpace{T,Degree}, chart) where
     return interpolate(nxfields, GWPCurlRefSpace{T,Degree}(), chart)
 end
 
-@testitem "divergence" begin
+@testitem "divergence - pointwise" begin
     using CompScienceMeshes, LinearAlgebra
 
     T = Float64
@@ -66,4 +68,72 @@ end
         # @show curl_num,  curl_ana, abs(curl_num - curl_ana)
         @test curl_num ≈ curl_ana atol=sqrt(eps(T))*100
     end
+end
+
+
+function divergence(localspace::GWPDivRefSpace, sh, ch)
+    fns = divergence(localspace, sh.cellid, ch)
+    α = sh.coeff
+    S = typeof(sh)
+    return S[S(s.cellid, s.refid, α*s.coeff) for s in fns[sh.refid]]
+end
+
+
+function divergence(localspace::GWPDivRefSpace, cellid::Int, ch)
+    divergence(localspace, cellid, ch, dimtype(localspace, domain(ch)))
+end
+
+
+function divergence(localspace::GWPDivRefSpace{T,D}, cellid::Int, ch, ::Type{Val{N}}) where {N,D,T}
+    function fields(p)
+        map(localspace(p)) do x
+            x.divergence
+        end
+    end
+    atol = sqrt(eps(T))
+    Dim = 2
+    NFout = div((Dim+1)*(Dim+2), 2)
+    lag = LagrangeRefSpace{T,D,Dim+1,NFout}()
+    coeffs = interpolate(fields, lag, ch)
+    S = BEAST.Shape{T}
+    A = Vector{Vector{S}}(undef, size(coeffs,1))
+    for r in axes(coeffs,1)
+        A[r] = collect(S(cellid, c, coeffs[r,c]) for c in axes(coeffs,2) if abs(c) > atol)
+    end
+    return SVector{N}(A)
+end
+
+
+@testitem "divergence - chartwise" begin
+    using CompScienceMeshes
+    const CSM = CompScienceMeshes
+
+    T = Float64
+    D = 4
+    NF = binomial(2+D, 2)
+    gwp = BEAST.GWPDivRefSpace{T,D}()
+    lgx = BEAST.LagrangeRefSpace{T,D,3,10}()
+
+    ch = CSM.simplex(
+        point(1,0,0),
+        point(0,1,0),
+        point(0,0,0))
+
+    divfns = BEAST.divergence(gwp, 1, ch)
+
+    p = neighborhood(ch, (0.2, 0.6))
+    gwp_vals = gwp(p)
+    lgx_vals = lgx(p)
+
+    err = similar(Vector{T}, axes(gwp_vals))
+    for i in eachindex(gwp_vals)
+        val1 = gwp_vals[i].divergence
+        val2 = zero(T)
+        for sh in divfns[i]
+            val2 += sh.coeff * lgx_vals[sh.refid].value
+        end
+        err[i] = abs(val1 - val2)
+    end
+    atol = sqrt(eps(T))
+    @test all(err .< atol)
 end

src/bases/local/laglocal.jl

@@ -386,7 +386,7 @@ function interpolate(fields, interpolant::LagrangeRefSpace{T,Degree,3}, chart) w
     s = range(0,1,length=Degree+1)
     Is = zip(I,s)
     idx = 1
-    vals = []
+    vals = Vector{Vector{T}}()
     for (i,ui) in Is
         for (j,vj) in Is
             for (k,wk) in Is
@@ -397,6 +397,10 @@ function interpolate(fields, interpolant::LagrangeRefSpace{T,Degree,3}, chart) w
                 idx += 1
     end end end
 
-    Q = hcat(vals...)
+    # Q = hcat(vals...)
+    Q = Matrix{T}(undef, length(vals[1]), length(vals))
+    for i in eachindex(vals)
+        Q[:,i] .= vals[i]
+    end
     return Q
 end

src/bases/local/ndlcdlocal.jl

@@ -78,7 +78,7 @@ function ttrace(x::NDLCDRefSpace, el, q, fc)
     return t
 end
 
-divergence(ref::NDLCDRefSpace, sh, el) = Shape(sh.cellid, 1, sh.coeff/volume(el))
+divergence(ref::NDLCDRefSpace, sh, el) = [Shape(sh.cellid, 1, sh.coeff/volume(el))]
 
 
 function restrict(ϕ::NDLCDRefSpace{T}, dom1, dom2) where {T}

src/bases/local/rtlocal.jl

@@ -28,7 +28,7 @@ end
 
 numfunctions(x::RTRefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 3
 
-divergence(ref::RTRefSpace, sh, el) = Shape(sh.cellid, 1, sh.coeff/volume(el))
+divergence(ref::RTRefSpace, sh, el) = [Shape(sh.cellid, 1, sh.coeff/volume(el))]
 
 """
     ntrace(refspace, element, localindex, face)