From 1813773501794f45b33dac89f4cb9b734e3c6b21 Mon Sep 17 00:00:00 2001 From: Kristof Cools Date: Mon, 7 Oct 2024 15:24:07 +0200 Subject: [PATCH] divergence tested and efie seems to work --- examples/efie_gwp2.jl | 31 ++++++ examples/mfie_gwp2.jl | 30 ++++++ src/BEAST.jl | 4 + src/bases/global/gwpdivglobal.jl | 39 +++++++ src/bases/global/gwpglobal.jl | 179 +++++++++++++++++++++++++++++++ src/bases/lagrange.jl | 3 +- src/bases/local/gwpdivlocal.jl | 69 ++++++++++++ src/bases/local/gwplocal.jl | 82 +++++++++++--- src/bases/local/rtqlocal.jl | 1 + src/bases/rtqspace.jl | 3 +- src/excitation.jl | 2 +- src/gridfunction.jl | 5 +- src/postproc.jl | 6 +- src/utils/lagpolys.jl | 6 +- test/test_gwp.jl | 48 +++++++-- 15 files changed, 478 insertions(+), 30 deletions(-) create mode 100644 examples/efie_gwp2.jl create mode 100644 examples/mfie_gwp2.jl create mode 100644 src/bases/global/gwpdivglobal.jl create mode 100644 src/bases/global/gwpglobal.jl create mode 100644 src/bases/local/gwpdivlocal.jl diff --git a/examples/efie_gwp2.jl b/examples/efie_gwp2.jl new file mode 100644 index 00000000..b343e991 --- /dev/null +++ b/examples/efie_gwp2.jl @@ -0,0 +1,31 @@ +using CompScienceMeshes +using BEAST + +Γ = readmesh(joinpath(dirname(pathof(BEAST)),"../examples/sphere2.in")) +# Γ = meshsphere(radius=1.0, h=0.1) +@show length(Γ) +# X = raviartthomas(Γ) +X = BEAST.gwpdiv(Γ; order=2) + +κ, η = 1.0, 1.0 +t = Maxwell3D.singlelayer(wavenumber=κ) +E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ) +e = (n × E) × n + +BEAST.@defaultquadstrat (t,X,X) BEAST.DoubleNumSauterQstrat(7,8,6,6,6,6) +BEAST.@defaultquadstrat (e,X) BEAST.SingleNumQStrat(12) + +@hilbertspace j +@hilbertspace k +# efie = @discretise t[k,j]==e[k] j∈X k∈X +# u, ch = BEAST.gmres_ch(efie; restart=1500) + +A = assemble(t[k,j], X, X) +b = assemble(e[k], X) +Ai = BEAST.GMRESSolver(A; reltol=1e-8, restart=1500) +u = Ai * b + +include("utils/postproc.jl") +include("utils/plotresults.jl") + +# Id = BEAST.Identity() \ No newline at end of file diff --git a/examples/mfie_gwp2.jl b/examples/mfie_gwp2.jl new file mode 100644 index 00000000..e235f38b --- /dev/null +++ b/examples/mfie_gwp2.jl @@ -0,0 +1,30 @@ +using CompScienceMeshes, BEAST + +# Γ = readmesh(joinpath(dirname(pathof(BEAST)),"../examples/sphere2.in")) +Γ = meshsphere(radius=1.0, h=0.1) +# X, Y = raviartthomas(Γ), buffachristiansen(Γ) +# X = brezzidouglasmarini(Γ) +# Y = brezzidouglasmarini(Γ) +# X = raviartthomas(Γ) +# Y = raviartthomas(Γ) +X = BEAST.gwpdiv(Γ; order=2) +Y = BEAST.gwpdiv(Γ; order=2) + +ϵ, μ, ω = 1.0, 1.0, 1.0; κ = ω * √(ϵ*μ) +# κ = 3.0 +NK, Id = BEAST.DoubleLayerRotatedMW3D(1.0, im*κ), Identity() +E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ) +H = -1/(im*μ*ω)*curl(E) +h = n × H + +BEAST.@defaultquadstrat (NK,X,X) BEAST.DoubleNumSauterQstrat(7,8,6,6,6,6) +BEAST.@defaultquadstrat (Id,X,X) BEAST.SingleNumQStrat(7) +BEAST.@defaultquadstrat (h,X) BEAST.SingleNumQStrat(12) + +@hilbertspace j +@hilbertspace m +mfie = @discretise (NK-0.5Id)[m,j] == h[m] j∈X m∈Y +u = gmres(mfie) + +include("utils/postproc.jl") +include("utils/plotresults.jl") diff --git a/src/BEAST.jl b/src/BEAST.jl index 96d4fc5d..879a8e94 100644 --- a/src/BEAST.jl +++ b/src/BEAST.jl @@ -156,6 +156,8 @@ include("bases/local/ndlcdlocal.jl") include("bases/local/bdm3dlocal.jl") include("bases/local/rtqlocal.jl") include("bases/local/gwplocal.jl") +include("bases/local/gwpdivlocal.jl") + include("bases/basis.jl") include("bases/lincomb.jl") @@ -163,6 +165,8 @@ include("bases/trace.jl") include("bases/restrict.jl") include("bases/divergence.jl") include("bases/global/globaldofs.jl") +include("bases/global/gwpglobal.jl") +include("bases/global/gwpdivglobal.jl") include("bases/lagrange.jl") include("bases/rtspace.jl") diff --git a/src/bases/global/gwpdivglobal.jl b/src/bases/global/gwpdivglobal.jl new file mode 100644 index 00000000..1b97f43d --- /dev/null +++ b/src/bases/global/gwpdivglobal.jl @@ -0,0 +1,39 @@ +struct GWPDivSpace{T,M,P} <: Space{T} + geo::M + fns::Vector{Vector{Shape{T}}} + pos::P + degree::Int +end + +function refspace(s::GWPDivSpace{T}) where {T} GWPDivRefSpace{T,s.degree}() end +function subset(s::S,I) where {S<:GWPDivSpace} S(s.geo, s.fns[I], s.pos[I], s.degree) end + +function gwpdiv(mesh, edges=nothing; order) + + T = coordtype(mesh) + S = Shape{T} + + space = BEAST.gwpcurl(mesh, edges; order) + fns = Vector{Vector{S}}(undef, length(space.fns)) + for (i,fn) in enumerate(space.fns) + fns[i] = [S(s.cellid, s.refid, -s.coeff) for s in fn] + end + + return GWPDivSpace(space.geo, fns, space.pos, order) +end + +@testitem "GWPcurl global: numfunctions" begin + using CompScienceMeshes + + h = 0.5 + mesh = meshrectangle(1.0, 1.0, 0.5) + edges = setminus(skeleton(mesh,1), boundary(mesh)) + + order = 2 + gwp = BEAST.gwpdiv(mesh, edges; order=order) + + ne = order+1 + nf = order * (order+1) + Nt = length(edges)*ne + length(mesh)*nf + @test numfunctions(gwp) == Nt +end \ No newline at end of file diff --git a/src/bases/global/gwpglobal.jl b/src/bases/global/gwpglobal.jl new file mode 100644 index 00000000..3270e4da --- /dev/null +++ b/src/bases/global/gwpglobal.jl @@ -0,0 +1,179 @@ +struct GWPCurlSpace{T,M,P} <: Space{T} + geo::M + fns::Vector{Vector{Shape{T}}} + pos::P + degree::Int +end + +function refspace(s::GWPCurlSpace{T}) where {T} GWPCurlRefSpace{T,s.degree}() end +function subset(s::S,I) where {S<:GWPCurlSpace} S(s.geo, s.fns[I], s.pos[I], s.degree) end + +struct GWPGlobalEdgeDoFs order::Int end +struct GWPGlobalFaceDoFs order::Int end + +function globaldofs(edge_ch, face_ch, localspace, dof::GWPGlobalEdgeDoFs) + f = trace(edge_ch, face_ch, localspace) + T = coordtype(edge_ch) + ds = one(T) / (dof.order+2) + return stack(range(ds, step=ds, length=dof.order+1)) do s + p = neighborhood(edge_ch, s) + t = tangents(p, 1) + [dot(t,x.value) for x in f(s)] + end +end + +@testitem "GWP edge global dofs" begin + using CompScienceMeshes + T = Float64 + face_ch = simplex( + point(3,0,0), + point(2,0,-1), + point(0,0,1), + ) + edge_ch = simplex( + point(0,0,1), + point(2,0,-1), + ) + + degree = 3 + localspace = BEAST.GWPCurlRefSpace{T,degree}() + dofs = BEAST.GWPGlobalEdgeDoFs(degree) + A = BEAST.globaldofs(edge_ch, face_ch, localspace, dofs) + display(round.(A, digits=3)) +end + +function globaldofs(edge_ch, face_ch, localspace, dof::GWPGlobalFaceDoFs) + + d = dof.order + T = coordtype(edge_ch) + ds = one(T) / (d+2) + + f = trace(edge_ch, face_ch, localspace) + Q = zeros(T, numfunctions(localspace, domain(face_ch)), d*(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) + # return stack(S) do (i,j,k) + idx = 1 + for (i,j,k) in S + u = (i*ds,j*ds) + p = neighborhood(edge_ch, u) + t1 = tangents(p,1) + t2 = tangents(p,2) + vals = f(u) + q1 = [dot(t1,x.value) for x in vals] + q2 = [dot(t2,x.value) for x in vals] + Q[:,idx] = q1; idx += 1 + Q[:,idx] = q2; idx += 1 + end + + return Q +end + + +@testitem "GWP face global dofs" begin + using CompScienceMeshes + T = Float64 + face_ch = simplex( + point(3,0,0), + point(2,0,-1), + point(0,0,1), + ) + # edge_ch = simplex( + # point(0,0,1), + # point(2,0,-1), + # point(3,0,0) + # ) + edge_ch = simplex( + point(3,0,0), + point(2,0,-1), + point(0,0,1), + ) + + degree = 3 + localspace = BEAST.GWPCurlRefSpace{T,degree}() + dofs = BEAST.GWPGlobalFaceDoFs(degree) + A = BEAST.globaldofs(edge_ch, face_ch, localspace, dofs) + display(round.(A, digits=3)) +end + +function _addshapes!(fns, cell, gids, lids, β) + T = eltype(β) + atol = sqrt(eps(T)) + for i in axes(β,1) + fn = fns[gids[i]] + for j in axes(β,2) + isapprox(β[i,j], 0; atol) && continue + push!(fn, Shape{T}(cell, lids[j], β[i,j])) + end + end +end + + +function gwpcurl(mesh, edges=nothing; order) + + T = coordtype(mesh) + atol = sqrt(eps(T)) + + if edges == nothing + edges = setminus(skeleton(mesh, 1), boundary(mesh)) + end + + C12 = connectivity(mesh, edges, abs) + R12 = rowvals(C12) + V12 = nonzeros(C12) + + ne = order+1 + nf = order*(order+1) + + Ne = ne * length(edges) + Nf = nf * length(mesh) + Nt = Ne + Nf + + localspace = GWPCurlRefSpace{T,order}() + refdom = domain(chart(mesh, first(mesh))) + localdim = numfunctions(localspace, refdom) + + fns = [Shape{T}[] for n in 1:Nt] + pos = fill(point(0,0,0), Nt) + for cell in mesh + cell_ch = chart(mesh, cell) + for k in nzrange(C12, cell) + e = R12[k] + i = V12[k] + gids = (e-1)*ne .+ (1:ne) + lids = localindices(localspace, refdom, Val{1}, i) + edge_ch = chart(edges, e) + vals = globaldofs(edge_ch, cell_ch, localspace, GWPGlobalEdgeDoFs(order)) + α = vals[lids,:] + _addshapes!(fns, cell, gids, lids, inv(α')) + end + + order < 2 && continue + face_ch = chart(mesh, cell) + gids = Ne + (cell-1)*nf .+ (1:nf) + lids = localindices(localspace, refdom, Val{2}, 1) + vals = globaldofs(face_ch, face_ch, localspace, GWPGlobalFaceDoFs(order)) + α = vals[lids,:] + _addshapes!(fns, cell, gids, lids, inv(α')) + end + + for e in edges pos[(e-1)*ne .+ (1:ne)] .= Ref(cartesian(center(chart(edges, e)))) end + for f in mesh pos[Ne+(f-1)*nf .+ (1: nf)] .= Ref(cartesian(center(chart(mesh, f)))) end + + return GWPCurlSpace(mesh, fns, pos, order) +end + +@testitem "GWPcurl global: numfunctions" begin + using CompScienceMeshes + + h = 0.5 + mesh = meshrectangle(1.0, 1.0, 0.5) + edges = setminus(skeleton(mesh,1), boundary(mesh)) + + order = 2 + gwp = BEAST.gwpcurl(mesh, edges; order=order) + + ne = order+1 + nf = order * (order+1) + Nt = length(edges)*ne + length(mesh)*nf + @test numfunctions(gwp) == Nt +end \ No newline at end of file diff --git a/src/bases/lagrange.jl b/src/bases/lagrange.jl index a0bbb336..77297b5d 100644 --- a/src/bases/lagrange.jl +++ b/src/bases/lagrange.jl @@ -977,7 +977,8 @@ function lagrangec0(mesh::CompScienceMeshes.AbstractMesh{<:Any,3}; order) N = nV + nE + nF localspace = LagrangeRefSpace{T,order,3,binomial(2+order,2)}() - localdim = numfunctions(localspace) + dom = domain(chart(mesh, first(mesh))) + localdim = numfunctions(localspace, dom) d = order fns = [S[] for n in 1:(nV+nE+nF)] diff --git a/src/bases/local/gwpdivlocal.jl b/src/bases/local/gwpdivlocal.jl new file mode 100644 index 00000000..8c21b9c4 --- /dev/null +++ b/src/bases/local/gwpdivlocal.jl @@ -0,0 +1,69 @@ +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 (ϕ::GWPDivRefSpace{T,Degree})(p) where {T,Degree} + ψ = GWPCurlRefSpace{T,Degree}() + dom = domain(chart(p)) + u = parametric(p) + vals = ψ(dom, u) + vals = pushforwardcurl(vals, p) + n = normal(p) + map(vals) do v + (value = -n × v.value, divergence=v.curl) + end +end + +function interpolate(fields, interpolant::GWPDivRefSpace{T,Degree}, chart) where {T,Degree} + + function nxfields(p) + vals = fields(p) + n = normal(p) + return map(v -> n×v, vals) + end + + return interpolate(nxfields, GWPCurlRefSpace{T,Degree}(), chart) +end + +@testitem "divergence" begin + using CompScienceMeshes, LinearAlgebra + + T = Float64 + s = simplex( + point(3,0,0), + point(0,2,1), + point(-1,-1,-1), + ) + + order = 2 + ϕ = BEAST.GWPDivRefSpace{T,order}() + + u = T(0.2432); dx = sqrt(eps(T)) + v = T(0.5786); dy = sqrt(eps(T)) + + # p00 = neighborhood(s, (u,v)) + # p10 = neighborhood(s, (u+du,v)) + # p01 = neighborhood(s, (u, v+dv)) + + p00 = neighborhood(s, (u,v)) + t1 = normalize(tangents(p00,1)) + t2 = normal(p00) × t1 + + p10 = neighborhood(s, carttobary(s, cartesian(p00) + dx*t1 + 0*t2)) + p01 = neighborhood(s, carttobary(s, cartesian(p00) + 0*t1 + dy*t2)) + + ϕ00 = ϕ(p00) + ϕ10 = ϕ(p10) + ϕ01 = ϕ(p01) + + + for (f00, f10, f01) in zip(ϕ00, ϕ10, ϕ01) + df1dx = dot(t1, f10.value - f00.value) / dx + df2dy = dot(t2, f01.value - f00.value) / dy + curl_num = df1dx + df2dy + curl_ana = f00.divergence + # @show curl_num, curl_ana, abs(curl_num - curl_ana) + @test curl_num ≈ curl_ana atol=sqrt(eps(T))*100 + end +end \ No newline at end of file diff --git a/src/bases/local/gwplocal.jl b/src/bases/local/gwplocal.jl index 6dfa4c7c..1dc4dd4d 100644 --- a/src/bases/local/gwplocal.jl +++ b/src/bases/local/gwplocal.jl @@ -1,5 +1,8 @@ struct GWPCurlRefSpace{T,Degree} <: RefSpace{T} end +function numfunctions(x::GWPCurlRefSpace{<:Any,D}, + dom::CompScienceMeshes.ReferenceSimplex{2}) where{D} (D+1)*(D+3) end + function (ϕ::GWPCurlRefSpace{T,Degree})(p) where {T,Degree} dom = domain(chart(p)) u = parametric(p) @@ -14,18 +17,15 @@ function (ϕ::GWPCurlRefSpace{T,Deg})(dom::CompScienceMeshes.ReferenceSimplex{Di 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) + s = range(zero(T), one(T), length=d+3) + 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[] + crls = T[] i = 0 for j in 1:d+1 @@ -35,7 +35,13 @@ function (ϕ::GWPCurlRefSpace{T,Deg})(dom::CompScienceMeshes.ReferenceSimplex{Di Rₖ = _sylpoly_shift(s, k+1, w) push!(vals, Rᵢ*Rⱼ*Rₖ*nd1) - # du = _sylpoly_shift_diff(s, i+1, u)*Rⱼ*Rₖ + dRᵢ = _sylpoly_diff(s, i+1, u) + dRⱼ = _sylpoly_shift_diff(s, j+1, v) + dRₖ = _sylpoly_shift_diff(s, k+1, w) + du = dRᵢ*Rⱼ*Rₖ - Rᵢ*Rⱼ*dRₖ + dv = Rᵢ*dRⱼ*Rₖ - Rᵢ*Rⱼ*dRₖ + curl = (du*nd1[2] - dv*nd1[1]) + 2*Rᵢ*Rⱼ*Rₖ + push!(crls, curl) end for i in 1:d+1 @@ -45,6 +51,14 @@ function (ϕ::GWPCurlRefSpace{T,Deg})(dom::CompScienceMeshes.ReferenceSimplex{Di Rⱼ = _sylpoly(s, j+1, v) Rₖ = _sylpoly_shift(s, k+1, w) push!(vals, Rᵢ*Rⱼ*Rₖ*nd2) + + dRᵢ = _sylpoly_shift_diff(s, i+1, u) + dRⱼ = _sylpoly_diff(s, j+1, v) + dRₖ = _sylpoly_shift_diff(s, k+1, w) + du = dRᵢ*Rⱼ*Rₖ - Rᵢ*Rⱼ*dRₖ + dv = Rᵢ*dRⱼ*Rₖ - Rᵢ*Rⱼ*dRₖ + curl = (du*nd2[2] - dv*nd2[1]) + 2*Rᵢ*Rⱼ*Rₖ + push!(crls, curl) end for i in 1:d+1 @@ -54,6 +68,16 @@ function (ϕ::GWPCurlRefSpace{T,Deg})(dom::CompScienceMeshes.ReferenceSimplex{Di Rⱼ = _sylpoly_shift(s, j+1, v) Rₖ = _sylpoly(s, k+1, w) push!(vals, Rᵢ*Rⱼ*Rₖ*nd3) + + dRᵢ = _sylpoly_shift_diff(s, i+1, u) + dRⱼ = _sylpoly_shift_diff(s, j+1, v) + dRₖ = _sylpoly_diff(s, k+1, w) + + du = dRᵢ*Rⱼ*Rₖ - Rᵢ*Rⱼ*dRₖ + dv = Rᵢ*dRⱼ*Rₖ - Rᵢ*Rⱼ*dRₖ + curl = (du*nd3[2] - dv*nd3[1]) + 2*Rᵢ*Rⱼ*Rₖ + + push!(crls, curl) end for i in 1:d+1 @@ -69,15 +93,34 @@ function (ϕ::GWPCurlRefSpace{T,Deg})(dom::CompScienceMeshes.ReferenceSimplex{Di 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) + + dRsᵢ = _sylpoly_shift_diff(s, i+1, u) + dRsⱼ = _sylpoly_shift_diff(s, j+1, v) + dRsₖ = _sylpoly_shift_diff(s, k+1, w) + dRᵢ = _sylpoly_diff(s, i+1, u) + dRⱼ = _sylpoly_diff(s, j+1, v) + dRₖ = _sylpoly_diff(s, k+1, w) + + du = dRᵢ*Rsⱼ*Rsₖ - Rᵢ*Rsⱼ*dRsₖ + dv = Rᵢ*dRsⱼ*Rsₖ - Rᵢ*Rsⱼ*dRsₖ + curlS1 = du*nd1[2] - dv*nd1[1] + 2*Rᵢ*Rsⱼ*Rsₖ + + du = dRsᵢ*Rⱼ*Rsₖ - Rsᵢ*Rⱼ*dRsₖ + dv = Rsᵢ*dRⱼ*Rsₖ - Rsᵢ*Rⱼ*dRsₖ + curlS2 = du*nd2[2] - dv*nd2[1] + 2*Rsᵢ*Rⱼ*Rsₖ + + du = dRsᵢ*Rsⱼ*Rₖ - Rsᵢ*Rsⱼ*dRₖ + dv = Rsᵢ*dRsⱼ*Rₖ - Rsᵢ*Rsⱼ*dRₖ + curlS3 = du*nd3[2] - dv*nd3[1] + 2*Rsᵢ*Rsⱼ*Rₖ + + push!(crls, curlS2 - curlS3) + push!(crls, curlS3 - curlS1) end end NF = length(vals) - SVector{NF}([(value=f, curl=zero(T)) for f in vals]) + SVector{NF}([(value=f, curl=c) for (f,c) in zip(vals, crls)]) end @@ -141,4 +184,19 @@ function interpolate(fields, interpolant::GWPCurlRefSpace{T,Degree}, chart) wher end end return Q -end \ No newline at end of file +end + +function localindices(localspace::GWPCurlRefSpace{<:Any,Degree}, domain, + dim::Type{Val{1}}, i) where {Degree} + + ne = Degree+1 + (i-1)*ne .+ (1:ne) +end + +function localindices(localspace::GWPCurlRefSpace{<:Any,Degree}, domain, + dim::Type{Val{2}}, i) where {Degree} + + ne = Degree+1 + nf = Degree * (Degree + 1) + 3*ne .+ (1:nf) +end diff --git a/src/bases/local/rtqlocal.jl b/src/bases/local/rtqlocal.jl index f7e38347..f1472d52 100644 --- a/src/bases/local/rtqlocal.jl +++ b/src/bases/local/rtqlocal.jl @@ -22,6 +22,7 @@ function (ϕ::RTQuadRefSpace{T})(p) where {T} (value=D*f.value/j, divergence=f.divergence/j) end end +function numfunctions(ϕ::RTQuadRefSpace, dom::CompScienceMeshes.RefQuadrilateral) 4 end function interpolate(fields, interpolant::RTQuadRefSpace{T}, chart) where {T} diff --git a/src/bases/rtqspace.jl b/src/bases/rtqspace.jl index 2935c223..7f701ee6 100644 --- a/src/bases/rtqspace.jl +++ b/src/bases/rtqspace.jl @@ -97,7 +97,8 @@ end num_cells = numcells(m) num_bfs = numfunctions(s) r = refspace(s) - num_refs = numfunctions(r) + dom = domain(chart(m, first(m))) + num_refs = numfunctions(r, dom) celltonum = BEAST.make_celltonum(num_cells, num_refs, num_bfs, s) els, ad, a2g = BEAST.assemblydata(s) diff --git a/src/excitation.jl b/src/excitation.jl index f194a86c..17bc17f0 100644 --- a/src/excitation.jl +++ b/src/excitation.jl @@ -86,7 +86,7 @@ function assemble!(field::Functional, tfs::subdBasis, store; end -function celltestvalues(tshs::RefSpace{T}, tcell, field, qr) where {T,NF} +function celltestvalues(tshs::RefSpace{T}, tcell, field, qr) where {T} num_tshs = numfunctions(tshs, domain(tcell)) interactions = zeros(Complex{T}, num_tshs) diff --git a/src/gridfunction.jl b/src/gridfunction.jl index 05da678f..56931843 100644 --- a/src/gridfunction.jl +++ b/src/gridfunction.jl @@ -279,7 +279,7 @@ function cellquadsamplingvaluesvalues!( qr ) - num_tshs = numfunctions(tshs) + # num_tshs = numfunctions(tshs) num_oqp = length(qr) @@ -289,7 +289,8 @@ function cellquadsamplingvaluesvalues!( tvals = qr[p].value - for m in 1 : num_tshs + for m in 1 : length(tvals) + # for m in 1 : num_tshs tval = tvals[m] for k = 1:size(data, 1) diff --git a/src/postproc.jl b/src/postproc.jl index 398f4c8f..b720a9b5 100644 --- a/src/postproc.jl +++ b/src/postproc.jl @@ -10,12 +10,14 @@ function facecurrents(coeffs, basis) T = eltype(coeffs) RT = real(T) + mesh = geometry(basis) + dom = domain(chart(mesh, first(mesh))) + refs = refspace(basis) - numrefs = numfunctions(refs) + numrefs = numfunctions(refs, dom) cells, tad, a2g = assemblydata(basis) - mesh = geometry(basis) D = dimension(mesh) # U = D+1 U = 3 diff --git a/src/utils/lagpolys.jl b/src/utils/lagpolys.jl index 65d3c924..60620706 100644 --- a/src/utils/lagpolys.jl +++ b/src/utils/lagpolys.jl @@ -11,15 +11,17 @@ function _lagpoly(nodes, i, s, i0=1, i1=length(nodes)) end function _lagpoly_diff(nodes, i, s, i0=1, i1=length(nodes)) + T = typeof(s) r = zero(T) si = nodes[i] for p in i0:i1 p == i && continue + sp = nodes[p] rp = one(T) for j in i0:i1 j == i && continue j == p && continue - sj = nodees[j] + sj = nodes[j] rp *= (s - sj) / (si - sj) end rp *= 1 / (si - sp) @@ -41,5 +43,5 @@ function _sylpoly_shift(nodes, i, s) end function _sylpoly_shift_diff(nodes, i, s) - _lagpoly_shift(nodes, i, s, 2, i) + _lagpoly_diff(nodes, i, s, 2, i) end \ No newline at end of file diff --git a/test/test_gwp.jl b/test/test_gwp.jl index bc47712f..ca2210e8 100644 --- a/test/test_gwp.jl +++ b/test/test_gwp.jl @@ -17,6 +17,7 @@ end @testitem "refspace: self-interpolate" begin using CompScienceMeshes + using LinearAlgebra T, Dim, Nverts, Degree = Float64, 2, 3, 2 dom = CompScienceMeshes.ReferenceSimplex{Dim,T,Nverts}() @@ -26,16 +27,45 @@ end # fields(p) = [ϕ(dom,p)[1].value] coeffs = BEAST.interpolate(fields, ϕ, dom) + nf = numfunctions(ϕ, dom) + @test coeffs ≈ Matrix{T}(I,nf,nf) atol=sqrt(eps(T)) display(round.(coeffs, digits=3)) - # pts = [ - # point(T, 0.3, 0.1), - # point(T, 0.1, 0.3), - # point(T, 0.5, 0.0)] +end + + +@testitem "confspace: interpolate lowest order" begin + using CompScienceMeshes + using LinearAlgebra + + T, Dim, Nverts, Degree = Float64, 2, 3, 2 + supp = simplex( + point(T,3,0,0), + point(T,0,2,0), + point(T,0,0,-1)) + + ϕ = BEAST.GWPCurlRefSpace{T,Degree}() + ψ = BEAST.NDRefSpace{T}() + + fields(p) = [v.value for v in ψ(p)] + coeffs = BEAST.interpolate(fields, ϕ, supp) + + nf1 = numfunctions(ψ, domain(supp)) + nf2 = numfunctions(ϕ, domain(supp)) + @test size(coeffs) == (nf1, nf2) + 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 + for (i,u) in enumerate(pts) + p = neighborhood(supp, u) + basisvals = ϕ(p) + fieldvals = ψ(p) + r = sum(c * v.value for (c,v) in zip(coeffs[i,:], basisvals)) + s = fieldvals[i].value + @test r ≈ s atol=sqrt(eps(T)) + end end \ No newline at end of file