From eb7589de612e7a20540b446d9298d22839f65c71 Mon Sep 17 00:00:00 2001 From: Kristof Cools Date: Wed, 2 Oct 2024 17:16:41 +0200 Subject: [PATCH] self-interpolation verified --- src/bases/global/globaldofs.jl | 4 +- src/bases/local/gwplocal.jl | 90 +++++++++++++++++++++++++++++----- test/test_gwp.jl | 31 +++++++++++- 3 files changed, 111 insertions(+), 14 deletions(-) diff --git a/src/bases/global/globaldofs.jl b/src/bases/global/globaldofs.jl index b736bb1c..ad7908bd 100644 --- a/src/bases/global/globaldofs.jl +++ b/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 diff --git a/src/bases/local/gwplocal.jl b/src/bases/local/gwplocal.jl index a0788f4a..6dfa4c7c 100644 --- a/src/bases/local/gwplocal.jl +++ b/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,11 +62,18 @@ 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) @@ -74,3 +81,64 @@ function (ϕ::GWPCurlRefSpace{T,Deg})(dom::CompScienceMeshes.ReferenceSimplex{Di 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 \ No newline at end of file diff --git a/test/test_gwp.jl b/test/test_gwp.jl index 7d1b83b7..bc47712f 100644 --- a/test/test_gwp.jl +++ b/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 \ No newline at end of file