From fca8afe999497f922ec6adb634023f3d1e1c4622 Mon Sep 17 00:00:00 2001 From: Kristof Cools Date: Fri, 11 Oct 2024 17:18:40 +0200 Subject: [PATCH] support to take the divergence of GWP space --- examples/divgwp_expansion.jl | 66 +++++++++++++++++++++++++++++ src/bases/divergence.jl | 5 ++- src/bases/global/gwpdivglobal.jl | 58 +++++++++++++++++++++++++ src/bases/lagrange.jl | 1 + src/bases/local/bdmlocal.jl | 2 +- src/bases/local/gwpdivlocal.jl | 72 +++++++++++++++++++++++++++++++- src/bases/local/laglocal.jl | 8 +++- src/bases/local/ndlcdlocal.jl | 2 +- src/bases/local/rtlocal.jl | 2 +- 9 files changed, 208 insertions(+), 8 deletions(-) create mode 100644 examples/divgwp_expansion.jl diff --git a/examples/divgwp_expansion.jl b/examples/divgwp_expansion.jl new file mode 100644 index 00000000..c6ed125d --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/bases/divergence.jl b/src/bases/divergence.jl index 12943376..a4b32add 100644 --- a/src/bases/divergence.jl +++ b/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) diff --git a/src/bases/global/gwpdivglobal.jl b/src/bases/global/gwpdivglobal.jl index 1b97f43d..f707b805 100644 --- a/src/bases/global/gwpdivglobal.jl +++ b/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 \ No newline at end of file diff --git a/src/bases/lagrange.jl b/src/bases/lagrange.jl index 77297b5d..cbc3a617 100644 --- a/src/bases/lagrange.jl +++ b/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}}} diff --git a/src/bases/local/bdmlocal.jl b/src/bases/local/bdmlocal.jl index 1eb173a4..1c83f5f7 100644 --- a/src/bases/local/bdmlocal.jl +++ b/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), diff --git a/src/bases/local/gwpdivlocal.jl b/src/bases/local/gwpdivlocal.jl index 8c21b9c4..342c531a 100644 --- a/src/bases/local/gwpdivlocal.jl +++ b/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 \ No newline at end of file diff --git a/src/bases/local/laglocal.jl b/src/bases/local/laglocal.jl index cd90f2c3..21355e0f 100644 --- a/src/bases/local/laglocal.jl +++ b/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 \ No newline at end of file diff --git a/src/bases/local/ndlcdlocal.jl b/src/bases/local/ndlcdlocal.jl index 1ecc2ad6..8b924607 100644 --- a/src/bases/local/ndlcdlocal.jl +++ b/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} diff --git a/src/bases/local/rtlocal.jl b/src/bases/local/rtlocal.jl index 7c16d6f8..c9350014 100644 --- a/src/bases/local/rtlocal.jl +++ b/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)